The Messinian Salinity Crisis as a trigger for high pore pressure development in the Western Mediterranean

Evaporites are typically described as impermeable seals that create some of the world's highest reservoir pressures beneath the salt seal. However, several laboratory studies demonstrate that evaporites can retain open pore spaces that hydraulically connect the sediments above and below them in sedimentary basins. During the Messinian Salinity Crisis (5.97–5.33 Ma), up to 2,400 m thickness of evaporites were rapidly deposited in the Western Mediterranean, which may have generated high pore fluid overpressure in the basin sediments. Here we use one‐dimensional numerical modelling to quantify the temporal evolution of overpressure at two distinct locations of the Western Mediterranean, the Liguro‐Provençal and Algero‐Balearic basins, from the Miocene to Present. We reconstruct the sedimentation history of the basin, considering disequilibrium compaction as an overpressure mechanism and constraining model parameters (such as permeability and porosity) using laboratory experiments and the literature. In the Liguro‐Provençal basin the highest overpressure of 11.2 MPa occurs within the halite during deposition of Pliocene to Quaternary sediment, while in the Algero‐Balearic basin at the base of the Emile Baudot Escarpment, the highest overpressure of 3.1 MPa also occurs within the halite but during stage 3 of the Messinian Salinity Crisis (5.55–5.33 Ma). In the Algero‐Balearic basin an overpressure of 3.1 MPa could have been sufficient to hydro fracture the sediments, which agrees with the development of fluid escape features observed on seismic reflection profiles. In general, our models with evaporite deposition rates above 20 m kyr−1 and permeabilities below 10–18 m2, suggest that high overpressure, approaching lithostatic, can be generated in salt basins.


| INTRODUCTION
The Messinian Salinity Crisis (MSC) has been described as an ecological crisis, generated by geodynamic and climate drivers (Roveri et al., 2014) including processes from plate convergence associated crustal deformation, mantle-resisted slab dragging and tearing, to isostatic responses of salt loads, possibly causing the Atlantic-Mediterranean gateway closure (Capella et al., 2020). The MSC led to the rapid sedimentation of thick layers of halite and other evaporite minerals in the Mediterranean. Evaporites are impermeable seals that create some of the world's highest reservoir pressures beneath the salt seal (Warren, 2016). However, several laboratory studies demonstrate that evaporites can have porosities of 0.5 to <10% (e.g. Casas & Lowenstein, 1989;Kröhn et al., 2015Kröhn et al., , 2017 this study) and that pore fluid flow with permeabilities from 10 -11 to < 10 -21 m 2 can occur through them by cracks and/or dilatancy of grain boundaries (e.g. Popp et al., 2001;Urai et al., 2019;Zhang et al., 2020;this study). At basin scale, this laboratory observation suggests that, despite their low permeability, evaporites are able to transmit pore fluid pressure through them. Hence, evaporite sedimentation can potentially generate overpressure within the evaporites and in the sediments below them, ultimately affecting the mechanical properties and pore fluid flow of sediments during the geological evolution of a basin.
This work contributes to the Marie Skłodowska Curie Innovative Training Network 'SALTGIANT' which aims to understand the Mediterranean Salt Giant, one of the largest salt deposits on Earth (https://www.saltg iant-etn.com/). Here we provide insights into the pore pressure evolution in the Western Mediterranean (WM) basin, where up to 1,000 m of thick well-preserved halite were deposited over a period of less than 50 kyr during the MSC (Dal Cin et al., 2016). Fluid flow and overpressure has been previously studied in WM sediments (e.g. Arab et al., 2016;Bertoni & Cartwright, 2015), although the impact of the overpressure on the hydrodynamics of the basin has primarily been addressed in Pliocene to Pleistocene sediment or areas (e.g., West Alboran Basin) where evaporites are absent (e.g. Fernandez-Ibanez & Soto, 2017;Lafuerza et al., 2009;Revil et al., 1999). Previous studies in the WM show that overpressure associated with the presence of methane gas can exist in unconsolidated shallow sediment (depths <350 m below seabed), with fluid overpressure observed to return towards hydrostatic below the overpressure zones (e.g. in ODP Site 975; Revil et al., 1999). At these shallow depths, the likely cause of overpressure is in-situ microbial degradation of organic matter that generate free gas, gas exsolution during sea-level lowering and disequilibrium compaction (Lafuerza et al., 2009;Revil et al., 1999). In contrast, studies in the Eastern Mediterranean (EM) basin have focused on fluid flow where evaporites are present, in a remnant area of the Neo-Tethyan oceanic basin that opened in the Early Mesozoic, and in an area known for being a prolific gas province (e.g. Al-Balushi et al., 2016;Bertoni & Cartwright, 2015;Eruteya et al., 2015). Focused dominantly on pipe structures in the Levant Basin, Eruteya et al., (2015) proposed their formation from (a) dissolution of Messinian evaporites (western group pipes) that predates deformation of the overburden, and (b) differential loading during late Pliocene deformation that elevated pressure within MSC evaporites (eastern group pipes). Other modelling on the petroleum system of the Levant Basin also suggest that both instantaneous drop in sea-level and evaporite loading impacted subsurface pressures (Al-Balushi et al., 2016). Quantification of overpressure from basin inception to present day and estimates of overpressure magnitude triggering fluid expulsion events during the Messinian has not been the dominant focus of previous studies in the WM. We use numerical modelling to quantify and assess the time evolution and role of pore fluid overpressure in two WM basins, the Liguro-Provençal and Algero-Balearic basins ( Figure 1). We propose that in the WM evaporite deposition during the MSC caused high overpressure that likely continues to exist within the MSC evaporites and pre-Messinian sediments, and may explain some of the fluid escape features observed on seismic data.

Mediterranean
During the Late Cretaceous, convergence of the African and Eurasian plates commenced (Olivet, 1996), with convergence of plate boundaries between the northern margin of the African plate and the Iberian Peninsula from the Late Eocene to Early Oligocene (ca. 35-30 Ma) (Jolivet et al., 2006). Following subduction of the Tethyan oceanic lithosphere and 2204 | EAGE DALE Et AL.
roll-back of the Apennines-Maghrebides subducting plate towards the north-east, south-east and south, extensional tectonics commenced in back-arc basins around the Eocene to Oligocene boundary (ca. 34 Ma; Carminati et al., 2012). The roll-back created microplate movements in the WM, developing a clockwise rotation of the Balearic promontory relative to Iberia that opened the Valencia Trough, and a counter-clockwise rotation of Corsica and Sardinia relative to Eurasia, leading to rifting of the Balearic and Ligurian extensional centres (Schettino and Turco, 2006).
The Liguro-Provencal basin comprises present day areas of the Gulf of Lion, Ligurian Sea and Mediterranean Sea between Corsica and Sardinia to the east of Menorca (Carminati et al., 2012). Continental rifting commenced during the latest Eocene to Early Oligocene with active extension in the oceanic portion of the basin continuing until the late Aquitanian to late Burdigalian (ca. 21-16 Ma) (Carminati et al., 2012). The origin and age of the Algerian basin is poorly constrained, with ages from about 25-10 Ma proposed (Carminati et al., 2012). One proposal for tectonic evolution of the area is extension terminating in the Liguro-Provençal basin and beginning in the Algerian basin during the Langhian (Mauffret et al., 2004), supported by Alger-1 well chronostratigraphy on the Algerian margin (Burollet et al., 1978). For all WM basins (excluding the Tyrrhenian basin), basin extension ended by the Late Serravallian to Tortonian (Jolivet et al., 2006).
The onset of the MSC in the Mediterranean basin at 5.97 Ma (CIESM, 2008;Roveri et al., 2014) was initiated by tectonic and glacio-eustatic processes progressively isolating the Mediterranean Sea from the world ocean. During this period, basin water volume decreased, partial desiccation occurred and evaporite minerals were precipitated | 2205 EAGE DALE Et AL. (Krijgsman et al., 1999;Lozar et al., 2018). MSC events of the Late Miocene are grouped into three stages: Stage 1 (5.97-5.6 Ma) is the first evaporitic stage; Stage 2 (5.6-5.55 Ma) includes the peak of the crisis and evaporite precipitation in deep depocentres (Roveri et al., 2014);and Stage 3 (5.55-5.33 Ma) is the final evaporitic stage and Zanclean flooding (5.33 Ma) -the end of the MSC and return to marine conditions in the Mediterranean basin (Roveri et al., 2014). Since the Pliocene, the WM basin is reconnected to the world ocean through the strait of Gibraltar (Jolivet et al., 2006). During deposition of Messinian and post-Messinian sediment in the Liguro-Provençal basin, gravity sliding, and sediment deformation occurs in the deep basin, caused by differential compaction of overburden in areas of basement structures (Maillard et al., 2003). In the Algero-Balearic basin, gravity sliding and deformation also occur in the deep basin from the late Tortonian to present, coincident with tectonic compression/uplift on the Algerian margin (Mourad et al., 2014).

| Stratigraphic framework
Depositional environments and stratigraphic lithologies of the WM have been established using borehole, outcrop and seismic facies analysis, with limited stratigraphic correlation between onshore successions and deep offshore basins (e.g. Driussi et al., 2015). The stratigraphic framework includes Oligocene to Miocene pre-Messinian successions, the three stage stratigraphic model for the Messinian, followed by Pliocene to Quaternary successions.
The Oligocene to Miocene deposits show significant facies variability from continental to brackish to marine environments (Cherchi et al., 2008). In the Gulf of Lion margin, drilled successions comprise shallow-water limestone to clastic deposits (Cherchi et al., 2008). Syn-rift Oligocene to Miocene successions from the Sardinia graben comprise similar lithologies of shallow-water limestones to clastics and hemipelagic marlstone deposits, interbedded locally with volcanic deposits, while post-rift Miocene successions comprise hemipelagic marly-silt with turbidite deposits (Cherchi et al., 2008). In the WM, Oligocene to Miocene pre-Messinian successions are characterized in general by transparent, non-reflective acoustic facies (Carminati et al., 2012).
The first stage of the MSC from Sicily for instance, is characterized by deposits of marine marlstone, alternating with diatomites and evaporites of limestone, gypsum and halite, interpreted as a deep peripheral basin (Krijgsman et al., 1999;Roveri et al., 2014), while deep basins without deep well calibration are inferred to contain deposits of organic shale and/or dolostone (Manzi et al., 2013). The second stage follows widespread desiccation of the WM basin and erosion in marginal basins, leading to deposition of primary halite, clastic deposits and resedimented evaporites in deep basins (Manzi et al., 2013). The third stage is characterized by variable evaporite deposition from primary evaporitic facies (selenite, laminar gypsum and halite cumulate) to clastic evaporitic facies (gypsrudites, gypsarenites and gypsum siltites), as well as fresh to brackish water deposits of the Lago Mare event (Krijgsman & Meijer, 2008;Roveri et al., 2014). From the Pliocene to Quaternary (P-Q), overburden sediments are dominated by deposition of marlstone with variable amounts of claystone and siltstone, intercalated locally with sandstone and volcanic deposits (Burollet et al., 1978;Hsü et al., 1978aHsü et al., , 1978bLeroux et al., 2017;Ryan et al., 1973aRyan et al., , 1973b. In the WM, the Lower Pliocene successions are characterized by semi-transparent reflections, becoming more reflective in the Upper Pliocene to Quaternary (Dal Cin et al., 2016).
In the deep basin of the WM, the MSC has also been described as a trilogy of seismic units defined as the Lower Unit (LU), Mobile Unit (MU) and Upper Unit (UU) (Roveri et al., 2014). On seismic data, the LU is characterized in general by high amplitude reflections, the MU is characterized by transparent acoustic facies of halite, while the UU is characterized by high amplitude reflectors of gypsum alternating with transparent layers of halite (Roveri et al., 2014). The stratigraphic model applied in this study integrates both these seismic stratigraphic units and the three stages of the MSC; Stage 1 corresponds to the Lower Unit (LU; 5.97-5.6 Ma), Stage 2 to the Mobile Unit (MU; 5.6-5.55 Ma), and Stage 3 to the Upper Unit (UU; 5.55-5.33 Ma; Figure 2).

| Boreholes, samples and seismic sections
To constrain the lithology of our modelled units, we reviewed seismic data and six boreholes in the WM. Lithology for the pre-Messinian succession was determined using data from boreholes Alger-1, on the Algerian shelf, and GLP-2, on the Gulf of Lion mid-slope (Burollet et al., 1978;Leroux et al., 2017). Late Miocene Messinian Upper Unit and the Pliocene to Quaternary were determined using Deep Sea Drilling Project (DSDP) Sites 122, 134, 371, 372 and Ocean Drilling Program (ODP) Site 975 (Comas et al., 1996;Hsü et al., 1978aHsü et al., , 1978bRyan et al., 1973a). Further details of primary lithologies assigned to each modelled unit are given in Section 4.3. A limited amount of geophysical log and core-based physical property data were available from offset wells GLP-2, (DSDP) Sites 371, 372 and (ODP) Site 975, with no EAGE DALE Et AL. data available in the basin-centre. Sonic log estimates of claystone porosity with depth were available from GLP-2 (Leroux et al., 2017), with density data and porosity estimates available at shallow depths for the Pliocene to Quaternary in (DSDP) Sites 371, 372 and (ODP) Site 975 (Hsü et al., 1978a(Hsü et al., , 1978b. Data for the IODP expeditions can be accessed from Expedition Science Operators at (http://web.iodp.tamu.edu/OVERV IEW/) and (https:// mlp.ldeo.colum bia.edu/logdb/).
We used evaporite core samples in our laboratory experiments to constrain some of the hydrogeological properties of evaporites in the Mediterranean and North Sea basins, prior to lithology assignment of our modelled units.  Table 1.
Seismic profiles MS-39 and E12-SF 03 ( Figure 3) were examined to ascertain thicknesses of pre-salt, Messinian and supra-salt units, and to identify locations where present-day sediment thicknesses may represent the thickest deposition in the ancient basin prior to any effects of lateral deformation, essential to conform with our 1-D vertical fluid flow modelling assumption described in Section 4 Modelling approach. The MS-39 multichannel seismic reflection data were acquired by the Italian National Institute of Oceanography and Experimental Geophysics (OGS) in 1972 as part of a regional exploration project to understand the Mediterranean tectono-stratigraphy and characterize crustal settings (Finetti & Del Ben, 2005). The SALTFLU multichannel seismic reflection data (including profile E12-SF 03) were acquired in 2012 to provide detailed preand post-stack time migration data and RMS velocity data over the continental slope, particularly the Emile Baudot escarpment, and deformed sequences of the Algero-Balearic abyssal plain (Wardell et al., 2014). See Figures S1 and S2 for regionally extensive and uninterpreted seismic profiles MS-39 and E12-SF 03.
We used well and seismic velocity data from the literature and project company GALSI S.p.A to convert seismic time thickness units to depth. 2D ultra-high resolution multichannel seismic velocity data were available from GALSI S.p.A providing high quality velocities for Pliocene to Quaternary (P-Q) and MSC Stage 3 (Upper Unit) seismic units along the GALSI (Gas Pipeline -Algeria via Sardinia to Italy) route acquired from 2007 to 2008. Below the Upper Unit (UU), seismic velocities are considered to be of poor quality. Refer to Section 4.3 for information on velocities data used in this study.  Garrison et al. (1978); (2) Lugli et al. (1999). Note that connected porosity estimates were determined using mercury injection porosimetry. P c , confining pressure.

| Laboratory experiments
Seventeen sediment cores were used, and from each of these, a total of twenty smaller core plugs were cored, and cut and their ends ground flat. This resulted in nineteen discrete samples for porosity and permeability determinations, and one for X-ray CT-scan image analysis. Cores were selected to represent a range of evaporite lithologies and avoid impurities (e.g. claystone), while core plugs were selected on texture variation and to avoid fractures where possible.
3.2.1 | Porosity and permeability determination A set of 1.1-2 cm height and 2.5 cm diameter cores samples were used for high pressure (17.2 and 1.38-4.83 MPa of confining and pore pressure) permeability to helium and Hginjection porosimetry (Micromeritics AutoporeTM IV 9520 system) determinations at the University of Leeds. A second set of ca. 2 cm height, 5 cm diameter samples were used for absolute porosity estimates with helium pycnometer and absolute permeability with N 2 at the National Oceanography Centre (NOC) in Southampton. Further details of primary lithologies assigned to each modelled unit are given in Section 4.3. Porosity and permeability were measured at room temperature (ca. 20°C), at atmospheric pressure conditions for porosity and under a minimum hydrostatic confining pressure of 1.5 MPa for permeability to ensure rig sealing during gas flow-through. The gas permeability to N 2 was estimated using steady state flow (SSF) and pore pressure transmission (PPT) methods (e.g. Falcon-Suarez et al., 2017), depending on the sample permeability. For all the samples, we first attempted to measure permeability using the SSF, based on Darcy's law, the most widespread method for high to moderate permeability media (above 10 -16 m 2 ). For those samples with permeability below 10 -17 m 2 , we used the PPT method instead, an alternative based on transient states of the pore pressure. The PPT method was proposed by Metwally and Sondergeld (2011) based on the pulse decay method introduced by Brace et al. (1968), which consists of inducing pore pressure disequilibria in the rock and determining the permeability through the evolution of pore pressure-time decay curves towards the original steady state. For further details on the SSF and PPT methods refer to, for example, Metwally andSondergeld (2011) andFernandez-Ibanez andSoto (2017).
In all the permeability determinations, Klinkemberg's correction was applied to correct the deviation resulting from slippage effect of the gas (Klinkenberg, 1941), and transform the apparent permeability into absolute permeability.

| Mercury injection porosimetry
The porosity of a selection of the core plug samples was also analysed using mercury injection porosimetry with the Micromeritics AutoporeTM IV 9520 system. This model has four low pressure ports and two high pressure chambers. As mercury is a non-wetting fluid, the pressure must build up before mercury intrudes into a certain pore size and the interface crosses the throat between pore bodies. The balance between internal and external forces or pressures acting on an interface can be described by the Young-Laplace equation. The samples are cut into suitable size depending on their porosity and the penetrometer to be used. Clean and dry samples, of known weight, are then loaded into a penetrometer and evacuated. The penetrometer is automatically backfilled with mercury. The pressure is then increased to 25 psi (0.17 MPa) in the low pressure port, and up to 60,000 psi (413 MPa) in the high pressure chamber following pre-selected pressures. After reaching each pressure increment the volume of mercury intruded is recorded. Each penetrometer has been individually calibrated, therefore, the volume of mercury needed to fill the penetrometer at ambient conditions is used to calculate the bulk volume of the sample. The total volume of mercury injected is recorded assuming that at 60,000 psi all the pore volume has been filled. The grain volume is the difference between the sample bulk volume and mercury injected volume. Then the grain density can also be obtained. The pore throat size distribution and other properties can be calculated from this information (ASTM D4404-84, 2004). If necessary, a manual volume conformance and other corrections can be applied during data interpretation. All results presented here have been conformance corrected.

| X-ray computed tomography
X-ray micro-CT (XCT) imaging was carried out on one halite core plug to fully understand anomalous results of permeability. To improve the signal to noise ratio, the core plug was cored and cut down to a core size of 14 mm diameter with 20 mm height. A scan image to 10.1 μm voxel resolution was achieved. The scan was conducted using a micro-focus Custom Nikon HMX ST Scanner at the University of Southampton (e.g. Callow et al., 2018). The settings used on the HMX are as follows: a source to object of 40.4 mm, source to detector of 797.9 mm, 200 kVp peak voltage, no pre-filtration of the beam, 134 ms exposure time, 3,142 projections (2 frames per projection) and voxel size of 0.01 mm.

Mediterranean and North Sea basins
Tables 1 and S1 and Figure 4 show the results of permeability and porosity measurements on 19 evaporite samples. Although, the measured permeability range is wide (10 -13 to 10 -20 m 2 ), most of the anhydrite sample permeability are between 10 -17 and 10 -18 m 2 and have similar absolute porosity of 2.4%-2.93% at hydrostatic confining pressure (P c ) of 1.5 MPa. The permeability of the anhydrite is stress dependant, F I G U R E 4 Physical property data compilation for evaporites. (a) Global permeability ranges of evaporites including this study's laboratory results of permeability obtained for Permian and Miocene anhydrite, Miocene gypsum and fractured Miocene halite. Boundary of undisturbed/ undamaged subsurface halite <10 -21 m 2 (Stormont, 1997;Warren, 2016) with disturbed halite permeability taken from Stormont, 1997. (b) Permeability and connected porosity measurements for Miocene evaporites from this study. (c) Grain density and connected porosity measurements for Miocene evaporites from this study. decreasing to ca. 10 -21 m 2 when P c increases above 6.0 MPa. Two anhydrite samples show anomalously high permeability (about 10 -13 m 2 ), and absolute porosity (6.1% and 13.3%), likely caused by their irregular speckled outer surface texture having prevented adequate rig sealing. The five gypsum samples show similar permeabilities in the range of 10 -17 to 10 -18 m 2 and connected porosity from 1.5% to 3.1%, within the P c range 1.5-17.2 MPa, indicating low stress dependence for both properties. In contrast, the three halite samples of similar origin showed anomalously high permeability of up to 10 -13 m 2 at P c of 1.5 MPa, and high stress dependence, as this value decreased to 10 -16 m 2 at 17.2 MPa. X-ray computed tomography on the halite shows the presence of fractures and isolated pore spaces ( Figure 6). All halite samples show low connected porosity of 1.0%-2.0%. Testing on the samples of kainite and polyhalite showed low connected porosity of 0.5% and 3.6% respectively. We used the results of dry density and porosity for gypsum and halite, and permeability of gypsum as input parameters in our modelling. The measured halite permeabilities were disregarded from modelling as they are significantly lower than most values reported in the literature (Figure 4), likely because of pre-existing micro fracturing in the samples. Hence, undisturbed halite permeabilities from literature were used (e.g. Beauheim et al. 1991;Brodsky, 1994). The measured anhydrite permeability and porosity results were also disregarded from modelling, as anhydrite was considered unlikely across our model areas. Further details on gypsum dehydration to anhydrite results are provided in Section 5.

| 1-D disequilibrium compaction model
Pore fluid overpressure generation in our 1-D models considers brine as the only pore fluid and overpressure due to the disequilibrium compaction mechanism. Sea level changes do not affect overpressure in sediments saturated with a near incompressible fluid such as water (e.g. Liu & Flemings, 2009). Hence, the >1,500 m sea level fall in the WM (e.g. Hsü et al., 1977) is not considered here. Overpressure due to disequilibrium compaction occurs by an imbalance between increasing compressive stress and ability of the sediment to expel fluid (Tingay et al., 2009). During slow burial (normal compaction) as vertical load increases, pore-volume decreases and pore fluid is expelled from the sediment allowing an equilibrium between overburden and reducing pore-volume to be maintained (Osborne & Swarbrick, 1997). However, during rapid burial as vertical load increases, if accompanied by fluid that cannot be expelled rapidly, part of the load will be supported by the pore fluid resulting in pore fluid pressure increasing above hydrostatic (a process called disequilibrium compaction; for example, Osborne & Swarbrick, 1997). The magnitude and time evolution of overpressure depends on the balance between sediment loading and compressibility, pore fluid dissipation controlled by permeability, and drainage (dissipation) distance. Note that in this work the porosity includes any type of connected void such as intergranular pores or micro fractures along grain boundaries. Based on seismic data interpretation, we apply our 1-D models in areas with sufficient laterally extensive horizontal layers and limited tectonic compression, so horizontal fluid flow is assumed to be negligible. We account for water viscosity and density changes with variations in temperature, pore pressure and salinity.
The detailed description of the mathematical and numerical models are given in Marín-Moreno et al. (2013a, 2013b, and here we only provide the main equations (Table 2). To describe the mechanical compaction of sediments we consider that the change in porosity is a function of vertical effective stress (Equation 1), where depth change is controlled by Equation (2). The change in lithostatic pressure with time is expressed in terms of sediment thickness h (Equation 3). The stress compaction factor β in Equation (1) can be related to the depth compaction factor (Sclater & Christie, 1980) using Equation (4). Here we assume the empirical compaction factor β is equivalent to the bulk compressibility of the saturated sediment, as described by Hart et al. (1995). To describe fluid flow, we use Darcy's relationship (Equations 5 and 6) and assume that changes in permeability depend on changes in porosity caused by changes in effective stress (Equation 7). Finally, combining the above equations the disequilibrium compaction model is given in Equation 8.

| Modelling strategy and scenarios
Our modelling strategy ( Figure 5) commences with a set of rock hydrogeological properties for each unit from laboratory experiments performed in this study and the literature (Table 1 and Figure 4). We then run our 1-D disequilibrium compaction model using these rock properties and estimates of sedimentation rate from pre-compacted thickness and sedimentation time for each unit. Pre-compacted thicknesses are determined applying a percentage increase above present-day thicknesses estimated from seismic data. If the present-day modelled compacted thicknesses and present-day seismicderived thicknesses are similar, within a 5% tolerance, we assume the calculated present-day pore pressure, bulk density, porosity, compressibility and permeability depth profiles represent those in-situ. Otherwise, we re-evaluate input parameters, considering their inherent uncertainties, and rerun the model until the observed and calculated thicknesses are within tolerance. A corollary of this approach for model validation is the assumption that the hydrodynamic and (1) β = Stress compaction factor γ = Parameter controlling the change in permeability with porosity µ = Viscosity ρ = Sediment density ρ f = Fluid density σ' zz = Vertical effective stress ϕ = Porosity ϕ 0 = Initial porosity at seabed conditions K = Permeability tensor K zz = Vertical permeability k i = Permeability k i0 = Initial permeability at seabed conditions c = Depth compaction factor g = Gravitational acceleration h = Sediment thickness at seabed conditions P = Total pore pressure P* = Overpressure Time evolution in lithostatic pressure with sediment thickness Depth compaction factor to stress compaction factor conversion (4) Permeability versus porosity relationship compaction history generated by the model represent those of our study area. For both the Liguro-Provençal basin and the Algero-Balearic basin at the base of the Emile Baudot Escarpment, our approach has been to model three scenarios incorporating low, most likely and high estimates of overpressure, which cover the full range of possible variations in fracture limit and permeability. Here we define the fracture limit as the ratio of overpressure to vertical effective stress under hydrostatic conditions, sometimes defined in the geoscience literature as λ*, above which vertical fractures can occur. To represent changes in fluid flow due to the generation and propagation of vertical fractures, on a numerical cell-by-cell basis, we assume that once overpressure exceeds the fracture limit the permeability increases by two orders of magnitude. This increase in permeability is related to a threshold value above which permeability does not influence our results (see Section 5 for further discussion). In areas without significant tectonic compression and with sufficiently extensive horizontal strata, the major principal stress is vertical and the intermediate and minor stresses are in the horizontal plane. Hence, here we assume that fractures propagate vertically and open horizontally and with fracture limits of 0.8 ± 0.1 (e.g. Luo et al., 2017;Nikolinakou et al., 2014). In-situ fracture pressure measurements (e.g. traditional and extended leak-off test data) were only available from wells on the basin margin, limiting our ability to constrain fracture limits in the deeper basin-centre. As no data exists for the basin-centre halite, the simplest model we could assume was a horizontal to vertical effective stress ratio of 0.8, taken from initial stresses applied in other modelling projects near a salt diapir under hydrostatic conditions (Nikolinakou et al., 2014). The relatively high λ* implicitly accounts for the additional overpressure required to also overcome the tensile strength of the material. As we apply 1D modelling, the magnitude of the two horizontal stress components do not influence our results.
We positioned our Liguro-Provençal model in the south of the basin between the North Balearic and Catalan transverse fracture zones, an area also referred to in the literature as the North Balearic Basin (Figure 1; black box) where undeformed to mildly deformed sediment on seismic data (e.g. seismic profile SPBal-15 & SPBal-27) progressively deepens from the North Balearic fracture zone towards the Gulf of Lion (Maillard et al., 2003(Maillard et al., , 2020. In comparison, the Gulf of Lion slope between the Catalan and Arlesian transverse fracture zones (Figure 1) tilts towards the southeast with listric faults in the upslope, salt anticlines and translation in the mid-slope and contraction and diapirs in the downslope area (Maillard et al., 2003). An extensive region with relatively undeformed sediment also exists on seismic reflection data across the basin plain in the Gulf of Lion (Mianaekere & Adam, 2020). To the east of our Liguro-Provençal model, salt diapirs exist restricted to a northeast to southwest boundary of the deep basin where steps in top basement reside (Figure 1; Maillard et al., 2003). Two limitations of our method for the WM are that it is only applicable for sediments that are relatively undeformed laterally and that we do not account for the overprint in overpressure generated during formation of a diapir. Our Liguro-Provencal model is located ca. 6 km from the synclinal axis of a salt diapir. If we had considered repositioning the model ca. 8 km further west of its current location (total of 14 km from the synclinal axis of the diapir), unit thicknesses (pre-kinematic) would be similar and so overpressure estimates would remain similar to our 1-D model results.
The most likely scenario uses a fracture limit of 0.8 and permeability at seabed for gypsum and halite of 10 -18 and 10 -20 m 2 respectively. The low scenario uses a fracture limit of 0.7 and permeability of gypsum and halite of 10 -17 and 10 -19 m 2 respectively. The high scenario uses a fracture limit of 0.9 and permeability of gypsum and halite of 10 -19 and 10 -21 m 2 respectively. See Table 3 for a summary of the non-default parameters used in the low and high scenarios. Four additional scenarios are presented that evaluate the sensitivity of overpressure to common halite properties (porosity and permeability), to understand the impact of downward fluid migration on our models, T A B L E 3 Non-default physical property parameters used in uncertainty analysis of each model or sensitivity analysis of additional scenarios to ground truth whether mineral dehydration is plausible at our model locations in the Algero-Balearic basin, and to determine timing of fluid expulsion in the WM. Sensitivity of overpressure to halite properties was modelled using halite permeabilities of 10 -16 -10 -22 m 2 and initial seabed porosity of 0.1%-4.0%. For timing of fluid expulsion events, conservative values for fracture limit of 0.7, permeability for gypsum and halite of 10 -19 and 10 -21 m 2 respectively and a Lower Unit package thickness of gypsum of 1,405 m are used. The sensitivity of the model overpressure to downward fluid flow into basement rock, was modelled using ranges in pre-Messinian claystone permeabilities of 10 -17 -10 -22 m 2 at porosities of 2%-14% for a 4,000 m subsurface depth. To ground truth mineral dehydration as an overpressure mechanism in the region, heat flow of 80-120 m W/ m 2 (Carballo et al., 2014), thermal conductivities for marlstone of 1.5-3.0 W/m K/m (Erickson & Von Herzen, 1978), thermal conductivities for gypsum of 1.0-1.3 W/m K/m (Balkan et al., 2017; Robertson & Geological Survey (U.S.), 1988), thick basin-ward units and seabed temperature are used to estimate the thermal structure of the 1D sediment column. Combining these temperature data with pressure data we estimate the P-T conditions of the Algero-Balearic Upper Unit Gypsum relative to the boundaries of the dehydration reaction. Here we assume that present-day thicknesses adequately represent maximum burial depth of the sediment. Location of the basin-ward thicknesses on the lowermost slope of the continental rise is given in Section 6. See Table 3 for a summary of the non-default parameters used in each of the additional scenarios.

| Modelling parameters and boundary conditions
Our seismic stratigraphic model extends to 4 km below the seabed and comprises five seismic stratigraphic units, three of which represent the Messinian Salinity Crisis. Seismic units were interpreted on PSTM data with thicknesses in time for each unit converted to depth using a range of velocities from well and seismic velocity analysis over the Mediterranean (Table 4). Present-day thicknesses were selected from the

Unit Velocity range (m/s) Remarks
Pliocene to quaternary (P-Q) 2,000-3,150 • 2,000 m/s is reported from the Q5 to seafloor surface (top P-Q layer) around the Gulf of Lion using an average of velocities derived from various datasets (1). • In the deep basin of the Western Mediterranean, we expect an average velocity of 2,930 m/s using velocities from 2D Ultra-high resolution Multichannel Seismic data over the Algeria to Sardinia basin centre (2). • 3,150 m/s is reported from the P11 to PXX surface (base P-Q layer) around the Gulf of Lion mid-slope GLP-2 well (1).
Messinian upper unit (UU) 3,100-3,500 • 3,100 m/s is reported from seismic velocity analysis of profile MS-39 in the Western Mediterranean (3). • In the deep basin of the Western Mediterranean, we expect an average velocity of 3,300 m/s using velocities from 2D Ultra-high resolution Multichannel Seismic data over the Algeria to Sardinia basin centre (2).

|
mid-range of thicknesses calculated for each unit, except where high quality seismic velocity data existed from the GALSI project. We use the GALSI data to calculate the thicknesses of Pliocene to Quaternary (P-Q) and MSC Stage 3 (Upper Unit) of the Liguro-Provençal deep basin model. A single representative lithology per unit is selected, using seismic stratigraphy and literature sources (Table 5). For the five units, the primary lithologies were marlstone, claystone, halite and gypsum. Marlstone was used for the pre-Messinian and Pliocene to Quaternary (P-Q) units, claystone for the MSC Stage 1 (Lower Unit), halite for the MSC Stage 2 (Mobile Unit), and gypsum for the MSC Stage 3 (Upper Unit).
Duration of deposition of the modelled units from Miocene to Present is provided in Table 5. Duration of deposition in our models for the pre-Messinian unit ranges from 14 Myr along the basin edge to 10 Myr in the deep central oceanic location (Carminati et al., 2012). Average fluid and solid properties and other modelling constants are provided in Table 6. We assume fully watersaturated sediment for all scenarios with initial seabed density and viscosity for water of 1,028 kg m 3 and 0.0012 Pa s respectively. For marlstone units we use an initial seabed porosity of 30% from claystone porosity trends (Magara, 1980), an initial compaction factor of 0.4 km −1 reported in Marín-Moreno et al. (2013a, 2013b for similar sediments, an irreducible porosity of 10%, and a permeability at the seabed of 10 -17 m 2 from porosity and permeability trends for argillaceous material (Neuzil, 1994). For evaporite units, we use an initial seabed porosity of 2.0%-3.0%, an initial compaction factor of 0.1-0.2 km −1 , an irreducible porosity of 1.0%, and a permeability at seabed for the most likely scenario of 10 -20 m 2 for halite and 10 -18 m 2 for gypsum estimated from laboratory tests as part of this study.
We assume a seabed (top boundary) temperature of 13°C, which corresponds to the estimated temperature at water depths of 2,585-2,638 m (Manca et al., 2004) and an average geothermal gradient of 36°C km −1 (Erickson & Von Herzen, 1978). The temperature is only used to calculate changes in pore fluid density and viscosity with depth. We impose boundary conditions of zero overpressure at the top of the models representing the seabed, and zero flow at the base of the models.
The mathematical model (Equation 8) is solved in Matlab (R2017) using an implicit finite difference scheme with backward differences to approximate the time derivative and second-order centred differences in space, an harmonic average to estimate the permeability in the interface between cells, and a fully compacted coordinate system for the depth axis . The numerical model uses 400-800 cells in the z-direction and 800-1,600 time steps per unit. We run the default model with different mesh sizes to assess their influence on our results and select the mesh size from which further refinement produced negligible changes.

| Numerical model results
To
The effect of permeability on our scenarios show that the most likely and high scenarios give similar results in terms of pore pressure within the halite using permeabilities of 10 -20 and 10 -21 m 2 respectively (Figure 7). This means that for permeabilities lower than 10 -20 m 2 , differences in pore pressure are small.

Deposition from 16 to 5.3 Ma
After the peak of the MSC, MSC Stage 3 is characterized by deposition of 646 m of Upper Unit gypsum (5.55-5.3 Ma; sedimentation rate 2,936 m Myr −1 ) contributing to further loading of underlying evaporitic units with overpressure increasing to 11.2 and 12.4 MPa within the base of the MSC Stage 2 halite and pre-Messinian units respectively (red lines).

Deposition from 16 Ma to present day
Following the MSC, deposition of 1,480 m of marlstone during the Pliocene-Quaternary (5.3 Ma to present day; sedimentation rate 278 m Myr −1 ) allowed pore fluid to dissipate to near hydrostatic pressure within the MSC Stage 3 gypsum and the Pliocene to Quaternary (black lines). However, below the MSC Stage 2 halite seal with present-day λ* near hydro fracture conditions, overpressure up to 21.6 MPa is retained within the pre-Messinian sediment. The present-day overpressure that remains is located below 1,409 m where porosity deviates from the clay normal compaction trend (Figure 8) as is expected during disequilibrium compaction. In Alger-1, (DSDP) Sites 134, 371 and 372 and (ODP) Site 975, the Pliocene to Quaternary (P-Q) overburden sediments are dominated by deposition of marlstone with various mixtures of sand, silt and claystone, while GLP-2 is dominated solely by carbonated claystone. When comparison is made to the sedimentation rate versus fluid retention depth relationship for silt, silty claystone and claystone from global data (Swarbrick, 2012), assuming a sedimentation rate of 278 m Myr −1 , and 'silty' lithology, we would expect top of overpressure to begin near the base of our P-Q unit. This is consistent with our estimates and hydrostatic pressures maintained to 2,000 m depth below the seabed in wells like Andalucia-G1 (Fernandez-Ibanez & Soto, 2017). Applying the same sedimentation rate and an alternative 'silty shale' lithology, we would expect top of overpressure to occur at depths anywhere from ca. 900 m below seabed to near base of our P-Q unit based on global equivalent examples of sedimentation rate.

Deposition from 20 to 5.97 Ma
Commencing from the Early Miocene (Pre-Messinian; yellow lines), deposition of 579 m of marlstone (20-5.97 Ma; sedimentation rate of 41 m Myr −1 ) allowed pore fluid dissipation to hydrostatic pressure.

Deposition from 20 to 5.6 Ma
Assuming the Lower Unit of MSC Stage 1 is absent along the edge of the basin (1 m inferred, 5.97-5.6 Ma; sedimentation rate of 2.7 m Myr −1 ), hydrostatic pressure conditions persist to 5.6 Ma (green lines).

1-D models Reference A-B L-P
Acceleration of gravity m/s 2 9.80 9.80 Robinson et al. (1995)

Deposition from 20 Ma to present day
Following the Messinian Salinity Crisis, deposition of 286 m of marlstone during the Pliocene-Quaternary (5.3 Ma to present day; sedimentation rate of 54 m Myr −1 ) allowed pore fluid to dissipate by up to 2.6 MPa within the MSC Stage 2 halite and pre-Messinian sediment (black lines).

| Sensitivity of the model to common halite properties
We evaluated the impact of uncertainty in initial seabed porosity, permeability and sedimentation rate on overpressure development during halite deposition (Figure 12), as this is the primary unit contributing to the major increase in λ* (Figure 13). We considered halite thicknesses of 200-1,000 m, the latter based on thickness estimates We used halite permeabilities of 10 -16 -10 -22 m 2 based on global literature ranges, derived from a combination of laboratory tests, modelled and inferred values, and F I G U R E 7 (a) Present-day pressure and (b) overpressure from seabed estimated for the Liguro-Provençal model. Red lines are uncertainty ranges. Results were calculated applying variation in fracture limit from 0.7 to 0.9 and permeability of evaporites from 10 -17 to 10 -21 m 2 . The most likely scenario (red dotted line) uses a fracture limit of 0.8 and permeability of gypsum and halite of 10 -18 and 10 -20 m 2 respectively. The low value scenario (red dashed line) uses a fracture limit of 0.7 and permeability of gypsum and halite of 10 -17 and 10 -19 m 2 respectively. The high value scenario (red solid line) uses a fracture limit of 0.9 and permeability of gypsum and halite of 10 -19 and 10 -21 m 2 respectively. The column on the right side shows the five units modelled [Colour figure can be viewed at wileyonlinelibrary.com] F I G U R E 8 Present-day variations of density, porosity compressibility and permeability with depth for the Liguro-Provençal model. (a-d) Results were calculated applying variation in fracture limit from 0.7 to 0.9 and permeability of evaporites from 10 -17 to 10 -21 m 2 . The most likely scenario (red dotted line) uses a fracture limit of 0.8 and permeability of gypsum and halite of 10 -18 and 10 -20 m 2 respectively. The low value scenario (red dashed line) uses a fracture limit of 0.7 and permeability of gypsum and halite of 10 -17 and 10 -19 m 2 respectively. The high value scenario (red solid line) uses a fracture limit of 0.9 and permeability of gypsum and halite of 10 -19 and 10 -21 m 2 respectively. The column on the right side shows the five units modelled. our experimental data. Halite permeabilities above 10 -17 m 2 generate hydrostatic pressures. Hence, for the low overpressure scenario with a fracture limit of 0.7 and the highest halite permeability of 10 -19 m 2 , if the overpressure exceeds the fracture limit, our assumed permeability increase of two orders of magnitude results in a halite permeability of 10 -17 m 2 , which is the threshold above which permeability does not influence our results. When the permeability drops below a threshold of about 10 -18 m 2 , halite with thickness greater than 600 m develops overpressure above 1 MPa. In contrast, halite with thickness of 200 m requires permeability below 10 -20 m 2 , to generate and maintain the same overpressure magnitude. This two orders of magnitude difference in threshold permeability is related to the ability of permeability to dissipate overpressure for a given length scale and time scale.
In our 1-D models, for the same time scale, the thinner the layer the shorter the distance the fluid needs to travel to dissipate overpressure, and so the lower the permeability needed to generate and maintain the same amount of overpressure. For halite thicknesses of 600-1,000 m, a permeability of 10 -20 m 2 develops overpressure within the range 3.9-7.1 MPa. Below 10 -21 m 2 , that being the permeability of pristine, undamaged halite, overpressure for a 1,000 m halite remains high at 7.7-8.5 MPa. When comparison is made for permeability ranges of 10 -20 -10 -22 m 2 , minor variation in overpressure, up to 1.3 MPa, is obtained.
Shallow halite layers such as that of Quaternary halite in the Saline Valley, CA display low porosities (<10% at 10 m below ground level) and tightly cemented layers below a depth of 45 m (Casas & Lowenstein, 1989). In our study, lower connected porosities of 1.0%-2.7% were obtained from laboratory testing of shallow Messinian halite collected in Sicily. Integrating literature sources and our laboratory measurements of porosity, we tested the impact of uncertainty in initial halite seabed porosity of 0.1%-4.0% on overpressure. For an initial seabed porosity of 1.0%, a significant increase in overpressure up to 6.5 MPa is obtained for a 1,000 m thick halite. For initial seabed porosities above 1.0%, overpressure plateaus with only minor increase in overpressure by about 0.9 MPa.
There is considerable uncertainty concerning the stratigraphic framework for the MSC evaporites, as well as their absolute chronostratigraphy in the deep basins owing to limited well control, lack of chronostratigraphic constraint, studies not structured within a regional context and scientific debate on the origin of the evaporites (Hardie & Lowenstein, 2004;Krijgsman et al. 1999).
Accounting for the uncertainty in stratigraphic models for the MSC, we test the impact of sedimentation rate on overpressure, using halite thicknesses of 200-1,000 m and total duration of Messinian halite deposition of 50-90 kyr. For duration of deposition of 50-90 kyr and halite thickness of 1,000 m (sedimentation rates 11-20 m kyr −1 ), a minor difference in overpressure, up to 0.6 MPa, is obtained. Halite with a lower thickness of 200 m and the same duration of deposition (sedimentation rates 2-4 m kyr −1 ) show even lower magnitude difference in overpressure, of 0.25 MPa.

| Sensitivity of the model to downward fluid migration
We evaluated the impact of downward fluid migration from pre-Messinian sediment into basement rock and the effect of permeability variation of pre-Messinian sediment on this type of migration. To do this, we assumed a highly fractured basement rock by imposing a zero overpressure bottom boundary condition. Although the nature of the basement in the Mediterranean is variable, we expect there to be oceanic crustal igneous rock in the Liguro-Provençal basin where our model is located (Figure 1; Sàbat et al., 2018). If a boundary condition of zero overpressure is imposed at the base of the model representing full dissipation through the basement, the ability to retain overpressure within pre-Messinian units depends largely on its permeability which is poorly constrained. We tested pre-Messinian permeabilities of 10 -17 -10 -22 m 2 (Neuzil, 1994) reasonable values at porosities of 2%-14% based on claystone compaction trends at a depth of ca. 4,000 m (Allen & Allen, 2013). For downward flow and permeability of 10 -17 m 2 , present-day overpressure is near hydrostatic pressures within the pre-Messinian sediment. For permeability lower than 10 -19 m 2 overpressure develops, which increases mid-unit up to 32.2 MPa (λ* of 0.75) for a permeability of 10 -22 m 2 . Below this at the boundary between pre-Messinian and basement, a regression in overpressure to hydrostatic conditions is modelled (see supporting material Figure S3). F I G U R E 1 0 (a) Present-day pressure and (b) overpressure from seabed estimated for the Algero-Balearic model. Red lines are uncertainty ranges. Results were calculated applying variation in fracture limit from 0.7 to 0.9 and permeability of evaporites from 10 -17 to 10 -21 m 2 . The most likely scenario (red dotted line) uses a fracture limit of 0.8 and permeability of gypsum and halite of 10 -18 and 10 -20 m 2 respectively. The low value scenario (red dashed line) uses a fracture limit of 0.7 and permeability of gypsum and halite of 10 -17 and 10 -19 m 2 respectively. The high value scenario (red solid line) uses a fracture limit of 0.9 and permeability of gypsum and halite of 10 -19 and 10 -21 m 2 respectively. The column on the right side shows the four units modelled [Colour figure can be viewed at wileyonlinelibrary.com] | 2221 EAGE DALE Et AL.

| Gypsum dehydration to anhydrite
In the Algero-Balearic basin, polygonal faults have been interpreted within the upper evaporites of the MSC and lowermost Pliocene sequences, suggesting presence of past fluid expulsion and migration events (Bertoni & Cartwright, 2015). Lofi (2018) interpret the polygonal faulting to be generated by overpressure induced by fluid from the gypsum to anhydrite dehydration process. Anhydrite has been cored before from the Upper Unit (UU) of (DSDP) Site 371 in the Algero-Balearic Basin, however the well was drilled in a zone where numerous shallow magnetic anomalies and a thin veneer of evaporites are present (Figure 1; Hsü et al., 1978aHsü et al., , 1978b. To understand if gypsum dehydration occurs in the location of the SALTFLU seismic data, we first evaluate pressure and temperature conditions of the Upper Unit gypsum at the base of the Emile Baudot Escarpment (Figure 1; Figure 13a Model A location) relative to the boundaries of the dehydration reaction (Figure 13b). Using the parameters described above (Section 4), we show that the Upper Unit of gypsum at the base of the Emile Baudot Escarpment is unlikely to reach the pressure and temperature conditions required for gypsum-anhydrite transformation. Considering thicker basinward units on the lowermost slope of the continental rise (Figure 1 green star on strike direction and in close proximity with Figure 13a Model B location), fluid release from evaporitic dehydration is possible if heat flow exceeds 105 mW m −2 in combination with low thermal conductivities for marlstone and gypsum of 1.5 and 1.0 W/m K/m respectively. However, these low modelled thermal conductivities are inconsistent with higher values obtained during (DSDP) Leg 42A (Erickson & Von Herzen, 1978). Alternatively, our disequilibrium compaction models suggest that sediment loading over the 5.55-5.33 Ma period can cause sufficient overpressure to hydro fracture the underlying MSC Stage 2 (Mobile Unit) halite (Figure 13c), allowing fluid to migrate into the Upper Unit of gypsum and leading to development of a polygonal fault system.

DISCUSSION
Our sensitivity analysis of evaporite petrophysical properties show that permeability is the dominant parameter controlling the generation of pore fluid overpressure. However, a broad range of permeability values are reported in literature ( Figure 4) which reduces the predictive ability of overpressure from numerical models, as illustrated in our modelling of the MSC Stage 2 (Mobile Unit) halite. Laboratory measurements of permeability from high quality, undamaged evaporites from borehole cores is then essential to produce reliable predictions. When permeability measurements of representative evaporites under undisturbed conditions are not available, low (most likely) and high permeability overpressure modelling scenarios and the threshold above which permeability does not influence overpressure results should be provided (Figure 12). F I G U R E 1 1 Present-day variations of density, porosity compressibility and permeability with depth for the Algero-Balearic model. (a-d) Results were calculated applying variation in fracture limit from 0.7 to 0.9 and permeability of evaporites from 10 -17 to 10 -21 m 2 . The most likely scenario (red dotted line) uses a fracture limit of 0.8 and permeability of gypsum and halite of 10 -18 and 10 -20 m 2 respectively. The low value scenario (red dashed line) uses a fracture limit of 0.7 and permeability of gypsum and halite of 10 -17 and 10 -19 m 2 respectively. The high value scenario (red solid line) uses a fracture limit of 0.9 and permeability of gypsum and halite of 10 -19 and 10 -21 m 2 respectively. The column on the right side shows the four units modelled [Colour figure can be viewed at wileyonlinelibrary.com] 2222 | EAGE DALE Et AL.
Overpressure build-up up to hydro fracturing during the MSC has likely caused fluid expulsion events in the WM basin. Fluid expulsion related features are evident on seismic data with examples of mud volcanoes, pipes and polygonal faulting in sediments of the Central and Western Mediterranean (Bertoni & Cartwright, 2015). Using seismicbased evidence across the entire Mediterranean, a conceptual framework for timing of fluid expulsion during the MSC indicates three possible fluid flow stages, the first commencing in the early stage of the MSC before about 5.6 Ma, the second during deposition of MSC Stage 2 basin centre evaporites from 5.6 to 5.53 Ma, and the third during deposition of MSC Stage 3 (Upper Unit) evaporites from 5.53 to 5.33 Ma (Bertoni & Cartwright, 2015). To evaluate the role that evaporite deposition played on these three stages, we model overpressure applying a scenario of conservative values for fracture limit of 0.7, permeability for gypsum and halite of 10 -19 and 10 -21 m 2 respectively and a high Stage 1 (Lower Unit) thickness of 1,405 m, as observed on seismic data in the Gulf of Lion deep basin (Leroux et al., 2017) with an alternative low permeability evaporite scenario of gypsum. Our models show that sediment loading by this thickness of LU gypsum does not cause overpressure to increase above hydro fracturing in this first stage, from 5.97 to 5.6 Ma. From our modelling in the WM, we identify two possible timings of fluid expulsion events relative to the MSC, the first by sediment loading of Stage 2 (Mobile Unit) halite from ca. 5.58-5.55 Ma and the second by sediment loading of Stage 3 (Upper Unit) evaporites from ca. 5.55-5.33 Ma causing overpressure of the underlying MSC Stage 2 halite to increase above hydro fracturing. We therefore show that in the WM the fluid expulsion events triggered by an increase in overpressure above hydro fracturing likely started during and after deposition of Stage 2 halite (Mobile Unit). The former event timing during Stage 2 differs slightly from seismic observations of basincentre pockmarks in the Eastern Mediterranean, described to have formed in the early stages of the MSC related to sea-level drop (Bertoni & Cartwright, 2015), while the latter event timing is consistent with brecciated limestone in Central Mediterranean outcrop and seismic observations of polygonal faulting in the Western Mediterranean, formed at the late stage of the MSC up to the early Pliocene (Bertoni & Cartwright, 2015;Iadanza et al., 2013).
The distribution and thickness of the Stage 1 Lower Unit in the Algero-Balearic basin is not entirely known, due to complex structures of salt deformation, erosion, and seismic imaging effects (Dal Cin et al., 2016). The Lower Unit appears absent on seismic line E12-SF03 in the Algero-Balearic basin at the base of the Emile Baudot Escarpment. To understand the impact of seismic unit thickness on our modelled sediment loading and overpressure generation, we apply a thickness for the Lower Unit of 73 m, representing the maximum estimated threshold for vertical resolution. This threshold was determined using a relationship between frequency, velocity, wavelength and resolution (Liner & McGilvery, 2019). For gypsum rock with compressional velocities of 5,700-5,800 m/s and dominant frequency of the seismic signal between 50 and 20 Hz, the theoretical thickness that can be resolved is estimated at 29-73 m. Our models show that during deposition of 73 m of LU gypsum 5.97-5.6 Ma, only low overpressure, of 0.4 MPa, is generated. Even if we assume erosion removed part of the LU gypsum and that a maximum thickness of 300 m may have been deposited 5.97-5.6 Ma (Stefano et al., 2010), then only low overpressure, of 0.95 MPa, is generated. We therefore show that for the Algero-Balearic basin at the base of the Emile Baudot Escarpment, sediment loading of the MSC Stage 1 (Lower Unit) 5.97-5.6 Ma played no role in overpressure increase above hydro fracturing and fluid release in the area.

| CONCLUSIONS
Thick low permeability evaporites from the Messinian Salinity Crisis have generated high overpressure within the evaporites and throughout pre-and post-Messinian sequences. The high overpressure within the evaporites would have been sufficient to hydro fracture them and generate vertical fluid release features. This study generates for the first time quantitative estimates of the time evolution of overpressure.
We completed laboratory measurements to constrain properties of evaporite minerals as input to a series of 1-D F I G U R E 1 3 Comparison between gypsum and anhydrite reaction and disequilibrium compaction as possible mechanisms explaining observed fluid escape features in the Algero-Balearic basin. (a) Seismic profile E12-SF 03 showing location of 1-D overpressure models, and interpreted horizons and faults. Inset showing the location of Algero-Balearic seismic profiles (black lines) and 1-D overpressure models evaluated in this study. Note that the Lower Unit (LU) is absent in this location. Green star shows the location of possible evaporite diagenesis and a fluid flow feature from Bertoni and Cartwright (2015). (b) Pressure and temperature phase diagram for gypsum-bassanite-anhydrite with dehydration boundaries by Klimchouk (1996;  disequilibrium compaction models. We conclude for the physical properties of Mediterranean evaporites that: • Evaporite porosities lower than 3% can become connected by cracks and/or dilatancy of grain boundaries allowing fluid flow. • Permeabilities of anhydrite and gypsum at different confining pressures range from 10 -17 to 10 -21 m 2 .
We used a 1-D disequilibrium compaction model to reconstruct fluid flow through time and to quantify the magnitude of overpressure generated in the Western Mediterranean basin. For the Liguro-Provençal basin and Algero-Balearic basin we conclude: • Permeability lower than 10 -18 m 2 can cause overpressure within Messinian evaporites. • Rapid sediment loading of low permeability Messinian evaporites inhibited vertical fluid flow causing high overpressure within pre-Messinian and Messinian sequences. • Rapid sediment loading caused sufficient overpressure to hydro fracture MSC evaporites. Hydro fracturing may have occurred during Stage 2 deposition of halite (Mobile Unit) from about 5.58-5.55 Ma in the Liguro-Provençal basin, and during Stage 3 deposition of Upper Gypsum from 5.55 to 5.33 Ma in the Algero-Balearic basin. • Fluid release features observed in seismic reflection data near the Emile Baudot escarpment of the Algero-Balearic basin, previously interpreted to be caused by gypsumanhydrite transformation, can also be explained by disequilibrium compaction-related hydro fracturing.