3D MODELING OF DYNAMIC RUPTURE
PROPAGATION BASED ON A BOUNDARYINTEGRAL EQUATION METHOD FOR
EARTHQUAKES IN CHINA
Zhang LifenMEE10501

Supervisor:

Taku TadaBunichiro Shibazaki

ABSTRACTWe developed a code for dynamic rupture propagation on arbitrarily shaped faults based on the stress
Green’s function corresponding to a triangular mesh fault element. A boundary integral equation
method is adopted in this study. Starting from representative theorem, after regularization anddescretization of a basic elastodynamic equation relating stress with slip, a more applicableformulation used in dynamic rupture simulation is constructed. For a boundary condition, we utilize asimple but typical slip weakening friction law. In order to avoid numerical oscillation, according to
Courant-Friedrichs-Lewy (CFL) condition, time intervals used in the simulation is given as ∆t ≤2l 6CL (l is the element size, and CL denotes P wave velocity). And then we applied the code to aplanar fault and a curved fault.
For a simple planar fault, spontaneous rupture propagation initiates from a given circularasperity and then propagates symmetrically and bilaterally around the central point. A bump of slip
distribution near the center of the fault is the most characteristic feature of the fault that starts from afinite asperity. In addition, we considered a curved fault with different final inclination angles. The
curved fault with a small inclination angle does not have significant effects on the rupture propagation.
The rupture still can expand symmetrically and bilaterally beyond a bend. However, a large inclination
angle may cause the rupture to propagate asymmetrically. Initial conditions determine the rupture
propagation behavior. Given the same initial stress on the same fault, the smaller the critical slip
weakening displacement is, the faster and longer the rupture propagates. Initial stresses inside andoutside of the asperity (T asp, Te) also have great effects on the rupture propagation. If T e is too small,the rupture will stop soon after initiation. If T e becomes larger, the rupture can propagate longer andfaster with subshear velocity.
Keywords: Boundary integral equation method, Triangular fault mesh, Dynamic rupture propagation,
Slip weakening friction law

1. INTRODUCTION
It has often been observed that earthquake rupture propagates along pre-existing faults with complex
geometry in the continental crust. These fault complexities, sometimes referred to as nonplanar faults,
include fault bending, jogs, step-overs, branches and so on. In order to understand how these complex
fault geometries affect the rupture propagation, dynamic numerical simulation provides insight into
physical effects of nonplanar fault geometry on the earthquake rupture and slip process. In addition, aEarthquake Administration of Hubei Province, China.
Asahi Shimbun General Service, Japan.

Building Research Institute, Japan.
1physical understanding of the earthquake rupture process can improve our capability for predicting
ground motion, and therefore our assessment of seismic hazard (Duan and Oglesby, 2005).
Many numerical methods have been used for dynamic simulation. Finite difference method(FDM) and finite element method (FEM) are so-called field domain methods, which are applicable toarbitrary geologic structures and fault geometries, but require a field discretization within the domain
of interest (Goto and Bielak, 2008). In contrast, boundary integral equation method (BIEM) is more
powerful, which is a boundary based method. Because a problem is solved not in the whole volume
but only on the fault surfaces, it simplifies the problem and requires small resources in computation. In
addition, for fault rupture problem, it is easy to consider the singularity at the tip of crack using BIEM.
Therefore, in the present paper, we adopt this method to describe the dynamics of 3-D cracks with
arbitrary geometry. Compared with rectangle meshes, triangle meshes are more appropriate to describe
the complex fault geometry. However, few researchers model the dynamic rupture propagation using
Green’s function calculated with triangular fault meshes. Therefore, in the present study, we willdevelop a code to do the 3-D dynamic rupture simulation using triangular fault meshes.

2. THEORY AND METHODOLOGY
2.1. Boundary Integral Equation Method
Boundary integral equation starts from representative theorem (Aki and Richards, 1980), which states
the relationship between displacement at receiver location and slip across a fault. However, for a
dynamic problem, stress should be taken into consideration. Based on this theorem, a basic
elastodynamic equation relating stress with slip is constructed. For the 3-D case, slip is assumed tooccur along the x1- and x2-axis, so the stress along the x 1- and x2-axis can be represented as:t

T1 x, t= T10 − μ2dS ξdτ ∆u1 ξ, τ0Γ+ ∆u2 ξ, τ

∂2
∂2
∂2
G
+
2
G
+
G
∂x1 ∂x3 13 ∂x12 33∂x32 11

∂2
∂2
∂2
∂2
G
+
G
+
G
+
G
12
13
23
∂x2 ∂x3∂x1 ∂x3∂x1 ∂x2 33
∂x32

(1)
where Gij denotes Green’s function, Δu1(ξ,τ) and Δu2(ξ,τ) denotes the slip along the x1- and x2-axis, µ
is the rigidity, Γ is the fault surface, T10 is the initial stress along x1-axis. And T2(x, t) has symmetricalform with T1(x, t).
2.2. Regularization and Discretization
The traction boundary integral equations, thus derived, are hypersingular, and cannot be directly usedin numerical modeling. Tada et al. (2000) solved this difficulty by enlarging Fukuyama andMadariaga's (1998) integration by parts technique to derive a non-hypersingular time-domainboundary integral equation for the 3-D non-planar crack in a fully explicit form. However, theboundary equation is still too abstract to be utilized in numerical simulation. Therefore, discretizationis needed. When discretizing the fault surface, in general, the size of the mesh elements may be either
homogeneous or variable along the fault, but the time step intervals are usually taken to behomogeneous. Additionally, based on collocation method and piecewise constant approximation
method, the boundary integral equation can be discretized to be a more applicable form (Eq. (2)). For
the simulation, we put receiver at the collocation point of each element. And the time collocation point
had better be put at the end point of each time intervals.
𝑁𝑋𝑝

𝑇𝑚𝑝 = 𝑇0𝑚 +𝑡

𝑉𝑛𝑞
𝑛=1 𝑝=1𝑑𝑆 𝜉Γ𝑑𝜏 𝐾 𝑥𝑚 , 𝑡𝑝 − 𝜏; 𝜉, 0 H 𝜉 ∈ Γ𝑛 H 𝜏 − 𝑡𝑞 −1 − H 𝜏 − 𝑡𝑞
0

(2)

2

where Vnq is the discrete slip rate, H(𝜉 ∈ Γ𝑛 ) and H 𝜏 − 𝑡𝑞 −1 − H 𝜏 − 𝑡𝑞 denote that if slip occurs
within Γ𝑛 during a certain time intervals, slip rate is unity, otherwise it is zero. 𝐾 is the integration
kernel, which is a response function. During the simulation, the calculation of the kernel is of great
importance.
2.3. Numerical Stabilization
An exponentially growing numerical oscillation in BIEM is well known in the BIEM community(Tada and Madariaga, 2001). It has been recognized that the oscillation is coming from thespatiotemporal discretization, therefore, an appropriate and optimum combination of space grid andtime intervals is very important, which can suppress the stability to some extent. We adopt thecommon method of Courant-Fridrichs-Lewy condition. Due to the triangular fault meshes, the timeinterval (Δt) is chosen as ∆t ≤ 2𝑙 6CL (l denotes the element size, and C L denotes the P wave
velocity).

3. DATA AND VALIDATION
3.1. Coordinate Transform
For the planar fault, we can describe its geometry just using global Cartesian coordinate system.
However, it is a little difficult to define the nonplanar fault. Instead, Tada (2009) solved this difficultyby making use of a local Cartesian coordinate system.
We assume that slip occurs on a triangular fault element, with its vertices named A, B and C,
which is embedded in a 3D infinite, homogeneous and isotropic medium. For convenience, the local
coordinate system is defined according to the following rules:(1) The system of local Cartesian coordinates (x1 , x2 , x3 ) is right handed, and shares the origin with thesystem of global Cartesian coordinates (X1 , X2 , X3 ), which is also right handed.
(2) The x3 -axis is perpendicular to the fault element and the direction is upward. Vertices A, B and C
are aligned counter-clockwise in this order when seen in the +x3 -direction.
(3) The x1 -axis can be defined in any arbitrary direction as long as it is normal to the x3 -axis, but inour case, we define it parallel to the slip direction on the fault. Totally, there are five coordinate
systems used to calculate the stress Green’s function.
3.2. Validation for the Green’s function
After obtaining the corresponding local
coordinates, to validate the accuracy of the stresscomputation, we calculated the different
components of the Green’s function Lij/k caused bya constant slip rate of unit magnitude taking placein the xk-direction across triangles ABC or CDA
at receiver location R(0.6,0.3,0.2). By comparison
with Tada (2006), as shown in Figure 1, ourresults are in very good agreements. In thisvalidation, triangular faults ABC and CDA, withvertices located at A(0,0,0), B(1,0,0),C(1,1,0) andD(0,1,0) in the local x1x2x3-coordinate system.
The shear wave velocity is 1.0 km/s, and the P
wave velocity is 1.73 km/s.

Figure 1. Comparison of the different componentsof Green’s function computed using the 3-D
kernels of BIEM by this study (the lower figures)
and Tada (2006) (the upper figures).

3

4. RESULTS AND DISCUSSION
For dynamic rupture simulation, we need to evaluate slip, slip velocity and stress. As a boundary
condition, we adopt simple but typical slip weakening friction law. We apply the code to a planar anda nonplanar fault with different inclination angles.
4.1. Planar and Nonplanar Fault Model
Table 1. Parameters used in dynamic simulation
The fault geometries of planar andnonplanar fault are shown in Figure 2.
It is well known that in order to initiate
rupture in such a model, an initial highstress patch must be ready to break(Andrews, 1976).
We study spontaneous
rupture of a fault starting from acircular asperity that is ready tobreak. It is loaded with stress that islarger than the peak strength andsurrounded by a bounded fault surface
where the initial stress is smaller thanthe peak strength. The parameters arelisted in Table 1. In this table, theFigure 2. Schematic fault geometries of a planar and a mesh size is chosen as 1.0 km, andnonplanar fault. Tasp and Te denote the stress inside and considering the numerical stability, theoutside of the asperity. T p is the peak strength. Figure 2(a) time intervals is 0.045 second.
and Figure 2 (b) show the planar and nonplanar fault model.
4.2. Rupture Propagation on a Planar Fault
For a planar fault, the rupture initiates from the circular asperity and then propagates symmetricallyand bilaterally. Figure 3 shows the slip distributions on the fault plane as a function of time andposition along the x axis of the fault. We can observe the rupture starts from the localized asperity
around the center of the fault and then grows bilaterally at a slightly increasing slip rate. Comparing
my result with Madariaga et al. (1998), we can find the consistence of a bump at the center of the fault.
This localized excess slip is the most characteristic feature of faults that start from a finite asperity,
which is common phenomenon for most slip functions (Madariaga et al. 1998). However, if the initial
radius of asperity is too small, the bump is insignificant.

Figure 3. Comparison of spontaneous slip of the bounded rectangular fault.
(Left -Madariaga et al (1998), Right -Present study)44.3. Rupture Propagation on a Nonplanar Fault
For a nonplanar fault, the bending position is put at x=14.14 km, and we consider two cases withdifferent inclination angles of 10˚ and 30˚. The series of snapshots in Figure 4 (a) and Figure 4(b)show the evolution of slip velocity for curved fault model.
For both cases, the ruptures
expand symmetrically on the
primary fault until the rupture
front reach the bend. In the caseof the small bend with the final
inclination angle of 10˚ (Figure4(a)), the rupture can propagate
symmetrically beyond the bend,
which is similar with planar
fault rupture propagation. Whenthe final inclination angle
increases to 30˚ (Figure 4(b)),
at 2.25 seconds, the rupture
front reaches the bending
boundary of the fault. However,
the rupture cannot propagate
symmetrically any more beyond
the bend. It illustrates thatcurvature fault with big final
inclination angle hinders the
rupture propagation, which
Figure 4. Snapshots of slip velocity for the nonplanar fault with behaves like a barrier to someinclination angle of 10˚ and 30˚. (Unit of horizontal and vertical extent.
axis: km)4.4. Influence of Initial Conditions on Rupture Propagation
In the simulation, we givedifferent initial conditions to
examine their effects on therupture propagation. When theinitial stresses inside andoutside of the asperity are
given, it is easy to find that thesmaller the slip weakening
displacement (Dc) is, the fasterandlongertherupturepropagates (Figure 5(a)). If Dc
is too large, the rupture easily
terminates(Figure5(c)).
Another important factor is thestress outside of the asperity(Te), by comparison, it is
Figure 5. Snapshots of stress for the different initial conditions on illustrated that small Te easily
causes the rupture to stop
the rupture propagation. (Unit of horizontal and vertical axis: km)(Figure 5(b)).

55. CONCLUSIONSThe purpose of this study is to simulate the complex fault dynamic rupture propagation based on theboundary integral equation method. In order to attain the goal, we developed a code and validated it.
We applied the code to the spontaneous rupture on a rectangular planar fault and curved fault in an
infinite homogeneous medium, starting from a circular asperity under the slip weakening friction
criterion. By analyzing the different cases, we had some findings.
For the planar rectangular fault with homogeneous initial stress field, the rupture initiatesfrom the circular asperity and then propagates symmetrically and bilaterally. The bump of the slip
distribution (excessive slip) near the center of the fault is the most characteristic feature of the fault
that starts from a finite asperity.
For the curved rectangular fault with different inclination angles, the small bend angle does
not have great influence on the rupture propagation. Instead, the large bend angle may behave like abarrier.
The initial conditions have great influences on the rupture propagation. Given the same
initial stress on the same fault, critical slip weakening displacement (D c) determines the rupture
propagation. The smaller the D c is, the faster and longer the rupture propagates. The initial stresses
inside and outside of the asperity (Tasp, Te) also have great effects on the rupture propagation. If the T e
is too small, the rupture will stop soon after initiation. If the T e becomes larger, the rupture can
propagate longer and faster with sub-shear velocity.
6. RECOMMENDATIONThe boundary integral equation method is good at dealing with nonplanar fault problems. However, it
has only limited applicability to a full space and a homogeneous medium. But for real earthquake
researches, the free surface has some effects on the rupture propagation. In addition, for the dynamic
simulation, it takes much time and computer memory. Therefore, according to above limitations, werecommend the following suggestions:i.
In order to compensate the weakness of boundary integral equation method, it is essential toresort to the advantages of other numerical methods. For the strike-slip faults, we can resort to
the imaginary mirror sources to consider the influence of free surface. For other faulting types,
we also should find the solutions.
ii.
More appropriate algorithm and flexible program are required for modeling dynamic rupture
propagation.

REFERENCES
Aki, K. and Richards, P.G., 1980, W.H. Freeman, San Francisco.
Andrews, D. J., 1976, J. Geophys. Res., 81, 5679-5687.
Duan, B,C. and Oglesby, David D., 2005, J. Geophys. Res., 110, B03304
Fukuyama, E., and Madariaga, R., 1998, Bull. Seism. Soc.Am. 88, 1-17.
Goto, H., and Bielak, J., 2008, Geophys. J Int., 172, 1083-1103.
Madariaga, R., Olsen, K., and Archuleta, R., 1998, Bull. Seism. Soc. Am. 88, 1182-1197.
Tada, T., Fukuyama, E., and Madariaga, R., 2000, Computational Mechanics, 25, 613-626.
Tada, T., 2006, Geophysical Journal International, 653-669
Tada, T., 2009, International Geophysics, Elsevier Inc , USA.

6