Regional-scale paleofluid system across the Tuscan Nappe–Umbria–Marche Apennine Ridge (northern Apennines) as revealed by mesostructural and isotopic analyses of stylolite–vein networks

We report the results of a multiproxy study that combines structural analysis of a fracture–stylolite network and isotopic characterization of calcite vein cements and/or fault coating. Together with new paleopiezometric and radiometric constraints on burial evolution and deformation timing, these results provide a first-order picture of the regional fluid systems and pathways that were present during the main stages of contraction in the Tuscan Nappe and Umbria–Marche Apennine Ridge (northern Apennines). We reconstruct four steps of deformation at the scale of the belt: burial-related stylolitization, Apenninic-related layer-parallel shortening with a contraction trending NE–SW, local extension related to folding, and late-stage fold tightening under a contraction still striking NE–SW. We combine the paleopiezometric inversion of the roughness of sedimentary stylolites – that constrains the range of burial depth of strata prior to layer-parallel shortening – with burial models and U–Pb absolute dating of fault coatings in order to determine the timing of development of mesostructures. In the western part of the ridge, layer-parallel shortening started in Langhian time (∼15 Ma), and then folding started at Tortonian time (∼8 Ma); late-stage fold tightening started by the early Pliocene (∼5 Ma) and likely lasted until recent/modern extension occurred (∼3 Ma onward). The textural and geochemical (δ¹⁸O, δ¹³C, Δ₄₇CO₂ and ⁸⁷Sr∕⁸⁶Sr) study of calcite vein cements and fault coatings reveals that most of the fluids involved in the belt during deformation either are local or flowed laterally from the same reservoir. However, the western edge of the ridge recorded pulses of eastward migration of hydrothermal fluids (>140 ^∘C), driven by the tectonic contraction and by the difference in structural style of the subsurface between the eastern Tuscan Nappe and the Umbria–Marche Apennine Ridge.

1 Introduction

The upper crust is the locus of omnipresent fluid migrations that occur at all scales, leading to strain localization, earthquake triggering, and georesource generation, distribution and storage (e.g. Cartwright, 2007; Andresen, 2012; Bjørlykke, 1994, 1993; Lacombe and Rolland, 2016; Lacombe et al., 2014; Roure et al., 2005; Agosta et al., 2016). Carbonate rocks host an important part of the world's exploited hydrocarbons, strategic ores, and water resources (Agosta et al., 2010). It is thus a fundamental topic to depict the history of fluid migration in deformed carbonates. Such knowledge no only impacts both the prediction and monitoring of energy prospects and potential storage areas but also may help refine our understanding of what mechanisms facilitate fluid migrations during diagenesis of the sedimentary rocks, along with temporal and spatial scales of fluid flow.

Fluid migration and accumulation events are strongly controlled by tectonics, especially by the related development of large-scale faults and fracture networks. Indeed, structural studies established that fracture networks in folded reservoirs are not exclusively related to the local folding history (Stearns and Friedman, 1972) and can also be influenced by burial history (Becker et al., 2010; Laubach et al., 2010, 2019) and long-term regional deformation (Lacombe et al., 2011; Quintà and Tavani, 2012; Tavani and Cifelli, 2010; Tavani et al., 2015; Bellahsen et al., 2006; Bergbauer and Pollard, 2004; Ahmadhadi et al., 2008; Sassi et al., 2012; Beaudoin et al., 2012; Amrouch et al., 2010). In fold-and-thrust belts and orogenic forelands, the complex deformation history can be unravelled by studying the development of mesoscale structures such as faults, veins, and stylolites (Tavani et al., 2015). Several deformational stages often affect the strata in such settings, starting with extension related to foreland flexure and bulging layer-parallel shortening (so-called pre-folding if kinematically unrelated with folding, early-folding otherwise). Once contraction is accommodated by strata tilting syn-folding, strata curvature-related structures develop under local extension. After the tilting is over, continuous contraction leads to the late-stage fold tightening. Mesoscale structures developing during the so-called post-folding events are kinematically unrelated with folding.

In the past decades, a significant volume of work has been conducted in order to reconstruct past fluid migrations through either localized fault systems or distributed sub-seismic fracture networks from the scale of a single fold to that of the basin itself (Engelder, 1984; Reynolds and Lister, 1987; McCaig, 1988; Evans et al., 2010; Evans and Hobbs, 2003; Evans and Fischer, 2012; Forster and Evans, 1991; Cruset et al., 2018; Lacroix et al., 2011; Travé et al., 2000, 2007; Bjørlykke, 2010; Callot et al., 2017a, b; Roure et al., 2005, 2010; Van Geet et al., 2002; Vandeginste et al., 2012; Vilasi et al., 2009; Barbier et al., 2012; Beaudoin et al., 2011, 2013, 2014, 2015; Fischer et al., 2009; Lefticariu et al., 2005; Di Naccio et al., 2005). In some cases, fluid migration is stratigraphically compartmentalized and directed by compressive tectonic stress. In other cases, mineralization records an infill of meteoric fluids flowing downward or of hydrothermal fluids (i.e. hotter than the host rock they precipitated in) flowing upward, either from the basin or the basement rocks through large-scale faults or décollement (Roure et al., 2005; Vandeginste et al., 2012; Cruset et al., 2018; Lacroix et al., 2011; Travé et al., 2000, 2007; De Graaf et al., 2019; Callot et al., 2010; Beaudoin et al., 2011, 2014; Bertotti et al., 2017; Gonzalez et al., 2013; Lucca et al., 2019; Mozafari et al., 2019; Storti et al., 2018; Vannucchi et al., 2010).

This contribution reports an orogen-scale paleofluid flow study in the northern Apennines (Italy). The study builds upon the mesostructural and geochemical analysis of vein and stylolite networks within the competent Jurassic–Oligocene carbonate platform along a transect running across the Tuscan Nappe and the Umbria–Marche Apennine Ridge (UMAR) (Fig. 1a). The data collection was organized to cover a large area comprising several folds in order to be able to differentiate regional trends from local, fold-related ones. We focused on identifying and characterizing first-order pattern of mesostructures – faults, fractures, and stylolites – associated with layer-parallel shortening and with thrust-related folding, along with the isotopic measurements (δ¹⁸O, δ¹³C, ⁸⁷Sr∕⁸⁶Sr), clumped-isotope measurements (Δ₄₇CO₂), and U–Pb absolute dating of their calcite cements. Without an appraisal of which fracture trends are relevant to the large-scale (i.e. regional) tectonic evolution, there is a risk to otherwise capture mesostructural and geochemical signals of local meaning only. In order to discuss the local versus hydrothermal fluid origin, we also reconstructed burial curves using strata thickness we correct from physical and chemical compaction. Novel constraints are added to the timing and minimal depth of layer-parallel shortening-related deformation based on the study of the roughness of bedding-parallel stylolites, the inversion of which reliably returns the maximum depth at which compaction under a vertical maximum principal stress was prevailing in the strata. U–Pb absolute dating of calcite steps on mesoscale faults further constrains the timing of folding. Such a multiproxy approach – one that combines structural analysis of fracture–stylolite networks and isotopic characterization of cements/coatings, together with new constraints on burial evolution and deformation timing – provides a picture of the regional fluid pathways during the Apenninic contraction.

Figure 1 (A) Simplified geological map of the study area, with locations of the sampling and measurement sites. Frames relate to the fracture study areas used in Fig. 4. Exact location of measurement sites is reported as black and red points, and labelled black points also represent the sampling sites for geochemical analysis. (B) Stratigraphic column based on stratigraphic and well data from the central part of the UMAR, after Centamore et al. (1979). (C) Crustal-scale composite cross section based on published seismic data interpretations; A-A' modified after Carboni et al. (2020); B-B' and C-C' after Scisciani et al. (2014). Note that both tectonic style (thick skinned and thin skinned) are represented by question marks for the UMAR.

2 Geological setting

The Neogene-to-Quaternary Apennines fold-and-thrust belt formed during the convergence of Eurasia and Africa (Lavecchia, 1988; Elter et al., 2012). It is associated with the eastward retreating subduction of the Adriatic Plate under the European Plate. The Apennines extend from the Po Plain to the Calabrian Arc and are divided into two main arcs, the northern Apennines that extend down to the south of the UMAR and the Southern Apennines that cover the remaining area down to the Calabrian Arc (Carminati et al., 2010). The evolution of the Apennines is characterized by a roughly eastward migration of thrust fronts and associated foredeep basins, superimposed by post-orogenic extension at the rear of the eastward propagating orogenic belt (Cello et al., 1997; Tavani et al., 2012; Lavecchia, 1988; Ghisetti and Vezzani, 2002).

The study area, the Tuscan Nappe and the Umbria–Marche Apennines Ridge, comprises a succession of carbonate rocks, Upper Triassic to Oligocene in age, which corresponds to a carbonate platform (Lavecchia, 1988; Carminati et al., 2010). The Umbrian carbonate units overlie early Triassic evaporites that act as a décollement level and that are themselves unconformably overlying the crystalline basement rocks (Fig. 1b). Above the platform, Miocene turbidite deposits record the progressive eastward involvement of the platform into the fold-and-thrust belt (Calamita et al., 1994). In the western part of the area, the belt is a thin-skinned assembly of piggy-back duplex folds (Fig. 1c), the so-called Tuscan Nappe, the folding and thrusting of which started by the late Aquitanian (ca. 23–20.43 Ma) and lasted until the Langhian (ca. 16 to 13.8 Ma; Carboni et al., 2020). The UMAR is an arcuate ridge exhibiting an eastward convex shape, with a line connecting Perugia and Ancona separating a northern part where structural trends are oriented NW–SE, from a southern part where structure trends are oriented N–S (Calamita and Deiana, 1988). Burial models suggest that from Burdigalian to early Messinian times, the Tuscan Nappe was further buried under the allochthonous Ligurian thrust sheet, reaching locally up to 1 km in thickness (Caricchi et al., 2015). In the eastern part (now UMAR), the foreland was progressively folded and thrusted from the Lower Miocene in the westernmost part of the current ridge to the Messinian in the foreland of the ridge (Mazzoli et al., 2002). UMAR was considered a thin-skinned thrust belt where shortening was accommodated by stacking and duplexing of sedimentary units detached above a décollement level located in the Triassic evaporites (Conti and Gelmini, 1994; Carboni et al., 2020). The seismic profile of the CROsta Profonda (CROP) project led authors to interpret the UMAR as resulting from thick-skinned tectonics, where the basement is involved in shortening (Barchi et al., 1998) through the positive inversion of normal faults inherited from the Jurassic Tethyan rifting (Fig. 1c). Even if the interpretation of basement-involved shortening is more accepted now, the subsurface geometry is still debated, with some models involving shallow duplexes (Tavarnelli et al., 2004; Mirabella et al., 2008), while in more recent works surface folds are rather interpreted as related to high-angle thrusts that either sole within the mid-Triassic décollement or involve the basement (Scisciani et al., 2014, 2019; Butler et al., 2004) (Fig. 1c). For the latter interpretations, the style of deformation of the UMAR strongly contrasts with the style of deformation of the Tuscan Nappe where shortening is accommodated by allochthonous, far-travelled duplex nappes (Carboni et al., 2020) (Fig. 1c). The cross section in Fig. 1c also implies that at least part of the motion on the décollement level at the base of the Tuscan Nappe postdates the westernmost activation of steep thrusts of the UMAR, as the thrust at the base of the Tuscan Nappe cuts and offsets the west-verging basement fault in the area of Monte Subasio. Currently, the whole Tuscan Nappe–UMAR area is experiencing extension, with numerous active normal faults developing trenches, as the contraction front migrated toward the Adriatic Sea (d'Agostino et al., 2001).

Our sampling focused on the carbonate formations cropping out from W to E in the Cetona area located west from Perugia in the Tuscan Nappe – Monte Corona, Monte Subasio, the Gubbio area, Spoletto Area, Monte Nero, Monte Catria, Monte San Vicino, and Monte Cingoli in the UMAR, as well as Monte Conero, the youngest onshore anticline related to the Apenninic compression, located on the coast line (Fig. 1a). The sampled units comprise, following the stratigraphic order (Fig. 1b): the Triassic anhydrites and dolostones of the Anidridi di Burano Formation focusing on the limestone and marl intercalation at the top; (2) Liassic massive dolomites of the Calcare Massiccio Fm. (Hettangian to Sinemurian, ca. 203–195 Ma); (3) the grey Jurassic limestones with chert beds of the Corniola Fm. (Sinemurian–Pleisbachian, ca. 195–184 Ma); (4) the micritic limestones, marls, and cherts of the Calcare Diasprini Fm (Bajocian–Tithonian, ca. 160–140 Ma); (5) the white limestones with chert beds of the Maiolica Fm. (Tithonian–Aptian, ca. 140–110 Ma); (6) the marly limestones of the Fucoidi Fm. (Aptian–Albian, ca. 110–100 Ma); (7) the white marly limestones of the Scaglia Bianca Fm. (Albian–Turonian, ca. 100–90 Ma); (8) the pink marly limestones of the Scaglia Rossa Fm. (Turonian–Priabonian, ca. 90–37 Ma); and (9) the grey marly limestones of the Scaglia Cinerea Fm. (Priabonian–Aquitanian, ca. 37–21 Ma). Up to 3000 m of Miocene turbidites were deposited when the area of interest was the foredeep ahead of the advancing fold-and-thrust belt and during fold development, including clay-rich limestones and silts of Marnoso–Aranacea (Aquitanian–Tortonian, ca. 21–7 Ma); in the eastern part of the ridge (east of the Cingoli anticline), thicker foredeep deposits are Messinian to Pliocene (ca. 7 to 3 Ma) in age.

3 Methods

The approach adopted relies on structural and geochemical analyses, the results of which are combined to establish a scenario of fluid flow and fluid–rock interaction during deformation of the Umbria–Marche Apennine Ridge in the northern Apennines.

3.1 Mesostructural analysis of joints, veins, and striated fault planes

Approximately ∼1300 joint and vein orientations, along with tectonic stylolite orientations, were measured along a WSW–ENE transect going from Cetona in the Tuscan Nappe to Monte Conero on the coastline (Fig. 1a). For each measurement site, fractures and tectonic stylolites were measured. Tectonic stylolites considered are either bedding-perpendicular dissolution planes displaying horizontal peaks after unfolding or vertical dissolution planes displaying horizontal peaks in the current bed attitude. Crosscutting and abutment relationships between mesostructures were carefully observed in the field (Fig. 2) and checked in thin sections under the optical microscope when possible (Fig. 3). Poles to fractures and stylolite peaks were projected on Schmidt stereograms, lower hemisphere, in the current attitude of the strata (raw) and after unfolding (unfolded) (Fig. 4). Assuming the same mode of deformation (i.e. mode I opening) and consistent chronological relationships and orientation, we use the pole density obtained from Fisher statistical analysis to define the sets of joints and/or veins which are the most documented and representative at fold scale. Tectonic stylolite planes and peaks were measured, and we consider that the average orientation of the stylolite peaks at the fold scale represents the orientation of the horizontal maximum principal stress (σ₁), as peaks grow parallel to the main shortening direction (Koehn et al., 2007). To complement this mesostructural analysis, striated fault planes were measured (1) in the Langhian carbonates from the Camerino syncline west from Monte San Vicino and (2) in the forelimb of Monte Subasio, with one site in the Scaglia Cinerea and one site in the Scaglia Rossa e Bianca. At each site, paleostress orientations (local trend and plunge) and regimes (reverse, extensional, strike-slip) were calculated using inversion techniques (Angelier, 1984) discussed in Lacombe (2012). Published studies in the UMAR highlight the complexity of fracture patterns at the fold scale that record several phases of stress perturbation and stress/block rotation due to the local tectonics and structural inheritance (Tavani et al., 2008; Petracchini et al., 2012; Beaudoin et al., 2016; Díaz General et al., 2015). In order to capture the mesostructural and fluid flow evolution at the regional scale during layer-parallel shortening and folding, we gathered the most represented fractures by structure, regardless of the structural complexity in the individual folds and corrected them from the local bedding dip using an open-source stereodiagram rotation program (Grohmann and Campanha, 2010) to discriminate between early, syn-, and late folding features.

Figure 2 Field photographs showing chronological relationships between veins, joints, and stylolites. (a) Monte Nero and (b) Monte Cingoli sets are reported along with local chronological order between brackets, as deduced from crosscutting and abutment relationships displayed. Dotted lines depict the labelled fracture and tectonic stylolite orientation in black and blue, respectively. Bedding-parallel stylolites spacing from (c) the Maiolica Fm. (Monte Subasio) and (d) the Massiccio Fm. (Monte Cingoli). Intersections between bedding-parallel stylolites and scanline are pointed out with orange arrows.

Figure 3 (a, b) Photomicrographs of various veins in plain polarized light, (c, d) with corresponding view under cathodoluminescence; top one is a set J1 vein from the Scaglia Fm., while the bottom one is a set J2 from the Maiolica Fm.

3.2 Inversion of sedimentary stylolites

Bedding-parallel stylolites are rough dissolution surfaces that developed in carbonates in flat laying strata during burial at the time when σ₁ was vertical. As proposed by Schmittbuhl et al. (2004) and later developed by Koehn et al. (2012), Ebner et al. (2009b, 2010), Rolland et al. (2014) and Beaudoin et al. (2019, 2020), the 1D roughness of a track along the bedding-parallel stylolite (i.e. difference in height between two points along the track) results from a competition between roughening forces (i.e. pining on non-soluble particles in the rocks) and smoothing forces (i.e. the surface energy at scale typically <1 mm and the elastic energy at scale >1 mm). The stylolite growth model (Koehn et al., 2007; Ebner et al., 2009a; Rolland et al., 2012; Toussaint et al., 2018) predicts that the surface energy-controlled scale returns a steep slope characterized by a roughness exponent (so-called Hurst exponent) of 1.1±0.1, while the elastic energy-controlled scale returns a gentle slope with a roughness exponent of 0.5 to 0.6 (Fig. 5). The length at which the change in roughness exponent occurs, called the crossover length (Lc, in mm), is directly related to the magnitude of differential and mean stress (${\mathit{\sigma }}_{\mathrm{d}}={\mathit{\sigma }}_{\mathrm{1}}-{\mathit{\sigma }}_{\mathrm{3}}$ and ${\mathit{\sigma }}_{\mathrm{m}}=\frac{{\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{3}}}{\mathrm{3}}$, respectively, in Pa) prevailing in the strata at the time the stylolite stopped to be an active dissolution surface following

$$\begin{array}{}\text{(1)}& \mathrm{Lc}={\displaystyle \frac{\mathit{\gamma }E}{\mathit{\beta }{\mathit{\sigma }}_{\mathrm{m}}{\mathit{\sigma }}_{\mathrm{d}}}},\end{array}$$

where E is the Young modulus of the rock (in Pa), γ is the solid-fluid interfacial energy (in J m⁻²), and β=ν(1-2ν)∕π, a dimensionless constant with ν being the Poisson ratio. Samples of bedding-parallel stylolites of which peaks were perpendicular to the dissolution plane were collected in specific points of the study area, and several stylolites were inverted. The inversion process follows the method described in Ebner et al. (2009b). Samples were cut perpendicular to the stylolite, hand polished to enhance the visibility of the track while being cautious to not altering the peaks, and scanned at high resolution (12 800 pixel per inches), and the 1D track was hand drawn with a pixel-based software (GIMP). Each track was analysed as a periodic signal by using the average wavelet spectrum with Daubechies D4 wavelets (Fig. 5) (Ebner et al., 2009b; Simonsen et al., 1998). In the case of bedding-parallel stylolites related to compaction and burial, the method assumes isotropic horizontal stress (σ_v>>σ_h=σ_H). This assumption can be checked by considering different tracks of the same sample, as horizontal stress isotropy during pressure solution implies a constant value of the crossover length regardless of the track direction. This leads to simplify Eq. (1) (Schmittbuhl et al., 2004) as

$$\begin{array}{}\text{(2)}& {\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{2}}={\displaystyle \frac{\mathit{\gamma }E}{\mathit{\alpha }\mathrm{Lc}}},\end{array}$$

with

$$\begin{array}{}\text{(3)}& \mathit{\alpha }={\displaystyle \frac{\left(\mathrm{1}-\mathrm{2}\mathit{\nu }\right)\cdot (\mathrm{1}+\mathit{\nu }{)}^{\mathrm{2}}}{\mathrm{30}\mathit{\pi }(\mathrm{1}-\mathit{\nu }{)}^{\mathrm{2}}}}.\end{array}$$

According to the sampled formation, we used the solid-fluid interfacial energy γ of 0.24 J m⁻² for dolomite and of 0.32 J m⁻² for calcite (Wright et al., 2001). As an approximation for the material mechanical properties, we use a Poisson ratio of $\mathit{\nu }=\mathrm{0.25}\pm \mathrm{0.05}$ and the average Young modulus derived from the Jurassic–Eocene competent core of E=24.2 GPa (Beaudoin et al., 2016). It is important to note that because of the nonlinear regression method we use – and because of uncertainty on the mechanical parameters of the rock at the time it dissolved – the uncertainty on the stress has been calculated to be about 12 % (Rolland et al., 2014). As the dissolution occurs along a fluidic film (Koehn et al., 2012; Rolland et al., 2012; Toussaint et al., 2018), the stylolite roughness is unaffected by local fluid overpressure until the system is fluidized and hydro-fractures (Vass et al., 2014), meaning it is possible to translate vertical stress magnitude directly into depth if considering an average dry rock density for clastic and carbonate sediments (2400 g m⁻³; Manger, 1963), without any additional assumption on the past thermal gradient or fluid pressure (Beaudoin and Lacombe, 2018). This technique has already provided meaningful results in various settings (Bertotti et al., 2017; Rolland et al., 2014; Beaudoin et al., 2019, 2020).

Figure 4 Poles of measured joints and veins and stylolite peaks projected on Schmidt stereograms, lower hemisphere, for the different structures. Data are projected in the current attitude of the strata (left, R), with pole to bedding in red, and after unfolding (right, U). Red areas represent highest density according to Fischer statistical analysis using the software OpenStereo, and main fracture set average orientations are represented as blue planes. For tectonic stylolite peaks, the blue square represents highest density according to Fischer statistical analysis. Striated fault inversion results are reported in the current attitude of the strata (bedding as dashed line).

Figure 5 Examples of results of stylolite roughness inversion, with signal analysis by Average Wavelet in Monte Nero (a) and in Monte Subasio (b).

Sedimentary stylolites also yield quantitative information on the volume of dissolved rocks during burial (Toussaint et al., 2018), the minimum of which can be approached in 1D by measuring the amplitude of the highest peak along the stylolite track (Table 1) and multiplying this height value (in m) by the average density of stylolites (number per metre) derived from field spacing measurement (Fig. 2c–d; in m).

Table 1 Results of stylolite roughness inversion applied on bedding-parallel stylolites.

^* Crossover length given within 23 % uncertainty. Vertical stress σv given within 12 % uncertainty using for Young modulus E=23.2 GPa (Beaudoin et al., 2014), Poisson ratio ν=0.25, and interfacial energy γ=0.32 J m⁻². Depth calculated using dry density of rock d=2400 g m⁻³, acceleration of the gravitational field g=9.81 m s⁻².

3.3 O and C stable isotopes

Calcite cements in tectonic veins related either to layer-parallel shortening or to strata curvature at fold hinges were studied petrographically (Fig. 3). The vein textures were characterized in thin sections under an optical microscope, and possible post-cementation diagenesis such as dissolution or replacement were checked under cathodoluminescence, using a cathodoluminescence CITL CCL 8200 Mk4 operating under constant gun condition of 15 kV and 300 µA. To perform oxygen and carbon stable isotope analysis on the cements that were the most likely to record the conditions of fluid precipitation at the time the veins opened, we selected those veins that (1) show no obvious evidence of shear; (2) the texture of which was elongated blocky or fibrous (Fig. 3; Bons et al., 2012); and (3) show homogeneous cement under cathodoluminescence (Fig. 3), precluding any later diagenetic alteration.

To begin with, 0.5 mg of calcite powder was manually collected for each of 58 veins and 54 corresponding host rocks (sampled ∼2 cm away from veins) in various structures and formations along the transect, in both Tuscan Nappe and UMAR. Carbon and oxygen stable isotopes were analysed on an Analytical Precision AP2003 mass spectrometer equipped with a separate acid injector system. Samples (calcite or dolomite) were placed in glass vials to conduct a reaction with 105 % H₃PO₄ under a helium atmosphere at 90 ^∘C, overnight. Results are reported in Table 2, in permil relative to Vienna PeeDee Belemnite (‰ VPDB). Mean analytical reproducibility based on replicates of the SUERC laboratory standard MAB-2 (Carrara marble) was around ±0.2 ‰ for both carbon and oxygen. MAB-2 is an internal standard extracted from the same Carrara marble quarry, as is the IAEA-CO208 1 international standard. It is calibrated against IAEA-CO-1 and NBS-19 and has exactly has exactly the same C and O isotope values as IAEA-CO-1 (−2.5 % and 2.4 % VPDB, respectively).

Table 2 Results of stable isotopic analyses of O, C, and ⁸⁷Sr∕^86Sr.

^* Values were corrected to reflect the fact that host rocks is dolomite; HR stands for host rock, V stands for vein.

3.4 ⁸⁷Sr∕⁸⁶Sr measurements

Sr-isotope analysis was performed at the Geochronology and Tracers Facility, British Geological Survey. First, 2–10 mg of sample was weighed into 15 mL Savillex teflon beakers and dissolved in 1–2 mL of 10 % Romil uPA acetic acid. After evaporating to dryness, the samples were converted to chloride form using 2 mL of Teflon-distilled HCl. The samples were then dissolved in ca. 1 mL of calibrated 2.5 M HCl in preparation for column chemistry and centrifuged. Samples were pipetted onto quartz–glass columns containing 4 mL of AG50x8 cation exchange resin. Matrix elements were washed off the column using 48 mL of calibrated 2.5 M HCl and discarded. Sr was collected in 12 mL of 2.5 M HCl and evaporated to dryness.

Sr fractions were loaded onto outgassed single Re filaments using a TaO activator solution and analysed in a Thermo-Electron Triton mass spectrometer in multi-dynamic mode. Data are normalized to ⁸⁶Sr∕⁸⁸Sr =0.1194. Three analyses of the NBS987 standard run with the samples gave a value of 0.710250±0.000001 (1SD).

3.5 Carbonate clumped-isotope paleothermometry (Δ₄₇CO₂)

Clumped-isotope analyses were carried out in the Qatar Stable Isotope Laboratory at Imperial College London. The technique relies on the tendency for heavy isotopes (¹³C, ¹⁸O) to “clump” together in the same carbonate molecule, which varies only by temperature. Since the clumping of heavy isotopes within a molecule is a purely stochastic process at high temperature but is systematically over-represented (relative to randomly distributing isotopes among molecules) at low temperature, the “absolute” temperature of carbonate precipitation can be constrained using clumped-isotope abundances.

The clumped-isotopes laboratory methods at Imperial College follow the protocol of Dale et al. (2014) as adapted for the automated clumped-isotope measurement system IBEX (Imperial Batch EXtraction) system (Cruset et al., 2016). Typical sample size was 3.5 mg of calcite powder per replicate. Measurement of ¹³C−¹⁸O ordering in sample calcite was achieved by measurement of the relative abundance of the ¹³C¹⁸O¹⁶O isotopologues (mass 47) in acid-evolved CO₂ and is referred to in this paper as Δ₄₇CO₂. A single run on the IBEX comprises 40 analyses, 30 % of which are standards. Each analysis takes about 2 h. The process starts with 10 min of reaction of the carbonate powder in a common acid bath containing 105 % orthophosphoric acid at 90 ^∘C to liberate CO₂. The CO₂ gas is then captured in a water∕CO₂ trap maintained at liquid nitrogen temperature and then moved through a hydrocarbon trap filled with Porapak and a second water trap using helium as carrier gas. At the end of the cleaning process, the gas is transferred into a cold finger attached to the mass spectrometer and into the bellows of the mass spectrometer. Following transfer, analyte CO₂ was measured on a dual-inlet Thermo MAT 253 mass spectrometer (MS “Pinta”). The reference gas used is a high-purity CO₂, with the following reference values: −37.07 ‰ δ¹³C_VPDB and 8.9 ‰ δ¹⁸O_VSMOW. Measurements comprise eight acquisitions each with seven cycles with 26 s integration time. A typical acquisition time is 20 min, corresponding to a total analysis time of 2 s.

Data processing was carried out in the freely available stable isotope management software, Easotope (https://www.easotope.org, last access: 19 February 2020, John and Bowen, 2016). The raw Δ₄₇ CO₂ is corrected in three steps. First, mass spectrometer nonlinearity was corrected by applying a “pressure baseline correction” (Bernasconi et al., 2013). Next, the Δ₄₇ results were projected in the absolute reference frame or carbon dioxide equilibrated scale (CDES; Dennis et al., 2011) based on routinely measured ETH1, ETH2, ETH3, ETH4 and Carrara marble (ICM) carbonate standards (Meckler et al., 2014; Muller et al., 2017). The last correction to the raw Δ₄₇ was to add an acid correction factor of 0.082 ‰ to obtain a final Δ₄₇CO₂ value (Defliese et al., 2015). Temperatures of precipitation were then estimated using the equation of Davies and John (2019). The bulk isotopic value of δ¹⁸O was corrected for acid digestion at 90 ^∘C by multiplying the value by 1.0081 using the published fractionation factor valid for calcite (Kim et al., 2007). Contamination was monitored by observing the values on mass 48 and 49 from each measurement, using a Δ₄₈ offset value >0.5 ‰ and/or a 49 parameter values >0.3 as a threshold to exclude individual replicates from the analysis (Davies and John, 2019).

3.6 U–Pb absolute dating of veins and faults

The calcite U–Pb geochronology was conducted in two different ways (specific methodology of which is reported in the Supplement):

• LA-ICPMS trace elements and U–Pb isotope mapping were performed at the Geochronology and Tracers Facility, British Geological Survey, UK, on six vein samples. Data were generated using a Nu Instruments Attom single collector inductively coupled plasma mass spectrometer coupled to a NWR193UC laser ablation system fitted with a TV2 cell, following protocol reported previously (Roberts et al., 2017; Roberts and Walker, 2016). Laser parameters were 110 µm spots, ablated at 10 Hz for 30 s with a fluence of 7 J cm⁻². WC1 (Roberts et al., 2017) was used as a primary reference material for Pb∕U ratios and NIST614 for Pb∕Pb ratios; no secondary reference materials were run during the session. Additional constraints on U–Pb composition were calculated from the Pb and U masses measured during the trace element mapping. Baselines were subtracted in Iolite, and Pb∕Pb and Pb∕U ratios were calculated offline in Excel. No normalization was conducted, as the raw ratios are suitable accurate to assess

• LA-ICPMS U–Pb isotope mapping approach was undertaken at the Institut des Sciences Analytiques et de Physico-Chimie pour l'Environnement et les Matériaux (IPREM) Laboratory (Pau, France). All the 29 samples were analysed with a 257 nm femtosecond laser ablation system (Lambda3, Nexeya, Bordeaux, France) coupled to an HR-ICPMS Element XR (ThermoFisher Scientific, Bremen, Germany) fitted with the Jet Interface (Donard et al., 2015). The method is based on the construction of isotopic maps of the elements of interest for dating (U, Pb, Th) from ablation along lines, with ages calculated from the pixel values (Hoareau et al., 2020). The ablation was made in a helium atmosphere (600 mL min⁻¹), and 10 mL min⁻¹ of nitrogen was added to the helium flow before mixing with argon in the ICPMS. Measured wash out time of the ablation cell was ∼500 ms for helium gas. The LA-ICPMS coupling was tuned daily, and the additional Ar carrier gas flow rate, torch position, and power were adjusted so that the U∕Th ratio was close to 1±0.05 when ablating the glass SRM NIST612. Detector cross-calibration and mass bias calibration were checked daily. The laser and HR-ICPMS parameters used for U–Pb dating are detailed in the Supplement.

4 Results

4.1 Mesostructural analysis of joints, veins, and striated fault planes

Based on the average orientation and the angle to the local fold axis, we group veins and joints in two main sets labelled J (Fig. 4):

• Set J1 comprises joints and veins at high angle to bedding, of which direction is perpendicular to the local strike of fold axis, i.e. E–W to NE–SW with respect to the arcuate shape of the folds. The trend of this set J1 evolves as follows: E–W in the westernmost part (Cetona, Subasio), E-W to NE–SW in the central part (Catria, Nero), NE–SW in the eastern part of the chain (San Vicino, Cingoli), and ENE–WSW in the far foreland (Conero).

• Set J2 comprises joints and veins at high angle to bedding that strike parallel to the local trend of the fold hinge, i.e. NW–SE in the ridge to N–S in the outermost part of the belt, where the arcuate shape is more marked.

Note that as set J1 strikes perpendicular to the strike of local strata, it is impossible to infer a pre-, syn-, or post-tilting (post-tilting then called J3 hereinafter) origin for its development. In most cases though, abutment relationships establish a relative chronology with set J1 predating set J2 (Fig. 3). The veins of sets J1 and J2 show twinned calcite grains (Fig. 3) with mostly thin and rectilinear twins (thickness <5 µm; Fig. 3). Another set comprising joints striking N–S and oblique to the direction of the fold axis is documented in Monte Catria. It is also encountered at other locations but can be regarded as a second-order set at the regional scale on a statistical basis. This set could be tentatively related to lithospheric flexure (Mazzoli et al., 2002; Tavani et al., 2012), but as it is the least represented in our data, it will not be considered hereinafter.

Most tectonic stylolites have peaks striking NE–SW (Fig. 4), but they can be split in two sets labelled S based on the orientation of their planes with respect to the local bedding.

• Set S1 comprises bed-perpendicular, vertical stylolite planes containing horizontal peaks in the unfolded attitude of strata;

• Set S2 comprises: vertical stylolite planes containing horizontal peaks in the current attitude of strata (set S2).

S1 and S2 were not always easily distinguished when both occurred at the fold scale because (1) stylolite data were often collected in shallow dipping strata, (2) peaks are not always perfectly perpendicular to the stylolite planes, and (3) the orientation data are scattered with intermediate plunges of the peaks. Another set showing stylolite planes with N–S peaks parallel to bedding is documented at Monte Subasio only, thus will not be considered in the regional sequence of deformation.

Finally, inversion of in-plane striation of mesoscale reverse and strike-slip conjugate fault population reveals (1) a NE–SW contraction in the unfolded attitude of the strata (early folding set F1, flexural-slip related, bedding-parallel reverse faults) and (2) a NE–SW contraction in the current attitude of the strata (late folding set F2, strike-slip conjugate faults and reverse faults).

4.2 Inversion of sedimentary stylolites

The paleopiezometric study of 30 bedding-parallel stylolites returned a range of burial depths, across the UMAR, from W to E, reported in Table 1. To ensure that the stress on the horizontal plane was isotropic, the stylolite inversion technique was applied on two orthogonal tracks for each stylolite. Reported Lc correspond to stylolite where both Lc values are similar within uncertainty. Most data come from the western part of the UMAR: in the Subasio Anticline (n=7), the depth returned by the Scaglia Bianca and the lower part of the Scaglia Rossa fms. ranges from ca. 800±100 to ca. 1450±150 m. In Fiastra area (n=6), the depth returned for the Maiolica Fm. ranges from 800±100 to 1200±150 m. In the Gubbio fault area (n=4), the depth returned for the Jurassic Corniola Fm. ranges from 600±70 to 1450±150 m. In Monte Nero (n=11), the depth data published by Beaudoin et al. (2016) and updated here range from 750±100 to 1350±150 m in the Maiolica. Fewer data come from the western part of the UMAR: in Monte San Vicino (n=2), the depth returned for the Maiolica Fm. ranges from 1000±100 to 1050±100 m. Finally, the depth reconstructed for the lower part of the Scaglia Rossa is 650±70 m in the foreland at Conero Anticline (n=1). The maximum height of peaks along the studied stylolite tracks ranges from 0.6 to 8.5 mm (n=30) with a mean value of 2.6 mm. Spacing values for these stylolites were measured on outcrops (Fig. 2) and range from 1 to 2 cm, averaging the number of stylolite per metre to 70. Considering that dissolution is isotropic along the stylolite plane, the volume of rock loss in relation to the chemical compaction is ∼18 %.

4.3 O and C stable isotopes

At the scale of the study area, most formations cropping out were sampled (Table 2), and oxygen isotopic values of the vein cements and striated fault coatings range from −16.8 ‰ to 3.7 ‰ VPDB, while in the host rocks values range from −5.3 ‰ to 0.4 ‰ VPDB. Carbon isotopic values range from −9.7 ‰ to 2.7 ‰ VPDB and from 0.0 ‰ to 3.5 ‰ VPDB in the veins and in the host rock, respectively (Fig. 6a–b). Isotopic values are represented either according to the structure where they have been sampled in, irrespective of the structural position in the structure (i.e. limbs or hinge), or according to the set they belong to, differentiating the sets J1, J2, F1, and F2 (Fig. 6a). At the scale of the belt, isotopic values of host rocks are very similar, the only noteworthy point being that the Triassic carbonates have lower δ¹⁸O than the rest of the column (δ¹⁸O of −5.5 ‰ to −3.5 ‰ versus −3.2 ‰ to −1.0 ‰ VPDB; Fig. 6b). Considering the vein cements, the δ¹³C values are rather similar in all structures and in all sets, a vast majority of veins showing cements with values of 1.5±1.5 ‰ VPDB. In most structures, the δ¹³C values of the veins are similar to the δ¹³C values of the host rock, with the notable exception of the veins hosted in Triassic carbonates of the Cetona anticline, where the shift between δ¹³C values of veins and δ¹³C values of host rocks ranges from +4.0 ‰ VPDB to +7.5 ‰ VPDB (Fig. 6c). The δ¹⁸O values range from −6.0 ‰ VPDB to +3.7 ‰ VPDB in most of the structures and formation, irrespectively of vein set. However, veins sampled in Monte Subasio and Monte Corona return very negative δ¹⁸O values <−15.0 ‰ VPDB (Fig. 6a). The shift between the δ¹⁸O value of the vein and the δ¹⁸O value of the surrounding host rock (Fig. 6d) increases from the western part of the belt (down to −15.0 ‰ VPDB in Monte Corona and Monte Subasio) to the central and eastern part of the belt (up to +5.0 ‰ VPDB in Monte Nero, Monte San Vicino, and Monte Conero).

Figure 6 (a) Plot of δ¹³C vs. δ¹⁸O (‰ VPDB) of veins represented by structure. Frames represent the different type of fluid system. (b) Plot of δ¹³C vs δ¹⁸O (‰ VPDB) of host rocks represented by structure. (c) Plot of δ¹³C_vein vs. δ¹³C_host (‰ VPDB) of veins represented by structures. (d) Plot of the difference between δ¹⁸O of host rocks and δ¹⁸O of veins (‰ VPDB) vs. eastward distance from the Cetona Anticline towards the Adriatic Basin across the strike of the UMAR. Data are represented by tectonic sets. The proposed extension of the Ligurian nappe overthrust is reported after Caricchi et al. (2014); red frames and arrows represent the fluid systems.

4.4 ⁸⁷Sr∕⁸⁶Sr measurements

Analyses were carried out on seven veins and six corresponding host rocks, distributed over three structures of the UMAR (Monte Subasio, Monte Nero, and Monte San Vicino, from the hinterland to the foreland) and three formations (the Calcare Massiccio, the Maiolica, and the Scaglia fms.; Fig. 7, Table 2). Vein sets sampled are the J1, J2, and J3 sets described in the whole area. ⁸⁷Sr∕⁸⁶Sr values of host rocks differ according to the formations, being the highest in the Scaglia Rossa Fm. (⁸⁷Sr∕⁸⁶Sr ≈0.70780), intermediate in the Calcare Massicio Fm. (⁸⁷Sr∕⁸⁶Sr ≈0.70760), and the lowest in both the Scaglia Bianca and the Maiolica Fm. (⁸⁷Sr∕⁸⁶Sr ≈0.70730). ⁸⁷Sr∕⁸⁶Sr values of host rocks are in line with expected values for seawater at the time of their respective deposition (McArthur et al., 2001). ⁸⁷Sr∕⁸⁶Sr values of veins scatter from 0.70740 to 0.70770, with lower values in Monte Nero and in Monte San Vicino (0.707415 to 0.707646, irrespective of fracture set) and higher values in Monte Subasio (0.707644 to 0.707690, set J2). One vein cement of J3 in Monte Subasio returned a lower ⁸⁷Sr^∕86Sr value of 0.707437.

Figure 7 Plot of ⁸⁷Sr∕⁸⁶Sr values vs. longitude, with filling colour related to tectonic set and point shape related to host rock formation. Note that the host-rock values are equivalent to the related formation seawater values.

4.5 Carbonate clumped-isotope paleothermometry (Δ₄₇CO₂)

Sixteen samples were analysed, including cements of NE–SW (J1) and NW–SE (J2) pre-folding vein sets, along with coatings of early folding reverse (F1) and late folding strike-slip conjugate mesoscale faults (F2). Regardless of the structural position in the individual folds, Δ₄₇CO₂ values for veins (Table 3) range from 0.511±0.004 to 0.608±0.000 ‰ CDES in Monte Corona (n=3), from 0.468±0.032 to 0.574±0.000 ‰ CDES in Monte Subasio (n=5), from 0.662±0.004 to 0.685±0.052 ‰ CDES in Monte Nero (n=2), from 0.568±0.035 to 0.658±0.018 ‰ CDES in Monte San Vicino and the syncline to its west (n=4), from 0.601±0.032 to 0.637±0.013 ‰ CDES in Monte Catria (n=2), and is 0.643±0.044 ‰ CDES in Monte Conero. Corresponding δ¹⁸O and δ¹³C are reported in Table 2. Analysis of Δ₄₇CO₂ returns the precipitation temperature (T), and the oxygen isotopic values of the mineralizing fluid can be calculated using the δ¹⁸O of the mineral, the clumped-isotope temperature and the fractionation equation of Kim and O'Neil (1997) providing the fractionation coefficient as a function of the temperature (Fig. 8). Veins and faults belong to the Calcare Massiccio Fm., the Maiolica Fm., the Scaglia Fm., and the marls of the Langhian (Table 3). In the outermost structure studied of the UMAR (Monte Corona), the fractures of set J2 (n=2) yield consistent precipitation temperatures $T=\mathrm{106}\pm \mathrm{8}$ ^∘C and ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{fluids}}=\mathrm{0.0}\pm \mathrm{1.8}$ ‰ VSMOW; the sample of the set J1 yields $T=\mathrm{56}\pm \mathrm{16}$ ^∘C and ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{fluids}}=-\mathrm{1.1}\pm \mathrm{1.8}$ ‰ VSMOW; in the UMAR, at the Subasio anticline, set F1 (n=3) returns temperatures T ranging from 80±5 to 141±19 ^∘C and a corresponding δ¹⁸O_fluids ranging from 8.4±1.0 ‰ to 16.1±2.1 ‰ VSMOW, while the set J2 yields $T=\mathrm{71}\pm \mathrm{0}$ ^∘C and ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{fluids}}=-\mathrm{5.2}\pm \mathrm{0.0}$ ‰ VSMOW; in Monte Nero, set J1 (n=2) yields consistent $T=\mathrm{30}\pm \mathrm{15}$ ^∘C and δ¹⁸O_fluids= [2.7±2.4, 6.8±0.2] ‰ VSMOW; in the syncline on the west of Monte San Vicino, set F2 (n=2) returns $T=\mathrm{36}\pm \mathrm{4}$ to 70±7 ^∘C and ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{fluids}}=\mathrm{2.5}\pm \mathrm{0.7}$ ‰ to 8.3±1.2 ‰ VSMOW; in Monte San Vicino, set J1 yields $T=\mathrm{47}\pm \mathrm{5}$ ^∘C and ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{fluids}}=\mathrm{3}\pm \mathrm{1.1}$ ‰ VSMOW, while set J2 yields $T=\mathrm{74}\pm \mathrm{10}$ ^∘C and δ¹⁸O${}_{\mathrm{fluids}}=\mathrm{7.2}\pm \mathrm{1.6}$ ‰ VSMOW. In Monte Catria, the sample of set J1 was characterized by a fluid with δ¹⁸O${}_{\mathrm{fluids}}=\mathrm{8.6}\pm \mathrm{0.7}$ ‰ VSMOW precipitated at $T=\mathrm{44}\pm \mathrm{4}$ ^∘C, the sample of set J2 by a fluid with ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{fluids}}=\mathrm{11.1}\pm \mathrm{2.3}$ ‰ VSMOW precipitated at $T=\mathrm{59}\pm \mathrm{10}$ ^∘C. In the easternmost structure (Monte Conero), the sample of J1 was characterized by a fluid with ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{fluids}}=\mathrm{5.8}\pm \mathrm{2.4}$ ‰ VSMOW precipitated at $T=\mathrm{42}\pm \mathrm{12}$ ^∘C.

Table 3 Fluid precipitation temperature and oxygen isotopic values derived from clumped-isotope analysis results.

Figure 8 Plot of δ¹⁸O_fluid (‰ VSMOW) vs precipitation temperature (^∘C) obtained from clumped-isotope analyses along with uncertainties reported as greyed contours. Dotted oblique lines are the measured δ¹⁸O${}_{{\mathrm{CaCO}}_{\mathrm{3}}}$ of the vein cements (‰ VPDB). Shape of the points correspond to tectonic set (layer-parallel shortening being U1 and compatible faults and folding U2), while filling colour relates to structure. LPS stands for layer-parallel shortening; LSFT stands for late-stage fold tightening.

4.6 U–Pb absolute dating of veins and faults

All samples from veins, whatever the set they belong to, are revealed to have a U∕Pb ratio not high enough to return an age, with too low a U content and/or dominated by common lead (see the Supplement), which seems to be common in tectonic veins (Roberts et al., 2020). Of the 35 samples screened, two faults returned variable U∕Pb and Pb∕Pb ratios, allowing the calculation of accurate ages (FAB5 and FAB6) by applying the robust regression approach (Hoareau et al., 2020). In sample FAB5, the pixels with higher U∕Pb ratios made it possible to obtain identical ages within the limits of uncertainty for the different plots in spite of a majority of pixel values dominated by common lead (5.03±1.2, 4.92±1.3, and 5.28±0.95 Ma for the TW (Tera Wasserburg), the 86TW, and the isochron plot, respectively) (Fig. 9a). The rather large age uncertainties are consistent with the moderately high RSE values, but the d-MSWD values close to 1 indicate good alignment of discretized data (Fig. 9b). In sample FAB6, the mapping approach returned distinct ages according to the plot considered because of low U∕Pb ratios (from 2.17±1.4 to 6.53±2 Ma). Keeping in mind their low reliability, the ages obtained for this sample grossly point toward precipitation younger than ∼8 Ma.

Figure 9 (a) Tera–Wasserburg concordia plot obtained from LA-ICPMS U–Pb dating of FAB5 calcite sample. The age was obtained by robust regression through the U–Pb image pixel values. The scale bar corresponds to the weight of each pixel as determined by robust regression. (b) Same plot but with discretized data represented as ellipses (one ellipse = 60 pixels). The running mean (window = 60 pixels) is also shown as a green line.

5 Interpretation

5.1 Sequence of fracturing events and related regional compressional and extensional trends

The previously defined joint and vein, fault, and stylolite sets were compared and gathered in order to reconstruct the deformation history at the scale of the belt. We interpret the mesostructural network as resulting from three stages of regional deformation, supported by already published fold-scale fracture sequence (Tavani and Cifelli, 2010; Tavani et al., 2008; Petracchini et al., 2012; Beaudoin et al., 2016; Díaz General et al., 2015; Di Naccio et al., 2005; Vignaroli et al., 2013) and in line with the ones observed in most recent studies (see Evans and Fischer, 2012; Tavani et al., 2015a, for reviews):

Layer-parallel shortening stage: chronological relationships suggest that set J1 formed before set J2. Set J1 is kinematically consistent with set S1 that recorded bedding-parallel, NE–SW directed Apenninic contraction, except in some places where sets J1 and S1 rather formed under a slight local rotation/perturbations of the NE–SW directed compression as a result of structural inheritance and/or of the arcuate shape of the fold. Bedding-parallel reverse faults of set F1 also belong to this LPS stage as they are likely to develop at an early stage of fold growth (Tavani et al., 2015).

Folding stage: set J2 reflects local extension perpendicular to fold axis and associated with strata curvature at fold hinges. The extensional trend, i.e. the trend of J2 joints and/or veins, changes as a function of curvature of fold axes in map view. We also interpret the stylolite peaks of which orientation are intermediate between set S1 and S2 (Fig. 4) as related to the folding stage (Roure et al., 2005).

Late-stage fold tightening: some tectonic stylolites with horizontal peaks striking NE–SW (set S2) and some veins and/or joints (set J3) postdate strata tilting, like the strike-slip and reverse faults of set F2. All these mesostructures formed slightly after the fold has locked, still under a NE–SW contractional trend which is now oriented at a high angle to bedding. They mark a late stage of fold tightening, when shortening is no longer accommodated, by e.g. limb rotation.

5.2 Burial depth evolution and timing of contractional deformation

Stylolite roughness inversion applied to bedding-parallel stylolites provides access to the maximum depth experienced by the strata at the time vertical shortening was prevailing over horizontal shortening (σ₁ vertical) (Ebner et al., 2009b; Koehn et al., 2007; Beaudoin et al., 2016, 2019, 2020; Beaudoin and Lacombe, 2018; Rolland et al., 2014; Bertotti et al., 2017). In this study, we compare the depth range returned by the inversion of a population of bedding-parallel stylolites to a local burial model (Fig. 10) reconstructed from the strata thickness documented in wells located in the western–central part of the UMAR (Nero–Catria area) (Centamore et al., 1979; Tavani et al., 2008). The resulting burial curves were constructed from the present-day strata thicknesses corrected from (1) chemical compaction by increasing the thickness by an estimated 18 %, then from (2) physical compaction by using the open-source software backstrip (PetroMehas), considering initial porosity of 70 % for the carbonates and 40 % for the sandstones, and compaction coefficients of 0.58 and 0.30 derived from exponential decrease of porosity with increasing burial for the carbonates and sandstones, respectively (Watts, 2001). The timing of exhumation was further constrained by published paleogeothermometric studies and by the sedimentary record (Caricchi et al., 2014; Mazzoli et al., 2002). To the west, tectonic reconstructions and organic matter paleothermometry applied to the Tuscan Nappe (Caricchi et al., 2014) revealed that most of this unit locally underwent more burial because it was underthrusted below the Ligurian Nappe but that the western front of the Ligurian Nappe did not reach Monte Corona (Caricchi et al., 2014). We therefore consider a unique burial curve for the whole western UMAR, and we project the range of depth values at which individual bedding-parallel stylolite stopped being active on the burial curves of the formations hosting the bedding-parallel stylolites. Recent application of this technique, coupled with absolute dating of vein cements (Beaudoin et al., 2018), showed that the greatest depth that a population of bedding-parallel stylolites recorded was reached nearly at the time corresponding to the age of the oldest layer-parallel shortening-related veins, suggesting that it is possible to constrain the timing at which horizontal principal stress overcame the vertical principal stress, switching from burial-related stress regime (σ₁ vertical) to layer-parallel shortening (σ₁ horizontal) (Beaudoin et al., 2020). In the case of the UMAR, 800 m is the minimum depth at which dissolution stopped along bedding-parallel stylolite planes, regardless of the studied formations. This confirms that burial-related pressure solution (i.e. chemical vertical compaction) likely initiated at depth shallower than 800 m (Ebner et al., 2009b; Rolland et al., 2014; Beaudoin et al., 2019, 2020).

Figure 10 Burial model constructed considering thickness from well data (central UMAR) corrected from both chemical and physical compaction. The range of depths reconstructed from bedding-parallel stylolite roughness inversion (with uncertainty) are reported for each formation as grey shades. The derived corresponding timing and depth of active dissolution are reported on the x axis and left y axis, respectively. The timing of the deformation is reported on the righthand side in the inset. The onset of layer-parallel shortening is deduced from the latest bedding-parallel stylolite to have been active; the effect of the 12 % uncertainty is represented by dashed red lines; the onset of late-stage fold tightening is given by U–Pb dating of fault coating in this study. The timing of folding and post-orogenic extension are reported from published sedimentary data (see text for more detailed explanations). LPS stands for layer-parallel shortening and LSFT for late-stage fold tightening.

Figure 10 also shows that bedding-parallel stylolites were active mainly from the Cretaceous (age of deposition of the platform) until Langhian times (∼15 Ma), which suggests that σ₁ likely switched from vertical to horizontal at ca. 15 Ma. For the sake of simplicity, we will consider this age of 15 Ma for the onset of layer-parallel shortening, but one must keep in mind that taking into account a 12 % uncertainty on the magnitude of the maximum vertical stress derived from stylolite roughness inversion, hence ±12 % on the determined depth, yields a 19–12 Ma possible time span for the onset of layer-parallel shortening (from Burdigalian to Serravalian) (Fig. 10). Syn-folding sedimentary deposits pins the beginning of folding of the UMAR to the Tortonian (11–7.3 Ma) in the west and to the Messinian (7.3–5.3 Ma) in the east (onshore) (Calamita et al., 1994). Consequently, in the central and western parts of the UMAR, we propose that the layer-parallel shortening stage of Apennine contraction lasted about ∼7 Myr (from 15–8 Ma – Langhian to Tortonian) before folding occurred. Absolute dating of faults related to late-stage fold tightening in the central part of the UMAR further indicates that fold development was over by ∼5 Ma, i.e. by the beginning of the Pliocene (5.3–1.75 Ma). We can therefore estimate an average duration of folding in the western–central part of the UMAR of ∼3 Ma. Knowing the oldest record of post-orogenic extensional tectonics in the UMAR is mid-Pliocene (∼3 Ma) (Barchi, 2010), we can also estimate the duration of the late-stage fold tightening to ∼2 Ma. In total, the probable period when the compressive horizontal principal stress σ₁ was higher in magnitude than the vertical stress (i.e. until post-orogenic extension) lasted for 9 Myr in the western–central part of the UMAR. We propose an average duration of fold growth about 3 Ma, quite in accordance with previous attempts to constrain fold growth duration elsewhere using either syntectonic sedimentation (3–10 Ma; Holl and Anastasio, 1993; Anastasio et al., 2018) – up to 24 Ma with quiescent periods in between growth pulses (Masaferro et al., 2002) – or mechanical or kinematic modelling applied to natural cases (1–8 Ma; Suppe et al., 1992; Yamato et al., 2011). The combination of bedding-parallel stylolite inversion, burial models and U–Pb dating of vein cements/fault coatings yields a valuable insight into the timing of the different stages of contraction in a fold-and-thrust belt (Beaudoin et al., 2018).

Figure 11 Conceptual model representing fracture development and regional-scale fluid migration during the formation of the Tuscan Nappes and Umbria–Marche Apennine Ridge. Red areas represent the extent of eastward pulses of hydrothermal fluids. Blue areas represent closed fluid system at the scale of the carbonate reservoirs. Potential effects of flexural event reported during upper Burdigalian and Lower Messinian in the eastern part of the belt (Mazzoli et al., 2002; Tavani et al., 2012) have not been documented in our dataset and, therefore, are not considered in this scenario.

5.3 Paleofluid system evolution

The combined use of bedding-parallel stylolites inversion and burial curves constrains the absolute timing of layer-parallel shortening, folding, and late-stage fold tightening in the UMAR (Fig. 10). The further combination of this calendar with the knowledge of the past geothermal gradient as reconstructed from organic matter studies in the eastern part of the Tuscan Nappe (23 ^∘C km⁻¹; Caricchi et al., 2014) therefore yields the expected temperature within the various strata at the time vein sets J1, J2, and J3 and faults F1 and F2 formed. Then it is possible to identify whether fluids precipitated at thermal equilibrium or not during the Apenninic contraction. Overall, most calcite grains from vein cements show thin (thickness <5 µm) and rectilinear twins, suggesting deformation at temperatures below 170 ^∘C (Ferrill et al., 2004; Lacombe, 2010), in line with the maximum expected temperature reached by the Upper Triassic–Eocene carbonate reservoir (120 ^∘C; Fig. 10). The fact δ¹³C values of veins are very close to the δ¹³C values of the surrounding host rocks (Fig. 6c), while the δ¹⁸O values of the veins are different from the δ¹⁸O values of the host rocks (Fig. 6d), suggests against that the fluids' original isotopic signatures were lost due to rock buffering. The similarity of δ¹³C values of veins compared to local host rock, along with the ⁸⁷Sr∕⁸⁶Sr values of the veins, which are in accordance with the expected values of the host rocks (Fig. 7) (McArthur et al., 2001), points towards very limited exchange between stratigraphic reservoirs and rules out external fluid input into the system. Indeed, other potential fluid reservoirs such as lower Triassic evaporites seawater have ⁸⁷Sr∕⁸⁶Sr values significantly higher (0.70800–0.70820) than the highest ⁸⁷Sr∕⁸⁶Sr values documented in the UMAR (Monte Subasio: 0.70760–0.70769; Fig. 7). These characteristics indicate a closed fluid system in most of the UMAR, with formational fluids precipitating at thermal equilibrium, limited reservoir fluid–host rock interactions in the reservoirs, and limited cross-strata fluid migrations.

When considering δ¹⁸O and Δ₄₇CO₂ values, the folds in the westernmost part of the UMAR, Monte Corona, and Monte Subasio require a different interpretation from the other folds of the UMAR (Fig. 6d). Indeed, the δ¹⁸O_fluids values derived from Δ₄₇CO₂ measurements, which are higher in Monte Subasio (8 ‰ to 16 ‰ VSMOW; Fig. 8) than in the rest of the UMAR (from 0 ‰ to 8 ‰ VSMOW; Fig. 8), suggest a fractionation usually interpreted as the result of rock dissolution during fluid migration (Clayton et al., 1966; Hitchon and Friedman, 1969) if considering an environment with limited connection between reservoirs and no implication of external fluids, i.e. where fluids are sourced locally from the marine carbonates. Thus, δ¹⁸O_fluids values in the western part of the UMAR point towards a higher degree of reservoir fluid–rock interaction there. This interpretation based on limited Δ₄₇CO₂ measurements can be further extended by considering the δ¹⁸O values of the veins, significantly higher than the δ¹⁸O values of their host rocks in the westernmost part of the UMAR, hence supporting a higher degree of fluid–rock interaction in this area regardless of the considered fracture set (Fig. 6). However, δ¹⁸O_fluids values derived from the Δ₄₇CO₂ measurements in the set J2 vein of Monte Subasio return a value of δ¹⁸O_fluids consistent with fluid–rock interaction (11 ‰ VSMOW), while another one returns a significantly lower value of δ¹⁸O_fluids (−5 ‰ VSMOW) that can be interpreted as infiltration of meteoric water (Fig. 8). Such evidence for both formational fluids with higher degree of fluid–rock interactions and meteoric fluids at Monte Subasio lead us to propose that the dataset reflects an alternation of fluid source during folding in the area.

Temperatures of precipitation are consistent with the predicted temperatures of host rocks considering the formation and timing of fracture development in most of the UMAR at all times of deformation (Figs. 6, 7, 10). This is different in Monte Corona and Monte Subasio where veins precipitated from fluids at higher temperature than the predicted temperatures for the host rock. In the Maiolica Fm. at Monte Corona, J2 veins with high value of δ¹⁸O_fluids returned a temperature >100 ^∘C, while the maximum predicted temperature during folding is <90 ^∘C. In the Scaglia Fm. at Monte Subasio, the layer-parallel shortening-related veins (J1) and faults (F1) precipitated from fluids much hotter (105–140^∘) than the predicted temperature during layer-parallel shortening (<70 ^∘C). In contrast, the J2 vein in the Scaglia Fm. at Monte Subasio precipitated at 70 ^∘C, i.e. at thermal equilibrium if considering burial history (Fig. 10). A hydrothermal dolomitizing fluid migration event was documented in the southeastern part of the UMAR (Montagna dei Fiori) and interpreted as vertical fluid migration from deeper Jurassic reservoirs (Mozafari et al., 2019; Storti et al., 2018). But the reconstructed fluid temperatures (100 ^∘C) and δ¹⁸O_fluids (6 ‰ VSMOW) are still much lower than the ones reconstructed from the fluids involved in Monte Subasio (105–140 ^∘C; 9 ‰ to 15 ‰ VSMOW), supporting that a different fluid system was prevailing in the westernmost part of the belt during layer-parallel shortening, folding, and late-stage fold tightening compared to central–eastern UMAR.

6 Discussion

6.1 Fluid origin and engine of fluid migration in the westernmost Umbria–Marche Apennine Ridge

During layer-parallel shortening and folding, a temperature of fluid precipitation of up to 140 ^∘C, i.e. significantly higher than the local host-rock temperature, implies that fluids flowed from depth >4 km, while a high δ¹⁸O_fluids reflects a high degree of reservoir fluid–rock interaction. Considering the subsurface geometry of the UMAR (Fig. 1) and discarding any input of external fluids originated from the lower Triassic formations or from the basement on the basis of ⁸⁷Sr∕⁸⁶Sr values, we propose that the fluids originated from the westward lateral extension of the carbonate platform reservoir that was buried under the Tuscan and Ligurian Nappes (Caricchi et al., 2014; Carboni et al., 2020). The coexistence inside a single deformation stage (layer-parallel shortening or folding) of both local/meteoric fluids and hydrothermal brines that migrated from depths can be explained by transient flushes into the system of hydrothermal fluids flowing from deeply buried parts of the same, stratigraphically continuous, reservoir (Bachu, 1995; Garven, 1995; Machel and Cavell, 1999; Oliver, 1986).

We therefore propose that the fluid system prevailing at Monte Corona and at Monte Subasio reflects an eastward, tectonically driven, flow of hydrothermal fluids under the influence of the Tuscan Nappe (Fig. 11), where shortening was accommodated by stacked nappes detached above the Triassic décollement level. That contrasts with the closed fluid system documented in the remainder of the UMAR, where shortening is distributed on deep-rooted faults. We propose that the long-term migration engine was a difference in hydraulic head due to the lateral variation in the burial depth of the reservoir related to the stacked Tuscan and Ligurian nappes west from the UMAR (up to 4 km for the Scaglia Fm.; Caricchi et al., 2014), which does not affect the UMAR (burial up to 3 km for the Scaglia Fm.; Figs. 10, 11b). The resulting hydraulic gradient allowed for the eastward fluid migration within the reservoir, enhanced by layer-parallel shortening and related fracture development (Roure et al., 2005). As the paleodepth variation was related to the thickness of the nappes stacking rather than to a foreland-type slope, the UMAR would then have formed a “plateau” without any large-scale lateral fluid migration (Fig. 11b). The inferred pulses of hydrothermal fluids (“hot flashes” in the sense of Machel and Cavell, 1999) also imply a strong influence of foreland-ward propagation of contractional deformation on the eastward fluid expellation (Oliver, 1986).

6.2 Influence of tectonic style on fluid flow during deformation history

If considering a thin-skinned tectonic model for the UMAR with shallow, low-angle thrusts rooting into the Triassic evaporitic décollement (Fig. 1) (Bally et al., 1986), one would expect some influence of Triassic fluids signatures to be involved in the reservoir paleohydrology at the time faults were active or during folding, as illustrated in similar salt-detached fold systems in the Pyrenees, in the Appalachians, and in the Sierra Madre Oriental (Lacroix et al., 2011; Travé et al., 2000; Evans and Hobbs, 2003; Evans and Fischer, 2012; Fischer et al., 2009; Smith et al., 2012; Lefticariu et al., 2005). On the other hand, if considering a thick-skinned tectonic model with high-angle thrusts crossing the Triassic down to the basement, it becomes more likely that these thrusts did not act as efficient conduits for deep fluids (evaporitic fluids or basement fluids) as fault damage zones in calcium-sulfates evaporites of the area can remain non permeable, if the displacement along the faults is smaller than the nonpermeable layer thickness (De Paola et al., 2008). This contrasts with paleohydrological studies of basement cored folds, where high-angle thrusts allow hot flashes of hydrothermal fluids into the overlying cover (Beaudoin et al., 2011; Evans and Fischer, 2012) in the absence of evaporites. Thus, the lack of Triassic influence in our paleofluid dataset, especially with respect to the low ⁸⁷Sr∕⁸⁶Sr values, seems to support a thick-skinned tectonic style of deformation in the UMAR (Fig. 11c). This fluid flow model therefore outlines important differences between belts where shortening is localized and accommodated by shallow nappe stacking, typical from thin-skinned belts, and belts where shortening is instead distributed on basement-cored folds related to high-angle thrusts, as commonly encountered in thick-skinned belts (Lacombe and Bellahsen, 2016). Tectonically induced fluid flow during layer-parallel shortening in response to hydraulic gradient and lateral tectonic contraction has also been described in other thin-skinned belts, such as the Canadian Rocky Mountains (Vandeginste et al., 2012; Roure et al., 2010; Machel and Cavell, 1999; Qing and Mountjoy, 1992), or in Venezuela (Schneider et al., 2002, 2004; Roure et al., 2003) where lithospheric bulging was at the origin of the depth difference leading to hydraulic gradient-driven fluid migration. The presented case study shows how stacking of sedimentary units typical of thin-skinned tectonics strongly influences the fluid system beyond the morphological front of the belt and allows occurrence of large-scale fluid migration.

7 Conclusions

Our study of the vein–fault–tectonic stylolite populations distributed in Jurassic to Eocene limestone rocks at the scale of the thin-skinned Tuscan Nappe and presumably thick-skinned Umbria–Marche Apennine Ridge reveals the occurrence of several fracture/stylolite sets that support a three-stage (meso)structural evolution of the Apenninic contraction: (1) layer-parallel shortening is reconstructed by a set of joint/veins striking NE–SW to E–W, perpendicular to the local trend of the fold, alongside with stylolite peaks striking NE–SW and early folding bedding-parallel reverse faults; (2) the folding stage is recorded by fold-parallel mode I joints and veins; (3) the late-stage fold tightening is recorded by post-tilting, late folding stylolite peaks and joints and veins, and also mesoscale reverse and strike-slip faults.

Thanks to burial models coupled to bedding-parallel stylolite paleopiezometry, along with (unfortunately scarce) U–Pb absolute dating of strike-slip faults related to late-stage fold tightening, we were able to reconstruct the timing of the onset and the duration of the Apennine contraction. The layer-parallel shortening likely started by Langhian time (ca. 15 Ma) and lasted for ca. 7 Myr. Folding started by the Tortonian time (ca. 8 Ma, as constrained by the sedimentary record of growth strata) and lasted ca. 3 Myr. Indeed, absolute dating of fault coatings implies that late-stage fold tightening started by the beginning of Pliocene (ca. 5 Ma), itself lasting for 2 Myr before post-orogenic extension affected strata since mid-Pliocene (3 Ma).

Accessing the starting and ending time of deformation in the UMAR also allowed us to predict the depth and expected temperatures of the paleofluids during fracturing assuming fluids precipitated at thermal equilibrium. By characterizing the cements related to sets of veins and faults using δ¹⁸O and δ¹³C, ^87∕86Sr, and Δ₄₇CO₂ values, we show that different paleofluid systems occurred during layer-parallel shortening and folding from west to east of the section. In the westernmost folds of the UMAR located beyond the arrow of the Ligurian Nappe thrusting over the Tuscan Nappe, we highlight a local fluid system with transient large-scale lateral, stratigraphically compartmentalized hydrothermal fluid migration. In contrast, these pulses are not documented in the rest of the UMAR and its foreland, where the fluid system always remained closed. We tentatively relate this change in fluid system to a lateral change in tectonic style of deformation across the belt, from thin skinned in the Tuscan Nappe to rather thick skinned in the UMAR. Beyond regional implications, the combination of stylolite roughness inversion and burial history reconstruction, linked to reliable estimates of the past geothermal gradient, appears as a powerful tool to unravel coupled structural and fluid flow evolution in fold-and-thrust belts.