A 667 year record of coseismic and interseismic Coulomb stress changes in central Italy reveals the role of fault interaction in controlling irregular earthquake recurrence intervals

Current studies of fault interaction lack suf ﬁ ciently long earthquake records and measurements of fault slip rates over multiple seismic cycles to fully investigate the effects of interseismic loading and coseismic stress changes on the surrounding fault network. We model elastic interactions between 97 faults from 30 earthquakes since 1349 A.D. in central Italy to investigate the relative importance of co-seismic stress changes versus interseismic stress accumulation for earthquake occurrence and fault interaction. This region has an exceptionally long, 667 year record of historical earthquakes and detailed constraints on the locations and slip rates of its active normal faults. Of 21 earthquakes since 1654, 20 events occurred on faults where combined coseismic and interseismic loading stresses were positive even though ~20% of all faults are in “ stress shadows ” at any one time. Furthermore, the Coulomb stress on the faults that experience earthquakes is statistically different from a random sequence of earthquakes in the region. We show how coseismic Coulomb stress changes can alter earthquake interevent times by ~10 3 years, and fault length controls the intensity of this effect. Static Coulomb stress changes cause greater interevent perturbations on shorter faults in areas characterized by lower strain (or slip) rates. The exceptional duration and number of earthquakes we model enable us to demonstrate the importance of combining long earthquake records with detailed knowledge of fault geometries, slip rates, and kinematics to understand the impact of stress changes in complex networks


Introduction
Advances in high-resolution topography and Quaternary dating tools suggest that transient fault slip rates and irregular earthquake recurrence intervals over thousand year timescales are a common feature of extensional areas [e.g., Friedrich et al., 2003;Palumbo et al., 2004;Bull et al., 2006;Nicol et al., 2010;Begg and Mouslopoulou, 2010;Akçar et al., 2012;Jewell and Bruhn, 2013;Nicol et al., 2016;Cowie et al., 2017] as well as other tectonic environments [e.g., Weldon et al., 2004;Oskin et al., 2008;Dolan et al., 2016;Gold et al., 2017]. A number of different processes have been invoked to explain this transient behavior including static elastic stress transfer [Benedetti et al., 2013], dynamic (i.e., coseismic) stress changes [Brodsky and van der Elst, 2014], temporal variations in the strength of brittle faults and ductile shear zones [Dolan et al., 2007], fluid migration through fault zones [Oskin et al., 2008], interactions with surface processes [Hetzel and Hampel, 2005], episodic release of strain stored in a crustal stress "battery" [Gold et al., 2017], and energy dissipation and minimization of work in response to flexural bending of normal fault footwalls [Cowie et al., 2017]. However, the timescales upon which these processes can affect fault slip rates are poorly constrained, and it is not clear which of these processes, if any, are the dominant mechanism controlling transient fault slip rates and irregular earthquake recurrence times. In this paper, we examine the role of elastic Coulomb stress interactions in the variability of fault slip rates and earthquake recurrence times over timescales of 10 3 years.
Elastic interactions between active faults have been invoked to explain a range of observations of fault and earthquake behavior within distributed fault networks including temporal and spatial changes in slip rates over periods of 10 2 -10 5 years [Cowie et al., 2005;Nicol et al., 2010] and triggering of subsequent earthquakes . Numerical modeling of this process is often used to explore the effects of stress interactions WEDMORE ET AL.
FAULT STRESS INTERACTIONS, CENTRAL ITALY 5691 on earthquake clustering, variation in recurrence intervals, and time varying fault slip rates [Robinson, 2004;Marzocchi et al., 2009;Robinson et al., 2009;Cowie et al., 2012]. However, many of these models are purely theoretical [e.g., Cowie et al., 1993;Zöller and Hainzl, 2007] or include only a relatively short record of recent earthquakes [e.g., Stein et al., 1992;King et al., 1994]. Other models lack many of the structural characteristics of fault networks that are thought to be key factors in determining how faults interact such as detailed fault geometry and kinematics, variable slip rates and interseismic loading rates for different faults, and earthquake histories that contain multiple seismic cycles [e.g., Robinson, 2004;Marzocchi et al., 2009;Robinson et al., 2009;Marzocchi and Melini, 2014]. Consequently, the characteristics of how faults interact over periods of time longer than individual earthquakes and what affect these longer-term interactions have on current earthquake activity, fault behavior, and seismic hazard are poorly understood.
Over timescales of 10 4 -10 5 years fault interactions lead to an increase in length and displacement rates on faults that are optimally located and orientated within a diffuse network of extensional faults [Cowie, 1998]. Over shorter timescales (10 2 -10 3 years) the effects of fault interactions on fault slip rates and earthquake recurrence intervals are less well understood, but it appears that earthquake activity is clustered in both time and space. Earthquake clusters seen in high-resolution dated stratigraphic sequences in actively extending regions suggest that fault slip rates vary over timescales of 10 2 -10 3 years [e.g., Bull et al., 2006;Nicol et al., 2010]. Evidence from 36 Cl cosmogenic exposure dating of bedrock fault scarps produced by normal faulting earthquakes show accelerated fault slip rates over periods of a few thousand years that are presumably associated with temporal clusters of earthquakes, alongside periods of quiescence (e.  Stein et al., 1997], the role of Coulomb stress changes in causing the sort of regional changes in seismicity over 10 2 -10 3 years observed in 36 Cl cosmogenic exposure data and highresolution stratigraphy is less clear. We use the central Apennines, Italy, to explore how faults are interacting through elastic Coulomb stress perturbations. Spatially transient temporal clusters of seismicity in the region over 10 2 -10 3 years ( Figure 1) [Schlagenhauf et al., 2011;Benedetti et al., 2013;Tesson et al., 2016], and shorter-term clusters of earthquakes (e.g., the 2016 central Italy earthquake sequence) suggest that faults are interacting over a range of timescales. Furthermore, the region has exceptional constraints on (i) the long record of historical earthquakes ( Figure 1) which appears to be complete for events M > 5.8 since 1349 A.D. [Guidoboni et al., 2007], (ii) the constraints on the mechanism controlling fault slip rates over multiple seismic cycles where slip rates appear to be controlled by viscous deformation at depth in response to long wavelength topographic uplift [Faure Walker et al., 2012;Cowie et al., 2013], (iii) the extensive record of Holocene fault slip rates in the region where nearly every fault has a slip rate derived from offset geomorphic features dating from 15 ± 3 ka ( Figure 2a) [e.g., Faure Walker et al., 2010], and (iv) the extensive record of fault geometries and kinematics across faults in the region (Figure 2b) Faure Walker et al., 2009. We use these constraints and records of earthquake activity in the central Apennines to model interseismic and coseismic stress interactions to explore the characteristics of an actively interacting distributed fault network and to determine whether Coulomb stress changes can explain observations of transient fault slip rates and irregular earthquake recurrence intervals.

Tectonics of the Central Apennines
Extensional faulting in the central Apennines initiated~2-3 Ma as thrusting caused by subduction of the Adriatic plate beneath Italy slowed, and a slab window opened allowing mantle upwelling to become the dominant geodynamic process [D'Agostino et al., 2001;Faure Walker et al., 2012]. This has produced sets of subparallel extensional faults up to 40 km in length, distributed over the crest of a long-wavelength topographic high along the axis of the central Apennine mountain chain that has been uplifted from sea level to >1000 m elevation since the Pliocene [D'Agostino et al., 2001;Galadini et al., 2003a]. Geodetic extension rates across the central Apennines of 2.9 ± 1 mm yr À1 [D'Agostino et al., 2011] broadly match the~3 mm yr À1 maximum extension rate calculated across the region by summation of 15 ± 3 ka strain rates [Faure Walker et al., 2010]. Upper crustal brittle strain rates also correlate with areas of high topography, high free air

Faults in the Central Apennines
Individual normal faults in the central Apennines offset prerift Tertiary and Mesozoic carbonates by up to 2.5 km . Following initiation of extension in the Plio-Pleistocene, fault interaction and crustal-scale strain localization led to fault growth and caused an increase in throw rates on the centrally located faults from~0.7 Ma [Roberts et al., 2002;Whittaker et al., 2007]. Offsets of sediments and hillslopes that have been preserved since the Last Glacial Maximum (LGM) have been measured at 174 locations on 48 faults and yield fault throw rates in the central Apennines averaged over 15 ± 3 ka ( Figure 2a) Faure Walker et al., 2010;Tucker et al., 2011;Mildon et al., 2016]. Measurements of striated and corrugated fault planes at 278 locations (>14,000 individual measurements) constrain the kinematics of the individual faults ( Figure 2b) [Morewood and Roberts, 2000;Faure Walker et al., 2010;Mildon et al., 2016].

Historical Earthquakes
Since 9 September 1349 A.D., the day of three well-document large magnitude earthquakes, the historical catalogue of earthquakes in central Italy, which comprises of macroseismic shaking intensities for damaged towns and villages, is thought to be complete for all earthquakes with M e ≥ 5.8 [Michetti et al., 1996;Guidoboni et al., 2012]. From 1349 A.D., the historical records include 27 large-magnitude  Tables 1 and S1 in the supporting information), contemporary reports of surface ruptures [e.g., Blumetti, 1995;, seismology [e.g., Westaway et al., 1989], and geodetic records [e.g., Walters et al., 2009] to constrain the faults that ruptured in 30 major earthquakes since 1349 A.D. (Table 1).
Projecting the macroseismic shaking from historical earthquakes since 1349 A.D. onto a transect perpendicular to the strike of faults in the central Apennines shows that the majority of the activity during this time has  [Guidoboni et al., 2007]. The inferred causative fault of each of the events and the literature sources for this inference are also listed. b The magnitudes of earthquakes that occurred prior to 1979 A.D. are listed in the Catalogue di Forti Terremoti [Guidoboni et al., 2007] and are derived from the macroseismic shaking records using the Boxer code of Gasperini and Ferrari [2000] and as such are described as equivalent magnitudes (M e ). For more recent earthquakes (1979-2016 A. D.) the magnitudes described are from seismological sources. c These earthquakes occurred at the extreme north-west and south-east of the study region. They were included to ensure that faults are not affected by edge effects in the modeling of fault interactions; however, they are not included within the analysis in Figure Walker et al., 2010Walker et al., , 2012 are more evenly distributed across both sides of the central Apennines, across the crest of the regional long wavelength topographic bulge [Cowie et al., 2017]. To investigate the role of elastic fault interaction in causing this skewed regional seismicity since 1349, we model the historical earthquakes and interseismic loading of faults from 1349 to 2016 A.D. to investigate how faults are interacting in the central Apennines and investigate the effects of these interactions on perturbing earthquake interevent times.

Model Setup
We are interested in determining the elastic stress perturbations caused by earthquakes on surrounding faults in a complex normal fault network due to both coseismic static stress changes and interseismic loading.
For coseismic stress changes we calculate the change in Coulomb stress (ΔCFF) following each historical earthquake in the central Apennine fault network. The change in Coulomb stress is defined as where Δτ is the change in shear stress on a fault, Δσ n is the change in normal stress on a fault, and μ 0 is the effective coefficient of friction. Shear stress is resolved in the slip direction of each fault, and normal stresses are positive if that fault becomes unclamped. We use a value for μ 0 of 0.4 although our results are not sensitive to the value used within the standard range of 0.2-0.6. The ΔCFF calculations were performed using the Coulomb 3.3 code [Lin and Stein, 2004;Toda et al., 2005] that uses the elastic half-space code of Okada [1992].
We model 97 faults in the central Apennines that show evidence of activity since the demise of the LGM as receiver faults. The fault map is constructed from geological maps of the region with evidence for Holocene activity from field mapping of post-LGM geomorphic offsets [Morewood and Roberts, 2000;Cowie and Roberts, 2001;Roberts et al., 2002;Papanikolaou et al., 2005;Faure Walker et al., 2010Wilkinson et al., 2015]. This probably includes all major faults (i.e., with slip rate >~0.1 mm yr À1 ) in the region as (a) slip rates as low as 0.1 mm yr À1 would produce cumulative offsets of 100 m over 1 Ma that would be resolvable on geological maps, (b) if significant contributors to the overall strain were missing then correlations between geodetic and geologic measurements of strain and between strain rate and topography would not hold [e.g., D' Agostino et al., 2011;Cowie et al., 2013], and (c) no significant clusters of microseismicity have occurred outside areas surrounding mapped faults during the last 30 years [Boncio et al., 2009].

Modeling Faults, Fault Geometries, and Fault Kinematics
Faults are modeled as linear features within the elastic half-space model and partitioned into 3 km long sections along strike with five sections downdip to a vertical depth of 15 km, constrained by the maximum depth of seismicity measured in the region [Boncio et al., 2009;Chiarabba et al., 2009]. Three-kilometer-long sections were chosen as this is the approximate resolution along strike at which kinematic and throw-rate measurements have been collected along the faults (Figure 2a). We include the fault kinematics as they can alter the local stress field toward the tips of faults [Maniatis and Hampel, 2008]. As fault strike is a controlling factor on fault kinematics, we also correct for the differences between the strike of the mapped faults and the planar fault trace by preserving the principal direction of strain and the magnitude of strain rate (Text S1 and Aki and Richards [1980], Faure Walker et al. [2009], Kostrov [1974], Wilkinson et al. [2015]).

Interseismic Stress Accumulation
In the central Apennines, Cowie et al. [2013] showed that below the seismogenic zone, nonlinear viscous deformation dominates and promotes strain localization. Consequently, we load the faults interseismically by continuous creep on shear zones at the base of each fault, below the locked brittle (seismogenic) upper crust ( Figure 3). The agreement between modern GPS extension rates and summed 15 ± 3 ka strain rates measured across all faults in the central Apennines indicates that there is no overall stress accumulation on the faults at geological timescales [Faure Walker et al., 2010]. Therefore, we impose the slip rates measured at the surface from the offset of 15 ± 3 ka landforms (Figure 2a) onto our modeled shear zones (Figure 3b). We then resolve the stress changes caused by the creeping shear zones onto patches in the upper seismogenic sections of the faults within the elastic half-space model described above (Figures 3b and 3c). The advantages of using this method to load the faults are as follows: (i) fault-specific slip rates can be used as a constraint on the long-term loading rate in such a way that would not be possible if bulk measures of regional strain (e.g., GPS measurements) were used; (ii) the approach relates well to physical models that show that in the presence of a brittle-ductile transition downdip, the viscous layer deforms continuously whereas the brittle layer remains locked between earthquakes [e.g., Huc et al., 1998]; and (iii) the setup creates an increase in stressing rate with depth such that the whole fault would reach its failure stress at about the same time. An alternative approach to load the brittle fault uniformly downdip would result in disproportionately high levels of stress (and therefore unstable slip) near to the surface not seen in observations of the depth of large earthquakes [e.g., Das and Scholz, 1983]. A depth of 15-24 km is used for the creeping shear zones as above this depth deformation occurs through brittle deformation as evidenced by active seismicity and the lower extent falls within the depth range found by Cowie et al. [2013] for localized viscous flow.
Throw measurements over 15 ± 3 ka were converted to slip rates per year for each section along strike using the dip of the modeled faults as calculated above (see Text S1 for details of this calculation and Table S2 for fault-specific parameters). Slip rate within the shear zone was assumed to be constant downdip and zero at the fault tips. The fault kinematics were assumed to match the kinematics of the portion of the fault in the upper crust. The resultant stress transferred reflects the annual interseismic loading stress acting on the fault. An internal check to ensure that the stress loading rates are consistent with earthquake recurrence intervals observed in palaeoseismic trenches is shown in the supporting information (section S2 and Cinti et al. [1992Cinti et al. [ , 2011, Field et al. [1999], Galli [1999, 2003], Galadini et al. [1997Galadini et al. [ , 2003b, Galli et al. [2002Galli et al. [ , 2005Galli et al. [ , 2011, Giraudi and Frezzotti [1995], Michetti et al. [1996], Moro et al. [2002], Pace et al. [2006], Pantosti et al. [1996], Salvi et al. [2003], Scholz [1994Scholz [ , 2002).

Modeling Coseismic Stress Changes From Historical Earthquakes
Although the time and location of the historical earthquakes are relatively well known, the slip distributions are poorly constrained by measurements with debate even surrounding the well-studied L'Aquila earthquake in 2009 (see Wilkinson et al. [2015] for a summary). Although some earthquake slip distributions are asymmetrical [Manighetti et al., 2005], we use triangular-shaped slip distributions to model the historical earthquakes, as normalized finite displacement profiles in the central Apennines, which average the earthquake surface slip distributions over the last 3 Ma, show a more symmetrical, triangular-shaped displacement profile (Figure 2c). We calculate average slip across the triangular-shaped slip distribution from empirical relationships between magnitude and average slip [Wells and Coppersmith, 1994] with the maximum slip centered in the middle of the faults (following Nalbant et al. [2013]). Rupture length was calculated using empirical relationships between magnitude and subsurface rupture length [Wells and Coppersmith, 1994]. Where palaeoseismic evidence and contemporary reports suggest multiple fault ruptures (e.g., 14 January 1703 [Blumetti, 1995], 2 February 1703 [Galli et al., 2010], 13 January 1915 [Michetti et al., 1996], and 24 August 2016 [Livio et al., 2016]), the empirically calculated rupture lengths were divided among the faults responsible for the earthquake (Figure 2b).
The stress changes on faults at any point in time are the sum of the interseismic and coseismic stress changes up to that point, starting from immediately prior to the first of the earthquakes in 1349 A.D.. Stress changes were resolved on all faults in the region (all faults act as receiver faults). Coulomb stress on the portions of the faults that ruptured in each event were reset to zero.

Errors and Uncertainties
We apply errors to the magnitudes of earthquakes prior to 1979 as these magnitudes are derived from inversions of macroseismic shaking records [Gasperini and Ferrari, 2000]. were run varying the source fault of each event, as two possible causative faults were indistinguishable given the palaeoseismology and historical macroseismic records available. We included the uncertainty of the time constraint on the fault slip rate measurements by running interseismic loading models with slip rates averaged over 12 ka, 15 ka, and 18 ka, respectively. This results in 36 end-member modeling scenarios by varying (i) the slip rates used for interseismic loading (three different scenarios), (ii) the faults that ruptured in 1456 and 1654 (four different scenarios), and (iii) the magnitude and average slip of earthquakes (three different scenarios). These alternative results are shown in Figure S5 in the supporting information and summarized below.
There is uncertainty over whether mean or maximum stress is the best indicator of future large earthquakes (see the debate between the cascade and deterministic models of rupture propagation Allen, 2005, 2006;Rydelek and Horiuchi, 2006]). We base our decision to analyze the mean Coulomb stress changes on need to account for the stress across the whole fault including the inherent heterogeneous stress distribution created by previous ruptures, closely spaced faults, and complex structural barriers [Aki, 1984;Steacy et al., 2004;Bull et al., 2006;Noda et al., 2013]. The implications of using the mean stress in our analysis is that it will (i) produce a more generalized view of the long-term behavior of the fault system, (ii) give an impression of the relative likelihood of future major earthquakes rather than small events, (iii) mean our results are not overly biased by the slip distribution of the events modeled, and (iv) allow our results to be affected by the realistic structural and kinematic variability we have incorporated in our model.

Results
The maximum interseismic Coulomb stress loading rate on any individual fault patch in the upper crust is shown in Figure 3d. Due to the model setup, these loading rates are broadly consistent with the loading rates that may be inferred from palaeoseismic recurrence times (see the supporting information).  Figure 7a where the mean Coulomb stress on each fault prior to each earthquake is shown. The mean accumulated Coulomb stresses on the faults that subsequently rupture (red circles) are positive and between 0 and 1 MPa for all but one earthquake: the M = 6.0 earthquake in 1933 on the Maiella Fault. However, 5 of the 15 patches of the Maiella Fault that ruptured were experiencing positive stress. Thus, a pattern emerges where earthquakes occur preferentially on faults that remain positively stressed after hundreds of years of interseismic loading and coseismic loading and unloading, even though 20% of faults are in "stress shadows" at any one time. To test whether this pattern is significant, we calculated the percentile of the stress on the fault that failed in each event relative to the distribution of stresses on each fault. This analysis was performed for each earthquake from 1654 A.D. onward, and the mean percentile of the fault that failed was calculated. Faults were then randomly selected from the distribution of stresses on all faults for each earthquake from 1654 A.D. onward, and the average percentile for this random sequence was then calculated. This random selection was made 10,000 times, and the resultant distribution was compared with the mean percentile of failure stress of the actual sequence of events. A nonparametric Mann-Whitney U test was performed to test whether the actual sequence of events was statistically different from the 10,000 random draws. For the modeling scenario shown in Figures 4, 5, and 7a, the value of U was 0.039, indicating that the actual sequence of earthquakes is significantly different from a random sequence of earthquakes at a 95% confidence level. In essence, earthquakes are more likely to occur on faults above the mean Coulomb stress acting on all faults. Crucially, this result holds not just for one earthquake but for 21 earthquakes since 1654. Figure 7b shows the change in the number of faults through time where the combined accumulated coseismic and interseismic Coulomb stresses are positive or negative. When summed over time, interseismic stresses (~10 À2 -10 À3 MPa yr À1 ) might be expected to dominate the final stress map ( Figure 6) and to overcome any (negative) coseismic stress changes (approximately À10 À1 MPa) after 10-100 years. However, after running our model for 667 years, the coseismic stress changes have not been overcome, with approximately 20% of faults (20-25 out of 97 faults) being net negatively stressed at any time since 1654 A.D.. As these negative  stresses are only generated within coseismic stress shadows, we investigate the effects that the coseismic stress changes for all 30 earthquakes could have on fault slip and earthquake activity in the longer term.
We assessed the impact of the combined coseismic stress changes on the recurrence interval of each fault by calculating the notional advance or delay time of future earthquakes on each of the faults from the 30 earthquakes modeled from 1349 to 2016 A.D.. By analyzing the combined coseismic stresses for all 30 earthquakes alongside the total accumulated interseismic stresses from 1349 to 2016 (Figure 8a) we calculated the time advance or delay to the next earthquake, ΔT, at the end of the earthquake sequence using the equation: where ΔCFF n is the mean coseismic stress change on a fault following earthquake n and CFF i is the interseismic loading rate. The advance or delay of future earthquakes due to the coseismic stress changes for all  Figures 8b and 9. The summed effect of the coseismic stress changes for 30 earthquakes can advance or delay the expected timings of earthquakes by thousands of years ( Figure 9).
To investigate potential structural controls on the impact of fault interactions on short timescales, we analyzed whether fault length or fault orientation affects the advance or delay of future earthquakes (Figures 10 and S4). Figure 10 shows how shorter faults are more prone to advances or delays in future earthquakes times than longer faults. This appears to follow the relationship between fault loading rates and fault length that shows that longer faults have faster interseismic loading rates (Figure 3e). Moreover, the standard deviation of advance or delay times is also greater for shorter faults: recurrence intervals for shorter faults are more variable on both an individual fault basis and across the whole fault system. In contrast, there appears to be no systematic relationship between overall fault orientation and the advance or delay of future earthquakes in our model ( Figure S4).
The results shown in  Figure S5. While section 4 above and discussion below refer to the results and figures shown in the main body of this paper, the discussions and conclusions hold irrespective of the iteration of the earthquake location, fault slip rate, or historical earthquake magnitude variations used (Figures 10d and S5).

Discussion
We show that earthquake interevent times on 97 individual faults become highly variable over the course of 30 earthquakes across the region, despite imposing linear interseismic loading. We show, using the exceptional record of historical earthquakes in central Italy, that the combination of interseismic loading over 667 years and the coseismic Coulomb stress changes from 30 earthquakes produces a statistically significant correspondence between the location of previous earthquakes and subsequent seismic events. We demonstrate the importance of coseismic stress shadows on perturbing the earthquake recurrence, particularly, of short faults. We argue that our results show how fault length plays a key role in the effects of fault interaction. Of the 21 earthquakes from 1654, all but one event occurred on faults where the mean combined interseismic and coseismic stresses were positive (Figure 7). This is despite 21-26% of faults being negatively stressed at any one time. This positive result for 20 events confirms the strong association between positive Coulomb stress change and earthquake occurrence observed by other authors [e.g., Deng and Sykes, 1997;Stein et al., 1997]. Furthermore, the use of a statistical Mann-Whitney U test enables us to reject the (null) hypothesis that the stresses on the faults that failed fall within a random pattern within the distribution of stresses on all the faults.
Our approach of imposing continuous creep in deep shear zones to set long-term fault slip rates (Figures 3a-3c) leads to some faults having higher interseismic loading rates (Figure 3d). The two longest faults have the highest loading rates because these faults have the highest measured Holocene throws (and thus slip rates) across them and many of the shorter faults are loaded at lower rates (although the correlation is not strong R 2 = 0.25; Figure 3e). While this result is due to our model setup (Figure 3), a physical explanation for this difference is that ductile deformation is likely to be more localized downdip of large faults and less localized downdip of smaller faults. Such behavior is typical in a material with a nonlinear viscous rheology as has been inferred for the central Apennines [Cowie et al., 2013].
Given that the coseismic stress changes are only 1-2 orders of magnitude greater than the interseismic loading rates of~10 À2 -10 À3 MPa yr À1 , it would be expected that after~10-100 years, negative coseismic stresses would be eclipsed by interseismic loading stresses. However, throughout the timescale modeled,~20% of the faults are consistently experiencing negative stresses (considering both interseismic and coseismic stress changes; Figures 6 and 7b). This may be due to the apparent cluster of earthquakes in the area surrounding the L'Aquila valley during the last 667 years. Earthquakes occurred in this area in 1461, 1703 (2 February), and 2009, and the concentration of relatively short faults across strike from the L'Aquila valley results in these faults having their stress decreased following each of the events (Figure 6). While the importance of stress shadows has been noted before [e.g., Harris and Simpson, 1998], by resolving stress changes on individual faults, we show the importance of stress shadows on perturbing the recurrence times of short faults with low interseismic loading rates (Figures 9 and 10). It takes longer for the interseismic loading stresses to overcome the negative coseismic stresses on shorter faults given their generally lower interseismic stress loading rates (Figures 10 and 3e). In addition, due to the comparative size of the stress shadows, it appears that negative coseismic Coulomb stresses have a greater impact on extending fault recurrence intervals compared with the shortening of fault recurrence intervals by positive coseismic Coulomb stresses ( Figure 10). Therefore, although the role of stress triggering in earthquake clustering over short timescales (i.e., days tõ 10 years) is well documented [e.g., Stein et al., 1994], it does not appear that Coulomb stress changes alone can cause the significant accelerations in fault slip rate over thousands of years revealed by 36 Cl data long bedrock fault scarps in the central Apennines [e.g., Cowie et al., 2017]. However, Coulomb stress changes are more effective at delaying future earthquakes by thousands of years, especially on shorter faults (Figures 9 and 10).
The sensitivity of small faults to slip rate perturbations has been observed in both theoretical and real fault systems. Palaeoseismic data from normal faults in the Taupo Rift, New Zealand, show that faults with lower slip rates experience greater variations in displacement rate through time [Nicol et al., 2010] and numerical simulations of fault interactions suggest that the percentage of across strike strain accommodated by a fault (which is analogous to fault slip rate) affects the slip rate variability . Figure 10 shows how the stress perturbations modeled have a greater effect on both the recurrence times (Figure 10b), and the variance in recurrence times (Figure 10c) of shorter faults. The dependence between slip rate variability and fault length has only been observed in extensional regions, and therefore may not apply to localized, plate-boundary faults. However, the systematic variations in earthquake recurrence intervals with fault length that arise not only in our model of the last 667 years in central Italy but also in observations from New Zealand suggest that in extensional regions, Coulomb stress interactions play an important role in influencing irregular earthquake recurrence intervals over thousands of years. Elastic stress interactions therefore provide a plausible physical explanation for the mechanism controlling observations of greater fault slip rate variability on shorter faults that is not readily explained by other possible causes of slip rate variability.
The effects of Coulomb stress interactions on smaller faults only become apparent due to the detail of the fault map used in our model which includes all (97) major faults in the region. In contrast, other studies often simplify the major faults in the region (e.g., Robinson [2004] and Robinson et al. [2009] in the Taupo Rift and Marzocchi et al. [2009] and Marzocchi and Melini [2014] in central Italy) or only resolve stress changes on a small number of faults (e.g., Li et al. [2015] in the Shanxi Rift, China, and Verdecchia and Carena [2016] in eastern California and western Nevada). By resolving stress changes on all faults in the region, we were able to explore the characteristics of the complicated interactions in areas of distributed faulting.
Our modeling suggests that the large-scale fault orientation appears to have little or no systematic control on the magnitude of fault interaction ( Figure S4) over timescales of 10 1 -10 3 years (i.e., beyond aftershock sequences). It appears that at the current stage of development of the fault system (i.e., beyond incipient localization), fault location has become more important than fault orientation on controlling fault interactions [e.g., Cowie, 1998].
As fault location appears to be more important than fault orientation on the impacts of stress interactions, areas with a greater density of shorter faults may be more likely to be susceptible to the impacts of Coulomb stress changes. Nicol et al. [2010] suggested that the activity of the shortest faults in the Taupo Rift is controlled by the largest earthquakes in the region (i.e., the earthquakes that occurred on the largest faults). In contrast, in the offshore part of the Taupo Rift, the spatial distribution of short faults is varied, with localized deformation and associated long faults in the southwest of the Whakatane Graben, and distributed deformation and a higher density of short faults in the northeast of the region [Nixon et al., 2014]. In this example, the shorter faults that are most susceptible to Coulomb stress changes are not situated in the areas where the largest earthquakes could occur as they are not adjacent to the longest faults. Therefore, more moderate-sized earthquakes may cause greater perturbations in earthquake recurrence intervals as these earthquakes occur in areas surrounded by a higher density of shorter faults. This may explain why earthquakes in the central Apennines since 1349 have clustered on the northeastern side of the mountain range ( Figure 1c): there appears to be a higher density of faulting in this region, especially surrounding the L'Aquila valley, yet the historical earthquakes are of a more moderate magnitude (M6-M6.7; Figure 2a). In contrast, the southwest of the region has hosted the largest-magnitude earthquakes (e.g., the M e 7.0 1915 earthquake) and may have a higher proportion of longer faults (Figure 2a).

Conclusions
The central Apennines, similar to many other areas of extensional tectonics, show evidence of spatially variable earthquake clusters and temporally variable fault slip rates over tens to thousands of years. To investigate the role of elastic fault interactions in controlling these observations we modeled the coseismic and interseismic Coulomb stress changes for 30 earthquakes in central Italy over 667 years on all active faults in the region.
We demonstrate that with a simple model that uses linear inputs for interseismic stress accumulation based on real measurements of the rate of movements of faults, coseismic Coulomb stress changes caused by the 30 earthquakes lead to highly nonlinear outputs. Our model illustrates how Coulomb stress changes can account for variability in the earthquake recurrence intervals of faults of hundreds to thousands of years.
Despite the interseismic loading stresses being on the order of <10 À2 MPa yr À1 and coseismic stress changes being on the order of 10 À1 MPa for each event, coseismic stress shadows are a pervasive feature of the sequence of earthquakes modeled with approximately 20% of faults in negative stress shadows at any one time. We also find that the pattern of faults that fail in each earthquake is not random, rather earthquakes are more likely to occur on faults with higher levels of Coulomb stress.
We calculate the earthquake delay or advance due to coseismic Coulomb stress changes and find that earthquakes can be advanced or delayed by hundreds to thousands of years. Furthermore, the susceptibility of faults to delays or advances in recurrence times appears to vary systematically with fault length (Figure 10). Faults with lower interseismic loading rates, which are generally the shorter faults, are more susceptible to coseismic Coulomb stress changes when compared to the longest faults in the array. This behavior has been observed in other extensional tectonic settings around the world suggesting that Coulomb stress interactions play an important role in controlling variable earthquake recurrence intervals in these settings and the magnitude of this effect may be a function of fault length.