Coupling rare earth element analyses and high-resolution topography along fault scarps to investigate past earthquakes: A case study from the Southern Apennines (Italy)

Italy is among the most seismically active areas of the entire Mediterranean basin, and it hosts extensional, compressional, and strike-slip earthquakes (e.g., de Nardis et al., 2022; Di Luccio et al., 2005; Herrmann et al., 2011; Pondrelli et al., 2006; Scognamiglio et al., 2006). Major extensional earthquakes concentrate along the Apennine chain in the so-called extensional seismogenic province (sensu Lavecchia et al., 1994, 2021). This province extends with an average NW-SE direction along the axial zone of the Apennines (Galadini and Galli, 2000; Lavecchia, 1988) and continues with a N-S strike in northern Calabria (Andrenacci et al., 2023; Brozzetti et al., 2009, 2017a, 2017b; Cirillo et al., 2022b; Napolitano et al., 2021; Tortorici et al., 1995) and a NNE-SSW strike in southern Calabria (Cucci, 2022; Galli and Bosi, 2003; Neri et al., 2020). In this active geological environment, the faults often form rock scarps, facilitated by Bahamian-type Mesozoic–Cenozoic carbonate rocks. These well-preserved rock scarps may reach up to tens of meters in height and stretch along the edges of great intramontane Quaternary basins for hundreds to thousands of meters (Brozzetti, 2011; Brozzetti et al., 2019, 2020; Bucci et al., 2013; Civico et al., 2018; Galadini and Galli, 2000; Galli et al., 2006; Galli and Peronace, 2014; Gori et al., 2011; Mirabella et al., 2011; Roberts and Michetti, 2004; Schirripa Spagnolo et al., 2021; Sgambato et al., 2020; Villani et al., 2018; Villani and Pierdominici, 2010). In some cases, these structures run in the highest portions of the main ridges, sometimes even with trends oblique to the ridges, where they form narrow intramontane along-strike elongated basins, some hundreds of meters long and filled with late Quaternary clastic deposits. This is the case for the 25-km-long Caggiano fault system (Mt. San Giacomo and Timpa del Vento segments in Galli et al., 2006), which, according to Bello et al. (2022b), could be the northern segment of a much longer system (Caggiano-Montemurro fault). The Caggiano-Montemurro fault is an ~65-km-long, SSW-dipping trans-ridge fault crossing a portion of the Monti della Maddalena ridge (Fig. 1), connecting the Auletta and Val d’Agri basins and bordering the northernmost portion of the Vallo di Diano (Fig. 1). The Caggiano-Montemurro fault was recently hypothesized to be the main fault of the area, with a length capable of generating the 1857 CE (M_w 7.1) earthquake (Bello et al., 2022b). In turn, previous paleoseismological investigations across the Caggiano fault by Galli et al. (2006, 2008) accounted for several surface-faulting events in the past 6 k.y. The latest ones occurred in the last 2 k.y. B.P., probably during/after slope-debris deposition related to the Little Ice Age (ca. 1400–1800 A.D.). According to Castelli et al. (2008), both of the fault segments comprising the Caggiano fault were responsible for the two 1561 CE main shocks (31 July and 19 August 1561 CE; cumulative M_w = 6.7 in the Catalogo Parametrico dei Terremoti Italiani [CPTI15]), and, tentatively, also the northernmost shock of the devastating 1857 CE earthquake sequence (cumulative M_w = 7.0).

In recent years, many advances have been made to study coseismic deformation at the surface using high-resolution topography (e.g., Bello et al., 2021b; Bubeck et al., 2015; Cirillo, 2020; Johnson et al., 2014; Westoby et al., 2012; Wilkinson et al., 2015, 2016). In fact, high-resolution imagery and topography from light detection and ranging (LiDAR) and unmanned aerial vehicle (UAV) technologies provide excellent data sets for measuring fault scarps (e.g., Bello et al., 2022a; Scott et al., 2022; Stewart et al., 2018; Wilkinson et al., 2015; Wolfe et al., 2020). The high-precision measurements obtainable allow more in-depth studies, pushing forward the frontier of knowledge on old topics, like past earthquakes. Small UAVs can be equipped with onboard real-time kinematic global navigation satellite system (GNSS/RTK) antennas working with the post-processing kinematic (PPK) technique. The latter allows the acquisition of georeferenced photographs (centimeter-scale resolution) without placing ground-control points (Bolkas, 2019; Cirillo et al., 2022a; Cledat et al., 2020; Vichi et al., 2022; Zhang et al., 2019).

With the help of a MATLAB algorithm, a recently developed technique allows semi-automatic measurement of the vertical displacement along the active fault strike (Bello et al., 2021b; Scott et al., 2020). This code, which permits the measurement without losing control by the geoscientist, is valuable to constrain both the vertical displacement due to coseismic ruptures and the vertical displacement of long-term scarps. Although this approach is extremely useful because it provides an estimate of the scarp height potentially attributable to deformation since the Last Glacial Maximum (ca. 18 ka; post-LGM hereafter), it is not able to provide details of the number of earthquakes and the displacement likely associated with each of them.

A new approach that could allow researchers to quantitatively hypothesize the number of past earthquakes recorded by a fault scarp is to analyze, from a geochemical point of view, the concentration of certain chemical elements in the fault planes. In fact, carbonate fault planes react geochemically with the environment to which they are subjected, being enriched or depleted in certain chemical elements as a function of the exposure time. In particular, Benedetti et al. (2002, 2013), Palumbo et al. (2004), Schlagenhauf et al. (2011), and Pousse-Beltran et al. (2022) investigated the postglacial earthquake record by analyzing bedrock fault planes using ³⁶Cl cosmogenic dating. The methodology was applied on the Sparta fault (Greece) by Benedetti et al. (2002), on the Fucino basin fault (Italy) by Benedetti et al. (2013), on the Magnola fault (Italy) by Palumbo et al. (2004) and Schlagenhauf et al. (2011), and on the Mt. Vettore fault (Italy) by Pousse-Beltran et al. (2022). Recently, Iezzi et al. (2021) combined a geodetic survey with the ³⁶Cl methodology on active normal faults near Athens, Greece. Carcaillet et al. (2008) and Mouslopoulou et al. (2011) were among the first to analyze the content of rare earth elements (REEs) in limestone fault planes to constrain the number of earthquakes released in post-LGM times by the Magnola and the Spili fault (Greece), respectively. Finally, the ³⁶Cl dating and REE analyses were combined by Tesson et al. (2016) on the Pizzalto fault, in central Italy.

Both approaches (i.e., topographic displacement and geochemical analyses) are based on the same two principles: (1) Since the Last Glacial Maximum (ca. 18 ka; Galli et al., 2006), the effects of surface faulting on rocky fault scarps are preserved due to the higher rates of faulting compared to the effects of erosion/sedimentation and (2) any earthquake of such magnitude as to allow the rupture to propagate at the surface will exhume a new portion of the fault plane, adding it to the cumulative fault scarp. For this last reason, the fault plane is ideally divided into ribbons of variable height, each of which would represent a seismic event. By sampling portions of rock parallel to the slip direction, a “log” of the slip increments can be recorded, and geochemical analyses can then be carried out on this record. Of course, fault planes can also be exhumed by non-tectonic (erosional and landslide) processes, and thus a geological and paleoseismological check must first filter out false active faults from real ones.

In this study, we concentrated on the northern section of the Caggiano-Montemurro fault. This highly segmented “trans-ridge fault” shows a well-defined en échelon geometric pattern at multiple viewing scales. Following the most widely adopted criteria in the literature (i.e., by identifying geometric and structural complexities such as gaps, overlaps, underlaps, and sudden strike variations), Bello et al. (2022b) constrained four orders of fault surface segmentation for the Caggiano-Montemurro fault. The section on which we concentrated, called Timpa del Vento, exposes fault planes that are suitable for our analyses because of their lithology, the quality of the outcrops, and the logistical ease of reaching them with instruments (Figs. 1–3). We combined the methodology of the displacement discontinuity approach (measuring displacements with high-resolution topography) to obtain information on ~5 km of fault length along with a detailed geochemical study at one site defining the variations in the content of trace elements (e.g., REEs + Y) on the fault plane.

The goal was to assess the number of earthquakes potentially released by the Timpa del Vento fault section, and likely by the Caggiano-Montemurro fault, which is considered to be a segmented or branched-at-depth normal fault (Bello et al., 2022b; Cello et al., 2003; Soliva et al., 2008). Finally, we used empirical relationships to estimate earthquakes magnitude from fault slip (Leonard, 2010; Thingbaijam et al., 2017; Wells and Coppersmith, 1994; Wesnousky, 2008).

The study area is located in the axial sector of the late Miocene–Pliocene Southern Apennines fold-and-thrust belt, which was dissected in the Quaternary by extensional structures (Lavecchia et al., 1994, 2021, 2022; Brozzetti, 2011, and references therein). We refer the interested reader to the broad available literature (e.g., Mostardini and Merlini, 1986; Nicolai and Gambini, 2007; Scrocca et al., 2007; Vezzani et al., 2010) for a complete and detailed description of the paleogeographic domains involved in the construction of the mountain chain. The extensional fault system offsets preexisting compressive structures of the fold-and-thrust belt beginning in the early Pleistocene (e.g., Barchi et al., 2007; Brozzetti, 2011; Casciello et al., 2006; Di Naccio et al., 2013; Papanikolaou and Roberts, 2007; Schiattarella, 1998), with low to moderate slip rates (Galli, 2020; Galli et al., 2006; Papanikolaou and Roberts, 2007; Sgambato et al., 2020) generating intramontane basins filled with Quaternary clastic and fluvio-lacustrine deposits (Amicucci et al., 2008; Brozzetti et al., 2017b; Brozzetti and Salvatore, 2005; Bucci et al., 2020; Di Giulio et al., 2016; Robustelli et al., 2014; Villani and Pierdominici, 2010; Villani et al., 2019). In the study area, these are the Auletta, the Vallo di Diano, and, marginally, the Val d’Agri (Fig. 1) basins, which show southwestward-thickening tilted beds in the Auletta basin (Amicucci et al., 2008; Brozzetti, 2011) and northeastward-thickening tilted beds in the Vallo di Diano and Val d’Agri basins (Giano et al., 2000; Zembo et al., 2009). For the whole area, current stress field data from tectonic structures show NNE-SSW extension, with an approximately N032E-trending near-horizontal s3 axis (Bello et al., 2021a, 2022b), slightly deviating from the approximately NE-SW s3 axis highlighted further south (e.g., Brozzetti et al., 2009; Cirillo et al., 2022b) or further north (e.g., Castaldo et al., 2018; Ferrarini et al., 2015; Lavecchia et al., 2017; Mariucci and Montone, 2016) along the Apennine chain.

The study area framed in Figure 2 is located along the Monti della Maddalena ridge (Fig. 1B), which is the structural continuity of the Mt. Marzano carbonate massif, the latter of which is bordered by the Irpinia fault system (Fig. 1A). This ridge is made of Triassic to Cretaceous marine carbonate and dolomitic rocks and Cretaceous–Paleogenic limestones, capped by Miocene siliciclastic deposits and by arenaceous flyschoid formations (Pescatore et al., 1970). The study sites are characterized by the presence of a portion of the Caggiano fault, which was recognized to be active by Galli et al. (2006). At the footwall of the fault, beyond the cataclastic layer, the rock type is characterized by shelf-marginal facies consisting of massive or poorly bedded bioclastic rudstones and framestones. The rock fault scarp can be followed along the narrow (200–400 m) Timpa del Vento and Campo di Venere basins, where it also affects active decametric-scale alluvial fans (Galli et al., 2006). The carbonate fault plane of our sampling site is located uphill of the southeastern tip of the Timpa del Vento basin, a few hundred meters from the beginning of a bedrock saddle dividing the two basins (Fig. 2). The sampling site is characterized by a tectonically exhumed fault plane portion (higher portion of ~370 cm) plus a lower part that was recently (between 2006 and 2016) bulldozed for the construction of a dirt road. This allowed us to reach the fault plane with vehicles and to carry the necessary sampling instrumentation close to the plane. In this site, Holocene deposits cut and are in contact with the fault, and, thanks to recent excavation, it is possible to correlate the modern soil covering the fault plane until recently with the REE concentrations. Close to the site described here, the topographic depression deepens abruptly and begins to widen to reach, within 500–600 m, its widest portion. According to some boreholes described in Galli et al. (2006), the infilling deposits of the basins are represented by a few meters of clays and sands along with alternating layers of slope debris and volcanic ashes. These deposits represent the last episode of a marsh environment related to the damming of the basin controlled by the fault. Below these continental deposits, boreholes showed ~15 m of yellowish and blue-gray sandy clays, attributable to the marine Pliocene section by Lucchetti (1943) and covering the bedrock.

In April 2022, we acquired 2677 photographs along part of the Timpa del Vento fault using a DJI Mavic 2 Pro drone and a Phantom 4 Pro V2 drone that flew at ~50–100 m altitude above ground level (Fig. 2A, inset A′). For the best positioning of the photographs, we used an Emlid Reach RS2 GNSS/RTK L1, L2, L5 base station, positioned in two central portions of the study area, and an antenna rover (L1/L2 RTK/PPK) installed on board both drones. In total, we covered ~5 km along the fault strike of the Timpa del Vento segment with an average imaging width of ~500 m (Fig. 4A). The acquired photographs were reviewed, and then those that were low quality, blurred, or taken during takeoff and landing were eliminated. We processed the photographs to produce dense point clouds, orthomosaics, and digital elevation models (DEMs) with the Agisoft Metashape Pro image-based photogrammetric modeling software (version 1.8.4; e.g., Bello et al., 2021b; James and Robson, 2012; Johnson et al., 2014; Westoby et al., 2012). We aligned the 2619 selected photographs and produced the models with high-quality settings that were exported with the default recommended resolution. The obtained accuracy of the models was ~3 cm in the horizontal and vertical directions. We built hillshade maps in ArcMap (ESRI ArcMap© 10.8) with a resolution of ~3 cm/pixel. The DEM we produced ranged from ~1000 to ~1200 m altitude above sea level and covered an area of ~3.2 km².

We mapped the Timpa del Vento fault at a fixed scale of 1:500 in ArcMap©, using the DEMs and orthomosaics produced as base maps (Fig. 4A). The fault trace is observable at the border of the small (i.e., 1–3 km²) intramontane Quaternary basins for the entire extent of the DEM (Figs. 3 and 4A). As mentioned, the fault pattern is highly segmented, and, even at the large scale (i.e., 1–3 km), it is articulated in small portions that border small intramontane basins (Figs. 2A and 2B). Based on the mapped traces and the DEMs produced, we used a MATLAB algorithm developed by Scott et al. (2020) for use by Bello et al. (2021b) and investigated the topography to measure the vertical separation (VS) along 98 topographic profiles. This approach allowed us to measure the VS defined as the vertical distance between the linear surface projections at the hanging wall and footwall of the fault measured at the fault trace position. For normal to oblique-normal faults (i.e., dip-slip or dip-slip with a slight strike-slip component), this assumption is considered appropriate because the lateral slip is thought to have minimal influence on the total displacement (Bello et al., 2021b; DuRoss et al., 2019; Scott et al., 2022). We reference the reader to Bello et al. (2021b) and to the user guide in Scott et al. (2020) for a detailed explanation of the steps for preparing the input files for the measurement process and a complete user guide of the MATLAB code.

We generated the topographic profiles, orientated perpendicular to the average fault strike (Fig. 4A), from the DEM with a 50 m spacing and a 2 m averaging window, which minimized the impact of the topography. To acquire the measurements, we manually marked two points on the hanging wall and two points on the footwall of the fault to provide the surface projection. The four points for the linear surface projections must be selected on bare ground, avoiding vegetation, which was easily identifiable on dense clouds, orthomosaics, DEMs, and topographic profiles (see details in Bello et al., 2021b; Scott et al., 2020). The position of the fault along the topographic profile was also manually picked by the operator and generally identified in the steepest part of the scarp face. The algorithm automatically measured the vertical separation. An example of one of the topographic profiles with the projection lines is shown in Figure 4B. The measurement was accompanied by a quality factor (measure quality ranking [MQR]), which characterizes the quality of the measurement based on parameters such as vegetation, the angle between the two linear surface projections (more or less regularized slope), and the ease of identification of the fault trace position. In addition, the measurements were accompanied by a statistical error, which was automatically calculated by the code as in the following equation by Bello et al. (2021b):

where m_footwall and m_{hanging wall} represent the best-fit slopes interpolated at the footwall and at the hanging-wall lines of the fault, while b_footwall and b_{hanging wall} are the topographic intercepts. ΔFx is 25% of the scarp height, representing the error in fault positioning (for details on the best-fit slopes calculations, see Bello et al., 2021b).

In total, we acquired 86 measurements on 76 transects considering the fault scarp near field, plus a further six measurements on six profiles considering the fault far field (i.e., the entire slope) and assuming a rough slope rectification during the strong LGM rhexistasy phases (ca. 24 ± 3 ka; Fig. S1¹; Galli et al., 2012, 2017, 2022). Profiles and measurements will be described in the next sections.

Sedimentary limestones are generally REE-poor (e.g., Rosatelli et al., 2023, and references therein; Stoppa et al., 2021), while the pedogenic environment (soil) they are in contact with is generally enriched during its development. The contents of REEs and other high field strength elements (HFSEs) in soils are highly variable, depending on the properties of the parent material, the degree of weathering, the content of organic matter, biological activity, and the presence of clay minerals, Fe/Mn oxides, phosphates, and hydroxides (Cao et al., 2001; Kabata-Pendias, 2010; Tangari et al., 2018, 2021; Tyler, 2004). In the uppermost part of the soil, REE forms organic complexes (Pourret et al., 2007), preventing their leaching and transfer to the deeper part of the soil. Thus, REE + Y enriches the organic-richer soil portion (Carcaillet et al., 2008, and references therein). Primary (i.e., limestone) and secondary (i.e., cement) carbonates easily dissolve in the presence of any fluids acidified by dissolved CO₂. Fault planes developing in limestone produce breccias subject to dissolution by rain and circulating (vadose) waters. As explained in detail by Carcaillet et al. (2008), the dissolution of carbonate along an exposed fault plane produces REE + Y enrichment in the runoff waters. These waters reach the pedogenized colluvial wedge at the base of the fault scarp, where REE + Y forms organic complexes (Pourret et al., 2007) and/or is taken up by specific bacteria (Chen et al., 2000), enriching the topsoil and its fluids. These fluids produce reprecipitation of carbonates on the fault plane. Elements with similar ionic radius and valence similar to Ca²⁺ (e.g., Mg, Sr, P, Mn⁺²; Carcaillet et al., 2008, and references therein) may enter the new carbonate crystal lattice. Other elements, including REE + Y, however, are not readily incorporated into the cement (Rosatelli et al., 2010), but they are fixed by bacterial, algae, and root activity along grain boundaries (Carcaillet et al., 2008). Therefore, the REE + Y concentration peaks during quiescent earthquake periods in the fault scarp portion in contact with topsoil. The estimated time for the REE + Y enrichment in a limestone fault scarp is 0.7 k.y. (Carcaillet et al., 2008). Once the enrichment is established, it remains constant regardless of the length of the time of the soil–fault scarp contact (see residence time in Carcaillet et al., 2008). Therefore, there is not a linear correlation between residence time (earthquake interval) and REE + Y concentration unless the earthquake interval is less than 0.7 k.y. In addition, the REE + Y enrichment processes are linked to organic complexation and soil activity, so variations in environmental parameters affecting both weathering processes and soil organic production would affect the rate of REE + Y buildup. The exchange process stops at the time when the fault plane is exposed due to exhumation (slope erosion processes or surface faulting). The conceptual model in Figure 5 shows the principle of a carbonate bedrock fault plane enriched with REEs by modern soil.

The samples for the bulk chemical analyses were acquired from a carbonate fault plane outcropping along the Timpa del Vento fault (Figs. 3D, 6, and 7A). For sample extraction, we used a 900 W, 115 mm grinder on which we installed a 3-cm-diameter and 4-cm-long diamond coring cup powered by an electric-powered generator (Fig. 6). On the fault mirror, the kinematic indicators (i.e., slickenlines and calcite steps) show a sense of motion with a slight right-lateral component (rake = −116.3°, in Aki-Richard’s format; see the stereoplot in Fig. 2A), which we considered as the correct direction line for the sampling. We laid a wire tied to two nails, positioned one on the head and the other at the base of the scarp, parallel to the slip direction (considering the slickenside). Starting from the highest part, we cored the fault plane every 10 cm along the wire, in the best-preserved points (e.g., without pervasive fractures in the plane and vegetation), acquiring 37 limestone cores (Fig. 7B). We continued coring the fault plane every 20 cm in the lower portion, which had been recently bulldozed, acquiring further two samples (Fig. 7B). To have an external control over the subsequent analyses (samples as background to be compared with samples from the fault plane), we sampled three additional cores (hereinafter referred to as Timpa limestone background [LIM]) from an outcrop belonging to the same geological formation (Lower Cretaceous–Liassic carbonate platform deposits belonging to the Monte Marzano–Monti della Maddalena unit; sensu Pescatore et al., 1999), cropping out at the footwall of the fault (location in Fig. 3B) and therefore not in contact with any soil or colluvium. After sampling, we classified each sample, giving it an identifier, taking a photograph of the extraction location, and describing its characteristics such as the presence of more or less developed lichens on the exposed surface. The next steps took place in the laboratory.

The 39 limestone sample cores from the scarp and the three Timpa limestone background (LIM1, LIM2, LIM3) samples were prepared for inductively coupled plasma–mass spectrometry (ICP-MS) analysis. The first step was washing with 5 mL of 99.8% acetic acid to remove the dust from the drilling. We used a mini-drill to remove the surface patina affected by lichens and weathering. Then, the mini-drill was used to collect a 2 mm portion of carbonates from the external surface of the cores. Special care was taken to avoid as much as possible the limestone grains in the fault breccia and collect the carbonate cement among them. This powder was dried at 105 °C and weighed and then washed with 5 mL of ultrapure 130 vol% hydrogen peroxide for removal of organic matter, dried, and weighed again. The powder was dissolved with 5 mL of 2 N ultrapure hydrochloric acid at 80 °C, assuring the dissolution of calcite and dolomite. The dried solution was diluted with 5 mL of ultrapure nitric acid at 65%. The method adopted here was similar to that of Carcaillet et al. (2008) and Mouslopoulou et al. (2011) but adapted to our sample suite.

Each solution was subjected to analysis by ICP-MS, looking for trace-element concentrations. The instrument used was an Agilent 7900 ICP-MS (Agilent Technologies, Tokyo, Japan) in the laboratory of newborn screening, proteomics, and endocrinology of the Center for Advanced Studies and Technology (CAST), University of Chieti. Detailed operating conditions and instrumental parameters are given in Table 1. The fourth-generation Octopole Reaction System (ORS) was able to measure all elements in helium (He) mode, even though low-mass elements are normally measured in “no gas” mode due to the lack of interferences (Balcaen et al., 2015). The optimization of ICP-MS was carried out to obtain maximum signal intensities for ⁷Li, ⁸⁹Y, ¹⁴⁰Ce, and ²⁰⁵Tl using a 1 μg L⁻¹ tuning solution containing Li, Y, Co, Ce, Mg, and Tl (Agilent Technologies, Palo Alto, California), while keeping the formation of oxides ¹⁴⁰CeO⁺/¹⁴⁰Ce⁺ and doubly charged species Ce2⁺/Ce⁺ ratios below 1% and 3%, respectively. The sample introduction system was washed between different analyses with 2% HNO₃. Three multi-element mixtures at 10 μg mL⁻¹ were used in acid solution:

• (A) Ag, Ba, Be, Cd, Co, Cr, Cu, Fe, Mn, Ni, Pb, Rb, Se, Sr, Tl, U, V, and Zn in 5% HNO₃;

• (B) Ce, Dy, Er, Eu, Ga, Gd, Ho, In, La, Lu, Nb, Nd, Pr, Sm, Th, Tb, Tm, Y, and Yb in 5% HNO₃; and

• (C) Hf, Nb, Sn, Ta, and Zr in 5% HNO₃.

These mixtures were employed to prepare daily diluted calibration solutions, and three calibration curves were prepared using these multi-element mixtures. An internal standard correction was performed by online addition of an internal standard solution of Rh (50 mg L^–1) in a T-piece. Ultrapure water (18 MΩ cm⁻¹) was obtained from a Milli-Q system (Millipore, Bedford, Massachusetts). Nitric acid (69% v/v) and an internal standard solution of Rh were bought from Merk (Darmstadt, Germany) and were ultrapure grade. The full data set was recorded with Agilent MassHunter Data Acquisition software (v. 4.2) and processed with the Agilent MassHunter Data Analysis software (v. 4.2).

At our sampling site (Figs. 3D, 6, and 7B), the anthropic excavation aside the slickenside exhumed one brownish paleosol, truncated upward and mantled by a brownish pedogenized level, both developed on a volcanic layer. We sampled both levels (see Fig. S2) for radiocarbon dating to constrain the age by radiocarbon analyses. The two samples were named TDV18–01 (upper one) and TDV18–02 (lower one). The analyses were performed on bulk soil samples using the accelerator mass spectrometry (AMS) technique at the Beta Analytic (Miami, Florida) laboratories. Samples were physically pretreated to remove roots or macrofossils and acid washed to remove carbonates. The measured radiocarbon ages were calibrated using the software CALIB 8.2 (Stuiver et al., 2021); for further details and for standards and analytical protocols see details available at http://www.radiocarbon.com/.

The analysis of the topographic profiles through the MATLAB algorithm (see Methods section) allowed the measurement of the vertical displacement on 76 topographic profiles out of 98 (Figs. 4 and 8A). The 22 profiles on which we did not acquire measurements were affected by complexities including man-made structures, such as roads, irrigation canals, quarries, or small landslides, and overdeveloped vegetation. The topographic profiles were traced every 50 m, roughly perpendicular to the main fault. On the profiles, we acquired 86 measurements of VS (Table 2; Fig. 8B). To characterize the long-term deformation and the variation of the fault scarps mapped in the study area, we plotted the VS data on a section parallel to the fault strike (Figs. 8B, 8C, and 8D). In Figure 8B, we plotted all the acquired measurements, while in Figures 8C and 8D, we considered the MQRs assigned during the measurement phase, and thereby we did not plot the quality = 4 values (i.e., low quality). We then used statistical errors (see Methods section) to add error bars to each measurement on the graph. Following Bello et al. (2022a), we computed envelope curves using the “envelope” function of the Signal Processing Toolbox in MATLAB (http://www.mathworks.com), which traces the envelopes of the peaks of the input signal, interpolating over local maxima separated by a given sampling interval, for which the units are the units of the data input. In Figure 8D, the sampling intervals are 2 and 30 m. We plotted the data both by differentiating the measurements based on the continuity of the mapped scarps (i.e., not cumulative in overlapping portions; Fig. 8C) and by cumulative values (Fig. 8D).

The displacement profile is characterized by a complex along-strike sinuosity with variable displacement, which overall follows a fairly regular trend. The displacement values (only those with MQR = 1–3) reach a maximum of ~400 cm and a minimum of around 70 cm, with mean and median values of ~210 cm. There are two major peak areas, a milder one located between 700 and 1800 m from the origin, and a more pronounced one located between 3300 and 3800 m from the origin. The latter area matches a portion of the fault where two en échelon strands overlap for ~500 m. Overall, the data show a two-peak profile (dark red dashed line in Fig. 8D) with two large slightly asymmetrical enveloping “bells” that grow rapidly moving from NW to SE and gradually decrease from the peak area southeastward. Considering the topographic profiles shown in Figure S1, i.e., considering the entire slope from the far field with respect to the fault location, we observed post-LGM vertical displacements ranging between ~9 and ~18 m, with average values of ~14 m.

In Table 3, we reported the ICP-MS analyses for the 39 samples picked from the slickenside and the 3 LIM representing the marine bedrock. As mentioned in the Sampling Strategy section above, LIM are samples acquired at the footwall of the fault plane on the same geological formation to be used as background. They had an average ΣREE + Y content of 2.56 (±0.56 standard deviation) ppm, with a La/LuN ratio (normalized to C1 chondrite of Sun and McDonough, 1989) of 5.58 (±2.02) (Fig. 9). LIM showed a very slightly negative Ce anomaly, prominent positive Eu anomaly, and variable Y content (Fig. 10). The Y variability is typical of limestones and can be related to sedimentation of Mn crusts with phosphate layers in the basin (Auer et al., 2016). The small Ce negative anomalies and the enhanced Eu positive anomalies, which determine high Eu/Eu* ratios (= Eu/√[Sm + Gd]; 12.2 ± 3.3 on average; Table 3; Fig. 9), can be related to the oxidation state of the waters (Madhavaraju et al., 2016). The bedrock LIM samples had V/Cr ratios between 25.7 and 30.6 (Table 3), i.e., well above the value of 4.3, which is the upper limit of suboxic conditions. Thus, the V/Cr ratio of LIM indicates very reducing conditions in the limestone formation environment. The poorly oxygenated waters favored the reduction of Eu³⁺ to Eu²⁺. Eu²⁺ is more compatible than REE³⁺ and so more prone to enter into the forming carbonate lattice.

Carbonates sampled on the fault plane (from S1 to S37 in Fig. 7; Table 3) above the present soil had very variable REE + Y contents with an average of 14.51 (±9.15) ppm (Fig. 10) and very variable La/Lu_N values with an average of 17.87 (±13.55) ppm (Fig 9). However, the ratio of light to heavy REEs (LREE/HREE) in the fault carbonate cements was higher than that in the Timpa limestone. An Eu anomaly was present in many but not in all of the fault carbonate samples. Where it was present, it was less pronounced than in the bedrock limestones, and the Eu/Eu* ratio was 1.61 (±1.2) on average (Fig. 9). Carbonate samples on the fault plane covered by the active soil (S38 to S39; Fig. 7; Table 3) were characterized by high REE + Y contents with an average of 23.43 (±6.44) ppm (Fig. 10), high La/Lu_N values between 8.47 and 11.61 (Fig. 9), and a low Eu/Eu* ratio with an average of 2.66 (±0.72) (Fig. 9). The data indicate that the soil at the contact with the fault plane scavenges REE + Y from the limestone breccia and deposits carbonate cements retaining some of the LIM characteristics, such as a negative Eu anomaly. However, the more oxygenated soil environment produces oxidation of Eu²⁺ into Eu³⁺, which can be incorporated into the newly formed cements with the same compatibility as the other REE³⁺ elements, reducing the Eu/Eu* ratio in those cements (Fig. 9). In the carbonate cements, the higher fractionation between LREEs and HREEs due to epigenetic crystallization is also evident (Fig. 9).

Figure 11 shows the variability of REE-Y concentration as a function of fault scarp height (cm). The vertical axis represents Δi/Δm, where Δi is the difference in concentration between each element and the average concentration of that element along the entire scarp (Ci – Cm), and Δm = Σ|Δi|N, with N = number of measurements (39) (e.g., Carcaillet et al., 2008; Mouslopoulou et al., 2011). The interpretation and discussion of the possible earthquakes highlighted by the REE-Y variations shown in Figure 11 are given in the next section.

We show in Figure S3 the detrended distribution (i.e., the difference between each abundance distribution and its trend) of REE-Y contents for each element, performed using MATLAB functions. REEs do not show homogeneous behavior in terms of LREE/HREE. Using this representation, we highlight the fractionation of REEs in each sample compared to the general trend of each element along the fault scarp.

We investigated the northernmost portion of the Caggiano-Montemurro fault, a highly segmented active fault believed to be responsible, or partially responsible, for the 1857 CE (M_w 7.1) and 1561 CE (M_w 6.7) earthquakes (Bello et al., 2022b; Galli, 2020; Galli et al., 2006, 2008; Spina et al., 2008). Three trenches opened by Galli et al. (2006) (location in Figs. 2A and 3A) account for multiple surface-faulting events during the past 6.5 ka (cal. yr B.P.) and show that the latest of these events occurred in the past 2 ka (cal. yr B.P.), possibly during/after the deposition of the slope debris of the Little Ice Age (ca. 1400–1800 CE). The paleosol levels that we sampled for radiocarbon dating (see “Radiocarbon Dating” section) at the side of the fault plane shown in Figures 3D, 6, and 7B provided a calibrated age of ca. 23.7 ka (lower level, sample TDV18–02, 0.9 m below ground surface) and ca. 11.0 ka (upper level, sample TDV18–01, 0.7 m below ground surface) (Table 4). Considering these ages, the lower brownish paleosol was truncated upward by slope erosion at the onset of the LGM, whereas the upper brown level might have formed at the expense of the Neapolitan Yellow Tuff tephra, dated at 14.5 ± 0.4 ka by Galli et al. (2017).

The preservation of the rock-fault scarp mapped in continuity with an en échelon geometry and overlapping portions of the 5 km section investigated with high-resolution topography (Fig. 4) allowed us to reconstruct the deformation profile of the fault over the area. On average, the scarp height profiles in the study area are best fitted by an asymmetric triangular shape (i.e., major peaks shifted from the center of the profile; for detailed description, see Manighetti et al., 2009), suggesting fault-lateral propagation and along-strike variations in the fault evolution (sensu Manighetti et al., 2009). Considering that the structure has recently been mapped for several tens of kilometers (Bello et al., 2022b), further research could be done to investigate the evolution of the structure with the method of high-resolution topography for scarp analysis. Indeed, the same measurement carried out on other exposed portions of the fault would allow the along-strike evolution to be reconstructed as completely as possible. This would allow hypotheses to be formulated about the average propagation direction of long-term ruptures as a whole, that is, those that most frequently characterize the earthquakes’ entire structure. The propagation direction of earthquakes in this area, or at least those that have been associated with the Caggiano-Montemurro fault (i.e., 1561 CE and 1857 CE), is still a matter of discussion (Bello et al., 2022b). In the studied portion, the trend most frequently characterizing the earthquakes is that of propagation from north to south, as suggested by the two asymmetrical envelope lines reconstructed from our VS data. In general, this characteristic of slip profiles is considered as self-similar, indicating areas of maximum and minimum asperity, and a function of fault maturity (Manighetti et al., 2007, 2009). However, we are aware that since the portion of the fault investigated in this work was only 5 km long, our observations could represent a local trend and could fall into a larger-scale trend in a different way. In addition, as also suggested by Bello et al. (2022b), the remarkable segmentation of this fault, and its deviation from the general trends in the bordering large basins aligned instead in “transridge” alignment and generating small, suspended basins within the ridges, could be an index of lower maturity with respect to most of the known faults in the Apennines (e.g., Galli and Peronace, 2014). For all these reasons, subsequent studies are required.

The results obtained from the analysis of the 39 cores acquired here (Figs. 1B, 2A, and 7) on the slickenside cropping out along the Timpa del Vento segment allowed us to obtain the concentrations of REE-Y contents in their up-scarp variation. Our data showed an enrichment and depletion in the up-scarp evolution of the REE-Y concentrations, with clear peaks along the scarp. The analyses also showed that the present soil is enriching the fault plane in contact with the soil (see Fig. 10). We also observed fractionation between LREE/HREE contents, which indicates that the greatest enrichment is due to the behavior of LREEs, while HREEs remain almost immobile during the exchange process at the soil-fault plane. In Figure 11, we show these concentrations along the fault plane, and, following the approach adopted by Carcaillet et al. (2008), Manighetti et al. (2010), and Mouslopoulou et al. (2011), we identified sections that could be considered representative of seismic events or abrupt erosion phases that have occurred over time at least since the LGM. In particular, we considered at least two different scenarios, one with seven higher-magnitude earthquakes, and another with 11 lower-magnitude earthquakes. We underline that since both scenarios have in common the total displacement of the scarp, they should be considered equivalent in terms of total energy release; moreover, they just represent the extremes of what might have happened, whereas mixtures of the two are not excluded. The substantial difference is that, in the first case, we only considered the major enrichments as representing phases of contact between the soil and the fault plane, while in the second scenario, we highlighted each enrichment peak. Each earthquake is represented by the portion of the scarp that is enriched (gray areas in Fig. 11) plus the nonenriched portion (white area immediately before the gray one in Fig. 11) that has not been in contact with the soil due to the minimum and variable thickness of the soil itself. In fact, as shown in the conceptual model in Figure 5, each enriched fault plane portion should not be considered with respect to all the others, but only with respect to an immediately preceding non-enriched portion. The sum of them represents the thickness of the earthquake displacement. Although the methodology applied in this work might have potential perspective for the study of past earthquakes, the findings are subject to at least two assumptions that need to be acknowledged: (1) Each of the peaks obtained is coseismic and not due to slope erosion; (2) the earthquakes that can be counted represent the minimum number of earthquakes sourced by the fault, as highlighted by Mouslopoulou et al. (2011). This latter statement is based in turn on the following three reasons: (2.1) Low- to moderate-magnitude earthquakes are not recorded because they did not rupture at the surface; (2.2) some earthquakes could have occurred temporally so close to each other as to not allow the soil to enrich the fault plane in REEs; and (2.3) some earthquakes may have had slip values less than the thickness of the soil. As far as assumption 1, as the fault activity and capability have been demonstrated by paleoseismological investigations (Galli et al., 2006), we feel confident that peaks could really be coseismic.

We speculated on the magnitude of the earthquakes that could have generated the exhumations of the fault plane as seen in Figure 11, considering different empirical relationships, which needed the maximum and average displacement as input data (i.e., Leonard, 2010; Thingbaijam et al., 2017; Wells and Coppersmith, 1994; Wesnousky, 2008). As shown in Table 5, we calculated the magnitudes for both scenarios in Figure 11, considering the VS obtained from the REE analysis both as average displacement and as maximum displacement (i.e., both for the empirical laws of Wells and Coppersmith [1994] and Wesnousky [2008] and only the average displacement for Leonard [2010] and Thingbaijam et al. [2017]). Obviously, this is an assumption, since, depending on the earthquake that generated it, the coseismic rupture that occurred at the measurement site may not represent the maximum or the average of the deformation. The result of the magnitudes can, however, be considered as a minimum, and therefore worthy of note. According to empirical laws, the observed displacements would account for earthquakes of 6.1 ≤ M_w ≤ 7.0 for scenario 1 and for earthquakes of 5.5 ≤ M_w ≤ 6.7 for scenario 2 (excluding Eq. 7 and Eq. 11 from scenario 1 and 2, respectively, which has a minimum VS of 20 cm but for which the maximum is unknown). Assuming an average of all the empirical laws for each earthquake, we estimate that the magnitude was always 6.1 ≤ M_w ≤ 6.7 for scenario 1 and 6.1 ≤ M_w ≤ 6.5 for scenario 2. Finally, regression curves offer an estimate of the most recurrent earthquake of M_w = 6.6 in scenario 1 and M_w = 6.4 in the case of scenario 2.

From the REE analysis, it is not possible to assign a time to the exhumation/faulting events. However, it is possible to constrain the number of supposed earthquakes within the post-LGM period, locally suggested by the radio-carbon age of the truncated paleosol in the hanging wall (ca. 24 cal. k.y. B.P.; Table 4). Therefore, on one hand, an earthquake with an average M_w = 6.6 would occur every 3.4 to 2.2 k.y. in scenario 1, while, on the other hand, an earthquake with an average M_w = 6.4 would occur every ~2.2–1.5 k.y. Certainly, these calculations derive from the assumption of the maximum of the considered time (~24 k.y.), while it is likely that the observation time scale should be further restricted (less than 16 k.y., considering climatic flattening of the slope, which, according to Galli et al. [2006], ended ca. 16 ka, and the time needed for the soil to develop after the LGM). Nonetheless, many authors in the literature agree that faults may cluster earthquakes in frequent short periods, which then alternate with long periods of nonseismicity (e.g., Bull et al., 2006; Cowie et al., 2012, 2017; Galli et al., 2008; Mildon et al., 2022; Mouslopoulou et al., 2009; Romano et al., 2013; Wedmore et al., 2017; Zinke et al., 2018; among others). This would explain the supposed activation of the Caggiano-Montemurro fault in 1561 and again in 1857. A recurrence time between ~1.6 and ~2.5 k.y. is in the same order as the other known active faults of the Apenninic chain (Galli, 2000) as well as a slip rate, which, according to our results, ranges between ~0.6 and ~0.9 mm/yr. These values were obtained by measuring ~14 m (average) of cumulative offset across the fault (see previous sections and Fig. S1) that occurred after the climatic flattening of the slope (ended ca. 16 ka). Our results confirm previous findings, contribute to our understanding of faults and earthquakes from this portion of Italy, and provide an approach that should be useful worldwide.

We are grateful for the reviews by Editor Ken McCaffrey and an anonymous referee, which substantially improved the manuscript. This research was supported by PRIN 2017 (2017KT2MKE) funds from the Italian Ministry of University and Research for the project “Overtime Tectonic, Dynamic and Rheologic Control on Destructive Multiple Seismic Events—Special Italian Faults and Earthquakes: From Real 4D Cases to Models” (principal investigator Giusy Lavecchia) and by funds provided to Francesco Stoppa and Francesco Brozzetti. The analyses were carried out at the “uDa Analytical high-Tech laborAtory” (DATA) and at the laboratory of newborn screening, proteomics, and endocrinology of the Center for Advanced Studies and Technology (CAST), and at the Laboratory of Structural Geology Cartography and Geological Modeling, all at the G. d’Annunzio University of Chieti–Pescara.