Initial stresses, strains, fault slip and traction (use quasistatic output as initial condition for dynamic simulation)

Dear Pylith developers,

I’m modeling fault rupture in a 2D (plane strain) elastic block using FaultCohesiveDyn. I have a horizontal planar fault with three sections (middle - static friction, sides - slip-weakening or rate-and-state).
I would like to load my model quasistatically (lithostatic compression and shear) to accumulate slip on the middle section of the fault and overstress the two side section and then use the output to run the dynamic simulation.
I have my quasistatic results, I wrote separate databases for stresses and strains in the block, traction on the fault and slip on the fault.

Domain output:
So far, I managed to get my initial dynamic stresses look like the last quasistatic output, however, strains are zero in the domain.
I’m doing smth like this:
db_initial_stress = spatialdata.spatialdb.SimpleDB
db_initial_stress.label = Initial stress
db_initial_stress.iohandler.filename = spatialdb/initial_stress_from_prestep.spatialdb
#db_initial_stress.query_type = linear

db_initial_strain = spatialdata.spatialdb.SimpleDB
db_initial_strain.label = Initial strain
db_initial_strain.iohandler.filename = spatialdb/initial_stress_from_prestep.spatialdb
#db_initial_strain.query_type = linear

This file contains both stresses and strains.
Should value-names for strains be total-strain-xx or strain-xx? (Doesn’t work either way though).
Also, should I use initial zz stress? When I try to do this (with stress-zz-initial in the file):
#db_initial_state = spatialdata.spatialdb.SimpleDB
#db_initial_state.label = Initial state
#db_initial_state.iohandler.filename = spatialdb/initial_stress_from_prestep.spatialdb
#db_initial_state.query_type = linear
Pylith outputs an error:
mpinemesis: …/…/…/pylith-2.2.1/libsrc/pylith/materials/ virtual void pylith::materials::Material::initialize(const pylith::topology::Mesh&, pylith::feassemble::Quadrature*): Assertion `numDBStateVars > 0’ failed.
[0]0:Return code = 0, signaled with Aborted

Fault output:
This is what I was doing for tractions:

Initial tractions

traction_perturbation = pylith.faults.TractPerturbation

db_initial = spatialdata.spatialdb.SimpleDB
db_initial.label = Initial fault tractions
db_initial.iohandler.filename = spatialdb/initial_fault_tractions.spatialdb
db_initial.query_type = linear

Didn’t work, tractions on the fault are zero.
How should I prescribe initial tractions?
Is there a way to prescribe initial slip for FaultCohesiveDyn or is it only possible for the kinematic one?

I just need to have the initial state on my model the same as my quasistatic output. Could you please tell me what I’m doing wrong or if there’s any simpler way to achieve that (like with an elastic prestep)?
Please let me know if you need a sketch or any of the codes.

Thank you in advance!

It isn’t clear what you mean by a planer fault with three sections that have different fault constitutive models. PyLith doesn’t support a single fault with multiple fault constitutive models. If you break up a single fault into three fault objects, each with its own fault constitutive model, the combined fault will have pinned (zero slip) at the intersections between the segments. The recommended approach is to use a single fault and fault constitutive model with spatially variable parameters.

If you are using a linearly elastic bulk constitutive model, then you do not need initial stresses and strains. All you need is the initial tractions on the fault. See examples/2d/subduction/step04 for an example of specifying initial tractions.

Dear Brad,

Sorry, that was indeed unclear. It’s the same fault with spatial parameter variation: e.g. slip-weakening friction everywhere with mu_s =mu_d for the middle section.

I’m doing exactly as in examples/step04.
I can get initial tractions on the fault (thought at time step 0 everything is still zero and they appear only at time step 1. Is this how it’s meant to be?). However, if I don’t prescribe any initial stresses within the block, stresses in the block are zero initially and are close to zero during the following time steps. Same with strains. And the slip on the fault is also zero initially.
If I prescribe initial stress field (uniform or generated from my quasistatic output) along with tractions on the fault, at time step 0 I only see stresses in the block and at time step 1 stresses on the fault are summed up with stresses in the block and I get twice larger values than I need.
I guess I can subtract block stresses from my fault stresses and only prescribe the perturbation. This will give somewhat realistic stresses (both in the block and on the fault) at time step 1, but that still doesn’t change the initial strain field or displacements.
Is there a way I can add my displacements (both on the fault and in the block)?


I can get initial tractions on the fault (thought at time step 0 everything is still zero and they appear only at time step 1. Is this how it’s meant to be?). However, if I don’t prescribe any initial stresses within the block, stresses in the block are zero initially and are close to zero during the following time steps. Same with strains. And the slip on the fault is also zero initially.

I think this is what you want. I think it is likely (depending on your setup) that even though the output of the fault tractions is zero at time step 0, they are not actually zero; this can happen as a result of how we do output in v2.2.x.

In general, you either impose initial fault tractions or initial stresses not both. If you do impose initial stresses in the bulk material, then when they are resolved onto the fault, you will get the corresponding tractions. Thus, when you impose both, you get double the fault tractions.

When you impose the initial fault tractions, the deformation will be relative to this traction state, so you should see zero displacements if there is no fault slip and no external forcing.

Dear Brad,

Thank you!
I guess, if I prescribe my block stresses and they are projected on the fault, this will be the closest to what I want.
Is the calculation still physical with zero initial displacements? Of course, I can always add them during postprocessing, but I’m worried if the response of the model would be the same with no initial deformation and just the stress field.

As long as you are using a linearly elastic bulk rheology all of the deformation will be the perturbation from the deformation associated with that initial stress state. As you mentioned, you can add the perturbation to the original displacement field to get a “total” displacement.

Thank you, Brad.
This was very helpful, as always!

Dear Brad,

I tried what you suggested: prescribed initial block stresses and they got resolved onto the fault.

However, due to all my displacements being zero, the slip is not propagating in the direction I want, e.g. with slip weakening case (middle weaker section - both coefficients 0.6, and side sections mud = 0.6 and mus = 0.7) rupture doesn’t propagate outwards into stronger sections, but propagates into the middle weaker section to generate a parabolic slip profile and catch up with the stress distribution first, which makes sense, given zero initial displacements but which is not at all what I want.

Do you have any suggestions on how to resolve the issue?
Do you think prescribing my displacement profile instead of tractions will give a better result?


Dear Brad,

I’m still very interested making this happen, so if you have any advice, please let me know.

The issue of which faults rupture and under what conditions is really a research question and is beyond the scope of the support we can provide to users. We can provide support in terms of how to use the code and making it is working properly. You may want to collaborate with an experienced modeler to address these research issues.