Gamra: Simple Meshes for Complex Earthquakes

The static offsets caused by earthquakes are well described by elastostatic models with a discontinuity in the displacement along the fault. A traditional approach to model this discontinuity is to align the numerical mesh with the fault and solve the equations using finite elements. However, this distorted mesh can be difficult to generate and update. We present a new numerical method, inspired by the Immersed Interface Method, for solving the elastostatic equations with embedded discontinuities. This method has been carefully designed so that it can be used on parallel machines on an adapted finite difference grid. We have implemented this method in Gamra, a new code for earth modelling. We demonstrate the correctness of the method with analytic tests, and we demonstrate its practical performance by solving a realistic earthquake model to extremely high precision.

realistic faults, these singularities are smoothed out by non-linear processes at the fault tips that are on a scale that are many orders of magnitude smaller than the fault itself. These characteristics make it challenging to numerically model realistic fault networks.
In addition, elastostatics is only one piece of the puzzle when modeling the earthquake cycle. We want to incorporate an elastostatic solver into an overall algorithm for modeling the entire earthquake cycle [9]. We desire a unified method, using the same mesh, architecture, and boundaries, that can solve elliptic equations (for static offsets of earthquakes), parabolic equations (for poro-elastic and visco-elastic evolution between earthquakes), and hyperbolic equations (for dynamic rupture during an earthquake). Then we will have a powerful tool for self consistent models of the entire earthquake cycle.
At present, one relatively successful approach to building this kind of tool uses boundary integral methods [9,32,37,28,34,41,52,57,58]. However, boundary integral methods inevitably make simplifications in the geometry or the physics of the problem. Finite-element methods [1,26,44,51,33,31] provide a natural way to fully represent the geometry and the physics as long as the mesh conforms to the faults. Generating these conforming meshes can be quite challenging and time consuming, especially when the faults intersect. The extended finite element method [10,15,64] shows great promise in addressing this problem with mesh generation, though it has yet to be applied to realistic 3D earthquake models.
Finite difference methods, on the other hand, have not traditionally been used for this kind of problem. Straightforward implementations of finite differences require that the displacement be continuous and differentiable. This limitation spurred the development of the Immersed Interface Method (IIM) [38]. IIM explicitly models the discontinuous jump, resulting in a series of corrections to the ordinary finite difference stencils. IIM has spawned a number of variations, and some of these have been applied to various problems in elastostatics [55,56,12,65]. None of them have looked at models most relevant to earthquakes, where we prescribe the discontinuity in the displacement. More importantly, none of them have discussed how to handle the difficulties associated with the singularity at the fault tip. Finally, none of these methods have been implemented on adapted grids or parallel machines.
The purpose of this paper is to describe a new method, inspired by IIM, that naturally handles all of the difficulties associated with faults. This method was developed with an eye towards performance, so it naturally extends to the use of parallel machines and highly adapted grids. With this solver in place, we can then use the existing deep understanding of how to implement hyperbolic and parabolic solvers for the equations specific to earthquakes in a finite difference framework [17,19,20,3,18,25,49,21,22,16,35,45].
We first describe the equations of linear elasticity, how we treat internal dislocations, and how we solve these equations on an adapted mesh. Then we demonstrate the correctness of the method and our implementation with a series of analytic tests. Finally, we document the performance of our implementation with a simulation of the 1992 Mw 7. 3 Landers earthquake. The algorithm described in this paper is implemented in Gamra, a code available at https://bitbucket.org/wlandry/gamra. Gamra is a French acronym for Géodynamique Avec Maille Rafinée Adaptivement, meaning "geodynamics with adaptive mesh refinement".

Methods
We begin by describing the equations of linear elasticity (section 2.1) and the mesh we use for solving them (section 2.2). Then we describe the Gauss-Seidel smoother that we use as a component in our solvers (section 2.3). Then we describe the corrections we make to treat internal dislocations of arbitrary orientation in two and three dimensions (section 2.4). Then we describe how we implement boundary conditions (section 2.5). With these components, we have a stable, accurate solver for earthquake physics.
However, this will not be a fast solver without multigrid. To implement multigrid (section 2.6), we need coarsening (section 2.6.1) and refinement (section 2.6.2) operators. To implement adaptive multigrid, we also need to set boundary conditions at coarse-fine boundaries (section 2.6.3).

Governing Equations
We solve the Navier's equation for elastostatic deformation with the infinitesimal strain approximation where the stress components σ ji are defined using Hooke's law in terms of the displacement components v i , Lame's first parameter λ, and the shear modulus We use Einstein summation notation, where each index i, j, k is understood to x, y, and z in turn, repeated indices are summed, and commas (,) denote derivatives. For all of our test problems, the stress tensor will be symmetric (σ ij = σ ji ). In addition, the forcing term f i is zero for many of our test problems. But equivalent body forces can be used represent inelastic deformation in quasi-static deformation simulations [7,54,53]. Therefore the inclusion of body forces in Eq.
(1) is critical for modeling quasi-static deformation due to off-fault processes.

Staggered Grid
We discretize the equations on a staggered grid, with the displacement located at cell faces as shown in Figure 1. Our method requires the shear modulus (µ) at both the cell centers and cell corners. Since µ is a given function of space, we could compute it exactly at both cell centers and corners. We have found that we get larger reductions in the residuals for each multigrid V-cycle by using the given function to compute the cell centers, and then using the geometric mean to fill the value at the cell corners. Specifically, in 2D, for a reference cell where the bottom left corner is located at x = 0, y = 0, µ at that corner is µ| 0,0 = µ| δx/2,δy/2 µ| −δx/2,δy/2 µ| δx/2,−δy/2 µ| −δx/2,−δy/2 The subscripts | 0,δy/2 indicate the variable located at an offset of x = 0, y = δy/2 from the bottom left corner. So | 0,0 is the bottom left corner, | 0,δy/2 is the left face, and | δx/2,δy/2 is the cell center. The Lame parameter λ is only needed at cell centers, so there is no extra interpolation step.
We can specify µ and λ one of two ways: analytic expressions and tables. We use the muparser library [11] to evaluate analytic expressions. To compute the modulus at the boundary, we may need the modulus at a point outside the boundary. For analytic expressions, we evaluate the expression at that outside point . For moduli given by a table, we choose the closest point covered by the  table. For multigrid, the modulus on coarser levels is interpolated from finer levels, not directly computed. Using the interpolated values rather than the directly computed values results in larger reductions in the residuals for each multigrid V-cycle. The interpolation onto the cell centered modulus is a simple arithmetic average of all of the fine points in the coarse cell.
This treatment of the modulus works well for the moderate jumps in material properties seen in realistic models of earthquake regions. More extreme jumps would require a more sophisticated treatment, such as applying IIM to material interfaces as well as faults.

Gauss-Seidel Relaxation
The core of the solver is a red-black Gauss-Seidel relaxation. We first define the residual as the non-zero remnant of equation 1 We discretize the residual in the usual way with centered differences. To be explicit, in 2D, we write the x component as where, in the reference cell (v y,y λ) ,x 0,δy/2 and We then define the expression ∂r i /∂ v i | x,y as the derivative of the finite difference expression of r i with respect to v i | x,y . For example, the derivative of = λ| δx/2,δy/2 + 2 µ| δx/2,δy/2 − λ| −δx/2,δy/2 + 2 µ| −δx/2,δy/2 /δx 2 .
The Gauss-Seidel update is then given by first pass second pass We perform the update in-place in two separate passes as seen in Figure 2. Our discretisation allows us to update each point within a pass independently of each other. Parallelizing the method involves partitioning the mesh into regions that each belong to a different processor. Synchronization only happens before each pass, where each region gets updates to a single layer of ghost zones.

Theory
We define faults as a finite-sized internal surfaces where there is a displacement discontinuity called slip. Fault slip is often described in piece-wise fault segments where displacement is uniform [47,48,62,42,6,24,46], and we follow this convention. This means that a model of a realistic fault will be made up of hundreds of fault segments, each with their own slip. Internal dislocations can cause stress and displacement singularities at the edges of these segments [50,59,13]. These singularities do not manifest themselves in real earthquakes because the rock behaves nonlinearly beyond a certain stress by, for example, breaking. However, the nonlinear behavior occurs over a length scale that is orders of magnitude smaller than the rest of the model. So the stress can still get quite high, and these stress concentrations are key to understanding localized deformation. So modeling algorithms must not break down in the presence of these singularities.
To illustrate the method, consider the single faults in 2D in Figure 3. The slip − → s = (s x , s y ) on the faults is given as an input to the problem. To compute v x,x If v x is constant on each side (v right , v left ), then the slip s x is the difference between them s This goes to infinity as the resolution improves and h decreases. However, the true value of v x,x at that point is zero because v x is constant. The core idea of the original IIM paper [38] is to model these discontinuities explicitly. Then we compute corrections to apply when computing derivatives. In this case, we can compute the correct derivative by carefully subtracting away the divergent term s x /δx. Then the corrected expression is One important note is that this correction is only applied if the line between v| x+δx/2 and v| x−δx/2 crosses the fault. If it barely misses the fault as in the case at point E in Figure 3, there is no correction. This is a significant difference from other methods such as extended finite elements, which can have difficulties arising from small cell volumes or bad aspect ratios [15]. This also implies that the tip of the fault, as seen by these corrections, is only determined up to O(h).
When looking at terms with second derivatives, we build them out of first derivatives. Since the slip is constant along the fault element, there is no correction in the derivatives, only in the displacements. This means that we can build ∆ (v x,xx ), the correction for v x,xx , out of ∆ (v x,x ), the corrections for v x,x . In the reference cell, this is To be concrete, when applying this method to Eq. 5, the correction at point B in Figure 3 is The correction to Eq. 6 at point C is which is zero if the modulus λ is constant. In contrast, the correction at point D, near the tip of the fault, is because only the derivative v y,y | Dx+δx/2,Dy crosses the fault. Finally, the correction to Eq. 6 at point B is zero because each individual correction ∆ (v y,y ) is zero. Note that these corrections do not depend on the type of slip on the fault. For example, if the slip has a tensile opening component, the corrections would have the same form. The only restriction is that the two sides of the fault must be in contact. With that said, we have only tested slip along the faults, so we can only speak with certainty about that kind of slip, referred to as mode II and III in fracture mechanics.
Excluding the tips, these corrections are exact for the type of slip being modeled. This means that the stress is consistent and well behaved across the fault. We might also expect that it would lead to a scheme that converges as O h 2 . However, the method's uncertainty about the location of the tips introduces a global error that converges as O (h). At the fault tips themselves, the logarithmic singularity introduces a local error that does not converge.
The above treatment describes a single fault. Since the problem is linear, we can handle multiple faults, each made up of multiple fault segments, by adding all of the corrections from individual fault segments together. This includes the cases where fault segments intersect.

Implementation
These corrections do not depend on the computed displacement field. In that sense, they could be interpreted as body forces f i in equation 1. In 3D, this would only require 3 additional numbers per cell. However, that analogy breaks down when we consider the corrections needed when interpolating between coarse and fine levels for multigrid (Section 2.6). With that in mind, we precompute and store the jump in several different directions as shown in Figure  4. In 2D, we store the jump across a cell (∆ f ) and the jump to the corner (∆ e ). Then, for example, the correction in Eq. 9 becomes In 2D, this requires storing extra numbers per cell in addition to the 6 (v x , v y , λ, µ, f x , f y ) already needed. In 3D, we store the jump across the cell (∆ f ), from the cell face to the edge (∆ e ), and from the cell face to the corner (∆ c ). This requires extra numbers per cell in addition to the 9 already needed.

Boundary Conditions
We have implemented two different kinds of boundary conditions: Dirichlet, where the displacement is fixed to a certain value at the boundary, and stress, where the displacement is set so as to dictate what the stress is at a point. When imposing these conditions, it turns out that there is an ordering dependency among the conditions. We must first impose Dirichlet conditions. Then the shear stress conditions use values that were just set by the Dirichlet conditions. Finally, the normal stress conditions use values that were just set by the Dirichlet and shear stress conditions.

Dirichlet
The simplest boundary condition is Dirichlet conditions on the displacement normal to the boundary, as shown in Figure 5. In this case, the value at the boundary is simply set to the boundary value: For Dirichlet conditions on the displacement tangential to the boundary, as shown in Figure 5, the point outside is set so that the average of the inner and outer points equal to the boundary value: The correction ∆ ey− | x+δx/2,y is necessary to handle any faults between x+δx/2 and x. For simplicity, we define the faults to never extend out of the mesh.

Stress
A more complicated boundary condition is to set the stress rather than directly setting the displacement.

Shear Stress
The y component of the shear stress at an x boundary is We apply this condition by setting v y at an outside point This depends on v x | x,y+δy/2 and v x | x,y−δy/2 , so the normal Dirichlet condition must be applied before this condition.
Normal Stress For the normal stress in the x direction in 2D, the analytic We discretize this condition as This interpolates the derivative v y,y onto (x, y + δy/2). The moduli, λ BC and µ BC , are also interpolated there with the usual formula The condition in 3D has an additional term, v z,z , which is computed in a similar manner. This discretization depends on v y | x+δx/2,y , so the shear stress condition must be applied before this condition.

Refinement
To refine the face-centered variables, we use the stencil shown in Figure 7. We first compute a derivative of the coarse values, which in 2D is We only refine corrections to the displacement, not the displacement itself. So there is no need to add fault corrections. If we are at the boundary where one of the variables is not available, we use a one-sided derivative. For example, at y = y lower , the expression is The fine value is computed from the closest coarse value and this computed derivative v x | 0,δy/2 = V x | 0,δy − dV x | 0,δy .

Coarse-Fine Boundaries
At the interface between coarse and fine levels, we need to compute boundary conditions for the fine mesh given the coarse surrounding mesh. There are two cases of coarse-fine boundaries: vector normal to the interface (e.g., v x at an x=constant boundary), and vector tangent to the interface (e.g. , v x at a y=constant boundary). When computing these internal boundary conditions, we must use at least quadratic interpolation to keep the overall error second order [40].
Vector Normal to the Interface Figure 8 shows the stencil that is used to compute the fine boundary value on the coarse-fine interface for the component of a vector that is normal to the interface. The first step is to interpolate the coarse values to the point C. First, we define some variables where∆ are the corrections on the coarse grid. Then the coarse value at C is In 3D, the interpolation for coarse values is along diagonal directions as shown in Figure 9. That means that we can replace Eq. 10 with and then use Eq. 11 as is. Eq. 12 is only slightly modified for 3D v x | δx,δy/2,δz/2 = v x | 0,δy/2,δz/2 + ∆ f x | 0,δy/2,δz/2 If one of the coarse points is outside the physical domain, then we use a simpler interpolation. If V + is outside, then and if V − is outside then Eq. 14 is used unchanged.
Vector Tangent to the Interface Figure 10 shows the stencil used for refinement in 2D when the vector is tangential to the interface. For the case where the coarse and fine values are on the same coordinate axis, the interpolation is v x | 0,δy/2 = When the fine value does not lie along the coarse grid, we use a simple average of the neighboring coarse values

Generating the Adapted Mesh
The final part of the method is generating a mesh. Starting with a uniform grid at the coarsest resolution 1. Compute a solution on the current set of grids (section 2.6).
2. If the current number of levels is less than the maximum number of levels (a) Compute the maximum curvature at each cell center (x + δx/2, y + δy/2). The curvature in the x direction with fault corrections is At the boundaries, not all points are defined. For example, at an x = x lower Dirichlet boundary, v x | x−δx,y+δy/2 may not be defined. In these cases, we use a one-sided curvature We then compute the maximum curvature where is a fixed number, unless the maximum number of mesh refinements has been reached. Note that is an absolute rather than a relative error.
(c) Recurse back to step 1 with the new set of grids.
At fault tips, the displacement is singular and so can never be adequately resolved. However, at a finite distance from the singularity, AMR solutions can still converge [4].

Accuracy
When solving equation 1 in the presence of faults, there will always be inaccuracies because of the singularities at the tips of the faults. Away from the singularity, we expect O (h) convergence (Section 2.4.1). At the singularity, analysis becomes difficult because the Taylor series approximation breaks down. However, the scheme in Section 2.7 monitors this error and refines where needed. This means that, where the algorithm has stopped refining, the discretization error should be less than the error bound . In practice, the actual error will be larger because the local error gets integrated along the points from the boundaries and singularities. An additional source of error arises because we only approximately solve equation 1. If there is an error in the displacement ξ i , that will generate an error in the derivative v i,j of approximately ξ i /δx, where δx is the grid spacing. This implies that, for a given ξ i , the error in the stress will be at least where min (λ, µ) is the smallest value of λ or µ. The modulus does not, in our problem, vary wildly, so ∇µ µ/δx. This implies that the error in the divergence of the stress is approximately Using equation 4, we relate this to the size of the residual r i Errors in v i do not contribute to errors in f i , so that term can be neglected. Simplifying this gives an estimate for the size of the error ξ i in terms of the residual This error will become comparable to the discretization error when ξ i = , so we can turn this around to find the minimum resolution required to ensure that the solver error is smaller than the discretization error To be clear, this analysis only covers errors in solving 1 using fault segments. We do not claim to model all of the physical effects (e.g. non-linear rheologies, topography, curved faults).

Analytic Tests
We have implemented this method in the parallel, adaptive code Gamra. Gamra uses the SAMRAI framework [29,30] to handle the bookkeeping associated with multiple levels, multiple grids, and multiple parallel processes. SAMRAI is a mature, freely available, actively developed framework for large-scale parallel structured adaptive mesh refinement. SAMRAI uses MPI to coordinate work among the different processors. This has allowed us to run Gamra on a wide variety of parallel architectures: SMP nodes, traditional Linux clusters, a Blue Gene /Q, and the Intel Xeon Phi 5110p GPGPU.
In this section we perform a number of tests to ensure that the algorithm works as expected and that we implemented it correctly. We have verified that the code works in both 2D and 3D, but mostly discuss the 3D results for brevity. The tests are available from the Gamra repository 1 .
Level  Table 1: L 1 , L 2 , and L ∞ errors and L ∞ convergence rate in v x at different maximum refinement levels for the 3D expanding cylinder.

Expanding Cylinder in a Heterogeneous Medium
This is a simple test to ensure that we handle variable elastic modulus correctly.
In cylindrical symmetry, if we set the moduli and body forces to then the basis functions for solutions to Eq. 1 which are purely cylindrical with no rotation or vertical components are To make the test more rigorous, we rotate the solution by an angle θ around the y axis. Figure 12 shows a numerical solution and its associated adapted grid for a model with µ 0 = 1.4, v − = 1, v + = 0, and θ = 18 • . Table 1 shows the L 1 , L 2 , and L ∞ error in v x . While the L 1 and L 2 errors do converge, they do not converge as O h 2 . The error in the unrefined regions no longer decreases, because the mesh does not get smaller there. The integral of these small errors over the large unrefined volume is large enough to affect the overall convergence rate. This is in contrast to the L ∞ error, which converges uniformly at the expected O h 2 rate.

Internal dislocations
Okada [47,48] derived an analytic expression for the displacement due to a single fault in a homogeneous elastic half space. Figure 13 shows a solution computed by Gamra for an inclined, rotated fault. As the grid size gets more Figure 12: A cutout of the scaled displacement magnitude of a computed solution and its associated adapted mesh levels for an expanding cylinder in 3D. The axis of the cylinder is angled 18 degrees from the x axis. The model covers (−5, 1, 0) to (5,11,10). The offset is to avoid the singularity at the origin. The boundary conditions, set from the analytic solution, are Dirichlet for the normal components 2.5.1 and stress for the tangent components 2.5.2. The equivalent resolution is 128 × 128 × 128.
refined, the mesh places points closer and closer to the singularity. This means that the global L ∞ error does not shrink, but rather grows with finer resolution.
To get around this, we cut holes around the singularities and compute the L ∞ error on that region. Figure 14 shows the L ∞ error as a function of resolution. We see that the error scales as O (h) up to the point where the error becomes comparable to the criteria for adapting the mesh. Moreover, Figure 15 shows that, for a line crossing near the singularity in the displacement, the stress is well behaved. We have also run tests where we replace one of the normal Dirichlet conditions (Section 2.5.1) with a normal stress boundary condition (Section 2.5.2) set using the exact Okada stress. Similarly, we ran tests which replaced one of the shear stress conditions (Section 2.5.2) with a tangential Dirichlet condition (Section 2.5.1). All of these tests converge in a similar manner. Figure 16 shows the residual versus the number of multigrid V-cycles for 2D and 3D. In spite of the singularity at the fault tips, the solvers perform well, with the per-iteration reduction of the residual tending asymptotically to about 0.25 in 2D and 0.12 in 3D. The 3D solver uses 4 rather than 2 sweeps per multigrid level, so the absolute reduction in the residual is larger.
This gives us some confidence that all of the moving parts involved in computing the solution: smoothing (Section 2.3), boundary conditions (Section 2.5), multigrid (Section 2.6), and adaptivity (Section 2.7) are correct and implemented correctly.

Setup
We construct a realistic model of the 1992 Mw 7.3 Landers earthquake using the slip model from Fialko [23] and the material model from the Southern California Earthquake Center Community Velocity Model -Harvard (CVM-H) [60]. The slip model consists of 426 individual fault segments ( Figure 17). Figure 17 also shows the variation of Lame's first parameter, λ. The second Lame parameter, µ, has similar structure.
The boundaries are about 100-200 km away from the faults. The boundary conditions on the sides and bottom are free slip: zero shear stress (Section 2.5.2) and zero normal displacement (Section 2.5.1). The boundary condition on the top is free surface: zero shear and normal stress (Section 2.5.2). Since these boundary conditions are imperfect, the error due to the boundaries is about the size of the displacement at the boundary: 1 cm. Getting the error down to the current limits of GPS technology (about 0.5 mm [27,36,63]), would require moving the boundaries so far away such that other effects not accounted for (e.g. topography, curvature of the earth) would become significant.
During a multigrid V-cycle, we used 4 pre-and post-sweeps. On the coarsest level, we smoothed 32 times to get an approximate solution. We set the refinement criteria (Section 2.7) to our estimate of the boundary error: 1 cm. We continue multigrid V-cycles until the L ∞ norm of the residual (Eq. 4) is less than 10 −3 m GPa km −2 . From equation 16, this implies a minimum resolution of 7 · 0.01/10 −3 = 8.37 km, which in this case is satisfied when the refinement level is at least 3. The mesh is globally refined to level 3, so the error is always dominated by the discretization.

Results
Gamra automatically generated the highly adapted mesh in Figure 18. This mesh has 8.1 × 10 7 elements, while an equivalent non-adaptive mesh would require 2.2 × 10 12 elements. The computed solution in Figure 19 highlights the discontinuous nature of the solutions. We expect the error to be concentrated close to the faults, as in Figure 14. So even though the error may be larger near the faults, this would not translate to a large offset error farther from the faults. With that in mind, we expect that the error in displacement in the regions covered by levels 3-10 to be about 1 cm, or about 0.125% of the maximum displacement. Otherwise, the automatic refinement criteria would have marked those regions for refinement.

Performance
We computed this Landers earthquake solution on a Xen virtual machine running in a Dell R720 with 16 physical cores (Intel Xeon CPU E5-2670) and 256 GB of RAM using OpenMPI 1.8.8 and gcc 4.7.3. Figure 20 shows the time to solve as a function of resolution and number of cores. Altogether, the scaling is quite good at finer resolutions on this shared memory architecture.
Although it is difficult to see in the plot, we see superlinear scaling from 1 to 4 cores for finer resolution. This superlinear scaling does not persist for higher core counts. This is probably a quirk due to running inside a virtual machine. On different hardware without a virtual machine (8 physical core Intel Xeon    We can roughly fit the relation between time and grid spacing on the plot with a power law t ∝ h −1.85 . This is significantly better than a solver on a fixed three-dimensional grid. Even an optimal multigrid solver would scale as t ∝ h −3 .

Conclusion
Elastic deformation due to the displacement of faults can be modeled efficiently with parallel multigrid methods using adaptive meshes and embedded interfaces. The multigrid efficiency is commensurate with what is expected for the simpler Poisson's equation multigrid solvers [61], in spite of the added complexity brought by internal dislocations and mixed boundary conditions. The computational efficiency is improved by the mesh adaptivity, which reduces the number of nodes by orders of magnitudes compared with uniform meshes. A key advantage of the proposed method is the ability to simulate complex fault geometries without manual and labor-intensive meshing. Even in these complex models, we experienced no problems due to instabilities in the solver or excess sensitivity of the final solution to small changes in the input.
In addition, the method offers high precision in the near field of faults, even capturing the stress singularity asymptotically (Figure 15). This is important for evaluating stress and other dynamic variables. All of these features make the proposed approach optimal for generating stress and displacements kernels for inversions for fault slip [5], investigation of the surrounding elastic structure [8,14], and building stress and displacement kernels for simulations using the boundary-integral method.
This study presents an important building block of earthquake cycle simulations. A future major undertaking will be to incorporate rupture dynamics and quasi-static off-fault deformation. Fault dynamics will require modeling the propagation of seismic waves. The mesh adaptivity may then be exploited to implement spatially variable adaptive time steps [43]. Quasi-static timedependent problems with off-fault plasticity and visco-elastic or poro-elastic deformation may be treated with the same elliptic solver using equivalent body forces (per unit time), requiring only more book-keeping to handle explicit time steps. Many other effects may be incorporated to enable even more realistic models of earthquakes and Earth deformation, such as a spherical geometry for global-scale models and topography to improve calculation of local stress.

A Adaptive Multigrid
For completeness, we detail the exact adaptive multigrid algorithm we use. This is mostly a restatement of Section 4 of Martin & Cartwright [40].
First, we define a Gauss-Seidel operator GS v, f , N , where v is an initial guess, f is the forcing term, and N is the number of times to apply the smoother. The output of GS v, f , N is a correction δ v = GS v, f , N . (c) Coarsen r i to make R i (Section 2.6.1) R i = Coarsen (r i ) .
(e) Refine the correction δ V to the fine level (Section 2.6.2) and add it to the fine correction δ v δ v = δ v + Refine δ V .
(f) Apply the smoother N post times to get a final correction

Return δ v.
Given these functions, the driver routine is short.
1. Compute a composite residual r i (Equation 4). This includes applying all physical (Section 2.5) and coarse-fine (Section 2.6.3) boundary conditions.
2. While the L ∞ norm of the residual is less than the stopping tolerance stopping (a) Compute δ v = MGRelax (l max , r i ).