Nonlinear Processes in Geophysics A Barnes-Hut scheme for simulating fault slip

To account for natural spatial and temporal complexity, large-scale, long-duration calculations are required for simulations of seismicity in fault zones that host large earthquakes. Without advances in computational methods, the rate of progress in “earthquake simulator” models and associated earthquake forecasts is limited by the rates at which computer speed and storage increase. To explore improvements in computational efficiency we develop the first implementation of the Barnes-Hut algorithm (Barnes and Hut, 1986) to calculate elastic interactions in a fault model. The Barnes-Hut method is an efficient, numerical scheme that treats local forces exactly and distant forces approximately. The approach is illustrated in example simulations of nonlinear fault strength in plane strain. Rudimentary error analysis indicates that efficient calculations, where execution time scales with number of grid points ( N ) as N logN , can be conducted routinely with errors on the order of 0.1%. We expect the Barnes-Hut method to be well suited for conducting initial exploration of parameter space for fault simulations with non-linear constitutive equations, and for efficient calculations of stress interaction in complex fault systems.


Introduction
A developing field in earthquake hazard research is the construction of "earthquake simulators", deterministic models that are intended to capture the statistical properties of seismicity in specific natural fault zones.Earthquake simulators use as input a particular natural network of faults from the mapped surface and subsurface expression, along with associated recurrence time or long term slip rates as constrained by paleoseismic and geodetic data (Rundle, 1988a, b;Ward, 1992Ward, , 2000;;Robinson andBenites, 1995, 1996).These models produce synthetic seismic catalogs but they Correspondence to: N. M. Beeler (nbeeler@usgs.gov)differ from the broader class of non-continuum models that includes cellular automaton (Burridge and Knopoff, 1967;Carlson and Langer, 1989a, b) and other cellular elastic models (Ben-Zion and Rice, 1993Rice, , 1995) ) in that the simulator models consider specific faults.The objective of synthetic seismicity models is to reproduce the overall statistical properties of earthquakes, e.g., frequency magnitude relations, interevent times or aftershock-foreshock statistics.Thus, the class of models that we refer to as "earthquake simulators" in this paper, such as Rundle (1988a, b), Ward (1992Ward ( , 2000)), Robinson andBenites (1995, 1996), are deterministic models of the occurrence and statistical properties of particular synthetic earthquakes within a specific natural fault network.The objective of the simulators is to reproduce the overall statistical properties of specific natural large earthquakes in the fault system of interest.
A related goal, or at least a hope, of some earthquake researchers such as select participants in the Working Group for California Earthquake Probabilities (WGCEP, 2003(WGCEP, , 2008)), is that, if the earthquake simulators can reproduce the earthquake occurrence statistics of a natural fault system, simulated earthquake probabilities can be calculated for large events which are poorly represented in the historical and paleoseismic record.Particular examples are multi-segment and multi-fault ruptures which are presently either ignored or characterized in forecasts using ad-hoc assumptions and expert opinion.Using an earthquake simulator, a synthetic catalog of multi-fault and multi-segment ruptures determined by the specified fault geometry, long term slip rates, stress interactions and on-fault frictional properties can be assembled from simulations over 1000's of earthquake cycles.Simulations, however, will need to capture the temporal and spatial complexity of natural seismicity.At present it is not known what fault properties are necessary to include in the earthquake simulators, so using them in probabilistic hazard calculations or in earthquake rupture forecasts is contingent on advances in earthquake science.In addition, the breath of temporal and spatial complexity under consideration for Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.
N. M. Beeler and T. E. Tullis: A Barnes-Hut scheme for simulating fault slip inclusion in simulator models is huge, suggesting that advances in computational approach are also needed, as follows.
In the temporal domain, natural earthquake clustering (foreshocks, aftershocks) presumably reflects on-fault physical processes that influence event occurrence times.Laboratory observations of frictional fault slip and rock fracture can explain some properties of temporal clustering.The key aspect is a small instantaneous dependence of fault strength on sliding rate that is observed during slip on pre-existing faults (Dieterich, 1979) and intact rock failure (e.g., Scholz, 1968a).The consequences of this rate dependence are profound; failure occurs following a prolonged period of gradual slip acceleration and failure time is somewhat insensitive to stress changes; both effects are due entirely to damping of the slip rate by the instantaneous rate dependence.In particular, time of failure depends on the size of the stress change and the fault's initial temporal proximity to failure (e.g., Scholz, 1972;Dieterich, 1994), resulting in what is known as "static fatigue" in the rock failure literature, (e.g., Scholz, 1972;Kranz, 1980).When this behavior is attributed to fault populations in model calculations, earthquake occurrence rates are time-dependent, such that aftershocks (Mogi, 1962;Scholz, 1968b;Knopoff, 1972;Das and Scholz, 1981;Yamashita and Knopoff, 1987;Hirata, 1987;Reushle, 1990;Marcellini, 1995Marcellini, , 1997;;Dieterich, 1994, among others) and foreshocks (Das and Scholz, 1981;Dieterich, 1994) are predicted.This class of models does not explicitly consider specific faults or fault geometry.The seismicity rate formulation of Dieterich (1994) is the most flexible of this type of synthetic seismicity model.
To incorporate laboratory-observed time-dependent friction rigorously in earthquake simulator models requires numerical solution of differential equations.As typical large earthquake periodicities exceed 100 years and require multicycle simulations to study recurrence probability, the extent that friction constitutive equations can be applied in earthquake simulator models is limited by computer speed.Consequently, published results from these kind of models to date (Rundle, 1988a, b;Ward, 1992Ward, , 2000;;Robinson andBenites, 1995, 1996) use rudimentary relations for failure and fault strength (e.g., static-kinetic friction) that lack the inherent time-dependence.
One promising development that allows for timedependent earthquake occurrence in these earthquake simulations is a fast computational scheme proposed by Dieterich (1995).The algorithm replaces nonlinear friction constitutive equations with analytical approximations; for the delayed failure described above the scheme uses a simplification that underlies the Dieterich (1994) seismicity rate formula and the precursory slip that accompanies interseismic strength recovery (Dieterich, 1972) is ignored.Furthermore, seismic fault slip is represented as quasi-static with a limit on slip speed from theory of Brune (1970).This Dieterich (1995) algorithm eliminates the numerical solu-tion of differential equations and thereby reduces computation time by orders of magnitude (Dieterich and Richards-Dinger, 2010).Similar approximations could be developed for other source processes that may operate during dynamic slip, such as pore pressurization, and particularly those that operate during the interseismic period.
In addition to natural temporal effects and the inherent computational problems they present, the overwhelming complexity of natural seismicity is in the spatial domain (Ben Zion, 2008).Mature natural faults are non-planar and segmented, having finite width and heterogeneous material properties.They are embedded in a heterogeneous and damaged elastic crust.Natural fault systems or networks are geometrically intricate collections of fault zones, in which there are earthquakes with over 10 orders of variation in moment magnitude.Because of the breadth of scale and geometrical complexity, to date the earthquake simulators are coarse and simplified.
All of the simulators in current use (e.g., Rundle, 1988a, b;Ward, 1992Ward, , 2000;;Robinson andBenites, 1995, 1996) use analytical static solutions to slip on a dislocation segment (e.g., Chinnery, 1963;Okada, 1985Okada, , 1992) ) to calculate stress interactions between segments.Calculating the static interactions places limitations on computation size or segment density because in 1-D stress at any segment center τ i changes due to the sum of slips δ j of all source segments j = 1, N as where the K ij 's are linear elastic stiffness coefficients and N is the number of fault cells.Calculating stress at each segment requires N computations for a total of N 2 over the entire grid.Thus, there are practical limits on spatial dimension or segment density in simulations of fault systems due to the computation time and memory requirements increasing strongly with the number of fault segments used.In summary, without ignoring some aspects of natural spatial or temporal complexity, without using approximate solutions or other time saving approaches, the rate of progress in earthquake simulation and associated earthquake forecasting is limited by the rates in which computer speed and storage increase.
In this paper we explore a potentially useful numerical method that intends to increase computational efficiency in earthquake simulator models.Rather than focusing on an approximate scheme for a specific type of fault friction (e.g., Dieterich, 1995), we develop a method that can be used to reduce the time of computation of any simulator.As they all use dislocation segments, improved efficiency can be achieved by using a spatially approximate solution for the elastic interactions.In approximate schemes such as the Barnes-Hut method (abbreviated BH and BHM throughout) (Barnes and Hut, 1986) and the fast multipole method (FMM) (Greengaard and Rokhlin, 1987), the objective is to relax the N 2 scaling associated with Eq. ( 1) and the implied computer storage and computation time constraints.Without going into great detail, the strategy of a FMM for dislocation modeling of fault slip would be to combine fault segments and compute interactions with other combined segments which are sufficiently far away from each other using a series expansion of the elastic stiffness coefficients.Interactions with segments that are nearby are calculated essentially without approximation.The FMM reduces the number of necessary pair-wise interactions from N 2 to order N ( Greengard and Rohklin, 1987).A method for calculating fault slip based on Warren andSalmon (1993, 1994) has been developed by Tullis and co-workers (Tullis et al., 1999).Though this method is referred to in the literature as a FMM, the number of interactions scales as N log N, rather than as N and more accurately this is a "tree" algorithm (see below) and not strictly a FMM.
Like the FMM, the Barnes-Hut method is an approximate numerical algorithm that reduces the number of mathematical operations necessary to solve a field problem (gravity, electrical, stress, etc.).Also like the FMM the BHM was used initially to calculate gravitational attraction amongst a large number of bodies (the N-body problem), but it can be used in other problems where the influence of a force decreases with distance.We develop the first explicit implementation of this method to simulate fault slip.As with the FMM, stress contributions from slip on nearby (local) fault segments are treated using the analytical solution for a dislocation and the stress contribution from the slip of distant segments uses approximate relations.Successful approximation requires that the local region about each segment is selected to be large enough to include the fault segments with the predominant stress interactions.The implementation of the Barnes-Hut algorithm is very different from the FMM.It uses no series expansion of the coefficients of the stiffness matrix K in (1) and it uses a tree computational structure.It is somewhat less efficient than the FMM, reducing the number of pair-wise interactions to order N log N, rather than N.It is, however, considerably simpler to implement relative to the FMM and the implementation is simple in the absolute sense as we hope to convey in this report.
In the following we describe a BH approach and illustrate the method in two example simulations, earthquake nucleation and periodic seismic faulting.Our calculations are not done using an earthquake simulator model, rather we wish to demonstrate that the BH method is efficient relative to the standard approach for calculating stress interactions.For our calculations we use a single planar fault and compare calculation duration and results from the method with the standard method.The simulations are conducted on an equidimensional grid with homogeneous off-fault properties in plane strain using stress interactions calculated with dislocation segment, boundary elements (Fig. 1).Plane strain is used to make the algorithm structure more easily explained and because potential scales with distance as in the analogous classic BH N-body gravity problem.Rudimentary error analysis indicates that efficient calculations which scale with grid size as N log N, can be conducted with errors of the order of 0.1%.Based on our results we expect that the Barnes-Hut method is suitable for use in earthquake simulators with long-range stress interactions.The BHM may be useful for simulating very long duration sequences of seismicity in complex fault systems, or for conducting initial exploration of parameter space for simulated faults with sophisticated constitutive behavior.The example calculations illustrate some of the tradeoffs between efficiency and error for the method.The Barnes-Hut approximation relies on forces decreasing with distance d between points of application and interest, for example for gravity as 1/d 2 (e.g., Turcotte and Schubert, 1982).At low cost to accuracy, distant forces have small influence and can be treated less rigorously than the nearby sources with larger influence.Using a gravity example, rather than calculating the effect of the attraction of each object in a distant galaxy on the Earth, the force of the summed mass of the entire galaxy can be treated as a point source applied its center of mass.The effect of the galaxy on the many other objects distant to it can be calculated in the same way.
Dislocation based faulting programs require the stiffness matrix K ij in Eq. ( 1), consisting of one coefficient for each pairwise elastic interaction.These are static stiffness coefficients and are applied in earthquake simulators in a quasistatic time stepping routine that ignores elastodynamic wave propagation.The implications of using this quasi-static approximation in earthquake simulators are not well understood and this is a topic of active research (Tullis et al., 2009).In this study we also ignore elastodyanmic effects.The K ij are distance dependent; for plane strain K ij ∝ 1/d 2 (see Eq. 6 below), and for 3-D ∝ 1/d 3 .The K ij coefficients are also proportional to the source segment width.The construction of the K matrix, consisting of N 2 calculations, is overhead prior to the actual stepping of displacement and stress (1) through time.This overhead is discussed in more detail in Sect.3.1 below.At each time step in a standard program there are an additional N 2 operations, those associated with calculating Eq. ( 1).The BHM replaces K ij with a smaller matrix.The result is that at each time step rather than there being N 2 operations there are fewer.This is where the timesavings in the method resides.There is overhead associated with setting up the smaller K matrix, but this is time well spent.Accordingly, the Barnes-Hut method consists entirely of constructing a K matrix that is much smaller but that is sufficient to produce an adequate estimate of the stresses due to fault slip.To illustrate the approach we use plane strain where the fault segment influence coefficients decay with distance from the source as in the classic Barnes-Hut gravity problem (Barnes and Hut, 1986) and the same decisionmaking criterion can be used unambiguously.
The computational algorithm uses a tree structure (Fig. 1b).Tree levels are denoted from k=1 to nn+1, where nn is an integer.The top level of the tree (k=1) is a single cell or segment with length equal to the entire fault.The second level (k=2) consists of 2 cells of equal width; these are daughter cells of the top-level single parent cell.The third level is constructed by dividing each parent on level 2 into two daughters so that tree level k=3 has 2 (k−1) =4 cells.Subsequent levels follow the same procedure.This is a binary tree.The bottom level consists of the actual fault grid of N segments.For simplicity the actual fault grid size is uniform in this example and throughout.Also, to simplify the tree structure, we only consider actual fault grids of integer powers of two, N=2 nn , and nn can be used to specify the height of the tree (nn+1) as well as the number of actual fault cells.The actual fault segments are referred to throughout as external segments or cells.The segments at upper levels of the tree (k < nn + 1) are internal segments.
The smaller BH stiffness matrix K that is analogous to K ij in (1) is constructed by sequentially traversing the tree structure from the top down for each external segment and comparing the ratio Dx/L for each internal segment to see if it is larger than a chosen value, φ.Here Dx is length of the internal segment, and L is the distance between the center of the internal and the external segments.If Dx/L < φ the internal cell is deemed far enough from the external cell to be treated approximately (distant).An influence coefficient commensurate with the internal cell's dimension and distance is used, and all its external and internal daughter cells are excluded from the computation (Fig. 2).If Dx/L > φ the internal cell is too close to the external cell, it is excluded from the computation, and 2 or more of its daughter cells will be used instead.Using this algorithm, K is generated containing a row for each external segment that consists of the influence coefficients corresponding to other external cells that are local, and approximate influence coefficients for distant cells.The value φ = 0 corresponds to the full N ×N K ij matrix in (1), and larger values to smaller K and increasingly more approximate calculations.Stress at each external cell is the product of the appropriate row of the K matrix and the corresponding external segment local values of slip and representative slips associated with the distant internal segments.In this implementation, the internal representative slips are the average of their external daughter segment slips.As a practical matter, φ is greater than or equal 0. Furthermore it also should not be much larger than 1.The long-range elastic interactions that produce coherent nucleation and growth of a simulated earthquake in the calculations in this paper disappear as φ exceeds 1 (for φ = 1, L = Dx, so even adjacent cells would be excluded from the calculation).
Though the method is described here using an equidimensional and effectively 1-D grid, it also can be used for true 2-D faults, 3-D (non-planar) faults and with irregular grid spacing.For 2-D and 3-D faults the tree structure is a quadtree.That is, for a uniform grid spacing, level k in the tree has 4 (k−1) cells, and there are N = 4 nn external cells at the highest level.For 2-or 3-D elastic problems the stiffness matrix K is N × N. A provisional 2-D fault BH code has been developed as an extension of the plane strain implementation presented in this paper.The case of 3-D requires combining the slips of non-coplanar segments, and remains a topic for future research.

Expected speedup and scaling relations
How the number of non-zero entries in K increases with increasing grid size or density determines how the duration of calculations scale with N. Scaling can be understood qualitatively by considering a single cell (red) within the grid structure as the grid extent is increased (Fig. 2); in this example φ=0.5.As the grid is increased by a factor of two, the cells used (grey) to calculate stress change rather than double, increase by only four cells, specifically as 4(nn-1).The grid has N=2 nn cells so by assuming that our single example cell is representative, the size of K is approximately 4N(nn-1) = 4N(−1 + log 2 N ) = 13.3Nlog 10 N − 4N .So for large calculations the scaling is NlogN as in the original Barnes-Hut method (Barnes and Hut, 1986).The fractional increase in execution time of the Barnes-Hut method, the speedup, should scale then as N/logN .The trend in scaling will extend to the practical time and memory limits of available computers.And since speedup is an increasing function of N , time-savings for large scale simulations should be immense.

Implementation specific properties of the Barnes-Hut scheme
Before moving on to the actual calculations a couple of differences between the problem of fault slip and the classic implementation of the BHM are worth noting.A principal difference is that for the BH gravitation problem the analogous influence matrix K is dynamic.That is, space is divided up into a static grid populated by moving bodies; so, an influence coefficient changes with time depending on the particular mass contained within the source grid element.As a consequence, the K matrix has to be recalculated at every time step.In our fault implementation the fault is fixed in space and the off-fault elastic properties are constant so the K matrix is static.K can be set up once and stored, leading to a more efficient method.A downside to this is that in the classic implementation the errors are nearly random and accumulate essentially as a random walk (Barnes and Hut, 1986) whereas for faulting the errors correlate in sign from one time step to the next leading to accumulation and the approximate solution may drift from the actual more rapidly in time.

Definitions of error and error magnitude that are used in this study
To judge the performance of the BH algorithm we track the error E associated with the method by comparing a particular solution at a grid point x at time t using the BHM, p(x,t), to the full solution f (x,t), typically, the solution without the BH approximation.Throughout we define error at a grid point as the difference between the solution of the BHM and that of the full solution, f (x,t) is assumed to be the correct solution.That is, our desire is to characterize how good the BH method is at reproducing the answer that would be otherwise obtained with the full dislocation model.The magnitude of the error M is characterized using the absolute value of the fractional error, or the percent error, 100 * M. Furthermore, we will refer to the magnitude of the error as "an error of order. . .".For example an error M=0.004 would be an error of order 0.001 or an error on the order 0.1%.

Geometry and boundary conditions
We apply the method to a fault in plane strain within a whole space (e.g., Weertman, 1979;Dieterich, 1992) (Fig. 1a).The fault is divided into N edge dislocation segments and the stress on segment i is the sum of an initial stress τ 0 i , stress resulting from self slip and slip of all the other segments τ i , and a remote uniform (tectonic) stressing where t is time and τ is the stressing rate.In the complete solution τ i is given by (1).More generally for the BHM and n i is the number of external and internal segments necessary to approximate the change in stress at external segment i due to fault slip at the particular choice of φ.To establish the errors of the BHM we start by choosing a reference solution.For this we use a numerical routine Rubin and Ampuero (2005) that employs a very accurate Fourier transform method to calculate the elastic interactions.The Fourier method stipulates periodic boundary conditions at the ends of the fault; that is, the fault is replicated in both directions (Fig. 3) an infinite number of times.Note that this choice of boundary condition is a requirement of the Fouier method that we are comparing with, rather than of the BM method Fig. 2. Schematic construction following the system of Fig. 1b, illustrating scaling as the grid extent increases from 16 to 64 cells.Outlined in red is a particular cell of interest.Gray shaded cells are those whose displacements are used to estimate the red cell's stress in a BH calculation with φ = 0.5.As the grid extent doubles, the number of gray cells increases by just 4 cells.The Table 1 summarizes the scaling which is NlogN (see text).itself.Influence coefficients consist of two terms for the actual source segment and 4 additional terms for each pair of replicate segments, one added to each side of the fault, where m is the number of replicate pairs and the distances for a single replicate pair, X ij L1 , X ij −1L1 , and X ij R1 to X ij −1R1 are as shown in Fig. 3.
To insure that solutions have small errors, calculations were conducted using 4096 grid points and 500, 5 s time steps for the BHM and for the Fourier approach.Since stress decreases with distance, for this particular geometry the in-fluence of the replicate faults is negligible beyond a few fault lengths.For calculations without the BH scheme (K ij in (3) is N × N), solutions using 4 replicate pairs are identical to 0.001% to the solution using 100 pairs.Comparisons with solutions having true periodic boundaries calculated in the wavenumber domain using Rubin and Ampuero's (2005) method show that our solutions with four replicates are identical to 0.1%.Keeping in mind that though the wavenumber approach is very efficient, it cannot accommodate nonuniform grids or boundary conditions and complex geometries that arise in nature and in earthquake simulator models.The wavenumber approach is less flexible.For the remaining comparisons in this paper we use reference solutions with Eq. ( 4), four replicates and φ = 0 except where noted; these are referred to as calculations done with the "standard" method.

Observed scaling in a faulting implementation of the Barnes-Hut method
We have conducted calculations with the Barnes-Hut and standard methods using the plane strain stiffness coefficients with periodic boundary conditions (4).As mentioned previously, the construction of the complete K ij matrix consists of calculating N 2 coefficients.This is setup time prior to the actual time-stepping calculation of displacement and stress.The setup time in our implementation of the standard method reflects the expected N 2 scaling (Fig. 4).The overhead associated with setting up the smaller K matrix for the Barnes-Hut method is in practice about ten times shorter than for the standard calculation (Fig. 4).So, there is less overhead associated with the Barnes-Hut method, despite it having a somewhat complicated tree structure.
At each time step in the standard method there are an additional N 2 operations, those associated with calculating Eq. ( 1).In calculations with grids of N = 2 11 to 2 14 the execution time in our implementation of the standard method scales as N 2 (Fig. 5).Our calculations show the BHM scales as NlogN as expected (Fig. 5).For 4096 elements the BHM is more than 10 times faster and for 16 384 elements more than 30 times faster.

Example faulting calculations with the Barnes-Hut method
For this method to be useful in large scale faulting calculations requires it to adequately reproduce the elastic interactions that it approximates.We illustrate the method in two different faulting applications: earthquake nucleation and periodic rupture of an asperity.Our interest in this paper is to document that the BHM produces reliable approximations and to quantifying the associated errors.

Nucleation
Up until recently, earthquake nucleation as inferred from laboratory experiments was thought result from slip increasing monotonically within a patch of constant size, the size being determined by the asperity contact dimension, elastic properties, and the dependence of the fault strength on slip rate.Direct observations of nucleation (Okubo andDieterich, 1981, 1984), and some plane strain simulations (Dieterich, 1992)  support this simple view.The Dieterich (1992) simulations use the particular constitutive relations for friction (Ruina, 1983).Here σ e is the effective normal stress, f is the ratio of shear to normal stress, V is slip velocity, and θ is a state-variable, which allows the fault to strengthen at very low sliding rates.The state variable has a steady-state value d c /V 0 .V 0 and f 0 are reference values of velocity and friction, respectively, and d c is a characteristic displacement associated with changes in shear resistance during sliding.Due to the complicated non-linear friction equations and spatial variations in slip, faulting simulations using (5a) and (5b) are necessarily done numerically.In Dieterich's (1992) numerical calculations the resulting nucleation patch size is constant.Extrapolations to natural stressing rates using laboratory measured parameters suggest that it would be impossible to resolve nucleation using surface and space-based strain sensors, and would be extremely difficult to detect in the subsurface using the most sophisticated borehole strain meters (Dieterich, 1992).
However, numerical simulations similar to those of Dieterich (1992) (Rubin andAmpuero, 2005;Ampuero and Rubin, 2008), show that only for relatively large negative rate dependence and high loading rates is there a constant nucleation patch size.Instead, for natural tectonic loading rates and typical lab values of rate dependence, nucleation has the form of a continuously expanding crack (Rubin and Ampuero, 2005).In their simulations there is transition from the fixed length nucleation patch to slow crack-like propagation as the ratio of the coefficients a/b in Eqs. ( 5a) and (5b) approaches 1.For this previously unrecognized regime, nucleation patch sizes can be orders of magnitude larger.Were this behavior general, nucleation might be easier to detect in the earth, and, in general, departures from fixed patch length nucleation cast some doubts on what we think we know about earthquake occurrence from lab studies (Dieterich, 1986(Dieterich, , 1992)).
During nucleation the state evolution relation can be simplified from the slip, time and slip rate dependent (5b) to an approximate form (Dieterich, 1992) which is exponential in slip.To illustrate the BHM in a first example we simulate earthquake nucleation using the slip and slip rate dependent relation for friction that results from combining (5a) and ( 5c) We assume the slip velocity is initially uniform at velocity V start except in a small region in the center of the grid where velocity is perturbed by a sinusoid of amplitude 10 • V start (Fig. 6).
During nucleation (5b) has a stable characteristic nucleation patch size (see Dieterich, 1992;Rubin and Ampuero, 2005) for ratios of a/b <0.5, and expanding, crack-like behavior for a/b >0.5.Instead, for (6), our calculation with a/b = 0.66 (Fig. 6) shows growth of a fixed length patch and generally we find that fixed length nucleation occurs for (5c) at all values of a/b.Though Rubin and Ampuero (2005) did not consider the behavior of (5c), our observation of a fixed length patch is expected based on their results.Particularly, they find that healing is not important for (5b) for low values of a/b.This lack of healing leads to the fixed length patch.Since (5c) is (5b) with the healing term removed, our result can be easily reconciled with Rubin and Ampuero (2005).The time dependent calculations carried out for this study show identical lengths for patches growing from perturbations longer and shorter than the final patch length (Fig. 6) and that wavelength is as inferred by (Dieterich, 1992) λ = 3.4Gd c /bσ n .Shown for comparison in Fig. 6a are the solutions without the BH approximation.These are in black beneath the BH solutions and are difficult to see at this scale; the BH approximation is a satisfactory representation of the solution over the entire grid.

Errors
Since distant forces are treated approximately in the BHM, successful approximation depends on 1) how well stress contributions from distant slip are approximated, and 2) how large the stress contributions from distant slip are relative the contributions from local slip.Rudimentary analysis has been conducted for a range of φ = 0.1 to 0.65 using 8192 cells.For the calculations shown in Fig. 6, which are terminated when the slip velocity reaches 1 mm/s, there are many orders of magnitude variation in slip and slip velocity over the grid.We consider the error at each grid point as the difference between the Barnes-Hut and the standard solution at that point, normalized by the standard solution as described in Sect.2.3.
For this range of φ the maximum errors in slip range from 0.3% to 25% and average errors range from 0.06% to 4%.To achieve maximum errors of 1% or better, φ < 0.2 is required.
To begin to understand the errors in the solution and how they vary with φ, first consider errors associated with replacing external cells with a parent; these indicate how well distant forces are approximated.Ignoring contributions from the replicates, (4) can be rewritten as and we see that the influence coefficients increase in proportion to the cell size Dx.Recall that Dx is a multiple n of the external cell size dx.When an external (child) cell is represented as part of a larger internal (parent) cell, stresses due to slip on that external cell will be represented by 1/n of the stiffness coefficient of the internal cell K parent .In the full problem this same slip would be represented by K child .The associated fractional error in the calculated stress due to using the internal cell instead of the external cell is The actual maximum errors associated with the influence coefficients as calculated from ( 7) are shown in Fig. 7a.These errors E K for φ = 0.1 to 0.65 (thick grey line) are large, ranging from 10% to 1% (Fig. 7a).They are much larger than the actual maximum errors in slip observed in the calculation (blue, Fig. 7a), and are still larger than the average errors in slip (green, Fig. 7a).The reason for the difference between the errors in stiffness and in slip is that the actual errors in slip result only in part from the error explicitly associated with the approximation.The further error reduction results because K ij for the nearest cells are not approximated and the nearest cells have the have the largest influence.Thus, Eq. ( 7) represents an upper and unrealized bound on the errors in this BHM.
To get a better idea of the source of error in the BH approximation, consider the absolute value of the ratio of the maximum parent influence coefficient to the maximum influence coefficient.As described in the appendix, this ratio, R, can be calculated analytically for our uniformly spaced grid, R is a small number ranging between 0.001 and 0.06 as φ increases from 0.1 to 0.65 (Fig. 7a).The ratio signifies that the distant and approximate influence coefficients are small relative to the nearby exact coefficients.This explains qualitatively why the maximum error in the stiffness coefficients E k greatly over-estimates the observed errors.For this nucleation problem Eq. ( 8) shows nearly the same variation with φ as the maximum and average slip.Empirically we can use 3.6R as an estimate of the maximum error in slip and 0.62R for the average error.These estimates of the errors are the black curves in Fig. 7a.To insure an error of 1% in slip requires φ < 0.2.Failure time, defined in these calculations as the time slip velocity reaches 1 mm/s, has a smaller error than slip.The BH solution with φ = 0.1 reproduces failure to 0.01% and executes in 6.5% of the time.Since these particular calculations www.nonlin-processes-geophys.net/18/133/2011/ Nonlin.Processes Geophys., 18,[133][134][135][136][137][138][139][140][141][142][143][144][145][146]2011 142 N. M. Beeler and T. E. Tullis: A Barnes-Hut scheme for simulating fault slip are done with a modest number of segments (N = 2048) and the fractional increase in execution time of the BHM, the speedup, should scale as N /logN, we expect more significant time savings for larger scale simulations.1% error on failure time can be achieved even with φ = 0.5 (Fig. 7b).Estimates of failure time can be efficiently calculated, preserving the overall behavior despite being approximate in detail.Again, R (8) shows nearly the same variation with φ is failure time.
Empirically we can use 0.29R as an estimate of the error in failure time.This error estimate is shown as the black line in Fig. 7b.

Periodic asperity failure
High-resolution seismological studies of creeping faults have revealed small repeating earthquakes with short recurrence intervals (typically t r <5 yrs).These repeating earthquakes occur on different faults and in different tectonic settings, for instance on the Calaveras fault in the aftershock zone of 1984 M 6.1 Morgan Hill, California earthquake (Vidale et al., 1994), within the aftershocks of the 1989 M 6.9 Loma Prieta, California earthquake (Schaff et al., 1998), on the Parkfield, California segment of the San Andreas fault (Nadeau and Johnson, 1998) where one such event is the target of the SAFOD project, on the Hayward fault in Northern California (Burgmann et al., 2000) and in subduction zones (Matsuzawa et al., 2002;Uchida et al., 2005).Failure is thought to occur on a single asperity or fault patch embedded in an aseismically creeping fault plane (Vidale et al., 1994;Ellsworth, 1995;Nadeau and Johnson, 1998;Nadeau and McEvilly, 1999).In this case, recurrence is influenced by the rate of aseismic creep of the fault surrounding the patch, as geodetically measured fault creep correlates with cumulative earthquake moment release rate and with recurrence interval (e.g., Ellsworth, 1995).Study of small earthquakes with short recurrence may provide insight into the behavior of larger repeating earthquakes that are of greater interest because of their greater damage potential, but are more difficult to study because they have much longer recurrence.On the other hand, small repeating earthquakes may be controlled by different fault physics due to their small slip and lower frictional heating.In that case source parameters and recurrence statistics of small repeating earthquakes may not be directly related to the large and hazardous events.
Our second example for testing the Barnes-Hut approximation is periodic failure of a small earthquake embedded on an otherwise aseismically creeping fault plane.The geometry is similar to 3D simulations by Chen and Lapusta ( 2009), but our calculations are plane strain and fully quasi-static while Chen and Lapusta (2009) use fully elastodynamic equations during rupture propagation.We simulate periodic asperity failure using Eqs.(5a) and (5b).The fault is loaded by an imposed constant slip rate V plate at the ends of the grid.The asperity has a < b, whereas outside the asperity b = 0.This configuration results in continuous aseismic slip of the fault well outside the asperity and stick-slip within the asperity.To account for radiated energy resulting from rapid slip we use the radiation damping approach of Rice (1993).Accordingly, stress change due to slip at each segment i is where the effective shear impedance η = µ 2β, and µ and β are the shear modulus and wave speed, respectively.For a 10 m patch in the center of a 1000 m long fault with V plate = 7.927 × 10 −10 m/s, the fault slips periodically with recurrence t r = 27.9 days (Fig. 8a).The space/time history of slip in the region containing the patch is well approximated with the BH approximation over the entire grid (Fig. 8b).Nucleation is different from the fixed width nucleation case (Fig. 6).Loading at the asperity edges leads to creep near the inner edge of the patch at nearly the same rate as outside the patch.Nucleation proceeds initially as the higher slip rate eats in towards the center from the edges.Over time the center of the patch accelerates.The region of accelerating slip is initially nearly the full asperity dimension but narrows with time.In the last stages of nucleation the patch length is similar in width to that in the simple nucleation example.Despite the patch being loaded directly at the patch edge, nucleation of rupture ultimately occurs at the patch center and rupture propagates outward bilaterally to the edges.Rapid slip extends outside the patch into the rate strengthening surroundings.In these calculations, deceleration and arrest proceeds preferentially in the patch center while the edges continue to creep at an appreciable rate.This creep at the edges is essentially the early stages of afterslip in the surroundings.Note that aspects of propagation and arrest in nature may be influenced by elastodynamics not included in these calculations.

Errors
We calculated the errors for φ = 0.1 to 0.65 using 8192 cells.To examine errors, as in the nucleation case, we stop the simulation when the slip velocity reaches 1 mm/s.In this implementation slip is not calculated so we determine the maximum errors in stress (Fig. 9) which range from 0.8% to 18%.The errors in stress are not as systematic as the errors in slip in the nucleation example.To achieve maximum errors of 1% or better errors, φ <0.2 is required.As in the nucleation example, the actual maximum errors associated with the influence coefficients E k as calculated from (7) are shown for reference (Fig. 9, thick grey line) and they are larger than the actual maximum errors in stress observed in the calculation.Shown also for reference is the ratio R. As in the nucleation example E k and R bound the errors in stress.Failure time (Fig. 9, red), has a smaller error than stress.The BH solution with φ = 0.1 reproduces the recurrence on the order of 0.1% and executes in 2.3% of the time.As with the nucleation example, the fractional increase in execution time of the BHM, Nonlin.Processes Geophys., 18,[133][134][135][136][137][138][139][140][141][142][143][144][145][146]2011 www.nonlin-processes-geophys.net/18/133/2011/ the speedup, scales as N/logN, and we expect more significant time savings for larger scale simulations.1% error on failure time can be achieved even with φ = 0. 5 (Fig. 7b).Again as with the nucleation example, results can be calculated efficiently that preserve the overall behavior despite being approximate in detail.

Conclusions/future work
Our interest in the Barnes-Hut method is to improve computational efficiency for developing "earthquake simulators".
The simulators use natural fault network geometries, longterm slip rates and recurrence data along with particular assumptions about source physics and friction to produce synthetic seismicity catalogs.The objective of the simulators is to reproduce the overall statistical properties of earthquakes, e.g., frequency magnitude relations, interevent times or aftershock-foreshock statistics associated with particular synthetic earthquakes on specific faults within an extensive model of a natural fault system.Our proposed approach to improving the efficiency of earthquake simulators is to use a spatially approximate solution technique for calculating the elastic interactions between fault segments.Our Barnes-Hut scheme for calculating elastic interactions due to fault slip was implemented in test cases using non-linear fault constitutive relations that depend on time, slip and slip rate.The results are that quasi-static simulations of seismic and aseismic faulting can be done routinely with 0.1% errors in slip, stress and time and that calculation duration increases with increasing number of fault segments N as NlogN.The method is more than an order of magnitude faster than the standard approach for N > 4000.We have not used the Barnes-Hut method in an earthquake

Fig. 1 .
Fig. 1.Fault geometry, cell and tree structure of the plane strain Barnes-Hut implementation.(a) Schematic representation of a fault in a whole space in plane strain.This is the geometry that is used in all the calculations in this paper.The fault is in the plane containing the z and x axes; the y direction is the fault normal.It has a finite length in the x direction and extends infinitely in the positive and negative z directions.The red dashed line, shown for reference, is the intersection of the fault with the plane containing the y and x axes at z = 0. (b) Cell and tree representation of the fault.Left panel shows the fault trace from (a) in red, that is, the intersection of the fault with the plane containing the y and x axes, corresponding to the red dashed line in (a).The fault is divided into 2 3 equi-dimensioned external cells.Above the fault in the left panel are the 2 3 − 1 internal cells, labeled by level.Right panel shows the internal and external cells schematically as a tree structure, a choice of corresponding cell numbers, and lines connecting parent to daughter cells.
136 N. M. Beeler and T. E. Tullis: A Barnes-Hut scheme for simulating fault slip 2 The Barnes-Hut method

Fig. 3 .
Fig. 3. Example grid, likeDieterich (1992) andWeertman (1979) but with periodic boundary conditions.This grid is used in all calculations in the present study.The periodic boundary conditions are used to approximate solutions obtained with an accurate Fourier method (Rubin andAmpuero, 2005).

Fig. 4 .
Fig. 4. Scaling of set-up time for boundary element calculations with the number of cells (N) used in the calculation.This plot shows how the observed time necessary to construct the stiffness matrix (e.g., Eq. 1), (setup time) increases with the number of cells N .Red symbols are for the standard implementation.The red line has a slope of 2 corresponding to scaling of N 2 .The black symbols define the scaling of setup for the BHM.The absolute setup duration is smaller for the BHM than for the standard approach.

Fig. 5 .
Fig. 5. Scaling of execution time with number of elements for nucleation with the Barnes-Hut scheme for φ = 0.5 (black symbols), and with the standard method (red symbols) for 100 time steps.Shown for comparisons are lines following N 2 and N logN scaling.

Fig. 6 .
Fig. 6.Rupture initiation using slip and slip rate dependent fault strength (5a) and (5c) and the Barnes-Hut scheme with φ = 0.5 on a plane strain fault grid of 40 m length and periodic boundary conditions.In these calculations f 0 = 0.7 a = 0.008, b = 0.012, d c = 1 µm, σ n = 6 MPa τ = 0, G = 2.1 × 10 4 MPa, and ν = 0.25.(a) Initial slip rate (black) is uniform at V s = 1 × 10 −12 m/s outside a centered perturbation of 0.25 m width.The form of the perturbation is one half of 0.5 m wave length cosine function of peak to peak amplitude 10•V start .Different curves show increasing slip velocity within the patch with time.The solution without the Barnes-Hut approximation is also shown but differences between it and BH are not really visible at this scale.Shown for reference is the dominant wavelength as inferred by Dieterich (1992).(b) Same as in (a) except the width of the perturbation (black) is 8 m.

Fig. 7 .
Fig. 7.Errors in the BHM for nucleation using the Barnes-Hut scheme as φ is varied between 0.1 and 0.65.(a) Errors associated with the stiffness coefficients (E k , Eq. 7 thick grey line) and with fault slip (M, blue and green).The M are from the numerical simulations.Shown also is R (Eq. 8, grey) which is the ratio of the maximum parent stiffness coefficient to the maximum stiffness coefficient.The black curves are fits of the form of R to the observed maximum and average errors in slip.These correspond to 3.6R and 0.62R.(b).Errors associated with time of failure (M, red).These are from the numerical simulations.Shown also is R (Eq. 8, grey).The black curve is a fit of the form of R to the observed failure time error.

Fig. 8 .
Fig. 8. Simulation of repeated slip of a seismic asperity embedded in an otherwise aseismically creeping fault.Equations (5a) and (5b) are used in the Barnes-Hut scheme on a plane strain fault grid of 1000 m length.The asperity is 10 m long.The calculation approximates periodic boundary conditions at the fault ends, that is, the fault is repeated along strike in both directions.For a true periodic boundary, the fault is repeated infinitely in both directions; in this calculation there are 1000 repeats of the fault in each direction The fault is loaded at both ends at constant rate V plate = 7.927 × 10 −10 m/s. a = 0.008, d c = 10 µm, σ n = 60 MPa τ = 0, G = 2.1 × 10 4 MPa, β = 3 km/s, and ν = 0.25.Inside the patch f 0 = 0.7 and b = 0.012, outside the patch f 0 = 0.2 and b = 0.0.(a) Time history of slip velocity on the element in the center of the patch from the BH approximate (φ = 0.1, black) and the standard solution (φ = 0.0, grey).The slight differences are difficult to see at this scale.The recurrence is characteristic with period ∼27.9 days.(b) Time-space evolution of slip rate about the patch during a single recurrence.Contours are spaced by changes of an order of magnitude at the center element.For φ = 0.1, black curves denote behaviour during nucleation, the lower most black curve corresponds to the initial stages of nucleation.For φ = 0.1, the red curves show the stages of rupture arrest and afterslip.Also shown in grey are the corresponding solutions for φ = 0.0.

Fig. 9 .
Fig.9.Errors in the BHM for repeated slip of an asperity as φ is varied between 0.1 and 0.6.Equations (5a) and (5b) are used in the Barnes-Hut scheme on a plane strain fault grid of 40 m length.The asperity is 10 m long.The errors associated with the stiffness coefficients (E k , Eq. 7 thick grey line), with shear stress (M, blue) and with failure time (M, red) are shown.These are from the numerical simulations.Shown also is R (Eq. 8, grey).

Table 1 .
The summary of the scaling which is N logN.