Dual control of fault intersections on stop-start rupture in the 2016 Central Italy seismic sequence

Abstract Large continental earthquakes necessarily involve failure of multiple faults or segments. But these same critically-stressed systems sometimes fail in drawn-out sequences of smaller earthquakes over days or years instead. These two modes of failure have vastly different implications for seismic hazard and it is not known why fault systems sometimes fail in one mode or the other, or what controls the termination and reinitiation of slip in protracted seismic sequences. A paucity of modern observations of seismic sequences has hampered our understanding to-date, but a series of three M w > 6 earthquakes from August to November 2016 in Central Italy represents a uniquely well-observed example. Here we exploit a wealth of geodetic, seismological and field data to understand the spatio-temporal evolution of the sequence. Our results suggest that intersections between major and subsidiary faults controlled the extent and termination of rupture in each event in the sequence, and that fluid diffusion, channelled along these same fault intersections, may have also determined the timing of rupture reinitiation. This dual control of subsurface structure on the stop-start rupture in seismic sequences may be common; future efforts should focus on investigating its prevalence.


Introduction
In regions of distributed continental faulting, networks of active faults are commonly segmented on length scales of 10-25 km, approximately equal to the seismogenic thickness of the Earth's crust (Scholz, 1997;Stock and Smith, 2000;Klinger, 2010). This intrinsic maximum fault size limits the magnitude of continental earthquakes that rupture a single fault or segment to <M w ∼ 6-7 (Pacheco et al., 1992;Triep and Sykes, 1997), depending on local seismogenic thickness and fault geometry. Therefore, large continental earthquakes above this threshold (Scholz, 1997) such as the 2010 M7.2 El Mayor-Cucapah, Mexico, 2016 M7.8 Kaikoura, New Overview of epicentral region and fault geometry used in this study. (a) Regional tectonic map, showing mapped active normal faults (magenta, modified from Roberts and Michetti, 2004), up-dip surface projection of model faults displayed in Figs. 6-9 (coloured to match (c) and Figs. 3, 4 and 6), and bodywave focal mechanisms for each earthquake (see Fig. 2). White dashed line shows inferred east-dipping fault from Fig. 7, and relocated aftershocks from Chiaraluce et al. (2017) are shown in black. Locations of short-baseline GNSS instruments are shown by blue triangles. Black box shows extent of Fig. 3c. (b) Regional map showing the location of (a) and direction of regional crustal extension. (c) 3D cartoon of model fault geometry adopted in this study. Thick coloured lines show the surface projection of each fault and correspond to the coloured faults in (a) and Figs. 3, 4 and 6. (For interpretation of the colours in this figure and other figures, the reader is referred to the web version of this article.) ther dynamic and static stress transfer cause cascading failure of multiple critically-stressed faults or rupture is arrested before all these faults have failed. In the latter case the start of rupture in subsequent subevents determines the temporal evolution of the seismic sequence. Large earthquakes and seismic sequences have vastly different implications for seismic hazard: high hazard in a single event, or moderate hazard spanning years or potentially decades. But our understanding of what controls whether multifault rupture occurs over days to years or in seconds, and of what controls the spatio-temporal evolution of seismic sequences, has been severely limited by a paucity of high-resolution observations of modern seismic sequences.
Combined analysis of geodetic and seismological data can image stop-start rupture behaviour and address these questions, by disentangling the spatial pattern and temporal evolution of slip in seismic sequences at high resolution. A sequence of 3 M w > 6 earthquakes from August to November 2016 in the Central Apennine mountains, Italy (Fig. 1) presents a rare chance to investigate a seismic sequence with modern datasets and here we exploit seis-mological and field observations, as well as geodetic data, to image the kinematics of the sequence, and to understand structural and dynamic controls on its evolution. Our results suggest that structural complexity, namely the intersections between two sets of oblique faults, may have played an important dual role in the Central Italy seismic sequence: first by limiting the extent of individual ruptures and second by channelling fluid flow and controlling the timing of subsequent failure throughout the sequence.

Seismological constraint on earthquake source mechanisms
The Central Italy seismic sequence started with an M ∼ 6 earthquake on the 24th August 2016, and was followed by tens of thousands of aftershocks, including two large M > 6 events on the 26th and 30th October (Chiaraluce et al., 2017, Fig. 1). We refer to these three major earthquakes as the Amatrice, Visso and Norcia events respectively. The seismic sequence continued into 2017, with several earthquakes M < 5.7 on January 18th, but here we focus on the three largest events only. Left column shows the best-fit seismological focal mechanism (black) determined using MT5 software package (Zwick et al., 1994). For each of the Amatrice, Visso, and Norcia earthquakes, the composite geodetic mechanism from the geodetic solution is shown, as both the full moment tensor (dashed green line) and best double couple (solid green line). Best-fit seismological mechanism parameters and source time function are also given. Central panels show depth-misfit curves. Each point on the curve is determined by fixing the centroid depth to the given value, and inverting the waveform data to determine the best fit mechanism and source-time function at that depth. Vertical green bars show the centroid for the geodetically-derived slip distributions shown in Fig. 6 for the Amatrice, Visso and Norcia earthquakes, for comparison. Right panels show dip-misfit plots. On each, red indicates the SW-dipping plane, blue the NE-dipping plane. Waveforms and best-fit synthetics for each earthquake are shown in Supplementary Figs. S1-S3.
For each of these three earthquakes, we invert teleseismic longperiod body waves for the best-fit focal mechanism (Fig. 2, Supplementary Figs. S1-S3). We treat each earthquake as a finite-duration point-source centroid, with a moment-release function parameterised by a series of 1 s triangular elements. We invert P and S H waveforms to determine the moment-rate function, a centroid depth, total moment and a focal mechanism (strike, dip, rake), using a least-squares approach (e.g. see Walters et al., 2009).
We estimate a seismological M w of 6.2, 6.1 and 6.6 for the Amatrice, Visso and Norcia earthquakes respectively, and find similar normal-faulting mechanisms on NNW-SSE striking faults for each event. Our seismological results suggest shallow centroid depths of ∼4 km and relatively shallow dips for all three earthquakes (<∼45 • for the Visso and Norcia earthquakes, <∼50 • for the Amatrice earthquake; Fig. 2). The centroid depths estimated here agree well with previous seismological estimates (e.g., Chiaraluce et al., 2017;Pizzi et al., 2017). All these seismological results are consistent in indicating that the majority of moment release is concentrated within the upper ∼8 km of the crust, with centroids from 3-6 km depth for all three of the largest events. Comparison between the focal mechanisms estimated here, and the composite focal mechanisms resulting from finite-fault inversion of near source data shows a good agreement in strike and dip, with a slight (∼10 • ) difference in rake that likely results from collapsing a distributed pattern of slip across several faults into a simple composite mechanism (Pizzi et al., 2017).

Field measurement of surface ruptures
The normal faulting mechanism from our seismological results is consistent with the NE-SW extension that characterises the Apennine mountain belt (D'Agostino et al., 2011, Fig. 1b). The earthquakes occurred in the region of the SW-dipping Laga and Mt. Vettore (hereafter referred to as 'Vettore') normal faults. Of these two major faults, the Laga fault was thought to have last ruptured in a major earthquake in 1639 A.D. (Rovida et al., 2016), whilst the Vettore fault was only known to be active from palaeoseismic investigations (Galadini and Galli, 2003).
Following each of the three major events, we started mapping the surface ruptures the day after the earthquake, contributing to the compilation of measurements collated by the openEMERGEO international working group (Civico et al., 2018). Measurements of throw, heave, net slip and slip vector together with fault strike and dip were collected in the region of the Vettore and Laga faults, using a handheld compass clinometer and a ruler. For the Amatrice earthquake, measurements were collected over the 1 month following the earthquake. For the Visso earthquake, due to the short time between this event and the subsequent Norcia earthquake, the identification and mapping of the ruptures is likely incomplete and we were only able to collect measurements of the slip.
The Amatrice and Visso earthquakes each generated semicontinuous surface ruptures with ∼10-20 cm slip on pre-existing bedrock scarps, towards the southern and northern ends of the mapped Vettore fault respectively (Figs. 3,4). In contrast, the Norcia earthquake ruptured portions of the Vettore fault along its entire length, generating offsets up to 2.3 m along the central portions of the fault, and re-rupturing some of the same sections that failed in the earlier earthquakes (Figs. 3, 4). These sections re-ruptured with the same kinematics but approximately an order of magnitude greater slip. During the Norcia earthquake, numerous smaller faults in the hangingwall of the Vettore fault were also activated, including metre-scale displacement of an antithetic structure ∼2 km SW of the main fault (Fig. 4c). These field measurements are consistent with our seismological solutions and the regional extension direction.

Geodetic datasets
We measure the surface displacements during the seismic sequence with eighteen separate geodetic datasets: one regional GNSS dataset for each of the three earthquakes (INGV, 2016a(INGV, , 2016b, one short-baseline GNSS dataset for the Norcia earthquake (Wilkinson et al., 2017) and fourteen InSAR datasets from the Sentinel-1 and ALOS-2 satellites, each constraining the coseismic displacement fields from one or more of the three earthquakes ( Fig. 5g, Supplementary Fig. S4).
We processed Sentinel-1 Synthetic Aperture Radar interferograms using the GAMMA software package (http://www.gammars .ch) and unwrapped them using the MCF algorithm (Costantini, 1998). ALOS-2 interferograms were generated using the JPL/Caltech/Stanford ISCE package (https://winsar.unavco .org /isce .html) and unwrapped using the SNAPHU method (Chen and Zebker, 2002). Orbital effects were corrected using precise orbits from the European and Japanese Space Agencies respectively, and topographic effects were removed using 1-arcsec topographic data from the Shuttle Radar Topographic Mission (Farr et al., 2007). Unwrapping errors were manually checked and corrected. Interferograms were resampled in preparation for modelling using a nested uniform sampling approach (e.g. Floyd et al., 2016), with higher density in the nearfield and lower density in the farfield, to obtain about 1500 line-of-sight data-points per interferogram ( Supplementary Fig. S4).

Model fault geometry
In order to relate our geodetic measurements of surface displacement to slip on faults in the subsurface, we first define a simplified array of nine rectangular model faults (Fig. 1a, c, Supplementary Table 1). This level of complexity in source geometry is commonly required for geodetic modelling of multi-segment, moderate-magnitude events like the Norcia earthquake (e.g. the S. Napa, California and Darfield, NZ earthquakes; Floyd et al., 2016;Elliott et al., 2012), and requires that fault geometries are fixed prior to inversion, often on the basis of additional geological or geophysical constraints. For the Central Italy seismic sequence, we use a wealth of such additional information to define fault geometries at the surface and at depth, including: 1) discontinuities and low-coherence regions in our InSAR data (e.g. Fig. 5a, c, e), which are commonly indicative of surface or near-surface faulting; 2) our field mapping of surface ruptures ( Fig. 3; Civico et al., 2018); 3) relocated aftershock clouds (Chiaraluce et al., 2017; e.g. Supplementary Fig. S5) and; 4) our body-wave focal mechanisms.
Our model geometry primarily consists of four major fault segments: three segments for the main Vettore fault, and one for the northern Laga fault. The locations, strikes and dips of these four faults are well constrained by the above datasets.
We place three additional minor faults within the hangingwall of the Vettore fault; one antithetic and two synthetic (Fig. 1c). These minor faults are necessary to explain important near-fault complexity both in the geodetic displacements for the Norcia earthquake (e.g. the region between the central Vettore fault and the minor antithetic fault in Fig. 5c) and the complex array of decimetric surface ruptures we mapped in the field (Fig. 3). Tests showing the increased local misfit to these data when the minor faults are each removed from the model are shown in Supplementary Figs. S13-S24 and summarised in Supplementary Table 2. We note that whilst the geodetic and field data constrain the surface location of these three faults and slip in the shallow subsurface, the geodetic data are insensitive to their geometry at depths greater than 1-2 km due to strong trade-offs with slip on the main Vettore fault.
Two more faults represent: a 12-km long ENE dipping structure antithetic to the Vettore fault that we call the Norcia Antithetic fault; and a NE-SW striking 14-km long normal fault we call the Pian Piccolo fault that cross-cuts the Castelluccio plain between the Vettore fault and Norcia fault system to the SW (Figs. 1, 3c). The intersections between model faults at depth are supported by alignments of relocated aftershocks (from Chiaraluce et al., 2017) when projected onto our model fault planes (Figs. 6, 7). The Pian Piccolo fault is required to explain a strong NE-SW aligned signal in the InSAR data that covers the Norcia earthquake (e.g. Fig. 5e). In addition, the location and strike of this fault are supported by the geomorphology of the Castelluccio plain ( Fig. 3d), previous geological mapping (Coltorti and Farabollini, 1995) and by relocated aftershocks (Fig. 7), the latter of which also constrains the dip at ∼40 • . Tests removing this structure require major (>2 m) slip in the Norcia earthquake on the Vettore fault at depths greater than 10 km, making the geodetic centroid depth incompatible with the centroid depths obtained from body-wave seismology (Fig. 2, Supplementary Figs. S13, S15). The Norcia Antithetic fault is strongly delineated in the aftershock data (Chiaraluce et al., 2017) so we include this structure in our model geometry, with this fault truncated at its southern end by the Pian Piccolo fault. This structure is less well constrained by the geodetic data than the other eight faults, and we ran several tests with: the Pian Piccolo fault truncated instead by the Norcia Antithetic fault; the two faults crossing and neither truncating the other; and the Norcia Antithetic fault removed completely. Whilst our preferred geometry for the Norcia Antithetic fault (Fig. 1c) does improve the fit to both the InSAR and GNSS data for the Norcia earthquake, these alternative geometries only result in marginally higher misfit to the data. However, it is important to note that in all of these alternative geometries, the distribution of slip on the Vettore-Laga fault system, and therefore also the major findings of this study, remain the same.

Inversion for distribution of slip in the sequence
Each of the nine model faults are discretised into patches ∼1 km along strike × 1 km in depth (Supplementary Table 1). We solve for the distribution of slip across this fault array during four discrete intervals: three coseismic intervals associated with each of the three M > 6 earthquakes and one postseismic interval that follows the Amatrice earthquake and precedes the Visso earthquake (Fig. 5g, red stars and arrow). We jointly invert all geodetic data for slip in these intervals following the method employed by Floyd et al. (2016). Surface displacements are modelled as resulting from slip on rectangular dislocations in an elastic half-space (Okada, 1985), with shear modulus 3.23e10 and a Poisson's ratio of 0.25. We solve for two components of slip for each fault patch to allow spatially-varying rake, with a non-negative constraint on the inversion. We also force slip to be zero on certain masked regions of our model during the first three time intervals (hashed regions on Fig. 6). This is for two reasons: to prevent fitting of noise in geodetic data away from the earthquake of interest, and to prevent high slip on shallow model patches that have low temporal resolution (see below) and for which field mapping revealed no significant slip in the relevant time interval.
Relative to the InSAR data, the regional and short-baseline GNSS data are weighted by factors of 30 and 6 respectively, to take into account both the relative variance of the different datasets and the much larger number of InSAR measurements. GNSS uncertainties are those given as formal uncertainties, and InSAR covariance is estimated for each dataset by fitting an exponential function to the 1D radial autocovariance from a non-deforming region of the data (e.g. Funning et al., 2005). We tested variations in the relative weightings, and higher weightings of the InSAR data led to degradation in the fit to the GNSS data without significant improvement to the fit to the InSAR, essentially overfitting noise in the InSAR ( Supplementary Fig. S11).
We regularise our inversion using a Laplacian smoothing criterion, and choose a smoothing factor that represents a compromise between smoothness of the solution and goodness-of-fit to the geodetic data ( Supplementary Fig. S12). Whilst changing the smoothing factor within reasonable bounds changes the peak magnitude of slip, it does not affect the spatial pattern of slip.
We estimate uncertainties on our geodetic slip distributions using a Monte Carlo approach, (e.g. Funning et al., 2005, Supple-mentary Fig. S6), and estimate the spatial and temporal resolution of our slip model using a resolution spike test (Du et al., 1992;Supplementary Figs. S7 and S8).

Recovered distribution of slip
The recovered slip distributions for the three coseismic intervals are shown in Fig. 6, along with composite geodetic focal mechanisms. We compare these focal mechanisms to our body-wave solutions in Fig. 2, and find that the total moment release, the centroid depth, and the geometry of the SW-dipping nodal plane match extremely well in all cases. The slip vector (and hence auxiliary plane) in the seismological solutions shows some discrepancy with the composite geodetic moment tensors, with the geodetic solutions in each case incorporating a slight oblique component to the moment tensor, compared with the almost pure dip-slip seismological moment tensors. Fig. 6 shows that slip is confined to the top ∼6 km of the crust in all events, which supports the shallow, <4 km centroid depths from our body-wave solutions.
In addition, the magnitude and location of slip in the top km of the model agrees well with our independent estimates from the field (Figs. 6, 3), which we take as validation of our choice of model geometry, given the complexity of the mapped network of surface ruptures.
Our results show that slip from the three M > 6 earthquakes on the Vettore and Laga faults (Fig. 7) is spatially complementary and slip is restricted to <∼6 km depth, first-order observations consistent with the results of lower spatial-resolution seismological models (Chiaraluce et al., 2017;Pizzi et al., 2017). However, whilst previous geodetic models of the three events (Cheloni et al., 2017;Xu et al., 2017) also show this same general interdependence of slip, our results for the Visso and Norcia earthquakes show some differences to these studies. Namely, both these previous studies show significant (>60 cm) slip in the Norcia earthquake at depths 6-10 km on the southernmost segment of the Vettore fault, whereas both our geodetic and bodywave results suggest no significant slip took place below ∼6 km in the Norcia earthquake.
Our model also includes slip on several minor but important faults associated with the Vettore fault; the inclusion of slip on these faults affects the recovered distribution of slip on the main Vettore fault in the Norcia and Visso earthquakes (Supplementary Figs. S13-S24). Our Pian Piccolo and Norcia Antithetic faults are similar to two additional faults proposed by Cheloni et al. (2017). Our Pian Piccolo fault has a similar strike to the model fault that Cheloni et al. (2017) infer to be a reactivation of the Sibillini Thrust, but in our model this oblique structure has a steeper dip and projects to the surface ∼2 km further to the N (Fig. 3c). However, despite these differences in detail, both our studies recover the same oblique sense of slip on these two antithetic and oblique structures.
Differences between our results and those from previous studies arise due to several factors. As well as differences in geometry of the major faults, our inversion technique (Floyd et al., 2016), has enabled us to jointly invert all geodetic datasets throughout the sequence for mutliple slip events on a single fault geometry, remov-ing the need to discard geodetic datasets that contain combined coseismic signals from the Visso and Norcia earthquakes (e.g. Xu et al., 2017) or make assumptions regarding the spatial separation of these signals (e.g. Cheloni et al., 2017). This inversion approach has enabled us to take advantage of an extra 4-5 Sentinel-1 interferograms spanning either or both of the Visso and Norcia earthquakes than are not used in previous studies (Cheloni et al., 2017;Xu et al., 2017). In addition, these extra interferograms, along with additional unique short-baselines GNSS data and other geological and geophysical constraints on the sequence (field mapping, bodywave mechanisms) have meant we are able to reasonably add several smaller faults into our model, as these are required to explain the new near-field datasets.

Why does rupture stop?
The clear spatial interdependence of slip from the three M > 6 earthquakes on the Vettore and Laga faults (Fig. 7) has been used to suggest that structural complexity along the main Vettore fault strand may have influenced evolution of the sequence, with previous studies focusing on the termination of the Amatrice earthquake (Chiaraluce et al., 2017;Cheloni et al., 2017;Pizzi et al., 2017;Mildon et al., 2017). Here we use our new results to investigate this idea further and propose that the intersection of several oblique and potentially seismogenic faults with the Laga-Vettore system halted rupture for each earthquake in the sequence, therefore determining their respective magnitudes and preventing cascading rupture in a single earthquake of M w 6.7 or larger.
Active normal faults in the central Apennines are relatively young, having initiated within the past 2-3 million yr, and as a result the faults are segmented with lengths generally <20 km (e.g. Roberts and Michetti, 2004). Boundaries between neighbouring fault segments are known to commonly arrest through-going rupture in earthquakes (e.g. Biasi and Wesnousky, 2016).
The primary structural trend in the region of the 2016 earthquake sequence relates to the main NW-SE striking normal faults. However, there is also a secondary oblique structural trend represented by NNE-SSW to NE-SW striking faults, many of which are currently active as normal faults, that either formed in the current extensional tectonic phase or are reactivated thrusts as-sociated with the previous compressional regime (e.g. Pizzi and Galadini, 2009;Coltorti and Farabollini, 1995;Civico et al., 2017). These have been suggested to act as structural barriers to rupture on the main NW-SE striking normal faults they crosscut or intersect (e.g. Pizzi and Galadini, 2009), by forcing segment boundaries on the major faults.
We note linear features along this NNE-NE regional trend in both the relocated aftershocks (Chiaraluce et al., 2017, Fig. 3e) and geomorphology of the Vettore basin (Fig. 3c, d). The Piano Grande basin is rhomb shaped and bounded to the north and south by topographic ridges striking approximately 040 • (NE-SW). There are also sharp changes in slope observed at the northern and southern sides of the Piano Grande and Pian Piccolo basins, suggesting that oblique faults may have some control on the morphology of this elevated plain (Fig. 3d). This is supported by previous geological mapping which has identified normal faults bounding these basins to the north and south (Coltorti and Farabollini, 1995).
The relocated aftershocks show predominant alignment along a similar, but slightly different NNE-SSW 015 • trend (Fig. 3e). Projecting aftershocks of the Amatrice earthquake onto the model Vettore fault along NNE or NE trends, we find lineations of aftershocks separate the hypocentres and major rupture extents of all three earthquakes (Fig. 7a, c). These aftershock alignments predate the Visso and Norcia earthquakes, so we interpret them as intersections between the oblique faults and the main Vettore-Laga fault system.
The Amatrice earthquake initiated on the Laga fault and rupture propagated northwards onto the southern Vettore fault (Tinti et al., 2016), but the northern termination of significant slip in Fig. 8. Coulomb stress change on the Vettore-Laga fault system calculated prior to the Visso earthquake (a), prior to the Norcia earthquake (b) and following the Norcia earthquake (c). Warm colours indicate fault regions that have been brought closer to failure, cool colours indicate regions taken further from failure. White stars represent the hypocentres for the three main earthquakes as in Fig. 6. the Amatrice earthquake is bounded closely by the intersection of the NW dipping Pian Piccolo normal-oblique fault with the Vettore fault (Fig. 6a). Slip is needed on a structure with this geometry to explain the surface displacement during the Norcia earthquake (e.g. Fig. 5e) and we see evidence for it from relocated aftershocks (Chiaraluce et al., 2017, Figs. 6a, 7c) and in the large-scale geomorphology of the Piano Grande basin, as discussed above (Fig. 3d). The presence of a specific fault in this location from the surface geology is debated (Coltorti and Farabollini, 1995;Pierantoni et al., 2013), but we also note that some oblique faults in this region are poorly expressed at the surface and may only be revealed at depth by geophysical surveys (Civico et al., 2017). Some authors have suggested that the Sibillini Thrust, with a similar geometry to our Pian Piccolo fault but a much shallower dip (Fig. 3c, Pizzi and Galadini, 2009), is instead responsible for halting the northwards rupture of the Amatrice earthquake (Chiaraluce et al., 2017;Cheloni et al., 2017;Pizzi et al., 2017). We favour our interpretation of a steep structural barrier at depth over a shallowly dipping planar barrier (e.g. Cheloni et al., 2017, Fig. 7b dashed grey line) on the basis of the aftershock data, the geomorphology and our new Sentinel-1 interferograms. However, our Pian Piccolo fault may well join with a steepened (∼40 • dipping) lateral ramp of the Sibillini Thrust at depth (e.g. Pizzi et al., 2017, Fig. 7b solid grey line) and these two scenarios are likely indistinguishable in our data. In either case, the key interpretation remains the same: a similar steep structure is inferred to have stopped northwards rupture of the Amatrice earthquake.
Similarly, slip in the Visso earthquake is closely bounded up-dip and at its abrupt southern termination by two lineations of aftershocks, which we propose represent an unnamed fault with strike NNE and dip ∼55 • to the east (dashed white line on Fig. 1, dashed black line on Figs. 3, 5, 7, 9b) and a possible small conjugate fault with apparent dip to the NW (Fig. 7, dashed purple line).
These same suggested oblique faults also appear to constrain the overall pattern of major slip in the Norcia earthquake (Fig. 7a,  b). The structural control here is twofold: as before the faults may directly act as barriers to rupture, but in addition the limitations they place on slip in the two previous events leaves stress shadows (Fig. 8b), which also act to constrain the slip in this last event. In particular we suggest this indirect structural control plays an important role for the southern-termination of major slip in the Norcia earthquake; the Pian Piccolo fault appears to have not stopped rupture during the Norcia earthquake and the slip maxima in this event instead terminates against the slip maxima of the Amatrice earthquake, which acts as a stress-shadow (Figs. 8b, 7).
On a larger scale, it appears that the intersection of the Norcia Antithetic fault with the Vettore-Laga system may also have restricted slip in all three earthquakes to shallow depths of <∼6 km.
If this sequence ruptured only the shallower half of the Vettore-Laga system in ∼12-15 km thick seismogenic crust (Chiarabba and De Gori, 2016), this could suggest that the deeper portion (depths >6 km) of the fault system is able to fail independently, and crosscutting structures may result in depth segmentation (e.g. Elliott et al., 2011) as well as along-strike segmentation of seismic slip.

Why does rupture start again?
Using our model slip distributions, we calculate the changes in Coulomb stress (King et al., 1994;Lin and Stein, 2004) on our model faults throughout the 2016 seismic sequence (Fig. 8; Supplementary Fig. S9), caused by the slip in each of the 4 modelled time intervals (Fig. 5g). We also include Coulomb stress changes from recent events prior to 2016: the 1997 Colfiorito earthquakes and the 2009 L'Aquila earthquake. Slip distributions for the 1997 and 2009 events are constrained by an inversion of ERS data for the Colfiorito earthquakes ( Supplementary Fig. S25) and a previ- Fig. 9. Spatio-temporal evolution of aftershocks on the Vettore fault before the Visso earthquake. (a) Time evolution of aftershocks. Aftershocks following the Amatrice earthquake (Chiaraluce et al., 2017) are plotted in grey showing distance in the Vettore fault plane in the direction of the red arrow in (b). The coloured boxes contain the most distant 20% of aftershocks for progressive time-periods, and the triangles show the median distance within each box. The star shows the location and time of the Visso hypocentre, and the black dotted, solid and dashed curves correspond to diffusive models with f = 0.18 and K = 4.8, 3.8 and 2.8 m 2 /s respectively. (b) Postseismic slip estimated from our geodetic model on the Vettore-Laga fault system, shown as in Fig. 6. The red arrow shows the direction in which distance is calculated in (a), and the coloured aftershocks correspond to those shown in (a). Intersections of model and inferred faults are shown as black solid and dashed lines as in Fig. 7. ously published slip distribution for the 2009 L'Aquila earthquake (Walters et al., 2009). Coulomb stress change ( CFF) is defined as: where τ is the change in shear stress, σ n is change in normal stress and μ is the effective coefficient of friction. Stress changes were resolved in the direction of slip of each fault patch. Where a fault patch did not slip in one of the time intervals, stress changes were resolved onto a rake of −90 • . Coulomb stress change calculations were performed using the Coulomb 3.2 software (Lin and Stein, 2004) and a value of 0.4 was used for μ , with elastic parameters kept the same as for our geodetic and seismological inversions.
The Norcia hypocentre was brought closer to failure by 1.7 ± 0.18 MPa by all previous events, and given the short interval between the Visso and Norcia earthquakes, it is likely that static stress interactions brought the Norcia hypocentre to the brink of failure, precipitating its rupture 4 days later.
The Visso hypocentre was brought closer to failure by 1.2 ± 0.55 MPa by the Amatrice earthquake and subsequent afterslip (Fig. 8a). Static stress transfer alone might be able to explain the two-month delay between the Amatrice and Visso earthquakes, with delayed failure triggered by a rate-and-state frictional response (e.g. Kroll et al., 2017). However, we also find a northwards progression of aftershocks along the Vettore fault during these two months, which reaches the Visso hypocentre at the time of the earthquake (Fig. 9a, b). Northwards aftershock migration and triggering was inferred to be associated with fluid diffusion following the two largest previous earthquakes in this region (Miller et al., 2004;Malagnini et al., 2012), so we therefore also investigate this possibility in the following section.

Temporal migration of aftershocks
In order to investigate the spatio-temporal pattern of aftershocks in the interval between the Amatrice and Visso earthquakes, we take the earthquakes in this interval, projected onto the Vettore-Laga model fault (Fig. 7a) and first apply a time-varying filter to remove earthquakes with magnitude below the magnitude of completeness (estimated using the goodness-of-fit method, Fig. S5 in Chiaraluce et al., 2017).
We calculate the distance from the location of peak slip in the Amatrice earthquake to all aftershocks to the north, in the plane of the model Vettore fault and in a direction approximately perpendicular to the intersection of the Pian Piccolo fault with this plane. This is also approximately parallel to (and directed updip along) the eastward dipping lineation of aftershocks seen in Fig. 7a.
We plot this distance versus the timing of aftershocks following the Amatrice earthquake in Fig. 9a. We split these data into four successive time-intervals, each containing 150 earthquakes, and find the 20% 'most-distant' aftershocks for each interval. We consider these earthquakes to represent the leading-edge of any aftershock propagation, and see a clear temporal trend. We see no such pattern if we repeat the analysis to the south. This up-dip, northwards-only trend resembles that seen in studies following the nearby 2009 L'Aquila earthquake (Malagnini et al., 2012). Plotted spatially on the fault plane (Fig. 9b), the aftershocks appear to be propagating along the minor antithetic fault that ruptured in the Norcia earthquake, and the eastward dipping structure that we infer acted as a barrier to rupture in the Visso earthquake. This aftershock migration reaches the Visso hypocentre at the approximate time of the earthquake (Fig. 9a, b), suggesting a possible underlying triggering mechanism for the Visso earthquake.
Northwards aftershock migration following the two largest previous earthquakes in this region was suggested to be driven by diffusion of over-pressured fluids from the region of mainshock rupture (Miller et al., 2004;Malagnini et al., 2012). We find the temporal evolution of aftershocks is also consistent with a similar process. If we plot the median distance of the 20% subset of aftershocks for each time interval, these points are consistent with a diffusive-like temporal trend (Fig. 9a). We consider a simple 1D model of a steady-state source of overpressured pore fluid that diffuses along the fault plane following the Amatrice earthquake (equation (8) in Malagnini et al., 2012). If aftershocks were triggered when the pressure increased by a fraction f of the difference between the overpressured source and the background hydrostatic pressure, then the aftershock sequence should propagate according to: where x is distance, t is time, K is diffusivity and erfc −1 is the inverse complementary error function. Varying f and K , we find that forward models with K varying between 2.8 and 4.8 m 2 /s and f = 0.18, corresponding to an 18% increase in pressure above the background, reasonably fit the data (Fig. 9a). We therefore find that the temporal evolution of aftershocks is consistent with a diffusive process and appears to be spatially focussed along the fault intersections described in the previous section. This pattern could also arise from mechanisms other than fluid migration, such as the propagation of afterslip. However, we note that the magnitude of postseismic slip in this interval is predominantly zero within uncertainty bounds (Fig. 9b, Supplemen-tary Fig. S6), and, in addition, aseismic-slip driven migration is typically over two orders of magnitude faster than the rates found here (Roland and McGuire, 2009). We highlight that more detailed analysis of additional geodetic data should be undertaken to fully rule out this alternative hypothesis: slip in the postseismic interval has high uncertainty here as it is constrained by three interferograms only, all of which also include coseismic signals (Fig. 5g).
If we do consider this process as driven by diffusion of porepressure CO 2 or water along the fault intersections, then we can estimate a permeability for these regions using our diffusivity estimates and equation (3) in Townend and Zoback (2000), keeping the same values for lithological and fluid parameters suggested in Malagnini et al. (2012). We obtain a value of 2.4 to 4.1 ×10 −14 m 2 , which is consistent with fractured bedrock limestone and previous estimates of fluid permeabilities along nearby faults (Townend and Zoback, 2000;Miller et al., 2004;Malagnini et al., 2012).
Irrespective of its exact nature, we suggest that the process driving the northwards propagation in aftershock activity through time brought the Vettore fault significantly closer to failure as it traversed the ∼12 km towards the Visso hypocentre over an interval of approximately two months. However, this diffusive process evidently did not trigger failure of the intervening Norcia segment. We suggest that if the process was fluid-driven and constrained to the fault intersections, it could have bypassed the Norcia hypocentre (Fig. 9b). This may explain why the seismic activity in the sequence jumped from the southern to northern ends of the Vettore fault before finally rupturing the central section.

Implications for multi-fault failure, seismic sequences and seismic hazard
Structural complexity in fault networks (including gaps, bends, stepovers and intersections between faults) sometimes appears to halt rupture propagation during earthquakes and sometimes permits through-going rupture, allowing large multi-segment earthquakes (e.g. Biasi and Wesnousky, 2016). However, whilst numerical models of dynamic rupture support this role of structural barriers (e.g. Oglesby, 2008), palaeoseismological records cannot determine the relative importance of these features in halting real earthquake ruptures, with respect to the effects of the unknown distribution of pre-earthquake stress across fault networks.
Our results from the Central Italy seismic sequence provide important real-world constraint on this problem, simply because we can consider the sequence as a failed multi-segment earthquake. As the different fault segments in our case were necessarily all near-critically stressed at the beginning of the sequence, our study clearly demonstrates the importance of pre-existing structure in stopping small earthquakes from becoming larger ones. Since the static stresses involved in eventual triggering of the Norcia earthquake are significantly smaller than stresses at the crack tip during dynamic rupture, it is likely that the Laga-Vettore fault system would have ruptured in a single large earthquake if pre-existing structure had not arrested slip. However, it is also important to note that despite this clear structural control for most of the seismic sequence, the Norcia earthquake appears to have ruptured through a barrier that halted the Amatrice earthquake. Our study therefore highlights that structural barriers appear to play a vital but enigmatic role in determining whether a large earthquake or a seismic sequence occurs on a segmented, critically stressed fault system.
Our results also suggest that these same structural barriers may have controlled the order and timing of earthquakes throughout the subsequent seismic sequence. We suggest that not only did the migration of pressure-driven fluids determine the temporal delay between the Amatrice and Visso earthquakes, but that the channelling of fluids along fault intersections caused the 'out-of-sequence' failure of the northern Vettore fault before the central portion.
Fluid-driven migrating clusters of seismicity are thought to be a common tectonic process across a range of tectonic regimes and environments (Vidale and Shearer, 2006), with fluids preferentially migrating along relatively high-permeability faults, and triggering seismicity due to increased pore-pressure. This process is likely to be more common in normal-faulting regions such as Central Italy (Chen et al., 2012) than in other tectonic environments, and has been proposed to play a role in other recent seismic sequences in the region (e.g. Miller et al., 2004).
Channelling of fluids along fault intersections is also to be expected. Significant variability in permeability is likely even across a single fault plane, and fault step-overs and transfer zones, intersecting oblique and antithetic faults are all thought to act as higher-permeability conduits (Sibson, 1996), focussing flow in these regions that are typically highly fractured and critically stressed in the long-term. This is supported by an investigation of over 200 geothermal systems in the western US (Faulds et al., 2011); over 50% of such sites are located on structural complexities in a region of active normal faulting. Conversely, only 6% of sites are located in mid-segments or at the point of maximum longterm displacement on major faults, implying low-permeabilities in these regions, possibly due to thick layers of clay gouge. This may also explain why migrating fluids bypassed the region of the Norcia earthquake hypocentre, which was located in a fault mid-segment away from the inferred high-permeability fault intersections.
We suggest here that intersecting structures in fault networks may play an important dual role in controlling the mode and timing of multi-segment failure: first acting as barriers to rupture and determining whether multiple segments fail together as one big earthquake or separately as several; and second in controlling the timing and order of subsequent earthquakes if failure does occur in a protracted seismic sequence.
Finally, we should also consider the possibility that these same processes operate on much larger spatio-temporal scales. The 2016 Central Italy seismic sequence may be part of a multi-decadal sequence of clustered seismic activity along with nearby earthquakes in 1997 and 2009 (Salvi et al., 2000;Walters et al., 2009). It has been suggested that previous decadal 'super seismic sequences' occurred in the Apennines in the 15th and 18th centuries (Chiarabba et al., 2011;Wedmore et al., 2017). We therefore may need to reconsider traditional concepts of the earthquake cycle, to include protracted coseismic periods that span decades rather than seconds.

Conclusions
The 2016 Central Italy earthquake sequence highlights the influence that structural complexity in fault systems may have in controlling and segmenting the rupture of critically-stressed faults, particularly in continental fault networks. Intersecting faults can act to limit rupture, but it is unclear under what conditions this will occur, and our results reinforce the recent suggestion that the final magnitude of a complex-rupture earthquake cannot be determined until rupture has stopped (Wei et al., 2011). In addition, these same structural barriers may play a dual role, also controlling the timing of failure in seismic sequences by channelling pressuredriven fluid flow along fault planes.
The 2010 El Mayor-Cucapah earthquake in Baja California represents an important counterpoint to the Central Italy sequence. Despite the difference in tectonic region and style of faulting, both episodes of strain release comprised multiple sub-events, taking place on a complex network of fault segments (Wei et al., 2011) and both featured fluid-driven migration of earthquakes following the initial onset of seismicity (Ross et al., 2017). However, at El-Mayor Cucapah, the fault network failed as a single earthquake, whereas in Central Italy, failure occurred sequentially over several months. In both cases, complexity of fault structure is key to understanding the pattern and evolution of seismic strain release.
A better understanding of the factors that may halt rupture, and the range of dynamic processes that may lead to cascading rupture over timescales ranging from seconds to years is critical for improving our ability to predict whether faults will fail in large, complex, earthquakes, or in temporally-distributed seismic sequences of multiple large events -two end member scenarios with very different implications for seismic hazard.