Paleo‐thermal constraints on the origin of native diagenetic sulfur in the Messinian evaporites: The Northern Apennines foreland basin case study (Italy)

Recent studies on the genesis of sedimentary native sulfur deposits indicate diagenetic mid‐low temperature Bacterial Sulfate Reduction (BSR) as the main process, involving organic compounds (kerogen/hydrocarbons), bacterial colonies and gypsiferous rocks. In the peri‐Mediterranean area (Southern Spain, Sicily, Northern Apennines, Israel), the main sulfur accumulations are always associated with late Miocene sulfates and organic‐rich successions encompassing the Messinian salinity crisis (MSC). In particular, the Messinian successions of the Apennine‐Adriatic foreland basin system, due to a large amount of high‐resolution stratigraphic data, represent a perfect case study for understanding the diagenetic conditions controlling the development of the BSR process during sedimentary basin evolution. In this work, thermal models performed in three sub‐basins in a sector of the Northern Apennines comprised of the Sillaro and Marecchia rivers (Italy), calibrated by means of organic and inorganic geothermometers, indicate a general thermal immaturity of the studied successions attained as a result of a constant heat flow similar to the present day one (ca. 40 mW/m2) since Late Tortonian and lithostatic loads between 615 and 1,710 m depending on different sub‐basins. These results suggest that the MSC deposits experienced maximum temperatures between about 39°C and 65°C. Temperatures derived from thermal models have been used to constraint occurrence of the diagenetic BSR associated with evaporitic deposits providing thermal constraints in sulfur genesis as well as new useful thermal‐constraints for basin analysis studies.

The sector of the Northern Apennines comprised between the Sillaro and Marecchia rivers (Italy) is a key area to understand the mechanism and timing of BSR-related sulfur formation that occurred in MSC evaporites in different depositional settings (Figures 2 and 3) of the Apennine-Adriatic foreland basin system. The Vena del Gesso wedge-top marginal basin to the west, where primary shallow-water evaporitic deposits accumulated, is characterised by very limited sulfur occurrence. Conversely, the inner foredeep Romagna sub-basins to the east contain clastic evaporites (turbidites and mass-wasting deposits) and show the largest sulfur deposits in the Apennines.
In this work, we adopted a multidisciplinary approach in order to reconstruct the thermal conditions at which the native sulfur deposition developed, combining thermal data from classical vitrinite analyses with those derived from the gypsum-anhydrite transition. Our findings shed new lights on the temperature and timing conditions necessary to generate diagenetic sulfur that may help similar studies in other worldwide areas.

| THE MESSINIAN SALINIT Y CRISIS STRATIGRAPHY AND THE BSR SULFUR DEPOSITS
Between 5.97 and 5.33 Ma Van Couvering et al., 2000), the Mediterranean Sea experienced one of the most catastrophic events in its recent geological history, known as the Messinian salinity crisis (MSC; Hsü et al., 1973;Selli, 1973), leading to widespread evaporite 2502 | EAGE ROSSI et al.

STRATIGRAPHY AND EVOLUTION OF THE STUDIED AREA
Our work is focused on three (sub-)basins of the Apennine-Adriatic foreland basin system showing different geological settings (Figures 2 and 3 Manzi et al., 2005Manzi et al., , 2007Roveri et al., 1998Roveri et al., , 2001Roveri et al., , 2003 and references therein). These (sub-)basins are located along the north-eastern flank of the Apennine belt and their tectono-sedimentary evolution was strictly related to the post middle Miocene tectonic events (Rossi et al., 2002(Rossi et al., , 2015Roveri et al., 1998Roveri et al., , 2001Roveri et al., , 2003Thomson et al., 2010;Zattin et al., 2002), which differently developed in adjacent sectors according to the timing of the local growth of individual morphostructural elements (Ricci Lucchi, 1973, 1981, 1986).
through gravity flows in the adjacent deeper basins. In the Vena del Gesso basin, the whole PLG unit was involved in the sliding (Roveri et al., 2003). The PLG is sealed by clastic deposits (Cusercoli Fm.) that belong to the third stage (5.53-5.33 Ma) and show an extremely reduced thickness (a few tens of metres) accumulated in a time span of about 200 ka; Manzi et al., 2007;Roveri et al., 1998Roveri et al., , 2001Roveri et al., , 2003Roveri & Manzi, 2006). Normal marine conditions, suggested by microfossil-rich hemipelagic deposits of the Argille Azzurre Fm., are observed in the basal Zanclean (Pliocene, 5.33 Ma; Hsü et al., 1973;Manzi et al., 2005;Roveri et al., 2003). This phase is characterised by tectonic quiescence lasting till the MPL3-MPL4 biozones boundary (4 Ma) after which a new tectonic phase, related to Apennines deep-rooted movements, caused the north-eastward tilting of the whole Apennine belt (Roveri et al., 2003;Thomson et al., 2010;Zattin et al., 2002). The sedimentation remained dominated by hemipelagic settling until the middle Pleistocene (0.781 Ma), when the onset of coastal deposits took place (Sabbie di Imola Fm.; Colalongo et al., 1982;Cremonini et al., 1969). The tectonic emplacement of the Ligurian unit (Landuzzi, 1994), heralded by submarine landslides (Val Sellustra Olistostrome Fm.;CARG, 1988), including calcareous olistoliths in a shaley matrix (Castellarin & Pini, 1989), occurred between ca. 4.5 Ma and 0.781 Ma (CARG, 1988;Colalongo et al., 1982) only in the westernmost portion of the VdG basin. Since 4.0 Ma the whole VdG basin was continuously uplifted up to today (Roveri et al., 2003;Thomson et al., 2010;Zattin et al., 2002). The described geological history of the area is simplified in order to provide an overall description of the accumulation/loss of burdening sediment, and thus, the overall increase/decrease of load (temperature); please see Roveri et al. (2003Roveri et al. ( , 2014, Rossi and Rogledi (1988), Rossi et al. (2002Rossi et al. ( , 2015 and references therein for further details. In this area, evidence of gypsum-anhydrite transformation in the PLG unit are generally lacking except in its westernmost portion, close to the Sillaro thrust lateral ramp of the Ligurian nappe. Here, native sulfur is present in small quantities (Figures 4 and 5).

| Giaggiolo-Cella and Sapigno subbasins
The GCs and Ss sub-basins (Figures 3-5) are separated from the VdG basin by the Forlì thrust front -Riolo anticline (FTF in Figure 3) and represent an inner foredeep depositional setting (Manzi et al., 2005;Roveri et al., 2001Roveri et al., , 2003. Both subbasins share a similar stratigraphy (Figure 4). Starting from the late Tortonian, the succession is characterised by the turbiditic lobes of the Marnoso-arenacea Fm. (Fontanelice Mb.), followed by thin-bedded turbidites including submarine slides (Ghioli di letto Fm.; Lucente et al., 2002;Manzi et al., 2005Manzi et al., , 2007Roveri et al., 1998Roveri et al., , 2001Roveri et al., , 2003. The upper portion of this Formation is composed by a 60 m-thick interval of organic-rich shale encompassing the Tortonian-Messinian boundary (7.24 Ma) and the entire early Messinian. The MSC onset (5.971 Ma, Manzi et al., 2013) occurs in the uppermost part of this unit and is marked by the disappearance of foraminifera. This is the Foraminifera Barren Interval unit (FBI; Manzi et al., 2018;Gessoso-solfifera Fm.) representing the deep-water time equivalent of the PLG evaporites (Figures 2 and 4) with scattered or concentrated dolomite and a significant organic matter content (TOC up to 5%; Manzi et al., 2007;Roveri et al., 2016) indicative of strong anoxia and water stratification . The FBI unit is sharply overlain by clastic resedimented evaporites composed of gypsum-bearing turbidites and olistostromes deriving from the dismantlement of the PLG unit (RLG unit, Resedimented Lower Gypsum; Gessoso-solfifera Fm.; Manzi, 2001;Manzi et al., 2005Manzi et al., , 2007Parea & Ricci Lucchi, 1972;Ricci Lucchi, 1973;Roveri et al., 2001Roveri et al., , 2003; see also Schlager & Bolz, 1977). The RLG unit shows its maximum thickness in the Ss sub-basin (Manzi, 2001;Manzi et al., 2005Manzi et al., , 2007 and is floored by an unconformity surface (MES) or its correlative conformity (MES-CC). The succession of the third stage (5.54-5.33 Ma, Bassetti, 2000;Bassetti et al., 1994;Roveri et al., 1998Roveri et al., , 2001Roveri et al., , 2003, is characterised by a greater thickness (up to 700 m) with respect to that of the VdG and can be subdivided into two Formations (Roveri et al., 1998(Roveri et al., , 2001(Roveri et al., , 2003, separated by a minor unconformity. The first formation (Tetto Fm., 5.54-5.42 Ma) is composed of a succession of monotonous shales, interrupted by rare sandstone horizons, and containing a thin volcanoclastic horizon dated at 5.53 Ma (Cosentino et al., 2013;Odin et al., 1997). Differently, the overlying Cusercoli Fm. (5.42-5.33 Ma) is characterised by the alternation of 3 coarse-(conglomerate and pebbly-sandstone fluvio deltaic bodies) and fine-grained (shale and fine sandstone shelf deposits) couplets. Three lacustrine limestone beds (Colombacci) and two paleosoils are also present and the fauna is represented by the typical fresh-water paratethian species found in the Lago-Mare interval (see Bassetti et al., 2004). As in the VdG basin, the return of normal marine conditions at 5.33 Ma is marked by hemipelagic sediments (Argille Azzurre Fm.) that continued in the GCs sub-basin until the Middle Pleistocene (0.781 Ma). Conversely, in the Ss sub-basin the deposition was interrupted by the emplacement of the Ligurian sheet (Val Marecchia thrust lateral ramp; Elter & Trevisan, 1973;Merla, 1951), capped by the sandy-conglomeratic fluviodeltaic episutural deposit of the Perticara Sandstone Fm. (De Feyter, 1991;Manzi, 2001;Manzi et al., 2005Manzi et al., , 2007Zattin et al., 2002). The Ss sub-basin loading by the Ligurian units ended in the Middle Pleistocene (0.781 Ma). Since about 4.0 Ma both the Ss and GCs sub-basins experienced a phase of general uplift affecting the whole Northern Apennine mountain range and were subaerially exposed since the Middle Pleistocene (0.781 Ma; Roveri et al., 2003;Thomson et al., 2010;Zattin et al., 2002; for further details). In the two sub-basins, gypsum was transformed into anhydrite during burial and was rehydrated back to gypsum after exhumation by meteoric waters.

| MATERIALS AND METHODS
In the VdG basin, we reconstructed five stratigraphic logs spanning from the Late Tortonian to the Pliocene and located in shallowing paleo-settings from the west to the east, progressively closer to the Riolo paleo-high (Figures 3 and 5): Sassatello (S), Fontanelice (F), Borgo Tossignano (BT), Monte Tondo (MT) and Brisighella (B). In the GCs sub-basin, we measured two stratigraphic logs at Giaggiolo (G) and Campea (C) (Figures 3 and 5), including a succession spanning from the Late Tortonian to the Pliocene and located, respectively, in a depocentral area and in a paleostructural high ( Figure 5). In the Ss sub-basin, we measured a single section along the Fanantello stream (FS, Figures 3 and  5) spanning from the Late Tortonian to the Pliocene and located in the core of a syncline dipping south-eastward below the Ligurian Sheet units.
A total of 41 organic-rich fine sediments (shales and marls) samples have been collected: 27 derived from the VdG, six from the GCs and eight from the Ss sub-basins ( Figure 5). For each sample, we carried out vitrinite reflectance analyses following the standard methods for sample preparation and measurement (Bustin et al., 1990;Corrado et al., 2019;Lucca et al., 2018;Schito et al., 2016). Results are summarised in Table 1, Figure S1 and Table S3. The eight stratigraphic sections have been modelled using the software BasinMod2D® (Platte River Associates, Inc) based on Sweeney and Burnham's (1990) kinetics model. A general overview on the Vitrinite Reflectance method and Thermal Modelling methods, as well as more information on model inputs are described in the supplementary notes (please see also Figure S2).
For thermal maturity calibration, vitrinite reflectance data were used. Moreover, the presence or not of re-hydrated gypsum that experienced gypsum-anhydrite transformation was used as a further temperature constraint for the base of evaporitic successions. As matter of fact, petrographic data suggest that gypsum in the western part of VgG basin and in Sapigno and Giaggiolo basins was turned into anhydrite and then re-hydrated back to microcrystalline gypsum (alabastrine) at surface exposure by meteoric waters (Murray, 1964). The presence/absence of alabastrine gypsum allows to define if the original evaporite deposit were/were not heated over the transition temperature, which has been experimentally determined at 51.9°C (Jaworska, 2012;Jowett et al., 1993;Klimchouk, 1996;Murray, 1964).
Assumptions used to run modelling are: (a) decompaction of the burial curves according to Sclater and Christie (1980); (b) thrusting duration is considered instantaneous (Endignoux et al., 1990); (c) a sediment-water-interface temperature of 11°C. Please see Figure S2 for further information.  Figure S1). In general, from 20 to 50 measurements were carried out on each samples. Vitrinite is the most common maceral and generally shows a small dimension (<10 µm) and different shapes varying from elongated to stocky and from square to slightly rounded. The fragments preservation state is good, although oxidation can occur. Shaly lithologies such as Argille Azzurre or some levels in the Ligurian Units show a high content in amorphous organic matter (AOM), while inertinite is found amongst all over the samples. Reworked organic matter with higher reflectance values was also recognised (blue bars in histograms in Figure S1). Sulfur crystals and sulfides, in particular rhombohedral pyrite, are widespread. Vitrinite reflectance values tend to increase towards the deepest portion (Marnoso-arenacea Fm. samples) in all eight logs.

| Thermal model calibration
In the VdG basin (F, BT, MT and B logs), thermal model calibration curves match the measured thermal maturity data using a constant heat flow (HF) value of 40 mW/m 2 and an amount of eroded burial of Argille Azzurre Fm. estimated at around 500 m ( Figure 6). In all tested scenarios, the temperatures reached by the gypsum and related ORS units are never higher than 44°C (see Table 2).
The S log represents the westernmost area of the VdG basin. Here, thermal model calibration can fit R 0 % and T gat constraints, with a constant heat flow of 40 mW/m 2 and an eroded burial of about 1,100 m belonging almost exclusively to the Val Sellustra unit of the Ligurian thrust sheet ( Figure 6). Paleotemperatures for the PLG and ORS units overtook T gat (53-55°C) and endured above it for minimum of 20 ka (see Table 2).
In the GCs sub-basin (G and C logs), thermal model calibration is validated using an HF of 40 mW/m 2 and an eroded burial (Argille Azzurre Fm.) of about 470 m ( Figure 6). Paleo-temperatures for the RLG unit rests above T gat for minimum of 60 ka showing values within 54°C and 64°C, while the FBI unit reached the maximum temperature of 65°C (see Table 2). Native sulfur deposits developed along the northeastern flank of the basin (Borgopaglia anticline; Figure 5).
For log FS (Ss sub-basin), its thermal calibration matches measured R 0 % and T gat constraints with a constant heat flow of 40 mW/m 2 and an eroded burial of about 1,200 m belonging to Ligurian and Epiligurian units ( Figure 6). Paleotemperatures resided above T gat but for a minimum of 110 Ka, maximum temperature reached by the FBI and RLG units is 65°C and 57°C, respectively (see Table 2). Extensive native sulfur generates alongside the north-eastern flank of the subbasin (sulfiferous limestone).
Using these data in thermal model calibration we considered possible solutions provided by the interplay between HF and eroded burial variation for all the studied sections ( Figure 6), reaching the conclusions that a constant HF from late Tortonian to present day of 40 mW/m 2 for the whole area is the most reliable solution since it best fits with the geological history of the area. This value agrees with the extrapolation of Della Vedova et al. (1995Vedova et al. ( , 2001, with the surface heat flow maps of Italy proposed by CNR in 1994 and is further confirmed by a more recent study by Pauselli et al. (2019) based on new available thermal data from 174 wells in the Northern Apennine. Using such HF values, calculated the thicknesses of eroded burial of both Ligurian and Argille Azzurre units are quite consistent amongst all sections. In detail, calculated thicknesses of overthrusted Ligurian units' range between 1,100 m in the Sassatello section to 1,200 m in the Fanantello (Ss) ones. These values are consistent with Corrado et al. (2010) and Thomson et al. (2010), that, in the same area, found that the AHe (Apatite-Helium geothermometers) system has been not reset and surface Ro% data are generally lower than 0.5% implying surface sediments did not experience temperatures higher than about 60°C-80°C. As well as the reconstruction of the total thicknesses of Argille Azzurre Fm. before erosion never exceed 600m that corresponds to the maximum present-day thicknesses measured in the nearby areas by Conti et al. (2016).
According to thermal models, temperatures experienced by FBI/ORS at the base of the evaporite units in the Ss, GCs and westernmost portion of the VdG (sub-)basins experienced temperatures between 55°C and 65°C (Figure 7), being enough to start organic matter biodegradation (i.e. maturation according to Bjørlykke, 2010;Tissot & Welte, 1984) and, consequently, to develop organic compounds necessary for the BSR process (Machel, 2001 and references therein). Conversely, in the VdG basin where sulfur did not develope, all tested scenarios indicate that the temperatures reached by the ORS were never higher than 44°C.
These results suggest that native sulfur can form from bacterial reduction in sedimentary basins for temperatures higher than 55°C and not lower than 44°C. Notably, the sulfur found in the Perticara Mine (Ss sub-basin) is probably related to the presence of migrated oil deriving from deposits buried in the depocenter of the sub-basin as suggested by Roveri et al. (2016), thus, not only dependent on maximum temperatures recorded in the FS succession.
Likewise, the gypsum-anhydrite transformation developed in high burial condition, thus when temperatures were higher than 53°C (GCs, Ss and western VdG sub-basins) and it never occurs below 44°C (central and eastern VdG basin).
The temperature boundaries of the gypsum-anhydrite transformation have been usually derived from experimental  Note: Results of thermal modelling calibration using HF=40 mW/m 2 . "Preserved" may include Tetto, Cusercoli, Argille Azzurre Fm.s and Ligurian sheet unit depending on log (see Table S1); "eroded" includes Argille Azzurre Fm. or Ligurian sheet unit (see Table S1).
Abbreviations works while little constraints regarding its occurrence in geological environment are available. According to Jowett et al. (1993), Hardie (1967) and Bredehoeft and Hanshaw (1968), the transition is mainly driven by temperature and by the activity of water. In the presence of coeval NaCl precipitation, the transition occurs at temperatures of 18°C. Nevertheless, according to Jowett et al. (1993), if there is no evidence of salt precipitation and the surrounding lithologies are siliciclastic (sandstone or shales) it can be assumed that pore water reflects the original formation water. The activity of water at which gypsum precipitates is 0.93 and assuming such a number the transition occurs at 51.9°C according to Hardies's (1967) experimental data. In our work, through temperature assessment at the base of the evaporitic successions by mean of the thermal model, we found that in this case history gypsum-anhydrite transition never occurred below 53°C thus in agreement with Jowett et al. (1993) observations. Moreover, we can set a lower limit suggesting that in similar geologic condition (e.g. lithologies and fluids composition), the transition never occurs below 44°C. We also suggest that in a geological environment, the process seems to be very rapid but not strictly time-dependent since it occurred in a time span ranging from 20 to 270 ka. The gypsum-anhydrite transition released a significant amount of crystallisation water, which may carry significant amounts of dissolved sulfate in contact with the most mature organic-rich shale interfaces, enhancing the BSR process in the Ss, GCs and westernmost portion of the VdG (sub-)basins. Thus, in this work, we suggest a genetic model for sulfur by BSR process in a diagenetic environment that involve the gypsum-anhydrite transformation from original gypsum facies (primary, PLG, or clastic, RLG) and the presence of slight immature organic matter in a temperature window from around 44°C to about 80°C. The lower boundary is suggested by the absence of sulfur in the eastern part of the VdG basin, while the upper limit is the temperature above which sulfate-reducers bacteria do not survive (Jorgensen et al., 2019;Machel, 2001). The association between gypsiferous rocks, presence of organic matter and native sulfur generation was already suggested in the literature by many authors (please see Ziegenbalg et al., 2010; and references therein) on the base of isotopic data. Nevertheless, the great novelty of this work is that we constraint the temperature range in which the process occurs in a natural case history.
The defined temperature window for the BSR-related sulfur genesis may be used for basin analysis purposes, where sulfur of sedimentary origin is present, thus in the other peri-Mediterranean basins and in other worldwide settings containing evaporites. Moreover, the evidence that gypsumanhydrite transformation and BSR native sulfur formation took place in the same areas is noted even if a certain relationship amongst these two processes cannot be definitely stated here: further analyses are needed to clarify the relationships of the sulfur genesis with the primary gypsum facies and the supposed temperature-related sulfate transition.
Furthermore, from a geological point of view, the results obtained by the reconstruction of the burial and exhumation history of the three (sub-)basins seem in contrast with the idea of a possible paleo-continuity amongst the Sillaro and Val Marecchia tectonic structures as proposed by some authors F I G U R E 7 Temperatures evolution across different sections studied in this work based on the simplified structural map (cf. with the simplified geological map of Figure 3) with temperatures indicated in (see Zattin et al., 2002). Conversely, it can be confirmed the hypothesis that the Ligurian sheet did not develop as a single thrust front, but its allochthonous Northeast-ward movement was hindered by the innermost thrusts involving the oldest portion of the Marnoso-arenacea clastic fill between the Sillaro and Marecchia valleys (Lucente et al., 2002).

| CONCLUSIONS
Paleo-thermal data obtained from the thermal modelling of three (sub-)basins of the Northern Apennines provided qualitative and quantitative information on the triggering conditions needed to form native sulfur from evaporites through BSR in both local and regional contexts. A heat flow of 40 mW/m 2 active since the late Tortonian, coupled with a variable lithostatic load produced by about 1,080 to 1,710 m of sediments, generated paleo-temperatures compatible with BSR from 55°C to 64°C at the ORS/FBI and PLG/RLG interfaces. Moreover, considering the minimum paleo-temperature (55°C) reached by the ORS/FBI interface, where native sulfur developed (Giaggiolo-Cella, Sapigno and westernmost portion of Vena del Gesso sub-basins), these conditions are enough for organic matter degradation and generation of organic compounds necessary for BSR. A paleo-temperature slightly lower than 44°C has been reached by the ORS interval underlying the evaporitic unit in the area where sulfur deposits did not develop (eastern Vena del Gesso basin), suggesting that this may be the lower limit for BSR. The evaporites that generated native sulfur also underwent the gypsum-anhydrite transformation, and thus remained over the transition temperature (52°C) for minimum of 20 ka, as shown in the S log example. Hence, the process of anhydritisation appears to be very fast and can release crystallisation water that may provide a significant amount of dissolved sulfate, facilitating the BSR. In summary, the association of BSR-related native sulfur, organic matter maturation and secondary (alabastrine after anhydrite) gypsum deposits define a specific temperature window between about 44°C and 80°C. This temperature interval may be used as paleo-temperatures constraints in the thermal modelling of basins where evaporite and organic-rich sediments were involved in burial-exhumations processes (foredeep-foreland system development). Finally, thermal modelling results obtained in this work also strengthens the idea that both Sillaro and Val Marecchia tectonic elements may be correctly considered as thrust lateral ramps, thus, Ligurian thrust sheet did not develop as a single front.