Orogen-scale uplift in the central Italian Apennines drives episodic behaviour of earthquake faults

Many areas of the Earth’s crust deform by distributed extensional faulting and complex fault interactions are often observed. Geodetic data generally indicate a simpler picture of continuum deformation over decades but relating this behaviour to earthquake occurrence over centuries, given numerous potentially active faults, remains a global problem in hazard assessment. We address this challenge for an array of seismogenic faults in the central Italian Apennines, where crustal extension and devastating earthquakes occur in response to regional surface uplift. We constrain fault slip-rates since ~18 ka using variations in cosmogenic 36Cl measured on bedrock scarps, mapped using LiDAR and ground penetrating radar, and compare these rates to those inferred from geodesy. The 36Cl data reveal that individual faults typically accumulate meters of displacement relatively rapidly over several thousand years, separated by similar length time intervals when slip-rates are much lower, and activity shifts between faults across strike. Our rates agree with continuum deformation rates when averaged over long spatial or temporal scales (104 yr; 102 km) but over shorter timescales most of the deformation may be accommodated by <30% of the across-strike fault array. We attribute the shifts in activity to temporal variations in the mechanical work of faulting.


5) References 6) Data tables
All analytical data (AMS and sample chemistry) for all samples 1. Field surveying methods

Ground Penetrating Radar
For this study, we operated a Sensors and Software PE-100 GPR system in common offset profiling mode. For the acquisition of high-resolution data we used 200 MHz antenna with a separation distance of 0.5 m and a step size of 0.1 m. According to wave theory, the highest achievable vertical resolution of the survey is one quarter of the operating wavelength (Jol and Bristow, 2003). Using the setup described above we can calculate an ideal vertical resolution of 0.125 m (assuming an average pulse velocity of 0.1 m/ns and dry conditions). However, this value represents the best that can be achieved; in reality the resolution will be slightly less owing to the complexity of ground responses.
Raw radargrams were processed using a common workflow but with varying parameter values due to favourable conditions at all sites. Data processing included time-zero correction, de-wow filtering, bandpass filtering and automatic gain control application to boost signals at depth. To determine average wave velocities a common mid-point (CMP) survey was conducted perpendicular to the survey lines and parallel to the fault planes. The results from these surveys give an approximate average pulse velocity of 0.1 m/ns, comparable to values obtained from the profile data using the shape of diffraction hyperbolae. This value was used to apply the topographic correction to the processed radargrams. Further details on the applications and processing procedures of GPR are presented by Neal (2004), Schrott and Sass (2008) and Jol and Bristow (2003). Data were plotted using Ekko View Deluxe 42 (https://www.sensoft.ca/products/ekko-project/overview/).

Airborne Laser Scanning (ALS)
Airborne laser scanning (ALS) is an active remote sensing technology that acquires 3D coordinates from the ground surface that can be used to generate automated digital terrain and surface models (Ackermann, 1999). Information on range, location and altitude of the target is generated in a 3D domain. The range information is determined by measuring the return time of the laser pulse from the surface to the aircraft mounted sensor (Flood & Gutelius 1997). Similarly, information on location is determined from ground based permanent or campaign GNSS control, and in-flight differential GNSS and inertial navigation systems (Wehr and Lohr, 1999).
For the present study, ALS data were captured using a Leica ALS 50-II system on board a Natural Environment Research Council's Airborne Research and Survey Facility Dornier 228-101 from an altitude of 18,000 ft. The Leica ALS 50-II is a discrete return system operating at a wavelength of 1064nm recording the first, second, third and last return and their intensity. This system has the ability to capture 83,000 pulses per second with a scan frequency of 45Hz, scan angle between 22 and 25º and beam divergence of 0.22mrad. In this study, the point density ranged between 0.9 returns/m2 and 1.2 returns/m2.
The initial stage of creating an elevation model from ALS is reassembling the data into separate return data layers using LAStools™ in ArcGIS 10.2™(http://www.esri.com/). In theory, the last returns represent ground surface, however, in the target area, a significant number of last returns came from forest canopy and low vegetation. The last return data was further filtered to remove vegetation using a Matlab™ Lowest_points script. An elevation model with a grid density of 1m was created using the nearest neighbourhood interpolation method in ArcGIS™ software from the filtered last return data. Subsequently, the elevation model was visualised as shaded relief model with standard ArcGIS™ software settings.

Terrestrial Laser Scanning (TLS)
A Riegl LMS z420i (RIEGL Laser Measurement Systems GmbH, Horn, Austria) was used to collect 38 laser scans during the period 06 April 2008 -10 April, 2008 to characterise the surface topography at each of the 6 sample sites. The scanner has a nominal range of 800m and a specified range error of 5mm at 50m range (0.0001%). Each fault was scanned from an optimum viewing position located (147.8 to 661.1 m) from the fault on the hangingwall side. This gave a good definition of the footwall, scarp and hangingwall and enabled profiles to be created to determine the throw across the fault (e.g., Fig 2(c) in main text plotted using Riscan Pro version 1.2.1 b9 (http://www.riegl.com/products/software-packages/riscan-pro/). Point spacing is determined by the scanner step angle and varies uniformly with range. For the overview scans, point spacing varied from 3 mm closest to the scanner to 0.5 m at maximum range on the footwall with an average of 0.16 m on the fault scarp. Where possible, the faults were also scanned prior to sampling from a position close to each sample site (1-15m) to enable the quality of the fault plane to be recorded prior to the destructive sampling process. Here point spacing varied from 4mm at the base of the scarp to 20mm at the top of the highest scarp ( At some sites (see details of individual sites in Section 3 below), a ruler survey method and clinometer measurements along the sample ladder were combined with the terrestrial LiDAR data to improve the detailed characterisation of the site geometry at both the meter scale and 10-100 m length scale.

Colluvial density estimation
Whole-soil bulk-density of hangingwall colluvium was calculated following the methods of Vincent & Chadwick (1994). Subsequent to trenching, soil horizons were identified and individual bulk-densities were determined for each horizon. An integrated bulk-density was then calculated for the entire sequence and this value was used in the Matlab® modelling. Density values for each site are given in Table S 4.4.1.

Fault plane sampling
For fault plane samples, carbonate bedrock fault samples were collected in continuous sample 'ladders' from a single locality on each of the Pescasseroli, San Sebastiano, Gioia dei Marsi, Parasano, Frattura, Tre Monte and Fiamignano faults. Prior to sampling, each fault was examined along-strike for optimal characterisation, i.e. the exposed bedrock reflected only seismotectonic exhumation rather than modification through a geomorphic process such as erosion, sedimentation or gravitational sliding. Detailed site characteristics are provided in Section 3, with relevant surface and sub-surface survey methods outlined in Section 1.
Once a suitable fault plane site has been selected a trench is dug down to 1-2m and the hangingwall colluvium examined and whole-soil bulk-density calculated. Samples are then taken from the fault scarp by use of a hand-held angle grinder and diamond blade. A continuous ladder of rectangular blocks with dimensions of 15 cm wide, 5 cm high and 2.5 cm deep are cut from the base of the trench to the top of the preserved fault plane, following the down-dip direction of the striae.
Overstepping of sample ladders is utilised to avoid isolated erosive features and fracture infill containing secondary calcite. Once cut, the samples are then numbered, described and    n i i q y n R 1 1 5 photographed before removal by use of a hammer and chisel. A sample numbering protocol is followed whereby the surface-rock interface at the top of the trench is set as the zero referenceframe.
For the upper slope sample at FIAM, the selected site was chosen on an elevated site on the footwall, away from the fault scarp. The sample site had a similar planar dip to the upper slope (c. 33) and was approximately 0.5m above the footwall surface. Multiple 2.5cm deep samples were selected from the surface and analysed separately (see Section 3 and Table 6.1.8 for further data).
Finally, the sample site is characterised in terms of topographic shielding, elevation and geographic coordinates as described in Gosse and Phillips (2001).

36 Cl sampling procedures and analyses
Prior to chemical digestion in the laboratory, the all individual samples are examined in detail for textural and lithological variability such as colour, entrained clasts within the breccia, presence of fine-grained clay cement and secondary calcite and/or iron-magnesium oxides. This information, coupled with separate chemical analyses, is used in conjunction with the 36Cl data to ensure that each sample is providing a robust measure of seismotectonic exhumation.
Preparation and chemical digestion broadly follows the methods outlined by Stone et al. (1996) and more recently in Schlagenhauf et al. (2010). After the removal of secondary calcite and/or metal oxides using a rock-cutting saw and handheld rotary tool, samples are crushed and an aliquot removed for PGNAA analysis. The remaining sample is then sieved to produce a 250-500 µm fraction for chemical processing. Multiple rinsing in 18 MΩ·cm H 2 O and leaching in 0.33N nitric acid ensures removal of contamination by atmospheric chlorine within grain boundaries and defects. A 30 g fraction of leached sample is then used for further preparation. An aliquot preserved for whole rock chemical analysis via ICP-MS. After addition of an isotopically enriched 35/37 Cl carrier, samples are fully digested in 2N nitric acid in an ice bath and any remaining residues are weighed so that precise total dissolved mass can be calculated. Following extraction of an aliquot for corroborative major element determination by ICP-OES, silver nitrate is added to precipitate silver chloride in darkroom conditions. In order to minimise isobaric interference from 36 S during accelerator mass spectrometry (AMS), the silver chloride precipitate is redissolved using an ammonium hydroxide solution and a barium nitrate solution is added. The sample is left for ~48 hours to promote precipitation of barium sulphate crystals. Samples are then passed through a 10µm Anatop syringe filter and the silver chloride re-precipitated. Further purification is achieved using multiple dissolution and reprecipitation steps, followed by water rinses of the precipitate using 18 MΩ·cm H 2 O. After drying, the silver chloride precipitate is pressed into a silver bromide substrate within a copper cathode. Two reagent blanks are processed alongside each batch of 14 samples in order to trace Cl contamination during laboratory analysis. Typical blank 36/35 Cl ratios are on the order of 10 -15 .
Sample 36 Cl and natural-Cl concentrations are measured by acceleratory mass spectrometry (AMS) at the Scottish Universities Environmental Research Centre. AMS is ultrasensitive isotope-ratio mass spectrometry done at high ion energies to resolve molecule and atomic isobaric interference. The sputter ion source is operated for low source memory and gas-stripping is preferred for beam high brightness. 30 MeV 36 Cl 5+ separation from 36 S isobar is wholly by efficient active ion-stopping measurement with a gas ionisation detector rather than involving passive post-stripping. Sample 36 Cl

Site Characterisation methods and individual site characterisation data
It is critical that fault plane sites studied with 36 Cl are fully characterized to show that exhumation from the ground is solely due to fault slip and not due to erosion or sedimentation processes. This supplement documents the geology and geomorphology of the sites in this paper, giving details of the variety of methods used for characterization. We provide two pages of figures for each site including the following datasets. Site parameters are given in Table S  A) Air photos and slope maps derived from LiDAR. These data are needed to show that the fault scarp is continuous along strike for hundreds of metres along a geological fault that offsets mapped stratigraphy (Roberts and Michetti 2004). The fault scarp should be characterized by an exposed fault plane in bedrock limestones that offsets an upper slope and a lower slope that were originally continuous across the fault during the high erosion rate period of the last glacial maximum (LGM; Roberts and Michetti 2004). The images should demonstrate that no alluvial fans or colluvial fans exist that have been fed from incised gullies located above the chosen sample site because their presence would reveal exhumation/burial of the sample site by erosion/sedimentation processes. They should also demonstrate that the hangingwall and footwall cut-offs of the slope formed during the demise of the last glacial maximum are preserved as parallel lines, allowing the slip to be reconstructed back to that time if the slip vector is known (see below). B) Stereographic data recording the fault and slip geometry. The slip vector orientation is needed to demonstrate that the slip can be reconstructed to guide the chosen orientation of sample transects for 36 Cl and the orientation of scarp profiles used to measure the offset of the slopes. The slip vector orientation is defined by the strike and dip of the fault, along with the plunge, and plunge-direction of any frictional wear-striae or corrugations on the fault plane. The preservation of millimetre-scale striae on the fault plane are also used prove minimal erosion of the fault plane after exhumation. C) Structural mapping along the strike of the fault. Strike and dip data, along with plunge and plunge direction data for frictional wear striae and corrugations, must exhibit a consistent pattern along strike to show that the kinematics of the chosen sample site are not anomalous and hence un-representative of the exhumation of the fault scarp. These data 7 are displayed as graphs of these parameters as a function of distance along strike. The data show that the kinematics of the chosen sites are consistent along strike for hundreds of metres and hence the sample sites are representative of the scarps as a whole. D) Scarp profiles. The amount of slip and hence offset of the slopes that formed during the demise of the LGM must be constrained because this is an input parameter for modeling the 36 Cl data. We have used airborne and terrestrial LiDAR supported by geomorphic observation in the field to define these offsets with cm-scale precision. We used the slip vector orientation to define the orientation of scarp profiles. We used geomorphic observations, observations in trenches, and ground penetrating radar (see below) to define the locations of the hangingwall and footwall cut-offs of the slope formed during the demise of the LGM. These cut-offs, alongside the dip of the fault measured from field structural measurements, are used to define the throw, heave and displacement. E) Ground penetrating radar. The position of the hanging wall cut-off must be defined so that slip in the plane of the fault can be measured and used in 36 Cl modelling. We used geomorphic observations to try to identify sites where no hangingwall sedimentation or erosion have occurred to obscure or destroy the hanging wall cut-off (see above). We then checked our interpretations with ground penetrating radar to define the lateral continuity of sub-surface layers of sediment. Confirmation that the hangingwall cut-off is preserved and coincident with the ground surface defined by LiDAR is provided if the sedimentary layers prove to be parallel to the ground surface with no major discontinuities, continuing up to the position of the fault, as is the case with our sites. We check this further by excavating a trench and logging the trench walls to determine sub-surface layering (see below). The ground penetrating radar data also demonstrate the sites have not been affected by massmovement/landsliding (Bubeck et al. 2015). The fault plane dip measured at the surface is plotted as a red arrow on the following figures to indicate the position and dip of the fault. F) Photos of the fault plane and sample locations. We select sampling locations mainly based on the criteria listed above, but sites are prioritised if they also exhibit well-preserved fault planes. We select fault planes where we can prove minimal erosion after exhumation due to their smoothness and the preservation of millimeter-scale frictional wear striae produced at depth by frictional fault slip. We take samples from striated surfaces using a rock-saw, avoiding sites where the fault plane has been degraded by erosion (chemical dissolution, physical plucking -usually along fractures, and biological disturbance due to plants exploiting weaknesses in the fault planes). We also excavate trenches and sample fault plane beneath the ground to allow sample collection and hence measurements of the pre-exposure 36 Cl concentrations. The sub-surface fault planes are in places disturbed by fractures, but we avoid these whilst sampling, selecting surfaces with clear frictional wear striae where possible. G) Trenching excavation. We excavate trenches to (i) check the subsurface layering and hangingwall cut-off location from ground penetrating radar and LiDAR, (ii) expose the fault plane for sub-surface sampling, and (iii) measure the density of the hangingwall material as this shields samples prior to excavation. The sites typically show 10-20 cm of organic rich soil that we assume is Holocene. Deeper material is usually scree with a fine matrix that contains markedly-less organic material -we assume this is colluvium deposited during the last glacial maximum. We correct the location of the hangingwall cut-off using the identified location of the Holocene to LGM sediment transition. We measure the density of the excavated material in trench because this is an input parameter for modeling the 36 Cl.
In summary, we take great care with sample site selection and characterisation because the exposure history and hence slip history is depends on selecting samples with well-defined exhumation solely by fault slip. 1 4. Modelling

Methodology used for modelling fault plane samples
When fitting slip histories to the observed 36 Cl data for each fault site we consider three different approaches depending on the complexity apparent in the data and the presence of independent data from LiDAR, historical earthquake and paleoseismic records: (i) Optimization; (ii) Bayesian Markov Change Monte Carlo (MCMC) with fixed slip-rate change points; (iii) Bayesian MCMC with flexible slip-rate change points that can be iterated. In each approach we use a Matlab script that generates possible slip-histories and which calls another Matlab program, published by Schlagenhauf et al. (2010) for modelling 36 Cl concentrations on bedrock scarps, to quantify the fit of each of the slip histories we consider.  Schlagenhauf et al. (2010) suggest several measures of fit, including weighted Root Mean Squared Error weighted by the uncertainty of the measurements (RMSw), Akaike Information Criterion (AICc) and Chi-Squared. Here we use RMSw and AICc. RMSw is defined as:

Measures of fit to the data
where and are the observed and modelled 36 Cl concentrations, is the significance or error of the measurement and is the number of measurements. Aikake's Information Criterion is defined as where is the number of observations and is the number of parameters used in the model. In the (Schlagenhauf et al., 2010) approach is taken to equal the number slip events used to model the data. This version of AICc is a modification of the original definition of AIC and should be used when the ratio ⁄ is small (i.e. ≤40), which it is in this study.
For the Bayesian parameter estimation approach we need to define the likelihood ( | ), that is the probability of observing the data ( ) given the parameters ( ). A standard approach is to define the likelihood as: where is the standard deviation of the data. Note that this no longer includes the analytical error of the measurements that was included in the .

(i)
Optimization approach (see sites PARA and SSB; Figs. S 4.5.6& S 4.5.7) Optimization can be used in cases where the slip-history is simple and can be explained by few parameters. It involves systematic scanning of the parameter space and in practice works for the estimation of at most 2 parameters. This works well in cases where we assume the slip rate is constant as then we run the model for a number of scarp-ages within a possible range to find the scarp age (SA) that results in the minimum value of the Root Mean Squared weighted by the uncertainty of the measurements (RMSw) or the Akaike Information Criterion (AICc). The latter is strongly related to the RMSw but takes into account the number ofslip events used to model the data. In the current paper slip-histories are modelled as a series of small slip events with constant offsets and inter-event times (slip size estimated from Wells and Coppersmith, 1994;Table S 4.4.1). We also considered offsets (or inter-event times) drawn from a random distribution with a given mean and standard deviation; tests showed that this did not change our conclusions and is not used in the results presented here.
In the case where the slip rate is not constant we need to infer the ages at which changes in the slip rate occurs as well as the total scarp age (SA). Initially, we assume that the height of the change point is known, based on independent observations of fault scarp morphology (e.g. from LIDAR measurements of surface roughness), see section (ii) below. In this case we generate a set of slips as before and count the number of offsets required to generate the sections of the scarp above and below the roughness change point(s). We assume constant slip rate above and below the change point(s) so that the inter-event times are given by the age of the change point and the number of slip events required.

(ii) Bayesian MCMC with fixed change-point heights
In the Bayesian approach we implement a sampling scheme based on the Metropolis Hastings algorithm (Metropolis et al. 1953, Hastings 1970, Sambridge et al. 2006. For a rigorous description of the underlying theory as well as a practical description of the implementation we refer to Sambridge et al. (2006). In brief, the method works as follows: The sampling starts with an arbitrary parameter set for which the fit to the data is calculated using the likelihood function. For each iteration a small random change is proposed to one of the parameters and the resulting fit is either accepted or rejected based on the ratio of the likelihoods. The method also allows for the inclusion of prior knowledge which in practice means that parameters that fit better with the prior distributions have a higher probability of being accepted (prior knowledge refers to information obtained from 3 independent sources, see more below). The theory behind the method prescribes that if we run this for long enough the sampled parameter sets will eventually "converge" and can be interpreted to represent the so called posterior distribution of the parameters, which takes into account the data as well as the prior knowledge about the parameters (Fig. S 4.1.2 is example below for site FIAM).
In our implementation for the 36 Cl data on fault scarps, the offset for each slip event is considered constant. The slip history in the fixed change point approach is fully defined by the elapsed time, ET (= time elapsed since most recent accrued slip), the height and age of the total scarp and the heights and ages of any number of change-points, CP1 etc. In practice we limit the number to those for which we have independent evidence (at most 2 change points) using measures of fault surface roughness from terrestrial LiDAR data. In each iteration we make a small change to the elapsed time, the scarp age or the age of any of the change points. An illustration of the results from this approach for site FIAM is given in Figs . The prior distribution for the scarp age (SA) is a normal distribution with mean of 15000 yrs and standard deviation 2500 yrs (implying that we are 95% certain that the scarp age is between 10000 and 20000 yrs). This is based on the timing of the demise of the LGM in this area (Giraudi & Frezzoti, 1997), associated with a marked drop in hill slope erosion rates and the onset of scarp preservation (Tucker et al., 2011), and also the upper slope age that we obtained at this site (see Section 4.  (Green 1995, Sambridge et al. 2006. In this approach, the number and heights of change-points is not fixed. In each iteration we either add, or remove change-points and/or change the timing of existing ones in addition to changing the elapsed time (ET) and/or scarp age (SA). As before, the slip rate between change points is assumed to be constant and as define in method (ii) above. Although any number of change-points can in theory be considered, the method favors simple solutions with few change-points (Sambridge et al. 2006) and in our results the number of change-points never exceeded 6. Due to the fact that this is a very high dimensional parameter space in which each offset is a potential change point this method takes a long time to converge (cf. trace plots for fixed change points   (2012) only as a source of synthetic data sets for testing our Bayesian approach.
Firstly, the impact of uncertainties that arise from site geometry characterised by LiDAR and structural data gathered in the field (Table S 4.4.1) are compared to uncertainty in scarp age (SA) that we anticipate based on the timing of the demise of the LGM in this area (12 -18 ka). Figure S 4.2.1 shows that the range in SA dominates all of the other uncertainties related to site geometry and this is why we include it as prior information in the Bayesian modelling. Moreover, the range in SA also means that we can only resolve with confidence temporal variations in slip rate above a certain magnitude ( 2 is that where we have inferred SRV = 0.0 for the field data we cannot exclude some temporal variations in slip rate but where we estimate SRV > 0.3 we are confident that the rate variations have indeed been significant (i.e., periods of rapid slip interspersed with periods of relative quiescence). In addition to slip rate variations, a long elapsed time (several hundreds to thousands of years) can be resolved by a change in slope in the 36 Cl profile at the ground surface (  Figure 3 in the main part of the paper reflects variations in average slip rate and that deviations from the simple fan shape reflect temporal variations in slip rate (SRV) and/or long elapsed times (ETs).
We also use synthetic data to test the flexible change point Bayesian approach ( 2(c)) we generated a synthetic 36 Cl data set with the site geometry, sample spacing and sample chemistry at site PARA. These data were then treated in exactly the same way as our field data to demonstrate not only the success of our approach in recovering the main features of the slip history, but also to determine the window length that best characterizes the SRV both for the actual slip history and the slip history inferred from the 36 Cl profile (Fig. S 4.3c&d). We obtain the median posterior estimate of SA = 14.6 ka (+4.0/-3.2 ka at the 90% C.I., which overlaps the actual SA = 15.6 ka); the maximum likelihood estimate SA = 14.4 ka. The maximum likelihood slip history has an SRV = 0.45, compared to 0.44 for the actual slip history, both calculated using a 3000 year sliding window. Furthermore, the SRV for the five highest likelihood fits range from 0.2 to 0.6, i.e, all resolve that SRV > 0 for the synthetic data.

Methodology used for modelling upper slope sample at site FIAM (LGM inheritance)
Using the exposure age calculator of Schimmelfennig (2009)  Cl that accumulated hence why we interpret the measured upper slope 36 Cl concentration as indicative of slope stabilisation and use Schimmelfennig's (2009) age calculator for this sample. Correcting for the finite erosion rate estimated for this site, using Tucker et al.'s formula, we obtain a 'corrected' upper slope stabilisation age of 17.0 +1.7/-1.8 ka (here the uncertainty includes both analytical uncertainty and uncertainty in the erosion rate estimate). This 'corrected age' lies within the 12-18 ka age range for the demise of the LGM in this area (Giraudi and Frezzoti, 1997).
The prior distribution for Scarp Age (SA) used in our Bayesian modelling of the fault plane samples (see Section 4.1) implies that we are 95% certain that the SA is between 10000 and 20000 yrs. It is a normal distribution based on the information on the 12 -18 ka age range for the demise of the LGM from Giraudi and Frezzoti (1997). Any effect of inherited 36 Cl on SA from finite LGM erosion rates ( +/-2 kyrs, see calculation above) thus lies well within the range of our prior and is less than the C.I.'s that we quote on results in this paper. In other words, the effect of inherited 36 Cl on our overall results/conclusions is negligible because we have taken a conservative approach in our implementation of the Bayesian MCMC methodology. We demonstrate this, for example, in  1. Impact on 36 Cl profile shape of uncertainties that arise from site geometry characterised by LiDAR and structural data gathered in the field. All tests (black dashed lines) are based on a single typical field site (similar to PARA) with a 15 ka scarp age. In (b) the variations in scarp height derive from the combined measurement error on upper slope dip angle and fault plane dip. In (d) density range reflects the range of densities measured across all our sites. We use one value for each site but also test that our conclusions on SRV are not sensitive to the specific value. Effect of variations in slip rate (SRV) are considered in Figure S.4.2.2. Differences in scarp age (SA) for the standard site (see grey lines) reflect the age range that defines the demise of Last Glacial Maximum (LGM) = 12-18ka (Giraudi & Frezotti, 1997). The uncertainty in SA clearly dominates over other sources of uncertainty associated with the site geometry. We include prior information on SA in our Bayesian modelling approach (see Section 4.1 of Supplementary Material).
Sensitivity of 36 Cl profile shape to site specific parameters Constant rate slip histories are shown by grey lines and denote the range of rates associated with the 12 -18 k yr time window that defines the demise of the LGM in this area (Giraudi and Frezzotti, 1997  Comparison between theory and field data for 36 Cl concentrations measured at the top of the trench. On faults that have a lower average slip rate (lower SR ave ) the 36 Cl concentration is higher and the rate of increase in 36 Cl concentration with height on the scarp is greater. The real data show a similar correlation but with more scatter, partly due to temporal variations in slip rate (SRV) and long elapsed times (see theoretical cases A, B and D in (c)) but also due to variations in cosmogenic production rate between sites. Grey arrow in (c) indicates the range within which SRV > 0 is difficult to resolve with only trench samples because of the likely range in total scarp age (i.e., 12 -18 ka; Giraudi & Frezzotti, 1997)    The synthetic data set was generated assuming the same chemistry for each sample; the sample spacing and the analytical error bars are the same as that of a real site in this study (PARA). The modelling of the 36 Cl data used the known site geometry (α, β, γ), scarp height, site elevation, and slip size to solve for the timing of each slip event and thus the slip history (red line in (b)). Graphs (c) and (d) show the slip rate variability (SRV) calculated for both the actual and the modelled slip histories respectively. SRV is defined (Cowie et al., 2012) as the standard deviation of the slip rates (SR , SR , etc.) measured over a fixed time window divided by the average slip rate (SR ). SRV is calculated for different window lengths and using two different methods (a sliding window (black dots and lines) versus consecutive time windows (grey dots and lines in (c)). The sliding window method gives more stable results and shows that for window lengths > 2500 years the SRV is much less sensitive to window length. In this study we use a window of length 3000 years, consistent with previous published work in central Italy (Cowie et al., 2012        Grey shade in top row indicates Scarp Age (SA) estimated for each site using our modelling approach. The mean SA = 17.8 ± 4.3 ka (9 independent estimates) represents a regional estimate of the onset of bedrock scarp preservation based on our modelling of the fault plane samples. Our upper slope cosmogenic sample at one site (FIAM) gave a slope stabilization age of 17.0 +1.7/-1.8 ka (corrected for LGM erosion rate, see Section 4.1.1) and lies within 1 standard deviation of this regional estimate. Age information on the demise of the LGM (12 -18 ka) and the associated reduction in erosion rates that led to scarp preservation in this area (Tucker et al., 2011), is included as prior information on scarp age (SA) (blue filled pdf in inset in part (c)). Fits X and Y indicate fits that approximately correspond to the 90% credible intervals (C.I.) on ages for CP1, CP2 and SA. An abrupt change in the shape of the 36Cl profile at the base of the scarp at this site indicates that there is an elapsed time (ET) of several hundred years since the last significant accumulation of slip (see synthetic data Fig. S 4.2.3(b)). The non-zero elapsed time that we use to constrain the modelling for this site is the timing of the 1349 AD earthquake (we set ET = 665 yrs), which historical records strongly suggest ruptured this fault (Guerrieri et al., 2002;Galli and Naso, 2009). Our age estimate for CP1 is 1973 yrs ago (+404/-479 yrs., 90% C.I.), for CP2 is 4516 yrs ago (+2015/-1753 yrs., 90% C.I.) and for SA 13.70 ka (±0.32 ka., 90% C.I.). These age estimates indicate that between 1349 AD and approximately 2000 years ago (i.e., from Roman times to the end of the Middle Ages) displacement accumulation occurred very rapidly, and the slip rate deviated significantly from the Holocene-averaged rate (~1.8 mm/yr) on this fault but that this rapid phase was preceded by a period of lower than average rates of slip, particularly between ~4.5 ka and 13.7 ka (~0.6 mm/yr). SRV = 1.0 is estimated for the highest likelihood slip history at this site (Table S 4.4.2). Including the second change point (CP2) leads to a better fit to the topmost samples, i.e., above height = 18 m. If CP2 is not included then the RMSw increases from 17.3 to 19.32 and the SA we then infer is only 9.5 ka which lies outside the the 12 -18 ka age range of the LGM demise and also inconsistent with the age of the stablisation of the upper slope at this site (17.0 +1.7/-1.8 ka; see Section 4.1.1). In addition we ran models with ± 0.5 m variation in the heights of CP1 and CP2 and there was no significant difference in the results shown here. In Fig. S   Below CP1 the fault plane is consistently low in roughness (mean deviation < a few cm) and preservation is mostly 90-100%. Between CP1 and CP2 roughness changes abruptly (mean deviation of ~10 cm) but shows no systematic increase with height; the degree of fault plane preservation is more variable in this zone (see red and blue areas). Above CP2 roughness and % preservation both show a clear dependence on height consistent with progressive erosion of the scarp top. Areas of the fault plane where blocks have been recently plucked (e.g., bottom left) are easy to distinguish by angular morphology. Our single change point approach allows us to search for the 1st order features of the slip history rather than details that may be artefacts caused by hanging-wall sedimentation/erosion during the Holocene (see Galli et al. 2012). We estimate SA = 11.76 ka (+2.4/-2.6 kyrs 90% C.I) for this site, which is at the lower bound of the prior that we use for SA based on the demise of the LGM in this region (12 -18 ka; Giraudi and Frezzoti, 1997). Our estimate for the age of CP1 is 4914 yrs ago (+1545/-1465 years 90% C.I). Our results suggest only subtle variations in slip rate over the last ~12 ka, no significant elapsed time (ET) and SRV = 0.2 (Table S 4 (a)). The implied SA for a constant rate is 26 ka and, moreover, the fit to the data along the subaerial portion of the 36Cl profile is worse than the maximum likelihood variable rate model shown in (b), i.e., the grey circles overlap some of the analytical error bars but none of the data points. (c) The change point (CP) height inferred using the flexible change point method coincides with a distinct change in scarp morphology revealed by the LiDAR topographic profile through the sample site. The age of CP is ~7.7 ka. SRV for the highest likelihood fit to the data is 1.4. With these data there is no significant ET resolved at this site, i.e. ET < a few hundred years. A relative high slip variability (Table S 4.4.2) is consistent with the highly oblique orientation of this fault (NE-SW) relative to the regional strike (NW-SE) of the overall fault array in this area but SRV is not well constrained by these data. . This age is consistent with the prior that we use for SA based on the demise of the LGM in this region (12 -18 ka; Giraudi and Frezzoti, 1997). These results suggest no significant variations in slip rate have occurred over the last ~15 ka; the average rate of slip over this time interval is 0.54 mm/yr (±0.07 mm/yr 90% C.I.). (c) For comparison a simple optimisation approach returns a similar result, i.e., lowest RMSw = 9.45 for a constant rate of slip and a scarp age of 16.6 ka; thus SRV = 0 as slip rate is constant (Table S 4.4.2). This fault ruptured in the 1915 Fucino earthquake (e.g., Michetti et al., 1996) although the surface offset at this location was less than our sample spacing at this site. The SA implied by the lowest RMSw constant slip rate model is ~20 ka which agrees to within two standard deviations the mean scarp age of 15 ka used as prior information in the modelling of the other sites. The slip rate is 0.2 mm/yr and SRV = 0 as the slip rate is constant for the fit to these data ( Optimisation modelling -Site SSB (Section 3, Table S  We also conducted model runs with a uniform prior on (ET) and found that if ET is less than ≈ 700 ± 300 years then RMSw increases to >20. The prior on the ET was set at 6 ± 2 ka yrs ago. The maximum likelihood fit implies a scarp age (SA) = 22.5 ka. Inherited 36 Cl was included in the modelling at this site to take into account the relatively low erosion rate compared to fault slip rate during the LGM indicated by a less well-developed planar upper slope at this site (Section 3). Using the formula of Tucker et al. (2011) and the measured site geometry we estimated that the inherited 36 Cl is equivalent to 6000yrs of preexposure at this site. The ET for the highest likelihood fit = 5.7 kyrs and a conservative estimate of SRV = 0.9 (using a 3000 year sliding window; Table S 4.4.2). If inherited 36 Cl is not included the estimated SA = 29 ka, which signifinantly older than the demise of the LGM (12-18 ka) and onset of scarp preservation due to reduced Holocene erosion rates (Tucker et al. 2011), but similar results are obtained for ET = 6.7 kyrs and SRV = 1.3. Thus a long ET and high SRV are robust conclusions obtained from our modelling of this site: a constant slip rate model does not fit the data well (RMSw increases from 16 to 24) and implies a scarp age of 33 kyrs. Some erosion of the fault plane during the Holocene is evident at this site (see site charactersation photographs in Section 3) but we sampled from the least degraded areas. The measured profile of 36 Cl concentration increases with height, as at the other sites, whereas even a modest erosion rate (0.02mm/yr) would predict little increase or even a decreasing 36 Cl concentration with height (see blue circle in (a)). Finally, even though the subaerial portion of the 36 Cl profile may be somewhat modified by erosion the subsurface 36 Cl profile at this site indicates that this fault has a very long ET, i.e., several 1000 years, consistent with the lack of evidence for historical ruptures on this fault.     35/37 Cl spike added to sample prior to dissolution. Spike concentration: for samples c1746 -c2038 the mgCl/g solution = 6.1850, 37 at/ 35 at = 0.0510; for samples c2338 -c2500 the mgCl/g solution = 7.1777, 37 at/ 35 at = 0.0510.