The protocataclasite dilemma: in situ ³⁶Cl and REE-Y lessons from an impure limestone fault scarp at Sparta, Greece

Reconstructions of palaeoseismicity are useful for understanding and mitigating seismic hazard risks. We apply cosmogenic ³⁶Cl exposure-age dating and measurements of rare-earth elements and yttrium (REE-Y) concentrations to the palaeoseismic history of the Sparta Fault, Greece. Bayesian-inference Markov chain Monte Carlo (MCMC) modelling of ³⁶Cl concentrations along a 7.2 m long vertical profile on the Sparta Fault scarp at Anogia indicate an increase in the average slip rate of the scarp from 0.8–0.9 mm yr⁻¹ 6.5–7.7 kyr ago to 1.1–1.2 mm yr⁻¹ up to the devastating 464 BCE earthquake. The average exhumation of the entire scarp up to the present day is 0.7–0.8 mm yr⁻¹. Modelling does not indicate additional exhumation of the Sparta Fault after 464 BCE. The Sparta Fault scarp is composed of fault breccia, containing quartz and clay-lined pores, in addition to host-rock-derived clasts of calcite and microcrystalline calcite cement. The impurities control the distribution of REE-Y in the fault scarp surface and contribute spatial variation to ³⁶Cl concentrations, which precludes the identification of individual earthquakes that have exhumed the Sparta Fault scarp from either of these data sets. REE-Y may illustrate processes that localize slip to a discrete fault plane in the Earth's near-surface, but their potential use in palaeoseismicity would benefit from further evaluation.

1 Introduction

Seismic hazard risks are significant in many parts of the world, and studying the magnitude, recurrence, mechanisms, and impacts of past earthquakes helps form a basis for mitigating current and future risk. While historical earthquake records are a crucial archive (Gürpinar, 2005), their spatial distribution is patchy and the recurrence interval of large earthquakes on many faults exceeds historical records. Geologic-based inferences of earthquake timing, recurrence intervals, and the magnitudes of slip and shaking intensity are an essential component of seismic hazard risk mitigation (McCalpin and Nelson, 2009, p. 24). Topographic expressions of tectonic faults, the displacement of surficial sediments revealed in trenches, and geochemical alterations on subaerially exposed fault surfaces may each provide evidence useful to the study of palaeoseismicity (e.g. Benedetti et al., 2002; Dramis and Blumetti, 2005; Michetti et al., 2005; Carcaillet et al., 2008; Manighetti et al., 2010; Mouslopoulou et al., 2011; Cowie et al., 2017a; Mozafari et al., 2022). Here we apply concentrations of cosmic-ray-produced (cosmogenic) ³⁶Cl and rare-earth elements and yttrium (REE-Y) to study palaeoseismicity on the Sparta Fault at Anogia, Greece (Fig. 1a, b).

Figure 1 Study site. (a) The study site location in the Peloponnese, Greece. Key tectonic features are shown. The box indicates the location of panel (b). (b) The location of the Sparta Fault, separating the Taygetos mountains from the Sparta basin. The location of the Anogia field site used both in this study and in Benedetti et al. (2002) is shown. Benedetti et al. (2002) located a second sampling transect at Parori (also shown). The digital elevation model has a 24 m resolution and is derived from ASTER GDEM (GDEM2), which is a product of NASA and METI (Japan). (c) Schematic diagram of the Sparta Fault scarp at Anogia, showing the locations of our vertical ³⁶Cl and REE-Y sampling transects and the ³⁶Cl sampling transect of Benedetti et al. (2002). (d) Photograph showing the location of the Anogia A profile prior to sampling. The existing sample scar is from Benedetti et al. (2002). (e) Photograph showing the location of our REE-Y and drill core profiles after sampling and the Anogia B profile before slab sampling.

The Mediterranean is a densely populated, seismically active region subjected to hundreds of earthquakes of magnitude (M) > 4 every year (Godey et al., 2013; Meng et al., 2021; Ozkula et al., 2023). Within the Aegean tectonic plate (Fig. 1a) and around its margins, there were > 1450 earthquakes during 1998–2010 (Godey et al., 2013), 77 of which were M > 5. In central Greece, earthquakes are associated with normal faults, which occur because of the extension of the Aegean plate (Jolivet et al., 2013). In limestone, they may be identified by spectacular scarps, which form from the accumulation of bedrock slip that occurs during successive earthquakes. Holocene fault scarps can be well-preserved (Armijo et al., 1991), making them suitable targets for palaeoseismic studies.

The concentration of cosmogenic ³⁶Cl has been used to infer palaeoseismic activity in limestone normal faults (Zreda and Noller, 1998; Mitchell et al., 2001; Benedetti et al., 2002; Palumbo et al., 2004; Schlagenhauf et al., 2010; Tesson et al., 2016; Cowie et al., 2017a; Iezzi et al., 2021; Mozafari et al., 2022). This nuclide is produced from spallogenic and muonic reactions occurring in ⁴⁰Ca in limestone exposed at the surface and in the uppermost few metres of the subsurface. Following an earthquake, the newly exposed scarp segment accumulates ³⁶Cl, the concentration of which is dependent upon the duration of subaerial exposure plus the preceding duration of subsurface exposure, potentially allowing the earthquake to be dated. The higher the concentration of ³⁶Cl, the older the earthquake which has exhumed the sampled scarp segment. Because the highest part of a scarp was exhumed first, ages increase towards the top of a scarp.

Despite the potential for ³⁶Cl dating of earthquakes, a given ³⁶Cl concentration profile may not have a unique tectonic reconstruction (Goodall et al., 2021). Scarp-exhuming earthquakes may be temporally clustered, making it challenging to resolve individual earthquakes within the uncertainties of ³⁶Cl measurements (Bubeck et al., 2015). Also, ³⁶Cl concentrations along vertical profiles that deviate from the theoretically predicted patterns appear to be a ubiquitous feature of normal faults developed in limestone (e.g. Benedetti et al., 2002; Palumbo et al., 2004; Tesson et al., 2016; Cowie et al., 2017a; Goodall et al., 2021; Mozafari et al., 2022; Dawood et al., 2024). These challenges have helped motivate the development of probabilistic models for determining exhumation histories from ³⁶Cl concentration profiles (Tikhomirov et al., 2019). For example, Bayesian modelling incorporates prior geologic information (Cowie et al., 2017a; Beck et al., 2018; Tesson and Benedetti, 2019; Goodall et al., 2021, Iezzi et al., 2021) to identify the most probable exhumation history from complex ³⁶Cl data and make inferences on the seismogenic potential of a fault (Tesson et al., 2016).

Measurements of REE-Y have also been used to unravel palaeoseismic information on limestone fault scarps, frequently together with ³⁶Cl dating (Carcaillet et al., 2008; Manighetti et al., 2010; Mouslopoulou et al., 2011; Tesson et al., 2016; Bello et al., 2023; Moraetis et al., 2023). The vertical distribution of REE-Y along fault scarps may result from exchanges with former hanging-wall soil REE-Y before uplift. REE-Y would be leached from subaerially exposed scarp surfaces through calcite dissolution and accumulate in the surfaces of hanging-wall soil, where they form organic complexes (Carcaillet et al., 2008; Bello et al., 2023; Moraetis et al., 2023). Because of low pH, calcite dissolution is highest where the soil surface abuts the scarp, and the REE-Y becomes locally enriched in the adjacent scarp surface through soil-to-scarp REE-Y exchange during reprecipitation of calcite. Peaks in REE-Y on fault scarp surfaces that are now subaerially exposed may therefore represent former soil surfaces, which are now exposed to leaching and subsequent accumulation in the hanging-wall soil, thus completing a cycle. The spacing of these REE-Y peaks may permit identification of the number of slip events and earthquake magnitudes from the vertical displacement lengths. These inferences can be made independently of ³⁶Cl measurements. Using both techniques could provide robust palaeoseismic information for seismic risk assessment models.

The pioneering cosmogenic ³⁶Cl study of the Sparta Fault by Benedetti et al. (2002) motivated our studies. Benedetti et al. (2002) found evidence at Parori (Fig. 1b) for the historically recorded 464 BCE earthquake that destroyed Sparta (Armijo et al., 1991) and five older earthquakes. Interestingly, they were unable to substantiate a displacement from the historical 464 BCE earthquake at nearby Anogia. Our study objectives were to (i) study slip rates on the Sparta Fault at Anogia by taking advantage of recent advances in both the measurement of ³⁶Cl and earthquake modelling, accounting for all ³⁶Cl production pathways and shielding effects (Schlagenhauf et al., 2010), and Bayesian modelling (Goodall et al., 2021) and (ii) complement the ³⁶Cl exposure dating with measurements of REE-Y to best constrain the palaeoseismic history of this fault.

2 Geological setting

The Sparta Fault is a 64 km long, NNW–SSE-striking, range-bounding normal fault in the southern Peloponnese (Fig. 1a, b). It separates the eastern flank of the Taygetos mountains (maximum elevation of 2407 m a.s.l.) from the Sparta basin (Fig. 1b). The Sparta Fault is part of a larger normal fault system, which exceeds 150 km in length, and is matched on the western margin of the Taygetos mountains by the antithetical Kalamata fault and other similar faults located offshore of the Mani Peninsula (Fig. 1a; Armijo et al., 1991; Çal et al., 2024). The subaerially exposed scarp of the Sparta Fault is developed in late Senonian–Eocene limestones of the Ionian unit (Institute for Geology and Subsurface Research, 1969; Armijo et al., 1991). Folded and tilted Permian to early Triassic pelitic and psammitic sedimentary and metasedimentary units outcrop in the Taygetos mountains and are also offset by the Sparta Fault at depth (Institute for Geology and Subsurface Research, 1969; Armijo et al., 1991). Geomorphic evidence for Quaternary uplift along the eastern flank of the Taygetos mountains includes steep triangular facets (20–40°) that are hundreds of metres high along the central portion of the range and decrease in height towards the N and S, wineglass canyons, perched valleys, alluvial fans with up to 4 m of entrenchment near the Sparta Fault trace, and slope-break knickpoints (Armijo et al., 1991; Benedetti et al., 2002; Pope and Wilkinson, 2005; Papanikolaou et al., 2013; Çal et al., 2024). Collectively, the evidence indicates an environment that has been tectonically active during the Holocene.

The Sparta Fault scarp is nearly continuous along strike and it reaches a maximum height of 10–12 m in its central portion but tapers towards both ends. Hanging-wall erosion associated with stream incision can locally form higher scarp segments. The scarp has a 61–64° dip and, in all but a few locations, slickensides have been eroded away following exhumation. The slope of the hanging-wall ground surface matches that upslope of the footwall, which indicates a contiguous hillslope prior to formation of the present scarp and that sediment accumulation at the scarp base is generally minor. Some wedges of sediment are locally present on the hanging wall and in some places are welded to the scarp face, in positions now perched above the hanging wall (Fig. S1). These wedges may have been perched by earthquake-induced displacement on the Sparta Fault or are debris deposits from mass movements that have partly eroded. It is possible that other sediment wedges were also formerly attached to the scarp face but have since fallen off.

Our sampling site is located at Anogia, where the Sparta Fault scarp is 6.8 m high (Fig. 1c), sparsely fractured, and displays a smooth surface texture (Fig. 1d–e). Apparent slickensides are faintly visible as grooves that may have been widened and deepened by weathering, and the surface displays a black coating, like those commonly occurring on limestone and which contain higher concentrations of SiO₂ and Al₂O₃ than the underlying rock (Carcaillet et al., 2008). The scarp surface at Anogia also displays a spatially variable distribution of subaerial weathering features, such as rills and dissolution pits, which we avoided in our sampling. The lower-angle hillslopes on both the footwall (above the fault scarp) and the hanging wall display a patchy distribution of bedrock outcrops and an indurated allochthonous regolith composed of limestone clasts, with a matrix of red aeolian dust and calcite cement. An outcrop of limestone about 50 m upslope of the fault scarp reveals folded and tilted bedding. The bedding nearest to the scarp has a dip of 45–60° and a strike of 268–279°, which corresponds with those for the fault scarp, indicating that faulting appears to exploit these structural weaknesses in the bedrock. We neither observed scarps with a total offset of 2–3 m within tens of metres downslope of the Sparta Fault scarp (Benedetti et al., 2002), nor observed fault scarps within hundreds of metres upslope of the Sparta Fault scarp. If earthquakes, including in 464 BCE, bypassed the fault scarp at Anogia (Benedetti et al., 2002), they did not leave geomorphic expressions that we observed on field reconnaissance.

3 Methods

To study the palaeoseismicity of the Sparta Fault, we combined accelerated mass spectrometry (AMS) measurements of cosmogenic ³⁶Cl concentrations from samples collected from the Sparta Fault scarp with field and laboratory analyses of scarp composition and mineralogy and with field measurements of hanging-wall soil composition and pH. We took these measurements by sampling a vertical ³⁶Cl profile at Anogia, upwards from the ground surface and adjacent to the sampling transect of Benedetti et al. (2002) for direct comparison with the pre-existing record (Fig. 1c). We also took samples for ³⁶Cl and mineralogical analyses, including REE-Y, from a second vertical profile located about 50 m to the south (Fig. 1c). We chose this additional site for its smooth, non-fractured fault scarp surface, and sampling was completed from the top of the scarp to 80 cm below the present ground surface, following hand excavation of a pit.

3.1 ³⁶Cl concentrations

We sampled the first profile (Anogia A), adjacent to the southern margin of the Benedetti et al. (2002) profile, for ³⁶Cl by using an angle grinder to cut $\mathrm{10}\times \mathrm{20}\times \mathrm{3}$ cm ($h\times w\times d)$ slabs from the ground surface to a height of 3.9 m (Fig. 1c, d). Because of a crack in the fault scarp at 1.1 m above the ground, the transect was shifted laterally (towards the north) by 40 cm, thus duplicating the measurement at 1.1 m. A total of 37 samples from this profile were measured for ³⁶Cl concentration. We sampled the second profile (Anogia B), ∼ 50 m further to the south, for ³⁶Cl and mineralogical analyses by initially drilling 14 cores of 4 cm diameter to a depth of 3 cm into the scarp surface (Fig. 1c, e). Four of these cores were spaced at 20 cm intervals below the ground surface, and eight were spaced at 80 cm intervals above the ground surface to a height of 6.4 m, which is 0.4 m below the top of the scarp. These samples were augmented by another two drill core samples at 1.2 and at 6.0 m. Subsequently, we took 20 samples from this profile using an angle grinder to cut $\mathrm{10}\times \mathrm{20}\times \mathrm{3}$ cm ($h\times w\times d)$ slabs from the ground surface to a height of 2.0 m (Fig. 1c). A total of 71 samples from the three profiles were subjected to preparation chemistry for ³⁶Cl targets and measured using AMS.

For ³⁶Cl measurements, limestone samples were crushed to approximately 0.5 mm diameter, and the whole sample was used without removing any size fraction through sieving. Prior to partial dissolution, approximately 120 g of crushed material was washed with deionized water to remove fines. Following Stone et al. (1996), meteoric ³⁶Cl was removed using two cycles of partial dissolution with nitric acid to dissolve 5 % (by mass) of the carbonate each time. To prepare the AMS target, we used 30 g of dried sample spiked with 1 mg of ³⁵Cl-enriched sodium chloride carrier (source: Icon Isotopes, ³⁵Cl 99.635 at. %, ${}^{\mathrm{35}}\mathrm{Cl}{/}^{\mathrm{37}}$Cl = 273) to measure native chloride by isotope dilution. A slurry of the sample and 120 g of deionized water was slowly dissolved with 60 g of concentrated trace-metal-grade nitric acid. Post-dissolution, both liquids and undissolved solids were quantitatively transferred to a centrifuge bottle, where the solids were removed by centrifugation. The supernatant was decanted to another centrifuge bottle, and chloride was precipitated using one molar silver nitrate. After a settling period, the bottle was centrifuged to isolate the silver chloride, which was then washed, dissolved with ammonium hydroxide, and treated with barium nitrate to remove sulfate in preparation of further purification by chromatography. The solution was loaded onto 5 mL of Bio-Rad AG 1-X8 strong anion exchange resin, and chloride was moved through, with 0.50 mmol and then 0.150 mmol nitric acid. After re-precipitation with silver nitrate and a washing step, the silver chloride was dried and packed into silver-bromide-cored copper holders. AMS measurements were performed at the Purdue Rare Isotope Measurement Laboratory according to procedures in Muzikar et al. (2003); standards used for measurement are described in Sharma et al. (1990).

3.1.1 Bayesian modelling of ³⁶Cl concentrations

We use the “SimpleSlips” version of the Bayesian Markov chain Monte Carlo (MCMC) code from Cowie et al. (2017b) to identify slip rates from ³⁶Cl concentrations. SimpleSlips builds upon “modelscarp” from Schlagenhauf et al. (2010), which models the number of earthquakes, their ages, and resulting displacements from ³⁶Cl concentrations based on user-defined inputs. “Modelscarp” accounts for each ³⁶Cl production pathway in limestone (Table 1 in Schlagenhauf et al., 2010). SimpleSlips is available through https://github.com/lcgregory/SimpleSlips (last access: 5 November 2024) and we follow the procedures of Cowie et al. (2017a, b) and Goodall et al. (2021) to generate potential slip histories.

Table 1 Parameters for MCMC modelling of slip rate.

α is hanging-wall colluvial surface dip angle. β is scarp dip angle. γ is the dip angle of the hillslope above the fault scarp. ε is scarp erosion rate. Pre is pre-exposure duration. Scarp age is the initial estimate of exhumation of the oldest (highest) part of the scarp. Elapsed time is the estimated duration following the last earthquake. The ³⁶Cl production rate of 59.4 ± 4.3 atoms g⁻¹ yr⁻¹ is taken from Schlagenhauf et al. (2010), re-calculated from Lifton et al. (2005). When using the ³⁶Cl production rate of 48.8 ± 3.5 atoms g⁻¹ yr⁻¹ from Stone et al. (1996), Pre is 10.6 kyr; otherwise, all other parameters are fixed.

SimpleSlips models (i) scarp age, which is the timing of the earthquake to exhume the uppermost, and therefore oldest, part of the fault scarp; (ii) the time at which each subsequent earthquake occurred and the corresponding height of the exhumed scarp; and (iii) the time since the most recent earthquake exhumed the lowest part of the fault scarp (elapsed time). The exhumation of the entire scarp is attributed to a user-defined number of earthquakes that each exhumed the same vertical length of scarp. This constant length stipulation is a requirement of SimpleSlips and is acceptable because the MCMC code focuses on modelling slip rates rather than identifying individual earthquakes, which is an appropriate methodology for our ³⁶Cl data. The timing of each earthquake, apart from the first and last, is therefore dependent on the selected number of earthquakes. We follow Goodall et al. (2021) in using the flexible change point method of Cowie et al. (2017b), which allows variable slip rates between iterations.

We parametrized SimpleSlips as follows: (i) we defined the scarp age as a 1σ normally distributed prior of 8000 ± 1500 years. This selection is partly based on the record that contemporary fault scarps in the Mediterranean region are Holocene in age. They have been exhumed since the Last Glacial Maximum (LGM) because hillslope bedrock erosion and regolith transport rates were much higher during the LGM (e.g. more than 10 times higher for the Magnola fault in Italy; Tucker et al., 2011), preventing ruptured fault scarps from persisting as subaerially exposed features (e.g. Benedetti et al., 2002; Cowie et al., 2017a; Goodall et al., 2021). The adopted scarp age is refined through fitting a modelled ³⁶Cl concentration profile to the measured ³⁶Cl concentration profile using the “modelscarp” code of SimpleSlips. The scarp age is also balanced by the period of pre-exposure, which is the ³⁶Cl inventory that accumulated in the bedrock while it was mantled by up to a few metres of colluvium before initial post-LGM subaerial exposure. A wide Gaussian prior (5000–16 000 years) is assigned in our modelling to account for the uncertainty in scarp age. (ii) The elapsed time is defined as 2500 years, based on the youngest known earthquake on the Sparta Fault of 464 BCE. We assign a 1σ uncertainty of 1000 years to reflect uncertainty in the historical record. (iii) To further define the most likely slip rate history for the Sparta Fault, we completed multiple model runs with varying numbers of earthquakes (three to six) and ³⁶Cl spallation production rates (48.8 ± 3.5 to 59.4 ± 4.3 atoms g Ca⁻¹ yr⁻¹). These end member production rates are respectively from Stone et al. (1996) and a calculation using Lifton et al. (2005) by Schlagenhauf et al. (2010). All model runs used the temporally variable geomagnetic field of Lifton et al. (2005) to scale the ³⁶Cl spallation production rate, and spallation production rates for K, Ti, and Fe are as shown in Table 1 in Schlagenhauf et al. (2010). Scarp age and elapsed time are the priors in the MCMC model, the number of earthquakes defines the timing and location on the scarp of slip change points, and prior probabilities are as defined in SimpleSlips.

The SimpleSlips MCMC algorithm generates a slip history, using the input parameters conditioned on prior probability, to construct a forward model of ³⁶Cl concentrations for this slip history (Goodall et al., 2021). The quality of the slip history solution is then assessed by comparing modelled and measured ³⁶Cl concentration profiles. The algorithm iteratively adjusts a parameter defining the slip history and recalculates a new forward model solution. Acceptance of the new slip history hinges either on its likelihood surpassing that of the prior model or on the ratio of new to current likelihood exceeding a randomly selected value from a uniform distribution between zero and one. Otherwise, the new model solution is discarded, adhering to the principles of the Metropolis–Hastings algorithm (Metropolis et al., 1953; Hastings, 1970). We ran this process for 200 000 iterations, using the parameters in Table 1, and results were either assessed on 160 000 iterations, following exclusion of a burn-in phase of 40 000 iterations to mitigate the influence of initial parameters, or on the 10 000 most probable models.

3.2 Sparta Fault scarp composition

Fault scarp chemical composition and mineralogy were analysed from Anogia B as follows. An initial elemental analysis was done in the field on the Sparta Fault scarp surface using an Olympus Innov-X Delta (40 kV) handheld X-ray fluorescence (XRF) device. This instrument performs elemental analyses with a circular sample spot of 8 mm diameter and can measure elements heavier than Na. All elements lighter than Mg are reported as lighter elements (LEs). Of the elements that compose REE-Y, the instrument was only capable of measuring yttrium. Sampling was done at an interval of 5 cm (or less) over a 7.7 m vertical profile, beginning ∼ 0.9 m below the hanging-wall soil surface. This profile corresponds with the location of the drill core and 2.0 m long ³⁶Cl profiles at Anogia B but was measured before either drilling or slab sampling (Fig. 1c).

For more detailed analyses of elements, including REE-Y, a total of 39 cores (22 and 35 mm diameters and to depths of ∼ 3 cm) were collected every 20 cm from a vertical transect at Anogia B using a portable drill (Fig. 1d). The outermost 1 mm was removed from each core prior to crushing to avoid contamination from the black surface coating. The next 15 mm of each core was then rinsed with cold water, air-dried, and crushed using a grinder with a steel mortar to a grain size of < 100 µm. This crushing technique might supply additional REE-Y to samples (Hickson and Juras, 1986), but, if so, this likely occurs systematically across samples, and we are more interested in spatial trends, which we confirm independently using the handheld XRF device, than absolute abundances. The crushed samples were then analysed for major and trace elements using fusion inductively coupled plasma mass spectrometry (FUS-ICP-MS) at Activation Laboratories (Ontario, Canada).

We complemented the FUS-ICP-MS analyses with spot elemental analysis of one rock core from 1.1 m above the scarp base at Anogia B to make a high-resolution determination of any spatial variations in the scarp composition. This was done with an energy-dispersive X-ray spectroscope (EDS) attached to an environmental scanning electron microscope (ESEM). We used a Quanta 650 FEG with an Oxford INCA EDS, and the analysis took place in a high-vacuum environment at 20 kV. The technique is incapable of detecting REE-Y because their concentrations are too low. Photomicrographs and backscatter images of pore spaces were also taken using the ESEM. These analyses were completed at the Department of Geological Sciences, Stockholm University.

3.3 Sparta Fault scarp mineralogy

A modal analysis of mineral fractions was completed on thin sections taken from the remaining 38 core samples. This was done by counting 1000 points on each thin section (Hutchison, 1974) using a Pelcon automatic point counter attached to a Leica (DM LSP) optical microscope. This point counter comprises a stepping frame attached to a control box (power supply) and is also connected to a computer for statistical analyses using Pelcon software version 2. The point counting and mineral identification was made using an objective working distance of 1.52 mm. The line section pre-set step-length was 0.3 mm, and the line section distance was 1.5–2 mm. The point counting permitted a detailed quantitative analysis of the mineralogy of the Sparta Fault scarp surface. This detailed mineralogy was then compared with the chemical composition data to determine whether phases other than the host limestone were present.

3.4 Hanging-wall soil chemistry and pH

Soil chemistry and pH were measured in samples taken at ∼ 10 cm intervals to a depth of ∼ 90 cm in the pit excavated at the base of the Anogia B profile (Fig. 1c). The elemental analysis was again done with the handheld XRF device. Indicator strips were used to measure pH from mixtures of a 1:1 mass ratio of soil:distilled water and soil:1 M KCl (Sikora and Moore, 2014). These analyses help determine the vertical distribution of REE-Y in the soil (using yttrium as a proxy) and indicate how they might correlate with pH and the vertical distribution of REE-Y in the fault scarp below the soil surface.

4 Results

4.1 Sparta Fault ³⁶Cl concentrations

The cosmogenic nuclide ³⁶Cl concentrations from our three profiles (Table S1) and the original Benedetti et al. (2002) ³⁶Cl concentrations are compared in Fig. 2. The Anogia A and Anogia B profiles display corresponding trends of increasing ³⁶Cl concentrations with increasing height on the fault scarp. Only at 1.6 m do the trends strongly deviate from each other. The Anogia B profile indicates generally lower ³⁶Cl concentrations, including 6 of 19 points that do not overlap within uncertainty with data points at corresponding elevations on the Anogia A profile. Four of those points are located from 1.0 to 1.3 m. In comparison with Anogia A, the adjacent segment of the Benedetti et al. (2002) profile (0–4 m) shows ³⁶Cl concentrations that are on average 19 % higher. Uncertainties (1σ) for data points comprising each profile are almost identical, displaying a mean of 3.8 % for the Benedetti et al. (2002) profile versus 3.9 % for the Anogia A and Anogia B profiles. However, the Benedetti et al. (2002) profile displays more variation between adjacent sample points than is evident in our profiles. While concentrations differ between the three longest profiles, they show a consistent gradient up to ∼ 4 m on the scarp. Above 4 m on the scarp, both our Anogia B drill core profile and the Benedetti et al. (2002) profile display matching lower gradients. While differences in measured concentrations between our two profiles and the Benedetti et al. (2002) profile might be expected given technical advances between measurements, the Anogia A and B profiles do not mirror each other, despite them being horizontally separated by only ∼ 50 m. This inability to consistently replicate measurements along the two profiles justifies a focus on identifying slip rates using the Goodall et al. (2021) model, rather than individual earthquakes, also because up-scarp ³⁶Cl concentration gradients are more consistent between the profiles.

Figure 2 Sparta Fault ³⁶Cl concentration profiles. Error bars indicate 1σ measurement uncertainties.

Slip rates for the Sparta Fault are explored through comparing scarp exhumation generated by three, five, and six modelled earthquakes, where each earthquake exhumes 183, 122, and 104 cm, respectively (calculated by dividing the scarp height by the number of model earthquakes). We focus our analyses on the Anogia A profile supplemented with drill core samples from above 3.9 m on the scarp and from the scarp surface buried by colluvium. This combined profile was chosen for modelling both because the Anogia A profile was sampled at 10 cm intervals up to 3.9 m on the scarp, versus only 2.1 m for Anogia B, and because Anogia A is located adjacent to the Benedetti et al. (2002) profile. Furthermore, MCMC modelling of ³⁶Cl concentrations did not converge with measured concentrations for the full Anogia B profile (i.e. including the drill core samples above 2.1 m) but rather only for the intensively sampled lowermost 2.1 m plus subsurface drill core samples. Modelling only the lowermost 2.1 m plus subsurface drill core samples necessitated changes to scarp age and pre-exposure from those used for the Anogia A plus drill core sample profile because this lowermost part of the scarp has a younger age and because the 2.1 m profile length is indivisible into the 6.5 m length of the Anogia A plus drill core sample profile. These changes, especially to scarp age, invalidate comparisons of slip rates between the two profiles. We did not measure compositions for the Anogia A samples, so we use a mean scarp composition from Anogia B in our modelling. Results from the SimpleSlips model (Cowie et al., 2017b; Goodall et al., 2021) applied to the Anogia A plus the drill core profile are shown below and in Fig. S2, respectively, for end-member ³⁶Cl productions rates of 59.4 ± 4.3 and 48.8 ± 3.5 atoms g Ca⁻¹ yr⁻¹. Geochemical data for the fault scarp used in modelling are shown in Table S2. Modelling results from Anogia B (lowermost 2.1 m and subsurface drill core samples and the entire profile) are shown in Fig. S3.

The results of the Bayesian inference MCMC modelling of ³⁶Cl data from the Sparta Fault are shown in Figs. 3–5. The accepted scarp exhumation models (n=160 000) are shown in slip versus time histograms (Fig. 3a). The maximum a posteriori probability (MAP) model, shown by the red line, deviates slightly from the maximum model density (mean model, black line) for each slip segment but more so for the slip segment at 4.9–6.1 m on the scarp. It indicates three exhumation events between 2.4 and 6.1 m on the scarp, which are closely spaced in time at 5000–6000 years ago. The 95 % confidence intervals (Fig. 3b) illustrate little change in variance between model results from lower, younger parts of the scarp to older, higher parts of the scarp, although the MAP model deviates towards being younger than the mean model towards the top of the scarp. The range of accepted models fits the measured ³⁶Cl data well (Fig. 3c) but accommodates a broad range of corresponding slip histories along the entire vertical length of the scarp (Fig. 3d).

Figure 3 Markov chain Monte Carlo (MCMC) model fits to measured ³⁶Cl concentrations and model slip histories, Anogia A + drill core profile. Slip accumulation is shown for five model earthquakes that each exhume the same vertical length of scarp rather than reflecting the magnitude and timing of historical earthquakes. The red line in panels (a), (b), and (d) is the maximum a posteriori probability (MAP) estimation model, which is the maximum likelihood multiplied by the prior probability based on scarp age. Each panel includes 160 000 iterations, following removal of a burn-in of the first 40 000 iterations. (a) Histogram showing the distribution of accepted model slip histories in slip space versus time. The density of overlapping models increases from warm to cool colours. The mean model and 95 % confidence bounds are also shown. (b) The 95 % confidence bounds of the smoothed model distribution (black lines) calculated for age at each step in the slip. The mean (black line) and MAP (red line) slip histories are also plotted. (c) Model fits to measured ³⁶Cl concentrations (circles). The coloured lines represent a selection of 160 model fits from low (yellow) to high (blue) probability at equal intervals (1000) through the distribution. The black lines indicate 1σ measurement uncertainties. (d) Slip histories through five model earthquakes corresponding to MCMC fits shown in panel (c). Results are shown for a ³⁶Cl production rate of 59.4 ± 4.3 atoms g Ca⁻¹ yr⁻¹. Refer to Fig. S2 for equivalent results using a production rate of 48.8 ± 3.5 atoms g Ca⁻¹ yr⁻¹.

Statistics for how well the MCMC modelling fits the measured ³⁶Cl data and our initial estimates of scarp age (8000 years) and elapsed time (2500 years) are illustrated in Fig. 4 and summarized in Table S2. The posterior probability distribution function indicates that the elapsed time since the most recent earthquake is consistent with the 464 BCE earthquake (mean of 2501 ± 173 years; Fig. 4a, Table S2). In contrast, the time when the scarp started to form (scarp age), presumably through a decrease in hillslope erosion following the LGM, is indicted by the posterior probability distribution to have been longer than our initial estimate of 8000 years (mean of 8742 ± 502 years; Fig. 4b, Table  S2). Mean values of likelihood, weighted root-mean-square (RMS_w), and corrected Akaike information criterion (AICc) are 0.25–0.28, 13.9–14.6, and 863–893, respectively, across the range of the modelled number of slip events (Fig. 4d and e, Table S2), indicating that the number of earthquakes (change points) has a minor influence on modelling a fit to measured ³⁶Cl concentrations.

Figure 4 Statistical plots for Markov chain Monte Carlo (MCMC) modelling. Results are shown for three, five, and six model earthquakes. Solid vertical lines indicate the median of each distribution, whereas dashed vertical lines indicate 95 % confidence intervals, colour-coded according to the number of modelled earthquakes. Posterior probability distribution functions from all models for (a) elapsed time and (b) scarp age. Distributions of (c) likelihood, (d) weighted root-mean-square (RMS_w), and (e) corrected Akaike information criterion (AICc) of slip history calculated for modelled ³⁶Cl concentrations compared to the measured values. The most probable (top 6 %) of models are shown in panels (a) and (b), whereas panels (c), (d), and (e) show statistics for 160 000 models following removal of a burn-in of the first 40 000 iterations.

The slip rate for the Sparta Fault is calculated from the most probable of models (i.e. the top 6.25 % of fits to the ³⁶Cl data (n= 10 000); Fig. 5, Table 2). For the entire vertical length of the fault scarp and five modelled earthquakes, both the mean and MAP slip rates are 0.7–0.8 mm yr⁻¹ for end-member ³⁶Cl production rates, calculated up to the present day (Fig. 5a). For the same calculation but excluding the 2500 years since the most recent known earthquake at 464 BCE, the slip rates are higher, with mean and MAP values of 1.1 and 1.2 mm yr⁻¹, respectively (Fig. 5b). The lowest 3.7 m of the fault scarp is the most recently exhumed scarp segment and the most intensively sampled. It displays a steep ³⁶Cl concentration gradient, which indicates matching mean and MAP slip rates of 1.0 mm yr⁻¹, for five model earthquakes (Fig. 5c). The highest 2.5 m of the scarp displays a gentler ³⁶Cl concentration gradient relative to the bottom 3.9 m of the scarp as indicated by our drill core samples and the Benedetti et al. (2002) profile. The mean and MAP slip rates for this scarp segment are therefore lower, at 0.8–0.9 mm yr⁻¹ (Fig. 5d). Varying the number of earthquakes between three and six has a minor influence on the calculated slip rates (Table 2). An increase in mean slip rate occurred between 6.7 and 5.3 kyr (Fig. 5e).

Figure 5 Slip rates for the Sparta Fault at Anogia (Anogia A plus drill core profile) from Markov chain Monte Carlo (MCMC) modelling. Results are shown for three, five, and six model earthquakes. In each panel, the most probable (top 6 %) models calculated from the median scarp age and scarp height are shown. Solid and dashed vertical lines indicate the mean and maximum a posteriori probability (MAP) estimation for each distribution, respectively. Slip rates are shown for three, five, and six model earthquakes, using a ³⁶Cl production rate of 59.4 ± 4.3 atoms g Ca⁻¹ yr⁻¹. (a) The distribution of the most probable slip rate for the entire scarp calculated up to the present day. (b) The distribution of the most probable slip rate for the entire scarp calculated up to the last known earthquake at 464 BCE. (c) The distribution of the most probable slip rate for the lower segment of the scarp. (d) The distribution of the most probable slip rate for the uppermost segment of the fault scarp. (e) Mean slip rate over time. Slip rates using a ³⁶Cl production rate of 48.8 ± 3.5 atoms g Ca⁻¹ yr⁻¹ are shown in Fig. S2.

Table 2 Slip rates for the Sparta Fault at Anogia from the best Markov chain Monte Carlo models (n= 10 000) for end-member ³⁶Cl production rates and a varying number of model earthquakes.

MAP is maximum a posteriori probability

4.2 Granulometry of the Sparta Fault scarp surface

A first look at the Sparta Fault scarp surface yields a misleading impression of homogeneous limestone (Figs. 1, 6a), whereas close inspection of the core samples instead reveals a typical fault breccia (Fig. 6b–d). This breccia consists of angular to rounded limestone clasts with axes of 1–7 mm (in the two-dimensional view provided by thin sections) surrounded by matrix/cement in which clasts are < 0.1 mm in length. The fault breccia is defined as a protocataclasite, according to the classification of Woodcock and Mort (2008). The composition of the protocataclasite displays large spatial variations, with some portions containing abundant clasts (Fig. 6c), whereas others are dominated by fine matrix (Fig. 6d). The proportion of clasts > 2 mm ranges from 5 % to 20 % vertically along the fault scarp, and the proportion of matrix ranges from 5 % to 60 %. We did not measure the thickness of the protocataclasite, but it exceeds the 3 cm depth of our drill cores.

Figure 6 The heterogeneous fault breccia that comprises the Sparta Fault scarp surface. (a) The Sparta Fault scarp surface appears smooth and homogenous, as illustrated by this photograph of the scarp base at Anogia (the upper half of the dug trench in the foreground). (b) Fault breccia is revealed in a cut drill core, where clasts of host limestone are cemented in a fine matrix. (c) A photomicrograph shows limestone clasts (dotted outlines) comprising about 60 % of the thin-section area. (d) A photomicrograph shows fine matrix comprising about 60 % of the thin-section area. Arrows indicate quartz.

4.3 Sparta Fault scarp composition and mineralogy

In addition to a spatially variable granulometry, the fault scarp shows a spatially variable distribution of major and trace elements. The major component is, as expected for limestone, CaO (mean 52.22 %), but its concentration varies between 43.83 % and 56.64 % (Table S3), which exceeds spatial variations in CaO seen elsewhere in limestone normal fault scarps (Carcaillet et al., 2008; Tesson et al., 2016). Quartz (SiO₂) also occurs, and it also displays spatial variations (0.10 %–20.82 %), with broad peaks occurring at 0.5–0.4 m below the ground and 0.9–1.2, 4.6–4.8, and 6.0–6.2 m along the vertical fault scarp profile (Fig. 7, Table S3). An additional peak in SiO₂, but which is not seen in point counting of quartz, occurs at 6.6 m (Fig. 7, Tables S3 and S4). The point counting and geochemical measurements were done on different aliquots of the sample slab extracted from each level, thereby causing variation in quartz concentration between the two methods, which is prominent at 6.6 m. Mean concentrations of other major elements are low in bulk samples, including Al₂O₃ (0.21 %), MgO (0.16 %), Fe₂O₃ (0.09 %), P₂O₅ (0.07 %), and K₂O (0.05 %; Table S3). However, EDS measurements, such as those shown in Fig. 8a, reveal that the concentrations of some elements are frequently much higher in intergranular pores (Fig. 8c) than elsewhere in the fault scarp, including Si ≤ 38.3 %, Al ≤ 11.7 %, Fe ≤ 48.4 %, and K ≤ 7.1 % (Table S5). Furthermore, intergranular pores and quartz frequently occur together (Fig. 8b), and the concentration of Al₂O₃ covaries with the much more abundant quartz (SiO₂) (Fig. 7).

Figure 7 Concentrations of Al₂O₃ and SiO₂, and quartz abundances from point counting, along a vertical profile, Sparta Fault scarp, Anogia. The concentration of each element (C_i) is normalized to its mean concentration through the profile ($C_i/\stackrel{\mathrm{‾}}{C}$). The ³⁶Cl concentration profile is shown for reference.

Figure 8 Energy-dispersive X-ray spectroscope (EDS) elemental abundances and environmental scanning electron microscope (ESEM) and backscatter SEM imagery of a thin section of fault breccia comprising the Sparta Fault scarp surface at 1.1 m above the hanging wall. (a) Element abundances in a pore, the location of which is shown in panels (b) and (c). Si, Al, and K are abundant relative to Ca, which indicates that clay, e.g. illite, is lining the pore. (b) Quartz is an abundant constituent of the thin-section matrix. (c) Porosity, shown in black; note its spatial association with quartz. The location of the sample used in panel (a) is in a small pore, shown in the inset.

Quartz is revealed by microscopy to be present as randomly oriented rounded to angular grains that are < 50 µm in diameter (Figs. 6d, 8b). Quartz is a constituent of the protocataclasite fine matrix that is mostly comprised of microcrystalline calcite precipitates and which cements larger host-rock-derived CaCO₃ clasts (Figs. 6b–d, 8b, 9a). Point counting further reveals quartz modes ranging from 0.1 % to 15.4 % of the thin-section area (Table S4), with higher abundances correlating to higher abundances of fine matrix. The spatial correlations between SiO₂, quartz abundances on point counting, and fine matrix are further strengthened by EDS spot elemental analyses (Fig. 9). Here, the two selected spots in the fine matrix display Si abundances of 29.7 % and 28.9 %, which contrasts with 1.7 % and 0.9 % for the two spots located on clasts. CaO abundances display an inverse relationship with SiO₂ (33.7 % and 31.2 % for the clasts versus 4.8 % and 5.1 % for the fine matrix). SiO₂ is present largely as quartz, as evidenced by the strong spatial correlation between quartz and SiO₂ along the vertical profile (Fig. 7). Quartz can therefore be used as a proxy for fine-matrix abundances in the Sparta Fault scarp.

Figure 9 Concentrations of Si in the Sparta Fault breccia, 1.1 m above the scarp base at Anogia. (a) A cut drill core from the Sparta Fault scarp at Anogia showing limestone clasts cemented in fine matrix. The circled fine matrix is examined under high resolution in panels (b) to (d). (b) An ESEM image showing the sample location for spot elemental analysis (rectangle). (c) Sample points for elemental analysis using EDS, with values shown in the table. (d) The abundance of Si in the fine matrix illustrated in magenta for the circled part of the thin section shown in panel (a).

In addition to the spatial relationship between quartz and fine matrix, we observed in backscatter SEM images that pore spaces, which frequently harbour higher concentrations of Si, Al, K, and/or Fe than host-rock-derived clasts, are also more abundant in the fine matrix (Fig. 8c). These observations provide evidence that clay particles (< µm-scale) frequently coat pore spaces. The abundance of quartz therefore also provides a proxy for the abundance of clay-coated pore spaces.

Concentrations of REE-Y vary in a wave-like pattern along the vertical profile, with maxima occurring at −0.4, 0.8, 2.6, 4.0, and 6.4 m (Y = 1.2–11.1 ppm; Table S6, Fig. 10). These maxima do not systematically decrease with vertical distance above the hanging wall and are not highest in the soil-mantled portion of the scarp. Yttrium (mean 6.3 ppm), La (mean 5.04 ppm), Nd (mean 3.54 ppm), and Ce (mean 2.31 ppm) have the highest concentrations, whereas all other REE-Y are < 1 ppm (Table S6). The concentrations of REE-Y elements co-vary vertically along the scarp surface (R²= 0.95; Fig. 10a).

Figure 10 Vertical distribution of REE-Y elements on the Sparta Fault (Anogia B profile). (a) REE-Y concentrations. Mean values for all REE-Y elements at each sample point are shown in red dots, whereas individual REE-Y elements are shown in blue dots. The concentration of each element (C_i) is normalized to its mean concentration through the profile ($C_i/\stackrel{\mathrm{‾}}{C}$). Concentrations of Al₂O₃ are shown for reference. (b) LREE:HREE ratio. There are two measurements at −0.6 m.

There is no depletion of light (LREE) relative to heavy (HREE) rare-earth elements with increasing height on the subaerially exposed fault scarp, where it ranges between 3.9 and 5.1 (Figs. 10b, 11a, Table S6). However, there is a relative depletion of LREE on the scarp surface buried by soil ($\mathrm{LREE}/\mathrm{HREE}$ is 3.2 to 4.0; Figs. 10b, 11a), with the least depletion at 0.40 m depth and progressively larger LREE depletion with increasing depth. Peaks and troughs in the $\mathrm{LREE}/\mathrm{HREE}$ ratio along the vertical profile poorly match peaks and troughs in REE-Y concentrations (Fig. 10a, b), although local minima correspond at 3.0 and 5.2 m on the scarp. Accordingly, the correlation between $\mathrm{LREE}/\mathrm{HREE}$ and total REE-Y concentration is only weak (R²= 0.36; Fig. 11b).

Figure 11 REE-Y elements on the Sparta Fault (Anogia B profile). (a) Concentrations of rare-earth elements and yttrium (REE-Y) normalized to chondrite composition (McDonough and Sun, 1995). Each line shows a measured location on the scarp surface. The two low REE-Y outliers at 5.2 and 5.4 m also have exceptionally low SiO₂ and Al₂O₃. (b) LREE:HREE versus REE-Y sum. The R² value is for a linear fit. (c) LREE, HREE, and LREE:HREE versus Al₂O₃ (wt %). The R² values are for linear fits. In each panel, the six subsurface samples (≤ 0 m), including a replicate measurement at −0.6 m, are shown as filled circles.

REE-Y concentration maxima occur at locations that correspond closely with the Al₂O₃ maxima (Fig. 10a, Table S6). Accordingly, LREE, HREE, and total REE-Y are strongly correlated with Al₂O₃ (R²= 0.92; Figs. 11c, S3a). Spatial correlations between REE-Y and SiO₂ and K₂O are also observed (R²= 0.56 and 0.87, respectively; Fig. S4c, e). While REE-Y concentrations vary in wave-like pattern along the scarp, REE-Y is not enriched, and LREE is depleted relative to HREE, in the soil-covered scarp surface.

4.4 Hanging-wall soil chemistry and pH

The terra rosa soil mantling the hanging wall primarily comprises aeolian dust (Muhs et al., 2010) and carbonate clasts. At our sample site, the soil thickness at the base of the Sparta Fault scarp is 0.8 m, and this appears to be stable, at least over the timescale of scarp surface dissolution, as evidenced by a much smoother scarp surface texture below the soil surface compared with the subaerially exposed scarp. Below the organic horizon (∼ 0.1 m thick) the soil is welded, probably by calcite precipitates, and horizons are absent. Soil pH is, in general, slightly acidic along the excavated vertical profile, remaining within a 6.2 to 7.0 range (Fig. 12a, Table S7). An outlier occurs at −0.30 m, where the pH is 5.6 ± 0.2. Soil composition varies with depth (Fig. 12b, Table S8). Concentrations of Si, Al, and K are lower in the organic horizon (11 %, 0 %–5 %, and 0.4 %, respectively) compared with the remainder of the profile (18 %–30 %, 5 %–10 %, and 0.5 %–0.9 %, respectively), whereas the concentrations of LE, which includes C, are, as expected, higher in the organic horizon (75 %–80 %) than in the lower profile segment (51 %–64 %). The concentration of yttrium ranges from a maximum of 36–39 ppm at 0.5–0.6 m depth to a minimum of 11 ppm at 0.1 m depth, and its vertical distribution correlates positively with Si (R²= 0.71), Al (R²= 0.45), and K (R²= 0.54) and negatively with pH (R²= −0.52; Figs. 12c, S4b, d, f, Table S8).

Figure 12 Hanging-wall soil chemistry, adjacent to the Sparta Fault scarp at Anogia. (a) Soil pH along a vertical profile measured from soil mixed with distilled H₂O and 1M KCl. Uncertainty ranges show the ≤ 0.5 resolution of the indicator strips. (b) Concentrations of Si, Y, and elements too light to be measured using handheld XRF (LE, including C) along the vertical soil profile. Each element has been normalized through division by its mean concentration through the soil. (c) Y concentrations plotted against pH (measured from 1:1 KCl) and Si concentration at each measured depth interval beneath the soil surface.

5 Discussion

5.1 Slip rate on the Sparta Fault at Anogia

The average exhumation rate of the entire scarp up to the present day is 0.7–0.8 mm yr⁻¹ (Fig. 5a, Table 2). This compares with an exhumation rate of 1.1–1.2 mm yr⁻¹ up to the 464 BCE earthquake (if an earthquake occurred now, the rate up to the present day would increase). These values compare favourably with estimates of 0.8 and 1 mm yr⁻¹ from Papanikolou et al. (2013) and Benedetti et al. (2002), respectively, but are higher than the 0.5 mm yr⁻¹ exhumation rate estimated for the northernmost segment of the Sparta Fault (Çal et al., 2024). Our data show an increase in average slip rate during exhumation of the scarp from an initial 0.8–0.9 mm yr⁻¹ between 6.5 and 7.7 kyr ago to 1.0 mm yr⁻¹ between 3.0 and 6.0 kyr ago (Fig. 5e). These slip rates directly reflect the steeper ³⁶Cl gradient for the lower 4.0 m of the fault scarp compared with the gentler gradient from 4.0 to 6.5 m (Figs. 2 and 3c). Although the sampling density is highest over the lowermost 4 m, we have confidence in the lower inferred average slip rate for the higher, older part of the scarp because both our dispersed drill core samples and the Benedetti et al. (2002) profile indicate a lower ³⁶Cl concentration gradient (in trend, rather than absolute values) above 4 m. The MAP model (Fig. 3a) indicates that three scarp-exhuming earthquakes may have occurred 5000–6000 years ago (MAP average slip rate of 1.1 mm yr⁻¹), which is consistent with an increase in average slip rate during this period observed in the slip rate versus time plot (Fig. 5e). The lower rate of exhumation for the upper ∼ 2.5 m reflects an apparent quiescent period prior to these earthquakes. Neither the historical record nor the ³⁶Cl concentrations measured on the Anogia A and B profiles (Figs. 3 and 4a), supported by measurements of geometry and hanging-wall colluvium, and in Benedetti et al. (2002) provide evidence for large, scarp-exhuming earthquakes after the 464 BCE event. Fault scarps may be exhumed by earthquakes clustered within several thousands of years and then lie dormant for similar, or even longer, periods (Wallace, 1987; Friedrich et al., 2003; Benedetti et al., 2013, Cowie et al., 2017a). The recent 2.5 kyr period of quiescence is therefore not necessarily indicative that another earthquake is imminent.

Our data do not uniquely specify the number and timing of scarp exhumation events, and we have been unable to identify other faults along the eastern flank of the Taygetos mountains suitable for ³⁶Cl analyses that with the Sparta Fault may form part of a system, across which slip is distributed. We therefore limit our interpretations to averaged slip rates and the timing of changes in these rates for the Sparta Fault at Anogia, rather than attempting to identify individual earthquakes or draw conclusions on regional fault kinematics and associated seismic hazards.

5.2 Methodological and geological sources of uncertainty in the ³⁶Cl data

A feature of the ³⁶Cl data is that our Anogia A and B profiles display systematically lower concentrations than the Benedetti et al. (2002) profile (Fig. 2). The Benedetti et al. (2002) profile also displays variations between adjacent sample points that exceed those observed in our profiles. We interpret the systematic differences in ³⁶Cl concentration between our profiles and the Benedetti et al. (2002) profile as reflecting methodological differences related to advances in sample preparation chemistry at PRIME Lab, Purdue University. For this reason, we elect not to model the Benedetti et al. (2002) data using the MCMC methodology.

While our Anogia A and B profiles display corresponding trends with increasing elevation on the fault scarp, Anogia B samples generally have lower ³⁶Cl concentrations (Fig. 2). Indeed, 6 of its 19 ³⁶Cl concentrations do not overlap within uncertainty with concentrations of corresponding samples on the 3.9 m Anogia A profile, including four points located between 1.0 and 1.3 m. We interpret these differences as indicating that the fault scarp at Anogia B has either been partly shielded from cosmogenic radiation, has eroded more than the scarp surface at Anogia A, or contains a higher concentration of non-calcite impurities (Sect. 5.3). Of potential additional relevance is that the texture of the scarp surface at Anogia B is smoother than at the location of Anogia A. Because a similarly smooth texture also characterizes a portion of the scarp surface presently buried by colluvium mantling the hanging wall, the smooth texture at the location of Anogia B may indicate recent burial of the scarp surface by colluvium and/or CaCO₃ dissolution/reprecipitation occurring at a higher rate than at locations where the exposed scarp surface texture is rougher. If a smooth texture reflects erosion through CaCO₃ dissolution, there might be preferential flow, or seepage, of water from the hillslope above the scarp at the location of Anogia B. Observed lumps of colluvium cemented to the Sparta Fault scarp, at locations perched above the present hanging-wall surface (Fig. S1), partially shield the underlying scarp surface today. However, had this previously occurred at the location of Anogia B, an eroded colluvial lump would be evidenced in the hanging-wall sediments. On the contrary, there is no colluvial lump, but rather a sub-horizontal surface is present with an expression that differs little from the surface below the Anogia A profile. The inter-profile differences in ³⁶Cl concentrations illustrate the value in taking samples for ³⁶Cl measurements from more than one vertical profile at a particular location because ³⁶Cl concentrations can vary either through spatial variations in non-calcite impurities or past shielding by sediments or bedrock, which can otherwise be difficult to detect. Partial shielding may impact the interpretation of palaeoseismicity, including the timing, number and magnitudes of earthquakes, through locally lowered ³⁶Cl concentrations.

5.3 The effects of mineralogical impurities on ³⁶Cl concentrations

Mineralogical impurities embedded in the fault breccia that comprise the scarp surface appear to be a key geological reason for spatial variations in the concentration of ³⁶Cl. Measurements of chemistry and mineralogy at Anogia B indicate that SiO₂ comprises 0.1 wt %–20.8 wt % of the scarp. Because the concentration of CaCO₃ is inversely correlated with SiO₂ (largely quartz), peaks in SiO₂ might coincide with troughs in ³⁶Cl, although a simple relationship vertically along the scarp is obscured by the relationship between ³⁶Cl concentration and exposure duration. A local peak in SiO₂ of 12 wt %–15 wt % coincides with a local low in ³⁶Cl concentration at Anogia B between about 0.6 and 1.2 m on the scarp (Figs. 2 and 7, Tables S1 and S3). A distinct low in ³⁶Cl concentration at 1.6 m also corresponds with a local peak in SiO₂ of 9 wt %. However, the magnitudes of the variations are inconsistent between these two locations, such that a high peak in SiO₂ corresponds with a small reduction in ³⁶Cl at 0.6–1.2 m and vice versa at 1.6 m. Because ³⁶Cl is also produced by spallation on K (162 ± 24 atoms g⁻¹ yr⁻¹ at SLHL; Evans et al., 1997), Fe (1.3±0.1–1.9 ± 0.2 atoms g⁻¹ yr⁻¹ at SLHL; Stone et al., 2005; Moore and Granger, 2019), and Ti (13 ± 3 atoms g⁻¹ yr⁻¹ at SLHL; Fink et al., 2000), noise in the ³⁶Cl data might also partly reflect the relative abundances of these elements. However, this appears to be insignificant given that measured concentrations of these elements are extremely low (concentrations of K₂O, Fe₂O₃, and TiO₂ are 0 %–0.12 %, 0.03 %–0.24 %, and 0 %–0.02 %, respectively; Fig. S4, Table S3). Other elements, seemingly present as trace amounts of clay, lining pores in the fault breccia (Fig. 7, Table S3), are also an insignificant contributor to variations in ³⁶Cl concentrations. For the Sparta Fault at Anogia, quartz embedded in the fault breccia may be the key mineralogical impurity that is likely contributing variance to the ³⁶Cl concentrations, which in turn impacts our ability to reliably date individual earthquakes.

5.4 Interpretation of REE-Y distributions and implications for palaeoseismicity

REE-Y cannot be used to infer imprints of former soil profiles on the Sparta Fault at Anogia. Petrographic analyses indicate that the Sparta Fault scarp is composed of a protocataclasite consisting of calcite clasts derived from the host limestone, microcrystalline calcite cement, and quartz (Figs. 6, 7). Furthermore, EDS analysis indicates that trace amounts of clay, such as illite, are lining pores where microcrystalline calcite cement and quartz are located (Fig. 8; Carcaillet et al., 2008). We infer that REE-Y are adsorbed onto clay minerals lining pores in the fine-grained matrix of the fault breccia, as indicated by correlations between REE-Y and Al, K, Si, and Fe (R²= 0.92, 0.87, 0.56, and 0.47, respectively; Fig. S4a–c) and between Y and both Si and Al in the hanging-wall colluvium (R²= 0.71 and 0.45, respectively; Figs. 12c and S4b, c). Supplementary data from the Kaparelli fault (R²= 0.95 for Si; Figs. 1a and S5a) and Magnola fault hanging walls (R²= 0.98 for both Si and Al; Fig. S5b, c and electronic appendix to Manighetti et al., 2010) also indicate REE-Y may be adsorbed to clay embedded in limestone fault scarps. These correlations generally contrast with a weaker negative correlation between Y and pH (R²= 0.52) for the hanging-wall soil on the Sparta Fault (Fig. 12c). Soil pH does not appear to be the dominant control on REE-Y distributions in the Sparta Fault scarp, which differs to interpretations on other limestone fault scarps (Carcaillet et al., 2008; Bello et al., 2023).

We propose a causative relationship between the vertical distributions of REE-Y and clay on the Sparta Fault scarp. This reasoning is supported by the following observations:

• i.

  The Sparta Fault scarp REE-Y concentrations are equivalent to (Nuriel et al., 2012; Goodfellow et al., 2017) or higher than those measured elsewhere in platformal limestone (Carcaillet et al., 2008; Mouslopoulou et al., 2011), but Y concentrations are lower in the adjacent hanging-wall soil (REE were not measured in the soil; Tables S6, S8).

• ii.

  If REE-Y exchange between the soil and fault scarp occurs according to the Carcaillet et al. (2008) model, fractionation of LREE and HREE elements is expected. For example, LREE might be preferentially mobilized (Takahashi et al., 2005; Carcaillet et al., 2008), leading to an enrichment of LREE relative to HREE in the fault scarp, where there are peaks in total REE-Y. Conversely, LREE may be depleted relative to HREE where there are troughs in total REE-Y. However, the proportion of LREE to HREE remains confined to a constant range vertically along the subaerial section of the Sparta Fault scarp (Figs. 10b, 11a), is weakly correlated with total REE-Y (R²= 0.36; Fig. 11b), and is relatively depleted at all measured depths beneath the soil surface (Fig. 10b).

• iii.

  There is no systematic decrease with distance above the hanging wall in total REE-Y (Fig. 10a, b), in contrast to declining concentrations with distance above the hanging wall on the Magnola fault (Carcaillet et al., 2008).

Adsorption of REE-Y onto clay has been observed in regolith (Borst et al., 2020) but has not been previously discussed in the context of interpreting palaeoseismicity on limestone fault scarps.

Although we infer that adsorption of REE-Y onto clay minerals embedded in fault breccia dominates on the Sparta Fault, the approximate coincidence of the subsurface peak in scarp $\mathrm{LREE}/\mathrm{HREE}$ and total REE-Y with the mid-profile peaks in soil pH and Y (Figs. 10, 12a, b) provides evidence of REE-Y exchange between the scarp and the soil. However, the consequence is LREE depletion in the scarp, rather than enrichment (Fig. 10b), and it is unclear why this apparent depletion is not replicated on the subaerially exposed scarp. One possibility is that colluvium accumulation postdates the most recent earthquake, although, if so, low ³⁶Cl concentrations in the buried scarp surface indicate that the soil accumulation was co-seismic with the last earthquake or accumulated soon afterwards. It is also unclear why colluvium would accumulate only after the most recent earthquake. An alternative possibility is that a superficial LREE-depleted zone has been eroded from the subaerial scarp surface through dissolution. This would imply erosion of centimetres of scarp surface since the last known earthquake on the Sparta Fault at 464 BCE (an erosion rate of 0.01 mm yr⁻¹ over the past 2500 years would remove 2.5 cm of scarp surface). However, another possibility is that perhaps more time is required to increase LREE to concentrations seen on the subaerial scarp surface, but 2500 years have already passed since the most recent known earthquake, and maximum REE-Y enrichment has been inferred to occur within 500 years on the Spilli and Magnola faults (Manighetti et al., 2010; Mouslopoulou et al., 2011). Alternatively, LREE enrichment occurs after scarp exhumation, perhaps through exchange with aeolian dust fallout, as has been observed in Dead Sea halite (Censi et al., 2023). Such dust inputs may supply REE-Y (Yang et al., 2007), as indicated by the correlation between Y and Si in the hanging-wall colluvium (Fig. 12b, c); contribute fine-grained mineral soil to the hanging-wall colluvium; and lower soil pH through buffering locally sourced CaCO₃. However, given that inputs of Saharan dust are ubiquitous throughout the Mediterranean (Stuut et al., 2009) and can comprise a large component of soils in the region (Muhs et al., 2010; Styllas et al., 2023), similar patterns of LREE depletion in the soil-covered scarp surface relative to the subaerial scarp surface are expected to have been observed elsewhere, which is not the case (Carcaillet et al., 2008; Manighetti et al., 2010; Mouslopoulou et al., 2011; Tesson et al., 2016; Bello et al., 2023).

For the Sparta Fault scarp, the presence of clay likely relates to fault breccia formation at considerable depths beneath the Earth's surface, rather than subaerial weathering processes. The formation of protocataclasite occurs beneath the Earth's surface at depths that may range from metres to up to thousands of metres. A model for this involves fluids moving along the Sparta Fault, primarily associated with seismic events. These fluids dissolve CaCO₃ from the host limestones and potentially also silicate minerals from psammitic and pelitic (meta)sediments, where they are dissected by the fault. In association with variations in temperature and pressure along the fault, chemical saturation of these fluids results in precipitation of clay, quartz, and microcrystalline calcite, which cements clasts of host-rock-derived limestone into the fault breccia. Subsequent faulting re-fractures the breccia, and particle comminution over time produces quartz grains that are rounded to angular in shape, randomly oriented, and < 50 µm (Figs. 6, 8). The fault breccia may also have undergone multiple generations of microcrystalline calcite re-cementing from re-circulating fluids. As an alternative to a dissolution–precipitation model, clay and quartz emplacement may involve fluid entrainment of particles and grains from clay- and quartz-bearing sedimentary units during faulting, as has been observed elsewhere (e.g. Darwin, 1891; Roy, 1946; Brandon, 1972; Röshoff and Cosgrove, 2002). This process may also be accompanied by comminution of fault-zone quartz grains derived from psammitic rocks. We tentatively exclude a contemporary aeolian source for the clay and quartz because there is no documented mechanism to transport clay particles and quartz grains from the soil to centimetres into a fault scarp. We cannot distinguish soil to scarp clay and quartz migrations on the Sparta Fault, which have been observed, for example, at the micrometre scale in surface coatings on the Magnola fault, because that scarp is comprised of pure carbonate (Carcaillet et al., 2008). It is likely that limestone fault scarps are generally composed of fault breccias (Agosta and Aydin, 2006; Carcaillet et al., 2008; Nuriel et al., 2012) and that, where a fault intersects varying lithologies, chemical and mineralogical heterogeneities may occur in the fault breccia, as observed on the Sparta Fault. Where they occur, these heterogeneities may control the spatial distribution of REE-Y, independent of any spatial reorganization of REE-Y attributable to subaerial weathering and pedogenesis. If, as we infer, the spatial patterning of REE-Y, quartz, and clay is inherited from depth, the observed wave-like signal (Figs. 7, 10) may reflect sorting and cementing of breccia around surface asperities on the fault plane. The resulting infilling of depressions with fault gauge may create a successively more polished and localized fault plane along which friction is lowered, thereby permitting larger slip (i.e. larger earthquakes) along the fault (Sagy and Brodsky, 2009). While REE-Y concentrations do not appear to be a reliable indicator of Holocene palaeoseismicity of the Sparta Fault, they may instead reveal processes that localize slip to a discrete fault plane.

While the Sparta Fault displays concentrations of clay and quartz impurities that are much higher than on other reported limestone fault scarps, three general implications emerge for using REE-Y in making inferences on palaeoseismicity. Firstly, the potential control on REE-Y distributions of even trace amounts of non-calcite impurities in the breccia comprising fault scarps should be considered through analyses of thin sections in addition to scarp chemistry. Secondly, soil acidity and REE-Y enrichment, including any resulting exchange with the buried scarp, may peak some tens of centimetres below the colluvium surface. Peaks in REE-Y concentrations on subaerial fault scarp surfaces may therefore not reflect former soil surfaces, even if there is soil–scarp exchange of REE-Y. In addition, the Sparta Fault scarp REE-Y data indicate that it may be rewarding to focus on up-scarp variations in $\mathrm{LREE}/\mathrm{HREE}$ ratios, rather than on REE-Y concentrations, because these may be a sensitive indicator of REE-Y exchange processes occurring beneath soil covers (Fig. 10b). Lastly, relationships between REE-Y distributions and soil mineralogy should be more closely assessed, in addition to the commonly modelled and studied effects of pH (e.g. Carcaillet et al., 2008; Manighetti et al., 2010; Mouslopoulou et al., 2011; Moraetis et al., 2023; Tesson et al., 2016; Bello et al., 2023). Fine-grained mineral inputs through aeolian dust fallout comprise substantial volumes of Mediterranean soils (Muhs et al., 2010; Styllas et al., 2023). These inputs may be complemented in the eastern Mediterranean by Holocene tephra from the South Aegean Active Volcanic Arc or volcanic centres in Italy (Bourne et al., 2010; Smith et al., 2011; Koutrouli et al., 2018; Vougioukalakis et al., 2019). Decadal to millennial variations in Holocene dust and tephra fluxes may directly impact on REE-Y distributions in hanging-wall soils and potentially in scarp surfaces, in locations where soil–scarp REE-Y exchange is important. These fluctuations may contribute to REE-Y patterns in soils that are difficult to predict and in scarp surfaces reflect (volcanic, climatic, and pedogenic) processes that may complicate potential palaeoseismic inferences.

Moraetis et al. (2023) consider REE-Y analyses an established method in palaeoseismicity. Our detailed study errs towards caution; there remain important uncertainties regarding processes of REE-Y enrichment and depletion in limestone fault scarps. Indeed, we maintain that there is considerable uncertainty regarding how the resulting patterns should be interpreted with respect to palaeoseismicity. Fundamentally, it remains unclear how far into buried scarp surfaces the REE-Y can be adsorbed from soil or incorporated into calcite through dissolution–precipitation. A dissolution rate of 0.001 mm yr⁻¹ will erode 1 cm from a subaerially exposed scarp surface over 10 000 years, which is about the timescale considered to be relevant to assessing full seismic cycles and therefore making accurate assessments of palaeoseismicity (Mouslopoulou et al., 2012; Tesson et al., 2016). Even such a slow rate of subaerial scarp dissolution will therefore remove any REE-Y signals inherited from former soil cover unless that exchange extends to centimetres into the scarp.

6 Conclusion

Modelling of slip rates from ³⁶Cl data from the Sparta Fault at Anogia, Greece, indicates an increase in average slip rate during exhumation of the scarp from 0.8–0.9 mm yr⁻¹ between 7.7 and 6.5 kyr ago to 1.0 mm yr⁻¹ between 6.5 and 2.5 kyr ago (the timing of the 464 BCE earthquake). Average exhumation of the entire scarp is 0.7–0.8 mm yr⁻¹. Earthquake ages were not modelled from our data, but there is no indication from our analyses that earthquakes may have contributed to exhumation of the Sparta Fault since 464 BCE.

The Sparta Fault scarp is impure; it is composed of fault breccia, which contains quartz and clay-lined pores in addition to calcite. The vertical distribution of REE-Y is highly correlated with the pore clay and may indicate processes that localize slip to a discrete fault plane deep below the ground surface. The potential exchange of REE-Y between the hanging-wall colluvium and the adjacent footwall scarp is overwhelmed at this site by REE-Y attached to the pore clays inherited from depth. Because of this, Holocene earthquakes and their slip distances and magnitudes cannot be inferred for the Sparta Fault from REE-Y concentrations. While this is probably also true for similar impure limestone fault scarps elsewhere, other controls on REE-Y distributions, in addition to hanging-wall soil pH, should be evaluated in attempting palaeoseismic inferences more generally from normal fault scarps developed in limestone.