Burial-Deformation History of Folded Rocks Unraveled by Fracture Analysis, Stylolite Paleopiezometry and Vein Cement Geochemistry: A Case Study in the Cingoli Anticline (Umbria-Marche, Northern Apennines)

of by A Case Study in the Abstract: Unravelling the burial-deformation history of sedimentary rocks is prerequisite information to understand the regional tectonic, sedimentary, thermal, and ﬂuid-ﬂow evolution of foreland basins. We use a combination of microstructural analysis, stylolites paleopiezometry, and paleoﬂuid geochemistry to reconstruct the burial-deformation history of the Meso-Cenozoic carbonate sequence of the Cingoli Anticline (Northern Apennines, central Italy). Four major sets of mesostructures were linked to the regional deformation sequence: (i) pre-folding foreland ﬂexure/forebulge; (ii) fold-scale layer-parallel shortening under a N045 σ 1 ; (iii) syn-folding curvature of which the variable trend between the north and the south of the anticline is consistent with the arcuate shape of the anticline; (iv) the late stage of fold tightening. The maximum depth experienced by the strata prior to contraction, up to 1850 m, was quantiﬁed by sedimentary stylolite paleopiezometry and projected on the reconstructed burial curve to assess the timing of the contraction. As isotope geochemistry points towards ﬂuid precipitation at thermal equilibrium, the carbonate clumped isotope thermometry ( ∆ 47 ) considered for each fracture set yields the absolute timing of the development and exhumation of the Cingoli Anticline: layer-parallel shortening occurred from ~6.3 to 5.8 Ma, followed by fold growth that lasted from ~5.8 to 3.9 Ma.


Introduction
The reconstruction of the burial-deformation history of sedimentary rocks is a complex issue but an essential exercise to understand the tectonic and sedimentary history in fold-and-thrust belts and foreland basins, with numerous implications spanning from the evolution of the fluid-flow system and associated resources to the understanding of the long-term behavior of the upper crust [1][2][3][4][5][6][7]. Mesostructures observed in fold-and-thrust belts and related forelands, such as faults, veins, and stylolites provide essential information for understanding the deformation pattern [8]. Numerous studies have indeed linked the development of fractures to the large-scale long-term folding history and geometry, either through a descriptive field-based approach [9][10][11], paleostrain and paleostress reconstructions [12][13][14][15][16][17], mechanical simulation (e.g., [18,19]), or through geochemical approaches (e.g., [20]). Besides, studies of distributed subseismic fractures demonstrated by gently dipping limbs and its hinge is relatively flat. The fold has a spatial extension of about 15 km from north to south, and about 5 km from west to east and it culminates 400 m above the local land (770 m above sea level). An WNW-ESE oriented left-lateral fault and a NNE-SSW oriented right-lateral fault bound the Cingoli Anticline to the north and south, respectively ( Figure 1A). In addition, [101] and [63] suggest that inherited pre-contractional structures striking~N-S (i.e., the Jurassic rifting and the late Miocene foreland flexure) strongly controlled the subsequent contractional tectonic evolution of the area ( Figure 1A,C).

Sedimentary Succession
Where not specified, the lithostratigraphic units described below are traditional and validated units [102], and correspond with formation-rank units.
The deformed units comprise Mesozoic and Cenozoic marine deposits, consisting of evaporites and platform carbonates at the base overlain by pelagic carbonates [63,64,101,103] ( Figure 1B). The succession comprises: (1) Upper Triassic anhydrites and dolostones, grouped in the Anidridi di Burano, unconformably deposited above the continental deposits of the Verrucano and the Hercynian basement; at the top of the Anidridi di Burano, euxinic interstratified marls are present.

Fracture-Stylolite Network Characterization and Striated Fault Planes Analysis
The structural data collection and sampling sites are distributed along the whole anticline ( Figure 1). The characterization of the fracture-stylolite network is based on field observations and measurements and analyses of joints, veins, striated faults, and tectonic stylolites. The dataset comprises more than 2300 orientations of mesostructures from the Cingoli Anticline. Furthermore, abutment and crosscutting relationships were carefully observed and analysed to establish the relative chronology of fracture sets. Fracture orientation data were projected on Schmidt stereodiagrams (lower hemisphere), in the current attitude of the strata (raw) and after unfolding (unfolded). In addition, the major sets of joints and/or veins (i.e., the most documented and representative at fold scale) were grouped and averaged by a Fisher statistical analysis, based on the following assumptions: (i) similar orientation considering natural variability (i.e., within 20 • ); (ii) deformation mode (e.g., opening, shearing, contraction), defined through thin sections observed with optical microscope; and (iii) chronological relationships.
Considering that stylolite peaks grow parallel to the main shortening direction [107], with respect to the distribution of dissolution gradients (i.e., non-soluble particles) [108], the orientation of the horizontal maximum principal stress (σ 1 ) was inferred from the maximum density of the peak orientation of tectonic stylolites. The early, syn-, and late folding sets of mesostructures were discriminated with the OpenStereo software [109]. Data were corrected for bedding attitude, by rotation about a horizontal axis to remove the dip angle of the strata, in order to investigate the relationships between fracturing and folding. Thus, fractures were grouped according to their geographical and structural position (respectively north/south and forelimb/backlimb), and according to the main stages of deformation.
Approximately 40 mesoscale striated faults were also measured to complement this mesostructural analysis, in 3 sites of measurements located in the northern and southern ends of the anticline. We used Angelier's inversion technique [110] which, under specific assumptions [111], allowed us to calculate paleostress orientations (i.e., local trend and plunge of principal stress axes) and stress ratios for each site of measurement.

Rock Mechanical Properties
To calculate the mechanical properties of the studied rocks, we used the Schmidt rebound hammer technique, which is a non-destructive method used for the estimation of the uniaxial compressive strength and Young modulus of concrete and natural rocks (e.g., [112]). It implies the use of a spring-loaded piston (the Schmidt hammer), pressed orthogonally against a surface of rock. The energy created by the resistance of the surface to the impact enables the piston to rebound. The distance traveled by the piston after the rebound is called the rebound value R, which is considered to be a proxy of the surface hardness [112], itself used to quantify uniaxial compressive strength and Young modulus of the rock [113]. A Silver Schmidt OS8200 (manufactured by PROCEQ) was used on 12 sites in various localities and sedimentary units of the fold. Each site consists of 50 to 90 rebounds performed perpendicularly to the surface, in most case lying flat (i.e., the hammer being vertical), from which rebound values were averaged as a single rebound value valid for the site. In order to ensure that this value is free from any outlier due to local heterogeneity, we calculated a moving average incremental mean until the mean rebound value stabilizes (supplementary material, Figure S1, Table S1).

Sedimentary Stylolite Roughness Inversion
Sedimentary stylolites are pressure-solution surfaces usually developed parallel to bedding in sub-horizontal strata during burial, i.e., when σ 1 was vertical. They are frequently observed in sedimentary rocks, and especially in carbonates. The 1D roughness of a track along these dissolution surfaces, i.e., the difference in height between two points along the track (sensu [53]), results from a competition between two forces [53]: (i) roughening forces, i.e., pinning on non-soluble particles in the rocks, and (ii) smoothing forces, associated with the surface energy at scale <1 mm and the elastic energy at scale >1 mm. Two main scaling regimes are discriminated by the stylolite growth model [56,107,108,114] depending on the predominant energy and the associated Hurst exponent (i.e., roughness exponent): the surface energy-controlled scale, characterized by a steep-slope and a roughness exponent of 1.1 ± 0.1, and the elastic energy-controlled scale, associated with a gentle slope and a Hurst exponent between 0.5 and 0.6. At the transition of these scale regimes, the change in roughness exponent is associated with a crossover length, estimated in mm by the signal processing approach. The authors in [53] directly link this crossover length L c to the magnitude of prevalent mean stress σ m and differential stress σ d in the strata at the time the stylolite stopped to be an active dissolution surface, with Equation (1): where L c is the crossover length converted to m, E the Young modulus of the rock (in Pa), γ is the solid-fluid interfacial energy (in J·m −2 ), and β a dimensionless constant depending on the Poisson ratio (ν) and calculated with the relation: The mean stress and the differential stress, are defined in Pa according to the Equations (3) and (4): Several studies successfully applied this approach, establishing the spectral analysis of the roughness of sedimentary stylolites as a robust paleopiezometry tool [54,57,59,60,108,[115][116][117][118], especially because the final roughness is acquired quickly at the end of the stylolite's growth, thus depends only on stress and no longer on strain rate [108].
More than 100 sedimentary stylolites, with peaks perpendicular to the dissolution plane, were sampled in several localities of the Cingoli Anticline ( Figure 1A) in the Maiolica, Scaglia Rossa, and Scaglia Variegata. The samples were cut and polished to better analyze the stylolite track. Each cut was made perpendicular to the plane of the stylolite and scanned with a resolution of 12,800 pixels per inches; the resulting file is an image on which the 1D track was hand drawn at magnifications 200% and 400% for greater precision. Then, each track was analyzed as a periodic signal. Usual analyses involve the Fourier Power Spectrum (FPS) and Average Wavelet Coefficient (AWC) methods [115]. We chose the method of Average Wavelet spectrum with Daubechies D4 wavelets [115,119], which is proven to be more stable and less sensitive to resolution effects. We used a non-linear regression with fixed Hurst coefficients of 0.5 and 1.1, corresponding to elastic and surface regimes, respectively (please refer to supplementary material, Figure S2, for the plots). The uncertainty for this regression approach to estimate L c has been previously estimated to 23% [57]. To calculate the vertical stress magnitude, the horizontal stress was then considered as isotropic (i.e., uniaxial strain hypothesis, σ v > σ h = σ H ) in the case of sedimentary stylolites. Thus, the Equation (1) established by [53] is simplified as: α being a constant defined as follows: To verify this assumption, sections perpendicular to the stylolite plane were cut for several samples and treated following the same inversion method, the isotropy of the horizontal stress implying a constant value of the crossover length regardless of the track direction [53,115].
While the solid-fluid interfacial energy γ is well-known and stable, fixed to 0.24 J·m −2 for dolomite and of 0.32 J·m −2 for calcite [120], and while the Poisson ratio ν is rather stable in carbonates and can be approximated to 0.25 ± 0.05, a major source of uncertainty lies in the values of the Young modulus E [57,121]. With known E, ν, and γ, the uncertainty on the calculated stress is 12% [57]. As in previous works [57,60,117,118] the calculated vertical stress σ v = σ 1 was translated directly into the burial depth (z) of rocks using the relation with ρ the average dry density of overlying rocks and g the gravity acceleration. Indeed, the stylolite roughness can be considered as unaffected by local fluid overpressure because the dissolution is located along a fluidic film [56,108,116], an assumption that remains valid until the system is fluidized [122]. For this reasons, ρ was considered as the average dry rock density for clastic and carbonate sediments (evaluated mean value at 2400 kg·m −3 ), without any additional assumption on the past thermal gradient or fluid pressure [7], and g the gravity acceleration, fixed at 9.8 m·s −2 .

O-C Stable Isotopes
The analysis of stable oxygen and carbon isotopes, associated with the study of the diagenetic state of the rocks, allows for the identification and characterization of fluid generations at the origin of mineralization and vein filling within sedimentary rocks [20,[123][124][125][126][127][128][129]. We focus hereinafter on the calcite cements filling the tectonic veins related either to LPS or to strata curvature at fold hinges. The mineralogy of some host-rock was checked with X-ray diffraction (supplementary material, Figure S3) using a Bruker D2 Phaser diffractometer from the ISTeP laboratory (Sorbonne Université) with X-Ray wavelength of 1.54056 Å. The resulting X-ray patterns showing mostly pure calcite with minor amounts of quartz. The diagenetic state of calcite veins was checked under cathodoluminescence microscopy and were performed on a cold cathode Cathodyne platform (CITL CCL 8200 Mk4) at stable vacuum of 60 mThor, a voltage of 12 kV, and a current of 200 µA, corresponding to the ideal voltage-current conditions to activate the luminescence of the carbonates.
We selected vein cements where the vein texture [130] and the diagenetic state support a single phase of filling, occurring at the same time or soon after fracture development (e.g., elongated blocky, Figure 2).
The isotope ratios of oxygen ( 18 O/ 16 O) and carbon ( 13 C/ 12 C) of samples collected in several localities of the Cingoli Anticline were obtained by Isotope-Ratio Mass Spectrometry (IRMS). The spectrometric equipment couples an automatic sample preparation line (KIEL IV) and an analysis section (DELTA V advantage) from Thermo Fisher Scientific at the ISTeP laboratory (Paris). We first selected veins of which (i) the cement witnessed a single growth step, unaltered by later diagenetic events, checked with cathodoluminescence microscopy; and (ii) the cement was likely related to vein opening, considering antitaxial or elongated blocky textures [130]. Thirty (30) to 50 µg of powder were sampled from the vein and the surrounding host-rocks as well, using a computer assisted micromill drill. The powder was reacted with anhydrous orthophosphoric acid at 70 • C to extract CO 2 gas, itself ionized in the spectrometer. The measured isotopic ratio reported in permil relative to the Vienna Pee Dee Belemnite (‰VPDB) are with an accuracy of 0.05‰ and 0.1‰ for carbon and oxygen, respectively.

Carbonate Clumped-Isotope Paleothermometry (∆ 47 )
All analyses were performed at the Laboratoire des Sciences du Climat et de l'Environnement (LSCE, Gif sur Yvette). Eight samples of homogenized carbonate powder were converted to CO 2 by anhydrous phosphoric acid reaction at 90 • C in a common, stirred acid bath for 15 minutes. Initial phosphoric acid concentration was 103% (1.91 g/cm 3 ) and each batch of acid was used for 7 days. After the cryogenic removal of the water, the evolved CO 2 was helium-flushed at 25 mL/mn through a purification column packed with Porapak Q (50/80 mesh, 1 m length, 2.1 mm ID) and held at −20 • C, and then quantitatively recollected by cryogenic trapping and transferred into an Isoprime 100 dual-inlet mass spectrometer equipped with six Faraday collectors (m/z [44][45][46][47][48][49]. Each analysis took about 2.5 hours, during which the analyte gas and working reference gas were allowed to flow from matching, 10 mL reservoirs into the Nier-type ion source through deactivated fused silica capillaries (65 cm length, 110 µm ID). Every 20 minutes, gas pressures were adjusted to achieve m/z = 44 current of 80 nA, with differences between analyte gas and working gas generally below 0.1 nA. Pressure-dependent background current corrections were measured 12 times for each analysis. All background measurements from a given analytical session are then used to determine a mass-specific relationship linking background intensity (Zm), total m/z = 44 intensity (I44), and time (t): Background-corrected ion current ratios (δ 45 to δ 49 ) were converted to δ 13 C, δ 18 O, and "raw" ∆ 47 values as described by [131], using the IUPAC oxygen-17 correction parameters. The isotopic composition (δ 13 C, δ 18 O) of our working reference gas was computed based on the nominal isotopic composition of carbonate standard ETH-3 [132] and an oxygen-18 acid fractionation factor of 1.00813 [133]. Raw ∆ 47 values were then converted to the "absolute" ∆ 47 reference frame defined by the "ETH" carbonate standards [132] using regression methods detailed by [134]. Full analytical errors are derived from the external reproducibility of unknowns and standards (N f = 78) and conservatively account for the uncertainties in raw ∆ 47 measurements as well as those associated with the conversion to the "absolute" ∆ 47 reference frame. The precipitation temperature was calculated using the calibration proposed by [135] and updated by [132].

Burial Model
The burial model associated with the Cingoli Anticline area was constructed from the open access well data collection of the ViDEPI project selected around the Cingoli Anticline (i.e., Misa1, Rosora1, Burano1, Treia1 wells). Thicknesses were sequentially uncompacted to obtain burial depths of the strata of interest through time. Two kinds of corrections were therefore applied, considering the effects of both physical and chemical compaction. The computer interface used to produce these burial curves is the Backstrip software, which performs 1D backstripping of sedimentary strata [136] and involves several parameters for modeling.
First, the layer thicknesses were corrected for chemical compaction, in order to be referenced in the software. Thickness information was provided by stratigraphic studies or well data (i.e., current formation thicknesses). In this case, thicknesses were corrected for chemical compaction, considering spacing and amplitudes of bedding-parallel stylolites (BPS) for each formation studied. The average number of sedimentary stylolites per meter was estimated from outcrop data. The height of the highest tooth (i.e., the height from tooth to base line) associated with each analyzed stylolite was also computed from the samples. Chemical compaction was then deduced from these two parameters, calculated as their product, and expressed as a percentage of bed thickness.
The second step consisted of defining parameters needed to evaluate physical compaction undergone by the different layers, related to the weight of the sedimentary column and possibly of the water column (according to the type of basin considered). These parameters, such as the dry density ρ and the porosity coefficient c, were defined on the basis of the work by [137,138], as follows: (i) the dry density ρ was chosen at 2700 kg/m 3 for carbonate rocks and 1800 kg/m 3 for other lithologies; (ii) the porosity coefficient c was calculated with the following equations and the porosity-depth curves established by [137,138]. Φ is the porosity at depth y, while Φ 0 the surface porosity, both given by the curves.
Thus, c is equal to 0.58 for carbonate rocks and to 0.3 for other lithologies. Corrections related to the weight of sediments and water being significantly different [137,138], the type of basin was also defined (0 for a marine basin, 1 for a continental basin), in order not to introduce bias into the resulting burial curves.

Fracture-Stylolite Network Characterization and Striated Fault Planes Analysis
The fracture deformation mode was defined through thin sections observed with optical microscope, either by considering the texture (i.e., elongate blocky or crack-seal), or the object shift in the matrix (Figure 2).  [130]. Red arrows represent the direction of the opening (mode I), and green arrows the direction of calcite crystal growth. (D-F) observation of these different sets in cathodoluminescence, indicating a single crystallization phase, synchronous with the mode I opening.
Three major sets of fractures were discriminated on the basis of their average orientation (Figures 3 and 4A,B) whereas their chronological sequence was established through abutment and crosscutting relationships at different scales (from outcrop to thin section): set I gathers bedding-perpendicular joints oriented N180 to N020 (after unfolding). This set is observed over the entire anticline, predates sets II and III, because it is intersected and abutted by a set of stylolites with peaks oriented N045, which are themselves intersected and abutted by the joints/veins of set II ( Figure 4C). -set II gathers joints and veins with N045 ± 10 • orientation, present throughout the study area. They are perpendicular to bedding strike. This set postdates set I and predates set III. -set III comprises bedding-perpendicular N130 to N160-oriented joints parallel to bedding strike (i.e., N135-140 in the North and N160 in the South). They are mainly parallel to the axis of the anticline and crosscut or abut all other joint/vein sets ( Figure 4C). -another set of E-W fractures, poorly represented at the scale of the anticline (i.e., only in the northern part, in three sites of measurements), includes N070 to N110-oriented joints (after unfolding) and perpendicular to the bedding, developed after set II.
Because of the low number of measurements (i.e., 25 of 3000 fractures analyzed) and because they systematically developed near faults (Figure 3), this family of fractures is considered as minor and of local meaning only, and therefore not affiliated to a major set. Consequently, it will not be interpreted thereafter.    Table S2). Each measurement point is associated with two stereodiagrams (lower hemisphere), representing main fracture orientations in current (R) and unfolded attitude (U); on each stereodiagram, the bedding is reported as dashed lines, and fracture planes by solid-colored lines, each color relating to one of the three major fracture sets defined (green: set I, blue: set II, pink: set III). The inversion of striated faults for stress was carried out in few sites in the anticline ( Figure 5): in the northern backlimb, conjugate NW-SE trending reverse faults reveal a compressional stress regime with a σ 1 axis roughly oriented N045; -in the northern forelimb, N170-180-oriented normal faults indicate either an extensional regime with σ 3 oriented N045, or, more likely correspond to tilted oblique-slip reverse faults consistent with a pre-tilting N020 compression; -in the southern backlimb, the few fault-slip data preclude any reliable stress tensor calculation. The dataset is however consistent with a post-tilting σ 1 oriented N045 and σ 3 oriented N135.

ORIENTATION OF σ1
Orientations given by poles/peak's orientation on stereograms (unfolded attitude) Supplementary data (given by litterature) Figure 5. Location and plot of measured tectonic stylolites on the geological map. Each measurement point is associated with two stereodiagrams (lower hemisphere), representing main orientations of tectonic stylolites peaks measured (current and unfolded attitude). On each stereodiagram, the bedding is reported as dashed lines, and peaks orientation (i.e., σ 1 orientation) is given by high pole density zones. Stereodiagrams with fault-slip data and principal stress axes are also reported.

Young Modulus Estimate
Rock elastic properties were measured on flat homogeneous surfaces, for Maiolica, Scaglia Rossa and Scaglia Variegata, corresponding to 12 sites of measurement (n = 1063). For each site, the mean for rebound value R was represented as a function of the number of rebound incorporated in the mean calculation (supplementary material, Figure S1, Table S1); the stabilized R value (represented as a plateau on the graph) is then believed to be corrected from heterogenous effect and outliers, and so represent the rebound value R for the rock studied. Strikingly, the average of these representative R values is similar in the Maiolica, Scaglia Rossa and Scaglia Variegata formations, with values of 45 ± 8.4, 48 ± 5.8 and 46 ± 8.5, respectively. R values were further interpreted as Young moduli following the empiric relationship determined in [113] for sedimentary rocks and return a E value similar for the 3 formations at about 20 GPa, very similar to the one reconstructed from stylolite inversion by [121] of 23 GPa.

Sedimentary Stylolite Roughness Inversion
The stylolite roughness inversion method was applied on 112 BPS sampled in the northern, central and southern parts of the Cingoli Anticline ( Figure 1A), within the Cretaceous to Eocene carbonate formations. The inversion was successful (i.e., returning a value of crossover length L c ) on 77 BPS covering the anticline and distributed as follows: Maiolica (early Cretaceous, n = 56), Scaglia Rossa (late Cretaceous-early Eocene, n = 18) and Scaglia Variegata (middle to late Eocene, n = 3). For several stylolites, this paleopiezometric inversion was applied on two orthogonal tracks, in order to ensure that the stress on the horizontal plane was isotropic. L c values are summarized in Table 1, and reported as an interval for each formation studied, considering an uncertainty of 23%:  40 17 720 * Crossover length given within 23% uncertainty. Vertical stress σ v given within 12% uncertainty calculated according to Equation (5), considering a Young modulus E = 23 GPa ( [139], this study), a 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 gravity g = 9.81 m·s −2 .

Burial Model
The burial curves resulting from the backstripping process are presented in Figure 6. They were reconstructed for the Triassic to Pliocene formations in the Cingoli area, considering: (i) the chemical compaction calculated at 8% for Maiolica, and 3% for Scaglia Rossa and Scaglia Variegata considering spacing and amplitude of BPS (following [108]); (ii) physical compaction by using the open-source software BackStrip [136]. The temperatures linked to these depths were calculated by considering a geothermal gradient of 23 • C·km-1 reconstructed in the outermost western part of the UMAR from organic matter thermal maturity [65] and clay minerals [140] (Figure 6). These curves illustrate a first phase of increasing burial, corresponding to the deepening of the Umbria-Marche basin and a second phase of exhumation since early Pliocene. The maximum burial depths computed for the formations of interest, and equivalent temperatures, can be deduced from the left and right y-axis of Figure 6, respectively. These curves are consistent with models established for the inner part of the belt in the area of the Monte Tancia thrust [71].

Oxygen and Carbon Stable Isotopes
Twenty-eight (28) vein calcite cements and surrounding calcite host-rocks from the Scaglia Rossa were analyzed for δ 18 O and δ 13 C (Table 2, Figure 7).  In the host-rock (n = 11), the δ 18 O isotopic values range from −2.95 to −1.45‰VPDB while the δ 13 C isotopic values range from 1.08 to 3.08‰VPDB. In the calcite veins (n = 28), the δ 18 O isotopic values range from −1.56 to 2.59‰VPDB while the δ 13 C isotopic values range from 0.05 to 3.23‰VPDB. The vein cements show variable isotopic values: for the set I (N-S, n = 10), δ 18 O ratio ranges from −1.56 to 2.19‰VPDB while δ 13 C ratio ranges from 0.05 to 2.08‰VPDB; for the set II (N045, n = 11), δ 18 O ratio ranges from 0.58 to 2.59 ‰VPDB while δ 13 C ratio ranges from 1.81 to 2.97‰VPDB; for the set III N140 (n = 7), δ 18 O ratio ranges from −1.12 to 0.79‰VPDB while δ 13 C ratio ranges from 1.1 to 3.23‰VPDB; δ 18 O ratio ranges from −1.57 to 1.93‰VPDB while δ 13 C ratio ranges from 2.02 to 2.21‰VPDB (Table 2, Figure 7A). In order to account for possible rock buffering effect, the isotopic values of the veins were plotted against isotopic values of the surrounding host-rock, for both carbon ( Figure 7B) and oxygen ( Figure 7C). Results show that most veins have a δ 13 C value similar to their host rock, with a difference ranging from −0.50 to 0.25‰VPDB in all sets except in the set I where this difference reaches −1.05‰VPDB. Considering the difference in δ 18 O values, the results are more scattered, ranging from −0.12 to 5.14‰VPDB. Notably, the difference in the set II is higher than the one in the set III.

Carbonate Clumped-Isotope Paleothermometry (∆47)
Eight of the 9 samples presented in the Supplementary Material were selected as being unambiguously related to a major fracture set. Consequently, 7 samples of vein cements and 1 sample of striated coating of fault plane were selected for ∆47 clumped isotope measurements (Table 3), with ∆47 values ranging from 0.593 ± 0.006‰ to 0.630 ± 0.006‰, (1SE) corresponding to precipitation temperature (T47) ranging from 38.30 ± 1.9 • C to 51.4 ± 2.2 • C (1SE). Veins belonging to different tectonic sets appear to yield distinct temperatures of precipitation, with T47 ranging from 48.7 ± 2.1 • C to 51.4 ± 2.2 • C for the set II (n = 5) and related fault cements, while vein cements from sets I and III have T47 ranging from 38.8 ± 2.0 • C to 45.1 ± 2.1 • C.

Sequence of Mesostructures in Relation to Folding
The main stages of regional deformation, already described in the literature [29,[60][61][62][63]68,121,142], were associated with sets of fracture-stylolite network identified in this domain of the UMAR.
Set I (N-S to N020-oriented fractures) is the oldest set encountered in the Cingoli Anticline. Because of its orientation and opening mode, we propose to interpret it as an along-strike joint set related to the flexure stage associated with forebulge development (sensu [8]).
Vertical, bedding, and fold-axis perpendicular set II joints/veins, associated with early folding stylolites with N045-oriented peaks likely reflects a stage of LPS with σ 1 striking perpendicular to the northern part of the UMAR structure axes. This is confirmed by reverse faulting associated with a N045 σ 1 after unfolding ( Figure 5). Local complexities are interpreted as resulting from LPS related stress perturbation, resulting in a slight local stress rotation in the vicinity of local heterogeneities such as inherited faults (Figure 3, e.g., [23,143,144]). For instance, we interpret the N020 contraction in the northern part of the fold as a local rotation around the WNW-ESE fault. We also consider that stylolites with peaks-oriented E-W documented in the Calcare Massiccio relate to LPS perturbed by the reactivation of N-S striking inherited normal fault.
The joints/veins of set III postdate those of set II ( Figure 4C) and are beddingperpendicular and strike parallel to the local fold axis and bedding strike. We propose to relate this set to the folding stage, reflecting outer-arc extension associated to strata curvature at fold hinge. The~20 • variation of the orientation of this set between the north and south of the fold (N140 in the northern part and N160 in the southern part, Figure 3 and Figure 4) is consistent with the arcuate shape of the fold and then strengthen this interpretation.
Late folding, tectonic stylolites with horizontal peaks striking N045, along with posttilting strike-slip faults ( Figure 5) are interpreted to be related to horizontal NE-SW contraction affecting the strata after the fold was locked, corresponding to LSFT [8,15].
Our results therefore demonstrate that the N045 compression prevailed during the entire contractional history, i.e., from LPS to LSFT.

Evolution of the Burial Depth
The calculation of vertical stress involves the use of the following mechanical parameters: (i) crossover lengths L c values, calculated by considering an uncertainty of 23% and reported in Table 1; (ii) mechanical and chemical parameters, defined in the literature and given above (i.e., Young modulus E, Poisson ratio ν and the solid-fluid interfacial energy γ). Then, the burial depths were calculated from each value of the vertical stress using Equation (2) and rounded to the closest 10 m ( Table 1). The corresponding depth ranges for each formation are: -Maiolica: from 720 ± 85 m to 1840 ± 220 m; -Scaglia Rossa: from 880 ± 100 m to 1590 ± 190 m; -Scaglia Variegata: from 720 ± 85 m to 1100 ± 130 m. These ranges of burial depth correspond to the ranges of depth in which pressure solution along sedimentary stylolites was active, i.e., at the time vertical shortening (σ 1 vertical) was prevailing over horizontal shortening [60,145]. Figure 6 shows these ranges of depth reported on the burial model for comparison. The inversion and modeling data appear to be consistent as the maximum burial recorded by sedimentary stylolites never exceeds the maximum depth of the formation they belong to ( Figure 6). In addition, the largest range of depths is associated with the Maiolica (i.e., the oldest formation); for the Scaglia Rossa and Variegata, intervals are overall narrower, and the more recent the formation, the narrower the depth range and the shallower the depth returned. In the case of the Scaglia Variegata (i.e., the youngest), the maximum burial recorded by the stylolites does not exceed 1200 m. Thus, BPS would not develop between 1200 and 1750 m, i.e., at the maximum burial values associated with the Scaglia Variegata ( Figure 6). However, these data indicate that the burial was continuous from the Cretaceous to the late Miocene until the maximum burial depths of the sedimentary layers studied were reached. The reconstruction of the complete burial/exhumation history in relation to the deformation stages requires the combination of these data with isotope analyses.

Fluid System
The fluid system in the Cingoli Anticline can be partially characterized using stable O, C and clumped isotope dataset. The positive difference of isotopic values between vein cements and related host-rock, being low to null for δ 13 C values ( Figure 7B) yet significant for δ 18 O values ( Figure 7C), argues against rock buffered fluid precipitation as well as diagenesis related to burial. Instead, it strongly suggests that cements precipitated from local fluids originated from the studied sedimentary sequence (hence with identical δ 13 C signature), with various but limited degrees of fluid-rock interaction likely related to migration, leading to an increase of the δ 18 O ratios. The ∆ 47 results complement and support this interpretation as the combination of the temperature of precipitation T 47 with the δ 18 O value of the cement yields the δ 18 O values of the precipitating fluid ( Figure 7D) using the temperature dependent equation of fractionation of [141]. The reconstructed δ 18 O values of the precipitating fluids range from 4.80 to 10.50‰ SMOW, irrespective of the tectonic vein sets where fluids precipitated. Positive δ 18 O values points towards a more or less evolved brine origin for the fluid, which is consistent with a scenario involving a migration of Turonian-Lutetian marine fluids inside the host Scaglia Rossa and precipitating at thermal equilibrium within the host. It is noteworthy that these characteristics of the fluid system in the Cingoli Anticline are consistent with the fluid systems reconstructed in most of the other folds and thrusts of the UMAR, except for the Subasio Anticline [60] and Monte Tancia thrust [71]. Interestingly, our dataset does not document a rock buffering of external fluids as interpreted in the Monte Tancia thrust (southern part of the UMAR) by [71], and it does not reflect the late meteoric derived fluid infiltration documented there and related to the currently active extensional tectonics. When considering the tectonic vein sets, temperature of precipitation T 47 differs significantly between set I, set II, and set III, supporting the interpretation of a thermal equilibrium between local fluids and the host rocks throughout the burial history. Moreover, the difference in δ 18 O fluids values between set II, related to LPS, and set III, related to local curvature of the strata at fold hinge, suggest that the degree of lateral fluid-rock interaction was higher during LPS than during folding, as is the case elsewhere [139].

Discussion
The results of mesostructural and isotopic analyses, together with inversion of the roughness of the sedimentary stylolites for maximum burial depth, have been combined in order to unravel the history of deformation in the Cingoli area. Moreover, because the fluid system appears to be at thermal equilibrium during deformation and that each fracture set has a specific T 47 signature, it is possible to infer the timing of fracture development by comparing the range of precipitation temperatures T 47 measured in the vein cements from each set with burial curves and maximum depths of active pressure-solution reconstructed from inversion of stylolite roughness. The stages of deformation (Figure 8), together with their absolute timing and sequence were therefore characterized, and compared with existing data for this study area [60,71], as well as in other localities of the belt, i.e., the anticlines of Monte Nero [121], Monte Catria [66] and Monte Conero [68]; the ages determined using the isotopic data are given with an uncertainty of 0.2 Ma, due to uncertainties on temperature (± 2 • C) ( Figure 6): (i) pre-contractional stage, marked by burial, vertical compaction and dissolution along BPS recognized at the scale of the fold-and-thrust belt [60,61,66,68,121], under a vertical σ 1 . It lasted until the early Messinian (ca. 6.4 Ma); (ii) this pre-contractional stage is partly coeval with an E-W extension related to the flexure of the Adriatic foreland [63,101], the onset of which is set to early Burdigalian (ca. 21 Ma, as defined by the inflection point of the burial curves), and which ended by middle Messinian (ca. 6.3 Ma, Figure 8). This E-W extension would be at the origin of the development of a network of N-S fractures (set I). Pre-folding N-S striking joints have already been described in this anticline [61], and in other anticlines, Monte Nero [121] and Monte Catria [66], without being related to regional extension. In the Conero anticline, however, [68] related a set of N-S, high angle to bedding, joints and veins associated with normal faults to a flexural event. Further considerations provide the existence of polyphase syndepositional normal faulting: the Barremian [146][147][148] and late Cretaceous phases of stretching (e.g., [95,149,150], well known in the Umbria-Marche-Sabina area, could have developed set I joints in the Maiolica and Scaglia Rossa, as well as the development of stylolites in the Mesozoic rocks. (iii) LPS stage, a pre/early-folding compressional stage with σ 1 NE-SW-oriented related to the Apenninic contraction [66,[97][98][99]. The onset of this stage corresponds to the switch from formerly vertical to horizontal σ 1 , associated with a N045 compression marked by the fractures and stylolites of set II. This stage of deformation is consistent regionally as it has been documented in the Monte Nero [121], Conero [68], and Monte Catria [66] anticlines, and identified in numerous folds of the UMAR [52,54]. Based on T 47 precipitation temperature and burial history, the LPS stage started since the middle Messinian (ca. 6.3 Ma). It has been reported in other case studies that the fold growth and associated underlying ramp activation is likely to be responsible for the uplift we reconstructed ca. 5.8 Ma [60]. Thus, we consider the LPS stage to have lasted from 6.3 Ma to 5.8 Ma; (iv) fold growth stage, characterized by a compression parallel to regional shortening, i.e., NE-SW-oriented [29] and local extension perpendicular to fold axis and related with strata curvature at fold hinge [60]. Based on the dating of the LPS, fold growth started at 5.8 Ma yet the T 47 points towards a related N135 striking joints/veins development during the latest Neogene-early Pliocene (ca. 5.2 to 3.9 Ma, Figure 8). That difference in timing suggests a 0.6 My long fault activity and strata tilting before curvature became high enough to developed outer-arc extension fractures; (v) LSFT, post-dating the fold growth stage and still associated with a NE-SW contractional trend. At this stage shortening is no longer accommodated by e.g., limb rotation [60] and is associated with tectonic stylolites with N045-oriented peaks. The E-W fractures locally measured cannot be associated with this deformation stage, because no consistent chronological relationships with the syn-folding fracture sets were identified. Isotopic analyses (Figure 8) suggest the onset of the LSFT by the Pliocene (ca. 3.9 Ma). Despite the record of recent seismic activity in this northern part of the Apennines, linked to a post-orogenic NE-SW extension [151], the end of this deformation stage can precisely be determined neither from previous studies carried out in Cingoli [61] nor from data collected during this study.
This deformation scenario is in line to the one proposed by [61], that discriminates seven sets of stylolites, of which complexity can be related to more local effect of the fold evolution.
The fold growth duration in Cingoli as constrained by the above results (ca. 1.9 My, from 5.8 to 3.9 Ma) is consistent with that established by [62] using the age of foredeep deposits. The age of the end of the foreland flexure, evaluated at ca. 6.3 Ma in our study, is also consistent within uncertainties with the early Messinian age (ca. 7.2-6.5 Ma) derived from the sedimentary dating of the flexure-related normal faults [63,64]. The abrupt increase in the slope of the burial curves ( Figure 6) likely reflects the initiation of this flexure at about 21 Ma, which is older than the late Burdigalian (ca. 16 Ma) age previously proposed on the basis of the distribution and variations in thickness of the Schlier (Aquitanian-Serravallian) marking the deepening of the basin [63], and older than Serravallian (ca. 13.8 Ma) proposed at the west of the Cingoli Anticline, in the western part of the belt [79]. The horizontal isotropy of the compaction-related sedimentary stylolites, along which dissolution was active until 7 Ma, however, suggests a very limited imprint of the flexure on the magnitude of the horizontal stresses in the Cretaceous-Paleogene rocks until the early Messinian, when the flexure became important enough to cause fractures and large-scale normal faults [63].
The folding duration, estimated in our work at~2 My (from 5.8 to 3.9 Ma) is consistent with the results of previous studies carried out in the more internal areas of the belt, using absolute dating methods to reconstruct the duration and timing of the deformations [60,71]: the difference in timing is in line with [62], and the duration of the folding is consistent with [60,71], which respectively dated the folding from 8 to 5 Ma [60], and from 9 to 7 Ma [71] (i.e., a duration of 2 to 3 My).
This case study illustrates the potential of the combination of mesostructural, paleopiezometric, and isotopic analyses, which reveals a regionally consistent sequence and timing of deformation stages, despite multiple sources of uncertainties. Namely, the refining of the timing of the flexure expression on the sedimentary reservoir is an example of how the study of mesostructures can provide insights on the large-scale structures.
Another example of such upscaling lies in the interpretation of the arcuate shape of the Cingoli Anticline, that the distribution and timing of LPS related veins (set II) bounds to be a primary feature of the fold development, likely linked to the reactivation of an inherited N-S normal fault during folding (Figure 8). For each stage of deformation, the main principal stress σ 1 is represented as red arrows, as well as the associated mesostructures (green: set I, blue: set II, pink: set III) and burial depths recorded by the Scaglia Rossa (deduced from burial curves). Faults are represented by red lines in map view, dotted lines when inactive, and solid lines when active. The inherited normal fault is also represented on the 3D block, by a red plane in transparency when it is inactive.

Conclusions
This work, focused on the Cingoli Anticline in eastern UMAR, shows how the burialdeformation history of folded rocks can be unraveled using an original combination of ubiquitous features of carbonate rocks: fracture analysis, BPS paleopiezometry and vein cement geochemistry. The main conclusions are: different stages of deformation were recognized: (i) E-W extension related to foreland flexure (σ 1 vertical); (ii) N045 oriented LPS; (iii) fold growth; (iv) LSFT, under a horizontal N045 contraction. Mesostructural analyses also support that the arcuate geometry of the Cingoli Anticline is a primary feature, probably linked to the oblique reactivation of a N-S inherited normal fault. -the burial history of strata was reconstructed with high resolution using roughness inversion applied to sedimentary stylolites. Our results highlight that this paleopiezometric technique yields consistent maximum depth estimates down to 2500 m, in agreement with previous studies in the western part of the UMAR. -the timing of deformation, and particularly the duration of the Apenninic contractional stages, was reconstructed from combined paleopiezometric, isotopic and mesostructral data. Following foreland flexure (ca. 21.2 to 6.3 Ma), LPS was dated from middle Messinian to early Pliocene (ca. 6.3 to 5.8 Ma) and fold growth occurred between early and middle Pliocene (ca. 5.8 to 3.9 Ma). The precise duration of LSFT remains out of reach. The duration of the fold growth phase is in line with previous estimates based on other proxies such as K-Ar and U-Pb absolute dating [71]. -the O and C stable isotope signatures and clumped isotopes of ∆ 47 of vein cements imply that the paleofluid system that prevailed during LPS and folding in this structure involve marine local fluids with limited interaction with the host rock, in agreement with earlier findings in the eastern UMAR.
Beyond regional implications, this study demonstrates the high potential of our new approach combining paleopiezometric, isotopic, and mesostructural data to reconstruct the sequence and to constrain the timing not only of local mesoscale deformation, but also of regional tectonic events in an orogenic system. Our results further confirm that the paleopiezometric inversion of the roughness of sedimentary stylolites for the vertical stresses is a reliable and powerful tool to unravel the amount and timing of burial without any assumption about the past geothermal gradient.
Supplementary Materials: The following are available online at https://www.mdpi.com/2076-3 263/11/3/135/s1, Figure S1: graphical representation of the incremental mobile average value of the elastic rebound as a function of the number of measurements incorporated in the calculation, for each measurement site, Table S1: stabilized average rebound value R for the Maiolica, Scaglia Rossa and Scaglia Variegata. Standard deviation and number of measures are detailed for each site of measurement, Figure S2: Average Wavelet analysis of the stylolite roughness for all studied samples, Table S2: GPS coordinates of measurement and sampling sites, Figure S3: interpretated X-ray diffractometry spectrum of the Scaglia Rossa, Figure S4: ∆ 47 analysis report, detailing analysis and results. Funding: This research was funded by the Sorbonne Université, grant number C14313, the Région Nouvelle-Aquitaine, and the Agence Nationale de la Recherche (Projet Investissement d'Avenir (programme E2S)).