1
3D hp

Adaptative Discontinuous
Galerkin
Method
for
Modeling
Earthquake
Dynamics
Josué
Tago
1
, V
íctor
M. Cruz

Atienza
1
,
J
ean
Virieux
2
,
V
incent
Etienne
3
and F
rancisco
J.
Sánchez

Sesma
4
1
Instituto de Geofísica, Universidad Nacional Autónoma de México, Mexico
2
Institut de Sciences de la Tèrre, Université Joseph Fourier, France
3
UMR Géoazur, Université de Nice
–
Sophia Antipolis, France
4
Instituto de Ingeniería, Universidad
Nacional Autónoma de México, Mexico
Submitted to Journal of Geophysical Research
March 14, 2012
2
Abstract
1
2
We introduce a novel scheme
, the DGCrack,
to simulate the dynamic rupture of earthquakes in
3
three dimensions
based on an hp

adapt
at
ive
d
iscontinuous Galerkin
method
. We solve the
4
velocity

stress weak formulation of elastodynamic equations in an unstructured tetrahedral
5
mesh
with arbitrary mesh refinements
(h

adaptivity)
and local approximation or
ders (p

6
adaptivity).
O
ur scheme considers second

order fault elements (P2) where dynamic

rupture
7
boundary conditions are verified through
ad hoc
fluxes
across the fault
.
T
o
model
the
Coulomb
8
slip

dependent
friction law, we introduce a predictor

corrector
scheme
for
estimating
shear
fault
9
tractions
,
in addition to
a special treatment of
the
normal traction
s
that
guarantee
s
the
10
continuity of
fault normal
velocit
ies
. We verify
the
DGCrack
by compari
son
with
several
methods
11
for tw
o spontaneous rupture test
s
and find excellent results
(i.e.
rupture times
RMS errors
12
smaller than
1.0
%)
provided
that
one
or
more fault elements
resolve
the
fault
cohesi
ve
zone.
To
13
make this comparison clear, we introduce a methodology based on cross

correlation
14
measurements that p
rovide a simple way to
quantify
the similarity between
solutions.
Our
15
verification tests include
a
6
0
o
dipping normal fault reaching the free surface
.
T
he
DGCrack
16
method reveals
power convergence rates
close
to those of well

established methods and a
17
numerical efficiency
significantly higher than that of
similar
discontinuous
Galerkin
approaches
.
18
W
e
appl
y
the method
to
the
1992 Landers

earthquake fault
system
in a layered medium,
19
considering
heterogeneous initial
stress
conditions and
mesh adaptivites
.
Our results
show that
20
previously proposed dynamic models
for
the Landers earthquake require a reevaluation in
21
terms of the initial stress conditions.
22
23
3
Key points
: Dynamic rupture, discontinuous Galerkin, unstructured mesh, Landers earthquake
24
25
1.
Introduction
26
27
The availability
of high

quality near

field
records
of
large
subduction
earthquakes in the
last
few
28
years
makes it possible
,
like
never before,
to test
and
validate
physic
s

based rupture models.
The
29
d
evelopment
of
sophisticated
models
to expl
ain such
aggregate
of observations is now
largely
30
justified
. Huge efforts have been made by the seismological community in the last ten years to
31
overcome technical limitations preventing most dynamic
source
models
,
at the end of
the
90s
,
32
from integrating
the
effect
of
fault geometry
in the spontaneous rupture of earthquakes
.
Because
33
of both its simplicity and efficiency, t
he finite difference
(FD)
method has been one of the first
34
and
most
persistent approaches
to simulate rupture dynamics
along planar faul
ts
(e.g. Andrews,
35
1976; Madariaga, 1976
;
Miyatake, 1980
; Day, 1982; Virieux and Madariaga, 1982
;
Harris and
36
Day, 1993;
Madariaga et al., 1998;
Peyrat et al., 2001;
Day et al., 2005;
Dalguer and Day, 2007
).
37
Although
different
strategies have been proposed
in
recent
years
to integrate com
plex fault
38
geometries into this
kind of method (Cruz

Atienza and Vi
r
ie
ux,
2004;
Kase and Day, 2006;
Cruz

39
Atienza et al., 2007
;
Kozdon et al.,
2011)
today
, mo
st common approaches handle
numerical
40
lattices
(meshes)
that are
ad
aptable to the
problem geometry (i.e.
fault geometry
)
.
One group
of
41
methods is based on
boundary integral
equations
(BIE)
and
have
also
been used
for longtime
42
(
e.g.
Das and Aki, 1977; Andrews, 1985;
Cochard and Madariaga, 1994
; Kame and Yamashita,
43
1999;
Aochi et al., 2000; Lapusta et al., 2000; Hok and Fukuyama, 2011
)
.
However, s
ince these
44
methods
discretize only boundaries
and require
semi

analytical approximations of Green
45
functions, they have difficulties
integrating
heterogeneities
of the bulk propert
ies
where the fault
46
4
is embedded
.
The other group
consists of
domain methods based on
weak formulation
s
of the
47
elastodynamic equations
,
and
can be separated in
to
two subgroups
depending on the way
the
48
lattice
boundaries are treated
.
On one hand t
he
continuo
us
finite element
methods
(FEM)
,
whose
49
formulations
require
continuity
between the
mesh
elements
except where special treatments of
50
boundary conditions are imposed
(
e.g. Oglesby and Day, 2000;
Ampuero, 2002; Festa and Vilotte,
51
2006;
Ma
and Archuleta, 2006;
Kaneko et al., 2008
; Ely et al., 2009
; Barall
, 2009
) and
,
on the
52
other
,
the d
iscontinuous
finite element methods,
better
known as the discontinuous
Galerkin
53
(DG)
methods,
which
only consider fluxes between elements and, therefore,
do
not impose any
54
field
continuity across the
ir
boundaries
.
When studying
the
earthquakes
source
physics
, the
55
discontinuity
produced
across the fault
by the rupture
process
must be accurately treated
,
so
56
that the DG strategy is naturally
suitable
for
tackl
ing
this problem.
57
58
The
fi
rst dynamic rupture model
based on a DG approach
was introduced
in two dimensions
59
(2D)
by Benjemaa et al. (2007)
for low

order (P0) interpolation functions
.
In this case,
where the
60
basis functions are constants,
the DG schemes are also known as finite vo
lume (FV) methods
61
(Le
V
eque
, 2002
) and provide
c
omputational
ly
efficient
algorithms
as fast as
second order FD
62
schemes
(i.e. they are
equivalent in efficiency
on rectangular meshes)
.
However, t
he extension to
63
three dimensions (3D) of this
model
(Benjemaa et al., 2009)
revealed to
have convergence
64
problems
for
unstructured tetrahedral grids
(e.g. non

planar faults)
(Tago et al., 2010)
.
In these
65
irregular grids, P0 elements
have
zero

order convergence for wave propagation modeling
due to
66
the center
ed flux approximation
(
Brossier et al, 2009; Remaki et al, 2011
)
, so increasing the
67
element interpolation order to achieve a proper numerical convergence of wave propagation
68
with a DG scheme is mandatory.
Nevertheless,
it is
also
know
n
that
the
convergence rate for
the
69
5
dynamic

rupture
numerical problem
is
in
sensitive to
high
interpolation
order
s
(i.e. 4th order or
70
higher)
, and that second order
interpolation
methods
are often the
most accurate
and efficient
71
approximation
s
for applying
the corresp
onding fault
boundary conditions
(
Cruz

Atienza et al.,
72
2007; Moczo et al., 2007;
Rojas et al., 2009;
Kozdon et al., 2011
).
A
notable
case where the
73
convergence rate is essentially not sensitive to increments in the
interpolation
order is the
74
ADER

DG discon
tinuous Galerkin method
for
2D and 3D
geometries
by de
l
a
Puente et al. (2009)
75
and Pelties et al. (2012)
, respectively
,
despite
its
spectral convergence for the wave propagation
76
problem (Dumbser and Kaser, 2006). The ADER

DG is based on a modal interpolation
77
formulation
,
instead of the nodal interpolation we
consider
here. Both
formulations
are
78
mathematically equivalent but computatio
nally different
(Hesthaven and Warburton, 2008)
.
O
ur
79
choice
of using
the nodal
approximation
essentially relies in the fact that
the evaluation
of fluxes
80
requires
fewer computations
than
in
a modal
scheme
, as we shall explain on Section 4.1
.
81
82
In this work we introduce a
novel
discontinuous Galerkin approach
, namely the
DGCrack
83
method,
to model the dynamic rupture of earthquakes in 3D
geometries
.
The numerical platform
84
of o
ur model
is the G
eo
DG3D
code
(Etienne et al., 2009), so that it
accounts
for multitask
85
computing using the Message Passing Interface (MPI) with
~
85
% scalability
, free surface
86
boundary conditions
with arbitrary geometry
and
Convolutional
Perfectly Matching Layer
87
(
C
PML) absorbing boundary conditions
along the external edges of t
he computation domain
88
(Etien
n
e et al., 2009
and references there
in
)
.
In this approach,
intrinsic
attenuation
via the rock
89
quality factors Q
has been introduced as well
,
but will not be discussed further in this work
90
(Tago et al., 2010)
.
To maximize both
the efficiency and the accuracy
of the scheme
depending on
91
the model properties and geometry
,
t
he method
handles
unstructured mesh
refine
ments
(i.e.
h

92
6
adaptiv
ity)
and
locally
adapts the order of the
nodal
interpolations
(i.e.
with
in every grid
93
element;
p

a
daptivity
)
(Etien
n
e et al., 2009)
.
We first introduce the mathematical and
94
computational
concepts
for the 3D dynamic rupture problem
,
and then assess
both
its accuracy
95
and convergence rate
by comparing
calculated
solutions with those yielded by
finite d
ifference
96
(
DFM
)
, finite element (FEM)
, spectral boundary integral (MDSBI)
and
spectral element
97
(
SPECFEM3D
) methods for two spontaneous rupture benchmark
tests
(Harris et al., 2009
).
Since
98
one of our
major goal
s
in the near future
is
the investigation of dy
namic
rupture propagation
99
along
realistic (nonplanar) fault geometries, we take special care to verify
the accuracy of the
100
highly singular normal stress field across the fault during rupture propagation, as the fault
101
normal tractions
strongly determine
the
radiated energy throughout the Coulomb
failure
102
criterion
. We finally illustrate
the capabilities of the
DGCrack
method
through
a spontaneous
103
rupture
simulation
along the 1992 Landers earthquake fault system, which is a geometrically
104
intricate and physical
ly interesting study case.
105
106
2.
Elastodynamic equations
107
108
V
elocit
ies
and stresses
induced by the propagation of waves
in a homogeneous elastic medium
109
can be modeled
with
a first order hyperbolic system
(Madariaga, 1976)
.
F
ollow
ing
the
110
transformation proposed by Ben
j
emaa et al. (2009)
, a
pseudo

conservative form
of the system is
111
given by
112
113
∑
(
)
{
}
7
∑
(
)
{
}
(
1
)
114
115
where
(
)
is the velocity vector and
(
)
is the stress
116
vector.
T
o avoid physical properties in the right hand side of
(
1
)
,
a change of variables
is
done in
117
the stress vector leading to its first three components
118
119
(
)
(
)
(
)
,
(chg_var)
120
121
which involve the
mean
and the deviatoric stresses.
122
123
In this model, t
he physical properties of the medium are the density
,
,
and
the
medium matrix
124
(
)
,
which is composed by the Lamé parameters
and
. The
125
external forces and the initial stresses are defined as
and
,
respectively
,
while
the terms
126
and
are constant real matrices
whose
definition can be
found
in Benjemaa et al. (2009).
127
128
2.1 hp

adaptive Discontinuous Galerkin scheme
129
130
8
To make the local approximation of the hyperbolic system
(
1
)
, we first need to
decompose the
131
domain
into
elements
, so that
132
133
∑
(
2
)
134
135
w
here each
is a straight

sided tetrahedron and the mesh
(i.e. decomposed domain)
is
136
assumed to be geometrically conforming.
By applying a
Discontinuous Galerkin approach to t
he
137
weak formulation of
(
1
)
,
as
proposed by Etienne et al. (2010)
,
we obtain
the following
velocity

138
stress
iterative
scheme
in
every

tetrahedron
,
139
140
(
)
∑
(
)
{
}
∑
[
(
)
(
)
]
(3
)
141
(
)
∑
(
)
{
}
∑
[
(
)
(
)
]
(4
)
142
143
where
is the group of adjacent elements to the

tetrahedron
and
represents the tensor
144
product
. T
he mass matrix
,
, the stiffness matrices
,
,
for all
{
}
,
the flux matrices
,
145
an
d
,
and
the auxiliary flux matrices,
and
,
are defined
by
Etienne et al. (2010). The
146
size of these matrices depends on the order
of the
polynomial basis (
e.g.
P0, P1
,
P2
, …, Pk
) used
147
9
for the nodal interpolation.
Staggered t
ime integration
is done through a
second order explicit
148
leap

frog scheme
,
which allow
s
the alternation
of
velocities
and stress
es
computation.
149
150
One of the main features of this scheme is the h

adaptivity
,
which allow
s
work
ing
with
151
unstructured tetrahedral meshes and
thus to adapt the mesh geometry
to
both
the physical
152
properties of the medium
and the problem geometry
(e.g. mesh refinement)
.
Furthermore,
the p

153
adaptivity
is possible
thanks to
the fluxes
between
neighbor
ing elem
ents, which
are such that
154
two
adjacent tetrahedra
may have
different interpolation orders.
The f
luxes
are computed via a
155
non

dissipative
centered scheme
that allows
choos
ing
different
degrees of freedom (DOF)
in
156
every tetrahedron.
As shown by Etienne et al. (2010), this
property is a powerful tool to optimize
157
the domain discretization
by adapting the element order to the medium
waves speed
.
Our
158
scheme includes
finite volume
approximation orders,
P0
(i.e. constant functions)
, linear
159
interpolation functions
,
P1
,
and quadratic interpolation functions
,
P2.
Etienne et al. (2010) have
160
shown that
this scheme is
efficient and
accurate enough
for
modeling
wave propagation
in
large
161
and
highly heterogeneous
elastic
media
.
The distribution of
in
terpolation
order
s
in the
162
computational domain is such that
P0 elements
describe
the CPML
slab while
elements with
163
both P1 and P2 approximations di
s
cretize the physical domain depending on their sizes. This
164
enhances the accuracy of the scheme and minimize
s
the computational load.
Current
165
computational developments will allow us in the near future to consider higher interpolation
166
orders away from the rupture zone for large

distance wave propagation.
167
168
3.
Dynamic rupture model
169
170
10
Earthquakes are highly nonlinear phenomena produced by sliding instabilities along geological
171
faults. The stability of the rupture system depends on several physical factors, like the initial
172
state of stresses in the
earth
, the material properties
and r
h
e
o
logy, the sliding rate,
the fault
173
geometry, and
the constitutive relationship governing the
mechanics of the rupture surface.
174
During rupture propagation, fault tractions evolve dynamically depending on all these factors
,
175
and the accuracy of an earthquake m
odel strongly depends on the correctness of boundary
176
conditions applied
to these tractions
in accordance with
the
fault
constitutive relationship (i.e.
177
friction law)
, which in turn depends on the accurate energy transportation through elastic waves
178
in the
medium.
Insuring the accuracy
of
both
boundary conditions (local feature) and wave
179
propagation (global feature)
has longly been a difficult task for many seismologists
.
In the
next
180
sections we will introduce
both
the fault boundary conditions and the frict
ion law
used
in our
181
dynamic source model
and
then its
formulation
in
to
the
DG
scheme
.
In spite of the attention they
182
deserve,
we shall
not discussed
here t
he features of wave propagation
since they
have been
183
previously
analyzed
by Etienne et al (2010)
and Tago et al. (2010)
.
184
185
3.1
Boundary conditions and friction law
186
187
The
fault
,
,
is
a
piecewise discretized
surface
with
a
normal vector
,
,
defined
on each
surface
188
element
.
Slip
and stresses over
are related through a friction law
so that
the
fault
tangential
189
tractions evolve
according to
that
law
,
which in turn
depends, for instance, on some fault
190
kinematic parameter
and wave propagation around the rupture tip
.
O
n the other hand,
the strain
191
field
is
accommodated
,
in
the
elastic medium, through displacement
(i.e. velocity)
discontinuit
ies
192
11
across
.
It is thus convenient to
split our domain into
a
plus

and
a
minus

side
with respect to
193
the fault
(Figure 1)
and
express
the
limiting
velocity field
,
,
over
as
194
195
(
)
(
(
)
)
196
We
define the vector
as
the velocity
discontinuity
across
,
where
its fault tangential
197
component, namely the fault slip rate, is
198
199
⟦
⟧
,
(
5
)
200
201
and
its fault normal component is
202
203
⟦
⟧
.
(
6
)
204
205
T
he slip magnitude at
any
time
is thus defined as the integral of
the modulus of
over
time
:
206
207
(
)
∫
‖
(
)
‖
.
(
7
)
208
209
The
dynamic rupture
boundary conditions
on
the
fault
are the
following two
jump conditions
,
210
involving the tangential fields
(Day et al., 2005)
,
211
212
‖
‖
(
8
)
213
‖
‖
,
(
9
)
214
12
215
and a third jump condition
applied
to (
6
)
216
217
,
(
10
)
218
219
where
is the
fault
tangential traction vector
,
and
is the
fault
frictional
strength.
The fault
220
strength
is
determined
by the Coulomb friction law, which depends on the
fault normal stress
,
221
,
the friction coefficient
,
,
and the fault cohesion,
,
as
222
223
.
(
11
)
224
225
C
ondition (
8
)
upper bounds
the magnitude of
to the
fault
strength
,
,
that always acts
226
opposite to sliding on
. The second condition (
9
)
forces
the slip rate
to be
parallel to
,
the
227
tangential traction.
Another implication of
these
jump conditions is that, whenever the
inequality
228
of (
8
) holds, the
slip rate vector is zero, which
can be
easily
seen through the modulus of (
9
).
229
Finally, condition (
10
)
implies
that there is neither fault opening nor
mass interpenetration
230
across the fault
during rupture propagation
.
231
232
In
nature, t
he friction coefficient
depends on the
fault
slip, slip
rate
and state variables
233
accounting for
the
sliding history and fault age
(
e.g.
Dieterich, 197
9
;
Ruina, 1983
)
.
However,
we
234
shall assume
a simple slip weakening law
(Ida, 1972
; Palmer and Rice
, 1973
),
which make
s
235
linearly
depend on the
slip as
236
237
13
(
)
(
)
(
)
(
)
,
(
12
)
238
239
where
and
are the static and dynamic friction coefficients
, respectively,
and
is the
240
critical slip weakening distance.
241
242
A series of studies have tried to estimate
from historical earthquakes based on indirect source
243
observations through the inversion of strong motion seismograms (e.g. Ide and Takeo, 1997;
244
Mikumo and Yagi, 2003). However, due to the limited bandwidth of
the seismograms,
the slip
245
weakening distance
was poorly resolved in these studies. Moreover, dynamic models based on
246
such indirectly
inferred
values may be biased and not able to resolve the stress breakdown
247
process over the fault (Guatteri and Spudich, 2000; Spudich and Guatteri, 2004).
Direct
248
observation of
from near

field data is
seldom
possible
(
Cruz

Atienza et al., 2009
);
except for
249
so
me isolated cases where rupture propagated with supershear speeds (Cruz

Atienza and Olsen,
250
2010).
251
252
The manner
by which
we incorporate
t
he
fault
model
given
by
both
the
jump conditions (
8
)
to
253
(
10
)
,
and the
friction law (
11
) and (
12
),
into our DG scheme is
presented next
.
254
255
3.2 Discrete
source
model
256
257
The domain decomposition
introduced
by
(
2
)
should account
for the presence of
,
so that
the
258
fault surface
is
discretized with triangles
lying
on
the face
s
of
adjacent
tetrahedra
(Figure 1)
, i.e.
259
14
we
preclude
to be
embedded in
any tetrahedron
.
The physical domain,
, is then
260
decomposed
as follow
261
262
∑
{
(
)
}
263
264
w
here each
is a straight

sided tetrahedron with surface
and the union of all
elements
265
describe
s
a
geometrically conforming
mesh. The order
of the
polynomial basis
chosen
in our
266
method corresponds to P2
quadratic functions because
higher
approximation
orders
does not
267
significantly improve neither the a
ccurac
y of dynamic

rupture
numerical schemes
nor the
ir
268
convergence rate (
e.g.
Moczo et al., 2007;
Rojas et al., 2009
; Pelties et al., 2012
).
As we shall see,
269
k
eeping
a
low
approximation
order
(i.e. P2 interpolation functions)
provide
s
both
good
accuracy
270
and efficiency to
our numerical scheme.
271
272
Since
every
tetrahedron
has
its
own nodes
in the DG method
(i.e.
ten
independent
nodes
for P2
273
elements
, Figure 1
)
,
a
fault node is
then
composed
of
two co

located nodes
(i.e.
a split

node
)
.
O
ne
274
o
f them belong
s
to the

tetrahedron
within
the plus

side
of the domain
,
and the other to the
275
adjacent

tetrahedron in
the
minus

side
(
see Figure 1
)
.
This
means that
each
split

node
lie
s
276
between
two tetrahedr
a
sh
aring
a
fault element
.
Furthermore, since the DG method does not ask
277
for
field
continuity
over the element faces
,
the
dislocation produced by the rupture may be
278
handled
naturally
over the
fault
,
through the
discontinuit
y
of
the tangential velocities
279
{
}
. However, this should be
treated carefully
because
most of the elastic fields
must
280
remain co
nt
i
nuous across
.
281
282
15
(Figure 1)
283
284
S
ystem
(
3
)
and
(
4
)
is solved
everywhere inside
except over
,
where jump conditions (
8
)
to
285
(
10
) must be verified.
However, as pointed out by Benjemaa et al. (2009),
before treating the
286
fault fluxes
accordingly
we
should notice
that the system is not symmetric because
of
287
288
(
)
,
289
290
which
is
due to the variable transformation
(chg_var)
required
to group all the medium
291
properties on the left

hand side.
Al
though
our
method
do
es
n
o
t require the system
to be
292
symmetric
,
that condition is
convenient because
this way
our
scheme for P0 elements
essentially
293
reduces
to
the one
proposed by Benjemaa et al. (i.e. finite volume
approach
)
,
which was derived
294
from energy balance consideration across the fault
.
295
296
To obtain a symmetric system, we
first
multiply
(
4
)
by the symmetric
al
positive definite matrix
297
,
defined in Benjemaa et al
. (2009),
and
then isolate
the updated field values
to
get the
298
equivalent symmetric system
299
300
(
)
[
∑
(
)
{
}
∑
[
(
)
(
)
]
]
(
13
)
301
16
(
)
[
∑
(
)
{
}
∑
[
(
)
(
)
]
(
14
)
302
303
where
,
and
.
304
305
Since the flux
es
across
the fault
elements
{
}
must satisfy
the jump
306
conditions (
8
)
to
(
10
)
,
we
cannot
simply
use the centered scheme
proposed
by
Etienne et al.
307
(
2010
)
for
fluxes between regular elements
.
Introducing
the
fault vector
fluxes
and
for the
308
velocity and stress schemes, respectively, the system (
13
) and (
14
) m
a
y be rewritten as:
309
310
(
)
[
∑
(
)
{
}
∑
[
(
)
(
)
]
]
(
)
(
15
)
311
(
)
[
∑
(
)
{
}
∑
[
(
)
(
)
]
(
)
(
16
)
312
313
17
where
is a Kronecker delta which is 1 if
and
0
otherwise.
314
315
Since the fault boundary conditions (
8
) and (
9
) must be applied to the traction vector
,
316
following Benjemaa et al. (2009
)
we
notice that the flux
in the
velocity scheme (
13
)
th
r
ough any
317
element surface
may be expressed in term
s
of
tractions
as
318
319
[
(
)
(
)
]
[
(
)
(
)
]
.
320
321
By i
mposing
continuity of
over the fault
element
(i.e.
)
and
assuming the
322
same approximation order in the two
tetrahedra sh
aring
that
element
(i.e.
)
,
we set the
323
flux vector across the fault,
,
to
the unique traction vector
,
and
define the flux across the
324
fault as
325
326
(
)
(
)
.
(
17
)
327
328
B
y substituting equation (
17
) into
the velocity scheme (
15
)
and regrouping the terms
excluding
329
the flux across the fault
on
,
we get
330
331
(
)
(
)
,
(
18
)
332
333
where
is
334
335
18
(
)
[
∑
(
)
{
}
∑
[
(
)
(
)
]
]
336
Because the fault normal stress
determines the frictional
strength via (
11
) and boundary
337
conditions are applied to the shear tractions,
the fault traction vector
has to be decomposed
338
into its tangential
,
,
and normal
,
,
components
to
finally rewrite the
velocity scheme (
18
)
339
as
340
341
(
)
(
)
(
)
.
(
19
)
342
343
3.2.1
F
lux
es
across the fault
for the velocity scheme (
19
)
344
345
To
update velocities on
we need
to specify
the
flux across the
fault
(
17
)
in
(
19
)
.
To this
346
purpose
we
require
such that the jump conditions (
8
) and (
9
)
are fulfilled
,
and
so that
347
the
continuity of the
fault
normal velocity
is
also
guaranteed
(condition (
10
))
.
All
traction
348
conditions
are
verified
at every fault node and
for
every time step.
349
350
To compute the fault tangential tractions,
, we first notice that the inequality
of condition
(
8
)
351
along with
the modulus of (
9
)
impl
ies
that
352
353
,
(
20
)
354
355
19
which
means that
,
when
ever
‖
‖
remain
s
below
the frictional strength
,
,
the tangential
356
traction must
ensure
the continuity of the tangential velocities
at
every fault node shared by the
357

and

tetrahedral.
358
359
Let
us assume
that
so
in (
19
)
,
where
the

and

tetrahedral
lie
at the
plus

side
360
and the minus

side of the fault
,
respectively.
Since we only deal with fault nodes, we thus
361
construct the matrix
, which is the inverse mass matrix whose components depend
362
exclusively on those nodes. It is simply constructed from
by eliminating its rows and
363
columns associated with the off

fault nodes.
364
365
T
o compute
verifying
condition (
20
) we
require
both
tetrahedra
sharing a
fault element
to
366
have the same
order
, so that the nodes in both sides of the fault match to each other
(see Figure
367
1)
.
Besides this,
for the specific contribution of
,
we need to compute the volume and surface
368
integrals
of
and
, respectively,
in a standard element (Zienkiewicz et al.
,
2005) such that
369
370
and
(
21
)
371
372
wher
e
, the
i

tetrahedron volume, and
, the
i

tetrahedron fault surface, are the
373
corresponding Jacobians.
Then by substituting (
18
) into definition (
5
) and using (
21
), we
374
express the slip rate vector as
375
376
20
(
)
(
)
(
)
.
(
22
)
377
378
For getting
the tangential traction
,
we
finally
use (
20
) into (
22
)
,
which leaves us the following
379
expression
380
381
(
(
)
)
(
)
(
)
(
)
.
(
23
)
382
383
This
procedure
ensures
the continuity of the tangential velocity across the fault. However
,
if
the
384
time

dependent
frictional strength
,
,
is overcome by
the
modulus of the
tangential traction,
385
‖
‖
,
rupture
must
occur
and the tangential velocity
should
not longer
be
continuous across
386
the associated fault node
. In that case
,
the equality
of
condition (
8
)
hold
s
so that
387
388
‖
‖
‖
‖
.
389
390
Therefore, to compute the slip rate (
22
) at
every fault n
ode
and for every time step,
the
391
tangential traction
is
adjusted according to
the following
criterion
, which depends on whether
392
the fault point has broken or not
:
393
394
{
‖
‖
‖
‖
‖
‖
(
24
)
395
396
21
Since the nodes
within a
fault element
are coupled
through
the flux matrix
,
,
and the mass
397
matrix
,
,
when
rupture
happens in a given node and
the
condition
(
24
)
is imposed, tractions
398
in the remaining
nodes change
for
the same time step
.
Thus,
to accurately
and simultaneously
399
satisfy (
24
)
in
every fault node, i.e. to allow rupture propagation inside the fault elements,
we
400
use
an iterative
predictor

corrector
(PC)
scheme.
If
is the number of nodes in a given fault
401
element (i.e. six for P2 elements) and
is the number of nodes that have broken, then the
402
predictor

corrector
scheme operates only if
.
The PC
scheme will basically
403
find
the tangential
traction
s,
,
in the
unbroken
nodes,
given the boundary condition applied
404
in the
broken
nodes
for the same time step
. Thus,
the procedure will only influence the two
405
interacting elements sharing
the
same
breakable
portion of the fault surface.
406
407
Our
predictor

corrector
scheme is simple and converges fast
:
when a fault node
breaks
in a given
408
element
, the modulus of
its tangential traction
is set to
(condition (
24
))
.
Once a
pplied this
409
condition, the
pr
edicted
tangential traction
s
in
the unbroken nodes
of the same element
410
must be recomputed
accordingly
via
(
23
)
.
To this purpose, a new mass matrix
,
,
must
411
be constructed
considering
only the unbroken
fault nodes
while
for
the broken nodes
the
412
tangential
tr
action condition is set
. If the
magnitude of the new

predicted traction
s
overcomes
413
the fault strength
, then
it is
corrected by setting
it
to
.
This
updating
cycle continue
s
iteratively
414
through n
ew predictions and corrections
until no other
node breaks
inside
the element
after the
415
last correction.
The maximum number of possible iterations is given by the order of the
416
approximation, i.e. DOF, used in the fault elements
,
and will always be smaller than
(i.e.
417
five or less iterations
for P2 elements).
The PC procedure
has to be performed locally
in
each
418
piece of fault surface
. It is
an efficient iterative verification of boundary conditions and
419
22
represents
an ad
ditional reason to preserve low interpolation order
s
(
i.e.
P2 specifically)
in the
420
tetrahedra
sharing a fault segment.
421
422
Since
the
DG scheme
s
do not enforce
continuity of
the fields between two adjacent tetrahedra
,
423
and the
accuracy
of our
method
depends on
the special treatment of
velocit
ies
over the fault
, we
424
must take care of the fault normal tractions as well in the same manner we
do
for the tangential
425
components
.
We now deduce a formula to compute the normal traction,
, which
ensures
426
continuity of the fault normal velocity field
.
This
model constrain
is
given by the jump condition
427
(
10
)
.
428
429
As done for the tangential slip rate (equation
(
22
)
),
definition
(
6
)
may be developed
as
430
431
432
(
)
(
)
(
)
.
433
434
Using condition (
10
)
to force continuity of the
normal velocity,
we then define the
fault
normal
435
traction
436
437
(
(
)
)
(
)
(
)
(
)
,
(
25
)
438
439
from which the fault normal stress is given by
440
23
441
.
442
443
T
he frictional strength
,
,
can now be
computed in every fault node using (
11
)
as a function of
444
both
and the
friction coefficient,
, which
depends
on the fault
slip
,
(equation (
7
)),
445
through the slip weakening law (
12
).
446
447
D
efinitions
given for
the
normal (equation (
25
)) and tangential (equation (
23
))
traction
448
components
finally
allow us to update the velocity field
in
every
fault
node
via
equation
(
19
)
.
449
450
3.2.2
Fluxes across the fault
for the stress scheme (
16
)
451
452
For the stress scheme
and within the
i

tetrahedron
,
the flux
across the fault,
,
is
computed
453
using only the velocity field
in that
tetrahedron
, so that the
flux is given by
454
455
.
(
27
)
456
457
Th
e simplicity of this flux
comes from
the computation of a unique traction vector
on
the fault
458
that guarantees either the
continuity or
discontinuity of the velocity field depending on whether
459
the fault has broken or not
.
T
his
fault

flux
approximation is equivalent to the
one
proposed by
460
Benjemaa et al. (2009)
,
where the flux
estimation is
based on
an energy balance
consideration
461
across the fault
.
The simplicity of d
efinition (
27
)
is
due to the continuity of the
fault normal
462
velocity
implicit
in the computation of
through
(
25
).
463
24
464
4.
Rupture m
odel verification
,
convergence
and efficiency
465
466
Verification of dynamic rupture models
is a particularly difficult task.
Since no analytical solution
467
exist
s
for the spontaneous rupture
problem
(i.e. closed form equations for the resulting
468
motions)
, the only possible way to be confident
of a given approach i
s
the comparison of
results
469
for a well

posed rupture problem
between various
numerical
techniques based on different
470
approximations.
This kind of exercise has been systematically
performed
in the last years
by an
471
international group of
modelers
(Harris et al., 2009). In this section we present results for two
472
benchmark tests proposed by this group
,
i.e. TPV3 and TPV10
(
see:
473
http://scecdata.usc.edu/cvws/index.html)
, and compare them with those
obtained with
finite
474
difference, finite elem
ent
,
spectral element
, discontinuous
Galerkin and
spectral
boundary
475
integral
methods.
Based on
these
comparisons
,
we assess the numerical convergence
rate
of
our
476
method
, its efficiency
and determine numerical criteria to guarantee its accuracy
.
477
478
4.1 T
he Problem
V
ersion
3
(TPV3)
479
480
Consider the
spontaneous rupture
of
a vertical
right

lateral
strike

slip fault
embedded
in a
481
homogeneous fullspace
with
P

and S

waves speeds of 6000
m/s
and 3
434 km/s
, respectively,
482
and
density of 2670 kg/m
3
.
The fault is rectangular and measures 30 km long by 15 km
wide
483
(
Figure 2
)
.
Rupture n
ucleation
happens
in a
3 km by 3 km
square region, centered
both along

484
strike and along

dip, because the
initial shear stress
there is higher than the fault strength
.
The
485
friction law is linear slip

weakening
with zero cohesion
(Equations (
11
) and (
12
))
, and both
486
25
static and dynamic friction coefficients are constant over the fault.
The initial fault normal
487
traction is also constant
,
as
are
the static and dynamic fault stren
gths. Values for all source
488
parameters
,
i.e.
initial stress conditions and friction parameters,
are shown in Table 1.
Results
489
for this problem are compared with those
yielded by
the DFM finite difference scheme (Day et
490
al., 2005)
, the ADER

DG
discontinuous Galerkin scheme (Pelties et al., 2012)
and
the spectral
491
boundary integral equation method
by Geubelle and Rice
(
1995) with the implementation
of
492
E.M. Dunham (MDSBI: Multidimensional
Spectral Boundary I
ntegral, version 3.9.10, 2008
)
; all of
493
th
em
for a
n equivalent
grid size of 50 m.
494
495
All DGCrack solutions presented in this section
were
calculated in the same
100
x
110
x
95
km
3
496
physical domain discretized with unstructured h

adaptive meshes so that the element
497
characteristic lengths go from 1 km in the CPML slab to the desired length over the fault plane
498
(i.e. 1.0, 0.8, 0.5, 0.4, 0.3, 0.2 and 0.1 km).
499
500
(Figure
2
)
501
502
Our
first comp
arison corresponds to the rupture times
on
the fault plane
with the DFM method
503
(Figure
3
).
W
e have used
a
characteristic elements size of 100 m
over the fault (i.e. an effective
504
grid size
(internode distance)
of
about
50 m
in
our P2 elements approximation)
.
The fit
between
505
both solutions is
almost perfect
.
No significant difference may be seen in this comparison.
Figure
506
4
compare DGCrack
seismograms at fault points PI
(pure in

plane deformation)
and PA
(pure
507
anti

plane deformation)
(
see
Figure 2
) for the
slip
(
4
a)
,
shear traction
(
4
b)
and
slip rate
(
4
c and
508
4
d)
fields
with those yielded by the DFM, ADER

DG and MDSBI methods.
Except for weak
509
26
oscillations, the comparison is
also
excellent.
Besides the stress build

up
s
, which are
nicely
510
resolved at both observational points before failure,
let us
n
otice how the friction law is
well
511
resolved
as compared
to
the
other solutions
,
with
stress overshoots around 7
s
and 8 s at PI
,
and
512
8.5
s
and 10.5 s at PA. The associated slip reactivati
on
s
are
also
well
modeled
and can be seen in
513
the slip rate functions at both points.
Stopping phases from the fault edges strongly determine
514
the slip rate and are
sharply resolved
at 6.5 s and 7 s at PI, and at 4 s and 7.5 s at PA.
A closer
515
comparison of s
lip rates
at
both observational points
(Figure
4
d)
suggests
that
the closest two
516
solutions
to each other
correspond to
ours
and t
he one
generated
by the spectral boundary
517
integral method (MDSBI).
Despite weak oscillations present in the DFM, MDSBI and DGCrack
518
waveforms,
both
amplitude and phase
of the DGCrack and MDSBI solutions are
remarkably
519
similar.
520
521
(
Figure 3
)
522
523
(
Figure
4
)
524
525
Since no analytic solution exists for this problem, quantitative compar
isons between all the
526
approximat
ions
may give insights about their correctness.
Figure 5
b
presents a quantitative
527
comparison
of
all numerical solutions based on cross

correlation
(cc)
measurements of the slip
528
rate time series on both PI and PA observationa
l points.
Each colored square
of the Cross

529
Correlation Matrix (CCM)
corresponds to a cc

based metric between two solutions:
for a given
530
method, its metrics with respect to the other approaches are
those
corresponding to
its
531
associated
raw and colon
of the matrix
.
Values in the upper triangular part of CCM are
given by
532
27
533
,
534
535
where
is the
maximum
cross

correlation coefficient between the
i
and
j
solutions, and
536
subscript
min
read
s
for the
smallest
coefficient
of all
possible
combinations.
Values in the lower
537
triangular par
t
of CCM are given by
538
539
,
540
541
where
is
the delay in seconds between the
i
and
j
solutions for the maximum correlation
542
coefficient, and the
max
subscript means the maximum delay of all possible combinations.
Both
543
measures provide a quantitative mean to assess the similarity between solutions relative to the
544
worst comparison found betw
een all combinations. So they do not provide absolute cross

545
correlation information except for the auto

correlations along the CCM diagonal.
While the
546
ADER

DG solution is the closest to both the DGCrack and MDSBI solutions
in terms of correlation
547
coefficie
nts (see
in Figure 5b)
, the smallest phase error is found between the DGCrack
548
and MDSBI solutions
(see
)
.
As
it
may also be seen on Figure 4d
, the solution with the
549
lowest correlation metrics with respect to the
other ones
is that from DFM (i.e. first column and
550
first row). By averaging both metrics per solutions couple
,
we could better
see
which time series
551
are the closest to each other. Figure
6
presents
the results of this exercise
, where the two
552
d
iscontinuous Galerkin so
lutions reveal to be the more similar, although very close to the one
553
yielded by the Boundary Integral method.
Measures provided by this method should be
554
28
interpreted carefully since they do not account
for
the computational co
st required by each
555
method to achieve
its solution.
556
557
(Figure
5
)
558
559
One important issue in dynamic rupture model
ing
is
the
control
of
spurious oscillations
560
produced by the advance of the crack tip throughout the
discrete
lattice.
Using
either
artificial
561
viscosity
or intrinsic dissipation procedures
in rupture models
is a delicate matter because the
562
associated damping does not distinguish between numerical and physical
contributions
. In other
563
words, if bad
ly
handled,
dissipation
may
absorb
frequencies belonging to the physical

problem
564
bandwidth affecting, for instance, peak slip rates and rupture speed
s
, as shown by Knopoff and
565
Ni (2001).
This is probably why the ADER

DG
scheme (de la Puente, 2009), which uses
566
intrinsically dissipativ
e Go
dun
ov fluxes (
González

Casanova, 2006; LeVeque, 2002
), requires
567
high
interpolation
orders
to achieve good accuracy
(e.g. compare O2 with O3 or higher order
568
solutions in Pelties et al., 2012).
In the DGCrack scheme, in contrast, the centered flux scheme we
569
use
is conservative so that the energy is not intrinsically dissipated (
Hesthaven and Warburton,
570
2008) but, because it is dispersive, our method
presents spurious oscillations.
However, since
571
these oscillations remain reasonable small and we expect them to be
even smaller for physically
572
attenuating media or different friction laws (e.g. rate

and state

dependent),
we have decided not
573
to integrate an artificial viscosity.
574
575
(
Figure
6
)
576
577
29
To assess the convergence rate of the DGCrack scheme and to determine a
quantitative criterion
578
that guaranties its accuracy, Figure
s
5a and 5
c present
two
different
error
metrics
, defined by
579
Day et al. (2005)
as the relative
root mean
square difference between a given solution and the
580
reference one (i.e. the DFM solution for h
= 50 m)
,
as a function of the
characteristic
elements
581
size on the fault.
These
metric
s
corr
espond to the rupture times
over the fault
and the peak slip
582
rate
s at fault points PA and PI (Figure 2).
Both figures
reveal
a power convergence rate of the
583
DGCrack
method
with
regression
exponents
reported in Table 2 and
compared
to
th
ose for other
584
methods. We also report in that table the
exponent
for the final slip on both observational points
585
(not shown in the figure).
It is important to notice that all
DGCrack
solutions
presented
in the
586
manuscript correspond to unstructured meshes with refinements around the rupture surfaces.
587
As mentioned by de
la Puente et al. (2009
), t
his is
a
critical
issue
since the
accuracy
of solutions
588
significantly
depends on the quality
of the tetrahedral lattice
built
independently
using
standard
589
tools
(
in this work we have used the
Gmsh
software developed by
Geuzaine and Remacle, 2009)
,
590
as can been seen in the
error
dispersion
on
both Figures 5a and 5c
with respect to the regression
591
lin
es
.
592
593
The accuracy of
dynamic rupture models depends on the resolution of the cohesive zone, which
594
may vary during rupture evolution.
This zone is defined as
the fault region where the stress
595
breakdown process takes place. The amount of fault elements inside this region behind the
596
rupture front, Nc,
is thus the numerical criterion that guaranties a given accuracy level. F
igures
597
5a and 5c
also reveal that one
or
more
fault element
s
inside the cohesive
zone (
i.e.
)
is
598
enough to achieve
errors
smaller
than 1 % and 10 % for rupture times and peak slip rates,
599
30
respectively.
In the TPV3 test case, this condition implies fault element sizes smaller or equal to
600
~45
0 m.
601
602
We finally address
a fundamental question in computational sciences
:
the numerical
efficiency.
603
Pelties et al. (2012) have recently introduced a method to solve the
dynamic rupture problem
604
based o
n a discontinuous Galerkin scheme
that incorporates
a sophisticated strategy allowing
605
arbitrarily high order approximations
in space and time
, i.e. the ADER

DG method. Since both the
606
DGCrack and ADER

DG methods share many different
capabilities
linked to the DG
607
approximation
,
it is worthy to
see difference
s and to estimate
the computational cost of each
608
method to achieve
the same
accuracy level.
Since the ADER

DG method is based on a modal
609
approximation,
the fluxes across the element faces are sensitive to
all
modes across
the
element
;
610
while in a nodal appr
oximation
,
like in the DGCrack method, they only depend on the nodes lying
611
on the element face where the flux is computed. This difference
translates
into fewer
612
computations in the
nodal
approximations (Hesthaven and Warburton, 2008).
Figure 5d
613
presents a
quantitative comparison
of rupture times errors as a function of total computing
614
times (
i.e. the
CPU ti
me
, which is given by the duration of each simulation multiplied by the
615
number of
cores
) for both methods. The DGCrack simulations were run in Pohualli,
the
parallel
616
computer with 1
72
cores
(i.e. 2.33 GHz quad

core Xeon processors)
of the Department of
617
Seismology at UNAM
.
Since the computational cost is not reduced by the ADER

DG method
618
when using high approximation orders (Pelties et al., 2012), values r
eported in Figure 5d
619
correspond to the ADER

DG O3 solution. Both regression lines
have about the same slope;
so
620
similar
CPU time difference
s
between the methods
are
expected for any grid size.
If we take
the
621
same accuracy level we used to establish the
Nc
condition
in the
last paragraph (i.e.
1 %
error
for
622
31
rupture times)
then the CPU time of the ADER

DG O3 method i
s
about 30 times larger than the
623
one
required
by the DGCrack method
in our computing platform. Since
the
simulation
times
624
reported by Pelties et
al. were obtained using a 0.85 GHz
BlueGene
parallel computer
, the actual
625
CPU time difference would be reduced to a factor
of
~10 if both methods were run in the same
626
platform,
which
still
is a
significant
factor
,
especially if multiple large scale simula
tions
are
627
required for an inversion modeling purpose, for example
.
This comparison was made with the
628
available data reported by Pelties et al. (2012), however it would be interesting to design a
629
dynamic

rupture speed test that could better reflect the CPU
time required by both schemes.
630
Besides, it is important to notice that ADER

DG show
s
smooth time series
on the fault
and
631
accurate wave propagation away from the fault (
Dumbser and Kaser, 2006
). Theses
632
characteristics are not quantified here and may be
essential for long

range wave propagation
633
problems, for instance.
634
635
4.2
The Problem Version 10 (
TPV10
)
636
637
Our last verification test consist
s
of
a 60
o
dipping
normal fault reaching the free surface of a
638
homogeneous halfspace
(
Figure 7
) with P

and S

waves spe
eds of 5716 m/s and 3300 km/s,
639
respectively, and density of 2700 kg/m
3
.
The fault has the same dimensions
as
in TPV3 (i.e. 30
640
km long by 15 km width) but the
center of
its
3
x
3 km
2
nucleation patch
is
located deeper
,
at
12
641
km along dip.
Frictional and initial stress conditions
on the fault plane
are reported in Table 3,
642
where the cohesion term of equation
(
11
) is not zero, and
both
pre

stress conditions are
643
dependent on the along

dip distance, h
d
.
The unstructured h

adaptive tetrahedral m
esh used to
644
32
obtain the DGCrack solution in shown in
Figure 7
, which has a characteristic element size
of 1
00
645
m
over the fault
plane
.
646
647
(Figure
7
)
648
649
This test case has interesting features that are essential to verify a dynamic rupture model for
650
non

planar
faults. Since we deal with a dipping
normal fault reaching the E
arth
’s
surface,
651
reflected
waves are
bounced
back
to
the
source
,
inducing
transient
variations
of
the
fault
652
traction
vector
that
significantly affect rupture propagation via
the Coulomb failure
criterion
653
(e
quation
(
11
)
) (
Oglesby et al.,
1998
).
Figure 8
presents a comparison of rupture times over the
654
fault obtained with three different methods (Harris et al., 2009): a finite element (FaulMod;
655
Barall, 2010), a spectral element (SPECFEM3D; Kaneko e
t al., 2008) and the DGCrack methods.
656
The three solutions were computed with a fault elements size of 100 m.
Despite the complexity
657
of the rupture model, the match between all solutions is remarkably good. Dynamic effects on
658
rupture propagation due to the
presence of the free surface are clearly seen in
the
upper most 2
659
km, where
a secondary rupture front propagating down

dip is initiated about 2.5 s after
660
nucleation.
661
662
(Figure
8
)
663
664
Figure 9 shows on

and off

fault waveforms
yielded by the same three methods
Σχόλια 0
Συνδεθείτε για να κοινοποιήσετε σχόλιο