Regional-scale paleoﬂuid 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 ra-diometric constraints on burial evolution and deformation timing, these in foreland fold-and-thrust belts: a reappraisal of the Umbria–Marche Northern


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, 1994Bjørlykke, , 1993Lacombe and Rolland, 2016;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., , 2019 and long-term regional deformation 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-andthrust 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 (socalled 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 postfolding events are kinematically unrelated with folding.
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 stylolitesassociated with layer-parallel shortening and with thrustrelated folding, along with the isotopic measurements (δ 18 O, δ 13 C, 87 Sr/ 86 Sr), clumped-isotope measurements ( 47 CO 2 ), and U-Pb absolute dating of their calcite cements. Without an appraisal of which fracture trends are relevant to the largescale (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.

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  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. 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 inver-sion 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(Scisciani et al., , 2019Butler 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, fartravelled 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).

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.

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 beddingperpendicular 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 (σ 1 ), 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.

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 σ 1 was vertical. As proposed by Schmittbuhl et al. (2004) and later developed by Koehn et al. (2012), Ebner et al. (2009bEbner et al. ( , 2010, Rolland et al. (2014) and Beaudoin et al. (2019Beaudoin et al. ( , 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 (σ d = σ 1 − σ 3 and σ m = σ 1 +σ 2 +σ 3 3 , respectively, in Pa) prevailing in the strata at the time the stylolite stopped to be an active dissolution surface following where E is the Young modulus of the rock (in Pa), γ is the solid-fluid interfacial energy (in J m −2 ), 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 According to the sampled formation, we used the solidfluid interfacial energy γ of 0.24 J m −2 for dolomite and of 0.32 J m −2 for calcite (Wright et al., 2001). As an approximation for the material mechanical properties, we use a Poisson ratio of ν = 0.25 ± 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 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 −3 ; Manger, 1963), without any additional assumption on the past thermal gradient or fluid pressure . This technique has already provided meaningful results in various settings (Bertotti et al., 2017;Rolland et al., 2014;Beaudoin et al., 2019Beaudoin et al., , 2020. 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 in m).

O and C stable isotopes
Calcite cements in tectonic veins related either to layerparallel 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 (    (Beaudoin et al., 2014), Poisson ratio ν = 0.25, and interfacial energy γ = 0.32 J m −2 . Depth calculated using dry density of rock d = 2400 g m −3 , acceleration of the gravitational field g = 9.81 m s −2 . 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 3 PO 4 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 labora-tory 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).

87 Sr/ 86 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 chlo-   ride 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 86 Sr/ 88 Sr = 0.1194. Three analyses of the NBS987 standard run with the samples gave a value of 0.710250 ± 0.000001 (1SD).

Carbonate clumped-isotope paleothermometry ( 47 CO 2 )
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 ( 13 C, 18 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 13 C− 18 O ordering in sample calcite was achieved by measurement of the relative abundance of the 13 C 18 O 16 O isotopologues (mass 47) in acid-evolved CO 2 and is referred to in this paper as 47 CO 2 . 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 2 . The CO 2 gas is then captured in a water/CO 2 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 2 was measured on a dual-inlet Thermo MAT 253 mass spectrometer (MS "Pinta"). The reference gas used is a highpurity CO 2 , with the following reference values: −37.07 ‰ δ 13 C VPDB and 8.9 ‰ δ 18 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 47 CO 2 is corrected in three steps. First, mass spectrometer nonlinearity was corrected by applying a "pressure baseline correction" (Bernasconi et al., 2013). Next, the 47 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 47 was to add an acid correction factor of 0.082 ‰ to obtain a final 47 CO 2 value (Defliese et al., 2015). Temperatures of precipitation were then estimated using the equation of Davies and John (2019). The bulk isotopic value of δ 18 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 48 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).

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 sam- ples. 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 −2 . 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 −1 ), and 10 mL min −1 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.

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 posttilting (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;, 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).

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 %.

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 δ 18 O than the rest of the column (δ 18 O of −5.5 ‰ to −3.5 ‰ versus −3.2 ‰ to −1.0 ‰ VPDB; Fig. 6b). Considering the vein cements, the δ 13 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 δ 13 C values of the veins are similar to the δ 13 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 δ 13 C values of veins and δ 13 C values of host rocks ranges from +4.0 ‰ VPDB to +7.5 ‰ VPDB (Fig. 6c)

87 Sr/ 86 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. 87 Table 2. Analysis of 47 CO 2 returns the precipitation temperature (T ), and the oxygen isotopic values of the mineralizing fluid can be calculated using the δ 18 O of the mineral, the clumpedisotope 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)  7.2 ± 1.6 ‰ VSMOW. In Monte Catria, the sample of set J1 was characterized by a fluid with δ 18 O fluids = 8.6 ± 0.7 ‰ VSMOW precipitated at T = 44 ± 4 • C, the sample of set J2 by a fluid with δ 18 O fluids = 11.1 ± 2.3 ‰ VSMOW precipitated at T = 59 ± 10 • C. In the easternmost structure (Monte Conero), the sample of J1 was characterized by a fluid with δ 18 O fluids = 5.8 ± 2.4 ‰ VSMOW precipitated at T = 42 ± 12 • C.

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.

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 foldscale 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 beddingparallel, 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 strikeslip 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.

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 (σ 1 vertical) (Ebner et al., 2009b;Koehn et al., 2007;Beaudoin et al., 2016Beaudoin et al., , 2019Beaudoin et al., , 2020Rolland 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 , 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 (σ 1 vertical) to layerparallel shortening (σ 1 horizontal) .
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., 2019Beaudoin et al., , 2020. 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 σ 1 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 σ 1 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 beddingparallel 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-andthrust belt .

Paleofluid system evolution
The combined use of bedding-parallel stylolites inversion and burial curves constrains the absolute timing of layerparallel 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 −1 ; 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, suggest- 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.
ing 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 δ 13 C values of veins are very close to the δ 13 C values of the surrounding host rocks (Fig. 6c), while the δ 18 O values of the veins are different from the δ 18 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 δ 13 C values of veins compared to local host rock, along with the 87 Sr/ 86 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 87 Sr/ 86 Sr values significantly higher (0.70800-0.70820) than the highest 87 Sr/ 86 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 δ 18 O and 47 CO 2 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 δ 18 O fluids values derived from 47 CO 2 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, δ 18 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 47 CO 2 measurements can be further extended by considering the δ 18 O values of the veins, significantly higher than the δ 18 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 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. the considered fracture set (Fig. 6). However, δ 18 O fluids values derived from the 47 CO 2 measurements in the set J2 vein of Monte Subasio return a value of δ 18 O fluids consistent with fluid-rock interaction (11 ‰ VSMOW), while another one returns a significantly lower value of δ 18 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 δ 18 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 δ 18 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. 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 δ 18 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 87 Sr/ 86 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 (layerparallel 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).

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 saltdetached 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 thickskinned 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 87 Sr/ 86 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 highangle 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(Schneider et al., , 2004Roure 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.

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 threestage (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 δ 18 O and δ 13 C, 87/86 Sr, and 47 CO 2 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 foldand-thrust belts.
Data availability. All data are available in the tables and in the Supplement.
Author contributions. NEB, OL, DK, AnB, and JPC were involved in the overall writing of the paper led by NEB; NEB, OL, DK, and AnB collected structural data and rock samples in the field; NEB, AL, and OL conducted microstructural inversion; GH, AdB, CJ, MM, NR, ILM, FC, and CP designed experiments and collected the geochemical data and wrote the related parts of the paper and appendices. All authors critically reviewed the multiple drafts of the paper.
Competing interests. The authors declare that they have no conflict of interest.
Special issue statement. This article is part of the special issue "Faults, fractures, and fluid flow in the shallow crust". It is not associated with a conference.