Dear Pylith developers,

I am using Pylith V3.0.2 and my question is about the stresses calculated on a left-lateral, 30 deg fault, based on the Green’s Function example (see attachment).

Our goal is to use Pylith to calculate the stresses on the fault for a heterogeneous earth, which we will then use in our own BEM code to simulate earthquake cycles (we can only have faults in a homogeneous earth right now). Thus, as a simple test, I ran a simple Pylith simulation, and then copied the same settings (geometry, constant Vp, constant Vs, constant density, etc), into the set of analytical solutions for a fault in a half space that I got from Paul Segall’s ‘Earthquake and Volcano Deformation’ book (eq’s 3.77-3.79). I have attached a figure with the shear stress for a single source on the fault.

I already noticed that Pylith gives the stresses on the positive and negative side of the fault, and both are different. I guess this is due to some interpolation of stresses from the center of a triangle to the vertices on the fault? But then the shape is also different for the two sides, as in, the negative side shows a stress that I would expect due to a single dislocation, but the stress on the positive side is a bit different, why is this?

If I then compare both solutions with that of the analytical solution (which are computed in fault local coordinates), I do not see a lot of overlap.

I am probably missing something here, but I can not figure out what it is that gives these differences. Could you help me?

simulation_compare.zip (72.4 KB)