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.