**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)**

**Aurélie Labeur 1, Nicolas E. Beaudoin 1, Olivier Lacombe 2,\*, Laurent Emmanuel 2, Lorenzo Petracchini 3, Mathieu Daëron 4, Sebastian Klimowicz <sup>2</sup> and Jean-Paul Callot <sup>1</sup>**


**Abstract:** Unravelling the burial-deformation history of sedimentary rocks is prerequisite information to understand the regional tectonic, sedimentary, thermal, and fluid-flow evolution of foreland basins. We use a combination of microstructural analysis, stylolites paleopiezometry, and paleofluid 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 flexure/forebulge; (ii) foldscale 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 quantified by sedimentary stylolite paleopiezometry and projected on the reconstructed burial curve to assess the timing of the contraction. As isotope geochemistry points towards fluid 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.

**Keywords:** Apennines; fold-and-thrust belt; burial and tectonic history; fractures; stylolites; fluid flow; clumped isotope thermometry; paleopiezometry

#### **1. 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–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–11], paleostrain and paleostress reconstructions [12–17], mechanical simulation (e.g., [18,19]), or through geochemical approaches (e.g., [20]). Besides, studies of distributed subseismic fractures demonstrated

**Citation:** Labeur, A.; Beaudoin, N.E.; Lacombe, O.; Emmanuel, L.; Petracchini, L.; Daëron, M.; Klimowicz, S.; Callot, J.-P. 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). *Geosciences* **2021**, *11*, 135. https://doi.org/10.3390/ geosciences11030135

Academic Editors: Domenico Liotta, Giancarlo Molli, Angelo Cipriani and Jesus Martinez-Frias

Received: 30 January 2021 Accepted: 1 March 2021 Published: 13 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

that mesostructures are markers of the local deformation sequences and provide access to the evolution of the associated paleostresses [8–10,15,21–36]. When it comes to the burial history, however, the depth-time paths are reconstructed from temperature-dependent proxies, such as organic matter thermal maturity, temperature-dependent clay minerals (e.g., [37–39]), fluid inclusion microthermometry, or low temperature thermochronology (apatite fission tracks, U-Th/He on apatite crystals [40–46]). All these methods rely on the occurrence of specific features (organic matter, fluid inclusions, specific minerals like zircons, or Ar-rich clays) and return only partial information (e.g., maximum temperature, timing of reaching closure temperature ... ) while requiring assumptions about the past geothermal gradient. The latter is often argued over in fold-and-thrust belts and foreland basins where uplift and erosion are common events, thus it remains strongly challenging to reconstruct the burial at which the deformation occurred [7].

We propose to build on recent methodological advances to better constrain the timing of deformation in fold-and-thrust belts and foreland basins. The development of U-Pb dating on calcite cement enables absolute dating on the tectonic vein filling and synkinematic clay minerals (e.g., [47–49]), clumped isotopes can be reliably used to reconstruct the exact temperature of precipitation under 90–100 ◦C [50], and its combination to fluid inclusion microthermometry allows to reconstruct the temperature-pressure-time path in sedimentary sequence [51]. Moreover, the understanding of how stylolites can be used as stress gauges [52,53] enables reconstruction of the burial experienced by sedimentary rocks (e.g., [7,54–60]). These methods rely on the ubiquitous features of carbonate rocks, and we propose in this contribution to combine paleopiezometry and isotope-based thermometry as a new methodology widely applicable to carbonates in order to assess the timing and depth of burial/contraction during deformation history. We present an original application by reconstructing the burial-deformation history of the carbonate sedimentary sequence of the Cingoli Anticline, an arcuate fold in the Umbria-Marche Apennines Ridge (UMAR, Northern Apennines, central Italy; Figure 1). The Cingoli Anticline is an excellent case study for testing such new methodologies: (i) it is a rather simple symmetrical fold, mainly formed by carbonate rocks exhibiting pervasive fractures and stylolites [61]; (ii) there is an abundant literature that discusses the timing of folding, based on sedimentologic and tectonic studies [62–64]; (iii) the geothermal gradient is known [65]; and (iv) it belongs to a young mountain belt where the past mesostructures received relatively poor attention [60,61,66–70] compared to their active counterparts.

In this study, we aim at identifying, characterizing and dating the mesostructures related to the main stages of deformation, i.e., faults, fractures, and stylolites linked to layerparallel shortening (LPS), folding, and late-stage fold tightening (LSFT). For this purpose, we first carry out a classical field-based mesotructural study. We further reconstruct the burial history of strata by building burial curves using present strata thickness corrected for physical and chemical compaction and by applying the roughness inversion technique to sedimentary stylolites in order to quantify the maximal vertical stress, and hence the maximum burial depth experienced by the sedimentary sequence prior to contraction and during exhumation.

The timing of the deformation is further constrained by the temperatures derived from isotopic data (18O,13C and clumped isotopes), showing a thermal equilibrium with the host rocks considering previously published geothermal gradients, and compared to previous studies carried out in nearby folds [60,71]. This original and transferable approach reveals the burial evolution, the timing of deformation, and the fluid system in the Cingoli Anticline during the Apennine contraction in eastern UMAR.

#### **2. Geological Setting**

#### *2.1. The Umbria Marche Apennine Ridge (UMAR)*

The Apennines fold-and-thrust belt extends from the Po Plain to the Calabrian Arc over a distance of 1500 km, and it is the result of the Eurasian and African Plates convergence [72,73]. The Apennine belt accommodated significant orogenic contraction, estimated

up to 50% through cross section balancing [74], and recorded shortening rates ranging from 6 mm/year up to 15-50 mm/year [74,75]. From a structural point of view, the Apennine belt is characterized by a succession of asymmetrical anticlines with eastward vergence, separated by narrower and often asymmetrical synclines. The Apennines are commonly divided into two main arcs, the Northern and the Southern Apennines arcs, each associated with its own geological and structural characteristics [76–78]. Moreover, this curved belt, with eastward convexity, is increasingly younger from west to east (from Oligocene to Pleistocene; [72,79]). This is the result of a roughly eastward migration of the deformation front (and its associated successive foredeep basin), which is related to the eastward retreating subduction of the Adriatic Plate under the European Plate, superimposed by post-orogenic extension at the rear of the propagating orogenic belt [29,72,80,81].

The UMAR represents the central-southern part of the Northern Apennines arc developed during late Miocene-Pliocene when the area was involved in the Apennines build-up. It is about 450 km long and rises to over 2000 m above sea level [75]. The sedimentary succession observed in the UMAR can be divided into three main units as follows: (i) Upper Triassic evaporites, the thickness of which is estimated at 1000 m before deformations and considered as a décollement level. They unconformably overlie crystalline rocks of the basement, which is barely or not even observable at the outcrop [72]; (ii) the Umbria-Marche carbonate-dominated succession (~2500 m thick), divided in several formations deposited from the earliest Jurassic to Oligocene [72,82–84]; (iii) Miocene hemipelagites and turbidites, deposited above these carbonate rocks, which record the progressive eastward involvement of the Meso-Cenozoic succession into the fold-and-thrust belt [62,85]. Indeed, during the foredeep stage, more than 3000 m of turbidites were deposited ahead of the advancing fold-and-thrust belt.

Both thick- and thin-skinned structural styles of deformation have been proposed for the Apennines. The thin-skinned interpretation considers a disharmonic deformation of the crust with the sedimentary units detached along the Triassic evaporites [86–88]. The thin-skinned model is opposed to the thick-skinned which considers the involvement of the basement during compressional deformation (e.g., [89]) Furthermore, several studies suggest that many thrusts are rooted on inherited pre-orogenic structures, mostly preexisting normal faults formed either during the evolution of the Mesozoic passive margin or during the foreland flexure [90–96].

The UMAR undergoes the following main stages of regional deformation: (i) the forebulge stage consisting of the foreland flexuring [8], dated from late Oligocene-early Miocene in the western part of the ridge (i.e., eastern Tuscany-Monte Nero) and from middle Miocene in the eastern part of the ridge (i.e., Gubbio, San Vicino, and Cingoli areas) [79]; (ii) the LPS event, a pre/early-folding compressional stage NE-SW-oriented related to the Apenninic contraction [66,97–99], occurred by early Miocene to the west, and by middle Pliocene to the east [79]; (iii) the folding stage, started by early Miocene in the western part of the UMAR, and by middle Miocene in the eastern part. This stage is characterized by a maximum stress trending parallel to regional shortening, i.e., NE-SW-oriented [29], and local extension perpendicular to fold axis and associated with strata curvature at fold hinge [60]; (iv) the LSFT, associated with a NE-SW contractional trend. This stage corresponds to the moment when shortening is no longer accommodated, by e.g., limb rotation [60]; (v) a post-orogenic extension, starting by Miocene times in the western UMAR [99] and by early Pliocene times in the eastern part of the ridge [62,63], and continuing today. This extensional stage is associated with NNW-SSE trending related normal faults causing the downfaulting of the fold succession [86,99,100].

#### *2.2. The Cingoli Anticline*

#### 2.2.1. Structural Pattern

The Cingoli Anticline is located in the eastern part of the UMAR along the footwall of the Sibillini Mountains (Figure 1A). This anticline is characterized by an arcuate geometry; its main NW-SE axis evolves toward a N-S orientation in its southern part; it is characterized

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

**Figure 1.** (**A**) Location and simplified geological map of the Cingoli Anticline. Red points represent the sampling sites for sedimentary stylolites analysis, and black points represent the measurement sites. (**B**) Stratigraphic column, not to scale, of the Umbria-Marche area, with thicknesses valid for the Cingoli Anticline area (modified from [61]). (**C**) SW-NE geological cross section through the northern part of the Cingoli Anticline [63].

#### 2.2.2. 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:


#### **3. Materials and Methods**

#### *3.1. 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.

#### *3.2. 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).

#### *3.3. Sedimentary Stylolite Roughness Inversion*

Sedimentary stylolites are pressure-solution surfaces usually developed parallel to bedding in sub-horizontal strata during burial, i.e., when σ<sup>1</sup> 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 Lc to the magnitude of prevalent mean stress σ<sup>m</sup> and differential stress σ<sup>d</sup> in the strata at the time the stylolite stopped to be an active dissolution surface, with Equation (1):

$$\mathcal{L}\_c = \frac{\gamma \mathcal{E}}{\beta \sigma\_m \sigma\_d} \tag{1}$$

where Lc is the crossover length converted to m, E the Young modulus of the rock (in Pa), <sup>γ</sup> is the solid-fluid interfacial energy (in J·m<sup>−</sup>2), and <sup>β</sup> a dimensionless constant depending on the Poisson ratio (ν) and calculated with the relation:

$$
\beta = \frac{\mathbf{v}(1 - 2\mathbf{v})}{\pi} \tag{2}
$$

The mean stress and the differential stress, are defined in Pa according to the Equations (3) and (4):

$$
\sigma\_{\rm m} = \frac{\sigma\_1 + \sigma\_2 + \sigma\_3}{3} \tag{3}
$$

$$
\sigma\_{\rm d} = \sigma\_{\rm l} - \sigma\_{\rm 3} \tag{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–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 Lc 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> σ<sup>h</sup> = σH) in the case of sedimentary stylolites. Thus, the Equation (1) established by [53] is simplified as:

$$
\sigma\_{\rm v}{}^2 = \frac{\gamma \mathcal{E}}{\alpha \mathcal{L}\_{\rm c}} \tag{5}
$$

α being a constant defined as follows:

$$\alpha = \frac{(1 - 2\mathbf{v})(1 + \mathbf{v})^2}{30\pi(1 - \mathbf{v})^2} \tag{6}$$

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 <sup>γ</sup> is well-known and stable, fixed to 0.24 J·m−<sup>2</sup> for dolomite and of 0.32 J·m−<sup>2</sup> for calcite [120], and while the Poisson ratio <sup>ν</sup> 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 σ<sup>v</sup> = σ<sup>1</sup> was translated directly into the burial depth (z) of rocks using the relation

$$
\sigma\_{\mathbf{v}} = \varrho \mathbf{g} \mathbf{z} \tag{7}
$$

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<sup>−</sup>2.

#### *3.4. 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–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 (18O/16O) and carbon (13C/12C) 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 CO2 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.

#### *3.5. 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 CO2 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/cm3) and each batch of acid was used for 7 days. After the cryogenic removal of the water, the evolved CO2 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–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):

$$\text{Zm} = \text{a} + \text{bl44} + \text{ct} + \text{dt}^2 \tag{8}$$

Background-corrected ion current ratios (δ<sup>45</sup> to δ49) were converted to δ13C, δ18O, and "raw" Δ<sup>47</sup> values as described by [131], using the IUPAC oxygen-17 correction parameters. The isotopic composition (δ13C, δ18O) 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 Δ<sup>47</sup> values were then converted to the "absolute" Δ<sup>47</sup> 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 (Nf = 78) and conservatively account for the uncertainties in raw Δ<sup>47</sup> measurements as well as those associated with the conversion to the "absolute" Δ<sup>47</sup> reference frame. The precipitation temperature was calculated using the calibration proposed by [135] and updated by [132].

#### *3.6. 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/m3 for carbonate rocks and 1800 kg/m3 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 Φ<sup>0</sup> the surface porosity, both given by the curves.

$$\Phi = \Phi\_0 \mathbf{e}^{-\mathbf{c}\mathbf{y}} \text{ with } \mathbf{c} = \frac{1}{\mathbf{y}} \ln \left( \frac{\Phi}{\Phi\_0} \right) \tag{9}$$

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.

#### **4. Results**

#### *4.1. 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).

**Figure 2.** Observation of the different types of fractures in optical microscopy. (**A**) N-S fractures, (**B**) N045 fractures and (**C**) N135-160 fractures. Textures and shifts in the matrix were characterized, in order to verify the deformation mode, based on the classification of [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):


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.

**Figure 3.** Location of fracture planes measured on the geological map (GPS locations are provided in the supplementary material, 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).

**Figure 4.** Main fractures orientations (after unfolding) plotted on histograms, stereograms and rose diagrams for (**A**) the whole anticline and (**B**) northern and southern parts of the anticline, discriminating forelimb and backlimb measurements. Three major sets of fractures discriminated according these orientations are represented with their specific color (green: set I, blue: set II, pink: set III). (**C**) Chronological relationships between fractures and stylolites (i.e., abutment and crosscutting), observed at mesoscopic scale.

> Tectonic stylolite peaks are mostly oriented N045 (Figure 5). The continuous change in dip of the stylolite plane from vertical to oblique suggests that part of the tectonic stylolites developed before folding and other after folding. Few prefolding stylolites peaks are oriented N140 (WP CIN 37 and WP CIN 27–28) in the Maiolica and Scaglia in the southern backlimb, and N090 in Calcare Massiccio in the southern part of the anticline (WP CIN 12).

> The inversion of striated faults for stress was carried out in few sites in the anticline (Figure 5):


**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., σ<sup>1</sup> orientation) is given by high pole density zones. Stereodiagrams with fault-slip data and principal stress axes are also reported.

#### *4.2. 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.

#### *4.3. 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 Lc) 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. Lc values are summarized in Table 1, and reported as an interval for each formation studied, considering an uncertainty of 23%:


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



**Table 1.** *Cont.*

∗ Crossover length given within 23% uncertainty. Vertical stress σ<sup>v</sup> given within 12% uncertainty calculated according to Equation (5), considering a Young modulus E = 23 GPa ([139], this study), a Poisson ratio <sup>ν</sup> = 0.25, and interfacial energy <sup>γ</sup> = 0.32 J·m<sup>−</sup>2. Depth calculated using dry density of rock d = 2400 g·m−3, acceleration of gravity g = 9.81 m·s−2.

*4.4. Burial Model*

#### Calculated depths with uncertainty Burial-related pressure solution Ranges of depth from BPS Quat. 80 100 Age (Ma) Burial depth (km) Equivalent temperature (°C - 23°C.km ) -1 Plioc. Middle Triassic Jurassic Cretaceous Pal. 3 2 1 0 250 200 150 50 100 40 20 Upper Triassic Eocene Olig. Miocene L.Trias. 4 60 - -- - - <sup>0</sup> Maiolica Scaglia Variegata Scaglia Rossa 2.2 40 60 30 25 20 10 15 0 30 50 *Set III N140 Set II N45 FOLDING LPS FLEXURE* Scaglia Rossa Maiolica *BURIAL RELATED PRESSURE-SOLUTION* ? Age (Ma) Equivalent temperature (°C - 23°C.km ) -1 *Set I N-S*  **ı<sup>1</sup> horizontal** 2 1.4 1.6 1.8 Burial depth (km) *LSFT* 5 **ı<sup>1</sup> vertical** ? - - - - - -

The burial curves resulting from the backstripping process are presented in Figure 6.

**Figure 6.** Burial model constructed considering thickness from stratigraphic and well data corrected for chemical and physical compaction. The range of depths reconstructed from BPS roughness inversion (with uncertainty shaded in light grey) are reported for each formation as grey levels. The corresponding timing and depth of active dissolution are reported on the x axis and left y axis, respectively. The results of clumped isotope analysis (i.e., temperatures of precipitation of vein cements at thermal equilibrium with the host rock) are reported on the right y axis. The timing of the deformation is reported on the right-hand side in the insert. Onsets of LPS, folding stage and LSFT are deduced from the results of the roughness inversion process applied on sedimentary stylolites, as well as from clumped isotope data.

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

#### *4.5. Oxygen and Carbon Stable Isotopes*

Twenty-eight (28) vein calcite cements and surrounding calcite host-rocks from the Scaglia Rossa were analyzed for δ18O and δ13C (Table 2, Figure 7).


**Table 2.** Results of Stable Isotopic Analyses of Oxygen and Carbon Isotopes.

**Figure 7.** Isotopic data from tectonic veins, faults and host-rocks sampled in the Scaglia Rossa. (**A**) Oxygen versus Carbon stable isotopic ratio (‰VPDB) of host rocks (black squares), and vein cements according to the vein sets. Red dotted frame represents the range of isotopic values documented in the UMAR from the Hettangian to Aquitanian carbonates, black dotted frame represents the range of isotopic values documented in tectonic related fracture fillings at the scale of the range [60]. (**B**) δ13C values of vein cements versus δ13C values of the surrounding host rocks (‰VPDB), according to the vein sets. (**C**) δ18O values of vein cements versus δ18O values of the surrounding host rocks (‰VPDB), according to the vein sets. (**D**) Δ47CO2 measured temperature of precipitation (◦C) versus δ18O values of calcite cements (‰VPDB, oblique dotted lines) and corresponding δ18O values of the related fluids (‰SMOW) calculated from temperature-dependent fractionation equation CaCO3-H2O of [141]. For A-D, tectonic veins are reported as full circles of which color relates to the set they belong to (green: N-S, blue: N045, purple: N140). Note that the circle with red contour on D correspond to an LPS-related fault.

In the host-rock (n = 11), the <sup>δ</sup>18O isotopic values range from −2.95 to −1.45‰VPDB while the δ13C isotopic values range from 1.08 to 3.08‰VPDB. In the calcite veins (n = 28), the <sup>δ</sup>18O isotopic values range from −1.56 to 2.59‰VPDB while the <sup>δ</sup>13C 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), <sup>δ</sup>18O ratio ranges from −1.56 to 2.19‰VPDB while <sup>δ</sup>13C ratio ranges from 0.05 to 2.08‰VPDB; for the set II (N045, n = 11), δ18O ratio ranges from 0.58 to 2.59 ‰VPDB while δ13C ratio ranges from 1.81 to 2.97‰VPDB; for the set III N140 (n = 7), δ18O ratio ranges from −1.12 to 0.79‰VPDB while <sup>δ</sup>13C ratio ranges from 1.1 to 3.23‰VPDB; <sup>δ</sup>18O ratio ranges from −1.57 to 1.93‰VPDB while <sup>δ</sup>13C 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 δ13C 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 <sup>δ</sup>18O 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.

#### *4.6. 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.



#### **5. Interpretation of Results**

*5.1. Sequence of Mesostructures in Relation to Folding*

The main stages of regional deformation, already described in the literature [29,60–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 σ<sup>1</sup> striking perpendicular to the northern part of the UMAR structure axes. This is confirmed by reverse faulting associated with a N045 σ<sup>1</sup> 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.

#### *5.2. Evolution of the Burial Depth*

The calculation of vertical stress involves the use of the following mechanical parameters: (i) crossover lengths Lc 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:


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 (σ<sup>1</sup> 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.

#### *5.3. 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 δ13C values (Figure 7B) yet significant for δ18O 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 δ13C signature), with various but limited degrees of fluid-rock interaction likely related to migration, leading to an increase of the δ18O ratios. The Δ<sup>47</sup> results complement and support this interpretation as the combination of the temperature of precipitation T47 with the δ18O value of the cement yields the δ18O values of the precipitating fluid (Figure 7D) using the temperature dependent equation of fractionation of [141]. The reconstructed δ18O values of the precipitating fluids range from 4.80 to 10.50‰ SMOW, irrespective of the tectonic vein sets where fluids precipitated. Positive δ18O 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 T47 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 δ18Ofluids 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].

#### **6. 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 T47 signature, it is possible to infer the timing of fracture development by comparing the range of precipitation temperatures T47 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):


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 T47 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;


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

**Figure 8.** Interpretative model of the history of deformation in the Cingoli Anticline. The structural evolution of the area is represented in 3D and in map view. For each stage of deformation, the main principal stress σ<sup>1</sup> 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.

#### **7. 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:


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


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: Δ<sup>47</sup> analysis report, detailing analysis and results.

**Author Contributions:** Conceptualization, A.L., N.E.B., O.L., and J.-P.C.; methodology, A.L., N.E.B. and S.K.; formal analysis, A.L., N.E.B., M.D., S.K. and L.E.; fieldwork, O.L., N.E.B., J.-P.C. and L.P.; writing—original draft preparation, A.L.; writing—review and editing, A.L., N.E.B., O.L., L.E., L.P., M.D. and J.-P.C.; project administration, N.E.B.; funding acquisition, N.E.B., O.L., J.-P.C. All authors have read and agreed to the published version of the manuscript.

**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)).

**Acknowledgments:** A.L and N.E.B are funded through the ISITE programme E2S, supported by ANR PIA and Région Nouvelle-Aquitaine. Authors also thank two anonymous reviewers for insightful comments that improved the manuscript and Giancarlo Molli, Angelo Cipriani and Domenico Liotta, for their editorial work.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

