Simulate 2 earthquakes with overlap fault

Dear Developers,

I have some questions. I am sorry if this is a basic question. I would like to ask whether it is possible to simulate 2 earthquakes event in a shorter time (for example: following the example in 3d subduction step04). After the first earthquake, the second earthquake happens several years later. If it is possible, can we say that the second earthquake will be affected by the stress from the first earthquake from this simulation? And if we can do that, is it also possible to do the simulation if the fault from the first earthquake overlaps the fault from the second earthquake? I attach the model below:


The purple patch refers to the 1st earthquake fault patch. The yellow patch refers to the 2nd earthquake fault patch. And the green refers to the overlap fault patch which is actually the part of the first EQ but also the second EQ.
I tried to do the simulation using both faults but I got an error for the second fault patch that says:

Ambiguous fault surface. Cell 1174454 has all of its vertices on the fault

I have tried to run the simulation for each fault patch but it works fine.

Thank you in advance.
Alvina KK

PyLith does not allow overlapping faults. In this case, you should create a single large fault that encompasses everywhere slip will occur over all the earthquakes. If you are prescribing slip, then use spatial databases to limit the extent of slip to just the patch you want. For spontaneous rupture, you can confine slip by adjusting the friction, shear traction, and/or normal tractions using the spatial databases.

Dear @baagaard ,
I am currently trying to follow your suggestion but getting some errors.

{default}::
– pyre.inventory(error)
– timedependent.interfaces.faultcohesivekin.eq_srcs.eqkinsrc.stepslipfn.simpledb.simpleioascii.filename ← ‘’
– Filename for spatial database not specified.
{default}::
– pyre.inventory(error)
– timedependent.interfaces.faultcohesivekin.eq_srcs.eqkinsrc.stepslipfn.simpledb.simpleioascii.filename ← ‘’
– Filename for spatial database not specified.
step02_spontan.cfg:124:
– pyre.inventory(error)
– timedependent.interfaces.faultcohesivekin.eq_srcs.eqkinsrc.stepslipfn.simpledb.filename ← ‘spatialdb/fault_slabtop_coseismic_merc_hill.spatialdb’
– unrecognized property ‘timedependent.interfaces.faultcohesivekin.eq_srcs.eqkinsrc.stepslipfn.simpledb.filename’

I would like to try first to run the simulation with 2 different source slip in different time before I try the spontaneous rupture. I only change the input for the second earthquake but I get this error.
Do I miss something in this part?
I also attach my cfg file.step02_spontan.cfg (6.9 KB)

Thank you in advance.

Best Regards,
Alvina KK

The SimpleDB object has an iohandler (SimpleIOAscii is the default) for I/O, so you need to use

slip.iohandler.filename = spatialdb/fault_slabtop_coseismic_merc_hill.spatialdb

This is why you get these two error messages.

It works fine now. Thank you,Brad!

Dear @baagaard ,

I tried to use static-dynamic friction and slip weakening for the second fault but the result of the slip is zero. Is it because the initial traction resulted from the first fault is too small? I add this part for the second fault :

Initial fault tractions

traction_perturbation = pylith.faults.TractPerturbation
traction_perturbation.db_initial = spatialdata.spatialdb.SimpleDB
traction_perturbation.db_initial.label = Initial fault tractions
traction_perturbation.db_initial.iohandler.filename = spatialdb/fault_slabtop_tractions_gusman.spatialdb

Friction

friction = pylith.friction.SlipWeakening
friction.label = Slip weakening

Force healing after each time step, so weakening is confined to each

time step and does not carry over into subsequent time steps.

friction.force_healing = True

We must define the quadrature information for fault cells.

The fault cells are 2D (surface).

quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 2

friction.db_properties = spatialdata.spatialdb.SimpleDB
friction.db_properties.label = Slip weakening
friction.db_properties.iohandler.filename = spatialdb/fault_slabtop_slipweakening_hill.spatialdb

Thank you in advance.
Regards,
Alvina KK

You can examine the tractions on the fault to assess how close to failure it is. For example, you can compute the ratio of the magnitude of the shear traction to the magnitude of the normal traction and compare this to the coefficient of friction.

Dear @baagaard ,
Thank you for your explanation. I got it. I am also wondering if it is possible to use this examination to see the state of failure for branching fault (example: splay fault)? So the geometry of the second fault is not at the same plane as the first fault plane but has a different dip angle. So the fault is rooted in the first fault plane.

Thank you in advance.
Best Regards,
Alvina KK

Yes, this same methodology works for other faults as well.

Dear @baagaard ,
Since I also want to analyse the stress distribution during the interseismic period, I run the interseismic period using this geometry but with the entire fault slabtop as the fault interface. But I get this error:

– Initializing integrators.
Fatal error. Calling MPI_Abort() to abort PyLith application.
Traceback (most recent call last):
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/apps/PetscApplication.py”, line 74, in onComputeNodes
self.main(*args, **kwds)
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/apps/PyLithApp.py”, line 128, in main
self.problem.initialize()
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/problems/TimeDependent.py”, line 120, in initialize
self.formulation.initialize(self.dimension, self.normalizer)
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/problems/Implicit.py”, line 121, in initialize
self._initialize(dimension, normalizer)
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/problems/Formulation.py”, line 470, in _initialize
integrator.initialize(totalTime, numTimeSteps, normalizer)
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/feassemble/ElasticityImplicit.py”, line 56, in initialize
ModuleElasticityImplicit.initialize(self, self.mesh())
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/feassemble/feassemble.py”, line 359, in initialize
def initialize(self, *args): return _feassemble.IntegratorElasticity_initialize(self, *args)
RuntimeError: Error while projecting location.
tolerance condition error
projection: merc
units: m
proj options:
lon: -0.923157
lat: 1.5708

I once got this kind of error as well and fixed by changing the mesh size from smaller to bigger size. I am wondering if this issue really comes from the mesh size or not. Also, I don’t understand why pylith suggested me to use this projection option (lon: -0.923157 lat: 1.5708) because this point is not in my geometry. Do you have any experience with this?

Thank you in advance.

Best Regards,
Alvina KK

This looks like a mesh or coordinate system issue. The latitude value corresponds to pi/2 or 90 degrees north. Make sure all of the vertices in your mesh have coordinate values in the range you expect and that your coordinate system for the mesh is consistent with the the coordinates of the vertices in your mesh.

Dear @baagaard ,
I am sorry for asking you this but how can I check that the vertices in the mesh are consistent with my coordinate system? I don’t understand why the vertices can have a different coordinate system since all the surfaces used to create the geometry are in the same coordinate system. Thus the vertices inside the geometry should also be the same.

Alvina KK

You create your mesh outside of PyLith, so there is a chance that the coordinate system you used when creating your mesh is not the same one specified in the PyLith parameter file. The easiest way to check what coordinate system PyLith is using is to add the command line argument --mesh_generator.reader.coordsys.help-properties when you run PyLith. This will print out the parameters for the coordinate system.

If you generate an Exodus file for your mesh, then you can examine the contents using ncdump. For example, ncdump -v coordx,coordy,coordz mesh.exo.

Dear @baagaard ,

I obtained this result when I add the command line argument:

properties of ‘csgeoproj’:
datum_horiz=: Name of horizontal datum.
default value: ‘WGS84’
current value: ‘WGS84’, from {file=‘./pylithapp.cfg’, line=59, column=-1}
datum_vert=: Name of vertical datum.
default value: ‘ellipsoid’
current value: ‘mean sea level’, from {file=‘./pylithapp.cfg’, line=60, column=-1}
ellipsoid=: Name of ellipsoid.
default value: ‘WGS84’
current value: ‘WGS84’, from {default}
help=: prints a screen that describes my traits
default value: False
current value: False, from {default}
help-components=: prints a screen that describes my subcomponents
default value: False
current value: False, from {default}
help-persistence=: prints a screen that describes my persistent store
default value: False
current value: False, from {default}
help-properties=: prints a screen that describes my properties
default value: False
current value: True, from {command line}
is_geocentric=: Flag indicating geocentric coordinate system.
default value: False
current value: False, from {default}
origin_lat=: Latitude of local origin in degrees (WGS84).
default value: 0.0
current value: 0.0, from {default}
origin_lon=: Longitude of local origin in degrees (WGS84).
default value: 0.0
current value: 0.0, from {default}
rotation_angle=: Rotation angle in degrees CCW of local x-axis from east.
default value: 0.0
current value: 0.0, from {default}
space_dim=: Number of dimensions for coordinate system.
default value: 3
current value: 3, from {file=‘./pylithapp.cfg’, line=58, column=-1}
units=: Units of coordinates.
default value: ‘m’
current value: ‘m’, from {default}

Which is the same coordinate system when I created the geometry. My mesh is an exodus file. If for example, there are some vertices which are read by Pylith as outside the geometry, how to fix it? I have no problem when I used smaller fault patch as seen in the figure to run the coseismic simulation. But the problem comes up with the entire fault slabtop.

Thank you for assistance.
Alvina KK

One of your earlier meshes had NaN values for some vertex coordinates. Is this happening again?

In general, the best way to debug these types of issues is to start with something that works and change one parameter at a time until you get to the desired set of parameters. If you encounter a problem, closely examine what you changed compared to what worked.

Dear @baagaard ,

I have fixed the previous NaN issue by recreate the surface to be used in the mesh geometry and there is no NaN values. With this mesh, I was able to run the coseismic simulation with small fault patch. That is why I continue to try the interseismic simulation and use the entire fault slabtop. But this issue came up.

I just tried to change the mesh size.


Previously, I used 2.5km for the long blue line (y-direction) and also the first point of a line at x-direction and varied up to 10km when it reached long red line (in the middle part of geometry). So I change the size from 2.5km (finest) to 10km (coarser) become 5km (finest) to 10km (coarser). With 5km as the finest mesh size, it can pass the error. It seems that this mesh variation becomes the issue? Because when I had this problem previoulsy, I used mesh size 2.5km (finest) at long blue line and 5km at the long red line. I fixed it by changing the mesh size at long red line to be 10km.

Alvina KK

Based on your description it looks like this is a problem with the mesh. I suggest checking the mesh quality and using condition number smoothing to improve it, if possible. You may need to play around with how you specify the mesh size on the geometric entities to get a better quality mesh with smooth variations in the mesh size.

Dear @baagaard ,
Yes it seems that the problem in my case is because I have two big main fault ( Now I have changed it as 1 big main fault) and also one small splay fault at shallow part. The mesh variation between these two faults trigger the issue. If my mesh size is too small for that main fault (for example 2.5km for fault 150x200km), it becomes a problem to run interseismic and spontaneous earthquake simulation but it is fine for coseismic simulation. Meanwhile the splay fault needs the small mesh size. Now, I get some senses related to this issue.

Now, I try again to run the simulation and got this notification which I dont understand:

RuntimeError: basic_string::substr: __pos (which is 11) > this->size() (which is 10)

What does this refer to? So, I will know things I have to check.

Thank you.
Alvina KK

Can you post the entire output? This error indicates something is wrong with a string but we need more information.

Dear @baagaard ,

Here is the report :

/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/faults/FaultCohesiveKin.py:160:initialize
– faultcohesivekin(info)
– Initializing fault ‘fault_slabtop_patch’.
Fatal error. Calling MPI_Abort() to abort PyLith application.
Traceback (most recent call last):
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/apps/PetscApplication.py”, line 74, in onComputeNodes
self.main(*args, **kwds)
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/apps/PyLithApp.py”, line 128, in main
self.problem.initialize()
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/problems/TimeDependent.py”, line 120, in initialize
self.formulation.initialize(self.dimension, self.normalizer)
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/problems/Implicit.py”, line 121, in initialize
self._initialize(dimension, normalizer)
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/problems/Formulation.py”, line 470, in _initialize
integrator.initialize(totalTime, numTimeSteps, normalizer)
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/faults/FaultCohesiveKin.py”, line 166, in initialize
FaultCohesive.initialize(self, totalTime, numTimeSteps, normalizer)
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/faults/Fault.py”, line 173, in initialize
self.output.writeInfo()
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/meshio/OutputManager.py”, line 222, in writeInfo
field = self.dataProvider().getVertexField(name)
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/faults/FaultCohesiveKin.py”, line 177, in getVertexField
field = ModuleFaultCohesiveKin.vertexField(self, name)
File “/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/faults/faults.py”, line 390, in vertexField
def vertexField(self, *args): return _faults.FaultCohesiveKin_vertexField(self, *args)
RuntimeError: basic_string::substr: __pos (which is 11) > this->size() (which is 10)
application called MPI_Abort(MPI_COMM_WORLD, -1) - process 0
/scratch/akuncoro/2021/software/pylith-2.2.2-linux-x86_64/bin/nemesis: mpirun: exit 255

and here the json file :
step02_s-parameters.json.tar.gz (18.8 KB)

Thank you in advance.
Alvina K K