Reconstruction of the Temperature Conditions of Burial-Related Pressure Solution by Clumped Isotopes Validates the Analysis of Sedimentary Stylolites Roughness as a Reliable Depth Gauge

Abstract

Rough surfaces known as stylolites are common geological features that are developed by pressure solution, especially in carbonate rocks, where they are used as strain markers and as stress gauges. As applications are developing in various geological settings, questions arise regarding the uncertainties associated with quantitative estimates of paleostress using stylolite roughness. This contribution reports for the first time a measurement of the temperature at which pressure solution was active by applying clumped isotopes thermometry to calcite cement found in jogs linking the tips of the stylolites. This authigenic calcite formed as a redistribution of the surrounding dissolved material by the same dissolution processes that formed the extensive stylolite network. We compare the depth derived from these temperatures to the depth calculated from the vertical stress inversion of a bedding parallel stylolite population documented on a slab of the Calcare Massiccio formation (early Jurassic) formerly collected in the Umbria-Marches Arcuate Ridge (Northern Apennines, Italy). We further validate the coevality between the jog development and the pressure solution by simulating the stress field around the stylolite tip. Calcite clumped isotopes constrain crystallization to temperatures between 35 and 40 °C from a common fluid with a δ¹⁸O signature around −1.3‰ SMOW. Additional δ¹⁸O isotopes on numerous jogs allows the range of precipitation temperature to be extended to from 25 to 53 °C, corresponding to a depth range of 650 to 1900 m. This may be directly compared to the results of stylolite roughness inversion for stress, which predict a range of vertical stress from 14 to 46 MPa, corresponding to depths from 400 to 2000 m. The overall correlation between these two independent depth estimates suggests that sedimentary stylolites can reliably be used as a depth gauge, independently of the thermal gradient. Beyond the method validation, our study also reveals some mechanisms of pressure solution and the associated p,T conditions favouring their development in carbonates.

1. Introduction

Pressure solution is one of the most efficient mechanisms to accommodate volume reduction at low differential stress and at low temperature [1] and can be found in most sedimentary rocks. In carbonates, this mechanism is associated with the ubiquitous development of serrated surfaces called stylolites [2,3,4,5], which can also be found less frequently in silicates [6,7,8]. Stylolites are dissolution surfaces that accumulate non-soluble material while driving the dissolved carbonate away into the rock volume, locally redistributing material and changing the local porosity and permeability of the rock (see [9] for a review). As such, stylolites are puzzling features with debated effects on rock strength, large-scale fluid flow and resource distribution [10,11,12,13,14,15,16,17]. Stylolites are characterized by prominent teeth that develop normal to the direction of the highest applied stress, usually interpreted as the direction of the main principal stress σ₁ [18]. Considering stylolites where peaks strike parallel to the plane, that relationship to σ1 allows two origins for stylolites to be distinguished. The ones with peaks perpendicular to bedding (Bedding Parallel Stylolites, BPSs) relate to compaction under a vertical σ1 (typically during burial). The ones with peaks parallel or at low angle to bedding (tectonic stylolites) develop during tectonic contraction associated with orogens development.

Consequently, stylolites are primarily used in structural geology as markers of paleostress directions (e.g., [19,20,21]), as well as for estimating the amount of shortening experienced by the rock [22,23,24]. The morphology of stylolites has been studied to propose classifications [3,4,25] and to discuss the origin of their wavelength [26]. Specifically, analysing the roughness of stylolites (i.e., the difference in height between two consecutive points along the trace) led to the establishment of a new paleopiezometry technique [27,28], known as the Stylolite Roughness Inversion Technique (SRIT).

SRIT has been widely applied across various tectonic contexts, ranging from sedimentary basins [29,30,31] and offshore passive margins [32,33] to broken foreland basins [34] and fold-and-thrust belts [35,36,37,38,39,40,41], being, so far, efficient at reconstructing depth between 400 m and 4000 m. Despite a few attempts to apply SRIT to tectonic stylolites for reconstructing past tectonic stress magnitudes [35,41,42], it has predominantly been used for burial pressure stylolites (BPSs). The cessation of a single stylolite’s development can be attributed to several processes: (1) those inhibiting dissolution, such as local saturation by insoluble materials (e.g., [9]) or the clogging of surrounding pore space by material redistribution (e.g., [23,43]); and/or (2) changes in the stress field, such as local fluid overpressure (e.g., [33]) or large-scale alterations in the far-field stress (e.g., [44]). When applied to a population of BPSs, SRIT provides a range of depths at which pressure solution was active. This range can be linked to various geological events. In sedimentary, poorly deformed areas, the average vertical stress reconstructed often aligns with the maximum burial depth estimated using independent proxies, such as the thermal evolution of organic matter (e.g., [31]). In orogenic forelands, however, the highest stress values consistently recorded by stylolites often differ from the maximum burial depth experienced by the host rock [34,36,37,39,45]. In such cases, this suggests that the maximum depth recorded by BPS development reflects the point in burial history when either the studied formation became saturated and ceased developing stylolites, the temperature and differential stress conditions favoured an alternative deformation mechanism or the far-field stress changed, causing σ1 to become non-vertical. The latter scenario is typically favoured, enabling researchers to predict the timing of the transition from vertical to horizontal σ1 by comparing the depth range obtained through SRIT with burial-time models.

This approach provides a means to constrain the onset of contraction that later leads to layer-parallel shortening and thrusting/folding [45]. In some case studies, the predicted timing has been successfully compared with the age of the oldest deformation features obtained through U-Pb radiochronology on syntectonic calcite in veins and faults [34,36]. Using this workflow, SRIT is emerging as a valuable tool for reconstructing basin and deformation histories, allowing for precise identification of contraction onset before significant deformation of the strata occurs. However, the relationship between stylolite roughness, vertical stress application and depth of stylolite development has yet to be directly validated by comparing depth predictions from SRIT with thermometry data obtained from the same object.

We examine a slab of Jurassic limestone of the Massiccio Formation recollected in the Umbria-Marche Arcuate Ridge (UMAR) in the Northern Apennines [46]. This slab hosts a population of presumed sedimentary stylolites, forming a network of either isolated or anastomosed stylolites [47], which are parallel to one another (Figure 1A). The network also features cemented cracks and jogs that link the tips of the stylolite (Figure 1B and Figure 2a). Aharonov and Karcz [48] associated the development of these jogs and their calcite filling with stylolite growth, making the calcite filling a suitable candidate for estimating the local temperature at the time and depth that the pressure solution occurred.

The present study confronts the vertical stress and derived depth reconstructed from applying SRIT on the population of BPSs to the depth of stylolite growth estimated from calcite clumped isotope (Δ₄₇) thermometry on redistributed calcite filling these jogs. Our findings demonstrate a first-order consistency between the two methodologies, further reinforcing the use of SRIT on BPSs as a reliable tool for reconstructing the depth at which pressure solution was active. Beyond providing methodological validation of SRIT as a depth gauge, this unique dataset offers valuable insights into the pressure and temperature conditions under which stylolites preferentially develop in limestones.

2. Geological Context

The Apennines issued during the Neogene from the convergence between Eurasia and Africa [50,51] and are still active nowadays in the Adriatic basin, while undergoing post-orogenic extension in the hinterland. They developed in the context of the westward retreat of a subduction of the Adriatic plate under the European plate [52]. The Apennines are divided into two main arcs [53], the Northern Apennines and the Southern Apennines, with the transition located in central Italy, south from the Umbria Marche Arcuate Ridge (UMAR, Figure 1A). The UMAR is a fold-and-thrust belt comprising piggy-back duplex folds likely developed on a Triassic décollement level ([49,54], Figure 1A). It developed in an eastward progressive deformation initiated in the Tortonian time in the west (Subasio Anticline) and in the Messinian time in the eastern part onshore (Conero Anticline) [55], as recorded by syn-tilting sedimentary records across the belt ([56,57]) and by absolute U-Pb radiochronology of tectonic veins calcite [45].

Across these folds, the whole of the Jurassic–Neogene sedimentary column of the so-called Umbria-Marche series crops out (Figure 1C). Above the Triassic, evaporites deposited a thick carbonate platform that forms most of the sedimentary stack, starting with the object of this study, the Calcare Massiccio Fm. (Hettangian to Sinemurian), dolomitic limestones deposited as a shallow pelagic platform characterized by benthic fauna which was deposited between 30 m to 200 m depth under the sea level, and which underwent a very slow subsidence until the Cretaceous period [58] and partial dolomitization at a regional scale [59].

The Calcare Massiccio is overlain by the Corniola limestones (Sinemurian to Pleisbachien), the Calcare Diaspirini limestones with cherts and marls (Bajocian to Tithonian), the Maiolica white limestones with cherts (Tithonian to Aptian), the Fucoidi marly limestones (Aptian to Albian), the Scaglia Bianca white marly limestones (Cenomanian), the Scaglia Rossa pink marly limestones (Turonian–Priabonian) and the Scaglia Cinerea grey marly limestones (Priabonian–Aquitanian). This platform is ca. 2100 m thick and has been buried below 300 m of Schlier pre-orogenic shales of the Schlier (Burdigalian to Langhian) and below up to 3000 m of Marnoso–Aranacea clay-rich limestones and silts deposited during the Miocene (Aquitanian–Tortonian) when it was the foredeep of the developing fold-and-thrust belt. In the eastern part of the ridge (Marche area), thicker foredeep deposits are Messinian to Pliocene.

The burial history of the westernmost UMAR was investigated using organic matter thermal evolution by [60] before being extended to the whole belt by [36,37,39] (Figure 1D). Based on the tectonostratigraphy of the area, the authors reconstructed a peak burial for the top of the Calcare Massiccio Fm. ranging from 3000 to 4000 m occurring between 8 My and 4 My from the west to the east of the belt [39]. The roughness of bedding parallel stylolites was studied in numerous folds of the UMAR by means of the Stylolite Roughness Inversion Technique, namely the Subasio anticline [37], the Monte Nero anticline [35], the Cingoli anticline [36], the Monte San Vicino anticline [39,45] and the Conero anticline [39]. Being mainly carried on the Maiolica and Scaglia Rossa Fms., SRIT results highlight that pressure solution was mainly active between 600 m to 2200 m if considering the whole population documented at the UMAR, with discrepancies when considering folds (Figure 1D). Indeed, in the central part of the belt, the range of depth is shallower and smaller, with BPS development limited from 500 to 1600 m, while it is deeper and bigger in the innermost structure (Subasio anticline, 1000–2900 m) and in the outermost structure (Conero anticline, 1100–2200 m). These differences in burial depth, once confronted with the local burial models, reconstructed a timing of the onset of contraction consistent with the sedimentological and radiogenic records across the belt [39].

3. Materials and Methods

3.1. Slab of Calcare Massiccio Fm

The study was conducted on a 3 m long, 1.2 m wide and 2 cm thick slab of dolomitic limestones from the Lower Jurassic Calcare Massiccio Fm., from the Umbria-Marches Arcuate Ridge in Italy. The exact location of the sample being unknown, a survey of the active and inactive quarries of Calcare Massiccio in the UMAR reveals that this slab could have been collected in any of the nine locations, all located in the central part of the UMAR, i.e., between the Subasio Anticline in the West and the San Vicino Anticline in the East (Figure 1A). These locations also all belong to the core of thrust-related anticlines, where the Calcare Massiccio crops out. It is important to note that the burial and tectonic history of this area is roughly similar in amplitude and timing (Figure 1D, [39]). The slab shows vuggy pores filled with calcite, leaving no porosity visible under the microscope (Figure 2). The slab is characterized by meter-scale stylolites (Figure 2a,b). Based on the observed morphology, stylolites fall in the suture and sharp peak descriptive class [25], meaning they are poorly suited to reconstruct the accurate amount of compaction they accommodated but are the best suited to reconstruct the vertical stress they experienced if it was the main principal stress in magnitude [37,39]. A complete description of the shape and fractal properties of the stylolite population of the slab has been published in [46] and further by [47]. The mesostructures also include the following: (1) meter-long veins that crosscut stylolites and strike mostly perpendicular to the stylolite planes; (2) centimetre-scale jogs that developed mainly between stylolite tips and some tension gashes around some stylolite peaks (Figure 3). The relationship between the jogs and the stylolites has been studied by [48], who explained the local orientation as the result of local volume reduction at the tip of a stylolite. The slab of Calcare Massiccio Fm. consists of a grainstone facies made of micritized ooids that are completely dolomitized (Figure 2c).

3.2. Numerical Model of Stress Fields Between Stylolite Tips

In order to test the relation between jogs and stylolites further, we model the stress field between two stylolite tips using the microstructural model Elle [61,62] using a lattice spring approach and inserting the stylolites as rows of particles that have elastic spring constants that are 1/10th of that of the matrix (Figure 4a). In contrast to the dynamic roughening models of Koehn et al. ([22] and references therein), the presented approach is static and similar to that of [48], where the stylolites were modelled as edge dislocations, with strain produced due to (dissolved) missing volume on a half plane. The edge dislocation model of [48] is expected to produce similar strains, and thus similar stresses, as the present weak inclusion model. Our Elle model is deformed vertically under a specified finite strain, and the resulting stress fields are analysed under three distinct geometrical configurations of neighbouring stylolites: neighbouring stylolites that are offset but terminate at the same horizontal position, stylolites whose tips overlap and stylolites with a horizontal gap separating their tips (Figure 4). The mean and differential stress fields near the tips of the modelled stylolites are then compared to a real stylolite example with a vein between the tips.

3.3. Analytical Methods for Past Temperature Estimation

3.3.1. O, C Isotope Measurements

In veins and jogs, the growth of crystals and occurrence of fluid inclusions was checked under optical microscopy on polished thin-sections, and the diagenetic state of the calcite cement was checked using cathodoluminescence microscopy using a cathodoluminescence CITL CCL 8200 Mk4 operating under constant gun condition of 15 kV and 300 µA (Figure 2c–e). For isotope measurements, the material was hand-collected both in calcite veins and dolomite host rocks using a scalpel or a dentist drill to avoid any contamination between vein and host rocks. Sampling included 37 tip-veins and tension gashes and 17 surrounding host rock samples (

Table 1. δ¹³C and δ¹⁸O measurements and calculated temperatures of precipitation.

Sample Name	δ 13C (‰VPDB)	δ 18O (‰VPDB)	Structure	Calc. Temperature * (°C)	Corr. Depth ** (m)
s01a	−3.71	−6.35	Jog
s01b	−2.47	−5.99	Jog
s02a	−2.31	−8.35	Jog	40	1300
s02b	−3.28	−8.55	Jog	41	1350
s02c	−1.90	−8.36	Jog	40	1300
s06a	−1.29	−8.56	Jog	41	1350
s10a	−1.69	−9.66	Jog	48	1650
s10b	−2.72	−10.57	Jog	53	1870
s10c	−0.52	−7.75	Jog	36	1130
s13a	−1.49	−9.38	Jog	46	1570
s15a	−2.83	−5.70	Jog	25	650
s15b	−3.11	−6.68	Jog	30	870
s17a.ii	−0.75	−8.08	Jog	38	1220
s01a	−4.76	−7.14	Jog	33	1000
s01b	−3.80	−6.10	Jog	27	740
s20a	−1.07	−6.25	Jog	28	780
s23a	−2.15	−9.32	Jog	46	1570
s24a	−0.83	−5.42	Jog	23	570
s24b	−1.12	−6.23	Jog	28	780
s27a	−2.24	−6.31	Jog	28	780
s28	−4.94	−6.77	Jog
s28a	−5.29	−6.90	Jog	31	910
s29b	−4.82	−6.56	Jog	29	830
s33a	−3.68	−7.17	Jog	33	1000
s37a	−0.56	−6.61	Jog	30	870
s08a	−1.00	−6.18	Jog	27	740
s03a	−1.86	−9.62	Jog	47	1610
s03b	−1.97	−8.55	Jog	41	1350
s04a	−2.18	−8.29	Jog	39	1260
s09a	−1.18	−8.63	Jog	41	1350
s02v	−8.21	−6.82	vein
s17v	−0.80	−9.34	vein
s23v	1.41	−5.07	vein
s25v	−4.28	−10.26	vein
s06p	−1.47	−8.27	Vug
s32p	−3.60	−6.02	Vug
s03m	1.25	−4.32	Matrix
s12m	1.51	−3.78	Matrix
s13m	1.26	−4.56	Matrix
s17m	1.50	−3.88	Matrix
s01m	1.36	−3.91	Matrix
s23m	1.56	−3.82	Matrix
s25m	1.33	−4.22	Matrix
s33m	1.07	−4.29	matrix
s37m	1.28	−4.27	Matrix
s8m	1.16	−4.12	matrix
s04m	1.61	−3.52	Matrix
s06m	1.61	−3.60	Matrix
s09m	1.45	−3.75	Matrix

* Calculated temperature using a δ¹⁸O value of the fluid of −1.3‰ SMOW and using fractionation equation from [64]. ** Calculated depth using a geothermal gradient of 23 °C/km [60] and a surface temperature of 10 °C. Samples reported in italics are derived from clumped isotopic measurements of Δ₄₇CO₂. Note that all isotopic values of the matrix were corrected as dolomite following the correction after [63].

). Samples were placed in glass vials and reacted with dehydrated phosphoric acid under vacuum at 90 °C, necessitating a correction for dolomite samples according to [63]. Values are reported in permil (‰) relative to the Vienna Peedee belemnite (VPDB/PDB) for carbon and for oxygen with an accuracy of 0.05‰ and 0.1‰, respectively. In order to convert the O isotopic values to temperature of precipitation, we used the relationship linking the δ¹⁸O_fluids, temperature and δ¹⁸O_Calcite established by [64], on the grounds that it is the best suited relationship for abiotic calcite precipitation in non-equilibrium, a classic feature for veins (Figure 5c). To solve this equation, we use a representative O isotopic value of the parental fluid measured from Δ₄₇ (see below).

3.3.2. Carbonate Clumped-Isotope Paleothermometry

Clumped isotope analyses were performed at the clumped isotope facility of the Laboratoire des Sciences du Climat et de l’Environnement (Gif sur Yvettes, France). Four powdered samples comprising only calcite fillings of tension gashes and jogs at the tip of stylolites were powdered and analysed three times each, following a protocol consisting of converting 4 mg of CaCO₃ to CO₂ by reaction with phosphoric acid at 90 °C in a stirred common and stirred acid bath for 15 min. After cryogenic removal of water, CO₂ was passed through a purification column (1 m long, 2.1 mm internal diameter) filled with Porapak Q (50/80 mesh) in a flow of Helium at 25 mL/min, then transferred to an Isoprime 100 dual-inlet mass spectrometer equipped with six Faraday collectors (m/z 44–49), with high-gain amplification (1012 Ω) applied to the m/z 47 cup. Final Δ₄₇ values are reported relative to the I-CDES reference frame [65,66]. A complete description of the methods, raw inputs and data processing is provided as Supplementary Materials.

3.4. Stylolite Roughness Inversion

Stylolite roughness inversion is a recent paleopiezometry technique that is based on the ability of stylolite roughness (i.e., the difference in height of two points along a profile separated by a given horizontal scale) to record the stress magnitude at the time the dissolution ended [27]. This approach is because the growth of some stylolites is favoured by the local dissolution gradient, which creates roughness, while it is inhibited by local forces tending to smooth the dissolution surface. These forces are different according to the observation scale of the roughness and are primarily the surface energy at the small scale, controlled by the chemistry of dissolving material, so by mineralogy, and, at the larger scale (typically above 1 mm), the elastic energy of the material (Hooke’s law), controlled by the rheology of the material and especially the elastic parameters ν (Poisson ratio) and E (Young modulus). It is key to note that, when treated as a signal, stylolite roughness analysis displays two power laws, which display two different slopes in a log–log plot, each characterized by a roughness exponent called the Hurst exponent [67]. The value of the exponent is 1.1 for the power law related to surface energy and 0.5 for the power law related to elastic energy. The observation scale at which the controlling forces switch from surface energy to elastic energy, called crossover length (Lc), is defined as follows:

$$\mathrm{Lc}={\displaystyle \frac{\mathsf{\gamma }\mathrm{E}}{\mathsf{\beta }{\mathsf{\sigma }}_{\mathrm{m}}{\mathsf{\sigma }}_{\mathrm{d}}}}$$

(1)

where Lc is the crossover length (in mm), E the Young modulus of the rock (in Pa), γ the solid–fluid interfacial energy (in J·m⁻²), β a dimensionless constant depending on the Poisson ratio (ν) ($\mathsf{\beta }={\displaystyle \frac{\mathsf{\nu }(1-2\mathsf{\nu })}{\mathsf{\pi }}}$), σ_m the mean stress magnitude (${\mathsf{\sigma }}_{\mathrm{m}}={\displaystyle \frac{{\mathsf{\sigma }}_1+{\mathsf{\sigma }}_2+{\mathsf{\sigma }}_3}{3}})$ and σ_d the differential stress magnitude (${\mathsf{\sigma }}_{\mathrm{d}}={\mathsf{\sigma }}_1-{\mathsf{\sigma }}_3)$. In case of sedimentary stylolites, i.e., developed under a maximum principal stress vertical that is supposedly much higher than the magnitudes of both horizontal stresses, we assume that σ_V >> σ_H = σ_h (uniaxial strain hypothesis), allowing Equation (1) to be simplified as follows [41]:

$${{\mathsf{\sigma }}_{\mathrm{v}}}^2={\displaystyle \frac{\mathsf{\gamma }\mathrm{E}}{\mathsf{\alpha }{\mathrm{L}}_{\mathrm{c}}}}$$

(2)

with α being a constant defined as $\mathsf{\alpha }={\displaystyle \frac{\left(1-2\mathsf{\nu }\right){(1+\mathsf{\nu })}^2}{30\mathsf{\pi }{(1-\mathsf{\nu })}^2}}$.

To find the value of Lc, the idea is to analyse the roughness of the stylolite (in 2D) as a signal using either Fourier Transform (Fourier Power Spectrum, FPS) or Average Wavelet Coefficient (Daubechies D4 wavelets, AWC), two methods that have been proven reliable for fault surface roughness and stylolite roughness [27,28,41,68,69,70]. Signal analysis is conducted on an 8-bit high-resolution image (12,800 dpi) of the stylolite track over a length ranging from a few cm to 5 cm. Scanned samples are digitized either manually on pixel-based drawing software or automatically with colour and shape-based segmentation. Then, the signal analysis is performed using the Matlab (R2022b) scripts provided by [41], modified in [39]. Because the stylolites in the studied slab are red and the host rock is white, the contrast allowed for a computer-assisted segmentation to be performed. For each log–log plot, a non-linear regression with set slopes, in agreement to the expected Hurst Coefficient, is fitted through the binned data and highlights the breakout in the linearity of the regression as being the Lc value, with an uncertainty of 23% [30]. In the present study, Lc was given using AWS signal analysis, and the resulting σ_v were transformed into depth (z) using σ_v = ρgz, where g is the gravitational acceleration (9.81 m/s) and ρ is the density of the overlying sediments. It is worth noting that σ_v, and, consequently, the depth estimated, is given within 12% uncertainty [37].

4. Results

4.1. Characterization of Stylolites and Petrography of Mineralization

The Calcare Massiccio slab exhibits stylolites exclusively of the suture and sharp peak type [25]. Considering stylolite as initiated as a planar surface, the minimum range of compaction can be derived from the amplitude of the biggest tooth on each stylolite (upper limit of the range). Morphometric measurement on 31 stylolites with different amplitudes and lengths return an amplitude of the biggest tooth ranging from 0.6 mm to 11.1 mm, with an average maximum value for the population of 6.3 mm. On a typical 1 m long column of the slab, the stylolite density is of 0.5 (stylolite/cm). Thus, the minimum compaction accommodated by BPSs represents between 3% and 6%, when considering an average or a maximum tooth height, respectively.

Jogs are systematically defined as small cracks bounded on both ends by stylolite tips and completely filled with calcite. The calcite textures are a non-fibrous blocky type (Figure 2c–e and Figure 3A–D). Those textures relate to precipitation that occurred in an opened crack (blocky). This does not exclude that the precipitation occurred very soon after the crack opened. The natural luminescence of the veins is calcite-related bright-to-dull orange, exhibiting zonation at crystal-scale, whereas the host rock is dolomite-related bright red-orange (Figure 2c–e). Crystal growth shows either a syntaxial or blocky type texture showing a growth from the edge to the centre of the vein. Also, some vein cements are affected by tectonic stylolites (Figure 2e), probably related to later tectonic stress, which support that those jogs and gashes developed prior to the maximum stress rotating to the horizontal. The vein cement is rich in primary and secondary fluid inclusions, which do not exhibit vapour bubbles (Figure 3E) or that are stretched (Figure 3F), respectively, suggesting that fluids precipitated at a temperature lower than 80 °C [71] or that the veins were deformed.

Calcite mineralization can also be found in vugs, former pores that were totally filled, and also in veins that crosscut, striking perpendicular to the studied stylolites (Figure 2 and Figure 3A). That systematic relationship, along with the en échelon veins (Figure 3A) that crosscut stylolites, support further the interpretation of the stylolites being sedimentary and related to burial. Also, it sets a relative chronology between the jogs being synchronous to pressure solution and the veins postdating pressure solution. Based on microstructures, it is impossible to state if the calcite filling the pores is anterior, synchronous or posterior to the pressure solution.

4.2. Stylolite Tips and Jog Simulations

The simulations of stylolite tip interactions are illustrated in Figure 4, with three scenarios: perfect tip alignment, overlapping tips and tips separated by a horizontal distance. Plots of the mean stress illustrate that the stylolites act like a LVR (localized volume reduction) defect, similar to an edge dislocation [61,72], with a compressional stress concentration in front of the tip (negative stress) and a tensile (positive) stress concentration behind the tip, in accordance with analytical solutions by [48]. If two stylolite tips are in close proximity, their stress fields may interact. The strongest mutual stress enhancement occurs when the tips are still separated by a small horizontal distance (Figure 4g). In contrast, stress between stylolite tips is only weakly enhanced once the stylolites overlap (Figure 4e). The differential stress plots show similar patterns as an edge dislocation, with the highest differential stresses at angles of 45 to 60 degrees on each side of the stylolite’s tip. As with the mean stress scenario, the region between stylolite tips experiences the largest stress concentration enhancement when there is still a small horizontal gap between the tips (of the same order as the distance between the stylolites) or if the tips are horizontally aligned. The natural example of stylolite tip veins (Figure 3) shows tips that have no horizontal distance between them (scenario c and d in Figure 4) and a vein that is slightly curved. If we overlay the differential stress field from the simulation with the natural example (Figure 4b), the region experiencing the largest concentration of differential stresses is very similar, if not identical, to the position of the vein including its curvature. This similarity between the geometry and actual stress field may indicate that the vein formed because of the high differential stresses between the tips and is potentially a combination between an extension and a shear vein, which is also indicated by its curvature. The vein in Figure 4 is not unique in its agreement with theoretical predictions; refs [23,48] measured orientations of a large number of veins from this slab and showed they agree with analytical predictions of stress from interacting stylolite tips. Since the stress field in the model would develop during stylolite growth in the real system, one could argue that the veins form in this stress field and are therefore of a similar age as the stylolites.

4.3. Stable Isotope Results

The δ¹⁸O and δ¹³C values of 29 calcite filling samples from jogs at the tip of stylolites, of four calcite filling samples from veins cross cutting stylolites and of two calcite filling samples from in a pore (vug), along with δ¹⁸O and δ¹³C values of 17 local dolomite host rock, are reported on a δ¹⁸O versus δ¹³C plot (Figure 5a) and in Table 1. Host-rock dolomite exhibits oxygen isotope signatures scattered within a small range from −4.5‰ to −3.5‰ PDB and carbon isotope signatures ranging from 1.2‰ to 1.6‰ PDB. The tectonic veins (n = 4) crosscutting the BPS have very scattered δ¹⁸O and δ¹³C isotopic values, from −5.07 to −10.26‰ PDB and −8.2 to 1.4‰ PDB, respectively. The vugs’ (n = 2) isotopic ratios range from −3.6 to −1.5‰ PDB and −8.27 to −6.02‰ PDB for δ¹⁸O and δ¹³C, respectively. Jog cement ratios of δ¹⁸O and δ¹³C (n = 30) are scattered over a range of −10.57‰ to −5.7‰ PDB and −4.8‰ to −0.5‰ PDB, respectively. The degree of fluid–rock interactions is assessed by plotting the δ¹⁸O of the vein cements versus the δ¹⁸O of the surrounding host rock in the vicinity (cm scale), using the same HR δ¹⁸O for several calcites when they are close to one another (Figure 5b).

4.4. Clumped Isotope Thermometry

Among the calcite filling from jogs of which carbon and oxygen stable isotopy were analysed, three were selected based on their size to measure the Δ₄₇ (in purple on Figure 5a) along with one cement sample in a vug. The four Δ₄₇ values are very similar (

Table 2. Clumped isotopes analysis.

Sample Name	Type	δ13C (‰VPDB)	δ18O Calcite (‰VPDB)	Δ47-ICDES	95% Uncertainty	T47 (°C)	95% Uncertainty	δ18OH2O (‰VSMOW)	95% Uncertainty
S01a	Jog	−3.71	−6.35	0.5512	0.0103	40.1	4.0	−1.1	0.7
S01b	Jog	−2.47	−5.99	0.5639	0.0104	35.3	3.8	−1.6	0.7
S28	Jog	−4.94	−6.77	0.5511	0.0103	40.1	4.0	−1.5	0.7
S32C	Vug	−3.6	−6.02	0.5350	0.0103	46.6	4.2	0.4	0.7

), ranging from 0.5511‰ to 0.05639‰. Based on the OGLS23 calibration of [73], these values correspond to precipitation temperatures (T₄₇) ranging from 35.3 ± 3.8 °C to 46.6 ± 4.2 °C when including the full dataset and from 35.3 ± 3.8 °C to 40.1 ± 4.0 °C if the vug is excluded. More detailed information about the data is available in the dedicated Supplementary Materials. The fractionation equation of [64] allows the δ¹⁸O signature of the parental fluids from which calcite precipitated to be calculated. For all three jogs, we reconstruct a rather similar δ¹⁸O_fluids (Figure 5; Table 2), ranging from −1.6 ± 0.7‰ SMOW to −1.1 ± 0.7‰ SMOW, i.e., around −1.3 ± 0.7‰ SMOW (95% confidence limits).

4.5. Stylolite Roughness Inversion

Thirty-one (31) stylolites were subjected to SRIT, with both FPS and AWC analysis. The first result is that a majority of SRIT tests conducted with FPS returned invalid results (sensu [39]), i.e., two growth regimes could not be found. In contrast, 87% of the same population returned a valid inversion with AWC, i.e., presented two well defined slopes corresponding to the modelized values of the Hurst exponent (Figure 6). Considering only the results from AWC (

Table 3. Results of SRIT applying to BPS (using AWC) and interpretation in depth.

Sample Name	Inversion Quality	AWS Lc (mm)	2σ	Compaction (mm)	σ1 (MPa)	Uncertainty (MPa)	Depth (m)	Uncertainty (m)
S3	A	0.25	0.06	0.64	41	6	1750	250
S10-1	A	0.42	0.10	4.95	32	4	1350	200
S10-2	A	0.77	0.18	5.82	23	3	1000	100
S23-1	A	1.41	0.32	6.11	17	3	700	150
S23-2	A	1.40	0.32	9.98	17	3	700	150
S30-1	A	2.23	0.51	10.13	14	2	600	100
S37	A	1.09	0.25	6.75	20	2	850	100
S29	A	1.61	0.37	7.08	16	2	700	50
S6	A	1.03	0.24	4.73	20	3	850	150
S33-1	B	1.68	0.39	10.87	16	2	700	50
S33-2	B	0.75	0.17	5.48	24	3	1000	150
S1	B	0.38	0.09		33	4	1400	150
S2-2	A	0.65	0.15	5.59	25	4	1050	200
S20-1	A	0.91	0.21	2.39	21	3	900	100
S28-1	A	1.07	0.25	9.09	20	2	850	100
S8-1	A	0.35	0.08	5.65	35	4	1500	150
S26-2	A	0.23	0.05	2.67	42	6	1800	250
S32-a	A	1.06	0.24	11.17	20	2	850	100
S32-b	A	0.42	0.10	6.10	32	4	1350	200
S32-c	A	0.94	0.22	4.69	21	3	900	100
S32-d	A	0.20	0.05	5.35	46	6	1950	250
S26-1	A	0.36	0.08	9.66	34	5	1450	200
S2-3	A	0.26	0.06	4.37	40	5	1700	200
s9-1	A	0.26	0.06	6.42	40	5	1700	200
S19-1	A	0.36	0.08	6.14	34	5	1450	200
s26-1	A	0.31	0.07	9.00	37	5	1550	250
s16-1	A	1.08	0.25	6.88	20	2	850	100
s16-2	A	0.39	0.09	4.36	33	4	1400	150
s16-3-g	A	0.30	0.07	4.38	37	5	1550	250
s16-3-d	A	0.30	0.07	7.57	37	5	1550	250
s16-4	A	0.19	0.04	1.20	46	7	1950	300

Inversion conducted following Equation (2) with E = 24.2 GPa, ν = 0.2, ϒ = 0.247 J·m⁻². Conversion to depth using ρ = 2400 g·m⁻³. All stress and depth results are given rounded up.

), Lc values are distributed (Figure 7) as follows: Min: 0.19 mm; first quartile: 0.30 mm. Median: 0.42 mm; average: 0.73 mm; third quartile: 1.07 mm. Max: 2.23 mm. (All values given within 23% uncertainty.) Applying Equation (2), with a Young modulus of 24.2 GPa corresponding to the one used for the sedimentary column of the UMAR ([37]), a Poisson coefficient of 0.2 and a solid–fluid interfacial energy for Dolomite at 0.247 J·m⁻² [74], the vertical stress σv gives a gaussian distribution centred around a value of 29 MPa with a minimum of 14 MPa, a first quartile of 20 MPa, a median of 31 MPa, a third quartile of 37 MPa and a maximum of 46 MPa.

5. Interpretation

5.1. Depth of Calcite Precipitation in Jogs

Jogs have been filled as a single event, either coeval or close in time to their opening, as suggested by the petrographic analysis and the numerical modelling (Figure 2, Figure 3 and Figure 4). Jogs are a common feature in stylolite development (e.g., [75,76,77]) and appear concomitant to the pressure solution. Indeed, in the mechanical simulation we provide, the stress patterns overlap the observed jogs, in line with the study of [48] on the same material. Simulation predicts a stress field allowing for dilation and shear, in contrast with the idea of stylolites as “anti-cracks” [72,78], where compressive stress surrounds the tip from all directions. Then, we can safely assume that the calcite filling precipitated from a fluid at thermal equilibrium with the host rock, as 1—the jog developed during pressure solution; 2—the limited size of the jogs prevented large-scale fluid migrations; and 3—the drop of local pressure due to crack development can both pump the fluids in the cavity and trigger the precipitation of calcite.

In terms of isotopic signatures, two trends can be defined (Figure 5a): One reflects a O and C fractionation related to the temperature increase of the fluid as the host-rock temperature increases during burial. This is a classical trend in diagenesis, where the temperature drives a strong fractionation of O (e.g., [79]) while the fractionation of C remains limited [80]. The second trend consists of data with light δ¹⁸O and heavy δ¹³C values, irrespective of the amount of compaction along the stylolite (Table 3). These are harder to explain, but it can be related to biogeochemical process, as there are numerous processes that can produce CO₂ and any depth and temperature and during dissolution [81,82,83,84,85,86,87]. The produced CO₂ is then dissolved into the fluid that precipitates in the jogs. As the data with a heavy δ¹³C have a light δ¹⁸O, we speculate that these data reflect some shallow depth bacterial activity (e.g., [88,89]). In any case, as these trends are defined based on δ¹³C, they do not affect the estimate of the temperature and depth presented below as it is only related to the δ¹⁸O values.

To narrow down the range of temperature at which pressure solution was an effective mechanism for accommodating compaction, we use the temperature of precipitation for calcite in the jogs, because jogs have been proven to be related to stylolite growth ([48] and Figure 3), and so the likelihood of being filled by a local redistribution of calcite coming from the stylolite plane is very high. In order to check potential solid-state reordering [65,90,91], we performed modelling using the ClumpyCool script [92] with the time–temperature paths described on Figure 1D (using the model for Monte San Vicino). As expected, the potential reordering would alter the temperature by less than 3 °C, as the maximum temperature reached by the strata is below 100 °C and as the strata stay below 80 °C during most of the time. Hence, the temperature of precipitation has been directly constrained by clumped isotope thermometry for four calcite samples in three jogs and one vug. This limited dataset can be indirectly extended to the whole dataset of O, C isotopes from jogs, assuming that all are precipitated from a fluid with a common oxygen isotopic signature. This assumption is based on the reasonable consideration that jogs developed during stylolite growth and that they are not likely to be drains for external fluids.

Using the representative δ¹⁸O_fluids value of −1.3‰ SMOW, close to the one expected for seawater (0‰ SMOW), we calculate the corresponding temperature using the relationship established by [64] (Table 1). For all the samples, the corresponding precipitation temperature, reported in red on Figure 5c, ranges from 25 °C to 53 °C. To transform that temperature into depth, one needs to (1) assume that the fluid precipitation occurred at thermal equilibrium with the geotherm and (2) use a past geothermal gradient valid for the area at the time of deformation. The first point is a sound assumption considering that the calcite is locally redistributed from pressure solution during the growth. The second point was estimated by means of vitrinite reflectance in the southern part of the central Apennines to be around 23 °C/km (Tuscan Nappe, [60]). In order to predict a depth with such data, it is important to consider the surface temperature at the time of precipitation, i.e., the stable temperature of the water table. On land, this is very variable and dependant on the climate [93], while being much more stable at the bottom of the sea [94]. The UMAR carbonate platform developed in a warm-humid climate [95,96], so, here, we used a surface temperature of 10 °C, i.e., corresponding to the stable temperature at the range of depth (50–300 m) for current carbonate platform construction in a tropical climate (e.g., [97]). The predicted depth for the range of temperatures reconstructed for calcite precipitation ranges from 570 m to 1870 m, with a median value of 1175 m (Figure 7a).

5.2. Depth of Active Pressure Solution from SRIT

SRIT on the BPS population has been successfully applied in the UMAR part of the Northern Apennines where the slab comes from and in the Jurassic to Palaeocene sedimentary sequence in the Apennines that mainly comprises carbonate rocks [35,36,37,39]. The method has returned a range of depths at which pressure solution was active in several contexts. The previous studies set the input parameters for the inversion of Lc values to vertical stress magnitudes, and they also provide the density of the dry sedimentary column, which was averaged for the whole UMAR at around 2400 g/m³ [35]. We use the dry density because it has been established that, because it develops at a fluid rock interface, the stylolite surface does not see the effect of fluid pressure in the overlying sediments [9]. As a result, the depth distribution for the stylolite population ranges from 600 m to 1950 m, with a median value of 1350 m (Figure 7).

6. Discussion

6.1. Independent Validation of SRIT on BPSs as a Reliable Vertical Stress Gauge

One of the main objectives of this study is to validate the stress magnitudes and related paleodepths obtained from applying SRIT on a population of BPSs. The slab of Jurassic Massiccio Fm studied by [46,47,48] was a perfect candidate to do so thanks to the common occurrence of syn-stylolitization fracture development at the tip of the BPSs, with a size that was large enough to allow for geochemical analysis on powdered samples. To our knowledge, this was the first study that used these overlooked objects to reconstruct the temperature and then depth at which pressure solution was active. Also, the roughness of the stylolite population associated with jogs returned 27/31 stylolites of which the roughness follows the model proposed by [27,28], i.e., they return two growth regimes characterized by a set Hurst exponent and separated by a crossover length Lc. This contrasts with the work of [46] who worked on the same slab, highlighting that these stylolites are fractal, i.e., their roughness is governed by a single power-law showing self-affine properties at all scales of observation. However, that outcome was from the study of only five stylolites, of which two did show a transition between two self-affine growth regimes at a crossover length Lc (CMa and CMc, Figure 10 of [46]), with Lc values of ca. 0.9 mm and 0.75 mm. These are in line with crossover length values we find in the present study.

Recently, SRIT applied to BPSs enabled an independent estimate of the timing of contraction and deformation in fold-and-thrust belts [45,98], when the BPS population records a shallower depth than the maximum depth reached by the deformed unit, hinting at a time/depth at which σ₁ is no longer vertical. The agreement of the SRIT with the estimated timing of the layer parallel shortening itself, obtained by rather scarce absolute dating [39], gives an indirect validation to SRIT. The potential of this tool, which is seemingly independent from temperature and fluid pressure, calls for a direct demonstration of its validity. Here, we identify an opportunity to validate the SRIT method using independent depth calculations. The depth derived from palaeothermometry can be associated with calcite precipitated during the pressure solution event, with a well-constrained geothermal gradient. When comparing the distribution of the population of depth from SRIT (n = 27) to the population of depths derived from δ¹⁸O (n = 27) and T₄₇ (n = 3) (Figure 7a), it appears clearly that there is a strong overlap between both results. The depth reconstruction from stylolite roughness inversion is slightly higher than the depth derived from cement temperature, but it might reflect—assuming the shift is significant—the fact that SRIT records (near) the end of the activity of the stylolite, while jogs develop during the stylolite growth [48]. However, direct comparison of SRIT on a given stylolite and precipitation temperature of the calcite at its tip (Figure 7b) show that there is no good correlation at the individual scale. Indeed, two-thirds of the data show discrepancy between both reconstructed depths that is higher than the methodological uncertainty. Such behaviour can be explained when the temperature of precipitation corresponds to a slightly shallower depth than the depth of cessation of stylolite development, as the redistribution might be occurring during the pressure solution and not after. Nevertheless, this illustrates that this method is more reliable when using a population of stylolite, and that limiting its application to one stylolite must be used only as an indicator of the order of magnitude of the vertical stress. On top of this, it is important to bear in mind the sources of uncertainty in reconstructing depth with both approaches (elastic parameters for SRIT, geothermal gradient and surface temperature for T₄₇). Consequently, the first order distribution of SRIT appears entirely consistent with the expected paleodepth at which stylolites were active. These results demonstrate that SRIT can be used as a reliable vertical stress gauge even when applied to a relatively small number of stylolites sampled from the same bedding following the protocol developed in [39,41].

6.2. Timing of Stylolite Development Regarding Burial History and Favourable p,T Conditions

In the Calcare Massiccio formation in the Umbria Marche sequence, results from both thermometry and SRIT (Figure 7) suggest a distribution of depths ranging from 600 m to ca. 1900 m for the activity of pressure solution. It is interesting to note that the maximum burial depth of the top of the Massiccio Fm. is estimated to vary towards the inner parts of the belt (eastward) from 2400 m to 3600 m depth (Figure 1D, [39]). Hence, the reconstructed depth values are not the maximum depth that was likely experienced by the strata. This raises the question of the timing of BPS development with respect to the burial history of a strata, along with the timing of jug development associated with the stylolite growth. In this study, we illustrate that pressure solution is a major deformation mechanism below 2000 m depth (Figure 8), but the shallowest stylolite development is just below 600 m, well after the induration of sediments and, in that case, posterior to the dolomitization of the ooids. Indeed, the rock needs to be coherent mechanically to develop the jogs, and the filling is not affected by the extensive dolomitization. In terms of absolute timing, the reconstructed depth range [600–1900 m] is very similar to the ones reconstructed elsewhere in the belt from shallower strata (mainly the Maiolica and Scaglia Rossa Fms.), ranging from 700 to 1600 m in the centre of the UMAR (San Vicino, Cingoli) and being a bit deeper in the outermost and innermost folds ([1100–2300 m] in the Conero Anticline and [1100–2900 m] in the Subasio Anticline, respectively). In terms of absolute timing, if we consider the burial model of the central area of the UMAR, where most quarries can be found (Figure 1A), the range of depth reconstructed from the slab predict a timing of active pressure solution spanning from 170 My to roughly 20 My (15 My if considering the 12% uncertainty), so a total (continuous?) duration of 150 Ma. That would correspond a very long time during which the maximum principal stress remained mostly vertical. In any case, it is worth noting that the end of burial-related BPS development in the slab is in line with the regional halt estimated around 10 My ago [39].

For the first time, the investigation of jogs developed synchronously to pressure solution along with the study of the roughness of the stylolites allow the temperature and pressure condition at which stylolites can form in carbonates to be refined. In the present case, pressure solution seems to be more active between 15 and 45 MPa and between 25 and 55 °C. This range is in line with previously published SRIT studies where most of the sedimentary stylolites developed between 600 and 2000 m [31,33,34,36,37,38,39,40,41]. However, that requires further investigation to confirm this trend, as this range might be too limited with regards to some previous studies that reconstructed a higher depth of pressure solution with SRIT [32]. This study further supports the idea that pressure solution formation progresses gradually with depth, with the cessation of activity on one stylolite occurring at a different depth than on another within the same strata [31,33]. The reconstructed pressure solution activity depth range is similar to the clogging of the vuggy porosity, as supported by the isotopic geochemistry of the vugs, wherein the precipitation temperature has been measured in one vug to be 46.5 °C and calculated in the other one to be 47 °C. However, although there is not enough data to check whether these calcite crystals directly reprecipitated from the dissolved material along the stylolite, it is safe to state that this clogging is synchronous to the activity of the pressure solution. Indeed, the development of the available space (i.e., the jog) is related to the local volume reduction around the stylolite tips and is unlikely to create a network allowing large-scale fluid flow. In other words, whether the calcite precipitated from the fluid produced by the pressure solution or whether it precipitated from local fluids contained in the pores, it can safely be considered to be at thermal equilibrium with the host rock at the time pressure solution occurred. All in all, this range of depth might represent the depths at which local porosity can be clogged and so permeability destroyed by mineral redistribution in reservoirs, key information in reservoir modelling.

7. Conclusions

The development of mixed-mode (opening and shear) jogs at stylolite tips during pressure solution, associated with coeval calcite precipitation in these cracks, provides a unique opportunity to conduct a critical assessment of applying SRIT to a population of BPSs as a reliable vertical stress/depth gauge and to further explore the p-T conditions for pressure solution. In the present study, depths derived from SRIT are compared to calcite crystallization temperature estimates from clumped isotopes and extended to O isotopes on the calcite filling the jog. This study presents a novel comparison to another reliable marker on the exact same object. The results clearly indicate that depths reconstructed from SRIT are in the same range as those inferred from the coeval precipitation temperature of calcite in jogs, ranging from approximately 600 m to 1900 m based on temperature and from approximately 600 m to 1950 m based on SRIT, with SRIT yielding a slightly higher median value (1350 m vs. 1150 m). This consistency validates the application of SRIT to a population of BPSs of type 3 morphology for reliably reconstructing the range of depths at which pressure solution was active and prominent during the burial history.