Pure and Applied Geophysics manuscript No.
(will be inserted by the editor)
Parallel 3D Simulation of a Fault Gouge Using
The Lattice Solid Model
Shane Latham1,2 , Steffen Abe1,2 and Peter Mora2,1
1
2
ACcESS MNRF
The University of Queensland
Brisbane, QLD, Australia, 4072
e-mail: [slatham,steffen,morap]@access.uq.edu.au
Earth Systems Science Computational Centre
The University of Queensland
Brisbane, QLD, Australia, 4072
e-mail: [slatham,steffen,morap]@esscc.uq.edu.au
Received: / Revised version:
Abstract Despite the insight gained from 2D particle models, and given
that the dynamics of crustal faults occur in 3D space, the question remains,
how do the 3D fault gouge dynamics differ from those in 2D? Traditionally,
2D modeling has been preferred over 3D simulations because of the computational cost of solving 3D problems. However, modern high performance
computing architectures, combined with a parallel implementation of the Lattice Solid Model (LSM), provide the opportunity to explore 3D fault micromechanics and to progress understanding of effective constitutive relations of
fault gouge layers. In this paper, macroscopic friction values from 2D and 3D
LSM simulations, performed on an SGI Altix 3700 super-cluster, are compared. Two rectangular elastic blocks of bonded particles, with a rough fault
plane and separated by a region of randomly sized non-bonded gouge particles, are sheared in opposite directions by normally-loaded driving plates.
The results demonstrate that the gouge particles in the 3D models undergo
significant out-of-plane motion during shear. The 3D models also exhibit a
higher mean macroscopic friction than the 2D models for varying values of
interparticle friction. 2D LSM gouge models have previously been shown to
exhibit accelerating energy release in simulated earthquake cycles, supporting
the Critical Point hypothesis. The 3D models are shown to also display accelerating energy release and good fits of power law time-to-failure functions
to the cumulative energy release are obtained.
Key words lattice solid model – discrete element method – parallel simulation – granular shear – macroscopic friction – accelerating energy release
2
Shane Latham, Steffen Abe and Peter Mora
1 Introduction
Many naturally occurring faults and shear zones contain regions of granular
material. The presence of gouge particles plays a fundamental role in influencing the macroscopic behavior of these systems. In order to gain a greater
understanding of earthquake nucleation processes in fault gouge zones, it
is important to characterize the relationships between the microscopic and
macroscopic mechanics. Computational simulation has played an important
role in the analysis of complex granular materials, allowing researchers to
vary micro-mechanical parameters and observe the resulting influence on the
macro-mechanics of the computational model.
In particular, Distinct Element Methods (DEMs) [5, 13, 15] have successfully been employed to simulate the dynamics of fault shear regions [7, 14–20,
22]. To date, the majority of results have been produced from 2D simulations,
despite the fact that real fault zones are 3D entities. With the continual performance improvement of computing hardware, comes the opportunity to
study computationally demanding 3D models of complex granular systems.
Of specific interest is the influence of the out-of-plane particle dynamics on
the macroscopic behavior of 3D systems.
The Lattice Solid Model (LSM) [13,21] is a variant of the DEM which
has been used to model the dynamics of fault gouge processes. LSM simulations of 2D faults have yielded results which provide an explanation of gouge
weakness and the heat flow paradox [14,15] (HFP). Simulations with other
2D DEMs, using round particles, also demonstrate low coefficients of friction
due to rolling [18, 20]. Typical 2D fault gouge models, using the LSM, have
involved tens of thousands of particles. For comparable 3D problems, particle
numbers can readily increase into the millions. These large 3D problems have
remained intractable for serial implementations of the LSM. Parallel computing architectures, such as the SGI Altix 3700 super-cluster, provide the
opportunity to solve much larger computational problems than traditional
single processor systems. In order to take advantage of high performance systems, a Message Passing Interface version of the LSM has been implemented
[1]. Recent benchmarks demonstrated an 80% parallel efficiency for the parallel LSM on 128 processors of the SGI Altix 3700 [10]. These results, for
large 2D wave propagation problems, indicate the potential for the LSM to
simulate more computationally challenging 3D fault gouge dynamics.
This paper describes results from simulations of 2D and 3D gouge regions
using the parallel implementation of the LSM. Section 2 gives an overview
of the LSM and the contact laws imposed on interacting particles, which is
followed by a description of the fault gouge models in Section 3. A comparison
of the macroscopic friction (also termed the effective or fault friction) for a
2D model and varying thicknesses of 3D model are presented in Section 4.
Power-law time-to-failure fits for the cumulative energy release and an
evolution in size–frequency statistics in the lead up to large earthquake events
are evidence that the crust behaves as a Critical Point (CP) system [4, 9,
23]. If this is the case, then intermediate earthquake prediction is possible.
Parallel 3D Lattice Solid Fault Gouge Simulation
3
However, the circumstances under which fault systems behave as CP systems
is still an open question. 2D LSM fault gouge simulations have previously been
shown to exhibit accelerating energy release and stress correlation function
evolution in simulated earthquake cycles [16,17]. In Section 5 it is shown that
3D models also exhibit accelerating energy release and good fits of the power
law time-to-failure functions are obtained.
2 Lattice Solid Model
The LSM [13,21] is a particle based model similar to the DEM [5]. The model
consists of spherical particles which are characterized by their radius, mass,
position and velocity. The particles interact with their nearest neighbours
by imparting elastic and frictional forces. Section 2.1 describes the particle
interactions employed in the 2D and 3D models introduced in Section 3.
Section 2.2 briefly describes the parallel design of the LSM implementation.
2.1 Interactions
In the gouge simulations, particles were restricted to interact in one of two
ways. A pair of particles could be involved in either a bonded interaction
or in a frictional interaction with one another. A volume of bonded particles
simulates a linear elastic solid within the model. The force Fbi,j which particle
i exerts on a bonded neighbouring particle j is given by
di,j
,
(1)
Fbi,j = k b (|di,j | − (Ri + Rj ))
|di,j |
where k b is the spring constant, Ri is the radius of particle i and di,j = ri −rj
with ri being the centre position of the particle i. From this equation, it can
be seen that the bonded force is repulsive if |di,j | < (Ri + Rj ) and attractive
if |di,j | > (Ri + Rj ).
A pair of particles i and j which are not bonded and which come into
contact (di,j < Ri + Rj ) undergo a frictional interaction. In this case, the
force Ffi,j which particle i exerts on particle j can be expressed in the form
Ffi,j
=
Ffi,jr
(
Ffi,js ,
+
Ffi,jd ,
i and j in static friction contact,
i and j in dynamic friction contact
(2a)
where Ffi,jr is the elastic repulsive force i exerts on j, Ffi,js is the static frictional
force which i imparts on j and Ffi,jd is the dynamic or slipping frictional force
which i imparts on j. The expression for repulsive elastic force is identical in
form to (1):
di,j
Ffi,jr = k r (|di,j | − (Ri + Rj ))
,
(2b)
|di,j |
4
Shane Latham, Steffen Abe and Peter Mora
with k r being the spring constant. When particles i and j are in static contact
the shear force is given by
Ffi,js = −k s dsi,j eTi,j ,
(2c)
where k s is the shear spring constant, dsi,j is a circumferal shear displacement
and eTi,j a unit vector formed by projecting the the relative particle velocity
onto a plane with normal di,j . This tangent unit vector can be expressed as
eTi,j =
where
ti,j = ui,j −
ti,j
,
|ti,j |
ui,j · di,j
2
|di,j |
di,j ,
and ui,j = vi − vj is the velocity of particle j relative to particle i, with vi
the velocity of particle i. The dynamic friction force has the simple Coulomb
form
fd
Fi,j
= −µP Ffi,jr eTi,j ,
(2d)
where µP is the intrinsic or inter–particle friction coefficient.
The transition from static friction to dynamic friction occurs if
Ffi,js + k s ti,j ∆t > µP Ffi,jr .
where ∆t is a timestep size. The transition from dynamic to static friction
occurs if
Ffi,jd + k s ti,j ∆t ≤ µP Ffi,jr .
All forces are applied at the particle centre, consequently, per–particle
rotational dynamics are not modeled. With this constraint, particles slide past
one another, as opposed to rolling over one another. An artificial viscosity is
also present in the model to prevent the buildup of kinetic energy in the
closed system. The amount of viscous damping has been chosen such that
the rupture dynamics are not significantly influenced [13].
2.2 Parallel Implementation
The parallel implementation of the LSM follows a modified master-worker
process model with MPI (Message Passing Interface [12]) used as the interprocess communication layer. A master process provides a high level of control
and external communication I/O facilities. The worker or slave processes perform the particle interaction computations. In contrast to a pure master-slave
approach, direct communication between worker processes is used instead of
communication involving the master process whenever possible. This design
minimizes the computation performed by the master process and also greatly
reduces the amount of communication between master and worker processes.
Parallel 3D Lattice Solid Fault Gouge Simulation
5
Computation is shared among processes by assigning a subset of particles
to each worker. This distribution of particles is based on a regular-grid partitioning of the geometrical domain. In addition to the particles located in
a particular subregion, the data set assigned to each worker–processor also
contains all particles interacting with any particle in the subregion. Due to
the short spatial range of the interactions, the number of duplicated boundary particles is small in comparison to the total number of particles. The
duplication of particles creates a small computational overhead, as the forces
computed for an interaction involving a duplicated particle are computed
in each relavent worker process. This small amount of duplicated computation is justified because it reduces the cost of communication. More detailed
descriptions of the parallel implementation can be found in [1, 10].
This parallelisation strategy has proven effective for a range of problems.
In particular, benchmark 2D wave propagation models involving millions of
particles [1, 10] have been used to show good scaling behaviour on SMP machines such as the SGI Origin3800, Compaq Alphaserver SC40 and the SGI
Altix3700. Large 3D fault gouge models, similar to the models described in
Section 3, involving hundreds of thousands of particles have also demonstrated
good scaling on the SGI Altix3700.
3 Fault Gouge Model
In the simulations, a fault gouge is represented as two rectangular elastic
blocks of bonded particles, with a rough fault zone separated by a region of
randomly sized non-bonded gouge particles. The elastic blocks are sheared
in opposite directions by normally-loaded driving plates. Figure 1 illustrates
the particle model setup. The block particles are uniformly sized with radius
R0 = 500µm and bonded in a regular 3D lattice. The roughness particles
provide surface-roughness and range in radial size from 0.4R0 to R0 in the 3D
models and from 0.165R0 to R0 in the 2D models. Each roughness particle is
bonded to neighbouring block and roughness particles. The interaction forces
on bonded particles are computed using equation (1).
The gouge particles are non-bonded, and impart forces on neighbouring
gouge and roughness particles via the friction interaction described in equations (2). The gouge particles have the same radial size range as the roughness
particles. The boundary particles at the top and bottom of the elastic blocks,
are elastically bonded to walls (not shown in Figure 1) which are parallel to
the x-z plane and apply compressive forces in the y direction. These walls are
also sheared in opposite x directions at constant velocity, as indicated by the
arrows in the leftmost illustration of Figure 1. A circular boundary condition
occurs on y-z planes at the left (x=0) and right (x = 34R0 ) extents of the
particle domain (a particle exiting the right hand side of the model reappears
at the left side and vice–versa). In the 3D models, fixed frictionless elastic
confining walls in the x-y planes prevent particles from being “squeezed” in
the z direction out of the gouge region.
Simulation results from a single
6
Shane Latham, Steffen Abe and Peter Mora
Fig. 1 Example fault gouge model setup: (left) Full view of the Z8 model from the
positive z direction with shearing motion in the positive and negative x directions
and constant normal force applied in the negative y direction. Initial gouge-region
particle configurations for the Z0 model (top right) and Z8 model (bottom right).
Number of gouge
particles
Total number of
particles
No. of worker
processes
Runtime (hours) for
2000000 timesteps
2D
Z0
Z2
713
1169 2106 4123
Z4
3D
Z16
Z8
7994
Z24
Z32
12003 16195
1701 2881 5476 10690 20972 31441 41865
8
12
12
24
24
36
48
2
4
7
8
14
16
17
Table 1 Total number of particles, number of gouge particles, number of worker
processes and run–times for the macroscopic friction simulations in Section 4.
Parallel 3D Lattice Solid Fault Gouge Simulation
7
2D model geometry and multiple 3D model geometries are presented in the
following section. The 2D and 3D models have the same x dimension and similar y dimension, with the difference in y dimension due to the different circle
and sphere regular–lattice packings used in 2D and 3D, respectively. Table 1
gives the particle counts in each of the models, both total number of particles
and number of particles in the gouge region, along with number of processors
and approximate run-times for the macroscopic friction simulations. The 3D
models have identical x and y dimensions and differing z dimension sizes.
The 2D model is referred to as the Z0 model and the 3D models are referred
to as the Z2, Z4, Z8, etc models (2 regular z-layer packings, 4 regular z-layer
packings, 8 regular z-layer packings, etc). The same constant pressure is applied to all models by the uppermost wall (which moves in the y dimension
to maintain constant pressure) and the bottom wall remains fixed in the y
dimension. Table 2 summarises the parameters and properties of the gouge
models.
Model size
Width × Height
2D: Z0
3D: Z2-Z32
34R0 × 64.6R0
34R0 × 64.4R0
ZN model depth,
N = 2, 4, 8, ..., 32
R0
3N
√−2
3
+2
Gouge particle
radius Rj
0.165R0 ≤ Rj ≤ R0
0.4R0 ≤ Rj ≤ R0
Roughness particle
radius Rj
0.165R0 ≤ Rj ≤ R0
0.4R0 ≤ Rj ≤ R0
Block Young’s
modulus
4.33GPa
5.66GPa
Block Poisson ratio
0.25
Block particle
radius R0
500µm
Gouge region
height
≈ 10R0
Normal stress
15MPa
2600kg/m3
Density
Shear velocity
b
r
k =k =k
s
µP
Table 2 Model parameters.
0.001R0 /s, 0.0001R0 /s
735MN/m
0.2, 0.4, 0.6
8
Shane Latham, Steffen Abe and Peter Mora
4 Macroscopic Friction
The shear stress on the driving walls divided by the normal stress on these
walls is termed the instantaneous macroscopic friction and denoted µI . In the
macroscopic friction simulations, each of the driving walls was sheared at a
speed of 0.001R0 /s ≈ 0.001VP (where VP is the P-wave velocity). Figure 2
shows snapshots of the gouge layer in the Z0 and Z8 models at 0%, 200% and
400% shear strain. Some localisation can be seen in the Z0 model near the
top of the gouge region at 200% shear strain. By 400% shear strain the linear
bands indicate that the strain is evenly distributed throughout the layer.
Figures 3(a) and 3(b) plot the macroscopic friction and dilatency statistics, respectively, of the 2D and 3D models for three different values of the
intrinsic friction coefficient: µP = 0.2, µP = 0.4 and µP = 0.6. Here, the
dilatency h is simply measured as the y component of the top wall position
Fig. 2 Snapshots of the Z0 (left) and Z8 (right) gouge regions (for µP = 0.4).
Gouge particles are shaded in bands according to their original x coordinate. The
top row is the original configuration, middle row is 200% shear strain and bottom
row is 400% shear strain.
Parallel 3D Lattice Solid Fault Gouge Simulation
9
(the bottom wall remains in fixed position). When calculating the statistics,
the initial “load-up” period (10% strain) is ignored. Tables 3(a) and 3(b)
contain the data plotted in Figure 3. In the 2D and 3D models the mean,
mininum, maximum and standard deviation macroscopic friction values all
increased with increased intrinsic friction. The mean macroscopic friction values for the µP = 0.4 and µP = 0.6 cases, are unrealistically high. These large
values are due to the absence of rotational dynamics in the model. Similar
large macroscopic friction values were produced by [18, 20] in non-rotational
2D numerical simulations. Laboratory experiments with glass rods measure
a 2D macroscopic friction value of E[µI ] ≈ 0.3 and experiments with glass
beads measure a 3D macroscopic friction value of E[µI ] ≈ 0.45 [6, 11]. Numerical investigations with round rotational particles in 2D [8, 18, 20] and 3D [8]
produced results which agree closely with the glass rod/bead experimental
(a)
(b)
p
Fig. 3 (a) Macroscopic friction µI statistics: E[µp
V [µI ], min[µI ] and
I ], E[µI ] ±
max[µI ] and (b) Dilatency statistics: E[h], E[h] ± V [h], min[h] and max[h].
Shane Latham, Steffen Abe and Peter Mora
(a) Macroscopic friction statistics.
µP = 0.2
µP = 0.4
µP = 0.6
p
p
p
E[µI ] min[µI ] max[µI ] V [µI ] E[µI ] min[µI ] max[µI ] V [µI ] E[µI ] min[µI ] max[µI ] V [µI ]
2D Z0 0.49 0.22
0.66
0.05
0.82 0.31
1.02
0.09
1.26 0.67
1.62
0.14
Z2 0.55 0.39
0.69
0.04
0.89 0.54
1.12
0.10
1.43 0.65
1.95
0.18
Z4 0.57 0.43
0.66
0.03
0.94 0.67
1.11
0.07
1.37 0.85
1.71
0.15
3D Z8 0.59 0.51
0.64
0.02
1.00 0.78
1.11
0.04
1.46 0.93
1.67
0.09
Z16 0.59 0.51
0.64
0.01
0.96 0.75
1.04
0.03
1.42 0.93
1.60
0.08
Z24 0.59 0.52
0.62
0.01
0.96 0.75
1.05
0.03
1.42 0.94
1.61
0.08
Z32 0.59 0.52
0.62
0.01
0.98 0.75
1.05
0.03
1.44 0.94
1.59
0.07
(b) Normalised dilatency statistics.
µP = 0.2
µP = 0.6
q
q
V [ĥ] range[ĥ] min[ĥ] max[ĥ] V [ĥ] range[ĥ] min[ĥ] max[ĥ] V [ĥ]
2D Z0 0.43
-0.29
0.14
0.05
0.89
-0.63
0.25
0.14
1.31
-1.00
0.31
0.27
Z2 0.24
-0.14
0.10
0.04
0.77
-0.57
0.19
0.15
1.49
-1.04
0.45
0.37
Z4 0.17
3D Z8 0.15
-0.10
0.07
0.03
0.77
-0.56
0.21
0.14
1.47
-0.94
0.53
0.29
-0.11
0.05
0.03
0.69
-0.59
0.10
0.13
1.42
-1.10
0.32
0.34
Z16 0.14
-0.11
0.03
0.02
0.61
-0.49
0.12
0.12
1.09
-0.87
0.22
0.24
Z24 0.15
Z32 0.14
-0.12
-0.11
0.03
0.03
0.02
0.02
0.57
0.64
-0.48
-0.54
0.09
0.10
0.10
0.13
1.22
1.17
-0.98
-0.84
0.25
0.33
0.27
0.26
range[ĥ] min[ĥ] max[ĥ]
10
µP = 0.4
q
Table 3 Macroscopic friction and normalised dilatency statistics ĥ = h − E[h].
Parallel 3D Lattice Solid Fault Gouge Simulation
11
values. Laboratory experiments with angular sand [6,11] give a higher 3D
macroscopic friction value of E[µI ] ≈ 0.6. 2D numerical simulations [15,16,
22], which account for rotational dynamics through the use of aggregate angular particles, also give macroscopic friction values of E[µI ] ≈ 0.6. In particular, the numerical results from [15, 16, 22] supported that for a wide range of
interparticle friction values, the macroscopic friction remains at E[µI ] ≈ 0.6
when shearing is non-localised. For highly localised shear, the macroscopic
friction values can drop as low as E[µI ] ≈ 0.3.
The first notable difference between the 2D and 3D cases are the mean
values of the macroscopic friction - the 3D values are greater for each value
of intrinsic friction than the respective macroscopic friction values in the 2D
simulations. This has been observed in similar numerical experiments [8] and
concurs with laboratory rod and bead experiments where the macroscopic
friction was measured to have higher mean values in 3D than 2D [6,11]. The
minimum macroscopic friction value also differs significantly between the 2D
and 3D models where the 3D models show a smaller range of instantaneous
macroscopic friction variation. While the minimum macroscopic friction value
is greater in all the 3D models than the 2D model, the maximum macroscopic
friction values for the larger 3D models are comparable to the 2D models.
The standard deviation values, in Table 3(a), of the Z2 and Z4 models
are comparable with (and in some cases exceed) those of the 2D model. It
is expected that the variance in the macroscopic friction decreases as the
number of particles in the model increases and certainly this is the case for
the Z8–Z32 models. In the thin z–dimension Z2 and Z4 cases, the particles
are hampered in the z direction of motion due to the close proximity of the
frictionless confining walls. Hence, particle reorganisation in this dimension
is restricted, promoting gouge dilation and subsequent larger oscillations in
the instantaneous friction. This is supported by the larger dilation maximum
values for the Z2 and Z4 models in Table 3(b).
The effect of the frictionless confining walls in the thin 3D models becomes
apparent when the gouge-particle trajectories are examined. Figure 4 plots
gouge-particle displacement histograms for the Z0, Z2, Z4 and Z32 models in
each of the x, y and z spatial dimensions. The 2D model gouge deformation is,
obviously, accommodated in the x and y dimensions only, and subsequently
the Z0 model has the highest proportion of particles which move furthest in
these directions. In the 3D models, as the model size increases, there is a
reduction in the amount of displacement of gouge particles in the x and y
directions, and an increase in the displacement in the z direction. It is interesting to observe the dramatic difference between the Z0 and Z32 trajectory
histograms in the x (shear) direction. Despite relatively modest (but significant) motion in the z direction, where >90% of Z32 gouge particles travel
between 2R0 and 6R0 in the z direction, there is a significant change in the
trajectory distance in the shear direction. In the Z0 model, >99% of particles
travel between 20R0 and 40R0 , while in the Z32 model only 50% of gouge
particles travel between 20R0 and 40R0 with the remaining 50% travelling
between 10R0 and 20R0 in the x direction.
12
Shane Latham, Steffen Abe and Peter Mora
While the 2D vs 3D values of Table 3(a) agree qualitatively with experiment, laboratory mean macroscopic friction values for 3D glass bead experiments (0.45) exceed the 2D glass rod mean friction values (0.3) by more than
30% [6, 11]. The 3D mean values in Table 3(a) exceed the corresponding 2D
values by at most 22%. Again, this is likely to be an artifact of the absence
Fig. 4 Gouge particle displacement histograms in the x (top), y (middle) and z
(bottom) directions in the µP = 0.2 models.
Parallel 3D Lattice Solid Fault Gouge Simulation
13
of rotational dynamics within the models. Further numerical investigations,
which incorporate rotational effects, need to be conducted in order to resolve
the current quantitative differences between laboratory experiment and the
numerical simulation results.
5 Accelerating Energy Release
An increasing number of field observations suggest an evolutionary behavior
of earthquakes. These observations support the Critical Point (CP) hypothesis, which suggests that the largest earthquakes can only occur after the
system has reached a critical state (when stress correlations exist in the system at all length scales). A power law time-to-failure function has found to
fit well with the observed Benioff strain release in the lead-up to some large
events [3,9]. Therefore, it is of interest to determine whether this same power
law time-to-failure function fits cycles in simulated earthquakes. 2D LSM
gouge models have previously shown accelerating energy release and stress
correlation function evolution [16, 17].
The 3D models were also seen to exhibit the accelerating energy release
behaviour. Figure 5 shows two earthquake cycles, one from a Z4 model and
one from a Z8 model. In these models, the inter-particle friction was set as
µP = 0.6 and the shear velocity of the driving plates was slowed by a factor
of 10 to 0.0001VP . Each small circle in Figure 5 represents the occurrence of
simulated earthquake. As well as the cumulative energy release for the cycle,
this figure plots an associated fit of the power law function
c
Ec (t) = A + B(tf − t)
(3)
where A, B, c and tf are the fit parameters with tf representing the failure
time. The RMS error of the power law time-to-failure fit divided by the RMS
error of a linear fit provides a quality–measure for the time-to-failure fits [2].
The lower the quality-measure value, the more superior the power-law fit
when compared with the linear fit. The earthquake cycles of Figure 5, show
good fits of the power law function, with quality–measure values of 0.536 for
the Z4 cycle and 0.478 for the Z8 cycle.
These initial 3D results provide further support for the CP hypothesis.
Accelerating energy release is not present before all large earthquake events
(simulated or real). It seems possible, however, that the analysis of 3D numerical models will provide insight as to the conditions under which a fault
(or fault network) behaves as a CP system.
14
Shane Latham, Steffen Abe and Peter Mora
Fig. 5 Examples of accelerating energy release (Z4 model top and Z8 model bottom) with power law curve fit.
Parallel 3D Lattice Solid Fault Gouge Simulation
15
6 Conclusions
Using a parallel/MPI implementation of the LSM, simulations of a simplified
fault gouge produced higher mean macroscopic friction values in 3D models
than in 2D models. These simulations support the hypothesis that friction is
greater in 3D than in 2D due to the extra dimension of particle interaction.
For the thin 3D models, an end wall effect was observed which restricted
gouge particle motion in the z-dimension and resulted in higher maximum
macroscopic friction response than the thickest 3D models. Particle trajectories, in the 3D models, were shown to have a significant z displacement
(Figure 4), and this dramatically reduced particle displacement in the direction of shear. These results demonstrate the important influence of the third
dimension on both the microscopic and macroscopic dynamics of shear zones.
The 3D gouge simulations have also produced instances of accelerating
energy release in the lead-up to the larger simulated earthquake events. This
provides further support for the CP hypothesis.
The unrealisticly high values observed for the mean macroscopic friction
in both the 2D and 3D models is due to the lack of rotational dynamics in
the model. Future simulation studies will incorporate rotational dynamics
both directly (by assigning angular properties to individual particles) and
indirectly (through the use of bonded aggregates of particles in the gouge
region). By generating 3D models with bonded-aggregate grains, the effects
of grain angularity and grain disintegration on the macroscopic behaviour of
the 3D system can be studied and their importance to earthquake processes
evaluated.
ACKNOWLEDGMENTS
The authors would like to thank Julia Morgan for her review comments,
which improved the quality of the manuscript. Funding support for this work
is gratefully acknowledged. Project work is supported by the Australian Computational Earth Systems Simulator Major National Research Facility, The
University of Queensland and SGI. The ACcESS MNRF is funded by the
Australian Commonwealth Government and participating institutions (Uni.
of Queensland, Monash Uni., Melbourne Uni., VPAC, RMIT) and the Victorian State Government. Computations were made using the ACcESS MNRF
supercomputer, a 208 processor 1.1 TFlops SGI Altix 3700, which was funded
by the Queensland State Government Smart State Research Facility Fund and
SGI.
References
1. Abe, S., Place, D., , and Mora, P. (2004), A parallel implementation of the lattice
solid model for the simulation of rock mechanics and earthquake dynamics. Pure
Appl. Geophys. 161, 2265–2277.
16
Shane Latham, Steffen Abe and Peter Mora
2. Bowman, D. D., Ouillon, G., Sammis, C. G., Sornette, D., and Sornette, A.
(1998), An observational test of the critical earthquake concept. J. Geophys.
Res. 94, 15635–15637.
3. Bufe, C. G., Nishenko, S. P., and Varnes, D. J. (1994), Seismicity trends and
potential for large earthquakes in the Alaska-Aleutian region. Pure Appl. Geophys. 142, 83–89.
4. Bufe, C. G. and Varnes, D. (1993), Predictive modeling of the seismic cycle in
the Greater San Francisco Bay Region. J. Geophys. Res. 98, 9871–9983.
5. Cundall, P. A. and Strack, O. A. D. (1979), A discrete numerical model for
granular assemblies. Geótechnique. 29, 47–65.
6. Frye, K. M. and Marone, C. (2002), The effect of particle dimensionality on
granular friction in laboratory shear zones. Geophys. Res. Lett. 29.
7. Guo, Y. and Morgan, J. (2004), Influence of normal stress and grain shape on
granular friction: Results of discrete element simulations. J. Geoph. Res. 109.
doi:10.1029/2004JP003044.
8. Hazzard, J. F. and Mair, K. (2003), The importance of the third dimension in
granular shear. Geophys. Res. Lett. 30(13). 10.1029/2003GL017534.
9. Jaumé, S. C. and Sykes, L. R. (1999), Evolving towards a critical point: a
review of accelerating seismic moment/energy release prior to large and great
earthquakes. Pure Appl. Geophys. 155, 279–305.
10. Latham, S., Abe, S., and Davies, M. (2004), Scaling evaluation of the lattice
solid model on the SGI Altix 3700, In HPCAsia2004, Proceedings of the 7th
International Conference on High Performance Computing and Grid in the Asia
Pacific Region, pp. 226–233.
11. Mair, K., Frye, K. M., and Marone, C. (2002), Influence of grain characteristics
on the friction of granular shear zones. J. Geophys. Res. 107.
12. Message Passing Interface Forum. (1997), MPI-2: Extensions to the MessagePassing Interface.
13. Mora, P. and Place, D. (1994), Simulation of the frictional stick-slip instability.
Pure Appl. Geophys. 143, 61–87.
14. Mora, P. and Place, D. (1998), Numerical simulation of earthquake faults with
gouge: towards a comprehensive explanation for the heat flow paradox. J. Geophys. Res. 103, 21067–21089.
15. Mora, P. and Place, D. (1999), The weakness of earthquake faults. Geophys.
Res. Lett. 26(1), 123–126.
16. Mora, P., Place, D., Abe, S., and Jaumé, S. Lattice solid simulation of the
physics of earthquakes: the model, results and directions. In GeoComplexity and
the Physics of Earthquakes (ed. Rundle, J., Turcotte, D., and Klein, W.), Geophysical Monograph Series, no. 120. (American Geophys. Union, Washington
DC 2000), pp. 105–125.
17. Mora, P. and Place, D. (2002), Stress correlation function evolution in lattice
solid elasto-dynamic models of shear and fracture zones and earthquake prediction. Pure Appl. Geophys. 159, 2413–2427.
18. Morgan, J. K. (1999), Numerical simulations of granular shear zones using the
distinct element method: II. The effect of particle size distribution and interparticle friction on mechanical behaviour. J. Geoph. Res. 104, 2721–2732.
19. Morgan, J. K. and Boettcher, M. S. (1999), Numerical simulations of granular
shear zones using the distinct element method: I. Shear zone kinematics and the
micromechnics of localization. J. Geoph. Res. 104, 2703–2719.
Parallel 3D Lattice Solid Fault Gouge Simulation
17
20. Morgan, J. K. (2004), Particle dynamics simulations of rate- and statedependent frictional sliding of granular fault gouge. Pure Appl. Geophys. 161,
1877–1891.
21. Place, D. and Mora, P. (1999), The lattice solid model to simulate the physics
of rocks and earthquakes: incorporation of friction. J. Comp. Physics. 150,
332–372.
22. Place, D. and Mora, P. (2000), Numerical simulation of localisation phenomena
in a fault zone. Pure Appl. Geophys. 157, 1821–1845.
23. Sykes, L. R. and Jaumé, S. C. (1990), Seismic activity on neighboring faults
as a long-term precursor to large earthquakes in the San Francisco Bay region.
Nature. 348, 595–599.