Determining Histories of Slip on Normal Faults With Bedrock Scarps Using Cosmogenic Nuclide Exposure Data

Cosmogenic exposure data can be used to calculate time‐varying fault slip rates on normal faults with exposed bedrock scarps. The method relies on assumptions related to how the scarp is preserved, which should be consistent at multiple locations along the same fault. Previous work commonly relied on cosmogenic data from a single sample locality to determine the slip rate of a fault. Here we show that by applying strict sampling criteria and using geologically informed modeling parameters in a Bayesian‐inference Markov chain Monte Carlo method, similar patterns of slip rate changes can be modeled at multiple sites on the same fault. Consequently, cosmogenic data can be used to resolve along‐strike fault activity. We present cosmogenic 36Cl concentrations from seven sites on two faults in the Italian Apennines. The average slip rate varies between sites on the Campo Felice Fault (0.84 ± 0.23 to 1.61 ± 0.27 mm yr−1), and all sites experienced a period of higher than average slip rate between 0.5 and 2 ka and a period of lower than average slip rate before 3 ka. On the Roccapreturo fault, slip rate in the center of the fault is 0.55 ± 0.11 and 0.35 ± 0.05 mm yr−1 at the fault tip near a relay zone. The estimated time since the last earthquake is the same at each site along the same fault (631 ± 620 years at Campo Felice and 2,603 ± 1,355 years at Roccapreturo). These results highlight the potential for cosmogenic exposure data to reveal the detailed millennial history of earthquake slip on active normal faults.

targets for investigating along-strike slip variation because they record a more temporally detailed history of progressive fault exposure compared to displaced landforms (Akçar et al., 2012;Benedetti et al., 2002;Cowie et al., 2017;Mechernich et al., 2018;Schlagenhauf et al., 2010). In the Mediterranean, scarps in limestone are suggested to be preserved since the Last Glacial Maximum (LGM; Armijo et al., 1992;Tucker et al., 2011), accruing slip over multiple earthquake cycles.
The rate of exhumation of the fault plane by earthquakes can be determined with measurements of the cosmogenic isotope chlorine-36 ( 36 Cl), which primarily accumulates in the scarp as a result of progressive exposure to cosmic radiation and production from abundant calcium (Ca) present in the limestone footwall (Gosse & Phillips, 2001). A forward model is required to determine fault slip rates and the pattern of exhumation through time. Normal fault scarps have a complex exposure history that starts when the fault is buried several meters below the surface, and the same profile of 36 Cl concentrations can result from different earthquake time and displacement histories as a result of attenuated nuclide production while buried. Consequently, although many previous studies suggest that the timing of individual earthquakes can be determined by this technique (Akçar et al., 2012;Benedetti et al., 2002Benedetti et al., , 2013Schlagenhauf et al., 2010;Tesson & Benedetti, 2019;Tesson et al., 2016), our primary aim is to determine fault slip rates and slip rate variations.
Cumulative fault slip can vary along strike on an individual fault as a result of (1) the natural along-strike displacement profile (Cowie & Shipton, 1998), (2) complexity of fault structure such as overlapping segments (Peacock & Sanderson, 1991), or (3) due to problems in the long-term preservation of displacement as a result of slope instability. In the case of (1) or (2), we would expect the total displacement to vary at different localities, but the timing of major slip rate changes should be temporally correlated along-strike if earthquakes typically rupture the length of the fault. In the case of (3), slip rate changes would not be correlated along strike, and slip rate may appear accelerated if material is removed from the scarp by nontectonic processes. Only one study has attempted to document the synchronicity of along-strike fault slip using cosmogenic isotopes on bedrock fault scarps (Schlagenhauf et al., 2011). They were able to model the data from multiple sites with a similar earthquake history, but only by changing the total amount of time that the scarp had been partially exposed at each site by several thousand years (2.5 ka vs. 13.0 ka, termed "preexposure"). If this parameter is kept constant between the sites, the data from different sites cannot be modeled with a temporally correlated exposure history, suggesting that the preservation of their sampling sites has been modified (supporting information Figures S1-S2).
To demonstrate the reliability of bedrock scarps for preserving earthquake and tectonic processes, we present five new 36 Cl data sets from the Italian Apennines: three localities on the Campo Felice fault and two on the Roccapreturo fault ( Figure 1). We focus on the central Italian Apennines because limestone fault scarps are common in the region and the faults are well exposed, well mapped, and easily accessible. There are 19 published 36 Cl sample sites in the region (Benedetti et al., 2013;Cowie et al., 2017;Schlagenhauf et al., 2010;Tesson et al., 2016, Figure 1). We also remodel data published by Benedetti et al. (2013) from a site on the Campo Felice fault and data published by Schlagenhauf (2009) from a site on the Roccapreturo fault, in order to directly compare with our new data on the same faults. Our sites were selected on the basis of field reconnaissance, mapping, terrestrial laser scanning (TLS), and remote sensing surveys, in order to attest that the slip preserved at the surface is primarily the result of earthquake displacements and not affected by hillslope processes. We use a Markov chain Monte Carlo (MCMC) method to model the data at each site, which constrains the timing of slip rate changes to facilitate comparison of sample localities along strike.  Palumbo et al. (2004), Schlagenhauf (2009), Schlagenhauf et al. (2011), Benedetti et al. (2013), Tesson et al. (2016), Cowie et al. (2017), earthquake moment tensors are from www.globalcmt. org, and the fault map is modified from Roberts and Michetti (2004). The regional extension direction indicated is based on D' Agostino et al. (2011). The DEM elevations are from 1 arcsecond (30 m) SRTM (Satellite Radar Topography Mission) data.
Our modeling approach incorporates uniform parameters related to early 36 Cl production at different sites from the same fault, using the timing of global climatic change as a constraint on how long the fault scarp has been preserved at all of our sampling locations. We show how our approach can be used to determine spatial and temporal variation in earthquake displacement on normal faults.

Quaternary Faulting in the Central Apennines
Global Navigation Satellite System (GNSS) measurements indicate that the central Italian Apennines is extending at a rate of 2.7 ± 0.2 mm yr −1 in a NE-SW direction (D'Agostino et al., 2011, Figure 1). This extension has produced a series of NW-SE trending normal faults that host >M w 6 surface rupturing earthquakes, which are recorded in both the instrumental and historical records (Rovida et al., 2019). Average extension rate estimates for the region based on the offset of postglacial slopes by active faults since the LGM are 3.1 ± 0.7 mm yr −1 (Faure Walker et al., 2010;Roberts & Michetti, 2004), in agreement with GNSS rates. Time variable fault slip rates and spatio-temporal earthquake clusters have been inferred in the region based on models of 36 Cl cosmogenic data (Benedetti et al., 2013;Cowie et al., 2017;Schlagenhauf et al., 2011), and several spatially correlated (along-strike) sequences of large earthquakes have occurred in the modern record (Wedmore et al., 2017).
Planar limestone bedrock fault scarps have been preserved along normal faults since the demise of the LGM (between 10 and 20 ka), and eventual reduction in hillslope erosion rates, allowing the bedrock exhumation rate, normally as a result of fault displacement during earthquakes, to exceed the erosion rate of the fault scarp ( Figure 2; Bubeck et al., 2015;Galli et al., 2012;Giraudi, 1995;Giraudi et al., 2011;Tucker et al., 2011). In the central Italian Apennines, fault scarps are observed in Mesozoic limestone, but scarps are poorly preserved where the faults pass into other lithologies, such as siliciclastic turbidite deposits. The preferential formation and preservation of fault scarps is due to the strong erosional resistance of limestone fault surfaces, and is also well documented in Greece and western Turkey, which host lithologies similar to central Italy (Akçar et al., 2012;Goldsworthy & Jackson, 2000;Mozafari et al., 2019). Exhumation of bedrock fault scarps in the central Apennines is not always only due to fault slip in earthquakes. In many areas, the footwall and hangingwall are subject to erosional and depositional processes that are currently active or have been active since the demise of the LGM. Removal or deposition of material on the hangingwall and footwall can contribute to the exhumation history of the scarp (Figure 3a; Bubeck et al., 2015).

Fault Geomorphology and Site Descriptions
We compare slip histories from multiple sites on two faults: the Campo Felice and Roccapreturo faults (Figures 1 and 4). We sampled three new sites on the Campo Felice fault and two new sites on the Roccapreturo fault, in 2013, 2014, and 2017, and we also make use of previously published data from one site on each fault (Benedetti et al., 2013;Schlagenhauf, 2009). We describe how the sites were selected, the background literature, and geomorphology relating to both faults and all sample sites.

Sample Locality Selection
Sample localities are selected to minimize the impact of post-LGM depositional or aggradational processes acting to expose or bury the planar fault scarp (Figure 3). We follow the criteria set out by Bubeck et al. (2015) and Cowie et al. (2017) to identify suitable cosmogenic sampling sites that have a stable hangingwall and footwall slope. The new sites presented in this study fulfill the following five criteria: (1) the footwall and hangingwall slopes are intact, planar, and show no evidence of incision; (2) the hangingwall slope is free of post-LGM sediments (typically associated with actively agrading alluvial fans, colluvial wedges, sloping footwall-hangingwall contacts, and the edge of major drainages); (3) the site is located away from fault relay zones (site RP2 is an exception); (4) the fault plane surface is well preserved; and (5) the contact between the free-face (fault plane) and the hangingwall slope is horizontal, ruling out along-strike mass movement.
We identify areas that conform to the first three of these criteria by investigating the contacts between the footwall, the fault scarp, and the hangingwall. Horizontal contacts at a consistent height over a distance of 10 m or more indicate a lack of significant erosion or deposition since the demise of the LGM (Figure 3).
The footwall slope should be smooth and uninterrupted by major drainages in the vicinity of the sample locality. We identify appropriate areas for sampling using a combination of satellite image analysis (Google Earth), interpretation of TLS-derived point clouds, and fieldwork.

Campo Felice
The Campo Felice fault has a total length of ∼15 km. It is composed of two overlapping segments, an ∼6 km southern section striking on average 130° and an ∼8 km northern section with an average strike of 120°. Campo Felice is part of a larger fault system (27 km) comprising the Ovindoli-Pezza fault to the SE and the Cefalone and Colle Cerasitto faults to the NW, as well as a splay that bounds the Piano di Pezza basin (e.g., Galli et al., 2008). These strands may rupture together. The Campo Felice fault and the 8 by 3 km basin it bounds have been the focus of several studies investigating the fluvial-glacial history of the basin and the Campo Felice fault (Benedetti et al., 2013;Giaccio et al., 2003;Giraudi, 2012;Giraudi et al., 2011;Wilkinson et al., 2015). The basin hosts moraine and lacustrine deposits which correspond to cool periods during global glacial cycles, including extensive moraines associated with the LGM (Giraudi et al., 2011).
Paleoseismic results from this fault system are summarized in Galli et al. (2008). The Campo Felice fault may have ruptured during an event in 1300 AD (∼720 ybp) based on data from a paleoseismic trench on the Cefalone fault segment 5 km north of the main structure, which Salvi et al. (2003) propose is linked to GOODALL ET AL.  the Campo Felice fault. Pantosti et al. (1996) undertook a paleoseismic study along the Ovindoli-Pezza fault 5-6 km south east of the Campo Felice fault; the results of which suggest that earthquakes of M 6.5-7.0 occurred sometime between 700 and 1,130 years ago, likely around 3,900 years ago, and between 5,300 and 7,000 years ago, and they estimated an average slip rate of 0.9-2.5 mm yr −1 . A TLS data set was collected along the length of the Campo Felice fault by Wilkinson et al. (2015) to investigate the Quaternary activity and geomorphology, and we used these data paired with field reconnaissance to select appropriate sampling localities.
The footwall of the Campo Felice fault is characterized by a slope that is affected by the bedding of Upper Jurassic to Upper Cretaceous carbonates. The bedding dips subperpendicular to the slope dip, and sometimes forms prominent steps in the landscape, but the slope formed during glacial periods is distinct and dips toward the hangingwall basin. The footwall slope has active drainage channels and gullies between ∼1-100 m wide that feed debris fans and gullies in the hangingwall slope. Away from active drainage, the hangingwall and footwall slopes form smooth planar surfaces (Figure 2), similar to the idealized model shown in Figure 3a. The hangingwall slope is composed of an apron of well-cemented colluvium, typical of faults in the region. The bedrock fault scarp is generally well exposed as a planar limestone breccia surface that strikes at approximately the same height across the slope (∼60-200 m above the basin floor), with a consistent displacement along strike away from drainage gullies, fans, and the tips of the segment . Preservation of the planar fault scarp varies along strike, becoming more degraded near the fault tips.
GOODALL ET AL.  In this study we present results from three sites on the Campo Felice fault that all meet our preservation criteria. Site CF1 was sampled in April 2014 and Sites CF2 and CF3 were sampled in April 2017. All sites are located on the southern segment of the fault along a ∼1 km long section ( Figure 4). The distribution of sample sites was dictated by the geomorphology of the fault; we can only sample where the site criteria is acceptable and the bedrock scarp is well preserved. Further details of each site location and characterization can be found in the supporting information, Figures S3-S8 and Tables S1-S5.
There is one preexisting 36 Cl sample site on the Campo Felice fault that was published by Benedetti et al. (2013). We compare our new data to this site and refer to it as CF4. The geomorphology at this sample site is stable and the scarp is well-preserved, as the site characteristics satisfy the criteria used for site selection outlined in our methodology.

Roccapreturo Fault
The Roccapreturo fault is part of the Middle Aterno Valley Fault system (MAVF), which has a total length of 21 km (Galadini & Galli, 2000). The fault is composed of two segments: the southern segment is ∼8 km long, and the northern segment is ∼3 km long. A 1 km long relay zone separates the two segments, with the distance between the segments varying between 400 and 900 m (Figure 4 inset). The footwall is characterized by planar slopes incised by gullies up to ∼300 m wide. The hangingwall slope is composed of forested colluvium and the bedrock footwall has low density bushy vegetation.
A paleoseismic study of the Roccapreturo fault identified two events based on the offset of stratigraphic layers dated with radiocarbon techniques (Falcucci et al., 2015). The trenches of this study are located ∼400-500 m northwest along strike from site RP2 (Figure 4). The most recent event occurred between 1879 and 2009 BP and 3787-6055 BP and the penultimate event occurred between 3787 and 4055 BP and 7329-7499 BP (reported by Falcucci et al., 2015, as 2σ age ranges). Falcucci et al. (2015) used the offset of early Pleistocene breccias to calculate a slip rate on the Roccapreturo fault of between 0.23 and 0.34 mm yr −1 . The fault has been seismically inactive during the time period covered by the historical record (approximately the past 700 years, Galadini & Galli, 2000). Schlagenhauf (2009) sampled and modeled one 36 Cl site (herein referred to as RP3) on the Roccapreturo fault. They find that the scarp did not form in one event, but multiple events of unknown number and magnitude. They suggest that the most recent event occurred approximately 2.0-3.0 ka and that the entire scarp was exhumed between 2 and 6 ka. They report the total offset in the plane of the scarp as 10.2 m and have calculated the average slip rate during the period of exhumation to be 1.7 mm yr −1 . The geomorphology of this site does not meet the necessary site selection criteria described in this paper and Cowie et al. (2017), because it is located close to a gully that appears to have contributed to exhumation of the fault scarp (supporting information, Figure S9). We remodel their data using our approach in order to demonstrate the effect of enhanced erosion on cosmogenic data from bedrock fault scarps (results identified as site RP3), comparing it with new data from localities with acceptable morphological characteristics.
We present data from two new sites on the Roccapreturo fault, the first (RP1) located 180 m northwest along strike from site RP3 and the second (RP2) a further 2.5 km northwest along strike ( Figure 4). Samples were collected in September 2013 at Site RP1 and in April 2017 at Site RP2. Site RP2 is located within the relay zone of the two strands of the fault, and because some deformation may be shared between the overlapping parts of the fault we expect the resulting slip rates to be slower than the central portion of the main strand. However, the timing of changes in slip rate (if there are changes) should coincide, as we assume large earthquakes rupture both strands of the fault. Details of site location and characterization can be found in the supporting information, Figures S10-S13 and Tables S1, S6-S7.
GOODALL ET AL.  , and CF3 were sampled during this study, and site CF4 was sampled and processed by Benedetti et al. (2013). On the Roccapreturo fault (lower panel), sites RP1 and RP2 were sampled during this study, and site RP3 was sampled and processed by Schlagenhauf (2009). The inset panel shows the relay where two strands of the RP fault overlap. Imagery from Google Earth, 2018.

Sample Collection and Preparation
Limestone fault scarps are composed of fractured limestones with an increase in fracture density into the fault core where an indurated carbonate fault gouge is present. Where unaffected by erosion, the limestone scarps have planar surfaces with slickensides and striations commonly visible on the surface. We use these indicators to identify areas where the fault plane is well preserved, because erosion will destroy fault surface features. We avoid areas of fault plane that are intensely fractured or where the scarp is eroded, as well as areas with obvious secondary precipitation of calcite. We avoid fractures and secondary calcite in an attempt to sample fault rocks that are not contaminated with vadose carbonate cements that might contain cosmogenic 36 Cl produced in the atmosphere and circulated in groundwater (Dunai, 2010).
Sampling involves excavating a trench in the hangingwall against the fault scarp to a depth of 1-2 m. At most sites, the density of the excavated colluvium is measured using a simplified version of the method outlined by Muller and Hamilton (1992), because colluvial density is a shielding parameter in the cosmogenic modeling. Using an angle grinder, discrete samples that are 15 cm wide, 5 cm high, and 2.5 cm thick are cut from the exposed fault plane along a line parallel to the slip vector on the fault (parallel to dip direction at all sites). Some samples are horizontally offset from the main vertical sample line to avoid eroded parts of the fault plane. Photos of each site including the location of the samples on the scarp are shown in the supporting information. We collect a 3D point cloud data set using TLS at each sampling site and extract the geometry of the slip parallel profile of the slope using the Matlab® code crossint (Figure 2e; Cowie et al., 2017;Wilkinson et al., 2015).
Sample preparation and measurement is undertaken following standard methods described by Cowie et al. (2017). Chemical sample preparation is conducted at the Leeds University Cosmogenic Isotope Laboratory and prepared samples are measured with the accelerator mass spectrometer (AMS) at the Scottish Universities Environmental Research Centre (SUERC). We report the 36 Cl concentration in atoms g −1 . Reported 1σ uncertainties are AMS analytical errors and include propagation of uncertainty based on procedural blanks and standard material measurements. Bulk rock chemistry is constrained by inductively coupled plasma optical emission spectrometry at the University of Leeds. Notably, sample aliquots for Ca weight % measurements must be diluted to ∼1 ppm Ca for accurate and repeatable measurements. A more detailed description of the sampling and laboratory processes, alongside full results tables enabling recalculation of 36 Cl concentrations, can be found in the supporting information. Cowie et al. (2017) show that the relationship between 36 Cl and height on a fault scarp should approximately scale with the average fault slip rate, such that faster faults have a steeper slope in 36 Cl concentration versus height and slower faults have a shallow slope. The concentrations from two sites with a similar slip history (or rate) should approximately overlap, though this may not be precisely true if there is a difference in site geometries or target element abundance (e.g., calcium). 36 Cl concentrations from a bedrock scarp do not represent direct exposure ages, because a portion of the total 36 Cl in each sample is accumulated while the fault is partially buried in the shallow subsurface. In order to model slip rate and the pattern of exhumation through time, we use a modified version of the Bayesian MCMC approach developed by Cowie et al. (2017) to explore the age-slip relationships that adequately explain the observed 36 Cl measurements within uncertainties (further described later in this section, the supporting information, and available online, github.com/lcgregory/SimpleSlips). Bayesian statistical methods are widely applied in earth science and geochronology in order to incorporate prior information and calculate the posterior distribution for a set of parameters given quantitative measurements, using a mathematical model (Bronk Ramsey, 2009;Montoya-Noguera & Wang, 2017). Bayesian inversions can also be transdimensional, meaning that the number of model parameters ("unknowns") for which we solve is allowed to vary, increasing or decreasing the complexity of the model depending on what is required by the data (Amey et al., 2019;Bodin & Sambridge, 2009;Dettmer et al., 2010;Green, 1995;Sambridge et al., 2006). Changes in slip rate can be added or removed, and the number of changes in slip rate is a hyperparameter, because it varies the number of model parameters. The number of slip rate changes is limited by a reversible-jump algorithm that favors simple solutions (Sambridge et al., 2006). Bayesian techniques are often applied to deal with uncertainty associated with limited data (Amey et al., 2019;Bronk Ramsey, 2009;Montoya-Noguera & Wang, 2017). Several different Bayesian MCMC approaches have been developed for modeling cosmogenic data from fault scarps (Beck et al., 2018;Tesson & Benedetti, 2019;Tikhomirov et al., 2011). In this study, we prefer the modified version of the approach in Cowie et al. (2017) because our primary aim is to identify and compare first-order variations in fault slip rate. While the code used here does not change some of the factors affecting production (attenuation depth and colluvial density), this code has fewer parameters than other available codes (focusing on parameters that have the most impact), and does not attempt to identify individual earthquakes, which fits within the limitations of our data.

Modeling of the Data
The MCMC code relies on the modified version of the Matlab® code from Schlagenhauf et al. (2010) to forward model the 36 Cl concentration. The forward model simulates exhumation of a normal fault plane and calculates the resulting 36 Cl concentrations including corrections for parameters such as site geometry, sample composition, and cosmogenic production. We employ a time-varying cosmogenic particle flux (from Lifton et al., 2014, using the LSD flux-based scheme) derived for each site using the most recent cosmogenic calculator implemented in CRONUScalc described by Marrero et al. (2016). Further details on site-specific production rate scaling are included in the supporting information, Table S8. Previous studies use a constant value for colluvium density at each site , and we also use a mean value of 1.5 g cm −3 .
At each sampling locality, the height of the preserved fault scarp is known ( Figure 2). There is a gap between the highest sample and the top of the scarp, because the top part has been subject to weathering processes for longer than the base and is poorly preserved. The exposure history of the unsampled portion of the scarp is modeled by assuming that the scarp has been preserved since approximately the demise of the LGM (10-20 ka). The MCMC algorithm explores a trans-dimensional parameter space, solving for both slip rate as well as the number and timing of changes in slip rate. A slip history is generated with parameters conditioned on the prior probability, to calculate a forward model of 36 Cl values for this slip history. The likelihood of the proposed slip history is calculated given the comparison between the modeled 36 Cl values relative to the measured data. The algorithm then varies one of the parameters used to define the slip history and runs the forward model again. The new slip history is accepted if it has a higher likelihood than the previous model or if the ratio of new/current likelihood is higher than a random number drawn from a uniform distribution between 0 and 1, otherwise the new model is rejected, as per the Metropolis Hastings algorithm (Hastings, 1970;Metropolis et al., 1953). We run this for 500k iterations and remove a burn-in of 50,000 iterations (from 500k) to exclude models that may be affected by the starting parameters.
The model parameters for which we solve to define a slip history are: (1) scarp age (SA, time of the first event that produced preserved fault scarp, with a 1σ normally distributed prior of 15 ± 3 ka), (2) elapsed time (ET, time since last earthquake, no prior unless something is known about the most recent earthquake), (3) timing of change points (timing of change in slip rate), (4) height on fault scarp of a change point, and (5) a hyperparameter (the number of change points). The slip rate between change points is kept constant. The actual number of parameters can vary between each iteration, dependent on how many change points are defined. In each iteration we make a small change to one of the parameters. Further details, including synthetic tests and examples from other faults, can be found in Cowie et al. (2017) and online. We use the flexible change point method of Cowie et al. (2017) rather than the fixed change point model (where the change point height up the fault scarp is fixed) because we have no additional data such as fault roughness to fix the height of the changes in slip rate up the scarp. The flexible change point model allows timing and number of changes in slip rate to vary between iterations, while the reversible-jump transdimensional algorithm naturally favors simpler models with fewer change points, potentially resolving the issue of overfitting the data (Sambridge et al., 2006).
The Bayesian MCMC algorithm results in a distribution of possible slip histories and their likelihood and misfit to the data. We then calculate the posterior probability by multiplying the likelihood by the prior. We use a constant slip size of approximately 1 m to exhume the scarp incrementally in our modeling as we find that using a smaller constant slip size has little effect on the overall model results but does make the inversion process more computationally expensive (supporting information Figure S14). If we were attempting to model individual events, then the choice of slip size would have to be further considered. We also run a suite of models at different constant slip rates, to determine whether a simple exposure history can adequately explain the data and fit the LGM hypothesis (supporting information, Figures S15-S16).

36
Cl data are plotted as cosmogenic isotope concentration versus sample height on Figure 5 for the Campo Felice and Roccapreturo faults. At each site, the 36 Cl concentration increases gradually with increasing height, due to higher parts of the scarp being exposed for longer. In general, data from different sites on each fault overlap, and the slope of 36 Cl concentration versus height is related to the slip rate on the fault, such that a steeper gradient indicates faster slip rates. Overall lower 36 Cl concentrations at 0 indicate faster slip rates. A more pronounced decay of 36 Cl with depth can indicate a longer elapsed time, because samples have spent more time partially buried and thus the pattern of measured 36 Cl more strongly reflects the decay of cosmogenic production rate with depth. 36 Cl concentration profiles are similar for sites CF2, CF3, and CF4 while concentrations are greater at site CF1 for samples from the same height. All profiles show a change in gradient at ∼2-3 m on the Campo Felice fault. On the Roccapreturo fault, site RP3 has lower 36 Cl concentrations for the same height than at sites RP1 and RP2, and has a steeper gradient at the base of the scarp compared to sites RP1 and RP2. The gradient at site RP3 gradually reduces with height. Sites RP1 and RP2 have similar 36 Cl concentrations, but with minor differences in gradient and the concentration at height 0. Site RP3 samples a section of GOODALL ET AL.
10.1029/2020TC006457 9 of 22  Benedetti et al. (2013), and data from Roccapreturo Site RP3 are from Schlagenhauf (2009). Analytical 1σ uncertainties are plotted as black lines. These are raw 36 Cl values, and some variation in 36 Cl between sites and noise in sample data is related to different production rates resulting from variable Ca, which is corrected in the forward model. preserved scarp that has an offset of 10.2 m, compared to sites RP1 and RP2 which have offsets of 7.2 and 4.7 m respectively; this difference in heights is discussed in the context of the site geomorphology and strain partitioning in the discussion. Here we present the modeling results from each cosmogenic sampling site in the context of each fault and how the results can be compared between sites on the same fault.

Campo Felice Fault
Each site was modeled with 500k iterations using the Bayesian-inference MCMC code described in the methods section, with the site characteristics listed in Table S8, and modeling parameters in Table S7 (Cowie et al., 2017, https://github.com/lcgregory/SimpleSlips). Two dimensional histograms of all accepted exhumation models for each site are shown in Figure 6. The histograms show the modeled distribution of time at which the fault surface was first exposed to the surface. The slip is modeled in approximately 1 m increments in the slip direction and is binned into 200-year intervals in the histograms. In order to compare between sites, we plot the 95 percentile range of these same exhumation histories for all sites (Figure 7a).
The models are poorly constrained above 7-10 m due to the lack of samples on the degraded part of the scarp, demonstrated by the increased variance between model results higher on the scarp. The 36 Cl concentration in each sample does reflect exposure of the fault surface at least 2 m above the sample due to being partially exposed to cosmic radiation while residing below the ground surface. As such, exposure of the GOODALL ET AL.

10.1029/2020TC006457
10 of 22 unsampled portion of the fault is somewhat constrained by the cosmogenic data at the top of the sampled portion, as well as by the independent prior in our modeling dictating that the top part of the scarp was preserved following the demise of the LGM (15 ka with a 1σ standard deviation of 3 ka). However, the older portions of the slip models have more variability in exposure time, primarily due to the range in predicted ages for the demise of the LGM and preservation of the scarp (Figure 6).
We calculate average fault slip rate based on the more probable models (or "top" 10% of models, Figure 7e), taking into account the number of slip rate changes to make sure the top models have the same distribution of slip rate changes as the full suite of models, because we do not want to introduce a bias for more or less complicated models. The models are ranked by posterior probability, separated into groups by the number of change points, and we select the top 10% of each change point group. The fit of the top 10% of CF1 models to the data is shown in Figure 8a with the corresponding exposure histories. We also show the probability of events over time for the full model distribution in Figure 7c. Figure S17 shows the fit to the data of the full model distribution for each site on the Campo Felice fault. In general, the range of accepted models fit the data well (and all fit within the standard deviation of the data), though the highest samples are poorly fit to the analytical errors by the least likely models.
For each site, we calculate the time-varying fault slip rate in mm yr −1 for the models that are the top 10% most probable (Figure 7e). Each model is a relatively simple time-slip vector, with a constant slip rate between the GOODALL ET AL.  (Figure 7e). Because one of the modeled parameters is the time elapsed since the last earthquake, each model has a period of time between the present day and the last proposed earthquake during which the incremental slip rate is zero. If another earthquake occurred today, the mean slip rate between the present day and what would then be the penultimate earthquake would change to accommodate the "new" slip, but modeled slip rates previous to the penultimate event would remain the same. Therefore, the apparent drop to 0 mm yr −1 in our slip rate calculations reflects the modeled elapsed time and does not imply that the fault is inactive-an important consideration if time-varying fault slip rates are to be incorporated into earthquake hazard assessment. Because slip rate is calculated as the mean of all of the models, we only show the rate up to 10 ka; older than 10 ka the slip rate is poorly constrained where the models "end" at the scarp age and the mean is not representative.
In general, the modeling results show agreement in exposure histories between the sites (Figures 6, 7a ,7c, and 7e). Sites CF1 and CF3 are located within 150 m of each other, as are CF2 and CF4, and there is approximately 1 km between both sets of sites ( Figure 4). There is a difference in scarp height between the group of northwestern sites (CF1 and CF3, average height 14 m) and the southeastern sites (CF2 and CF4, average height 21.5 m), but the height change does not lead to a significant difference in slip rate and the timing of change in slip rate during the past 0-7 ka, which is the best constrained time interval, because the 36 Cl concentrations at each site are all fit by an increase in slip rate during the same time interval.
GOODALL ET AL.
10.1029/2020TC006457 12 of 22  Figures (a and c) show the fit to data of models used to calculate average fault slip rates, for sites CF1 and RP1. The circles and error bars represent the 36 Cl measurements and the standard deviation of the data (used in the likelihood calculation); each colored line represents a model, with dark to light colors representing highest to lowest probability models, regulated by the scarp age probability and the number of change points. We present 400 models for each site ranging from the highest (dark blue) to lowest (yellow) probability at equal intervals (100) through the distribution. Figures (b and d) represent the corresponding model slip histories. The full distributions from the inversion for each site can be found in the supporting information.

(a) (b) (c) (d)
Models of sites CF1-CF3 have peak slip rates of 1.5-3 mm yr −1 between 500 and 2000 years, with a reduced slip rate of <1.5 mm yr −1 before approximately 3-4 ka (Figure 7e). Models of site CF4 have a higher peak slip rate of just over 3 mm yr −1 occurring more recently than at sites 1-3 (around 0.5-1 ka), reflecting the lower 36 Cl concentration and steeper 36 Cl versus height at this site, requiring a faster slip rate more recently. Models of sites CF2 and CF4 have a second longer period of increased slip rate between the demise of the LGM and ∼8 ka, though this part of the exposure history is not well resolved, based on the spread of model results at the top of the scarp in Figures 6a-6d. The results from all sites on the Campo Felice fault (Figure 7c) indicate that the fault was relatively active between 1 and 4 ka and relatively less active between 4 and 8 ka. The fault at sites CF1 and CF3 likely has less total slip in each event, compared to at sites CF2 and CF4, because the total scarp height is lower. Despite having less total slip, the timing of peak slip rate and rate change is correlated between all sites (Figure 7c and 7e).

Roccapreturo Fault
The modeling results from the Roccapreturo fault, including 2D histograms of slip versus time and mean slip rate through time, are shown in Figure 9, and the fits to the data of the 10% most probable models are shown in Figure 8. The fault slip rates in Figure 7f are calculated based on the 10% most probable models, and the probability of events over time is calculated for the full model distribution (Figure 7d). The general pattern of exhumation is characterized by relatively faster slip rate between 2 and 6 ka, and slow to zero slip between 6 ka and the demise of the LGM, and between the present to 2 ka, which implies a long elapsed time. The maximum slip rates are 2.2, 0.7, and 2.5 mm yr −1 at sites RP1, RP2, and RP3, respectively. The difference in average slip rate between sites is primarily due to the difference in scarp height, as a larger scarp requires a faster rate of exhumation averaged over the time period. The decrease from fast to slow slip rate occurs at the same time at sites RP1 and RP2 (2.5-3 ka), but RP3 has younger and more rapid peak in slip that lasts until 1.5-2 ka, required by the much lower overall 36 Cl concentration and higher scarp height at the site (Figures 5 and 9). The modeling provides a good fit to the data at sites RP2 and RP3, but data just above the ground surface at RP1 are poorly fit (Figure 8 and Figure S18). The pattern of these outlier data is not systematic, and suggests that there is additional noise that is not accounted for in some samples.

Discussion
Fault slip rates derived from cosmogenic isotopes measured on bedrock fault scarps can contribute to our understanding of fault behavior over multiple earthquake cycles and should be considered when estimating seismic hazard (Akçar et al., 2012;Beck et al., 2018;Benedetti et al., 2002;Cowie et al., 2017;Schlagenhauf et al., 2010). However, until now, it has not been demonstrated that results are consistent at different sites along strike on the same fault. Here we show that the timing of slip rate changes is similar at different sites along strike on the Campo Felice and Roccapreturo faults, but there are some differences in slip rate and total displacement between sites due to multiple factors. We discuss how the modeling parameters can be compared between sites, and highlight the assumptions and limitations of the 36 Cl method. We outline how results from sites with acceptable indicators of morphological preservation can be used to infer that both spatial (e.g., along strike on the fault) and temporal (changes in slip rate) variability is resolved on millennial timescales. We also compare our results with those from paleoseismic trenching on the same fault, which further supports the ability of cosmogenic isotopes measured on bedrock fault scarps to provide a reliable measure of fault activity.

Along-Strike Comparison of Fault Activity
The magnitude of surface displacement in individual earthquakes can vary along-strike (Ando et al., 2017;P. O. Gold et al., 2013;Rockwell & Klinger, 2013;Wedmore et al., 2019;Wesnousky, 2008), and, as a result, the pattern of cumulative slip on a fault may be temporally synchronous but spatially variable in slip magnitude. Cosmogenic analyses on bedrock scarps provide constraints on both the time and cumulative displacement, so data from multiple sites can be used to isolate the spatial variation in slip along-strike over multiple earthquake cycles. Previous 36 Cl studies on bedrock normal scarps have concluded that significant temporal slip rate variations occur on thousand year time scales Schlagenhauf et al., 2010). Faults are demonstrated to have intervals of relatively fast slip rate or short earthquake recurrence intervals, interspersed with periods of relative quiescence. Changes in slip rate are maybe linked to elastic interactions or strain partitioning processes that are larger in scale than a single fault Dolan et al., 2016) and, therefore, the timing of significant slip rate changes are likely to be temporally correlated along one fault. The alternative, that slip rate changes are not correlated along the fault, may imply that more random geomorphic processes have exposed the fault surface at different rates along the strike.
We compare the posterior probabilities of the time of scarp preservation, the time since the most recent earthquake, the number of change points in each model run, the average fault slip rate, and timing of "events" between the different sites on the same fault (Figures 7 and 10). One of the greatest uncertainties in modeling the cosmogenic data is the timing of preservation of the fault scarps, associated with the demise of the LGM and the transition from relatively fast to slow erosion of the scarp. We assign a wide Gaussian prior in our modeling to account for the uncertainty in how long it takes for fault activity to outpace erosion, and that this transition may be different for different faults. Figures 10a and 10d show the similarity in the posterior probability for scarp age at each site. We modeled the results from the Schlagenhauf et al. (2011) study using the same approach, and the scarp age posterior probabilities do not overlap ( Figure S2), suggesting that morphological factors not associated with the LGM have affected the development of scarps at those sites.
Using the total displacement measured at each site and the median value for each scarp age posterior probability, we compute the Holocene average slip rate for each site (Table 1)

(a) (b) (c) (d) (e) (f)
of the distributions and is not affected as much by a small number of very large or small outliers as the mean value. We report the median and standard deviation of all sites based on the median of the combined posterior values for scarp age (divided by scarp height) of all sites. The median slip rate over the time since the demise of the LGM is 1.15 ± 0.36 mm yr −1 at Campo Felice based on all four sites and 0.42 ± 0.14 mm yr −1 at Roccapreturo, based on sites RP1 and RP2. We exclude site RP3 from the slip rate calculation because erosion of the hangingwall has likely contributed to the exposure of the scarp.
The elapsed time parameter only has a positivity prior value assigned in the modeling, such that the elapsed time cannot be negative, because there are no paleoseismic trench sites within 2 km of our sites on the Campo Felice fault and there is no historical seismicity associated with either fault. The Campo Felice sites have a similar posterior probability distribution favoring an elapsed time of less than 800 years, with a nonnormal distribution that is skewed toward younger values (median values range between 411 and 771 years, Figure 10b, Table 1). Paleoseismic data from trenches north and south of our site on adjacent fault strands indicate that the most recent surface rupturing event on the Campo Felice fault was 720 years (north segment) and between 700 and 1,140 years (southeast segment-the Ovindoli-Pezza fault; Pantosti et al., 1996;Salvi et al., 2003). These results are in agreement with our estimated elapsed time, and they suggest that large earthquakes on the Campo Felice fault may involve multiple fault strands rupturing in the same earthquake, or sequences of events that occur over a relatively short time period, similar to several modern sequences in the region such as the 2016 central Italian sequence (Chiaraluce et al., 2017;Villani et al., 2018;Walters et al., 2018). Based on the interpreted surface displacements of 2-3 m, Pantosti et al. (1996) suggest that the causative event was a M 6.5-7.0, which is reasonable for a combined fault length of up to 35 km if all of the sampled segments were involved in one event. Pantosti et al. (1996) estimated that the average slip rate of the Campo Felice fault is 0.9-2.5 mm yr −1 , on the basis of multiple events occurring over the past 5,300-7,000 years, which also fits well with our long-term average fault slip rate (1.15 ± 0.36 mm yr −1 ).
The Roccapreturo sites have a broad distribution of elapsed time values, with the median values at RP1 and RP2 in agreement (median = 2.6 ± 1.4 ka, Figure 10e, Table 1), and a younger value for site RP3 (median = 2.1 ± 1.1 ka). These results agree with paleoseismic data that suggest the most recent event on the Roccapreturo fault was between 2 and 6 ka, with another large event occurring between 3.8 and 7.5 ka (Falcucci et al., 2015). These dates agree with the rapid slip rate between 2 and 7 ka at site RP2, which is located approximately 500 m from the paleoseismic trenches (Figures 4, 7d, and 7f). While traditional paleoseismic data have been compared to 36 Cl slip histories in previous studies, these have either been on different fault strands at distances of >5 km (Tesson et al., 2016) or have suggested disparate results (Benedetti et al., 2003;Kokkalas et al., 2007). The agreement we find between the two techniques, which have been applied in such close proximity on the Roccapreturo fault, provides further evidence for the reliability of slip histories derived from modeling of 36 Cl on bedrock fault scarps, and the potential for these two techniques to be combined for more informed seismic hazard analysis.
There are many assumptions that must be taken into account when interpreting cosmogenic data from bedrock fault scarps. Because the top of the fault scarp is not well preserved and cannot be sampled, exposure histories older than ∼8-10 ka are poorly resolved and the modeling is reliant on estimates of how the scarp is preserved through the demise of the LGM (Figures 6 and 9; Beck et al., 2018;Benedetti et al., 2013;Mechernich et al., 2018;Schlagenhauf et al., 2011;Tesson & Benedetti, 2019). The trans-dimensional nature of the Bayesian inversion favors simple slip histories with the lowest number of changes in slip rate and we do not apply any weighting to the data other than the standard deviation of each data set. The sampling bias is a challenge for calculating fault slip rates, because of the higher sample density at the base of the scarp. Consequently, the inversion favors simple slip histories that fit the data well in the bottom section of the scarp where there is a higher density of data, and that may fit less densely sampled data further up the GOODALL ET AL.  scarp poorly. Models can fit the data with a more simple slip history by not fitting the top few data points as well as the data at the base of the scarp. If the oldest part of the exposure history can be better quantified, perhaps by incorporating more sophisticated geomorphological models and data constraining the timing of the LGM (e.g., Tucker et al., 2020), the entire slip history may be better determined.
At some sites, no models fit the data to within the analytical uncertainties because they have outliers or noisy data that are not fit by any model. Applying site averaged calcium values at the two sites where we did not collect the data ourselves reduces the ability of models to fit the data because small variations in Ca concentration have a significant effect on the production rate of 36 Cl in each sample. One challenge in interpreting the output of MCMC Bayesian modeling is that, while there is a single best fit or most likely model, there are commonly hundreds or thousands of models that fit the data almost as well ( Figure 8). All of the models fit within the standard deviation of each data set and can be incorporated when calculating average slip rates and making broad interpretations. Identifying higher frequency variations or individual slip events (earthquakes) is challenging because the data can be fit with a range of models, and is not possible using the data and modeling methods in this study. However, the first-order variations in slip rate including pulses of rapid slip rate, which may represent temporal clustering of earthquakes, are consistent features in results from multiple sites along the same fault. The Bayesian MCMC approach with minimal parameters ensures that cosmogenic data are not overfit, and the result is an acceptable range of exposure histories, rather than a non-unique earthquake history.
Based on the results presented here and the large time and financial costs associated with sample processing and 36 Cl measurements, future studies may benefit from sampling multiple sites with discrete spaced samples rather than a continuous sample ladder at one sample site. The multisite sampling approach also allows more information to be gained on along strike variability of slip rates. Future research can focus on testing how many sites are required to determine a reliable slip rate, how outliers can be definitively identified, and how many samples are required for an acceptable measure of fault slip rate and slip rate variability (e.g., Beck et al., 2018). The geomorphology of each sample site should be carefully understood and documented to demonstrate the tectonic origin of bedrock fault scarps. Sampling at regular intervals on the fault scarp limits sampling bias and can reduce the complexity of interpreting modeled slip rates. While the prior that scarps are preserved only since the demise of the LGM is strongly supported in the Central Italian Apennines Tucker et al., 2011), application of the method to other regions will require equally robust evidence to define the scarp age prior distribution. Combining other data sources with the 36 Cl data, such as historical records and estimates from other dating techniques, helps to support results from cosmogenic data.

Temporal Slip Rate Variability
Temporal slip rate variability is observed at all of the sites on the Campo Felice and Roccapreturo faults (Figure 10c and 10f). Both faults experience pulses of relatively fast slip rate (over thousands of years), peaking at 3 mm yr −1 at Campo Felice and 2 mm yr −1 at Roccapreturo, separated by intervals relative quiescence with slower average slip rate. Fault slip rate variability or discrepancies between geodetic, Quaternary, and geological slip rates are observed on many faults in various tectonic settings (Dolan et al., 2016;Faure Walker et al., 2010;Ferry et al., 2011;Oskin et al., 2007;Papanikolaou et al., 2005;Zinke et al., 2017), with several mechanisms invoked to explain the variability. Orogen-scale changes in erosion patterns or the kinematics, growth, and localization of faulting may affect the comparison of geological rates (>10 5 ) with geodetic and Quaternary rates (Hoth et al., 2006;Nicol et al., 2010). In the Italian Apennines, Cowie et al. (2017) suggest that time variable slip rates are primarily caused by large scale interaction across the whole fault network, in order to minimize the work done by faults. In this geodynamic model, different regions of faults are active at different times as a result of the change in gravitational potential energy acting on the uplifted footwall, inducing flexural bending of the normal fault footwall and time-varying fault strength. Coulomb stress changes due to earthquakes are suggested to play a role in causing clustering of earthquakes and variable fault slip rates (Dolan et al., 2016;Wedmore et al., 2017). Dolan and Meade (2017) indicate that there is not yet a single mechanism that can explain this behavior across different faults, and suggest that it is caused by the complex interaction of processes that may be controlled by properties of a particular fault as well as the fault system as a whole. We observe peak slip rates at different times on the Campo Felice and Roccapreturo faults (Figure 7e and 7f), suggesting that the relatively close faults may interact in a way that one suppresses slip on the other (Figure 1). However more data from neighboring faults are required to determine whether, and how, they interact.
In order to understand the mechanism behind slip rate variability on a single fault, it may also be informative to constrain the activity of faults in the rest of the network based on observations over multiple timescales. Probabilistic seismic hazard models currently use time averaged constant slip rates on faults (Valentini et al., 2017) and have limited temporal and spatial data coverage due to the sparsity of paleoseismic data sets (Dolan et al., 2016). What can be inferred from 36 Cl data on bedrock scarps is also limited in time, but we are able to capture two major changes in slip rate at some sites, helping to better understand the variability of earthquake recurrence on timescales that are important for understanding fundamental geological problems and seismic hazard. The method can be widely applied where scarps are preserved to reveal fault interaction on thousand year timescales and to determine how several faults contribute to the large scale pattern of deformation. Slip rate variability may also be captured by quantifying slip rates using alternative methods that have different spatial and temporal coverage and resolution. Faure Walker et al. (2012) show that slip rates averaged over the Holocene (based on fault scarp heights) match the geodetic deformation rates, when averaged over large spatial scales (10's-100's of km). Cowie et al. (2013) suggest that the 1-10 ka year strain rates are representative of long-term geological rates based on the correlation between high strain and high topography, suggesting that faulting is driven by viscous flow on localized shear zones in the lower crust. 36 Cl derived slip histories have the potential to fill some of these spatial and temporal gaps and will help to elucidate the timing and mechanisms responsible for earthquake clustering and fault interaction.

Spatial Fault Complexity
The agreement between results from the Campo Felice Fault demonstrates that 36 Cl data from multiple sites spaced ≤1 km on one fault can be modeled successfully with similar slip histories. The larger fault scarp and a period of additional slip between 7 and 15 ka only observed at the two southern most sites (sites CF2 and CF4, Figures 7c and 7e) may suggest that the fault does not always rupture continuously or uniformly along strike, which matches modern observations of faults in the region (Boncio et al., 2004;Villani et al., 2018;Walters et al., 2018). Sites CF1 and CF3 are closer to the overlap between the central and northwest Campo Felice fault strands than sites CF2 and CF4 (Figure 4), and displacement may decrease as strain is shared across the two fault strands. Although Benedetti et al. (2013) determined an exposure history at site CF4 similar to our results, with two earthquakes at 1.1 ka, and events at 3.4, 4.2, 4.4, and 9.4 ka, their solution is nonunique, as there are many other exposure histories that fit the data at CF4 equally well ( Figure 6). The continuous ladder at CF4 leads to tighter constraint on our modeled parameters compared to sites CF1-CF3 ( Figure 10). Models at site CF4 also include more change points than CF1-3 (Figure 10c), suggesting more complex models (but not precise earthquake timings) can be generated from densely spaced data. This agrees with synthetic tests in Beck et al. (2018), which show that continuous sampling of the fault scarp does not necessarily resolve better constraints on absolute slip rates and the timing of change in slip rate compared to discrete sampling every 25-50 cm.
Modeled slip histories from sites on the Roccapreturo fault are not as similar as results from the Campo Felice fault. Site RP3, sampled by Schlagenhauf (2009), has a significantly larger fault scarp and a longer period of fast slip rate reaching 2 mm yr −1 from the present until 7-8 ka. The larger scarp at site RP3 is most likely a result of fast erosion of an unstable hangingwall on the edge of a major gully that incises the hangingwall and footwall of the fault (Figure 4). The fault scarp has been subject to active net erosive slope processes that likely removed material from the hangingwall slope exposing the fault surface in the gully, resulting in the higher scarp and faster slip rate than at other sites on the fault. The difference in timing of peak slip at site RP2 suggests that site RP1 experienced a more recent or larger slip event, implying that the fault does not always rupture continuously or that there is a variation in surface slip in a single event along the fault. We suggest that the slower average slip rate and shorter scarp height at site RP2 compared to site RP1 is because strain is partitioned between the end of the strand sampled (at RP2) and an overlapping fault strand located 1 km west and across strike (Figure 4). Between the two fault strands, there is a ramp in the topography that slopes down toward the southeast, and a step across each fault segment, perpendicular to fault strike, which is typical of a classic relay ramp morphology, where the length of the relay ramp is approximately three times the width (Fossen & Rotevatn, 2016, Figure 4).
Interaction between closely spaced fault segments can reduce the total displacement across individual segments due to strain partitioning, including at fault splays (Cowie & Roberts, 2001;Manighetti et al., 2015;McLeod et al., 2000). Our observation of lower slip rate where two strands overlap (Figures 4 and 7) at Roccapreturo suggests that over thousand year time scales the overlapping fault segments do not become completely inactive, but instead each overlapping segment has slower average slip rates (or less slip per event) relative to the center of the main fault segment (compare RP1 and RP2 slip rates on Figure 7f). Quaternary slip rate variation along strike is not typically observed at this scale and temporal resolution, demonstrating that 36 Cl provides a unique ability to investigate fault segment interaction and strain partitioning over millennial timescales. Due to the relatively young age of the normal fault network in the central Italian Apennines (2-3 Ma) and low extension rates across the region (2.7 mm yr −1 ; D'Agostino et al., 2011), the fault system is immature, with a complex network of faults in the region that are highly interactive on relatively short timescales (including in earthquake sequences, e.g., Nixon et al., 2016). Individual faults in the central Apennines are still growing, and through the process of localization, splays like the ones observed along the Roccapreturo fault may eventually become hard linked through to the surface and be capable of larger earthquakes and faster slip rates.
On the Campo Felice and Roccapreturo faults, we can observe the cumulative effect of the complexity of earthquake surface ruptures and resulting variation in displacement along strike. Some complexity arising during individual earthquakes may cumulatively cancel out over multiple earthquake cycles, if it is localized or random, and may contribute relatively insignificant noise to calculated slip rates. However, if patches of high or low slip occur repeatedly in the same location on the fault, the displacement at any one site is not representative of either the fault rupture as a whole during that event, or that particular site over multiple earthquake cycles. We find that the variation in slip is consistent over multiple earthquakes cycles at some sites, such as site RP2 having lower slip than at RP1, and that sites CF2 and CF4 have higher slip than CF1 and CF3. By comparing ruptures from individual events with multiple offsets accumulated over longer timescales, it is possible to better understand along-strike variability (Brozzetti et al., 2019;Cinti et al., 2019).

Conclusions
We present 36 Cl cosmogenic isotope results and modeled exposure histories from four sites on the Campo Felice fault and three sites on the Roccapreturo fault. Unlike previous work, our modeling approach can be uniformly applied to all data in order to test whether they agree, without arbitrarily varying parameters related to the preservation of the fault scarps. Models from different sites on the same fault have the same long-term preservation age and elapsed time since the last earthquake, as well as similar long-term patterns of slip rate variability. The slip histories do not agree where samples are collected from unstable slopes, and there are parts of the faults that have slower average rates (likely caused by lower displacement in cumulative events), due to rupture complexity or strain partitioning between overlapping fault strands. Each fault experiences periods of time with much faster slip rates than the average Holocene rate, likely due to multiple earthquakes occurring over a few thousand years. The average slip rate is 1.15 ± 0.36 mm yr −1 on the Campo Felice fault, and 0.42 ± 0.14 mm yr −1 on the Roccapreturo fault, with peak slip rates of 3 and 2 mm yr −1 at each fault, respectively. The range of along-strike variability in slip rate means that one location may not represent the typical behavior or hazard of an entire fault, and while sampling faults multiple times along strike is not always feasible, it can improve confidence in results by elucidating the range of slip rate and the timing of changes in slip rate.