Dynamic rupture

Dynamic rupture#

SeisSol is verified for a wide range of dynamic rupture problems (Pelties et al. 2014): it is possible to use branched faults, dipping fault geometries and laboratory-derived constitutive laws such as the rate-and-state friction law.

Sign conventions#

The following definitions regarding fault and stress orientation are conventions. The initial stresses need to be given in the global coordinate system. The Matlab script ‘/preprocessing/science/stressrotation3D’ is available for correct rotation of initial stress values on arbitrary orientated faults.

Definitions#

  • Traction = (normal vector)* (stress tensor)

    • Shear traction in strike direction right (left) lateral: positive (negative)

    • Shear traction along dip (anti-dip) direction: positive (negative)

  • Normal stress = negative in compression

  • Cohesion = negative (acts on normal stress, encapsulates the effect of pore pressurization)

  • Reference point = defines dip direction, starting point of the normal vector n pointing towards the fault, at +y

Note

SCEC benchmark descriptions define normal stress as positive in compression, opposite to SeisSol 3D convention.

Fault geometry#

The right-handed coordinate system is required. Please, make sure to follow this convention during model/mesh generation!

  • Arbitrary fault orientation possible, except a fault in the xy-plane, i.e. parallel to the free surface

  • Fault Output is in fault coordinate system

    • n = ( nx, ny, nz ).T normal vector, from + y → - y

    • s = ( ny, - nx, 0).T * sqrt ( nx^2 + ny^2) along-strike vector

    • d = (-sy*nz, sx*nz, synx-nysx).T / || d || along-dip vector, pointing in -z direction (downwards)

Example 1: Initial stresses for a 60° dipping fault#

Following the SCEC benchmark description of TPV10 we need to transfer the following initial stresses given in fault coordinates in global coordinates for SeisSol:

INPUT (from tpv10)

normal stress = - 7378 Pa/m * (meter downdip-distance)
shear stress = 0.55 * (normal stress) = + 4057.9 Pa/m * (meter downdip-distance)

Assuming the fault plane in xz-plane, dipping 60° in -y-direction we can use the script ‘/preprocessing/science/stressrotation3D’ to rotate the stresses 30° clockwise at the x-axis (considering we look towards +x) to get the full stress tensor in global coordinates:

OUTPUT (to be used in SeisSol’s parameter file)

s_yy = -  2019.26 Pa/m * (meter downdip-distance)
s_zz = -  5358.74 Pa/m * (meter downdip-distance)
s_yz = + 5223.72 Pa/m * (meter downdip-distance)
s_xx = s_xy = s_xz = 0.0 Pa

Example 2: Initial stresses for a branched fault#

Normal and shear stresses in fault coordinate system for strike-slip fault aligned with the x- or y-axis are already in the global coordinates. But if we have additional branches the initial stresses need to be rotated again. Following for example the SCEC benchmark description for TPV14, the stresses for the right-lateral strike-slip fault in xz-plane are:

s_yy = - 120 MPa
s_xy = + 70 MPa
s_xy (nucleation patch) = + 81.6 MPa
s_xx = s_zz = s_xz = s_yz = 0.0 Pa

For the 30° branch in -y direction we rotate the given stresses (in that case, they are the same as on the main fault) by 330° counterclockwise (or 30° clockwise) around the z-axis using the script ‘/preprocessing/science/stressrotation3D’.

OUTPUT (to be used in SeisSol’s parameter file)

s_xx = + 30.62 MPa
s_yy = - 150.62 MPa
s_xy =  - 16.96 MPa
s_zz = s_xz = s_yz = 0.0 Pa

In case of a left-lateral strike slip fault (for example TPV15) the stresses on the main fault are:

s_yy = - 120 MPa
s_xy =  - 70 MPa
s_xy (nucleation patch) = - 81.6 MPa
s_xx = s_zz = s_xz = s_yz = 0.0 Pa

Projecting the state variable increment#

A sudden decrease of slip rate to very small values, as for example occurring at a rupture front about to be arrested, may cause numerical issues and pollute the solution - in the worst case even leading to re-activation of rupture. Such numerical issues are easy to diagnose in the fault output visualization: a checkerboard pattern with unphysical values for the slip rate and slip in the affected region will be visible. This is a sign of mesh resolution (h-refinement) being locally not sufficient. SeisSol mitigates such artifacts by projecting the state variable (e.g. cumulative slip for linear slip weakening) increment on the 2D fault interface basis function after evaluation of the friction law. This implementation lowers the required h-refinement across the fault in comparison to earlier implementations (cf. Pelties et al. , JGR 2012; GMD 2014)

Visualisation: SlipRateOutputType (default =1)#

By default, the fault output will be showing regions affected by numerical problems. However, the user may choose to smooth out such artifacts for visualization purposes. Switching SlipRateOutputType in the DynamicRupture namelist from the default value to 0, will evaluate the slip-rate from the difference between the velocity on both sides of the fault, rather than evaluating the slip-rate from the fault tractions and the failure criterion directly. Note that this fix only applies to the output, but does not suppress numerical problems in the simulation. Also note that that SlipRateOutputType=0 is slightly less accurate than the default SlipRateOutputType=1 without numerical problems.

Friction laws#

Linear slip-weakening friction (FL=6, FL=16)#

The linear slip-weakening friction is widely used for dynamic rupture simulations.

The fault strength is determined by

\[\tau = -C - \min\left(0, \sigma_n\right) \left( \mu_s - \frac{\mu_s - \mu_d}{d_c} \min\left(S, d_c\right)\right),\]

where \(S(t) = \int_0^t |V(s)| ds\) is the accumulated fault slip, and the other variables are parameters of the friction, detailed below.

Friction parameters:

symbol

quantity

SeisSol name

\(\mu_s(x)\)

static friction coefficient

mu_s

\(\mu_d(x)\)

dynamic friction coefficient

mu_d

\(d_c(x)\)

slip-weakening critical distance

d_c

\(C(x)\)

cohesion

cohesion

\(T(x)\)

forced rupture time

forced_rupture_time

\(t_0\)

characteristic nucleation time

t_0

\(v_0\)

threshold velocity

lsw_healingThreshold

Friction law 16 implements linear slip-weakening with a forced rupture time. If you are only interested in linear slip weakening friction without forced rupture time, do not supply the parameter forced_rupture_time in the fault yaml file.

Friction law 6 uses Prakash-Clifton regularization for bimaterial faults. For friction law 16, we resample the slip rate in every step to suppress spurious oscillations. In the case of Prakash-Clifton regularization, we do not resample the slip rate. If the slip rate \(V\) drops below the threshold velocity \(v_0\), we reset the friction parameter \(\mu = \mu_s\). Also, we reset the state variable \(S = 0\). The threshold \(v_0\) is set to \(-1.0\) by default, such that healing is disabled.

Examples of input files for the friction laws 6 and 16 are availbable in the cookbook.

Linear slip weakening can be seen as a special case of rate-and-state friction with

\[\begin{split}\begin{aligned} f(V, \psi) &= C - \left( \mu_s - \frac{\mu_s - \mu_d}{d_c}\right) \min\left(\psi, d_c\right), \\ g(V, \psi) &= V. \end{aligned}\end{split}\]

Now the state variable stores the accumulated slip.

TP proxy slip-weakening friction (FL=1058)#

The TP proxy slip-weakening friction has been proposed by Herrera et al. (2024), GJI, to approximate thermal pressurization in a computationally efficient way. The fault strength is determined by

\[\tau = -C - \min\left(0, \sigma_n\right) \left( \mu_d + \frac{{(\mu_s - \mu_d)}}{{\left(1 + \frac{S}{d_c}\right)^{\alpha}}} \right),\]

All variables are the same as defined in previous section for FL=16. The friction law also supports forced rupture time. You can modify the default value 1/3 of \(\alpha\), by adjusting the TpProxyExponent parameter in the main parameter file (namelist: DynamicRupture).

Rate-and-state friction#

Rate-and-state friction laws allow modeling the frictional shear strength variations as a function of slip rate and of the evolving properties of the contact population (Dieterich, 1979, 1981; Ruina, 1983). In SeisSol, we currently support 3 types of rate-and-state friction laws, which differ by the set of ordinary differential equations describing the evolution of the state variable. The type of rate-and-state friction law is set by the FL variable in the DynamicRupture namelist (parameters.par): Friction law 3 implements the aging law, friction law 4 implements the slip law, and friction law 103 implements a slip law with strong rate-weakening. More details about these friction laws can be found in the SCEC benchmarks descriptions (TPV101 to 105) or in Pelties et al. (2013, GMD).

Some parameters are considered homogeneous across the fault and defined in the main parameter file. Others can spatially vary (rs_a, rs_sl0 for FL=3,4 and 103 and rs_srW for FL=103) and are defined in the fault yaml file. Examples of input files for the aging law and for the rate and state friction with strong velocity weakening are available at the given links.

All rate-and-state friction laws are described by the following system of differential algebraic equations, which depend on the state variable \(\psi\) and the slip velocity \(V\).

\[\begin{split}\begin{aligned} \tau &= \sigma_n f(V,\psi) \\ \frac{\partial\psi}{\partial t} &= g(V,\psi) \end{aligned}\end{split}\]

Aging law (FL=3)#

Reference benchmarks: TVP101 and TPV102

Friction parameters:

symbol

quantity

SeisSol name

\(a(x)\)

frictional evolution coefficient

rs_a

\(b(x)\)

frictional state coefficient

rs_b

\(L(x)\)

characteristic slip scale

rs_sl0

\(V_0\)

reference slip velocity

rs_sr0

\(f_0(x)\)

reference friction coefficient

rs_f0

\(V_{ini1}\)

initial along-strike slip velocity

rs_inisliprate1

\(V_{ini2}\)

initial along-dip slip veloctiy

rs_inisliprate2

Here, \(b(x)\) and \(f_0(x)\) can either be varied or given in the parameter file as constants. Once they are given in the fault yaml file, they are assumed to be varied and their parameter file values are ignored. For the parameter file and the fault yaml file, the names are the same.

\[\begin{split}\begin{aligned} f(V, \psi) &= a \sinh^{-1}\left[\frac{V}{2V_0} \exp\left( \frac{f_0 + b \ln(V_0 \psi / L)}{a}\right) \right] \\ g(V, \psi) &= 1 - \frac{V \psi}{L} \end{aligned}\end{split}\]

Slip law (FL=4)#

The slip law has the same parameters as the Aging Law.

\[\begin{split}\begin{aligned} f(V, \psi) &= a \sinh^{-1}\left[\frac{V}{2V_0} \exp\left( \frac{f_0 + b \ln(V_0 \psi / L)}{a}\right) \right] \\ g(V, \psi) &= -V\frac{\psi}{L}\ln \left(V \frac{\psi}{L} \right) \end{aligned}\end{split}\]

Severe velocity weakening (FL=7)#

No reference benchmark.

The Severe Velocity Weakening Law has the same parameters as the Aging Law does.

Strong velocity weakening (FL=103)#

Reference TPV103 and TPV104

In addition to the Aging and the Slip Law, strong velocity weakening requires two more parameters:

symbol

quantity

seisSol name

\(V_w(x)\)

weakening slip velocity

rs_srW

\(\mu_w(x)\)

weakening friction coefficient

rs_muW

Here, \(\mu_w(x)\) behaves similarly to \(b(x)\) and \(f_0(x)\); i.e. it can be specified in either the parameter file as a global constant or in the fault yaml file to vary; the same rules as above apply here as well.

\[\begin{split}\begin{aligned} f(V, \psi) &= a \sinh^{-1}\left[\frac{V}{2V_0} \exp\left(\frac{\psi}{a}\right) \right] \\ g(V, \psi) &= - \frac{V}{L} \left(\psi - a \ln\left[ \frac{2V_0}{V} \sinh\left( \frac{\mu_{ss}(V)}{a} \right) \right] \right) \end{aligned}\end{split}\]

with

\[\begin{aligned} \mu_{ss}(V) = \mu_w + \frac{f_0 - (b-a) \ln\left( \frac{V}{V_0} \right) - \mu_W}{\left( 1 + \left[ \frac{V}{V_W}\right]^8\right)^{1/8}} \end{aligned}.\]

Note that from the merge of pull request #306 of March 17th, 2021 to the merge of pull request #752 of December 22nd, 2022, the state variable was enforced positive in this friction law. This enforcement aimed at avoiding the state variable getting negative because of Gibbs effects when projecting the state increment onto the modal basis functions (resampling matrix). Since then, we realized that the state variable can get negative due to other factors, and, therefore, reverted this change.

Thermal Pressurization#

Seissol can account for thermal pressurization (TP) of pore fluids. As deformation occurs within the fault gauge, frictional heating increases the temperature of the rock matrix and pore fluids. The pore fluids then pressurize, which weakens the fault. The evolution of the pore fluid pressure and temperature is governed by the diffusion of heat and fluid. TP can be activated using thermalPress in the DynamicRupture namelist. TP can be enabled with all rate-and-state friction laws (FL=3,4 and 103). The TP parameters for which no spatial dependence has been implemented are defined directly in the DynamicRupture namelist:

&DynamicRupture
thermalPress = 1                     ! Thermal pressurization 0: inactive; 1: active
tp_iniTemp = 483.15                  ! Initial temperature [K]
tp_iniPressure = -80.0e6             ! Initial pore pressure; have to be added to normal stress in your initial stress yaml file [Pa]
tp_thermalDiffusivity = 1.0e-6       ! Thermal diffusivity [m^2/s]
tp_heatCapacity = 2.7e6              ! Specific heat [Pa/K]
tp_undrainedTPResponse = 0.1e6       ! Pore pressure change per unit temperature [Pa/K]

Two additional thermal pressurization parameters are space-dependent and therefore have to be specified in the dynamic rupture yaml file:

!ConstantMap
map:
  tp_hydraulicDiffusivity: 1e-4   # Hydraulic diffusivity [m^2/s]
  tp_halfWidthShearZone: 0.01     # Half width of shearing zone [m]

TP generates 2 additional on-fault outputs: Pore pressure and temperature (see fault output).

Nucleation#

SeisSol provides several options to enforce artificial nucleation within a limited temporal and spatial extent. The additional nucleation stress applied to the fault plane needs to be specified using the parameters Tnuc_s (along strike), Tnuc_d (along dip), and Tnuc_n (normal stress). Alternatively, if the stress state is described by the full stress tensor, use parameters nuc_xx, nuc_yy, nuc_zz, nuc_yz, nuc_xz, and nuc_xy. The parameter s_0 defines the start time of artificial nucleation, while t_0 sets the characteristic nucleation time over which the additional stress is applied. By default, s_0 = 0 and t_0 = 0.

Multiple nucleations#

SeisSol supports multiple episodes of artificial nucleation by extending the nucleation parameters. For each episode, parameters can be specified with numbered suffixes:

  • Nucleation start times: s_0, s2_0, s3_0, …

  • Characteristic nucleation times: t_0, t2_0, t3_0, …

  • Projected nucleation stress components: Tnuc_s, Tnuc2_s, Tnuc3_s, …

  • Nucleation stress tensor components: nuc_xy, nuc2_xy, nuc3_xy, …

Forced rupture time#

If forced_rupture_time is specified in the fault yaml file, the friction coefficient is artificially reduced from its static value to its dynamic value at the given time (e.g., Gabriel et al., 2012). This method enables a smooth nucleation process without affecting the overall stress drop. In this case, the characteristic nucleation time t_0 defines the duration of the linear transition from static to dynamic friction, and s_0 is ignored.

Contents