Coulomb stress transfer and fault interaction over millennia on non-planar active normal faults: the Mw 6.5–5.0 seismic sequence of 2016–2017, central Italy

In order to investigate the importance of including strike-variable geometry and the knowledge of historical and palaeoseismic earthquakes when modelling static Coulomb stress transfer and rupture propagation, we have examined the August–October 2016 A.D. and January 2017 A.D. central Apennines seismic sequence (Mw 6.0, 5.9, 6.5 in 2016 A.D. (INGV) and Mw 5.1, 5.5, 5.4, 5.0 in 2017 A.D. (INGV)). We model both the coseismic loading (from historical and palaeoseismic earthquakes) and interseismic loading (derived from Holocene fault slip-rates) using strike-variable fault geometries constrained by fieldwork. The inclusion of the elapsed times from available historical and palaeoseismological earthquakes and on faults enables us to calculate the stress on the faults prior to the beginning of the seismic sequence. We take account the 1316–4155 yr elapsed time on the Mt. Vettore fault (that ruptured during the 2016 A.D. seismic sequence) implied by palaeoseismology, and the 377 and 313 yr elapsed times on the neighbouring Laga and Norcia faults respectively, indicated by the historical record. The stress changes through time are summed to show the state of stress on the Mt. Vettore, Laga and surrounding faults prior to and during the 2016–2017 A.D. sequence. We show that the build up of stress prior to 2016 A.D. on strike-variable fault geometries generated stress heterogeneities that correlate with the limits of the main-shock ruptures. Hence, we suggest that stress barriers appear to have control on the propagation and therefore the magnitudes of the main-shock ruptures.


I N T RO D U C T I O N
The 2016-2017 A.D. seismic sequence of central Italy began on the 24th August at 03:36 with an M w 6.0 earthquake (Event 1, magnitude from the INGV) occurring close to the town of Amatrice. This earthquake caused the destruction of many small towns in the region and 297 deaths. Numerous aftershocks were recorded in the following weeks, which followed an Omori-type decay (Peruzza et al. 2016). The decay was interrupted on the 26th October by two earthquakes that occurred in quick succession (19:10 and 21:18 (Event 2) local time, M w 5.4 and 5.9 respectively, magnitude from the INGV) to the north, centred on the town of Visso. No deaths were reported, but there was further damage to buildings in the region. The largest earthquake in the sequence occurred on the 30th October at 08:40, M w 6.5 (Event 3, magnitude from the INGV) between, and re-rupturing, the locations of the 24th August and 26th October events (Fig. 1). Three indirect deaths were reported, but buildings experienced further serious damage, including the collapse of the San Benedetto basilica in Norcia. The time since 18th January 2017 A.D. has been marked by the occurrence of further seismicity in the SE of the epicentral area including events of M w 5.1, 5.5, 5.4, 5.0 (magnitudes from the INGV) which may have indirectly triggered an avalanche that resulted in a further 29 deaths.
The two previous large (≥M w 6.0) earthquakes in this region prior to 2016 A.D. occurred to the northwest (1997 A.D. Umbria-Marche seismic sequence, Cello et al. 1998;Vittori et al. 2000) and the southeast (2009 A.D. L'Aquila earthquake, Alessio et al. 2010;Chiaraluce et al. 2011). The main shocks (with focal mechanisms) and aftershocks from these three events are shown in Fig. 1. These events both occurred along strike of the location of the 2016-2017 A.D. seismic sequence, leading to discussion of the possibility that the 2016 A.D. seismicity filled a seismic gap (Stein & Sevilgen 2016) with both the 1997 A.D. and 2009 A.D. events increasing the stress in the 2016 A.D. epicentral area (Nostro et al. 2005;Walters et al. 2009). If this region is a seismic gap, we are 1206 C The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.  For the examples shown in this paper, the specified length is 1 km. (2) Using an assigned projection direction that is perpendicular to the mean strike of the fault, the segment boundaries (shown as crosses) are projected to depth. interested in why the pattern of main shocks began in the south, jumped to the north, and then finally ruptured in between the two previous 2016 A.D. main shocks. In this paper we only consider the static stress transfer during the sequence, as the time between the earthquakes is probably too long to be due to dynamic stress transfer.
In particular, we have calculated the static Coulomb stress changes (Harris & Simpson 1992;Reasenberg & Simpson 1992) leading up to and during the 2016 A.D. seismic sequence on faults with strike-variable geometry. We add additional information to the discussion of the possible seismic gap: the elapsed times since previous large magnitude earthquakes on the Laga, Mt. Vettore and Norcia faults and modelling of all coseismic events and interseismic loading since 1349 A.D. We emphasise that the Mt. Vettore fault most likely had the longest elapsed time in the region, and hence, presumably, the highest accumulated stress at depth. We also emphasise that the accumulated stress on the Mt. Vettore fault prior to 2016 A.D. was most-likely non-uniform, due to the spatial distribution of historical events and its strike-variable geometry. We discuss whether such stress heterogeneity on the Mt. Vettore fault represented a stress barrier, such as those that have previously been hypothesised to limit the area of fault that ruptures during earthquakes (Steacy & McCloskey 1998;Rydelek & Sacks 1999). We suggest that stress heterogeneities may be used to explain the sequence and limits of the three main shocks that occurred during the 2016 A.D. seismic sequence.

M E T H O D
The surface traces of active normal fault scarps in Fig. 1 and subsequent figures are part of an active fault database compiled in Google Earth TM based on extensive fieldwork (>14 000 slip vector and fault orientation measurements at >800 sites, from Morewood & Roberts 2000;Roberts & Michetti 2004;Papanikolaou et al. 2005;Papanikolaou & Roberts 2007;Roberts 2008;Faure Walker et al. 2010Wilkinson et al. 2015;Mildon et al. 2016a), satellite imagery (Google Earth TM ), palaeoseismology studies (e.g. Galli et al. 2008) and geological maps. The locations of these faults are constrained to within a few tens of metres and the strike and dip of these fault are known to ∼±5 o accuracy. These uncertainties do not greatly affect the generation of strike-variable fault planes and subsequent modelling of Coulomb stress transfer. It has been shown that the faults that have been mapped are sufficient to accommodate most or all of the tectonic strain across the width of the Apennines because strains implied by fault scarps that have formed since 15 ± 3 ka match those from GPS when summed across the width of the Apennines (Faure Walker et al. 2010) and correlate with large-scale features such as orogen-wide topography, total throws (vertical component of the slip) across the faults and SKS splitting delay times in the underlying mantle (Faure Walker et al. 2012;Cowie et al. 2013). Other minor structures may also exist with little or no geomorphic expression; we do not include such features in our modelling. The traces of the active fault scarps we include are used to build 3-D fault planes comprised of numerous rectangular elements with strike variable geometry and a constant value of dip. This is a first-order approximation, because the fault geometry could change at depth, but few data constrain this. The dip that we use for each fault is determined from field measurements or from the literature (e.g. Boncio et al. 2004;Pace et al. 2006). The size of elements used is 1 km; this is a trade-off between observations and computational time. The strike-variable geometry is defined by field measurements (Roberts & Michetti 2004;Roberts 2007;Faure Walker et al. 2009, 2010Wilkinson et al. 2015;Mildon et al. 2016a). Faults are assumed to project to depth perpendicular to the overall strike of the fault, maintaining a cylindrical geometry, and short faults are assumed to maintain an aspect ratio of 1. Again, this is a first-order approximation, because the fault geometry could change at depth, but few data constrain this. The thickness of the seismogenic zone and brittle faulting is assumed to be 15 km (Cowie et al. 2013). This method is outlined in more detail in Mildon et al. (2016c). Ductile shear zones are modelled to be projections of the brittle faults from 15-24 km depth. Stress changes between the strike-variable faults generated using this method are calculated in Coulomb 3.4 (Lin & Stein 2004;Toda et al. 2005), which uses an elastic half-space model.

Interseismic modelling
To model the interseismic loading of the active normal faults in the upper crust, ductile shear zones are modelled directly downdip of the brittle portion of the faults (Wedmore et al. 2015Wedmore 2017). The shear zones are modelled as slipping incrementally and transferring Coulomb stress to the base of the brittle upper portions of the faults. The rate of slip on the shear zones equals the rate of post 15 ± 3 ka slip measured in the field across the active bedrock scarps exposed at the surface (Roberts & Michetti 2004;Papanikolaou et al. 2005;Papanikolaou & Roberts 2007;Faure Walker et al. 2009, 2010Wilkinson et al. 2015;Mildon et al. 2016a) and new throw measurements presented herein (Fig. 2), following the mechanical model advocated by Cowie et al. (2013). The post 15 ± 3 ka sliprates measured at the surface are assumed to be representative of the long-term rate as the scarps are formed over numerous seismic cycles. When summed across strike, the surface throws and the strain rates they imply correlate with total throws, SKS splitting in the mantle and long wavelength topography (Faure Walker et al. 2012;Cowie et al. 2013). A cross-section illustrating this model setup and the resulting stress pattern is shown in Fig. 2(b). The interseismic slip distribution is constructed by linearly interpolating between the values of slip measured at the surface and zero at the tips of the faults. The rake is assumed to be −90 • for the entire extent of the shear zone. It is assumed that the shear zones have the same dip, and are direct continuations of the upper brittle portions of the faults. Again, this is a first-order approximation, because the shear zone geometry could change at depth, but few data constrain this.
To illustrate how post 15 ± 3 ka slip rates are determined we present examples of profiles across fault scarps in Fig. 2. The profiles are constructed using chain-surveying techniques using steel rules and clinometers, and a number of sites are checked with barometric altimeters and tripod-based LiDAR (e.g. Wilkinson et al. 2015). In particular, we note during field work the position of the upper slope, the degraded free face, the preserved free face, the colluvial wedge and the lower slope (see Papanikolaou et al. (2005) for further information). This information is used to construct the profiles and to interpret the Holocene throw. Scarps show natural variability in offset along strike and our experience of multiple scarp profiles (e.g. Wilkinson et al. 2015) leads us to assign an error of ±20 per cent to the throw values. The throw since 15 ± 3 ka on the Mt. Vettore fault used herein is 5.20 ± 1.04 m close to the centre of the fault (Fig. 2c). The throw on the Laga fault used herein is 9.82 ± 1.96 m close to the centre of the fault (Fig. 2d). The expression of the Norcia fault is less clear in the field as there is little well developed and exposed Holocene bedrock scarp. Hence using the relationship between fault length and post 15 ± 3 ka throw from forty-eight other faults in the Apennines (Roberts & Michetti 2004;Papanikolaou et al. 2005;Papanikolaou & Roberts 2007;Faure Walker et al. 2009, 2010Wilkinson et al. 2015;Mildon et al. 2016a), a throw of 11.9 m is assigned to the Norcia fault. This value will have a higher error due to the spread in data; hence an error of 50 per cent has been assigned. This gives slip rates of 0.35 ± 0.1, 0.65 ± 0.13 and 0.79 ± 0.43 mm yr −1 for the Mt. Vettore, Laga and Norcia faults respectively. These values are within error of values reported from palaeoseismology (Galadini & Galli 2000).
With regard to the time over which interseismic loading should be applied to each fault, we utilise information from the historical record and palaeoseismology. The historical record is thought to be complete for events >M w 5.8 since 1349 A.D. (Michetti et al. 1996) due to shaking records made by scholars through history and compiled in the CFTI catalogue (Guidoboni et al. 2007). The last earthquake >M w 6.0 on the Norcia fault was in 1703 A.D. (Blumetti 1995;Galli et al. 2005), giving an elapsed time of 313 yr. On the Laga fault, the last earthquake was in 1639 A.D. (Pace et al. 2002;Guidoboni et al. 2007), giving an elapsed time of 377 yr. On the Mt. Vettore fault, there is no known earthquake in the historical catalogue that can be associated with this fault, even prior to 1349 A.D. The elapsed time is determined from palaeoseismic trenching to be in the range of 1316-4155 yr B.P. (Galadini & Galli 2003).

Coseismic modelling
We have modelled thirty-four earthquakes with >M w 5.5 from 1349 A.D. to 2016 A.D. and calculated the coseismic stress changes These profiles are constructed using a metre ruler and a clinometer and noting the location of the upper slope, degraded fault plane, preserved free face and lower slope. UTM coordinates and heights of the start and end points are given. (e) Graph to show the relationship between the lengths of active faults in the central Apennines and the maximum throw measured along each fault. This relationship is used to assign a throw to the Norcia fault as this fault is poorly exposed in the field, so topographic profiles of the Holocene scarp cannot be taken.
resolved onto all the surrounding active faults for the central Apennines. Detailed slip distribution information is not available for these pre-instrumental earthquakes hence we have used simple linear slip gradients to produce concentric slip distributions. Assuming a triangular slip profile at the surface (following from field data presented in Roberts (2007)) and a triangular slip profile with depths produces a concentric slip distribution. At the time of writing no non-preliminary slip-distributions have been published for the 2016-2017 A.D. sequence based on modelling InSAR so we have also used so again we use simple linear slip gradients (although see below for more detail). We also note that such modelling tends to use simplified planar fault geometries, so we would be unable to directly apply them onto the non-planar faults we have measured in the fault and modelled, another reason why we have used simple linear slip gradients. The rupture length used to define the dimensions of the slip distribution is calculated from the magnitude of the event in the CFTI catalogue, using the relationship presented in Wells & Coppersmith (1994). If the calculated rupture length is shorter than the assigned fault that ruptures, then the location of the slip along the responsible fault is inferred from the distribution of shaking (Guidoboni et al. 2007). Palaeoseismic records are also used to help constrain which fault and how much of that fault ruptures (e.g. Galadini & Galli 2003;Galli et al. 2008). The distribution and dates of historical ruptures since 1349 A.D. is shown in Fig. 1. In addition to the historical ruptures in the CFTI catalogue, the three largest earthquakes in 2016 A.D. and the four M w > 5.0 earthquakes that occurred on 18th January 2017 A.D. have been modelled to determine the evolution of Coulomb stress during the sequence. For Event 1 (24th August), a slip distribution based on one released by the INGV, that is determined through preliminary modelling of InSAR fringes and observations of surface ruptures with ∼20 cm maximum surface offset throw (Livio et al. 2016), has been used to model the coseismic stress changes (INGV Working Group 2016). For Event 2 (26th October), no slip distribution has been published at the time of writing, so a simple concentric slip distribution has been used. This may well be consistent with the relatively simple concentric structure of the InSAR fringes for this event (COMET 2016 we measured the distance covered by the InSAR fringes parallel to the mean fault strike; this distance was 13 km. A simple concentric slip distribution was used, with zero slip at the edges of the fault and at 11 km depth. The depth of the rupture is based on the width of the fringes, which were comparable to Event 1. A maximum slip of 0.58 m in the centre of the fault was chosen to agree with the published moment magnitude from the USGS. Surface faulting was observed during this event of a few centimetres to decimetres (A. Michetti, personal communication, 2016) hence we include this in the model such that surface slip is 10 per cent of the maximum slip at depth (based on relationships of slip at depth and at the surface from historical and recent earthquakes; Wells & Coppersmith (1994)). For Event 3 (30th October), no slip distribution has been published at the time of writing; hence a simplified concentric slip distribution is also used. The slip distribution has the location of maximum slip skewed to the southeast and located directly down dip of the observed location of maximum surface slip (Mildon et al. 2016b). To a first-order approximation, this is likely to be consistent with unpublished InSAR fringes (COMET 2016) which appear to be skewed and more numerous to the south-east. It is assumed that the whole of the Mt. Vettore fault ruptured during Event 3, based on the observations of surface ruptures and the extent of the InSAR fringes (Iezzi et al., in preparation). Again, the maximum slip of 2.05 m was chosen to agree with the published moment magnitude. We acknowledge that these are simplified slip distributions. However, published slip distributions for other recent earthquakes in Italy are non-unique (e.g. see Wilkinson et al. 2015) and are most commonly generated assuming a planar geometry. Hence we have taken the most prominent features of the available information (InSAR fringes and surface observations) to generate the simplified slip distributions we have used on strike-variable faults, and this could be improved in future iterations.
The four earthquakes that occurred on 18th January 2017 A.D. (M w 5.1, 5.5, 5.4, 5.0) can be related to the Laga fault, based on their location. No published slip distributions or INSAR results are available. The stress change associated with these events has been included, despite their small magnitude. Using the magnitude for each event, the rupture area was calculated based on Wells & Coppersmith (1994). A rectangular or square rupture area was used for each event and a concentric slip distribution assigned. The rupture dimensions and maximum slip used are 3 × 4 km and 0.35 m, 6 × 5 km and 0.6 m, 5 × 5 km and 0.41 m, and 3 × 3 km and 0.28 m, respectively. The location of the patches of the Laga fault that slips is assigned based on earthquake locations provided by the INGV.

Cumulative transferred stress
In order to determine the stress history on the faults of interest, the interseismic and coseismic Coulomb stress transfers onto the brittle portions of the active faults are summed together over time. This summation is done for each individual element of the gridded strike-variable faults. As a starting point, we assume that the stress on all the faults was zero in 1349 A.D., as this is the time from when the historical record for M w > 5.8 is considered complete. This approach has been adopted following previous studies (Deng & Sykes 1997;Nalbant et al. 1998). However we also calculate the stress for the Mt. Vettore, Norcia and Laga faults implied by the elapsed times from palaeoseismology and the historical record, which in the case of the Mt. Vettore fault could extend back to 4155 yr B.P. We assume that following an earthquake, the Coulomb stress on the portion of the fault that slips drops essentially to zero. While this is likely a simplification, the stress must drop to close to zero, especially in the largest earthquakes, otherwise there would be a large accumulation of stress over many seismic cycles, which is unrealistic. Some off-fault deformation may also occur to relieve stresses, however in contrast to the stress drop experienced by the fault that ruptures, the magnitude of this stress relieved is likely to be much smaller.

Interseismic loading
The Coulomb loading rate from interseismic slip across the shear zones at depth has been calculated for all the faults, and shown in detail for the Mt. Vettore, Laga and Norcia faults. The accumulated interseismic Coulomb stress is heterogeneous with depth and distance along the fault as shown in Fig. 4 which shows the total interseismic stress accumulated since the last earthquake on each individual fault. The highest magnitudes of interseismic Coulomb stress are seen at the base of the brittle faults, adjacent to the shear zones (which are not shown in Fig. 4, but are shown in Fig. 2). The Coulomb stress decreases updip of the brittle fault plane. Fig. 4 shows that while the accumulated interseismic Coulomb stress in the upper portions of the faults are comparable between the Mt. Vettore, Laga and Norcia faults, close to the bottom of the brittle faults, there are marked differences in the accumulated stress between the different faults, due to the differences in elapsed time. It is assumed that any other forms of crustal and/or mantle deformation is constant across the three faults, due to their proximity. An important observation is that the Mt. Vettore fault has the highest magnitude of stress from interseismic loading at the base of the brittle fault, whether the minimum (1316 yr, 60 bars) or maximum elapsed time (4155 yr, 200 bars) is considered.

Coseismic changes during the 2016-2017 seismic sequence
The modelled slip distributions and implied coseismic stress changes from the three largest earthquakes in the 2016-2017 A.D. sequence have been modelled and summed with the cumulative interseismic and coseismic stress transfer since 1349 A.D. (Fig. 5). We show interseismic loading for the relatively long elapsed time on the Mt. Vettore fault as insets for each time step. Considering only the stresses that developed post 1349 A.D., prior to the onset of the sequence (Fig. 5a) Fig. 1). The northern termination of the rupture on the Mt. Vettore fault occurs at a prominent bend in the fault (B on Fig. 5) that exhibited a strongly negative prior Coulomb stress pattern. The line of maximum curvature on the fault at this bend plunges to the SW, parallel to the slip-vector we measured at the surface from this earthquake (∼220 • -240 • , Livio et al. 2016). So this bend does not represent a kinematic barrier to slip associated with the earthquake.
Our results show that Coulomb stress is reduced by ∼2 bars at this bend compared to adjacent fault sections; this may be an example of a stress barrier to rupture propagation (e.g. Steacy & McCloskey 1998). Following the first earthquake, the northern portion of the Mt. Vettore fault was the most positively stressed region along the fault (Fig. 5b). This is in part due to the positive stress transfer from the 1997 A.D. Umbria-Marche seismic sequence, but also in part due to the long interseismic loading since the last earthquake (see the inset showing interseismic stress since 4155 yr B.P.). It is this highly stressed patch on the northern portion of the Mt. Vettore fault that ruptured in the 26th October earthquake (Fig. 5c); this portion is also called the Mt Bove fault (Pizzi & Galadini 2009). Note that this rupture terminated in the SE where Coulomb stress becomes negative across another bend in the fault (B on Fig. 5).
We have no slip vector data (from surface observations) for this bend, but if the slip vector plunged SW parallel to the axis of maximum fault curvature at this bend then again, this would not have represented a kinematic barrier to slip during the earthquake. This is a sensible assumption given the variation of slip vectors on other normal faults (Roberts 2007;Faure Walker et al. 2010). However, again, the marked decrease in stress makes this location another candidate for a stress barrier to rupture propagation. The 26th October earthquake increased the stress on the shallow portions of central section of the fault by 0.3-1 bars. Note the two thin stripes of positive stress located on two fault bends extending from 15 km depth almost to the surface in our model in the portion of the fault un-ruptured during Event 1 or 2 in the sequence (Fig. 5c, B and B ). These highly stressed patches plus the highly stressed base of the northern and southern portions of the Mt. Vettore fault were included in the ruptured fault surface in the 30th October M w 6.5 earthquake with surface slip of up to 2.3 m (Mildon et al. 2016b), as shown in Fig. 5(d). The 30th October event increased the stress on the Laga Fault, with the 18th January 2017 A.D. (M w 5.1, 5.5, 5.4, 5.0) events occurring on a part of the fault that mostly shows positive Coulomb stress.

Stress history on the Mt. Vettore fault
A history of the mean stress across the Mt. Vettore fault is plotted in Fig. 6 Note that the 1328 A.D. earthquake was excluded from previous calculations because it is before 1349 A.D. and probably less well-constrained, however for the purpose of highlighting the stress evolution of the Mt. Vettore fault, it has been included here. The absolute stress between 4155 yr B.P. and 1328 A.D. is assumed to accumulate at the mean rate of interseismic loading (as shown in Fig. 3) as there is no reliable information regarding coseismic stress transfer during earthquakes in this time period. Hence the absolute value of stress is subject to a number of assumptions as outlined herein and would be subject to variation if any of these assumptions were changed. The magnitude of the stress drops and increases are calculated independently. It is interesting to note the positive and negative stress increases on the Mt. Vettore fault are in the range of 0-10 bars, with interseismic Coulomb stress changes are in the range of 0-60 bars or even 0-160 bars given the elapsed time on the Mt. Vettore fault. So it is unlikely that the interseismic loading will have been overwhelmed by stress decreases produced by earthquakes on faults across strike from the Mt. Vettore fault. This also implies that the major driving force causing earthquakes to occur is the long-term tectonic loading, and to a lesser degree the short-term coseismic Coulomb stress changes (although see Wedmore (2017) who shows that smaller faults are more sensitive to coseismic changes due to their relatively low slip-rates and hence interseismic loading rates).

January 2017 earthquakes
All four earthquakes with M w > 5.5 that occurred on 18th January 2017 A.D. have been modelled. The Coulomb stress transfer from these events and the cumulative stress pattern as of 19th January 2017 A.D. is shown in Fig. 7. These earthquakes occurred in a region of mostly relatively high Coulomb stress (2 to >10 bars, see Fig. 7a) on the Laga fault, with the exception of the third event, which occurred in a region of negative Coulomb stress. However, the second earthquake that occurred on this day was the largest event and may have increased the Coulomb stress in region that ruptured in the third event (along-strike). Also, given the short time duration between the second and third earthquakes (11 min), there may be some component of dynamic triggering involved. The coseismic stress transferred from these four earthquakes on the Laga fault has been added to the cumulative (interseismic and coseismic) stress since 1349 A.D. on this fault and is shown in Fig. 7. The patches of zero stress are the regions of the fault that slipped; these should be considered as minimum values of stress as we do not know if all the stress was relieved on these patches. Hence most of the Laga fault could be positively stressed at the present time.

D I S C U S S I O N
Coulomb stress transfer is routinely calculated following large earthquakes and used to infer the change in probability of the locations of aftershocks and subsequent main shocks (e.g. Harris & Simpson 1992;Reasenberg & Simpson 1992;Toda et al. 2012;Pace et al. 2014) and has been used for earthquake hazard assessment (e.g. Toda & Enescu 2011). Many such studies do not include strikevariable faults, long earthquake histories from the historical record or the elapsed time since the last earthquake on receiver faults, however it has been shown for the 2009 A.D. L'Aquila earthquake that the magnitude of Coulomb stress transfer varies around bends in normal faults (Mildon et al. 2016c) consistent with our results This shows that the coseismic stress loading from the two most recent earthquakes are minor, relative to prior historical events.
herein. Although valuable results can be obtained from planar fault geometries, such as the 1st order regional patterns of stress transfer, especially if they include interseismic loading and a large number of events from the historical catalogue(e.g. Wedmore et al. 2017), we argue that extra insights may be gained if we include strike-variable fault geometries and the known elapsed time on the faults of interest. This extra insight may help to explain the details of complicated seismic sequences like that in 2016-17 using Coulomb stress transfer (Fig. 5).
We suggest that the sequence could not be fully explained if: (1) The whole historical sequence was not modelled, because stress heterogeneities arising from coseismic stress changes from historical earthquakes would not be evident on the Mt. Vettore fault; (2) The differences in elapsed time on the surrounding faults were not taken into account, as this helps to explain why the sequence began on the Laga fault and the continuation of the sequence on the Mt. Vettore fault;  The locations are coloured coded with a grey scale to indicate the order in which they occurred. The earthquake that occurred in the negatively stressed region occurred 11 min after the largest earthquake to occur on this day. Therefore there may have been additional positive static Coulomb stress transferred to the location of the third earthquake, or dynamic triggering may have occurred. (b) The slip distribution of the four M w > 5.0 earthquakes on 18th January. Note that two of the slip patches overlap. (c) Cumulative interseismic and coseismic stress from 1349 to 19th January on the Laga fault. The traces of active faults are marked in green and the Norcia (NOR), Mt Vettore (VET) and Laga (LAG) faults are marked. The regions that slipped during the four earthquakes on the 18th January are shown by black dashed lines, in these regions the stress is displayed as zero, however this is a minimum as it is not known whether all the stress has been relieved.
(3) The variable geometry was not modelled because it is this that creates a region of highly negative Coulomb stress in the bends on the fault and we hypothesise that these acted as a stress barriers to the ruptures of the 24th August and 26th October main shocks (Figs 5b and c).
We also note that the relative geometry of the Norcia and Mt. Vettore faults affects the build up of Coulomb stress on the latter. During interseismic loading, the lower portion of the Mt. Vettore fault experiences its highest increase in Coulomb stress. When an earthquake occurs on the Norcia fault (e.g. the 1328 A.D. and 1703 A.D. earthquakes; Fig. 6), the coseismic stress transferred onto the Mt. Vettore fault varies with depth. The upper portion of the fault experiences negative coseismic Coulomb stress transfer, whereas the lower portion experiences positive coseismic stress transfer. This positively reinforces the region of positive stress generated on the lower portion of the Mt. Vettore fault by interseismic loading even though the faults are located across strike. This is due to the horizontal across-strike distance between the two faults. In addition, the distance between the two faults is relatively constant in the southern portion of both faults, however in the northern region, the geometry of the Norcia fault changes (the strike changes from NNW-SSE to NW-SE) and the distance between the faults increases. This affects the stress transferred to the northern portion of the Mt. Vettore fault, as shown on Fig. 6 for coseismic stress transfer from the 1703 A.D. earthquake (which likely ruptured the whole of the Norcia fault). This positively stressed region at the northern end of the Mt. Vettore fault is clear in the cumulative stress pattern in Fig. 5(a) and shows the importance of the 1703 A.D. earthquake in contributing to the heterogeneous stress pattern on the Mt. Vettore fault. We hypothesise that this earthquake in particular caused the northern end of the Mt. Vettore fault to be more positively stressed and contributed to the stress pattern associated with the earthquakes of 26th October to occur on the northern portion of the Mt. Vettore fault, skipping the central section. This indicates that the relative geometry between Downloaded from https://academic.oup.com/gji/article-abstract/210/2/1206/3858476 by Birkbeck College, University of London user on 13 November 2017 faults is important, as well as the variations in geometry along single faults.
The southern 8 km of the Mt. Vettore fault ruptured during both Events 1 and 3 as evident from surface ruptures (Livio et al. 2016) and INSAR results (COMET 2016) and our own field observations (Iezzi et al., in preparation;Mildon et al. 2016b). At first sight this suggests that our assumption that the stress on a fault drops essentially to zero following an earthquake is incorrect, because the same patch ruptured twice within ∼2 months. However, note that our models show that the base of the Mt. Vettore fault in its southernmost 8 km is positively stressed after the relatively shallow slip associated with the 24th August 2016 M w 6.0 event (up to 30 bars if calculated from 1349 A.D.; up to 200 bars if calculated from 4155 yr B.P.). Furthermore, after the 26th October, a region of highly positive Coulomb stress existed along the full 28 km length of the base of the fault with two apophyses of positive stress projecting to the surface at two prominent fault bends. This stress pattern appears to have contributed to the preparation of the Mt. Vettore fault for an M w 6.5 earthquake with up to 2.3 m surface slip and a ∼28 km rupture length (Mildon et al. 2016b). It may be that slip initiated at depth and was able to propagate up to near the surface via the positively stressed bends close to the central portion of the fault, before spreading along strike to the SE and NE portions of the fault. These sections of the fault had already slipped in the 24th August and 26th October earthquakes, but the stress may or may not have dropped to zero during these events. This suggests that zero or low Coulomb stress patches on faults that have ruptured in recent events can be ruptured through subsequent earthquakes, if the pattern of positive stress, especially at the base of the fault, is extensive.
Overall this study shows that when the historical record, elapsed time and fault geometry are taken into account, the notion that the earthquakes occurred due to a seismic gap between the 1997 A.D. and 2009 A.D. earthquakes is perhaps correct, but oversimplified. As shown in Fig. 6, the stress transferred to the Mt. Vettore fault from these two prior earthquakes was very small (∼0.03 bars for the 1997 A.D. earthquake and ∼0.01 bars for the 2009 A.D. earthquake) compared to the total accumulated stress on the fault (tens to hundreds of bars). The stress changes due to the other earthquakes shown in Fig. 6 (1328 A.D., 1639 A.D. and 1703 A.D.) are of much greater magnitude (−0.93, 0.36 and −2.22 bars respectively) and therefore have more influence over the stress on the Mt. Vettore fault, despite being much longer ago (assuming the stresses have not dissipated over time). However, it is possible that some static stress may have dissipated over this short time period and therefore the stresses presented in this paper should be considered to be a maximum. We also note that the greatest magnitude of received coseismic stress change on the Mt. Vettore fault is negative and occurs when earthquakes occur on the Norcia fault, that is, across strike. This shows that Coulomb stress drops play a more dominant role in delaying rupture than positive stress changes in preparing a receiver fault for rupture [see Wedmore et al. (2015Wedmore et al. ( , 2017 and Wedmore (2017) who first described this effect]. Moreover, the interseismic Coulomb stress loading must also be considered due to the high stress imparted to the base of the brittle faults. This emphasises the importance of including knowledge of the elapsed time on faults receiving stress (which can be thousands of years), even if knowledge of the coseismic stress changes are limited to a few hundred years by historical and palaeoseismic records, as the interseismic stress may dominate over relatively small magnitude fluctuations produced by coseismic stress changes. Finally, the above calculations suggest that a large portion of the Laga fault is positively stressed after the 18th January earthquakes.

C O N C L U S I O N S
Using an approach to model Coulomb stress changes on strikevariable faults, we suggest that the central Apennines 2016-2017 A.D. seismic sequence on the Mt. Vettore and Laga faults can be explained by Coulomb stress transfer, by including the coseismic changes from historical earthquakes and the interseismic loading for the elapsed time on the faults in the epicentral region. The first main shock occurred on the Laga and Mt. Vettore faults due to regions of high Coulomb stress on both faults. The rupture length on the Laga fault appears to have been controlled by the extent of the aftershocks from the 2009 A.D. L'Aquila sequence (Fig. 1). On the Mt. Vettore fault, the rupture appears to have been controlled by the presence of a fault bend creating a stress barrier. Following Event 1, the northern portion of the Mt. Vettore fault was the most positive due to the historical Coulomb stress transfer; hence the sequence jumped the central section. Four days later, in Event 3, the central section ruptured and the northern and southern sections re-ruptured due to high values of Coulomb stress along the entire base of the fault and two patches of positive stress extending almost to the surface on fault bends. The mean stress history on the Mt. Vettore fault as shown in Fig. 6 highlights the importance of considering the full historical catalogue of large earthquakes, as the coseismic changes from the most recent earthquakes in this example are negligible compared to previous stress drops and increases. Without including fault geometry, the historical catalogue and the elapsed time, we suggest that that the sequence of rupture locations with the jump to the north missing the central patch of the Mt. Vettore fault would be difficult to explain. It appears that it is only by including the fault geometry and coseismic stress changes since 1349 A.D. that generates stress heterogeneities, particularly on the Mt. Vettore fault, which can explain the distribution and limits of the main-shock ruptures. Finally, our calculations suggest that a large portion of the Laga fault is positively stressed after the 18th January earthquakes. Where future events occur will test the extent to which this approach can help inform seismic hazard analysis.

A C K N O W L E D G E M E N T S
This study was funded by NERC Studentship NE/L501700/1, NERC Urgency Grant NE/P01660X/1 (EEFIT Reconnaissance Mission to the Amatrice, Italy, 24/09/2016 Earthquake) and NERC Standard Grant NE/I024127/1. The method to model strike-variable was developed during JSPS Short Term Fellowship (PE15776) undertaken by Zoe Mildon. Work was aided by a GBSF (Great Britain Saskawa Foundation) grant 4602 to JFW. We thank Ruth Harris and an anonymous reviewer for their comments which helped to improve the manuscript.