3D hp-Adaptative Discontinuous Galerkin Method

moanafternoonInternet και Εφαρμογές Web

11 Δεκ 2013 (πριν από 3 χρόνια και 10 μήνες)

190 εμφανίσεις



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