Evaluating models of Coulomb stress transfer: Is variable fault geometry important?

To determine the importance of receiver fault geometry in Coulomb stress calculations a new methodology is presented to model faults with variable geometry. Although most models use planar faults, it is known that these are inaccurate representations of faults observed in the field. The central Italian Apennines are chosen as a straightforward tectonic system with well‐exposed normal faults to investigate the effect of variable geometry. It is shown that the static Coulomb stress transfer is most sensitive to changes in strike of the receiver faults, rather than changes in dip and rake. Therefore, a novel methodology to generate strike‐variable faults composed of discrete rectangular patches is developed. Using the 2009 L'Aquila earthquake (Mw = 6.3) the calculated stress transferred to planar versus variable faults is assessed. The stress transferred to variable faults is sufficiently different when compared to planar cases to merit other earthquake sequences being reassessed with available fault geometry.


Introduction
Coulomb stress transfer [Harris and Simpson, 1992;Reasenberg and Simpson, 1992;King et al., 1994] is routinely calculated following large earthquakes (recent examples include Gorkha, M w = 7.8, 2015 [Diao et al., 2015]; L'Aquila, M w = 6.3 [Walters et al., 2009], and Emilia-Romagna, M w = 6.1, 2012 [Saraò and Peruzza, 2012]). In each of these examples, the Coulomb stress is calculated from a planar source fault and resolved onto planar receiver faults. However, the fault maps in these areas [Decelles et al., 2001;U.S. Geological Survey, and California Geological Survey, 2006;Basili et al., 2008] show that these surface fault traces are not straight; therefore, a planar model is an inaccurate representation of the fault geometry at depth. Efforts to include variable fault geometry have been made prior to this study, including using triangular dislocation elements [Parsons et al., 1999;Meade, 2007;Bie and Ryder, 2014]. In this paper we investigate the effect of including variable source and receiver fault geometry (as rectangular elements) on the Coulomb stress calculation by using Coulomb 3.4 [Lin and Stein, 2004;Toda et al., 2005].
Coulomb stress transfer may be calculated either by resolving onto receiver faults of a known location and orientation [e.g., Toda et al., 1998], onto generalized faults of a specified orientation [e.g., Lin and Stein, 2004, Figure 13] or onto optimally orientated planes [e.g., King et al., 1994]. In this paper, we only consider the first approach in order to compare the difference between planar and strike-variable faults. It has long been established that when performing stress calculations the mean strike and dip of the receiver faults are important. However, what has not been routinely included is the use of variable strikes and dips along a single fault within Coulomb modeling. This is demonstrated in Lin and Stein [2004] for a strike-slip earthquake and two types of receiver faults, strike-slip and thrust faults, in different orientations. This example uses planar faults, despite the published fault traces being variable: the Coalinga anticline [Namson and Davis, 1988], White Wolf reverse fault [Stein and Thatcher, 1981], Garlock strike-slip fault [Dokka and Travis, 1990] Jiang et al. [2013] and Bie and Ryder [2014]). These examples use a triangular mesh approach outlined in Meade [2007]. This approach has not been widely adopted, and a simple rectangular fault, even discretized into numerous small patches, is still the standard approach to calculate the Coulomb stress transfer. This paper presents a methodology to model variable fault strike that has been developed for rapid use within Coulomb and has been tested rigorously to ensure that the results calculated are robust. In order to assess the importance of fault geometry, the Italian Apennines is used as case study region. This area of continental crust is undergoing extension at~3mm/yr [Hunstad et al., 2003;Serpelloni et al., 2005;D'Agostino et al., 2009;Faure Walker et al., 2010] which is distributed across a soft-linked array of normal faults [Cowie and Roberts, 2001]. No other faulting styles are widely recognized in the field with Holocene activity. Hence, the central Apennines is a relatively uncomplicated tectonic setting with one dominant style of faulting, and we are able to assess the effect of variable geometry on Coulomb stress transfer without having to consider different faulting styles. Throughout the central Apennines, normal faults are exposed at the surface as limestone bedrock scarps. These have formed over numerous earthquake cycles since the end of the Last Glacial Maximum (LGM, 15 AE 3ka) due to a decrease in erosion rates relative to throw rates since the LGM [Roberts and Michetti, 2004;Tucker et al., 2011]. These scarps can be mapped, and the changes in the geometry (strike, dip, and rake) and throw can be measured from meter to kilometer resolution [Boncio et al., 1996;Piccardi et al., 1999;Bonini et al., 2003;Roberts and Michetti, 2004;Papanikolaou et al., 2005;Roberts, 2007;Faure Walker et al., 2009;Wilkinson et al., 2015;Mildon et al., 2016]. Hence, the geometry of the normal faults can be well constrained by using the gathered field data, and the effects of including variable fault geometry into Coulomb stress calculations can be assessed. The 2009 L'Aquila (M w = 6.3) earthquake is used as an example to calculate the static Coulomb stress transfer onto the surrounding faults and to compare the planar and variable fault geometry approach.

Coulomb Stress Sensitivity to Receiver Fault Geometry
There are three parameters that are used to describe the geometry of a fault: strike, dip, and rake (slip direction). Sensitivity tests have been performed for each parameter independently to identify the spatial sensitivity of Coulomb stress transfer around a source fault.
A simple pure dip-slip normal source fault is used for the test (180°/65°/À90°for the strike/dip/rake following the convention of Aki and Richards [1980], 10 km long with 2 m of uniform slip), and the Coulomb stress is calculated at a depth of 5 km. The coefficient of apparent friction used μ 0 [King et al., 1994] is 0.4. For the receiver faults, each geometry parameter is varied individually through a range of values (strike: 0 to 360°, dip: 30 to 80°, rake: À60 to À120°), while the other two parameters are kept constant and equal to the source fault. The limits of each parameter are based on the range of values seen in the Apennines for normal faults. The results of each test are shown in Figure 1. Other depths and values of μ 0 have been tested and are shown in

Methodology
Surface fault traces of normal faults in the central Apennines are drawn in Google Earth TM based on fieldwork (>10,000 measurements at >800 sites [from Morewood and Roberts, 2000;Roberts and Michetti, 2004;Papanikolaou et al., 2005;Papanikolaou and Roberts, 2007;Roberts, 2008;Faure Walker et al., 2010Wilkinson et al., 2015;Mildon et al., 2016], satellite imagery, and geological maps). The field measurements are used to check that the strike of the surface traces drawn in Google Earth TM is a good representation of the true fault strike (for faults where there is limited field data coverage); there is good agreement between the traces and the field data.
From the fault traces drawn in Google Earth TM , a projection direction is assigned to each fault; this is the direction perpendicular to the mean fault strike (because the mean rake is À90°for a normal fault) and controls the dip direction of the fault. For each fault, a single value of dip is used which is either the mean dip measured from field measurements (where available) or taken from the literature. If no information is available, a dip of 65°is assigned as the mean dip from all data [Roberts and Michetti, 2004]. A single value of rake is used (À90°), despite variations present in the field data. The variability in dip and rake is not modeled in full in order to isolate the effects of including the strike variability and because the sensitivity tests ( Figure 1) showed that the transferred Coulomb stress has low sensitivity to the changes in the dip and rake of the receiver faults. In addition, changes in strike can be resolved readily from geological maps; hence, the importance of variations in this parameter for other fault systems can be understood and easily applied. The fault traces are divided into segments of a specified length of 1 km, and the points that divide the segments are projected to depth to create rectangular elements to be used as the input for Coulomb. The specified length of 1 km is chosen based on a trade-off between the level of detail and computational time. If the grid size were changed, this would only change the resolution of the stress field calculated. The faults are projected to a depth of 15 km by using a constant value of dip, the assumed depth of the brittle-ductile transition based on heat-flow modeling [Cowie et al., 2013], and the termination of earthquakes [e.g., Valoroso et al., 2013]. This process is done for all faults in the central Apennines in the Marche, Umbria, Abruzzo, Lazio, and Molise regions, by gridding the surface traces to a size of 1 km. Faults that are shorter than the depth of the seismogenic zone are assumed to maintain an aspect ratio of 1. This produces a model of the region shown in Figure 2 which contains >17,000 individual elements. Regional stress [King et al., 1994] is not assigned for these models because it is unnecessary when the stress is being resolved onto known receiver faults. Changing the size of the elements will affect the area of the stress discrepancies that can be resolved, and decreasing the element size will rapidly increase the computation time. This code is written to generate input files for immediate use within Coulomb and is available on request.
The L'Aquila 2009 earthquake is used as a case study to test the importance of planar versus strike-variable receiver fault geometry. There are numerous slip distributions published for this earthquake, see Wilkinson et al. [2015] for a summary. The most prominent feature in all these slip distributions is that the highest slip occurs close to but southeast of the center of the fault plane and the slip forms a bulls-eye pattern around this point of maximum slip. All of these slip distributions are calculated for a planar source fault; this is not comparable to the strike-variable source fault used in the models presented in this paper. Hence, a simplified slip distribution of maximum slip skewed slightly to the southeast of the strike-variable fault plane, with zero slip at the ends and at 15 km depth is used here (Figures 2c and 3). The value of maximum slip was chosen to fit the moment magnitude of the event.
It is not geometrically possible to grid a smooth surface with variable strike using only rectangular elements (as allowed in Coulomb). Therefore, it is important to understand the effects of underlap and overlap between the elements on the generated stress field. Tests were undertaken to investigate the possible effects ( Figure S1). It was found that any spurious stress signals are confined to the region directly above the source fault (when the Coulomb stress is calculated for a specified depth) and does not affect the Coulomb stress resolved directly onto receiver faults. Hence, these effects do not affect the results of this paper.

Straight Versus Variable Fault Geometry
The results of Coulomb stress transfer for planar versus strike-variable faults are shown in Figure 3 for the L'Aquila 2009 earthquake. For the planar faults, the transferred stress is largely as expected for a normal fault, with receiver faults along strike from the source fault generally experiencing positive stress transfer and receiver faults across strike from the source fault experiencing negative stress. When the faults are modeled with strike-variable geometry, a similar pattern is evident, but there are notable differences between the planar and strike-variable cases. The difference in stress between the strike-variable and planar models is shown in Figure 3c. There are two main features that can be identified. The faults that show the most difference between the planar and strike-variable cases are located in the across strike from the source fault, in the footwall, and hanging wall. The areas on individual faults that show the largest differences in stress are located where there are prominent bends in the fault trace. Some examples have been highlighted in Figure 3c such as grid references 370/4670 (the tip of the Fucino fault), 395/4697 (Campo Imperatore East fault), and 385/4685 (Barisciano fault). The Velino Magnola fault has a smooth bend along its full length, and differences in stress are seen along its full length in Figure 3c. This indicates that the stress on these sections of the faults would be incorrectly resolved if modeled by the standard planar approach. Differences between the planar and strike-variable geometry cases are observed across most of Figure 3c; this shows that this is a problem across the majority of the receiver faults, especially close to the source fault and across strike.

Sensitivity Tests
The sensitivity tests shown in Figures 1 and S2-S4 show that of the parameters strike, dip, and rake, Coulomb stress transfer is most sensitive to changes in the strike of the receiver faults when the strike, dip, and rake are considered independently. Almost all regions around the source show that the stress can switch from positive to negative with changes in strike of 10-20°, except for the region directly in the hanging wall of the source fault, which is insensitive to strike. The areas that are 45°into the footwall from the tips of the fault (receiver faults in these regions would be described as being en echelon to the source fault) show especially high sensitivity where the stress switches from positive to negative around the strike of the source fault. This is important because many published models resolve the stress onto faults of exactly the same orientation as the source fault. If the actual receiver faults varied in strike by as little as 10°, then the Coulomb stress will be incorrectly calculated.

Discussion
The hypothesis that the Coulomb stress transferred to receiver faults is highly dependant on geometry has been reassessed and improved. Based on sensitivity tests, it is shown that different areas around the source fault are particularly sensitive to changes in different parameters of geometry. It has been shown that there are particular regions where the receiver fault geometry should either be well known or the Coulomb stress transferred should be calculated for a range of geometries. One limitation is that this method cannot be used robustly to investigate the localized stress field around the source fault because spurious stress signals may be generated. Previous examples of Coulomb stress calculations should be re-examined with this receiver fault sensitivity in mind.
Based on the high sensitivity of Coulomb stress to receiver strike, an approach to modeling normal faults with strike-variable geometry has been outlined, where the fault planes are modeled as a series of discrete rectangular patches. Using the central Italian Apennines and the L'Aquila 2009 earthquake as a case study, it has been shown that the stress resolved onto strike-variable geometry receiver faults is considerably different to the stress resolved onto discretized segmented planar faults (Figure 3c). This will have implications for the assessment of the stress transferred to the surrounding faults and the distribution of aftershocks following a major earthquake.
Although the results here are shown for a normal faulting earthquake in a relatively simple tectonic setting, similar effects are highly likely to be seen for any faulting style when the strike-variable geometry is included. Hence, any Coulomb modeling performed for any faulting style may be missing parts of the transferred stress pattern if strike-variable fault geometry is not considered where possible. Few places show such excellent fault exposure as the Italian Apennines, but this should not be considered a barrier to modeling other earthquakes and faults using this methodology. Existing fault maps from the literature, geological maps or satellite imagery can be utilized to generate variable fault planes by using the methodology outlined in this paper. The static Coulomb stress transfer of previous earthquakes and sequences should be re-examined by using this approach wherever variations in fault geometries are available.

Conclusions
A new methodology has been proposed to model fault planes with variable fault strike as planes composed of rectangular elements in Coulomb 3.4. This method has been tested rigorously and produces interesting and valid results and is available on request. This method is a powerful new approach to studying the stress transferred to receiver faults with known variable strike along their lengths. There are large discrepancies observed between a planar and a variable geometry model of faults in the Italian Apennines, with the 2009 L'Aquila earthquake shown as a case study. These discrepancies show that the transferred stress will not be uniform across a fault and regions of high stress may be generated on faults with a low average stress. These differences may better explain the locations and distributions of aftershocks and the probability and locality of the next major earthquake. Previous Coulomb stress calculations should be re-examined taking account of the receiver fault geometry wherever possible because, as demonstrated herein, the geometry has a great effect on the Coulomb stress transferred to differently orientated sections.