Skip to content

Latest commit

 

History

History
1429 lines (1140 loc) · 58.3 KB

index.rst

File metadata and controls

1429 lines (1140 loc) · 58.3 KB

Introduction

The software package SeisSol (http://www.seissol.org/) allows for realistic simulations of the three-dimensional seismic wave field propagating in complex Earth structures generated by a finite dynamic earthquake source governed by a constitutive law that describes the relationship between fault stress and slip along a geometrically complex fault. SeisSol is a high-order accurate Discontinuous Galerkin Finite Element solver, based on the ADER-DG method presented in:raw-latex:citep{kaser2006}, enabling precise modelling of on-fault frictional failure coupled to seismic waves travelling over large distances in terms of propagated wavelengths with minimal dispersion errors, whereas it is intrinsically dissipative and removes frequencies unresolved by the mesh without affecting longer and physically meaningful wavelengths.

The software has recently proven to be highly scalable on current and future HPC infrastructure. It reached multi-petaflop/s performance on some of the largest supercomputers worldwide, and was a Gordon Bell prize finalist application 2014 (Breuer et al., 2014; Heinecke et al., 2014) in a pioneering simulation of the 1992 M7.2 Landers earthquake. High detail rupture evolution and synthetic ground shaking in the engineering frequency band (0-10 Hz) were modeled on a non-planar earthquake fault structure. In early 2017, SeisSol performed the longest and largest dynamic rupture scenario to date, enabled by local time stepping :raw-latex:`\citep{uphoff2017}`, resolving the 2004 Sumatra-Andaman earthquake including complex splay fault geometries. The paper won the prestigous “Best Paper Award” of the International Supercomputing Conference (SC17). Recently, large scale dynamic rupture forward scenarios of the 1994 Northridge earthquake and the Husavik-­Flatey­ fault system in Northern Iceland have been developed shedding light on pressing questions of fault mechanics (Heinecke et al., 2014; Gabriel et al., 2014; Wollherr et al., 2015; Ulrich et al., 2016, Weingärtner et al., 2016). SeisSol results imply that acknowledging geometrical complexity, realistic fault properties and velocity models affect not only earthquake source dynamics but the synthetic ground shaking crucially. The software package is available to the community as an open source distribution (www.github.com/SeisSol/SeisSol).

This documentation is a collection of useful dynamic simulation examples to help users build models from the scratch with little efforts. Each example is demonstrated carefully with geometry building, parameter setup and result visualization. Users are suggested to repeat each example in order to get a comprehensive idea of how to set up dynamic simulation models with SeisSol. The updates of source software and news of SeisSol should be posted in SeisSol or .

Examples

Overview

SeisSol is a part of SCEC dynamic code validation project :raw-latex:`\citep{harris2018}`. Here we show several SCEC benchmarks for beginners to quickly catch up with SeisSol workflow. Each benchmark example is composed of a short problem description, a section of geometry, initial setups (stress, nucleation, friction, etc.), and simulation results.

| p2.4cm | p2.6cm | p2.4cm| p7.5cm | No. & Fault type & Difficulty & Description
TPV5 & strike-slip & beginner& slip-weakening and heterogeneous initial stress conditions
TPV6 & strike-slip & beginner& bi-material fault and, slip-weakening and heterogeneous initial stress conditions
TPV12 & normal fault & beginner & linear elastic and initial stress conditions are dependent on depth
TPV13 & normal fault & beginner & non-associative Drucker-Prager plastic with yielding in shear ad initial stress conditions are dependent on depth
TPV16 & strike-slip & intermediate & randomly-generated heterogeneous initial stress conditions
TPV24 & branching strike-slip & intermediate & a rightward branch forming a 30 degree angle. There are linear elastic material properties
TPV29 & strike-slip & difficult & stochastic roughness. Linear elastic material properties in a homogeneous half-space.
TPV104 & strike-slip & difficult & Rate-state friction, using a slip law with strong rate-weakening.
Point Source & strike-slip & intermediate & benchmark of SISMOWINE WP2_LOH1.
Kinematic & reverse fault & intermediate & Kinematic source of 1994 Mw6.7 Northridge earthquake.

Prerequisites

Before you begin any of the examples, you will need to install latest SeisSol from . The instruction can be found at . All geometry and tetrahedral meshing are generated using free software Gmsh () is needed. If you do not wish to create your own mesh at this time, the meshes are also provided as part of the example. The ParaView visualization package may be used to view simulation results. You may use other visualization software, but some adaption from what is described here will be necessary. Furthermore, you can complete a subset of the example using files provided (as described below), skipping the steps for which you do not have the proper software packages installed

Input file resources

The files needed to work through the examples are found in . Users may download freely from the repository.

SCEC TPV5

TPV5 is the first SCEC benchmark. It has spontaneous rupture on a vertical strike-slip fault in a homogeneous halfspace. There are slightly heterogeneous initial stress conditions. The earthquake rupture is artificially nucleated in a square zone at the center of the fault surface. The rupture then spontaneously propagates over the rest of the fault surface. As it propagates away from the nucleation zone, it encounters two square patches with initial stress conditions that are different from the rest of the fault surface.

Geometry

The fault within the three-dimensional medium is a vertical right-lateral strike-slip planar fault that resides in a half-space. The fault reaches the Earth’s surface. The rupture is allowed within a rectangular area that is 30000 m long \times 15000 m deep. The bottom boundary of and the right and left ends of the allowed 30000 m \times 15000 m rupture area are defined by a strength barrier. The nucleation point is centered both along-dip and along-strike of the 30000m \times 15000m rupture area, on the fault plane, at 15000m along-strike and 7500m depth.

The mesh is generated in GMSH. All the files that are needed for the simulation are provided in . The tpv5.geo file contains the geometry for the fault in a cubit region.

Then the .geo file can be meshed by using:
gmsh tpv5.geo -3 -optimize -o tpv5.msh
Then convert the .msh file to 3D Gambit neutral file
gmsh2gambit -i tpv5.msh -o tpv5.neu

The toolbox of gmsh2gambit can be found in .

The 3D Gambit file can be convert to PUML format for LTS in latest version of SeisSol by:
pumgen tpv5.neu tpv5

The compilation and usage of PUMGen can be found in

Parameters

Nucleation

occurs because the initial shear stress in a 3000 m \times 3000 m square nucleation patch is set to be higher than the initial static yield stress in that patch. Failure occurs everywhere on the fault plane, including in the nucleation patch, following a linear slip-weakening fracture criterion.

TPV5 uses a linear-slip weakening friction everywhere on the fault. There are ten parameters associated with the friction constitutive law and fault properties in the parameters.par. It can be found at .

Four friction constitutive parameters are: mu_s, mu_d, d_c and cohesion. Six stress parameters are: s_xx, s_yy, s_zz, s_xy, s_xz, and s_yz. All the parameters are homogeneous on the fault except for the nucleation patch in the center of the fault, where s_xy is larger compared with that elsewhere. The parameters in TPV5 are listed in Table [table:tpv5].

Parameter Description Value Unit
mu_s static friction coefficient 0.677 dimensionless
mu_d dynamic friction coefficient 0.525 dimensionless
d_c critical distance 0.40 m
cohesion friction cohesion 0.0 MPa
s_yy stress 120 MPa
s_xy stress 70 MPa
s_xx,s_zz,s_yz,s_sxz stress 0 MPa
s_xy shear stress 81.6 MPa

Table: Table of LSR parameters on the fault.

Results

All examples here can be illustrated in Paraview (Detailed instruction can be found at ). The output folder contains a series of files for fault dynamic rupture (netcdf), wave filed (netcdf), receiver (.dat) and off-fault receivers (.dat). The fault dynamic rupture and wave filed files can be loaded in Paraview directly. For example, open Paraview and then go through File >> import >>prefix-fault.xdmf.

In the wave filed output file (prefix.xdmf, prefix_vertex.h5 and prefix_cell.hf), the variables are shown in Table [table:wavefield]

Index Parameter Description
1 U displacement in x-axis
2 V displacement in y-axis
3 W displacement in z-axis
4 u particular velocity in x-axis
5 v particular velocity in y-axis
6 w particular velocity in z-axis

Table: Table of wave field output in SeisSol. Index denotes the position used in iOutputMask in SeisSol parameter file.

In the fault dynamics output file (prefix-fault.xdmf, prefix-fault_vertex,h5 and prefix-fault_cell,h5), the variables are shown in Table [table:faultout]

Index Parameter Description
1 SRs and SRd slip rates in strike and dip direction
2 T_s, T_d, P_n transient shear stress in strike and dip direction, transient normal stress
3 U_n normal velocity (note that there is no fault opening in SeisSol)
4 Mud, StV current friction and state variable in case of RS friction
5 Ts0,Td0,Pn0 total stress, including initial stress
6 Sls and Sld slip in strike and dip direction
7 Vr rupture velocity, computed from the spatial derivatives of the rupture time
8 ASl absolute slip
9 PSR peak slip rate
10 RT rupture time
11 DS only with LSW, time at which ASl > d_c

Table: Table of fault dynamic output in SeisSol. Index denotes the position used in iOutputMask in SeisSol parameter file.

SCEC TPV6

TPV6 is intended to reside in the “well-posed” regime for bimaterial problems and so uses a very high shear modulus (density*vs*vs) contrast. Material properties are homogeneous within each side of the fault, but change when one traverses to the other side of the fault. This is the bi-material problem.

| p3.5cm | p7cm | p2cm| p2cm | Parameter & Description & Value & Unit
Vp & P velocity & 3750 &m/s
Vs & S velocity & 2165 & m/s
\rho & density & 225 & kg/m3
Vp & P velocity& 6000 & m/s
Vs & S velocity& 3464 & m/s
\rho & density & 2670 &kg/m3

Geometry

TPV6 uses the same geometry as TPV5. The fault within the three-dimensional medium is a vertical right-lateral strike-slip planar fault that resides in a half-space. The fault reaches the Earth’s surface. The rupture is allowed within a rectangular area that is 30000 m long \times 15000 m deep. The bottom boundary of and the right and left ends of the allowed 30000 m \times 15000 m rupture area are defined by a strength barrier. The nucleation point is centered both along-dip and along-strike of the 30000m \times 15000m rupture area, on the fault plane, at 15000m along-strike and 7500m depth.

Parameters

TPV6 uses a similar parameter setup as TPV5 except for the bulk parameters.

| p3.5cm | p7cm | p2cm| p2cm | Parameter & Description & Value & Unit
\rho & density & 2225 & kg/m3
\lambda & Lame parameter& 10.4 & GPa
\mu & Lame parameter & 10.4 & GPa
\rho & density & 2670 &kg/m3
\lambda & Lame parameter& 32 & GPa
\mu & Lame parameter& 32& GPa

Results

Figure [fig:tpv6-4s] and [fig:tpv6-7s] show the fault slip rate at 4 s and 7 s, respectively. The slip front is asymmetric when compared with TPV5 (Figure [fig:tpv5-4s]). Figure [fig:tpv6_velocity] shows velocity recorded at two off-fault receivers. The wave picks arrives at the far-side receiver lower than those at the near-side receiver.

SCEC TPV12

TPV12 and 13 are recommended by SCEC for elastic/plastic wave propagation code validation. TPV 12 describes spontaneous rupture on a 60-degree dipping normal fault in a homogeneous half-space. Material properties are linear elastic. Initial stress conditions are dependent on depth. Strongly super-shear rupture conditions.

Geometry

The model volume is a half-space. The fault is a 60-degree dipping, planar, normal fault. The fault reaches the Earth’s surface. Rupture is allowed within a rectangular area measuring 30000 m along-strike and 15000 m down-dip.

Note that 15000 m down-dip corresponds to a depth of 12990.38 m. A node which lies exactly on the border of the 30000 m \times 15000 m rectangle is considered to be inside the rectangle, and so should be permitted to rupture.

The portions of the fault below, to the left of, and to the right of the 30000 m \times 15000 m rectangle are a strength barrier, within which the fault is not allowed to rupture.

The nucleation zone is a square measuring 3000 m × 3000 m. The center of the square is located 12000 m down-dip (at a depth of 10392.30 m), and is centered along-strike.

The geometry is generated with GMSH. All the files that are needed for the simulation are provided in

The geometry and mesh generation process is similar as TPV5. The planar-fault geometry is build with Gmsh (Figure [fig:tpv12geo]). All the files that are needed for the simulation are provided in .

Nucleation

In previous benchmarks, nucleation was achieved by imposing a higher initial shear stress within a nucleation zone. In TPV12 and TPV13, nucleation is achieved by selecting a lower static coefficient of friction within a nucleation zone, so that the initial shear stress (which is implied by the initial stress tensor) is greater than the yield stress.

Outside the 30000 m * 15000 m rectangular rupture area there is a strength barrier, where nodes are not allowed to slip. Some codes implement the strength barrier by setting the static coefficient of friction and frictional cohesion to very large values. Other codes implement the strength barrier in other ways.

Parameters

LSR parameters

TPV12 uses a linear slip weakening law on the fault with different parameters inside and outside the nucleation zone. The parameters are listed in Table [table:tpv12lsr].

| p3.5cm | p7cm | p2cm| p2cm | Parameter & Description & Value & Unit
mu_s & static friction coefficient & 0.54 &
mu_d & dynamic friction coefficient & 0.10 &
d_c & critical distance & 0.50 & m
cohesion & shear stress cohesion & -200 000 & Pa
mu_s & static friction coefficient & 0.70 &
mu_d & dynamic friction coefficient & 0.10 &
d_c & critical distance & 0.50 & m
cohesion & shear stress cohesion & -200 000 & Pa

Initial stress

The initial stress on the fault is depth-dependent in TPV12/13. In the shallower portion above 11951.15 m, the stress field is optimal orientated while the other is isotropic.

| p7cm | p7cm | Parameter & Value
\sigma_1 & 26460 Pa/m * H
\sigma_3 & 15624.3 Pa/m * H
\sigma_2 & (\sigma_1+\sigma_3)/2
P_f & 1000 kg/m^3 *9.8 m/s^2 *H
\sigma_1,\sigma_2,\sigma_3 & 2700 kg/m^3 *9.8 m/s^2 *H

Results

SeisSol output xdmf file that can be loaded in Paraview directly. The wave field and fault output files have the same format as in TPV5.

SCEC TPV13

TPV13 is similar with TPV12 except for that material properties are non-associative Drucker-Prager plastic. TPV13 is using non-associative Drucker-Prager plasticity with yielding in shear. The material is characterized by six constitutive parameters:
*Bulk friction = 0.85 Fluid pressure = 1000 kg/m3 * Based on the SCEC benchmark, the stress components \sigma_{ij} are defined as :
*Mean stress: \sigma_m = (\sigma_{11}+\sigma_{22}+\sigma_{33})/3 Stress deviator: s_{ij} = \sigma_{ij} - \sigma_m \delta_{ij} Second invariant of the stress deviator: J_2(\sigma) = 1/2 *\sum_{ij} s_{ij} s_{ji} Drucker-Prager yield stress: Y(\sigma) =\max(0,c\cos \phi - (\sigma_m +P_f)\sin \phi) Drucker-Prager yield function: F(\sigma)=\sqrt{J_s(\sigma)-Y(\sigma)} *
The Drucker-Prager material is required to satisfy the yield equation: F(\sigma)\leq 0
When F(\sigma) < 0 , the material behaves like a linear isotropic elastic material, with Lame paramters \lambda and \mu.

Wen F(\sigma) = 0, if the material is subjected to a strain that tends to cause an increase in F(\sigma), then the material yields. For TPV13, we assume that the material yields in shear. Yielding in shear means that when the material yields, the stress tensor \sigma_{ij} changes by an amount proportional to the stress deviator s_{ij}, so as to preserve the condition F(\sigma) with no change in mean stress \sigma_m .

Nucleation

TPV13 uses the same nucleation method as TPV12

Plasticity parameters

To turn on plasticity in SeisSol, add the following lines in parameter.par:

&SourceType

Plasticity = 1 ! default = 0 Tv = 0.03 ! Plastic relaxation /

In the material.yaml, add plasticity parameters:

!Switch

[rho,mu,lambda,plastCo,bulkFriction]: !ConstantMap
map:
rho: 2700
mu: 2.9403e10
lambda: 2.9403e10
plastCo: 5.0e6
bulkFriction: 0.85.

The format of yaml file can be found at

Results

Figure [fig:tpv13compare] shows the comparison between TPV12 (elastic) and TPV13 (plastic). The peak of slip rate in TPV12 is higher than TPV13. This difference attributes to the response of the off-fault plasticity. Refer to :raw-latex:`\citep{wollherr2018}` for detailed discussions.

SCEC TPV16/17

TPV16/17 has spontaneous rupture on a vertical, right-lateral, strike-slip fault in a homogeneous half-space with randomly-generated heterogeneous initial stress conditions. The earthquake rupture is artificially nucleated in a circular zone on the fault surface. The rupture then spontaneously propagates outward on the fault surface and encounters heterogeneous stochastic initial stress conditions,some of which prevent it from propagating into certain regions on the fault surface.

Geometry

The fault is a vertical, planar, right-lateral, strike-slip fault. The fault reaches the Earth’s surface. Rupture is allowed within a rectangular area measuring 48000 m along-strike and 19500 m down-dip. A node which lies exactly on the border of the 48000 m \times 19500 m rectangle is considered to be inside the rectangle, and so should be permitted to slip.

The portions of the fault below, to the left of, and to the right of the 48000 m \times 19500 m rectangle are a strength barrier, within which the fault is not allowed to rupture.

In the example, a vertical fault is generated with Gmsh in Figure [fig:tpv16mesh]. All the files that are needed for the simulation are provided in .

Material parameter

Rock properties are taken to be linear elastic throughout the 3D model volume. The problem description can be found at . Table [table:tpv16material] lists all the material parameters.

Parameter Description Value Unit
\lambda Lame’s first parameter 3.2044e10 Pa
\mu shear module 3.2038e10 Pa
\rho density 2670 kg/m^{3}
Q_p P-wave attenuation 69.3  
Q_s S-wave attenuation 155.9  
h_{edge} element edge length 200 m

Table: Table of bulk and material parameters in TPV16/17.

Mapview of fault randomly-generated initial stress in TPV16.

Mapview of fault randomly-generated initial stress in TPV16.

Nucleation parameters

Initial stress (Ts0) is randomly-generated in TPV16/17 (Figure [fig:tpv16ts]).

In TPV16/17, a two-stage nucleation method is used. The first stage is a circular zone of forced rupture which surrounds the hypocenter. Its radius is approximately 1 km (the exact radius is determined as part of the stochastic method that generates the initial stresses). At the hypocenter, the value of then increases with distance from the hypocenter, which creates an expanding circular region of forced rupture. The forced rupture expands at a speed of for 80% of the way, and then for the remaining 20% of the way to the edge of the zone. Outside the zone of forced rupture, is equal to 1.0E9, which means that forced rupture does not occur outside the zone.

The second stage is a circular zone of reduced which surrounds the hypocenter. Its radius is approximately 4 km (the exact radius is determined as part of the stochastic method that generates the initial stresses). In the innermost 10% of the zone, equals 0.04 m. The value of then increases linearly with distance from the hypocenter, and reaches its final value of 0.4 m at the edge of the zone. Outside the zone, equals 0.4 m. The effect is to create a circular region of reduced fracture energy surrounding the hypocenter, which helps the rupture to expand during the early part of the simulation.

Results

The earthquake nucleates and the rupture propagates on the fault surface due to the heterogenous stress ratio on the fault. Figure [fig:tpv16slip] shows the fault slip rate along strike-direction at T=5.5 s.

Mapview of fault slip rate along strike-direction.

Mapview of fault slip rate along strike-direction.

There are several receivers on the fault surface. Figure [fig:tpv16fault] shows slip rate along the strike- and downdip-direction on the fault at point (15 km, 0 km, -9 km).

Fault slip along strike- (left)and downdip- (right) direction.

Fault slip along strike- (left)and downdip- (right) direction.

SCEC TPV24

TPV24 is designed to illustrate dynamic rupture in a fault branching system. TPV24 contains two vertical, planar strike-slip faults; a main fault and a branch fault intersecting at an angle of 30 degree (Figure [fig:tpv24]). The earthquake rupture is artificially nucleated in a circular zone on the main fault surface and then spontaneously propagates to the branching fault.

Geometry

There are two faults, called the main fault and the branch fault (Figure [fig:tpv24]). The two faults are vertical, planar, strike-slip faults. The faults reach the earth’s surface.

The main fault is a rectangle measuring 28 000 m along-strike and 15 000 m deep. The branch fault is a rectangle measuring 12 000 m along-strike and 15 000 m deep. There is a junction point. It is located 12 000 m from the right edge of the main fault, and the main fault passes through it.

The branch fault makes an angle of 30 degrees to the main fault. The branch fault ends at the junction point.

The hypocenter is centered along-strike at a depth of 10 km in the left side of the main fault. That is, the hypocenter is 8000 m from the junction point, and 10 000 m deep.

Figure [fig:tpv24mesh] shows the fault model generated in Gmsh. The mesh file can be found at . The mesh can be generated following the detailed process in Section [sec:tpv5].

Initial stress

The initial stress condition is depth-dependent at the depth above 15600 m. Table [table:tpv24] summarizes the initial stress contidions in TPV24.

| p2cm | p4cm | p5cm| p2cm | Parameter & Description & Value & Unit
s_zz & :math:`sigma_{zz} `& -2670*9.8*depth & Pa
Pf & fluid pressure& 1000*9.8*depth & Pa
s_xz,s_yz & \sigma_{xz}, \sigma_{yz}& 0 & Pa
s_yy & & b33*(\sigma_{zz} + P_f) - P_f & Pa
s_xx & & b22*(\sigma_{zz} + P_f) - P_f & Pa
s_xy & & :math:`b23*(sigma_{zz} + P_f) `& Pa
s_yy & & \sigma_{zz} & Pa
s_xx & & :math:`sigma_{zz} ` & Pa
s_xy & & 0 & Pa
Note that the effective stress tensor is :

\bar{\sigma}_{effective}=
\begin{bmatrix}
&\sigma_{xx} + P_f , & \sigma_{xy} ,& \sigma_{xz}  \\
&\sigma_{xy}, &\sigma_{yy} +P_f , &\sigma_{yz} \\
&\sigma_{xz} ,&\sigma_{yz} ,  &\sigma_{zz} +P_f
\end{bmatrix}

Nucleation parameters

Nucleation is performed by forcing the fault to rupture, within a circular zone surrounding the hypocenter. Forced rupture is achieved by artificially reducing the friction coefficient, beginning at a specified time . The parameter specifies how long it takes for the friction coefficient to be artificially reduced from its static value to its dynamic value. So, the friction coefficient reaches its dynamic value at time . We reduce the friction coefficient gradually, over an interval of time, in order to smooth the nucleation process and reduce unwanted oscillations.

T = \left\{
\begin{array}{lr}
& \frac{r}{0.7Vr} + \frac{0.081*r_{crit} }{0.7Vr} (\frac{1}{1-(r/r_{crit})^2} - 1),  r \leq  r_{crit} \\
& 1E+09, r > r_{crit}\\
\end{array}
\right.

The cohesion zone is defined as :

C_0 = \left\{
\begin{array}{lr}
& 0.3 + 0.000675 * (4000 - depth),  depth < 4000 m \\
& 0.3 MPa, depth \geq 4000 m\\
\end{array}
\right.

Note that the frictional cohesion is 3.00 MPa at the earth’s surface. It is 0.30 MPa at depths greater than 4000 m, and its value is linearly tapered in the uppermost 4000 m.

The friction parameters are listed in Table [table:tpv24fric].

Parameter Description Value Unit
mu_s static friction coefficient 0.12  
mu_d dynamic friction coefficient 0.18  
d_c critical distance 0.30 m
C_0 fault cohesion   Pa
T forced rupture time   s
t_0 forced rupture delay time 0.5 s

Table: Table of LSR parameters on the fault in TPV24.

Results

The model is run for 12.0 seconds after nucleation. The earthquake rupture is artificially nucleated in a circular zone on the main fault surface. The rupture then spontaneously propagates on the main fault and encounters a branching fault. The branching fault continues to rupture as well as the rest main fault. The fault slip rate is shown in Figure [fig:tpv24result1].

|Snapshot of slip rate in branching fault system. Top: slip rate at 2 s. Bottom: slip rate at 3.5 s. | |Snapshot of slip rate in branching fault system. Top: slip rate at 2 s. Bottom: slip rate at 3.5 s. |

SCEC TPV29

TPV 29 constains a vertical, right-lateral fault with rough fault interface (Figure [fig:tpv29]). The fault surface has 3D stochastic geometrical roughness (blue and red colors). In TPV 29, the surrounding rocks respond elastically.

Geometry

The roughed fault interface model is generated with Gmsh is complicated than planar faults in previous sections. There are 5 steps to generate the model.

1. Download fault topography data from SCEC. There are 2001 nodes along the strike and 1201 nodes along the downdip. The node files should contains:
*Line 1: nx, ny Line 2 to nx: positions of nodes along the strike (in meters) Line nx+3 to ny+nx+3: positions of nodes along the downdip (in meters) Line to the end: fault topography of each nodes (nx*ny, in meters) *

Save this file as mytopo_tpv29, which can be found in .

2. Make a model with plane fault as Figure [fig:tpv29geo]. The Gmsh tpv29.geo file can be found at

3. Use gmsh_plane2topo.f90 and gmsh_tpv29.in to shift the planar fault according to positions given in mytopo_tpv29.

$ ./gmsh_plane2topo gmsh_tpv29.in
This will generate a step1_modified.msh file which containing rough fault surface.

4. Make a new step2.geo file that contains the new rough fault and mesh following general Gmsh process.

Material parameters

In TPV29, the entire model volume is a linear elastic material, with the following parameters listed in Table [table:tpv29material].

Parameter Description Value Unit
\rho density 2670 kg/m^{3}
\lambda Lame’s first parameter 3.2044e10 Pa
\mu shear module 3.2038e10 Pa
h_{edge} element edge length 200 m

[table:tpv29material]

Initial stress

The initial stress are listed in Table [table:tpv29fault].

Parameter Description Value Unit
mu_s static friction coefficient 0.12  
mu_d dynamic friction coefficient 0.18  
d_c critical distance 0.30 m
s_zz :math:`sigma_{zz} ` -2670*9.8*depth Pa
Pf fluid pressure 1000*9.8*depth Pa
s_xz,s_yz \sigma_{xz}, \sigma_{yz} 0 Pa
s_yy   :math:`Omega * b33*(sigma_{zz} + P_f) - P_f ` Pa
s_xx   :math:`Omega * b11*(sigma_{zz} + P_f) - P_f ` Pa
s_xy   :math:`Omega * b13*(sigma_{zz} + P_f) ` Pa

Table: Table of initial stress in TPV 29. b11, b33,b13 are 1.025837, 0.974162, −0.158649, respectively.

Note that the effective stress tensor is :

\bar{\sigma}_{effective}=
\begin{bmatrix}
&\sigma_{xx} + P_f , & \sigma_{xy} ,& \sigma_{xz}  \\
&\sigma_{xy}, &\sigma_{yy} +P_f , &\sigma_{yz} \\
&\sigma_{xz} ,&\sigma_{yz} ,  &\sigma_{zz} +P_f
\end{bmatrix}

where \Omega is defined as:

\Omega = \left\{
\begin{array}{lr}
&1, depth \leq 17000 m \\
& (22000 - depth)/5000 m,  17000 < depth < 22000 m \\
& 0, depth \geq 22000 m\\
\end{array}
\right.

Nucleation parameters

TPV29 use the similar strategy for dynamic rupture nucleation.

T = \left\{
\begin{array}{lr}
& \frac{r}{0.7Vr} + \frac{0.081*r_{crit} }{0.7Vr} (\frac{1}{1-(r/r_{crit})^2} - 1),  r \leq  r_{crit} \\
& 1E+09, r > r_{crit}\\
\end{array}
\right.

The cohesion zone is defined as :

C_0 = \left\{
\begin{array}{lr}
& 0.4 MPa + 0.000675 MPa * (4000- depth),  depth < 4000 m \\
& 0.4 MPa, depth \geq 4000 m\\
\end{array}
\right.

The friction parameters on the fault are listed in Table [table:tpv29fric].

Parameter Description Value Unit
mu_s static friction coefficient 0.12  
mu_d dynamic friction coefficient 0.18  
d_c critical distance 0.30 m
t_0 forced rupture delay time 0.5 s

Table: Table of friction parameters in TPV 29.

Results

The earthquake rupture is artificially nucleated in a circular zone on the fault surface.

SCEC TPV104

In this example, we illustrate how to implement rate-state friction law using a slip law with strong rate weakening (RS-SL-SRW) and setup parameters in SeisSol.

TPV104 has a planar rectangular vertical strike-slip fault with the main rupture region of velocity-weakening friction, a zone on the fault surface with transitional friction surrounds the main fault rupture region, and the outer regions on the fault surface have velocity-strengthening friction (Figure [fig:tpv104]).

Geometry

TPV104 uses a vertical fault same as TPV5. We use the mesh file of TPV5 directly.

RSL parameters

TPV104 uses rate-state friction where shear stress follows:

\begin{aligned}
   \tau = f(V,\psi) \sigma\end{aligned}

The friction coefficient is a function of slip rate V and state \psi:

\begin{aligned}
    f(V,\psi) = a *  arcsinh [\frac{V}{2V_0} \exp(\frac{\psi}{a})]\end{aligned}

The state variable evolves according to the equation:

  \begin{aligned}
  \frac{d \psi}{dt} = - \frac{V}{L}[\psi - \psi_{ss}(V)]\end{aligned}
and
\begin{aligned}
   \psi_{ss}(V) = a \ln [\frac{2V_0}{V} \sinh (\frac{f_{ss}(V)}{a})]\end{aligned}

f_{ss}(V) is the stead state friction coefficient that depends on V and the friction parameters f_0, V_0, a, b, f_w and V_w.

\begin{aligned}
f_{ss}(V) = f_w + \frac{f_{LV}(V) - f_w}{[1+(V/V_w)^8]^{1/8}}\end{aligned}

with a low-velocity steady state friction coefficient:

\begin{aligned}
f_{LV}(V) = f_0 + (b-a) * \ln (V/V_0)\end{aligned}

In SeisSol input file, Rate-state friction law can be used by choosing FL=103 in parameter.par (Section ). The friction parameters of RS-SL-SRW are shown in Table [table:tpv104rsl].

| p3.5cm | p7cm | p2cm| p2cm | Parameter & Description & Value & Unit
RS_f0 & reference friction coefficient & 0.6 &
RS_a & direct effect & 0.01 &
RS_b & evolution effect & 0.014 &
RS_sr0 & reference velocity scale & 1d-6 & m/s
RS_srW & weakening slide rate & 0.1 & m/s
RS_sl0 & critical slip length & 0.4 & m
Mu_W & weakening friction coefficient & 0.2 &
RS_iniSlipRate1 & initial sliding velocity & 1d-16 & m/s
RS_iniSlipRate2 & initial sliding velocity & 0 & m/s
t_0 & forced rupture decay time & 1 & s
Rs_a & direct effect & 0.02 &
RS_srW & weakening slide rate & 1.0 & m/s

To stop the rupture, the friction law changes from velocity-weakening in the rectangular interior region of the fault to velocity-strengthening sufficiently far outside this region. The transition occurs smoothly within a transition layer of width w = 3 km. Outside the transition layer, the fault is made velocity-strengthening by increasing a by \triangle a= 0.01 and V_w by \triangle V_{w0} = 0.9 . The exact format can be referred to .

Results

The earthquake nucleates in the velocity-weakening zone spontaneously. The rupture propagates through the transition zone into the velocity-strengthening region, where it smoothly and spontaneously arrests. Nucleation is done by imposing additional shear stress in a circular patch surrounding the hypocenter.

Figure [fig:tpv104sr] shows the slip rate on the fault along the downdip direction at T=5s.

|Slip rate along-strike on the fault at 2 s(top) and 5 s (bottom) of TPV 104. | |Slip rate along-strike on the fault at 2 s(top) and 5 s (bottom) of TPV 104. |

Point Source

SISMOWINE WP2_LOH1

SISMOWINE is intended as a long-term interactive web interface for verifying numerical modeling methods in seismology. Numerical-method developers and numerical modelers may compare their solutions with other solutions. SISMOWINE is a continuation of the original SPICE Code Validation interface established within the 6th Framework Programme project .

LOH1 is used as an example here to illustrate the implementation of source point for earthquake nucleation in SeisSol. The details of LOH1 model can also be found at .

The model uses Right-handed Cartesian, x positive North, y positive East, z positive downward, all coordinates in meters. The source is buried at 2000 m in a half-space Earth (Figure [fig:loh1]. The top layer is 1000 m thick and the bottom layer is 33000 m. The material parameters are listed in Table [table:loh1].

  Vp (m/s) Vs(m/s) density Qp Qs
layer 4000 2000 2600 Inf Inf
half-space 6000 3464 2700 Inf Inf

Table: Material properties in LOH1 .

Geometry of LOH1 .

Geometry of LOH1 .

Geometry

The mesh is generate using Gmsh.

Point source input

The point source needs to be turned on in parameter.par file.

&SourceType Type = 50 FileName=’LOH1_source.dat’ /

The source input file can be found at . Duration of the source is 4 seconds.

Results

The comparison with solution is shown in Figure [fig:compare_loh1].

Kinematic source example - 1994 Northridge earthquake

We use this earthquake to demonstrate how to setup dynamic rupture model with kinematic rupture source in SeisSol.

The 1994 Northridge earthquake occurred on January 17, at 4:30:55 a.m. PST and had its epicenter in Reseda, a neighborhood in the north-central San Fernando Valley region of Los Angeles, California, USA. It had a duration of approximately 10–20 seconds. The blind thrust earthquake had a magnitude of 6.7 (Mw). This is a typical reverse-slip earthquake. The fault orients to N122^\circE and dips at 40^\circ. The simulation can be used to build similar model with moderate modifications.

Geometry

The fault geometry is made in Gmsh. Fault: plane fault 20 km*25 km dipping at 40-degree.

Region: 100 km*100 km *60 km.

Kinematic rupture Source

The kinematic source of the earthquake can be found at . The standard rupture format can be used directly in SeisSol, with the following lines in parameter.par file.

&SourceType Type = 42 FileName=’northridge.nrf’ /

Download standard rupture format file (northridge.srf) can be found in . Please note that the SCEC units are different with SeisSol units in some aspect.

The fault are divided in to 20 grids along the strike and 25 grids along the dip. The source time function (STF) of each rectangular elements is given in the file , whose format looks like the following:
verison (1.0) PLANE 1 ELON ELAT NSTK NDIP LEN WID STK DIP DTOP SHYP DHYP POINTS 500 LON LAT DEP STK DIP AREA TINIT DT RAKE SLIP1 NT1 SLIP2 NT2 SLIP3 NT3 SR1[1] SR1[2] SR1[3] . . . SR1[NT1] SR2[1] SR2[2] SR2[3] . . . SR2[NT3] SR3[1] SR3[2] SR3[3] . . . SR3[NT3] ...
Explanations:
Line 1: version
Line 2: Number of fault planes
Line 3:
ELON top center longitude
ELAT top center latitude
NSTK number of point sources (subfaults) along strike
NDIP number of point sources (subfaults) down-dip
LEN segment length (km)
WID segment width (km)
STK segment strike
DIP segment dip
DTOP depth to top of fault segment (km)
SHYP along strike location (from top center) of hypocenter for this segment (km)
DHYP down-dip location (from top edge) of hypocenter for this segment (km)
Line 4: Number of points per fault plane
Line 5-9:
LON: longitude of subfault center
LAT: latitude of subfault center
DEP: depth of subfault center (km)
STK: strike
DIP: dip
AREA: area of subfault (cm2)
TINIT: initiation time when rupture reaches subfault center (sec)
DT: time step in slip velocity function (sec)
RAKE: direction of u1 axis (rake direction)
SLIP1: total slip in u1 direction (cm)
NT1: number of time points in slip rate function for u1 direction
SLIP2: total slip in u2 direction (cm)
NT2: number of time points in slip rate function for u2 direction
SLIP3: total slip in u3 (surface normal) direction (cm)
NT3: number of time points in slip rate function for u3 direction
SR1[1],…,SR1[NT1] slip rate at each time step for u1 direction (cm/sec)
SR2[1],…,SR2[NT2] slip rate at each time step for u2 direction (cm/sec)
SR3[1],…,SR3[NT3] slip rate at each time step for u3 direction (cm/sec)
Project geographic coordinates

The geographic coordinates of source model is projected to Cartesian coordinates wit the pre-processing tool rconv.

rconv -i northridge.srf -o northridge.nrf -m “+proj=merc +lon_0=-118 +y_0=-4050981.42 +x_0=57329.54 +units=m +axis=enu” -x visualization.xdmf

To find the center of fault, use cs2cs in proj.4 to convert the cooridinates:

echo -118.5150 34.3440 0.0 | cs2cs +proj=lonlat +axis=enu +units=m +to +proj=merc +lon_0=-118 +axis=enu +units=m

This cooperation will project the coordinates and shift the center of fault to the origin (0,0) in Cartesian coordinates.

Results

Source rupture starts at 7.0 s and propagates in the domain. A snapshot of velocity is show in Figure [fig:northridge1]. The surface velocity output is refined by subdividing each triangle into 4 subtriangles while the domain output is not.

Acknowledgements

Copyrights

The copyrights belong to seismology group @LMU