Gamra: Simple Meshing for Complex Earthquakes
Abstract
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 [38], 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.
1 Motivation
A common feature of many earthquakes is a complex network of intersecting faults. Accurately modeling the static offsets and associated large scale deformation due to this fault geometry is crucial to a reliable understanding of seismic hazards [39]. The behavior of these faults is relatively well described by the equations of variable modulus elastostatics. However, for realistic faults, the displacement does not gradually taper off, but rather ends abruptly. This abrupt termination gives rise to a logarithmic singularity in the displacement [48]. In realistic faults, these singularities are smoothed out by nonlinear 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 poroelastic and viscoelastic 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. Finiteelement 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”.
2 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 GaussSeidel 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 coarsefine boundaries (section 2.6.3).
2.1 Governing Equations
We solve the Navier’s equation for elastostatic deformation with the infinitesimal strain approximation
(1) 
where the stress components are defined using Hooke’s law in terms of the displacement components , Lame’s first parameter , and the shear modulus
(2) 
We use Einstein summation notation, where each index , , is understood to , , and in turn, repeated indices are summed, and commas (,) denote derivatives.
For all of our test problems, the stress tensor will be symmetric . In addition, the forcing term is zero for many of our test problems. But equivalent body forces can be used represent inelastic deformation in quasistatic deformation simulations [7, 54, 53]. Therefore the inclusion of body forces in Eq. (1) is critical for modeling quasistatic deformation due to offfault processes.
2.2 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 Vcycle 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 , , at that corner is
(3) 
The subscripts indicate the variable located at an offset of , from the bottom left corner. So is the bottom left corner, is the left face, and 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 Vcycle. 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.
2.3 GaussSeidel Relaxation
The core of the solver is a redblack GaussSeidel relaxation. We first define the residual as the nonzero remnant of equation 1
(4) 
We discretize the residual in the usual way with centered differences. To be explicit, in 2D, we write the component as
where, in the reference cell
(5)  
(6)  
and
We then define the expression as the derivative of the finite difference expression of with respect to . For example, the derivative of is
The GaussSeidel update is then given by
(7) 
We perform the update inplace 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.
2.4 Treatment of Internal Dislocations
2.4.1 Theory
We define faults as a finitesized internal surfaces where there is a displacement discontinuity called slip. Fault slip is often described in piecewise 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 on the faults is given as an input to the problem. To compute at point , we would ordinarily write the finite difference expression
If is constant on each side , then the slip is the difference between them . The finite difference then becomes
This goes to infinity as the resolution improves and decreases. However, the true value of at that point is zero because 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 . Then the corrected expression is
One important note is that this correction is only applied if the line between and 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 .
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 , the correction for , out of , the corrections for . In the reference cell, this is
(8) 
To be concrete, when applying this method to Eq. 5, the correction at point B in Figure 3 is
(9) 
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
crosses the fault. Finally, the correction to Eq. 6 at point B is zero because each individual correction 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 . However, the method’s uncertainty about the location of the tips introduces a global error that converges as . 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.
2.4.2 Implementation
These corrections do not depend on the computed displacement field. In that sense, they could be interpreted as body forces 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 () and the jump to the corner (). Then, for example, the correction in Eq. 9 becomes
In 2D, this requires storing
extra numbers per cell in addition to the 6 (, , , , , ) already needed. In 3D, we store the jump across the cell (), from the cell face to the edge (), and from the cell face to the corner (). This requires
extra numbers per cell in addition to the 9 already needed.
2.5 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.
2.5.1 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 is necessary to handle any faults between and . For simplicity, we define the faults to never extend out of the mesh.
2.5.2 Stress
A more complicated boundary condition is to set the stress rather than directly setting the displacement.
Shear Stress
The component of the shear stress at an boundary is
We apply this condition by setting at an outside point
This depends on and , so the normal Dirichlet condition must be applied before this condition.
Normal Stress
For the normal stress in the direction in 2D, the analytic condition is , which implies
We discretize this condition as
This interpolates the derivative onto . The moduli, and , are also interpolated there with the usual formula
The condition in 3D has an additional term, , which is computed in a similar manner. This discretization depends on , so the shear stress condition must be applied before this condition.
2.6 Multigrid on an Adapted Mesh
With a smoother (Section 2.3), corrections for faults (Section 2.4), and boundary conditons (Section 2.5), we can compute highly accurate solutions to Eq. 1 on a single grid. This will, however, be very slow. To shorten the time to solution, we implement adaptive multigrid (Appendix A). This is essentially an enhancement of the multigrid method for adapted grids. To implement this, we must first implement coarsening, refinement, and coarsefine boundary operators.
2.6.1 Coarsening
Following Albers [2] we use weighted arithmetic averages to coarsen the face centered displacement and residuals. Figure 6 shows the fine values used to compute the coarse value for . The corresponding expression in the reference cell is
The expression in 3D is a straightforward extension
At physical boundaries where not all of the values are available, we average only over the face. In 2D, the expression is
and in 3D it is
2.6.2 Refinement
To refine the facecentered 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 onesided derivative. For example, at , the expression is
The fine value is computed from the closest coarse value and this computed derivative
In 3D, the expressions look very similar although now we interpolate along diagonals. For a fine variable on a coarse face, the derivative is
and the fine value is
For fine variables in between coarse faces, we average the fine values on each coarse face:
2.6.3 CoarseFine 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 coarsefine boundaries: vector normal to the interface (e.g., at an x=constant boundary), and vector tangent to the interface (e.g., 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 coarsefine 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
(10) 
where are the corrections on the coarse grid. Then the coarse value at C is
(11) 
The final step is to interpolate along a line to get the fine value at F
(12)  
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
(13) 
and then use Eq. 11 as is. Eq. 12 is only slightly modified for 3D
(14)  
If one of the coarse points is outside the physical domain, then we use a simpler interpolation. If is outside, then
and if 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
(15)  
When the fine value does not lie along the coarse grid, we use a simple average of the neighboring coarse values
and the interpolation becomes
At the or corner, some of the fine corrections (e.g. ) are not necessarily defined. For the boundary, we work around this by correcting the coarse value at to first, and then using the same correction from to . With this, the interpolation becomes
Figure 11 shows the points used for refinement in 3D when the coarse and fine values are on the same coordinate axis. Defining
2.7 Generating the Adapted Mesh
The final part of the method is generating a mesh. Starting with a uniform grid at the coarsest resolution

Compute a solution on the current set of grids (section 2.6).

If the current number of levels is less than the maximum number of levels

Compute the maximum curvature at each cell center . The curvature in the direction with fault corrections is
At the boundaries, not all points are defined. For example, at an Dirichlet boundary, may not be defined. In these cases, we use a onesided curvature
