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 laboratoryderived constitutive laws such as the rateandstate 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 (antidip) direction: positive (negative)
Normal stress = negative in compression
Cohesion = negative (acts on shear stress components, 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 righthanded coordinate system is required. Please, make sure to follow this convention during model/mesh generation!
Arbitrary fault orientation possible, except a fault in the xyplane, 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) alongstrike vector
d = (sy*nz, sx*nz, synxnysx).T /  d  alongdip 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 downdipdistance)
shear stress = 0.55 * (normal stress) = + 4057.9 Pa/m * (meter downdipdistance)
Assuming the fault plane in xzplane, dipping 60° in ydirection we can use the script ‘/preprocessing/science/stressrotation3D’ to rotate the stresses 30° clockwise at the xaxis (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 downdipdistance)
s_zz =  5358.74 Pa/m * (meter downdipdistance)
s_yz = + 5223.72 Pa/m * (meter downdipdistance)
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 strikeslip fault aligned with the x or yaxis 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 rightlateral strikeslip fault in xzplane 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 zaxis 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 leftlateral 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 reactivation 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 (hrefinement) 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 hrefinement 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 sliprate from the difference between the velocity on both sides of the fault, rather than evaluating the sliprate 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 slipweakening friction (FL=6
, FL=16
)
The linear slipweakening friction is widely used for dynamic rupture simulations.
The fault strength is determined by
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_d(x)\) 
dynamic friction coefficient 

\(d_c(x)\) 
slipweakening critical distance 

\(C(x)\) 
cohesion 

\(T(x)\) 
forced rupture time 

\(v_0\) 
threshold velocity 

Friction law 16
implements linear slipweakening 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 PrakashClifton regularization for bimaterial faults.
For friction law 16
, we resample the slip rate in every step to suppress spurious oscillations.
In the case of PrakashClifton 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 rateandstate friction with
Now the state variable stores the accumulated slip.
Rateandstate friction
Rateandstate 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 rateandstate friction laws, which differ by the set of ordinary differential equations describing the evolution of the state variable.
The type of rateandstate friction law is set by the FL variable in the DynamicRupture namelist (parameters.par):
Friction law 3
implements the ageing law, friction law 4
implements the slip law, and friction law 103
implements a slip law with strong rateweakening.
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 ageing law
and for the rate and state friction with strong velocity weakening
are available at the given links.
All rateandstate friction laws are described by the following system of differential algebraic equations, which depend on the state variable \(\psi\) and the slip velocity \(V\).
Ageing law (FL=3
)
Reference benchmarks: TVP101 and TPV102
Friction parameters:
symbol 
quantity 
seisSol name 

\(a(x)\) 
frictional evolution coefficient 

\(b\) 
frictional state coefficient 

\(L(x)\) 
characteristic slip scale 

\(V_0\) 
reference slip velocity 

\(f_0\) 
reference friction coefficient 

Slip law (FL=4
)
The slip law has the same parameters as the Ageing Law.
Strong velocity weakening (FL=103
)
Reference TPV103 and TPV104
In addition to the ageing and the slip Law, strong velocity weakening requires two more parameters:
symbol 
quantity 
seisSol name 

\(V_w(x)\) 
weakening slip velocity 

\(\mu_w\) 
weakening friction coefficient 

with
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 rateandstate 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.0e6 ! 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 spacedependent and therefore have to be specified in the dynamic rupture yaml file:
!ConstantMap
map:
tp_hydraulicDiffusivity: 1e4 # Hydraulic diffusivity [m^2/s]
tp_halfWidthShearZone: 0.01 # Half width of shearing zone [m]
TP generates 2 additional onfault outputs: Pore pressure and temperature (see fault output).