Spatial migration of temporal earthquake clusters driven by the transfer of differential stress between neighbouring fault/shear-zone structures

Uncertainty concerning the processes responsible for slip-rate fluctuations associated with temporal clustering of surface faulting earthquakes is a fundamental, unresolved issue in tectonics, because strain-rates accommodated by fault/shear-zone structures are the key to understanding the viscosity structure of the crust and seismic hazard. We constrain the timing and amplitude of slip-rate fluctuations that occurred on three active normal faults in central Italy over a time period of 20 – 30 kyrs, using in situ 36 Cl cosmogenic dating of fault planes. We identify five periods of rapid slip on individual faults lasting a few millennia, separated time periods of up to 10 millennia with low or zero slip-rate. The rapid slip pulses migrated across the strike between the faults in two waves from SW to NE. We replicate this migration with a model where rapid slip induces changes in differential stress that drive changes in strain-rate on viscous shear zones that drive slip-rate variability on overlying brittle faults. Earthquakes increase the differential stress and strain-rate on underlying shear zones, which in turn accumulate strain, re-loading stress onto the overlying brittle fault. This positive feedback produces high strain-rate episodes containing several large magnitude surface faulting earthquakes (earthquake clusters), but also reduce the differential stress on the viscous portions of neighbouring fault/shear-zones slowing the occurrence of large-magnitude surface faulting earthquakes (earthquake anticlusters). Shear-zones on faults experiencing anticlusters continue to accumulate viscous strain at a lowered rate, and eventually this loads the overlying brittle fault to failure, initiating a period of rapid slip through the positive feedback process described above, and inducing lowered strain-rates onto neighbouring fault/shear-zones. We show that these patterns of differential stress change can replicate the measured earthquake clustering implied by the 36 Cl data. The stress changes are related to the fault geometry in terms of distance and azimuth from the slipping structure, implying that (a) strain-rate and viscosity fluctuations for studies of continental rheology, and (b) slip-rates for seismic hazard purposes are to an extent predictable given knowledge of the fault system geometry.

Uncertainty concerning the processes responsible for slip-rate fluctuations associated with temporal clustering of surface faulting earthquakes is a fundamental, unresolved issue in tectonics, because strain-rates accommodated by fault/shear-zone structures are the key to understanding the viscosity structure of the crust and seismic hazard.We constrain the timing and amplitude of slip-rate fluctuations that occurred on three active normal faults in central Italy over a time period of 20-30 kyrs, using in situ 36 Cl cosmogenic dating of fault planes.We identify five periods of rapid slip on individual faults lasting a few millennia, separated time periods of up to 10 millennia with low or zero slip-rate.The rapid slip pulses migrated across the strike between the faults in two waves from SW to NE.We replicate this migration with a model where rapid slip induces changes in differential stress that drive changes in strain-rate on viscous shear zones that drive slip-rate variability on overlying brittle faults.Earthquakes increase the differential stress and strain-rate on underlying shear zones, which in turn accumulate strain, re-loading stress onto the overlying brittle fault.This positive feedback produces high strainrate episodes containing several large magnitude surface faulting earthquakes (earthquake clusters), but also reduce the differential stress on the viscous portions of neighbouring fault/shear-zones slowing the occurrence of large-magnitude surface faulting earthquakes (earthquake anticlusters).Shear-zones on faults experiencing anticlusters continue to accumulate viscous strain at a lowered rate, and eventually this loads the overlying brittle fault to failure, initiating a period of rapid slip through the positive feedback process described above, and inducing lowered strain-rates onto neighbouring fault/shear-zones.We show that these patterns of differential stress change can replicate the measured earthquake clustering implied by the 36 Cl data.The stress changes are related to the fault geometry in terms of distance and azimuth from the slipping structure, implying that (a) strain-rate and viscosity fluctuations for studies of continental rheology, and (b) slip-rates for seismic hazard purposes are to an extent predictable given knowledge of the fault system geometry.

Introduction
Slip-rates on brittle faults in the upper crust fluctuate through time, but why this happens is a fundamental, unresolved issues in tectonics (Dolan et al., 2016).A growing body of results demonstrate the occurrence of temporal clusters of surface faulting earthquakes on individual faults, where several large magnitude earthquakes occur closely spaced in time, raising the slip-rate averaged across the cluster duration above that of the background, long-term slip-rate (e.g., Sieh et al., 1989;Marco et al., 1996;Rockwell et al., 2000;Klinger et al., 2000;Dawson et al., 2003;Dolan et al., 2007;Oskin et al., 2008;Ganev et al., 2010;Scholz 2010;Gold and Cowgill 2011;Grapes and Holdgate 2014;Clark et al., 2017;Klinger et al., 2015;Dolan et al., 2016;DuRoss et al., 2016;Lefevre et al., 2018;Verdecchia et al., 2018;Wechsler et al., 2018;Angster et al., 2019;Cooke et al., 2020;Hatem et al., 2021;Iezzi et al., 2021;Elston et al., 2022;Mildon et al., 2022;Dolan et al., 2024).These studies also demonstrate the converse, that is temporal anticlusters, where a lack of large magnitude surface-faulting earthquakes in given time periods produces a slip-rate that falls below that of the background, long-term slip-rate.Such fluctuations are, in general, poorly explained, with causal processes suggested to include stress interactions between neighbouring faults and between faults and underlying shear-zones, microstructural and rheological evolution, and fluid involvement (e.g see Dolan and Meade, 2017 for a review).They are poorly explained because we lack examples where measurements of the slip-rate fluctuations have been replicated by quantitative models of the causal processes.Thus, uncertainty concerning the causal processes remains, leaving a fundamental, unresolved issue in tectonics.This is important because (a) the strain-rates accommodated by combined fault/shear-zone structures that link between the brittle upper crust and viscous lower crust are the key to understanding the viscosity structure of the overall crust and hence crustal rheology and evolution (e.g.Moore and Parsons 2015), and (b) the strain-rates, and their variability in space and time, control the spatial and temporal occurrence of hazardous earthquakes (Cowie et al., 2012;Pace et al., 2014).
A theme emerging to explain slip-rate fluctuations across temporal clusters and anticlusters is that the viscous roots of major faults, that is mid-and lower crustal shear zones, may evolve through time, changing their resistance to deformation (Dolan et al., 2007(Dolan et al., , 2016; Dolan and Meade, 2017;Mildon et al., 2022).The idea is that changes in the viscous roots drive fluctuations in the slip-rate of overlying faults.The questions that arise are as follows.(a) If such fluctuations exist, what causes them, and can knowledge of the processes inform our understanding of tectonics?(b) Are the processes quantifiable to an extent that knowledge of the processes can be used in seismic hazard assessment?
For example, Dolan et al. (2007) describe an example with "out-of-phase" high and slow slip-rate pulses on faults in the Los Angeles region versus the faults in the eastern California shear zone.They note the likely presence of ductile shear zones underneath the faults, and speculate that periods of high slip-rate at the surface correlate with high strain-rate on underlying shear zones.They suggest that the high ductile strain-rates lead to strain-hardening that ultimately terminates a period of high slip-rate; slip is then transferred to a neighbouring shear zone that has been experiencing relatively low slip-rate and hence has had time to anneal and weaken.In other words, the deformation migrates temporally and spatially, seeking out the weakest portion of the system to activate during the ongoing regional deformation.This explanation is appealing, because Dolan et al. (2007) show that the relatively large distance between the faults in question (>~100 km) rules out static Coulomb stress transfer as an explanation.
However, it has been noted that for other regions of active faulting, "out-of-phase" fault activity occurs where the across-strike fault spacing is smaller, and stress transfer has a role to play.Mildon et al. (2022) provide an example of "out-of-phase" high and slow slip-rate pulses on two separate Holocene-active normal faults in central Italy.The faults are only ~30 km apart across strike, so Coulomb stress transfer is significant between these structures.Mildon et al. (2022) suggest that stress transfer has changed the differential stress on viscous shear zones beneath the brittle faults, promoting changes in strain-rate due to the power-law relationship for dislocation creep, where strain-rate, ε, and differential stress, σ, are related by ε ∝ σ n , with n = ~3 (Hirth et al., 2001).In this scenario, the deformation migrates temporally and spatially, away from fault/shear-zones whose differential stress has been reduced by stress interactions.
The above describes two appealing explanations for the clustering of surface faulting earthquakes involving fluctuations in the properties of underlying shear zones.It may be that these explanations are in some way linked, or instead act separately, or at different length-scales, but this is not presently known.Indeed, this is challenging to ascertain due to the difficulty in obtaining samples of active shear zones from the viscous crust.We need further case studies that allow the role of one or other of the two appealing explanations to be quantified.
This paper presents a case study where we have been able to quantify the timings and magnitudes of slip-rate variations that we ascribe to clustering for three faults in the central Apennines, Italy.We model the transfer of differential stress between the structures coupled to calculations of implied viscous strain-rate fluctuations.Our approach is to recover slip-histories for the brittle active normal faults at the surface using in situ 36 Cl cosmogenic dating.Using the recovered slip-histories, we calculate the changes in differential stress on neighbouring fault/ shear-zone structures through time.We then calculate the implied strain-rate changes on underlying shear zones through time by using a flow law for viscous deformation (Hirth et al., 2001).We investigate whether the strain-rate fluctuations are of sufficient magnitude to drive and hence explain the fluctuations in slip-rate on the overlying brittle faults.We have no data on rates of strain-hardening and annealing, or rates or fluxes of fluid involvement or microstructural changes, so it is challenging to include these in our calculations; our approach is to see if we can explain the slip-rate fluctuations without including these variables, thus testing whether stress interactions are the dominant control.
In particular, the study area in the high topography of the Apennine mountains in central Italy is prone to damaging large magnitude normal faulting earthquakes (Mw 5.8-6.9), with normal faults spaced ~10-20 km across strike and with fault lengths of ~20-40 km (Roberts and Michetti 2004, Figs. 1 and 2).Using our 36 Cl results, we document the systematic spatial migration across strike of five periods of rapid surface slip-rate between three sub-parallel active normal faults, the Pescasseroli, Scanno and Maiella faults that are respectively ~21 km, ~31 km and ~23 km in length.The migration occurs over the time period since ~20 ka, with each period of rapid slip lasting ~1.3-4.5 kyrs.We interpret the periods of rapid slip-rate, involving 3.5-8.0m of slip, as temporal clusters of surface faulting earthquakes because (a) surface slip is only produced by earthquakes with magnitudes >~Mw 6.0 in the central Apennines (Galli et al., 2008), and (b) magnitudes of slip of 3.5-8.0m are too large to explain with single earthquakes on faults whose length are only ~21-31 km in length (Wells and Coppersmith 1994).In other words, each period of rapid slip must have contained several > Ms 6.0 to < Ms 7.0 earthquakes.
Furthermore, evidence exists that the slip-rates averaged over multiple surface faulting earthquakes on the brittle faults in central Italy are controlled by underlying viscous processes.The region exhibits a spatial pattern of long-term slip-rates measured across brittle faults, averaged over 15 ± 3 kyrs, that correlates with topographic elevation (Faure Walker et al., 2012, Fig. 2).The correlation shows a power-law relationship where the stress exponent, n, is equal to ~3.26 (Cowie et al., 2013, Fig. 2b).This relationship between driving stresses and strain-rates resembles the power-law relationship for dislocation creep described above (Hirth et al., 2001).The relationship has been interpreted to suggest that a fault and its underlying shear zone behave as a paired system, with the strain-rate on the underlying shear-zone driving the slip-rate on the overlying brittle fault, and with the strain-rate on the shear zone related to the value for differential stress (Cowie et al., 2013).
The approach in this paper exploits this relationship between strain-rate and elevation (which contributes to the vertical driving stress) by investigating what happens if the differential stress on shear zones is perturbed due to interaction between one fault/shear-zone pair and neighbouring fault/shear-zone pairs.Our results show that implied changes in differential stress predict strain-rate and slip-rate changes that mimic and hence explain the firstorder pattern of slip-rate fluctuations measured on the faults with in situ 36 Cl cosmogenic dating of the fault scarps.This does not rule out the action of strain-hardening and annealing, microstructural evolution or  Chiaraluce et al., 2011;Pace et al., 2002).The mainshocks and aftershocks provide information on the geometries of the ruptured faults at depth and the approximate position of the brittle-viscous transition with most seismicity limit to above 10 km, but with tens of earthquakes at depths of up to 16-18 km (see Fig. 12 of Chiaraluce et al., 2011).(c) Schematic cross-section showing the interpreted geometries adopted in this paper to perform stress modelling calculations.Fault dips are based on surface measurements of outcropping fault planes not the aftershock data from 1984 to 2009.(d) Projections to 15 km depth of the surface fault traces, with fault traces located in the inset.Colours show annual interseismic loading rates on the brittle faults from Mildon et al. (2019) given continuous creep on viscous shear-zones at depth at the rates defined by surface offsets of landforms or deposits formed since the last glacial maximum at 15 ± 3 ka (see Fig. 2a for justification of this; the 15 ± 3 ka slip-rates are from measurements in Cowie et al., 2013).These time averaged loading rates will be modulated by stress changes produced by fault interactions (see Fig. 2c and d, and results in this paper).(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)fluid involvement (see Dolan and Meade, 2017), but does suggest that the primary character of the slip-rate fluctuations and clustering can be explained without them, so they are not the dominant control.Moreover, we argue that because the stress changes produced by interactions are related to the geometry of the fault system, the occurrence of slip-rate fluctuations and hence earthquake clusters is not a random Poissonian process.Instead, the occurrence of slip-rate variations has a systematic underlying character defined by dynamic interactions produced by slip in both the brittle and viscous domains.Our overall interpretation is that because stress interaction is influenced by fault/shear-zone geometry, which can be measured, interaction has aspects that are inherently quantifiable and hence we should be able to use this to improve probabilistic seismic hazard assessment.

Details of the model to be explored
As introduced above, our understanding of the relationship between brittle faults in the upper crust and their underlying shear zones has recently been clarified by a study of strain-rates and topographic elevation in the Apennines (Faure Walker et al., 2012;Cowie et al., 2013, Fig. 2).Although it has long been considered that brittle faults are underlain by viscous shear zones (e.g.Savage and Burford 1973;Sibson Fig. 2. The relationship between strain-rate and topographic elevation in the Apennines and implications for stress modelling of earthquake clustering.(a) Strainrates averaged since 15 ka from offsets across fault scarps exposed at the surface.Values are summed across all faults within 5 × 90 km boxes traversing the Apennines mountains.(b) Strain rates from (a) are plotted against average elevations from SRTM data in 5 × 90 km boxes.A power law relationship exists, with an exponent of 3.26.This resembles the power-law relationship for dislocation creep, where strain-rate is proportional to differential stress raised to the power n, with n = ~3 (Hirth et al., 2001).Topographic elevation contributes to the vertical stress, implying that it contributes to the differential stress at depth.Thus, it is implied that the slip-rates on the brittle faults are driven by the viscous strain-rates on underlying shear zones, because the strain-rates are controlled by differences in elevation and hence differential stress in the viscous layer.(c) Brittle-elastic strain in the upper crust (spring and frictional ratchet) is matched by visco-elastic strain at depth (spring and dash-pot) when summed over numerous earthquakes cycles.Strain in the upper crust, and how it accumulates over many millennia can be constrained, for example with 36 Cl data.(d) The background differential stress and viscosity changes with depth for a given depth versus temperature regime (after Allison and Dunham 2018).(e) 1 m of coseismic slip on a brittle fault increases the differential stress on its underlying shear zone (e.g.Ellis and Stöckhert, 2004), and produces heterogenous changes of differential stress on neighbouring fault/shear-zones.The uneven across-strike fault spacing shown mimics the spacing for the faults studied in this paper.(f) 1 mm of slip on the underlying shear zone, which occurs every year if the slip-rate is 1 mm/yr on the overlying fault (recurrence interval for 1m coseismic slip-increments is 1000 years, assuming no temporal earthquake clustering), increases the differential stress on its overlying fault, and induces heterogenous changes of differential stress on neighbouring fault/shear-zones.The 1 mm/yr of slip on the underlying shear zone represents the long-term rate of interseismic slip, but the slip-rate in the immediate post-seismic phase is relatively rapid, decaying through time, yet this strain-rate transient lasts a full seismic cycle (1000 yrs).(g) Hypothetical time versus slip showing clustering of surface faulting earthquakes.10 earthquakes with 1m of slip occur in 10,000 years, but the slip is clustered in time.The relationship between slip-histories of this kind, covering many millennia, and the stress histories implied by panels (a) to (f) are poorly known and form the subject of this paper.
1977), we have lacked examples with quantitative constraints on the timing, duration and magnitudes of fault slip and strain at depth, measured over timescales including multiple surface faulting earthquakes.However, in the study of Faure Walker et al. (2012), strain-rates were calculated by averaging the slip across faults scarps that have formed since the demise of the high periglacial erosion rates during the last glacial maximum (15 kyrs ± 3 kyrs; Piccardi et al., 1999;Roberts and Michetti 2004;Papanikolaou et al., 2005;Faure Walker et al., 2010) (e.g. Figs. 3-5).This time period includes multiple surface faulting earthquakes (e.g.Fig. 1) with surface slip of up to a few tens of metres accumulating on individual faults (e.g.Galadini and Galli 1999;Galadini and Galli 2000;Roberts and Michetti 2004;Galli et al., 2008).The strain-rates calculated by Faure Walker et al. (2010) and Faure Walker et al. (2012) have been summed between faults located across strike The over-steepened slopes have dips that are a few degrees less than the dips of the preserved fault planes that occur at their bases.The triangular facets are the eroded remnants of parts of the fault plane that are older than those parts sampled for 36 Cl, and erosion has produced a lower dip angle.The hangingwall is formed of colluvial slope deposits sourced from the eroding triangular facets.The sample sites all fall within un-eroded portions of the hangingwall lower slope, indicated by the lack of eroded gullies within the red dashed lines.This ensures that sampled fault planes have been exhumed by fault slip rather than by erosion.(d), (e) and (f) Scarp profiles at the 36 Cl sites showing offsets since the stabilsation of the periglacial slopes.The scarp profiles are sited at the pink dots in (a), (b) and (c).These slope values are input values for the modelling to recover slip histories.The profiles were constructed by measuring sequential changes in dip angle using a 1 m rule, combined with information from terrestrial LiDAR scanning.(g) to (l) Photos of the preserved fault planes at the sample sites.The faults planes retain their planar geometry and are relatively free from erosion indicated by preserved mm-scale striations on the slickensides.The positions of trenches are indicated and samples for 36 Cl were taken from within and directly above these trenches.(m), (n), (o) and (p) show density measurements from the colluvium in the trenches and densities recovered from MCMC modelling of the 36 Cl data.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)from each other within 5 × 90 km boxes (Fig. 2a), and compared with SRTM topographic data from the same areas.Cowie et al. (2013) averaged the strain-rates from Faure Walker et al. (2012) for 50 km along strike in 5 km along-strike steps, and plotted this against topographic data averaged in the same manner (Fig. 2b).As described above, the resultant correlation between strain-rates and topographic elevation exhibits a power-law relationship where the elevation exponent, n, is equal to ~3.26 (with similar exponents derived if data were averaged for various values between 5 and 60 km along strike; Cowie et al., 2013).Note that elevation differences contribute to the vertical stress and  Cowie et al. (2017).Zero indicates the location of the base of an organic-rich layer identified in the trench which corresponds to the base of deposits from the warmer climatic phase following the demise of the periglacial conditions during the last glacial maximum (LGM).Also shown are the least squares solution and median and 90% confidence from the full posterior distribution derived from the MCMC modelling.(d), (e) and (f) Slip histories derived through modelling of the 36 Cl data.The least squares solution is shown (blue line) along with density plots for the top 2000 least squares solutions.Red lines show the 90% confidence from the posterior distribution.Note that the total slip is relatively small for the Maiella site, and we have less samples compared to the Scanno and Pescasseroli sites; to counteract the trade-off between sample number and slip magnitude that can impair visualisation from such 2d histograms, the inset labelled (d) i shows the top 100 least squares solutions for the Maiella site with an alternative set of pixel sizes and a peak density of 60 models per pixel; this visualisation improves visualisation of a slip-rate change.The vertical arrow labelled "surface samples" shows the extent of the exposed fault sampled for 36 Cl (see text for how this affects the derived slip history).More details of how the slip histories were recovered, and related interpretations are given in Fig. 5 and the text.For this figure, our interpretation is that the slip-rates appears to be variable through time at all 3 sites, with 2 periods of higher slip-rate on the Pescasseroli and Scanno fault and a single period of high slip-rate on the Maiella fault.Periods of high slip-rate occurred at different times on different faults.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)hence differential stress at depth, hence as described above (Cowie et al., 2013), this resembles the power-law relationship for dislocation creep, where strain-rate, ε, and differential stress, σ, are related by ε ∝ σ n , with n = ~3 (Hirth et al., 2001).Thus, it is implied that the slip-rates on the brittle faults are driven by the viscous strain-rates on underlying shear zones, because the strain-rates appear to be controlled by differences in elevation and hence differential stress.
In detail, slip on the brittle faults occurs episodically, dominated in terms of strain accumulation by damaging earthquakes (Mw 5.8-6.9;Fig. 1a).In contrast, at depth the strain is presumably accommodated by quasi-continuous viscous processes due to the increases in pressure and temperature with depth, but with the viscous deformation driven by elastic deformation of the deforming material.This visco-elastic deformation is widely accepted and replicated in relatively high-temperature rock deformation experiments, where material is heated, deformed by hydraulic rams, and then deforms in a viscous manner to dissipate the elastic strain and associated stresses (e.g.Rutter and Brodie 2004).In other words, given relatively high temperatures, stress increases related to elastic strains imposed on rocks at depth are dissipated through time by viscous flow.The depth of the brittle to viscous transition in central Italy is debated and may even fluctuate during the seismic cycle (Rolandone et al., 2004), but probably occurs at a depth of ~12-15 km (Fig. 1b).The geometries of the faults and shear zones at depth is uncertain (Boncio et al., 1998;Barchi et al., 2003;Finetti et al., 2001;Collettini et al., 2006), but mainshock and aftershock data suggest that at least some faults are relatively planar to depths of ~10-15 km (Chiaraluce et al., 2011), and underlying shear zones are likely to share similar planar geometries, at least for a few kilometres beneath the brittle-viscous transition, as seen in other well-imaged regions of normal faulting (e.g.Reston 1990;Phillips et al., 2016).For this paper we have assumed that in the depth range of ~15-20 km shear-zones are likely to exhibit similar dips to their overlying brittle faults (Fig. 1c).
Earthquake slip on brittle faults produces changes in Coulomb stress and hence differential stress on neighbouring faults and their underlying shear zones (Fig. 2e), (e.g.Reasenberg and Simpson 1992;King et al., 1994).Such stress changes can induce changes in strain-rate on underlying shear zones because the differential stress is related to strain-rate (e.g.Hirth et al., 2001;Ellis and Stöckhert 2004;Trepmann et al., 2007;Nüchter and Ellis, 2010).For example, shear zones underlying earthquakes experience a positive change in differential stress (Ellis and Stöckhert 2004), whilst neighbouring shear zones can receive a negative change in differential stress (Mildon et al., 2022; compare shear zones A, B and C in Fig. 2e); thus the implied change in strain-rate can be an increase or a decrease depending on location.These changes in shear zone strain-rate will be transient because post-seismic viscous strain accumulation works to dissipate accumulated stress (Verdecchia et al., 2018).These strain-rate transients have durations similar to the recurrence time on the overlying brittle fault, and the finite viscous strain at depth matches the finite brittle strain in the upper crust when summed over timescales involving multiple large magnitude earthquakes (e.g.1c).Thus, strains produced by the deformation in shear zones, over timescales involving multiple large magnitude earthquakes on their overlying faults, will also induce differential stress changes on neighbouring structures, again both positive and negative (Fig. 2f).This will change the timing of failure on brittle faults, bringing the time of the next earthquake closer or further away in time depending on location, which in the longer term may produce fluctuations in brittle slip-rates and earthquake clustering (Mildon et al., 2022).It also produces either positive or negative changes in differential stress on neighbouring shear zones depending on their geometries and location (Mildon et al., 2022).The cumulative magnitude of these differential stress changes and strain-rate changes on shear zones, integrated over entire seismic cycles, or even across multiple seismic cycles, will be large (e.g.differential stress reductions of − 1.78 to − 3.42 MPa at 15-16 km depths on receiver shear zones during periods of rapid slip lasting several millennia on neighbouring structures imply reductions of strain-rate from 1.46 x 10 − 16 s − 1 to 7.71 x 10 − 17 s − 1 and 3.73 x 10 − 17 s − 1 ; Mildon et al., 2022).Furthermore, although such stresses will be dissipated by viscous processes, the associated finite strains will be proportionately large and permanent (e.g. 36Cl data reveal slip of 6.15 m in 13.5 kyrs at 0.46 mm/yr, and 8.0 m in 8 kyrs at 1.0 mm/yr, on faults that remained active, representing slip-rate enhancement factors of 1.22 and 2.13 compared to slip-rates over ~17 kyrs of 0.37 and 0.47 mm/yr respectively; Mildon et al., 2022).
In summary, there is a complex feedback process involving interactions between paired fault/shear-zone structures, that is, faults and their underlying shear zones, but also interactions between neighbouring fault/shear-zones; these interactions will produce changes in differential stress that will produce large, permanent strain-patterns whose temporal development can be studied if time and slip constraints can be gained (e.g.Mildon et al., 2022).We argue that the 36 Cl data presented in this paper, collected from fault planes located across strike from each other, provide an additional example to that in Mildon et al. (2022).
The panels in Fig. 2a-f raise several questions that it may be possible to answer using the constraints provided by the 36 Cl data in this paper.
(1) Is the strain-rate and hence slip-rate on shear-zone A in Fig. 2e likely to be constant or fluctuating through time during fluctuations in sliprate over time periods including multiple large magnitude earthquakes on its overlying fault?(2) If it remains constant, what levels of differential stress are induced on the overlying fault during times when the slip-rate on the fault is below the longer-term rate and the inter-seismic periods are relatively long, e.g.Period x in panel (g), and are such values of differential stress plausible?(3) If the strain-rate fluctuates, how does the heterogenous transfer of differential stress onto neighbouring fault/ shear zones accumulate through time and how does this influence the strain-rates and slip-rates on those structures (e.g. is the slip clustered in time as in Fig. 2g)? (4) If the neighbouring fault/shear-zones B and C also undergo clustered slip, is the magnitude of the differential stress changes induced by clustered slip on fault/shear-zone A sufficient to explain the strain-rate and slip-rate changes on fault/shear-zones B and C? The answer to question ( 4) is particularly interesting because if the answer is positive it provides new insights into the causes of temporal earthquake clustering that presently makes probabilistic seismic hazard forecasting very challenging.What is needed to answer these questions are datasets that constrain the timings, rates and geometries of the types of interactions shown in Fig. 2. As suggested above, this paper investigates these questions using an example from the Italian Apennines.

Geological background
The central Italian Apennines are undergoing north-east/south-west orientated extension at a rate of ~ 3 mm/yr based on evidence from geodetic data, Holocene fault scarps and historical earthquakes (Selvaggi, 1998;D'Agostino et al., 2011;Faure Walker et al., 2010).Extension is accommodated by a network of normal faults, which strike north-west/south-east and predominantly dip to the south-west (Fig. 1a).The present-day extension has been ongoing since 2-3 Ma (Cavinato and De Celles, 1999;Roberts and Michetti, 2004) and the normal faults cross-cut the pre-existing Alpine fold and thrust structures (Doglioni, 1993;Nijman, 1971).
What it is known is that these three faults all show offsets of pre-rift strata of ~2000 m as shown on published geological maps and in a review of such fault offsets (e.g.Roberts and Michetti 2004).This study suggests that these three faults show evidence for late Quaternary and Holocene activity (see Fig. 3) in that they all exhibit (a) long-term fault escarpments with triangular facets hundreds of metres high representing fault planes degraded by erosion (e.g Tesson et al., 2021;Tucker et al., 2020), (b) an over-steepened base to the triangular facets representing more recent fault plane emergence at the surface and hence relatively less degradation by erosion (Tucker et al., 2011), and (c) offsets of rectilinear periglacial slopes across emergent bedrock fault scarps with exposed bedrock carbonate fault planes which are generally interpreted as signs of fault activity after the demise the high erosion rate that acted during last glacial episode in the Apennines (e.g.Piccardi et al., 1999;Galadini and Galli 2000;Morewood and Roberts 2000;Roberts et al., 2002;Roberts and Michetti 2004;Papanikolaou et al., 2005;Faure Walker et al., 2010;Boncio et al., 2016).Thus, for (c), these carbonate bedrock fault scarps have been envisaged to have formed since the end of the LGM (15 ± 3 ka) due to a reduction in erosion rate relative to fault throw rate, thus allowing coseismic exhumation and offsets to accumulate and be preserved (Roberts and Michetti, 2004).The age of such faults scarps in central Italy is also constrained by in situ 36 Cl cosmogenic dating of the exposed fault planes (e.g.Palumbo et al., 2004;Schlagenhauf et al., 2010Schlagenhauf et al., , 2011;;Benedetti et al., 2013;Cowie et al., 2017;Goodall et al., 2021;Mildon et al., 2022), which in general, confirms the age of 15 ± 3 ka, with some examples older by a few millennia.
Debate exists about fault scarps in central Italy, with some suggesting that they form due to "non-tectonic exposure" due to "gravitational and landsliding movements aided by weathering and slope degradation processes" (Kastelic et al., 2017).However, rates of exhumation of the fault planes measured over a 3.4 year observation period of 2.9-25.6 mm/yr (Kastelic et al., 2017) do not explain the concentrations of 36 Cl measured on the Fiamignano, Tre Monti and San Sebastiano fault planes (Cowie et al., 2017).These concentrations lie in the range of 0.5-2.3x10 5 atoms/g at the bases of the fault planes and 4.0-13.5 x10 5 atoms/g at the tops of the sampled fault planes (Cowie et al., 2017;Benedetti et al., 2013).To produce such 36 Cl concentrations given surface 36 Cl production rates of ~48 atoms/g per year (e.g.Schlagenhauf et al., 2010), exhumation of the fault planes must occur over many millennia at rates much lower than that suggested by Kastelic et al. (2017).For example, Electronic Supplement 1 suggests exhumation since 12.4-19.7 ka when appropriate slope dips, rock densities and production rates are applied in the 36 Cl modelling code produced by Schlagenhauf et al. (2010), implying slip-rates in the range of 0.2-2.18mm/yr over these timescales.These values can be explained by tectonic slip given the ~3 mm/yr extension rate across the region measured with GPS (D'Agostino et al., 2011).On the other hand, it is challenging to explain why the offsets across the bedrock scarps are only a few metres or more (see Cowie et al., 2017;Benedetti et al., 2013) if rates as high as 2.9-25.6G.P. Roberts et al. mm/yr from Kastelic et al. (2017) are applied over many millennia.An explanation may be that we note that the data from Kastelic et al. (2017) were not collected at the exact same sites as the 36 Cl data in all cases, so "non-tectonic" exhumation may apply to their sites due to localised erosion.This is consistent with the ideas of Cowie et al. (2017) who state that 36 Cl sample sites must be picked with extreme caution to exclude non-tectonic exhumation.This is because sites that do record tectonic slip over many millennia may neighbour sites only a few tens of metres away that have experienced significant erosion; the sites in Kastelic et al. (2017) may fall in the latter category.Moreover, as will be shown later in this paper, periods of rapid slip on the 3 faults we study do not coincide in time with periods of increased erosion and landslide activity measured in Europe using cosmogenic isotopes and correlated with well-known period of heavy rainfall in the Holocene at 4.2 ka (see references in Zerathe et al., 2014), so non-tectonic exposure due to gravitational and landsliding movements, as suggested by Kastelic et al. (2017), appear unlikely for the sites we study.Overall, although we agree with Kastelic et al. (2017) that local erosion can spoil attempts to extract long-term slip-rates on faults, like many others (e.g.Bosi 1975;Armijo et al., 1991;Piccardi et al., 1999;Galadini and Galli 2000;Morewood and Roberts 2000;Roberts et al., 2002;Roberts and Michetti 2004;Papanikolaou et al., 2005;Faure Walker et al., 2010;Boncio et al., 2016;Villani et al., 2018;Iezzi et al., 2018;Faure Walker et al., 2021), we conclude that the offsets across fault scarps along faults can be extracted if suitable sites are chosen, and these sites will record slip produced by earthquakes, with slip-rates averaged since the demise of the LGM.
Neither the Pescasseroli or Scanno faults have any large magnitude historical earthquakes associated with them, based on the distribution of historical shaking records (Guidoboni et al., 2019), although relatively small ~5 km faults in the area at the SE extremity of the Pescasseroli fault ruptured in the 1984 earthquake sequence (Ms 5.8; Pace et al., 2002).Normal faults around the Maiella mountain may have ruptured in 1706 (M 6.8) and in 1933 (M 5.7) (Scisciani et al., 2000;Pizzi et al., 2010) based on the distribution of historical shaking records, but there is some debate about the exact fault that ruptured (e.g.Pizzi et al., 2010 suggest the local Mt.Porrara or Palena fault exhibit late Pleistocene-Holocene activity versus the Caramanico fault, but were unable to conclude on which of these faults ruptured in 1706 and 1933; see their Fig.3), with others suggesting the 1706 earthquake may be due to thrusting (de Nardis et al., 2011).No paleoseismic trenching has been conducted across the three primary faults of interest, although trenching has been undertaken on the Sulmona fault (Galli et al., 2015) and Cinque Miglia fault (D'Addezio et al., 2001).On these faults, the paleoseismic record extends back to approximately 9 kyrs BP and 7 kyrs BP respectively.

Selection of a suitable across strike transect across the Apennines
Our selection of a transect to study was driven by the availability of existing 36 Cl data sets and suitable fault outcrops to sample for 36 Cl, but also by the desire to study a relatively simple geometry with no antithetic faults and a small total number of faults.Datasets on 36 Cl have been published for the Pescasseroli and Scanno faults and herein we have re-modelled these data (see Cowie et al., 2017; the Scanno fault is termed the Frattura fault in that paper).A new 36 Cl dataset is presented herein for the Maiella fault.

Field data and sample collection
Once the three faults for study had been selected, the specific sample sites on outcropping fault planes for 36 Cl were assessed in an attempt to gain representative structural information about the faulting, and to exclude the possibility that erosion or sedimentation has influenced the surface exposure of the fault plane.The sample sites were selected on parts of the faults that have similar strike to the overall fault trace.Holocene activity is indicated by offsets of periglacial slopes that formed during the demise of the last glacial maximum (Fig. 3 d, e and f).The fault planes at the surface have suffered low erosion since exhumation by fault slip, evidenced by the preservation of millimetre scale striations on the slickenside surfaces.The contacts between the lower slopes and the fault planes (i.e. the hangingwall cut-offs), and the contacts between the fault planes and the upper slopes (the footwall cut-offs) are subhorizontal, rising or falling in elevation by less than a few decimetres over tens metres along strike, and are therefore sub-parallel.This geometry of cut-offs is consistent with exhumation due to fault slip and not erosion, because erosion would produce incised gullies, due to streamflow or mass wasting, that would be obvious in the field; such features are not present at the sample sites (Fig. 3).
Samples for 36 Cl analyses were collected from transects up fault planes exposed at the surface and within shallow (~1 m) trenches exposing the fault planes in the sub-surface. 36Cl production is ~48 atoms/g/yr at the Earth's surface, resulting mainly from interactions between cosmic radiation and the atmosphere, and hence neutrons, muons, and calcium atoms in the bedrock.Variation from this value occurs due elevation and hence atmospheric thickness and latitude, and production decreases rapidly over ~2 m depth in the subsurface to <5 atoms/g/yr due to shielding by surrounding bedrock or colluvium, depending on the density and chemistry of the materials involved and other less predominant production pathways (see Schlagenhauf et al., 2010 for full explanation of production pathways).Fault plane samples from the sub-surface will have 36 Cl concentrations solely due to sub-surface production whereas surface samples will have had an early history of sub-surface production followed by surface exposure due to fault slip and surface production.Earthquake slip histories responsible for exhuming the samples through the sub-surface 36 Cl production zone to the surface and then upwards on growing surface fault planes can be retrieved by modelling the production through time to replicate the measured 36 Cl data (e.g Beck et al., 2018, Section 4.4).Note that the overall slip history is constrained by the overall shape of the 36 Cl profile, the positions of samples from the trench, and the positions of samples from the fault plane exposed at the surface, but the younger part of the slip history is likely to be better constrained because on some faults the upper portions of the fault planes have been eroded and cannot be sampled (Fig. 3d, e and f).
Shallow trenches (~1 m deep) were excavated to collect sub-surface samples, but also to make sure that no signs of any infilled colluvial channels exist at the sample sites; no such infilled channels were encountered.Trenching also allows the density and chemistry of the hangingwall colluvial deposits to be measured (Fig. 3 m, n, o and p).Densities were measured from each distinct lithological layer in the trenches by weighing all the material from each layer, and using rulers to measure the sequential change of volume of the excavated trench.The measured colluvial densities are comparable to those derived from the Monte Carlo Markov Chain (MCMC) iteration of this variable presented later in the paper (Fig. 3m and o).
Samples of the limestone fault planes were taken with a rocksaw.Samples were 2 x 5 × 20 cm in size and typically have a density of 2.7 g/ cm 3 .
We also conducted structural and geomorphic mapping for distances of several hundred metres around the samples sites, augmenting these where possible with tripod-based LiDAR surveys and ground penetrating radar (see NERC data repository).These data have been interpreted to confirm that the fault planes at sample sites have been exhumed by surface faulting and not by erosion.

Processing of 36 Cl samples
The samples were prepared following the approach of Stone et al. (1996) and Cowie et al. (2017).Samples underwent initial crushing and clean-room chemistry at the University of Leeds, and were then analysed with AMS to determine the concentrations of 36 Cl in each sample at the SUERC laboratory in East Kilbride.

Modelling of 36 Cl histories
The history of surface exposure for fault planes produced by earthquake slip can be retrieved from in situ cosmogenic 36 Cl data through modelling that aims to replicate measured 36 Cl concentrations (Mitchell et al., 2001;Benedetti et al., 2002Benedetti et al., , 2003;;Palumbo et al., 2004;Schlagenhauf et al., 2010).We retrieved slip-histories for the fault scarps using the Bayesian Monte Carlo Markov Chain (MCMC) reverse jump approach of Beck et al. (2018).This code proposes earthquake slip histories, calculates the implied 36 Cl concentrations for each slip history, comparing these with the measured 36 Cl concentrations, iterating until the measured concentrations are mimicked.The code requires input data files on the rock chemistry, colluvial chemistry, palaeomagnetic field, and data on the geometry and character of the scarp including the elevation and latitude, dips of the upper slope, fault plane and lower slope, the offset across the scarp, the depth of the trench, and densities of bedrock and colluvium.The code iterates the factors affecting 36 Cl production and the slip histories, drawing slip histories to model and compare with the measured 36 Cl data from a Brownian Passage Time model of earthquake recurrence.The code ranks outcomes using both highest likelihood outcomes from the posterior distribution of results, taking full account of uncertainties on model inputs and model error, as well as using least squares comparisons.The code assesses convergence onto the final posterior distribution by running at least two parallel Markov chains, comparing them using a standard test described by Gelman and Rubin (1992).The code runs relatively rapidly, allowing several million model runs to be achieved within a few days using a standard laptop computer.Results for the 3 studied faults are shown in Fig. 4 and in more detail in Electronic Supplement 2.
A number of codes exist to model 36 Cl data from fault planes (Schlagenhauf et al., 2010;Cowie et al., 2017;Beck et al., 2018;Tikhomirov et al., 2019;Tesson and Benedetti 2019).The reasons why we have chosen the code from Beck et al. (2018) are explained in detail in Electronic Supplement 3 and summarised below.(1) The Beck et al. (2018) code deals with early 36 Cl production prior to the demise of high glacial erosion rates by proposing and iterating slip histories chosen from a Brownian Passage Time model of earthquake recurrence (Matthews et al., 2002), thus allowing the possibilities of anticlustering, clustering or constant slip-rate in the pre-Holocene time period, rather than assuming a single periglacial slip-rate or single value for pre-exposure (compare with Schlagenhauf et al., 2010;Tesson and Benedetti 2019).(2) The code iterates slip per earthquake rather than using values for slip interpreted and input by the user (compare with Schlagenhauf et al., 2010;Tikhomirov et al., 2019), or assuming constant slip magnitudes (compare with Cowie et al., 2017).(3) The code iterates a large number of factors affecting 36 Cl production, including timing of the demise of the high erosion rates from the LGM, colluvial density, spallogenic production rate, spallogenic attenuation length, muonic production rate, and muonic attenuation length (compare with Schlagenhauf et al., 2010;Cowie et al., 2017;Tikhomirov et al., 2019;Tesson and Benedetti 2019).(4) The code addresses uncertainty in muonic production at depth by assigning a prior uniform distribution between 1300 and 1700 g/cm 2 for the attenuation length for muonic production, iterating these values for each model run (compare with Tesson and Benedetti 2019).(5) As mentioned above, the code quantifies convergence using a standard test from Gelman and Rubin (1992).

Stress modelling
After the earthquake slip histories were retrieved through modelling to replicate the measured 36 Cl concentrations, we modelled the stress transfer implied by the slip histories.The stress modelling is not restricted to modelling single earthquakes as is common for Coulomb stress studies (e.g.King et al., 1994;Stein 2003).Instead, following the approach described by Mildon et al. (2022), we model the stress history that accumulated during surface slip that accumulated over numerous surface faulting earthquakes on individual faults, involving time periods of 20-30 kyrs.We model periods of rapid slip (clusters) and slower slip (anticlusters) that last several thousand years and involve surface slip of several metres or more that are too large to explain with single earthquakes.Note that the slip that we model includes numerous surface faulting earthquakes and will therefore include coseismic slip, and any post-seismic or inter-seismic slip that accumulated between earthquakes.
The codes from Mildon et al. (2016Mildon et al. ( , 2022) ) model stress on (a) non-planar faults with along-strike corrugations and bends, (b) underlying shear zones, and (c) include a code that allows the stress changes on receiver fault/shear-zones to be summed for coseismic slip on neighbouring faults, plus post-seismic and inter-seismic increments of slip on underlying shear-zones and neighbouring fault/shear-zones.The horizontal strain-rates for the shear zones are assumed to match that measured across the brittle faults at the surface, an assumption that is necessary to allow strain compatibility between the viscous and brittle crust (see Fig. 2c).We model the stress transfer to the viscous shear zones in the depth range of 15-24 km, assuming that the shear zones have the same dips as the overlying faults in that depth range.

Viscous flow calculations
For each increment of the stress history we convert changes in differential stress on shear zones into strain-rate changes using knowledge of viscous flow, and use this to drive the slip-rates on the overlying faults.In particular, we study the horizontal strain-rate between two points on either side of the fault/shear-zones, which we term the farfield horizontal strain-rate.We do this because it does not require us to know the exact thickness of the deforming shear zone.The viscous flow modelling implements the quartz flow law from Hirth et al. (2001), where strain-rate, ε, and differential stress, σ, are related by ε ∝ σ n , with n = ~3.In detail, we use the quartz flow law for dislocation creep in quartz shown in equation ( 1), where, ε is strain rate, A is a material parameter, fH 2 O is water fugacity, m is the water fugacity exponent, σ is the differential stress, n is the stress exponent, Q is the activation energy, R is the ideal gas constant, and T is absolute temperature.
Values for variables in equation ( 1) are described in detail in Mildon et al. (2022) and implemented in Electronic Supplement 4. Stress changes from the codes provided by Mildon et al. (2022) are retrieved and input into Electronic Supplement 4 to calculate implied changes in the far-field horizontal strain-rate.Through comparison with far-field strain-rates implied from the scarp offsets since ~15 ka, Electronic Supplement 4 calculates the "difference factor" between the long-term ~15 kyrs strain-rate and those implied by the changes in differential stress.This "difference factor" is then used to perturb the long-term slip-rates on the three faults in question for the next period of slip, to gain the implied slip-rates histories.
The differential stress changes are not uniform across the modelled 3D fault/shear-zones (e.g.Figs. 6 and 7) so the user must choose appropriate value(s) to input into the calculations of the viscous strainrates.We test two hypotheses for the most appropriate choice of value (s).To model periods of reduced slip-rate, we follow the approach Mildon et al. (2022) who studied fault/shear-zone interactions constrained with 36 Cl data further north in the Apennines.They input the minimum value of differential stress at 15-16 km on the underlying shear zone, that is, the greatest stress decrease, as their hypothesis is that this element of the 3D model will have the highest induced increase in viscosity (see Fig. 2d), and hence be the rate-limiting element in the system.Mildon et al. (2022) showed results that are consistent with this hypothesis, and here we re-test it by using the results in this paper.To model periods of rapid slip, we test the hypothesis that the most appropriate value is the average differential stress change along the entire fault/shear-zone length at 15-16 km produced by two earthquakes plus the intervening interseismic shear zone creep.This value may be appropriate because it includes the influence of the entire length and displacement gradients on the fault/shear-zone, allowing the inclusion of elements with stress changes that promote both more rapid slip and less rapid slip, and also includes coseismic, postseismic and interseismic strains and stresses.

Overall modelling workflow
The overall model workflow is described by Mildon et al. (2022), and can be summarised as follows.(1) We derive the history of slip on each of the three faults by attempting to replicate the measured 36 Cl concentration profile up the dip of the fault plane using the Beck et al. (2018) code; we term this the "measured slip-history".(2) We derive the history of stress changes that are implied by the "measured slip history".
(3) We extract the stress changes that occurred at the start and end of periods of rapid or slow slip.(4) We input the stress changes into the viscous flow law in Equation (1), using Electronic Supplement 4, to derive implied strain-rate changes for shear zones underlying the brittle faults.( 5) We use strain-rate changes in the viscous shear zones to drive strain-rate changes on the overlying brittle faults for the next period of slip, comparing the predicted slip-rate histories with the "measured slip-history" from 36 Cl.The overall workflow tests whether changes in differential stress on underlying shear zones imply strain-rate changes that are sufficient to explain the slip-rate changes and hence earthquake clustering on the overlying brittle faults.

Observations derived from 36 Cl data
The results of MCMC modelling the 36 Cl data to recover slip histories are shown in Figs. 4 and 5, revealing the ages of scarps and fluctuations in slip-rate.
The results show that the concentrations of 36 Cl increase from a minimum within samples excavated from the trenches up the dip of the sampled fault planes, but with different concentration profiles for each site (Fig. 4a, b and c).The results of modelling to recover the slip histories from the 36 Cl profiles are shown in Fig. 4d-f and Fig. 5, detailed in Electronic Supplement 2, and summarised below.
Firstly, the MCMC modelling indicates that evidence for slip may be preserved from as far back as at least ~30,000 yrs BP for the Pescasseroli and Maiella faults, and ~20,000 years for the Scanno fault (Fig. 4d-f).The age of periglacial slopes that stabilised during the demise of the glaciation in central Italy, and hence offset by post-LGM scarps, has in the past been suggested to be 15 ± 3 ka by other authors, so the ages of scarps are expected to be no older than this age (Roberts and Michetti 2004;Palumbo et al., 2004;Faure Walker et al., 2010, 2012, 2021;Schlagenhauf et al., 2010Schlagenhauf et al., , 2011;;Tucker et al., 2011;Benedetti et al., 2013;Cowie et al., 2017;Tesson and Benedetti 2019;Goodall et al., 2021).However, our results suggest that some evidence for the slip before 15 ± 3 ka has been preserved on the three faults we study.Preservation of fault scarps and rock outcrops in general is thought to be controlled by the interplay between periglacial erosion rates and periglacial fault slip-rates, and it is generally accepted that the former outpaced the latter in the high mountains of the central Apennines (e.g.Bosi 1975;Piccardi et al., 1999;Galadini and Galli 2000; Morewood and Fig. 6.Map view down onto the faults and shear zones in the sub-surface.Implied stress transfer to brittle faults overlying viscous shear zones assuming 10,000 years of interseismic slip on the viscous shear zones with no earthquakes on overlying faults.The value of 10,000 years is derived from time periods of zero or little slip inferred through MCMC modelling of the 36 Cl data in this paper.Differential stresses in excess of 100 MPa would accumulate on the bases of the brittle faults in 10,000 years, so it is suggested that this is an unlikely scenario.Instead,it is suggested that the rates of interseismic on underlying shear zones may vary through time, with reduced rates during anticlusters, and the magnitudes of possible rate changes and their effects on slip-rates on overlying brittle faults are explored in the following figures.The rates of interseismic slip in following figures are varied depending on the values of differential stress imposed on them by interactions between faults and their underlying shear zones and interactions between neighbouring fault/shear-zone structures (see text for more explanation).The full stress history is shown in a movie included as Electronic Supplement 5. A combination of stress changes from coseismic slip on faults and interseismic slip on shear-zones is used to calculate strain-rate changes on neighbouring shear-zones and slip-rates on their overlying faults (shown in Fig. 8).Roberts 2000;Roberts et al., 2002;Roberts and Michetti 2004;Papanikolaou et al., 2005;Tucker et al., 2011).However, it is reasonable to assume that some bedrock outcrops may have existed in places along fault scarps in the periglacial conditions that prevailed in the central Apennines, especially if they experienced a pulse of rapid slip as erosion rates started to decline during the demise of the LGM, or if frozen screes or permanent snow patches protected the scarps from erosion.This may be the explanation for the relatively old ages of the three faults we have studied.However, it is challenging to constrain this with any certainty because it is only possible to sample for 36 Cl close to the top of the scarp where the fault plane itself is preserved (see Fig. 4e-g and 5d-f).To take this into account, we concentrate our interpretations of the slip and stress modelling on the slip since ~20 ka after which all 3 of the scarps we study show 36 Cl derived evidence for the accumulation of surface fault offsets (Fig. 4d-f) and hence earthquake surface rupturing.
Secondly, the results indicate that slip-rates on the faults since ~20 ka were not constant through time (Fig. 4e-g and 5a-c; Cowie et al., 2017).There are periods of time lasting up to 10,000 years when little or no slip accumulated, interspersed with periods of time lasting a few thousand years when up to 8 m of slip accumulated.Increases and decreases in slip-rate are indicated by ( 1) steepening and shallowing staircase patterns to the least squares solution that provides the best fit to the data, (2) concave-up and concave-down inflections in the 90% confidence lines derived from the full posterior distribution, and (3) the density of models from, for example, the top 10000, 2000 and/or top 100 least squares solutions.As it is the choice of the user to give their own weights to these three indicators, interpreting the durations and amounts of slip associated with the slip-rate fluctuations is subjective, and the reader should keep this in mind.However, note that our interpretations of the slip-rate changes considered all 3 indicators.In Fig. 5a-c, blue rectangles indicate our interpretations of the durations of the phases of high slip-rate, with vertical double-headed arrows and values in metres indicating the amount of slip.
Thirdly, our interpretations suggest that periods of high slip-rate occurred at different times on different faults, and that two waves of high slip-rate migrated from SW to NE across the fault system (see numbered blue rectangles in Fig. 5).For the single least squares solutions, which are the best fits we have achieved to the 36 Cl data, a relatively high slip-rate is recorded first on the Pescasseroli fault at ~20.5 ka, then on the Scanno fault at ~17 ka with the slip-rate on the Pescasseroli fault decreasing at approximately that time.A relatively high slip-rate then occurred on the Maiella fault at ~13.5 ka with the slip-rate on the Scanno fault decreasing at approximately that time.Relatively high slip-rates occurred again on the Pescasseroli fault at ~6.5 ka by which time the high slip-rate on the Maiella fault had ceased.Relatively high slip-rates then occurred on the Scanno fault at ~3.75 ka just after the high slip-rate on the Pescasseroli fault had ceased.In addition to the least squares solution for each fault, the ensembles of the top 100-2000 least squares solutions (Fig. 5a-c) in general support the above interpretations, but suggest some uncertainty on the exact timings.Similar uncertainty in the exact timings are also suggested by 90% confidence lines from the posterior solutions (Fig. 4d-f; see Electronic Supplement 2 for more detail), which nonetheless show upward or downward inflections that are mostly consistent with the least squares interpretations.Note that the 90% confidence lines from the posterior solutions include proposed slip histories from as early as 50% of the total model runs (i.e. after a 50% model burn-in), including those after reverse jumps in the Markov chains that occurred after 50% burn-in.In other words, the 90% confidence results include results from well before the final convergence of the code onto the final solutions (indicated by a Gelman-Rubin test that compares the results of parallel Markov chains), and include results where the code has intentionally degraded the fit to the data by performing a reverse jump.It may be that the ensemble of least squares solutions provide better constraints on the timing than the posterior solutions, but of course the uncertainty illustrated by the posterior solutions is also important to keep in mind.
With the above uncertainties in mind, our interpretation (see the blue rectangles in Fig. 5a-c) is that two waves of high slip-rate moved across the strike of the Apennines, with Wave 1 from the SW (Pescasseroli fault), to the centre (Scanno fault), to the NE of the Apennines (Maiella fault), and with Wave 2 again starting on the Pescasseroli fault and then onto the Scanno fault.Below we investigate whether such waves can be replicated through stress modelling.

Stress modelling results and implied strain-rate and slip-rate histories
The approach we have used for stress modelling is to consider the effects of the entire history of slip back to ~20 ka.Initially we exclude the effect of slip-rate fluctuations measured with 36 Cl (Fig. 6; Section 5.2.1) to assess whether differential stress values implied by constant shear-zone strain rates are plausible.Then, we include the effects of strain-rate fluctuations indicated by the 36 Cl results (Figs. 7 and 8; Section 5.2.2).

Investigating whether constant shear-zone strains rates are plausible given earthquake anticlusters on overlying faults
Here we investigate whether it is plausible that the fluctuations in slip-rate on the brittle faults constrained with 36 Cl can be accompanied by constant shear-zone strain-rates, or instead whether shear-zone strain-rates are likely to fluctuate through time.This addresses questions raised in Section 2.0, such as (1) Is the slip-rate on a shear-zone likely to be constant or varying through time during fluctuations in slip-rate on the overlying brittle fault over time periods including multiple large magnitude earthquakes?(2) If it remains constant, what values of differential stress are implied on the overlying fault, and are these values plausible, during times when the interseismic period is relatively long, for example the periods with little or no slip lasting up to 10,000 years shown in Fig. 5a-c?To explore the effects of long-periods of zero or relatively low slip-rate, Fig. 6 shows the effect of 10,000 years of interseismic slip on all the supposed viscous shear zones in the region, with no earthquakes on overlying faults, using the fault/shear-zone model shown in Fig. 2. The value of 10,000 years mimics that derived from time periods of zero or slow slip inferred through MCMC modelling of the 36 Cl data (e.g.Fig. 5 a, b and c).The longer-term slip-rates on shear zones in Fig. 6 not constrained with 36 Cl in this paper are from data in Cowie et al. (2013), who present measurements of offsets across fault scarps that are interpreted to constrain slip-rates averaged over 15 ± 3 ka.The results in Fig. 6 show that differential stresses in excess of 100 MPa would accumulate on the bases of the brittle faults in 10,000 years, which seems to us to be an implausibly high value.Stresses in the upper crust may be dissipated to an extent by processes related to diffusive mass transfer and perhaps small-scale fracturing, because storage of stresses of ~100 MPa is probably unrealistic given that differential stresses are thought to be ~175-200 MPa at 15 km for 25-30 • C/km geothermal gradients (Allison and Dunham 2018).Alternatively, stress build-up may be less than ~100 MPa if the strain-rates on underlying shear zones are reduced at times when the overlying faults exhibit reduced or zero slip-rate.To differentiate between these alternatives, what is needed is a quantitative study of the implied stress accumulation given information on slip-rate fluctuations on overlying faults such as that provided in this paper from 36 Cl data.This prompted the modelling described below (Figs. 7 and 8) that examines whether implied changes in differential stress (see Fig. 2) are sufficient to explain the observed changes in slip-rate from 36 Cl data.

Investigating the effects of time-varying slip-rates on shear zone stress
Here we investigate the stress histories for ~20 kyrs, including strain-rate fluctuations implied by the 36 Cl results and the stress modelling results in Section 5.2.1.We model the shear-zone strain-rate fluctuations implied by the slip-rate fluctuations on the overlying brittle faults implied by the 36 Cl data.Uncertainty in the slip-histories is described in Section 5.1, so here we attempt to replicate a simplified  e) and (f) Periods of rapid slip interpreted from the 36 Cl data, along with slip predicted by perturbing the long-term slip-rate (averaged since either 30, 20.5 or 18 ka) using stress interactions.The stress interaction model was driven initially by inputting the slip implied by 36 Cl data for the Pescasseroli fault between 20500 and 16900 years BP and calculating the implied stress changes on the other fault/shear-zones and their implied strain-rates and hence slip-rates via the quartz-flow calculator (Electronic Supplement 4).Subsequently, earthquake slip increments implied by the interpretation of the 36 Cl data, only for faults found to be slipping, were input into the stress model, and the values of differential stress increase were recorded on shear zones beneath that actively slipping fault, and on neighbouring fault/shear-zones that were not slipping according to the 36 Cl data.The stress changes were used to calculate strain-rate and slip-rate changes on all fault/shear-zones.(g) to (l) Values of differential stress change at 15-16 km depth on the underlying shear zones for fault/shear zones that are not slipping.Values are not shown for times when a fault is slipping, or when no faults are slipping.9 out of the 10 examples shown exhibit decreases in differential stress during the time when a fault/shear-zone across strike is slipping.(k) Measured time of the next period of rapid slip versus predicted time of the next period of rapid slip from the stress transfer model combined with the quartz flow calculator (Electronic Supplement 4).The red line indicates x = y for comparison with the regression line (dashed blue line).The overall positive correlation with predicted earthquakes falling within ±3404 years of the central time of each cluster interpreted from the 36 Cl, suggests that clusters initiate in the stress/quartz-flow modelling at approximately the correct times, giving some confidence in the modelling approach.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)version of the slip-history shown in Fig. 8 d, e and f (more detail is given in Electronic Supplement 5).We perturb the long-term slip rates using calculated stress changes, showing uncertainty in our results by perturbing the long-term slip-rates calculated using a range of values for the scarp preservation age (18 ka, 20.5 ka and 30 ka).We calculate the components of stress transfer (Fig. 7) implied by both (a) slip on underlying viscous shear zones, and (b) slip on overlying brittle faults.To show the cumulative effects, we sum the stress transfer from source fault/shear-zones over several metres of slip.We do this because although stresses may be dissipated on receiver shear zones by viscous deformation, this dissipation produces strain associated with the viscous shear-zones that will be accompanied by slip on overlying brittle faults.It is this slip recorded by the 36 Cl data that we attempt to replicate, so stress dissipation by viscous slip is included in our modelling.As mentioned above, we do this by assuming that the strain in the viscous crust matches that in the brittle crust to maintain strain compatibility (see Shimamoto and Noda 2014; Allison and Dunham 2018; see Fig. 2c).
Fig. 8a-c shows the "measured slip-histories" constrained with 36 Cl, whilst Fig. 8d-f shows slip-histories implied by the stress transfer model derived by sequentially inputting the developing stress into the viscous flow law as values of differential stress change through time as shown in Fig. 8g-l.
Firstly, the results show that it is not adequate to only model coseismic stress transfer (Fig. 7).For example, beneath our Pescasseroli 36 Cl sample site, interseismic viscous slip on the Maiella shear zone produces an increase in differential stress on the Pescasseroli brittle fault, but a negative differential stress change at depth on the Pescasseroli shear-zone (Fig. 7e).In contrast, slip on the brittle portion of the Maiella fault/shear-zone (Fig. 7f) produces negative differential stress change across most of the Pescasseroli fault/shear-zone, with the exception of the top 3-4 km where the stress change is slightly positive.Similar complexity is evident for the Scanno and Pescasseroli faults during slip on fault/shear-zone structures across strike (Fig. 7 a, b, c and  d).This complexity shows that it is insufficient to only model stress changes from coseismic stress transfer, but instead it is desirable to gain the integrated value of stress change implied by both brittle and viscous slip, because this integrated value for stress change will be related to the cumulative strain that we constrain with the 36 Cl data.Below we described these integrated values through time, comparing the implied slip histories with the slip histories constrained with 36 Cl.
Secondly, the results show that the magnitudes of measured slip can be replicated by our stress modelling approach (Fig. 8 d,e,f), and this is considered to be the most significant result in this paper.In detail, the stress interaction model was driven initially by inputting the slip implied by 36 Cl data for the Pescasseroli fault between 20,500 and 16,900 years BP and calculating the implied stress changes on the other fault/shearzones and their implied strain-rates and hence slip-rates via the quartz-flow calculator.This is why the implied and measured slip-rates are identical for this time period in Fig. 8d.Subsequently, earthquake slip increments implied by the interpretation of the 36 Cl data, were input into the stress model for faults found to be slipping at specific times, and the values of differential stress change were recorded on shear zones beneath that actively slipping fault, and on neighbouring fault/shearzones that were not slipping according to the 36 Cl data.The stress changes were used to calculate strain-rate and slip-rate changes on all fault/shear-zones.We also show values of differential stress change at 15-16 km depth on the underlying shear zones (Fig. 8g-l; Electronic Supplement 6 shows more detail).Values are not shown for times when a fault is slipping, or when no faults are slipping.

Key results from the stress modelling
The key results from the stress modelling are as follows.
(1) The measured slip is mimicked in terms of slip magnitude and hence slip-rate by the slip predicted by the modelling of stress, viscous flow and sliprate (compare Fig. 8a-c with 8d-f).The predicted slip at the end of the time period (i.e. at the present-day) is within ~29% of the measured slip from the ~20 ka and 18 ka results that are probably better constrained by the 36 Cl data.The absolute differences in metres of slip increase through time in the model run, probably reflecting the incremental accumulation of errors in accumulated slip introduced by uncertainty on input values and the simplified nature of the modelling, including the values of stress change input into the viscous flow calculations (see the two hypotheses for these input values described in Section 4.5.2). (2) In general, the value for differential stress change decrease becomes increasingly negative through time when a fault across strike is slipping (Fig. 8g to l).The only exception is the period 3700-2500 years BP on the shear zone beneath the Maiella fault.However, for all cases where faults enter a period of reduced slip-rate after a period of high slip-rate, those times are marked by decreases in the differential stress on their underlying shear-zone in the range of 0.3-4.2MPa (Fig. 8g-l).(3) We can assess how well the model explains the start times of periods of rapid slip.We do this by taking the strain-rate on a shear zone at the start of an anticluster and forward-projecting to the time when enough shear zone slip is implied to produce a 1m slip event on the overlying brittle fault, equivalent to a ~Mw 6.5 earthquake, and then we compare this with when 36 Cl data imply that a new period of rapid slip initiated on that fault.Fig. 8k shows there is a positive correlation between measured and predicted times of renewed slip (R 2 = 0.88), and a Spearman Rank study shows that there is a significant positive correlation (p value of 0.019 is less than the α value of 0.05, so we can reject the null hypothesis, and accept the alternative hypothesis that the observed correlation is statistically significant; see Electronic Supplement 7 for details).The predicted values fall within ±3404 years of the central time of each cluster interpreted from the 36 Cl.Given the fact that renewed surface slip could be with coseismic slip of ~0.2-2.0 m given the lengths of the faults and uncertainty in scaling between coseismic slip and fault length (Wells and Coppersmith 1994), yet we impose 1.0 m in this calculation, it is not surprising that our predictions of the timing of renewed slip does not match the measured values precisely.However, it is encouraging to note that the model does not predict initiation times for new periods of rapid slip that are radically different to those measured with 36 Cl (mean value of predicted cluster initiation times is within 3404 years of the value implied by the 36 Cl data, with a standard deviation of 2580 years), giving further support to the hypotheses underlying our modelling approach (see Section 4.5.2).
Our results where the modelled slip history mimics the measured slip history (compare Fig. 8a-c with 8d-e), suggests that our study appears to be capturing the first order processes involved, and hence goes some way to constraining why periods of rapid slip and slow/no slip occur; this is the main finding in this paper.The first order control on slip-rate appears to be the decrease or increase in differential stress experienced by shear zones underlying brittle faults produced by interactions between neighbouring fault/shear zones.
The results also provide insights into how strain is partitioned between different structures to accommodate the regional extension, and how this partitioning may change through time.Fig. 9 shows a summary of the migration of rapid slip pulses, plus a calculation of the contribution of each fault to the regional heave-rate at different times since ~20 ka.The regional heave rate averaged over 20 kyrs along the transect studied is 0.68 mm/yr, constrained by the timing and total slip constrained with 36 Cl, and summed across the 3 faults, which are the only known major faults on this transect.During the 5 clusters that occurred since 20 ka, slip on individual faults was sufficient to approach or exceed the regional heave-rate (between 44% and 496% of the long-term regional extension rate).At other times no slip occurred on the transect, suggesting that the regional heave-rate on the major faults was zero (e.g. during the time periods of 10.0-6.5 ka and 2.5-0.0 ka).Together the (1) short periods of time when the heave-rate exceeded the regional heave-rate averaged since ~20 ka, and (2) the time periods with no measured heave, combine to produce the heave-rate averaged since ~20 ka, implying time-varying regional extension.Compared to extension rates measured with GPS over the last few decades, between 4000 and 2500 years BP, a single fault (e.g. the Scanno fault) produced a heaverate equivalent to the geodetic extension rate, which lies in the range of 2-4 mm/yr (Anzidei et al., 2001;D'Agostino et al., 2008;D'Agostino et al., 2011;Serpelloni et al., 2005;Galvani et al., 2012;Devoti et al., 2011;Cenni et al., 2012;Farolfi and Del Ventisette 2016).Again, this implies time varying extension rates.As strain-rates are used to infer both the mechanical process involved in continental deformation (e.g.England and Molnar 1997), and underpin calculations of seismic hazard (Console et al., 2008), these results emphasise how important it is to recognise and understand the causes of earthquake clustering, because our results suggest that it is a simplification to assume that the present-day spatial distribution of strain-rate constrained by geodesy is stable through time.

Discussion
Our results show that interactions between neighbouring fault/ shear-zones produce changes in differential stress that in turn produce changes in strain-rate on underlying shear zone and hence slip-rate on overlying brittle faults.The results show that the start times and slip-rates within measured clusters and anticlusters can be replicated by calculating the effects of these interactions.This result is consistent with the results of Mildon et al. (2022), but that paper only attempted to explain the slip-rates and start times of anticlusters, whereas herein we explain slip-rates and start times of anticlusters and clusters.
Changes in the deformation on shear zones has been previously suggested to explain why rapid slip associated with an earthquake cluster ceases, in that strain hardening due to microstructural evolution on the shear zone at depth is suggested to outpace annealing (Dolan et al., 2007(Dolan et al., , 2016;;Dolan and Meade, 2017).This explanation suggests that the location that is affected by a new period of rapid slip will be the location with the most annealed and hence weakest shear zone.Given the complexity of the hardening and annealing/weakening processes envisaged, the new location may be difficult to predict, but the advantage of this explanation is that it does not rely explicitly on interaction between neighbouring fault/shear-zones, allowing potential explanation of why earthquake clustering migrates between structures many tens to hundreds of kilometres apart.Our results augment the suggestions by Dolan et al. (2007Dolan et al. ( , 2016) ) and Dolan and Meade (2017) adding further insights into clustering on faults with relatively close spacing.It Fig. 9. Summary of the migration of periods of rapid slip across the strike of the Apennines, and comparison with strain-rates implied by GPS data, the summed rate averaged over ~20 ka constrained by 36 Cl data, and the summed rate assuming the total offset across the scarps accumulated in 15 ± 3 ka.(a), (b) and (c) Five periods of rapid slip have been measured and replicated via modelling.(d) Locations of the periods of rapid slip through time and the approximate width across strike that is implied to experience ~0.5 g horizontal accelerations, assuming this extends 20 km into the hangingwall and 10 km into the footwall of a fault that ruptures in an earthquake.Fault dips are from field measurements of the surface outcrop of the fault plane.(e) Slip-rates within each period of rapid slip-rate converted into heave rates using the measured fault dip.Assuming a 4 mm/yr regional extension rate the faults have individually taken up between 8% and 84% of the regional extension.Assuming a 2 mm/yr regional extension rate the faults have individually taken up between 15% and 169% of the regional extension rate.Assuming a 0.68 mm/yr extension rate, constrained by 36 Cl, the faults have individually taken up between 44% and 496% of the regional extension rate.may be that the interaction-induced changes in differential stress and strain-rate we envisage may well work in addition to the competing rates of hardening and annealing envisaged by Dolan et al. (2007Dolan et al. ( , 2016) ) and Dolan and Meade (2017), perhaps with annealing promoted by the lowered strain-rates that accompany reductions in differential stress on shear zones.
Our suggestion is similar to that advanced to explain post-seismic slip following a single earthquake (e.g.Ellis and Stöckhert, 2004;Trepmann et al., 2007;Nüchter and Ellis, 2010), but here we use the stress interactions between faults and shear zones not only for single earthquakes, but also to explain slip-rate changes that have been measured with 36 Cl over millennial timescales containing multiple earthquakes, and between neighbouring fault/shear-zone structures.
The implication of our suggestion is that, at least for neighbouring structures that are close enough to interact through static stress transfer, the timings and locations of earthquake clusters should be quantifiable to an extent, given knowledge of the fault geometries, the conditions operating across the entire system of neighbouring fault/shear zones, and the long-term background deformation rates.For example, Fig. 10 shows that the spatial changes in differential stress induced by slip on a single structure are systematic with distance and azimuth, and hence quantifiable.We have modelled the differential stress transfer across a hypothetical system of thirty-five normal fault/shear-zones, each with a dip of 60 • and length of 25 km, with the brittle viscous transition at 15 km depth.We show a model with four Mw 6.6 earthquakes on fault eighteen, each with 1 m of coseismic slip at the surface, with three intervening periods of interseismic slip lasting 1000 years, mimicking a period of relatively rapid slip averaged over the four earthquakes.We include stress changes induced by coseismic slip on the brittle fault and interseismic slip on the underlying viscous shear zone.We show the spatial variations in differential stress change induced on neighbouring faults, revealing how values vary with both distance and azimuth relative to the ruptured fault (Fig. 10 ci-v).We have also calculated the same for up to eight earthquakes, again with interseismic slip and loading lasting 1000 years between each earthquake and shown the resulting change in differential stress for one point at 15.5 km depth at the centre of the trace for fault/shear-zone nineteen.The change in differential stress is a linear function of the number of earthquakes (Fig. 10di).The differential stress change of − 3.5 MPa on the shear zone beneath fault nineteen induced by eight earthquakes on fault/shear-zone eighteen implies that the slip-rate on fault nineteen decreases by a factor of 0.25, given the calculations of strain-rate variations implied by changes in differential stress implemented in Electronic Supplement 4. In other words, overall, Fig. 10 implies that the slip-rate reduction on any of the receiver fault/shear-zones due to a cluster on a neighbouring fault/ shear-zone can be calculated given inputs on the kinematics, magnitudes of slip, and geometry of the fault system; we term this distance and azimuth-related change in differential stress the "differential stress transfer function".This is similar to the so-called "system effect" proposed by Visini and Pace (2014), who calculate the effect of Coulomb stress changes from the neighbouring faults in a system, but here we include the effects of slip on underlying shear-zones in addition to the effect of the brittle faults.The "differential stress transfer function" changes for different slip magnitudes in a cluster on the fault undergoing rapid slip (Fig. 10d), and will also change for different values for dip of the fault/shear-zones, depth of the brittle-viscous transition, and other changes in fault geometry such as non-planar along strike bends in faults.However, given information on such variables, we argue that the differential stress transfer function will be systematic and hence it ought to be possible to calculate how the slip-rate on any particular fault/shear-zone in a region will be perturbed by a period of rapid slip on another fault/shear-zone within that region.In other words, the magnitudes and timings of slip-rate changes on faults ought to be quantifiable to some extent (Fig. 11), and it should be possible to utilise this information when assigning probabilities of rupture to faults in a region for seismic hazard purposes even if slip-rates vary through time due to temporal clustering (e.g.compare with Console et al., 2008;Pace et al., 2014).This may well be much less challenging than attempting to assign rates of fluid influx to fault zones, shear-heating, frictional parameters, or rates of strain-hardening and annealing to shear zones, if those are the processes that control temporal earthquake clustering (e.g.see Dolan and Meade 2017 for a review).In summary (Fig. 11), the aspects of clustering that appear to be quantifiable given the approach advocated in this paper appear to be (1) the slip-rate that applies across clusters and anticlusters, and (2) the duration of clusters and anticlusters.Aspects that appear unpredictable are the inter-earthquake times within a cluster, as this will clearly be influenced by the complex nature of the interacting frictional-elastic system (McCloskey and Bean 1994).
Nonetheless, a problem with our assertion that the magnitudes and timings of slip-rate changes on faults ought to be quantifiable to some extent, is that we have only considered a sub-set of all possible variables in our study of the Pescasseroli, Scanno and Maiella faults.(1) We have not included stress transfer from faults located in both directions along strike due to a lack of data covering the last ~20 ka.(2) The role of fluids penetrating the fault zones has not been explored in this paper and it has been suggested this may be a mechanism that promotes/inhibits temporal earthquake clustering (e.g.Dolan and Meade, 2017).( 3) Other factors that we have not explored include, but are not limited to, shear heating and its effects on fluid expulsion (e.g.Dolan and Meade, 2017) and changes in frictional properties due to microstructural and mechanical evolution (Biemiller and Lavier 2017;Collettini et al., 2019).The fact that we have not explored these processes, yet we are able to explain, at least to a 1st order, our observations of temporal earthquake clustering made with 36 Cl data, with a quantified model where differential stress changes alter the strain-rate, suggests these other factors may not dominate such systems, and may be secondary processes, but clearly this needs further exploration.Furthermore, the patterns of long-term slip-rate are challenging to explain using fluids, shear heating, and changes in frictional properties because those processes would need to produce regional patterns where the fault system has (a) scaling relationships between fault length and fault slip, and fault system length and fault system slip that adhere to global values measured across many fault systems (see Cowie and Roberts 2001;Roberts et al., 2002;Roberts and Michetti 2004), (b) a quantified relationship between elevation differences and strain-rate (e.g.Fig. 2; see Faure Walker et al., 2012;Cowie et al., 2013), and (c) time-correlated changes in slip-rate on different fault/shear-zones, such as the "out-of-phase" slip-rate changes documented in this paper.It is challenging to envisage, and remains to be explained, how fluids, shear heating and frictional changes would be spatially organised across a region in a way that would produce such patterns.
Overall, our work suggests that temporal anticlusters of earthquakes occur because rapid slip on a neighbouring fault/shear-zone reduces the differential stress on the fault/shear-zone experiencing the anticluster.Temporal earthquake clusters start because the reduced deformationrate in the preceding anticluster eventually produces enough stress to cause earthquake rupture of a fault, analogous to the tortoise catching up with the hare, which in turn loads its underlying shear zone, starting a new phase of positive feedback and rapid slip.The relationship between these processes and shear zone annealing/hardening, fluid effects, frictional changes and shear heating are unclear, and should be explored.Unlike annealing/hardening, fluid effects, frictional changes and shear heating, the effect of differential stress changes described herein may be less challenging to quantify given knowledge of the geometries (Fig. 10), providing insights into continental rheology and seismic hazard.

Conclusions
We have measured slip rate fluctuations on three ~20-30 km long faults located across strike from each other in the central Italian Apennines, constrained with in situ 36 Cl cosmogenic isotopes collected from (caption on next page) G.P. Roberts et al. exposed fault planes.We identify five periods of rapid slip within the last 20,000 years, each lasting between 4500 and 1250 years, and displaying slip in the range of 8.0-2.2 m.We interpret these periods of rapid slip as clusters of surface faulting earthquakes because the magnitudes of slip are too large to be explained by single earthquakes.The clusters are separated by time periods of between 10,000 and 8500 years within which surface slip of < 1m occurred.We interpret these periods as anticlusters.We have modelled the incremental history of transfer of stress implied by the anticlusters and clusters, calculating the implied changes in strain-rate on underlying shear zones using a viscous flow law.We use the implied strain-rate changes to derive the implied sliprate changes on overlying brittle faults.We assume that the horizontal strain-rate in the viscous crust matches that in the brittle crust.We show that the predicted slip-histories for the brittle faults closely resembles the measured slip-histories constrained with 36 Cl.This match implies that a relatively simple model, where strain-rates are controlled by stress changes in the viscous crust, captures the dominant process responsible for the clustering of surface faulting earthquakes, without the need to include fluid interactions, shear heating, microstructural hardening and annealing.These latter processes may be secondary effects promoted by or independent to the stress changes.Overall, because the dominant control on stress transfer is the geometry of the interacting fault system, we conclude that fluctuations in strain-rate, crustal viscosity and seismic hazard associated with clustering of surface faulting earthquakes may be predictable to an extent.earthquakes on fault eighteen, again with 1 m slip at the surface and a recurrence interval of 1000 years within the earthquake cluster.(ii) Implied changes in sliprate given the stress changes in (i), showing that eight earthquakes in the cluster on fault eighteen implies a slip-rate reduction on fault nineteen by a factor of 0.25.Overall, this figure implies that the slip-rate reduction on any of the faults due to a cluster can be calculated given inputs on the kinematics, magnitudes of slip, and geometry of the fault system.Fig. 11.Aspects of clustering that appear to be predictable or unpredictable.A = Clusters containing non-uniform inter-earthquake times and non-uniform slip magnitudes.Earthquake timings and slip magnitudes are unpredictable due to the complexity of the frictional/elastic system (e.g.McCloskey and Bean 1994).B = Intercluster times that may be predictable given the conclusions of this paper.The start of a cluster is driven by the occurrence of an earthquake on that fault which loads its underlying shear zone inducing positive fault-to-shear-zone feedback, and hence a cascade of earthquakes; the time of this initial earthquake can be calculated from the strain-rate on the underlying shear zone produced by the differential stress history.The end of a cluster is determined by the occurrence of an earthquake on a neighbouring fault whose geometric relationship (see Fig. 10) reduces the differential stress on the initial shear zone.

Fig. 1 .
Fig. 1.(a) Location map for the central Apennines, Italy, showing locations of historical earthquakes, and faults that offset landforms or deposits formed since the last glacial maximum at 15 ± 3 ka, and are hence considered to be active faults.Historical surface faulting is shown alongside the date of the earthquake.The three 36 Cl sample sites are shown.Only three active faults are crossed by the transect A-B drawn through the three 36 Cl sample sites.(b) Mainshocks and aftershocks from the 1984 and 2009 earthquake sequences (data from within rectangles X and Y) projected along-strike onto transect A-B (data from Chiaraluce et al., 2011; Paceet al., 2002).The mainshocks and aftershocks provide information on the geometries of the ruptured faults at depth and the approximate position of the brittle-viscous transition with most seismicity limit to above 10 km, but with tens of earthquakes at depths of up to 16-18 km (see Fig.12ofChiaraluce et al., 2011).(c) Schematic cross-section showing the interpreted geometries adopted in this paper to perform stress modelling calculations.Fault dips are based on surface measurements of outcropping fault planes not the aftershock data from 1984 to 2009.(d) Projections to 15 km depth of the surface fault traces, with fault traces located in the inset.Colours show annual interseismic loading rates on the brittle faults fromMildon et al. (2019) given continuous creep on viscous shear-zones at depth at the rates defined by surface offsets of landforms or deposits formed since the last glacial maximum at 15 ± 3 ka (see Fig.2afor justification of this; the 15 ± 3 ka slip-rates are from measurements inCowie et al., 2013).These time averaged loading rates will be modulated by stress changes produced by fault interactions (see Fig.2c and d, and results in this paper).(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 3 .
Fig. 3. (a), (b) and (c) Views from Google Earth showing the geomorphology of the fault escarpments at the 36 Cl sample sites.The sample sites are indicated with pink dots.The sample sites share similar geomorphology, situated on over-steepened slopes located at the base of triangular facets formed in Mesozoic carbonates.The over-steepened slopes have dips that are a few degrees less than the dips of the preserved fault planes that occur at their bases.The triangular facets are the eroded remnants of parts of the fault plane that are older than those parts sampled for 36 Cl, and erosion has produced a lower dip angle.The hangingwall is formed of colluvial slope deposits sourced from the eroding triangular facets.The sample sites all fall within un-eroded portions of the hangingwall lower slope, indicated by the lack of eroded gullies within the red dashed lines.This ensures that sampled fault planes have been exhumed by fault slip rather than by erosion.(d), (e) and (f) Scarp profiles at the 36 Cl sites showing offsets since the stabilsation of the periglacial slopes.The scarp profiles are sited at the pink dots in (a), (b) and (c).These slope values are input values for the modelling to recover slip histories.The profiles were constructed by measuring sequential changes in dip angle using a 1 m rule, combined with information from terrestrial LiDAR scanning.(g) to (l) Photos of the preserved fault planes at the sample sites.The faults planes retain their planar geometry and are relatively free from erosion indicated by preserved mm-scale striations on the slickensides.The positions of trenches are indicated and samples for 36 Cl were taken from within and directly above these trenches.(m), (n), (o) and (p) show density measurements from the colluvium in the trenches and densities

Fig. 4 .
Fig. 4. 36 Cl measurements and slip versus time for the 3 faults derived from modelling.(a), (b) and (c) 36 Cl concentrations versus distance in the plane of the fault.Data from (b) and (c) from Cowie et al. (2017).Zero indicates the location of the base of an organic-rich layer identified in the trench which corresponds to the base of deposits from the warmer climatic phase following the demise of the periglacial conditions during the last glacial maximum (LGM).Also shown are the least squares solution and median and 90% confidence from the full posterior distribution derived from the MCMC modelling.(d), (e) and (f) Slip histories derived through modelling of the 36 Cl data.The least squares solution is shown (blue line) along with density plots for the top 2000 least squares solutions.Red lines show the 90% confidence from the posterior distribution.Note that the total slip is relatively small for the Maiella site, and we have less samples compared to the Scanno and Pescasseroli sites; to counteract the trade-off between sample number and slip magnitude that can impair visualisation from such 2d histograms, the inset labelled (d) i shows the top 100 least squares solutions for the Maiella site with an alternative set of pixel sizes and a peak density of 60 models per pixel; this visualisation improves visualisation of a slip-rate change.The vertical arrow labelled "surface samples" shows the extent of the exposed fault sampled for 36 Cl (see text for how this affects the derived slip history).More details of how the slip histories were recovered, and related interpretations are given in Fig.5and the text.For this figure, our interpretation is that the slip-rates appears to be variable through time at all 3 sites, with 2 periods of higher slip-rate on the Pescasseroli and Scanno fault and a single period of high slip-rate on the Maiella fault.Periods of high slip-rate occurred at different times on different faults.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 5 .
Fig. 5. Details of how the slip-histories were recovered and interpreted.(a), (b) and (c) Slip in the plane of the fault indicated by the least squares solutions derived from MCMC modelling of the 36 Cl data.(d), (e) and (f) Medians and 90% confidence from the posterior distributions for colluvial density, spallogenic production, slow muon capture and age from which production is preserved (Tinit), and the values for the least squares solution.In (a), (b) and (c) blue rectangles show an interpretation of periods of rapid slip.Numbers and arrows show our interpretation of howperiods of rapid slip migrated between the 3 faults through time.The exact values for the time duration of the periods of rapid slip and the value of slip are subjective (see also Fig. 4), and this is discussed in the text.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 7 .
Fig. 7. Examples of differential stress change distributions produced by either interseismic slip on underlying shear zones or coseismic slip on overlying faults.(a) and (b) Interseismic and coseismic stress changes induced by deformation on the Pescasseroli fault/shear zone.(c) and (d) Interseismic and coseismic stress changes induced by deformation on the Scanno fault/shear zone.(e) and (f) Interseismic and coseismic stress changes induced by deformation on the Maiella fault/shear zone.The full stress history is shown in a movie included as Electronic Supplement 5. A combination of stress changes from coseismic slip on faults and interseismic slip on shear-zones is used to calculate strain-rate changes on neighbouring shear-zones and slip-rates on their overlying faults (shown in Fig.8).

Fig. 8 .
Fig. 8. Plots showing the relationships between the measured and modelled slip-rate on the fault/shear-zone that is slipping and stress changes on fault/shear-zones across strike that are not slipping.1-5 indicate the sequential times of periods of rapid slip.(a) and (b) -Top 2000 least squares solutions and the least squares solution derived from the 36 Cl data.(c) -Top 100 least squares solutions and the least squares solution derived from the 36 Cl data.(d), (e) and (f) Periods of rapid slip interpreted from the 36 Cl data, along with slip predicted by perturbing the long-term slip-rate (averaged since either 30, 20.5 or 18 ka) using stress interactions.The stress interaction model was driven initially by inputting the slip implied by36 Cl data for the Pescasseroli fault between 20500 and 16900 years BP and calculating the implied stress changes on the other fault/shear-zones and their implied strain-rates and hence slip-rates via the quartz-flow calculator (Electronic Supplement 4).Subsequently, earthquake slip increments implied by the interpretation of the 36 Cl data, only for faults found to be slipping, were input into the stress model, and the values of differential stress increase were recorded on shear zones beneath that actively slipping fault, and on neighbouring fault/shear-zones that were not slipping according to the 36 Cl data.The stress changes were used to calculate strain-rate and slip-rate changes on all fault/shear-zones.(g) to (l) Values of differential stress change at 15-16 km depth on the underlying shear zones for fault/shear zones that are not slipping.Values are not shown for times when a fault is slipping, or when no faults are slipping.9 out of the 10 examples shown exhibit decreases in differential stress during the time when a fault/shear-zone across strike is slipping.(k) Measured time of the next period of rapid slip versus predicted time of the next period of rapid slip from the stress transfer model combined with the quartz flow calculator (Electronic Supplement 4).The red line indicates x = y for comparison with the regression line (dashed blue line).The overall positive correlation with predicted earthquakes falling within ±3404 years of the central time of each cluster interpreted from the 36 Cl, suggests that clusters initiate in the stress/quartz-flow modelling at approximately the correct times, giving some confidence in the modelling approach.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 10 .
Fig.10.Cumulative differential stress changes across a grid of 35 normal fault/shear-zones, each with a dip of 60 • , given (a) four Mw 6.6 earthquakes on fault eighteen (each with 1 m coseismic slip at the surface) with a recurrence interval of 1000 years, and (b) interseismic loading from the shear-zone beneath fault eighteen at rate that would produce 1 mm/yr of slip on the overlying fault.(a) Map view.(b) Oblique view looking horizontally towards the NNW.(c) Graphs sampled at a depth of 15.5 km showing the stress change at the centre of the along strike fault trace and the average stress change at 15.5 km depth along Transects T1 through T6.(d) (i) Calculated values of differential stress change at 15.5 km depth for the centre of the fault trace of fault nineteen, given between 1 and 8 earthquakes on fault eighteen, again with 1 m slip at the surface and a recurrence interval of 1000 years within the earthquake cluster.(ii) Implied changes in sliprate given the stress changes in (i), showing that eight earthquakes in the cluster on fault eighteen implies a slip-rate reduction on fault nineteen by a factor of 0.25.Overall, this figure implies that the slip-rate reduction on any of the faults due to a cluster can be calculated given inputs on the kinematics, magnitudes of slip, and geometry of the fault system.