Are fluid inclusions in gypsum reliable paleoenvironmental indicators? An assessment of the evidence from the Messinian evaporites

Fluid inclusions (FIs) in marine minerals such as calcite, gypsum, and halite are a very useful tool for the understanding of the depositional environment because they represent microsamples of ancient seawater from which the minerals precipitated. This information can be obtained by microthermometric analyses unless FIs were modified after trapping. This approach has been applied to the primary bottom-grown selenite gypsum accumulated in the Mediterranean Sea during stages 1 and 3 of the Messinian salinity crisis (MSC; CIESM, 2008; Roveri et al., 2014a; Fig. S1 in the Supplemental Material¹). However, the salinity of the water body from which the evaporites precipitated is still debated. Several FI studies (Attia et al., 2004; Natalicchio et al., 2014; Evans et al., 2015; Costanzo et al., 2019) led to the idea that the Messinian evaporites precipitated from low-salinity waters, well below the gypsum precipitation field (GPF) for normal marine evaporites. The hypotheses explaining the origin of these low-salinity waters are complex. Natalicchio et al. (2014) suggested that the source of ions came from dissolution of preexisting gypsum by continental fresh water, but no evidence of previous marginal gypsum deposits has ever been found. Grothe et al. (2020) theorized that the Mediterranean-Paratethys water exchanges at the beginning of the MSC caused the formation of a low-salinity surface layer rich in calcium and sulfate, triggering gypsum precipitation at concentrations as low as 40 wt‰ NaCl equivalent (eq.). Also in this case, the source for the solutes remains unexplained.

Contrasting evidence for the continental origin of the evaporites comes from the isotope geochemistry and the biological content, which point to a significant marine signature (Roveri et al., 2014b). From the MSC onset, the Mediterranean began to develop its own hydrology due to the reduced exchanges with the Atlantic Ocean, but during the arid periods at (sub-)precessional scale, the ⁸⁷Sr/⁸⁶Sr returned to equilibrium with the global ocean (Reghizzi et al., 2018). Only at stage 3 of the MSC is a signal clearly distinct from the global ocean recorded by evaporites, limestone, molluscs, and ostracods (Fig. S2).

The δ³⁴S and δ¹⁸O from MSC stage 1 gypsum fall very close to those of Miocene seawater, suggesting precipitation from a mostly marine water body (Müller and Mueller, 1991; Lu and Meyers, 2003; Evans et al., 2015; García-Veigas et al., 2018). These two proxies recorded the periodic marine ingressions in the initial (Lugli et al., 2007) and final stages of the crisis (Longinelli, 1979; Manzi et al., 2009), outlining precessional arid-wet cycles.

In this discussion about the salinity of the parent brine from which gypsum precipitated, some other important aspects appear to be commonly overlooked:

1. Marine faunal assemblages are found in the shale interbedded in the stage 1 gypsum deposits (Neraudeau et al., 2002; Roveri et al., 2020).

2. Clear evidence of a Paratethyan water influx is found only in the stage 3 deposits with the development of the “Lago-Mare” brackish-water fauna (Roveri et al., 2008, 2014a).

3. Huge volumes of evaporites from stage 1 and especially stage 2 (halite) suggest that the only viable source of ions could have been the Atlantic Ocean.

4. The δ³⁴S of the gypsum of stage 1 commonly points to a marine origin.

Our study explores a new methodology for assessing the reliability of previous FI analyses. We systematically measured FIs never described before within gypsum crystals. Our revised methodological approach demonstrates the existence of high-salinity FIs that were not previously reported.

We analyzed 12 samples from 9 stratigraphic sections in different geological contexts throughout the entire Mediterranean Basin from the bottom-grown primary gypsum accumulated during MSC stages 1 and 3 (Fig. S3). Ten samples are from the Primary Lower Gypsum (PLG; stage 1; 5.97–5.60 Ma; Lugli et al., 2010; Manzi et al., 2013) of Spain (Sorbas and Almeria-Nijar Basins), northern Italy (Piedmont, Vena del Gesso, and Adriatic Basins), Sicily (Caltanissetta Basin), Egypt (Al-Barqan area), and Cyprus (Pissouri Basin). Two samples are from the Upper Gypsum (UG; stage 3; 5.55–5.33 Ma; Manzi et al., 2009) collected in Sicily (Caltanissetta Basin) and Cyprus (Tokhni Basin).

We distinguished two primary time-equivalent FI types aligned along different growth bands of the gypsum crystals merging at the oblique boundary between the dark core and the clear portion of the crystal (Figs. 1A, 1B, and 1F):

1. Type A: Pyramidal-shaped FIs 10–200 µm across, aligned along the crystal growth surface within the reentrant angle of the swallowtail twin (Fig. 1C), which were analyzed in previous studies.

2. Type B: Mainly tabular, hexagonal FIs 5–30 µm across (Figs. 1D and 1F) and minor pyramidal and triangular inclusion 5–200 µm across, aligned along the vertical-oriented growth bands parallel to the twin plane; these FIs have not been described in the literature to date.

A total of 55 millimeter-sized fragments were obtained by cleaving crystals along the 010 cleavage plane with a razor blade. Microthermometry was conducted on a total of 593 FIs to observe the last ice melting temperature, T_m(ice), from primary FIs within lateral growth bands and in the reentrant angle of the twin. The analytical method differs significantly from those of previous studies (Attia et al., 2004; Natalicchio et al., 2014; Evans et al., 2015; Costanzo et al., 2019) that allowed the measurement of only a small portion of the FI population. With respect to these previous works, we ran temperature cycles with larger excursions (from −100 °C to +120 °C) at a faster rate (50 °C/min) to enhance the mechanical stretch to the nucleation of a stable bubble. Performing as many as nine cycles, we dramatically increased the number of measurable FIs compared to previous method, using a stress cycle from −90 °C to +30 °C at rate of 30 °C/min (e.g., Attia et al., 1995). This is because the presence of the bubble is the prerequisite for the microthermometric measurements, otherwise the phase changes cannot be observed (Roedder, 1984). More details on our procedure are provided in the Supplemental Material. After the extreme stretching procedure, we performed the measurement run following the same criterion described by Attia et al. (1995).

The salinity values obtained are shown in Figure 2 and in Tables S2 and S3 (in the Supplemental Material). Overall, the salinities vary from 9‰ to 235‰ (wt‰ NaCl eq.) corresponding to T_m(ice) from −0.5 °C to −21.7 °C. Nine samples have average salinity values falling within the GPF and ranging from 123‰ to 160‰. Only three samples (RA16-3, SE11-15, and PISS-2) have average salinities (93‰, 94‰, and 85‰, respectively) below the gypsum saturation point. All samples show a wide range of values (Fig. 2; Figs. S5 and S6). Sample EMOG-13 shows the smallest salinity range (108‰–195‰; Figs. 2B and 4), with all FIs except one having marine values and with the highest frequency around 150‰. Except in sample EMOG-13, the values are concentrated in multiple high-frequency points on the normal distribution curve (Fig. 2B; Fig. S6). However, in samples AL16-4, BQ1+BQ3, MT14_4-3-10, SAB6, and TO-1, two main frequency peaks are clearly distinguished, one within the GPF and the other below (Fig. 4; Fig. S6).

Our results from FIs, both the new type (type B) and the type usually measured (type A), show a wide salinity range not matched by previous studies that provided only very low salinity values. The main explanation until now for these low values has been that gypsum precipitated from a water body with a strong mixing of seawater and fresh water enriched in Ca²⁺ and SO₄²⁻. An alternative process could be an unusual precipitation of gypsum triggered by microbiological activity (Natalicchio et al., 2014). An unexplored, more simple explanation is that FIs were modified after trapping. As shown in Figure 2A, the low-salinity FIs are statistically larger, and many of them show clear signs of modification due to post-trapping processes (leakage, necking-down, recrystallization, etc.; Roedder, 1984), whereas high-salinity FIs are statistically smaller. Fluid inclusions with salinity below the gypsum saturation point (110‰) range in size from 12 to ~12,000 µm², whereas those falling in the GPF are much smaller, from 12 to 720 µm²; however, 99.7% of FIs within the GPF are <600 µm² (Fig. 3). The FIs falling within the GPF actually represent the original parent brine, while those with lower salinity have probably been modified by post-depositional processes. This is because gypsum, unlike other evaporitic minerals, presents a perfect cleavage along the 010 plane, which makes its crystals structurally weak by load-unload cycles or tectonic stress. This weakness allows multiple crack-and-seal cycles, which can promote the introduction of low-salinity late fluids. Larger FIs have a higher probability of being intercepted by crack-and-seal processes along 010 planes than the small ones and thus can be more easily involved in mixing by secondary fluids. The result is that we find low salinities especially in large FIs that show clear signs of modification such as partial dissolution and recrystallization (highlighted in red in Fig. 2). The size of 720 µm² (Fig. 3) appears to be the boundary beyond which it is no longer possible to find salinities falling within the expected range of normal seawater precipitation. The size boundary is probably controlled by the spacing of the crack-and-seal discontinuities. One of the major indicators of post-depositional modification is the presence of FIs with a large salinity range (~10‰–150‰) within a single growth band (e.g., Fig. 1D). There would be no reason for FIs along the same growth band to display different salinities unless they were modified in different proportion according to their size by crack-and-seal processes. The crack-and-seal mechanism may not be identified under the microscope because of the strong tendency of gypsum to grow syntaxially without leaving visible traces.

A comparison between previous studies and this work (Fig. 4) shows how the new data are distributed over a wider salinity range reaching higher values. However, the frequency is not entirely random. Especially in the case of the Piedmont Basin (Fig. 4A), two frequency peaks are clearly recognizable. The peak at lower salinities, below the gypsum saturation point, partially overlaps with the previous data set, while that at the highest salinities falls entirely within the GPF and covers salinities never obtained before. The largest FIs belong to the first peak and could represent the preferentially modified ones, whereas the second peak consists only of the smallest FIs, which preserved the parent brine. The methodology used in the previous works allows measurement of only the low-salinity FIs and therefore does not provide representative data. A more representative methodology is only possible by applying a more intense mechanical stress that can allow enlargement of the population of measurable FIs. Finally, the most surprising result is from the Eraclea Minoa (Sicily) Upper Gypsum sample, from which all the values except one fall within the GPF. This is very impressive because the ⁸⁷Sr/⁸⁶Sr data from MSC stage 3 show the largest proportion of nonmarine waters mixing with seawater. This because ⁸⁷Sr/⁸⁶Sr is independent of salinity. Our results are in agreement with the scenario of seasonal input of marine waters into the Mediterranean not only during stage 1 but also during stage 3, as suggested by Manzi et al. (2009) and Vasiliev et al. (2017).

To summarize, our interpretation suggests that Messinian gypsum precipitated from a mostly marine water body and that post-depositional processes (crack and seal) introduced secondary low-salinity fluids that modified the parent brine composition inside the primary FIs. This appears to be the most simple and convincing explanation that fits the evolutionary model of the MSC and does not need to introduce complex additional processes, including microbially mediated reactions. To address the timing of the introduction of new fluids, we need further studies. However, our data indicate that the crystals subjected to intra-Messinian mass-wasting gravitational processes and exhumation (Vena del Gesso Basin: Roveri et al., 2006a; Sorbas and Almeria Basins: ROmodeo Salé et al., 2012; Roveri et al., 2019; Sicily: Roveri et al., 2006b, 2008; Manzi et al., 2021; Piedmont Basin: Dela Pierre et al., 2007; Cyprus: Manzi et al., 2016) show the largest FI size and yielded the lowest salinity values. Conversely, the strata still buried (Sabbioncello, Italy; Manzi et al., 2020) or exhumed during the Quaternary (Eraclea Minoa, Sicily) (Manzi et al., 2009) show the smallest pristine inclusions and the largest proportion of salinity values falling within the GPF. We cannot exclude multiple crack-and-seal phases.

Each analyzed sample shows a wide salinity range, from low to very high values. The high values fall within the expected range of the evaporation curve from seawater and have not been measured before. The large differences in salinity for FIs lying along the same growth band could only be interpreted via variable post-depositional modification of the trapped primary fluid. This study reveals that gypsum is a mineral unable to best preserve salinity information throughout geological time due to its remarkable crystallographic weakness.

The concept that Messinian gypsum precipitated from a low-salinity water body is probably the result of a strong bias in FI selection and measurement. The data obtained previously do not cover the entire range of FIs originally trapped within the crystals, but only those capable of nucleating a stable bubble after a weak stretching, which are only those having extremely low salinity. However, the low salinity was acquired via post-depositional crack-and-seal processes. The true witness of the depositional environment, i.e., normal precipitation from seawater, are the small FIs, which can be measured only after extreme stretching.

If this hypothesis is correct, then another possible implication can be inferred. The definition of the depositional environment of gypsum based on the geochemistry (δ¹⁸O and δD) of the crystallization water could also be affected. This is because the crack-and-seal mechanism may introduce new variable aliquots of water due to the syntaxial growth of gypsum during burial and exhumation.

We thank E. Salvioli Mariani (University of Parma, Italy) for technical support during microthermometric analysis on fluid inclusions, and Ahmed Sadek Mansour (Alexandria University, Egypt) who gave us the sample from Egypt. D.A. Vanko, two anonymous reviewers, and editor G.R. Dickens are acknowledged for helpful suggestions and comments on an early version of the manuscript. This work has benefited from the equipment and framework of the COMP-HUB Initiative, funded by the Departments of Excellence program of the Italian Ministry for Education, University and Research (MIUR, 2018–2022).