# Seismological Research Letters

- © 2016 by the Seismological Society of America

## ABSTRACT

It has recently become practical to use triangular fault elements in earthquake simulators instead of the usual rectangular fault elements. A simulator spends most of its time computing how slip on some fault elements affects the stresses on other fault elements. We explore whether rectangles or triangles yield more accurate stress values on curved fault surfaces. One might expect triangles to perform better, because triangles can follow the shape of a curved surface, whereas rectangles leave gaps and overlaps between adjacent fault elements. Our test results show that, contrary to expectations, rectangles overall perform as well as or better than triangles when computing stresses on curved fault surfaces. We also find that one triangulation may perform significantly better than another triangulation.

*Online Material: *Figures and tables describing a larger set of results than in the main article, and step‐by‐step algorithms for constructing fault discretizations.

## INTRODUCTION

An earthquake simulator is a computer program that generates a synthetic earthquake catalog spanning thousands of years or longer. Most of the computational effort in an earthquake simulator goes into computing how slip on one part of a fault affects stresses on other parts of the fault and on other faults. The computation is done by discretizing the fault system into a large number of fault elements and using Green’s functions to determine how a pattern of slip on some fault elements affects the stresses on all the fault elements (Pollitz, 2012; Richards‐Dinger and Dieterich, 2012; Sachs *et al.*, 2012; Tullis *et al.*, 2012; Ward, 2012).

Traditionally, earthquake simulators have used rectangular fault elements, chosen so that the Okada Green’s functions can be used (Okada, 1992). Recently, because of the development of new Green’s functions for triangular dislocations, it has become practical to use triangular fault elements (Meade, 2007; Gimbutas *et al.*, 2012). The purpose of this project is to assess the accuracy of stress calculations performed with triangular fault elements, as compared to the accuracy of the same calculations done with rectangular fault elements.

For planar faults, rectangles and triangles can be expected to give the same results (except possibly if the fault has an irregular boundary, which triangles can follow more accurately than rectangles); however, for curved faults, rectangles and triangles give different results. When a fault is curved, partitioning it into rectangular fault elements will necessarily create gaps and overlaps between adjacent elements. (It should be noted that the Okada Green’s functions do not work with arbitrary rectangles, but instead require rectangles for which the upper and lower edges are horizontal. Meeting this constraint further exacerbates the creation of gaps and overlaps.) In contrast, partitioning a curved fault into triangular fault elements can be done using a triangular mesh, which has no gaps or overlaps between adjacent elements.

Because triangles can represent curved fault geometry more accurately than rectangles, one intuitively expects that stress calculations performed with triangles should be more accurate than stress calculations done with rectangles. However, our results are contrary to the intuitive expectation. In our tests, triangles are not superior to rectangles. One or the other may be superior in a particular case, but, overall, rectangles perform as well as or better than triangles. Another unexpected result is that one triangulation of a fault surface may perform significantly better than another triangulation with a different pattern of triangles.

## METHODS

We performed stress computations on three curved fault surfaces: one with negative curvature, one with positive curvature, and one in the shape of a cylinder (Fig. 1).

Our negatively curved surface is motivated by the shape of the San Andreas fault in southern California (Fuis *et al.*, 2012). We ran our tests on a strike‐slip fault in the shape of a helicoid, with dip angles ranging from +45° (at one end of the fault) to −45° (at the other end of the fault), passing through 90° (at the center of the fault). The fault trace is a 38,400‐m‐long straight line, and the fault extends from the Earth’s surface to a depth of 19,200 m. The rake angle on the helicoidal fault is a constant 180°. We use the conventions of Aki and Richards (2002) for strike, dip, and rake.

Our positively curved surface is motivated by the shape of the Cascadia subduction fault (McCrory *et al.*, 2012). We ran our tests on a thrust fault in the shape of an ellipsoid, with dip angles ranging from 10° (at the Earth’s surface) to 30° (at the base of the fault). The fault trace is a 38,400‐m‐long arc, with strike angles ranging from +30° to −30°. The fault extends from the Earth’s surface to a depth of ∼11,065 m so that the fault depth is 38,400 m, as measured along dip on the curved surface. The rake angle on the ellipsoidal fault is a constant 90°.

For the helicoidal and ellipsoidal surfaces, we used greater curvature than the natural fault to emphasize the effects of the curvature.

Our cylindrical surface is the actual shape of the North Frontal West fault (NFWF), located in southern California ∼95 km east‐northeast of Los Angeles. We obtained the fault geometry from the third Uniform California Earthquake Rupture Forecast (UCERF3), fault model v.3.2 (see Data and Resources) (Field *et al.*, 2014). We chose the NFWF because of its unusually irregular shape, making it one of the more difficult faults to tile with rectangles. The shape of this fault is a cylinder with an irregular cross section, which dips at an angle of 49°. (Mathematically, a cylinder is considered to have zero curvature because it can be unrolled into a plane without distorting the surface.) Although we use the actual shape of the NFWF, we magnify it by a factor of ∼1.846 so that we can use the same‐sized fault elements as for our other two model faults. After magnification, the fault trace is an irregular curve with an arc length of 90,666 m. It extends from the Earth’s surface to a depth of ∼28,981 m, so that the depth of the fault is 38,400 m as measured parallel to the cylinder axis along dip. UCERF3 assigns the NFWF a rake angle of 90°. We interpret this to mean that fault slip is parallel to the cylinder axis. Because of the fault’s irregular shape, the rake angle must vary over the fault surface to maintain slip parallel to the cylinder axis. Accordingly in our tests, we use a rake angle that varies between ∼52° and 127°.

Our strategy is to impose uniform slip on part of the fault surface and then compute the induced shear and normal stresses elsewhere on the fault surface. We examined the induced stresses in two places: immediately adjacent to the slipping portion of the fault and 4.8 km away from the slipping portion. The induced stresses immediately adjacent to the slipping area are relevant for the simulation of rupture propagation during an earthquake, because the earthquake rupture proceeds from the currently slipping elements to their neighboring elements. The stresses 4.8 km away from the slipping area are relevant for examining the accuracy of stress transfer over larger distances, particularly for simulators that employ rate–state friction, in which stress changes influence the evolution of fault state even before the onset of seismic slip. We chose a separation of 4.8 km because it is roughly the maximum size of a fault‐to‐fault jump (Wesnousky, 2006).

To assess the accuracy of the calculation, we discretize the fault surfaces using either rectangles or triangles, each in six different sizes: 300, 600, 1200, 2400, 4800, and 9600 m. The 300 m results are taken as the reference solution, and the accuracy of the other five sets of results is assessed by how well they approximate the reference value. We now describe in detail how this is done.

### Accuracy Assessments

Earthquake simulators use fault elements in two ways: as sources and as targets. For each pair of fault elements S and T, the simulator computes the shear and normal stress acting on the centroid of the target element T, due to fault slip occurring on the source element S. Each fault element acts as both a source and a target. (Some simulators use both shear and normal stress, whereas other simulators use only the shear stress.)

When the fault is curved, the fault elements can only approximate the curved surface. When computing stresses, there is some error due to the fact that the source element does not lie precisely in the surface, and there is some error due to the fact that the target element does not lie precisely in the surface. We perform two tests to separately evaluate these two causes of error and determine whether triangular or rectangular elements are better:

A source test evaluates the error originating at the source elements. This test is performed using a source region and a target region that are separated by 4800 m. We vary the discretization of the source region, using either triangles or rectangles with sizes ranging from 600 to 9600 m. The discretization of the target region is held fixed. In this article, the target region is discretized with 600 m rectangles; other discretizations yield similar results. In each case, a reference solution is calculated by discretizing the source region with 300 m elements.

A target test evaluates the error originating at the target elements. This test is also performed using a source region and target region that are separated by 4800 m. The discretization of the source region is held fixed. In this article, the source region is discretized with 300 m rectangles. We vary the discretization of the target region, using either triangles or rectangles with sizes ranging from 600 to 9600 m. In each case, a reference solution is calculated using 300 m elements in the target region, placed so that their centroids align with the centroids of the larger elements. (There will be large gaps between the 300 m elements, because the number of 300 m elements is the same as the number of larger elements.) In order to have a valid comparison, it is necessary to preserve the centroid locations, because simulators evaluate the induced shear and normal stress at the centroids of the target elements.

When an earthquake rupture is propagating along a fault surface, the rupture expands from the region that is currently slipping into the immediately adjacent fault elements, so it is important for a simulator to be able to accurately compute stresses in the fault elements adjacent to a slipping region. We use a third test to evaluate whether rectangles or triangles are better for this important case:

A propagation test evaluates errors occurring during rupture propagation. This test is performed using a target region that borders on the source region, with no separation. We vary the discretization of both the source region and the target region, using either triangles or rectangles with sizes ranging from 600 to 9600 m. The source and target regions are discretized in the same way (e.g., if one is 1200 m triangles, then so is the other), as would be the case in a simulator. In each case, a reference solution is calculated by discretizing both the source and target regions with 300 m elements. In the target region, the centroids of the 300 m elements are aligned with the centroids of the larger elements, as in the target test.

In all our tests, we calculate stresses using the Green’s functions of Okada (1992) if the source region consists of rectangular elements or the Green’s functions of Gimbutas *et al.* (2012) if the source region consists of triangular elements.

### Error Metric

A calculation produces one stress value for each target fault element. If there are *N* target elements, then the stress values can be expressed as an *N*‐dimensional vector (*s*_{1},…,*s*_{N}). The components *s*_{i} can be either shear stresses or normal stresses. We handle shear stresses and normal stresses separately, so that the vector (*s*_{1},…,*s*_{N}) consists of either all shear stresses or all normal stresses. We also have available a reference solution computed using 300‐m‐fault elements of the same shape and configuration that were used to compute the vector (*s*_{1},…,*s*_{N}).

In order to produce a numerical measure of accuracy, we use a modified version of the *Q* metric developed by the Southern California Earthquake Center Dynamic Rupture Code Verification Project (Barall and Harris, 2015). The *Q* metric value can range from 0% (for perfect agreement) to 200% (for vectors pointing in opposite directions); lower values indicate better agreement. Suppose **a**=(*a*_{1},…,*a*_{N}) is an *N*‐dimensional vector. Define the *L*^{1} norm of **a** to be

Then, given two *N*‐dimensional vectors **a** and **b**, we compare them using the *Q* metric, defined as

In order to express *Q* as a percentage, the above value is multiplied by 200. All our results use the *Q* metric, expressed as a percentage, to measure the accuracy of a stress calculation.

## FAULT GEOMETRY AND DISCRETIZATION

In order to discretize a fault at multiple resolutions, we need a mathematical description of the fault surface geometry and a procedure for partitioning the surface into fault elements. In order to perform our tests, we also need to specify the location of the source and target regions.

### Mathematical Description of Fault Geometry

The fault is embedded in a 3D space with (*x*, *y*, *z*) coordinates. The *z* coordinate is vertical and points upward, so negative *z* values are underground. We introduce a fourth coordinate *d*, which represents distance along strike (Barall, 2012). Every point on the fault surface is identified by the coordinate pair (*z*, *d*). In order to specify the 3D shape of the fault, we must specify the mapping (*z*,*d*)→(*x*,*y*,*z*). In other words, we must specify *x* and *y* as functions of *z* and *d*. The electronic supplement, available to this article, contains the equations and parameter values that define these mappings for our three fault surfaces, as well as additional figures.

Figure 1 shows the source and target regions for the tests reported in this article. For the helicoidal fault, the source region is the left half of the fault, and the target region is a vertical strip one‐element wide, extending from the top to the bottom of the fault. For the other two faults, the source region is the bottom half of the fault (as measured along dip), and the target region is a strip one element high, extending from the left edge to the right edge of the fault. In all cases, the target region is immediately adjacent to the source region when performing the propagation test and is 4.8 km away from the source region when performing the source and target tests.

### Fault Discretization Patterns

We construct fault elements in several different patterns, so that we can compare the performance of rectangular and triangular fault elements, and also compare the performance of different triangulation patterns. Each pattern consists of either rectangles or triangles, and is produced with element sizes of 300, 600, 1200, 2400, 4800, and 9600 m. The discretization patterns are illustrated in Figure 2. The electronic supplement contains detailed step‐by‐step procedures for producing these patterns.

For the negatively curved (helicoidal) surface, we use three discretization patterns called rectangle, triangle‐2, and triangle‐4. For the positively curved (ellipsoidal) surface, we use three discretization patterns called rectangle, triangle‐1, and triangle‐3. For the cylindrical surface, we use four discretization patterns called rectangle, triangle‐1, triangle‐3, and triangle‐4. All the rectangular discretizations produce rectangles with upper and lower edges that are horizontal, which is a prerequisite for using the Okada formulas.

## TEST RESULTS

Figures 3–5 show some of our test results. Additional test results are available in the electronic supplement. The graphs show how the accuracy of the stress calculation varies as the size of the fault elements is varied. Each graph is a log–log plot, in which the approximation error (*Q* metric), in percentage, is plotted against the number of fault elements. For the source and propagation tests, it is the number of elements in the source region. For the target test, it is the square of the number of elements in the target region, to account for the fact that the size of the target region varies with the element size. (On a logarithmic plot, it makes little difference if the number of elements is squared or not, because squaring just changes the horizontal scale by a factor of 2.) The number of fault elements is the right thing to use on the horizontal axis, because the number of fault elements is what determines the computational effort that a simulator must exert to calculate the stress transfer.

(The Green’s functions for triangular fault elements are much more complicated than the Green’s functions for rectangular fault elements so take much longer to evaluate. But the time required to evaluate the Green’s functions is not significant, because an earthquake simulator computes and stores the Green’s function values at the start of the simulation and, from then on, just uses the stored values. So the total computational effort is dominated by the effort required to perform stress transfer calculations using precomputed Green’s function values, which is determined by the number of elements.)

There are three or four curves on each plot. One curve shows results for rectangular fault elements. The other curves show results for triangular fault elements, with each curve corresponding to a different method for constructing the triangulation. Each plot shows results for either shear stress or normal stress, as stated in the figure caption.

In eight of the nine cases shown, the rectangular fault elements perform as well as or better than the triangular fault elements. The exception is Figure 3b (source test for normal stress on the helicoidal fault), in which the triangle‐4 pattern performs better than the rectangular elements, although the rectangular elements still perform better than the triangle‐2 pattern.

In Figure 3c (target test for shear stress on the helicoidal fault), the rectangular fault elements greatly outperform both triangle‐2 and triangle‐4. We attribute the poor performance of triangles to the fact that, on a negatively curved surface, the triangulation is slightly corrugated. That is, differently oriented triangles have somewhat different strike and dip angles, which generate errors when the stress tensor is resolved onto the triangular fault elements. In contrast, our algorithm for producing rectangular fault elements attempts to preserve the strike and dip angles, which leads to improved accuracy in resolving the stress tensor.

An unexpected result is that sometimes different triangulation patterns can produce different levels of accuracy. Figure 3a,b,d shows cases in which the triangle‐4 pattern performs better than the triangle‐2 pattern. In contrast, Figures 4a–c and 5a,b show that the triangle‐1 and triangle‐3 patterns yield almost identical performance. Figure 5b (propagation test for shear stress on the cylindrical fault) shows a case in which triangle‐4 performs worse than triangle‐1 or triangle‐3.

Figures 3b, 4b, and 5a show normal stress, whereas the others show shear stress. Normal stress has larger errors than shear stress. This is due to the fact that, for source and target elements on the same fault, normal stresses tend to be smaller than shear stresses, so a small absolute error in normal stress can equate to a large percentage error. Bear in mind that when source and target elements are coplanar, the normal stress is zero except for free‐surface effects (and is exactly zero if the elements are vertical). So normal stress percentage errors are particularly large for the propagation test, where the source and target elements are adjacent and so are nearly coplanar.

## CONCLUSIONS

We performed accuracy tests using both rectangular and triangular fault elements, on three curved fault surfaces: a negatively curved fault in the shape of a helicoid, a positively curved fault in the shape of an ellipsoid, and a fault which has the cylindrical shape of the NFWF. We performed three kinds of tests: a source test that measures the errors when a fault element acts as a dislocation source, a target test that measures the errors when a fault element acts as a target, and a propagation test that measures the errors in propagating a rupture from a slipping region into the immediately adjacent fault elements.

In eight of the nine test results presented in this article, rectangles performed as well as or better than triangles. Including the results presented in the electronic supplement, this totals to 23 out of 28 tests. We also observed that one triangulation may perform better than another in an unpredictable way. For example, among the triangulations tested, triangle‐4 is the best on the negatively curved fault but the worst on the cylindrical fault. Although beyond the scope of this article, a systematic study of the performance of a larger variety of triangle configurations is warranted by our results.

Although either rectangles or triangles may be superior in a particular situation, we conclude that, in an earthquake simulator, overall, rectangles should perform as well as triangles. So the decision of whether to use rectangles or triangles need not be driven by concerns about accuracy, but instead should be based on other factors such as the capabilities of the simulator code and the discretization used in available fault models.

## DATA AND RESOURCES

Data for the North Frontal West fault geometry was obtained from the third Uniform California Earthquake Rupture Forecast (UCERF3) fault model v.3.2 (http://www.wgcep.org/components-deformation_model_3x, last accessed November 2015). All figures were made using Mathematica (http://www.wolfram.com/mathematica, last accessed November 2015).

## ACKNOWLEDGMENTS

This research was supported by the Southern California Earthquake Center (Contribution Number 6015) and by National Science Foundation Award EAR‐1135455. Southern California Earthquake Center (SCEC) is funded by National Science Foundation (NSF) Cooperative Agreement EAR‐1033462 and U.S. Geological Survey Cooperative Agreement G12AC20038. Thank you to Editor‐in‐Chief Zhigang Peng, Associate Editor J. Douglas Zechar, and three anonymous reviewers for many helpful suggestions.