Stresses on fault computed by Pylith compared to the analytical solution

PyLith
1

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?

pylvsanalytical
pylvsanalytical1848×961 44.6 KB

simulation_compare.zip (72.4 KB)

2

If you use a basis order of 1 for the displacement field, then the stresses in the cells will be uniform (basis order of 0). Thus, the stress field is not continuous and you will see issues like you describe. If you use a basis order of 2 for the displacement field, then the stress field will have a basis order of 1 (linear within cells) and nearly continuous acriss cells (in the finite-element formulation, we do not impose C1 continuity across cells).

If you want the tractions on the fault, then you should output the Lagrange multiplier (lagrange_multiplier_fault) as it directly corresponds to the tractions on the fault in the model coordinate system (not the fault coordinate system). This feature is available in the PyLith main branch. See tests/fullscale/linearelasticity/faults-2d/pylithapp.cfg for an example of outputting the Lagrange multiplier.

In the next PyLith release, we plan to have be able to output the fault traction in the fault coordinate system.

1 Like