Investigation of the Inherent Variability of the Mediterranean Sea Under Contrasting Extreme Climatic Conditions

Sampatakaki, Angeliki; Zervakis, Vassilis; Mamoutos, Ioannis; Tragou, Elina; Gogou, Alexandra; Triantaphyllou, Maria; Skliris, Nikolaos

doi:10.3389/fmars.2021.656737

ORIGINAL RESEARCH article

Front. Mar. Sci., 24 August 2021

Sec. Physical Oceanography

Volume 8 - 2021 | https://doi.org/10.3389/fmars.2021.656737

This article is part of the Research Topic Impact of Deep Oceanic Processes on Circulation and Climate Variability: Examples from the Mediterranean Sea and the Global Ocean View all 17 articles

Investigation of the Inherent Variability of the Mediterranean Sea Under Contrasting Extreme Climatic Conditions

$\nAngeliki Sampatakaki,$ Angeliki Sampatakaki^1,2^*

Maria Triantaphyllou³

Nikolaos Skliris⁴

¹Department of Marine Sciences, University of the Aegean, Mytilene, Greece
²Hellenic Centre for Marine Research, Institute of Oceanography, Anavyssos, Greece
³Faculty of Geology and Geoenvironment, National and Kapodistrian University of Athens, Athens, Greece
⁴Ocean and Earth Science, University of Southampton, Southampton, United Kingdom

The internal variability of the thermohaline circulation of the Mediterranean Sea is examined under contrasting extreme thermal and mass atmospheric forcing conditions. Three millennium-long numerical simulation experiments were performed under: (a) the current climatology, (b) a strong buoyancy forcing (SBF) scenario due to cold and dry conditions resembling the Younger Dryas event, and (c) a weak buoyancy forcing (WBF) scenario due to S1a sapropel deposition-like conditions (warm and wet). To isolate the inherent variability of the system, independent of interannual atmospheric forcing variability, the latter was defined as a perpetual year pertinent to each experiment. Self-diagnosed heat and salt fluxes, consistent to sea-surface characteristics of the above periods, forced three millenium-long, relaxation-free numerical experiments. These simulations were preceded by initial spin-up periods. The inherent spatiotemporal variability of the Mediterranean Sea was analyzed using the empirical orthogonal function (EOF) and spectral analysis on the simulated density fields. Our results revealed that the Mediterranean Sea exhibits high sensitivity to climatic conditions, allowing its circulation to change from anti-estuarine (for the SBF scenario, leading to a buoyancy loss to the atmosphere) to estuarine (for the WBF scenario, corresponding to a buoyancy gain from the atmosphere). In all three experiments, the interannual and decennial variabilities dominate in upper layers, and the decennial variability dominates in the Gibraltar and Sicily Straits. Under current climatic conditions the first two EOF modes express only 60% of the density variability in the deep layers. This contribution exceeds 90% under more extreme conditions. Moreover, the first EOF modes correspond to a basin-wide in-phase variability of the deep layers under the reference and WBF conditions. During SBF conditions the first modes reveal a vertical buoyancy exchange between upper and deeper layers. The second EOF mode of deep waters under both extreme scenarios showed that the western and eastern basins exchange buoyancy in decennial (for the cold/dry) and interdecennial (for the warm/humid) timescales. The residence time of the Eastern Mediterranean deep water was diagnosed to be centennial, semicentennial, and intercentennial for the cases of current period, SBF, and WBF, respectively.

Introduction

The Mediterranean Sea is characterized as a concentration basin (Nielsen, 1912; Sverdrup et al., 1942). Atlantic water enters as a surface current through the Gibraltar Strait and travels eastward along the southern shores of the Mediterranean Sea, undergoing the transformation of its characteristics, whereas the discharge of the Mediterranean water to the Atlantic occurs via the intermediate and deepest layers. The overturning circulation pattern of the Mediterranean Sea is powered by the formation of deep waters in the Gulf of Lions (Stommel, 1972) for the western subbasin and in the Adriatic (Malanotte-Rizzoli et al., 1997) and Aegean Seas (Zervakis et al., 2000; Nittis et al., 2003) for the eastern one, while the Levantine intermediate water (LIW) has a key role to these formations (Lascaratos et al., 1993; Theocharis, 2009). Numerous three-dimensional (3D) numerical models have been developed to simulate the circulation characteristics of the Mediterranean Sea and their evolution throughout the instrumental period (e.g., Wu and Haines, 1996; Skliris et al., 2007; Somot et al., 2008; Oddo et al., 2009; Beuvier et al., 2010; Pinardi et al., 2015).

Nevertheless, the overturning circulation of the Mediterranean Sea was not always the same as is recorded today. In the past the Mediterranean basin occasionally exhibited periods of dramatic reduction of deep-water ventilation, which can be identified in the sedimentary record by the presence of sapropels. These dark organic-rich layers are found throughout the Eastern Mediterranean, the Tyrrhenian Sea, and some parts of the Western Mediterranean and demonstrate wet and warm climatic conditions dominating the basin. The sedimentary record shows their reoccurrence in the past 7 million years while, in the last 450,000 years (ka), a minimum of 11 sapropel layers were formed (e.g., Rohling, 1994; Rohling et al., 2015; Andersen et al., 2018). The warm periods of sapropel deposition have been linked to Milankovitch cycles and precession minima (e.g., Rossignol-Strick, 1985; Lourens et al., 1996). The last sapropel deposition S1 was recorded during the Holocene Climatic Optimum (~10–6 ka BP). Both S1 sapropel sublayers (S1a and S1b) correspond to the warmest and most humid Holocene intervals for the Mediterranean, respectively (Gogou et al., 2007; Triantaphyllou et al., 2009a,b, 2016; Geraga et al., 2010; Triantaphyllou, 2014).

In contrast, the cool event of Younger Dryas was recorded just before the onset of the Holocene. Younger Dryas (12.9–11.7 ka BP) was the driest interval of the past 20,000 years (Alley, 2000; Carlson, 2013), and its impact has been recorded all over the Mediterranean Sea (e.g., Kotthoff et al., 2008).

The abovementioned changes recorded in the Mediterranean Sea sediments are linked with atmospheric forcing variability. Climatic teleconnections have been proposed to be linked with the Younger Dryas event (Cacho et al., 2001; Gogou et al., 2007; Marino et al., 2009). For the sapropel deposition intervals, data from the different sites indicate that, during these warm periods, the Eastern Mediterranean Sea was influenced by increased river discharges due to an enhanced monsoon activity over North Africa (Rossignol-Strick, 1985; Almogi-Labin et al., 2009; Athanasiou et al., 2015, 2017; Rohling et al., 2015) and a higher amount of rainfall around the Northern borderlands (Gogou et al., 2007; Kotthoff et al., 2008; Triantaphyllou et al., 2016).

In the last postglacial period, sediment records indicate that climatic alterations over the Mediterranean Sea between the warm and cold intervals may have modified the sea circulation in the subbasins and the Mediterranean outflow as well (Llave et al., 2006; Toucanne et al., 2012; El Frihmat et al., 2014; Jiménez-Espejo et al., 2015; Rohling et al., 2015; Tachikawa et al., 2015; Dubois-Dauphin et al., 2017). In the postglacial period, most 3D numerical paleo-modeling simulations referred to the S1 sapropel deposition period (Myers et al., 1998; Myers and Rohling, 2000; Myers, 2002; Meijer and Tuenter, 2007; Adloff et al., 2011; Grimm et al., 2015). Several simulations of the Mediterranean circulation during sapropel formation period converge to the finding that the basin continued to function as an overall concentration basin, at least regarding the Gibraltar exchange (Myers et al., 1998; Myers, 2002; Meijer and Tuenter, 2007; Adloff et al., 2011), a conclusion that is not necessarily supported by the sediment record (e.g., Thunell and Williams, 1989). We are not aware of efforts to numerically simulate the 3D Mediterranean circulation during the Younger Dryas period. In any case, all paleoceanographic simulations of the Mediterranean aim to reconstruct the long-term mean circulation and hydrographic characteristics of the basin. Of course, the lack of instrumental climatological information during the paleo periods, combined with the low sedimentation rate in the basin, constitutes a limiting factor to the simulations, especially in regard to the interannual to interdecennial variability.

However, both the observational record and the simulations of the basin-wide circulation converge to the conclusion that, in the last 50 years, the Mediterranean has exhibited a very energetic decennial variability, even in its deep and bottom layers (e.g., Robinson et al., 2001; Roether et al., 2007; Schroeder et al., 2016). In the 1980s, the interannual variability of the Mediterranean deep waters was first evidenced as a long-term linear trend in the western basin (Lacombe et al., 1985). The lack of frequent measurements in the Eastern Mediterranean until the mid-1980s delayed the revelation of any variability until the discovery of a major shift of dense water formation from the Adriatic to the Aegean Sea that took place between 1987 and 1993 and is known since then as the Eastern Mediterranean Transient (EMT) (Roether et al., 2007). Recent numerical hindcasts of the Mediterranean circulation spanning several decades have successfully reproduced a significant interdecennial variability of the deeper layers and dense-water formation processes (Pinardi et al., 2015; Waldman et al., 2018).

The interdecennial variability observed during the instrumental period has been attributed mainly to atmospheric forcing (e.g., Béranger et al., 2010; Incarbona et al., 2016; Somot et al., 2018); however, a significant part of the variability has been attributed to the internal dynamics of the basin. Waldman et al. (2018), analyzing a 31-year-long hindcast, showed that about 28% of the deep convection events can be attributed to the chaotic behavior of the system. Several mechanisms have been identified as a contributor to this chaotic behavior, of which the most prominent one—for the eastern Mediterranean—may be the Adriatic-Ionian Bimodal Oscillating System (BiOS; Gačić et al., 2010; Theocharis et al., 2014). Several theoretical and experimental studies have shown that internal dynamics may be responsible for the observed variability of the system, attributing the interannual–interdecennial variability of the overturning circulation and observed alteration of deep-water formation sites to internal exchange processes, causing a stochastic resonance-type behavior (Ashkenazy et al., 2012; Crisciani and Mosetti, 2016; Rubino et al., 2020).

Thus, several Mediterranean circulation analyses converge to the conclusion that interior dynamics are responsible for the interannual–interdecennial variability of its deeper layers in the instrumental period. However, as the mean overturning circulation and the corresponding density stratification of the Mediterranean Sea can vary over geological timescales, how would this affect the inherent variability of the basin in the decennial to centennial periods?

This study aimed to investigate the internal/inherent variability of the Mediterranean Sea under extreme cases of thermohaline functioning. To achieve this aim, three long-term numerical experiments were designed and performed to simulate the circulation of the basin under the climatic conditions typical of (a) the current (instrumental) period, hereafter called as the “reference” period, (b) the Younger Dryas-like cold and dry period, hereafter called as the “strong buoyancy forcing” (SBF) period, and (c) the warm and humid period of the S1a sapropel deposition-like interval, hereafter called as the “weak buoyancy forcing” (WBF). As shown in the following sections, SBF corresponds to a strong buoyancy loss from the sea to the atmosphere, whereas WBF corresponds to a buoyancy gain. Long-term simulations, forced by a perpetual year, are performed until “equilibria” states are reached. The heat and freshwater fluxes are diagnosed and then applied to the model experiments for a period of 1,000 more years. The internal/inherent variability of the basin under the three climatic scenarios was examined using the empirical orthogonal function (EOF) and spectral analyses.

Materials and Methods

Experimental Setup of the Model

Baseline Model (Reference Period)

The Princeton Ocean Model (POM, Blumberg and Mellor, 1987) was used with a configuration for the Mediterranean Sea similar to that carried out in the study by Skliris et al. (2007). The model uses 25 sigma layers and its horizontal resolution is 1/4° × 1/4°. The western boundary of the model is located in the Atlantic Ocean, and the inflow through the boundary, temperature, and salinity are prescribed from the Mediterranean Oceanic Database (MODB) seasonal climatology (see details in Skliris and Lascaratos, 2004; Skliris et al., 2007). The model uses artificial surface fluxes based on the bulk formulae as mentioned in Skliris et al. (2007), i.e., Rosati and Miyakoda (1988), exploiting the empirical formula of Reed (1977) for cloud attenuation and of Bignami et al. (1995) for long-wave radiation, while the calculation of sensible and latent heat fluxes follows the neutral Budyko (1963) scheme. The surface boundary conditions of the heat and salt fluxes are initially relaxed to the MODB (Brasseur et al., 1996) seasonal climatology (Equations 1 and 2). Precipitation is taken from the monthly data set as proposed by Jaeger (1976). The model is driven by a perpetual year atmospheric forcing of the monthly averages obtained, from 1973 to 1993, using 6-h European Center for Medium Range Weather Forecasts (ECMWF) reanalysis data and is initialized from the MODB spring temperature and salinity fields (for more details about the model, see Skliris et al., 2007).

The surface boundary conditions of the heat flux F_H and the salt flux F_S at the air–sea interface for the spin-up are given by

\begin{array}{l} F_{H} = Q_{N} + C_{T} (T^{c} - T_{1}) & (1) \end{array}

\begin{array}{l} F_{s} = (E - P) S_{1} + C_{S} (S^{c} - S_{1}) & (2) \end{array}

where Q_N is the net heat flux at the air–sea interface, P is the precipitation and E is the evaporation, T^c and S^c are the seasonally averaged climatological sea surface temperature (SST) and salinity from the MODB (Brasseur et al., 1996), and T₁ and S₁ are the temperature and salinity of the top-level of the model, respectively.

As we would like to investigate long-term internal dynamics without damping terms, the model was set up to run with no internal relaxation, and only the following modifications were applied to the model configuration mentioned in the study of Skliris et al. (2007): the outflow of River Po into the Adriatic was added to the original model as a source point (according to Skliris et al., 2007 river runoff setup), with the average discharges of river being equal to 2,000 m³s⁻¹, corresponding to the average discharges of rivers into the Adriatic (Ludwig et al., 2009) and a sinusoidal variation with a range of 500 m³s⁻¹ and two peaks a year in spring and autumn. For the Dardanelles exchange, the salinity of the inflow water and the outflow water was set to 29.6 and 38.9 throughout the year, respectively (Ünlüata et al., 1990). For the Nile River, discharges prior to the Aswan Dam construction are considered as described by Skliris et al. (2007), in all the experiments.

Then, the correction coefficients of the surface heat and salinity fluxes are estimated and redesignated as follows: the root-mean-square error (RMS) index was used to compare the surface seasonal values of temperature and salinity fields of the Mediterranean Sea, as produced by the model runs with the modified correction coefficients after a 40 year spin-up, with the respective climatological values. The RMS index was applied again on the average seasonal values produced by the model to an average depth of 300 m with the corresponding climatology values. The average of the results of the indicator on temperature and salinity in the surface and at an average depth of 300 m gave the final choice of correction coefficients (C_T, C_S).

The baseline model was spun-up for 320 years and reached its steady state from the first 160 years run.

Baseline Model (Reference Period): Application of Diagnosed Fluxes

According to previous studies (Myers and Haines, 2000; Pisacane et al., 2006), the diagnosed fluxes are applied to the model and the internal variability of the basin is examined.

The diagnosed monthly heat and freshwater fluxes were estimated after the spin-up period, and the model ran for another 1,000 years with these fixed fluxes.

The monthly diagnosed freshwater fluxes F_F were calculated by

\begin{array}{l} F_{F} = E - P + \frac{C_{S} (S^{c} - S_{m})}{S_{m}} D & (3) \end{array}

where S_m is the monthly averaged surface salinity as calculated from the last year of the spin-up of the model and D is the thickness of the first σ-layer of the model.

The monthly diagnosed surface heat fluxes $F_{H}^{^{'}}$ were calculated by

\begin{array}{l} F_{H}^{^{'}} = Q_{N} + C_{T} (T^{c} - T_{m}) & (4) \end{array}

where T_m is the monthly average temperature of the top-level of the model as calculated from the final year of the spin-up period.

The application of the diagnosed heat and salt fluxes as boundary conditions allowed the model to run free of any relaxation terms, leading to realistic SST and salinity fields. The relaxation terms of the model were removed to ensure that any observed internal variability of the basin would not be a product of artificial nudging (Weaver and Hughes, 1992).

Baseline (Reference Period): Model Validation

The reference experiment is evaluated based on the comparison of hydrographic characteristics in the interior of the basin between a recent climatology [from the available temperature and salinity profiles from 1900 to 2015 (Iona et al., 2018)] and the last 500 years of the simulation.

The evaluation of the simulations has been conducted via a comparison of the mean vertical profiles and the combined θ/S diagrams (Figure 1). Overall, the simulation exhibits a faithful reproduction of the vertical structure of the water column and an impressive agreement in the θ/S diagrams despite the presence of significant biases in the deeper layers. The simulation possibly underestimates the temperature by about 1°C and the salinity by about 0.3 psu. Such biases are not acceptable for operational simulations employing data assimilation or process studies using relaxation within the water column; however, we consider them acceptable in the framework of the present simulation, especially considering the agreement in the water masses as displayed in the θ/S diagram.

FIGURE 1

Figure 1. Comparison of simulated (red) and observed (black) salinity (top) and temperature (middle) profiles for the Western (A,C) and Eastern (B,D) basins. Both mean profiles (solid lines) and mean ± one standard deviation (STD) (dashed lines) over the complete data sets are presented. The comparison of the mean values is extended to the corresponding θ/S diagrams in the bottom panels (E,F).

An indirect way to evaluate the overall performance of the simulation regarding deep water formation processes is via the overturning stream function, which is presented in the Results section of this study, in comparison to the results of the extreme forcing simulations.

Setting Up the Model for Extreme Climatic Conditions

To determine the extreme conditions required for the two experiments, we selected to simulate two contrasting cases in regard to the Mediterranean hydrographic conditions, at least as witnessed immediately before or within the Holocene period in the geological record. The period selected as a representative of the least stratified conditions was the Younger Dryas period (12,900–11,600 years BP, Cheng et al., 2020) during which the basin suffered a great buoyancy loss due to predominant cold and dry conditions. The period chosen as a representative of the highest stratification in the Mediterranean was the S1a sapropel deposition period (10,200–8,000 years BP, Triantaphyllou et al., 2016). To approach such conditions in the simulations, the forcing fields of the baseline model were modified as described below, based on the reconstructed sea surface data of marine cores and the insolation of each simulating period.

Determination of the Paleo-Modeling Adjustments

First of all, spatially averaged temperature (at) and salinity (as) ratios were used to assess the heat and freshwater fluxes from the SST and salinity fields. The coefficients at and as were calculated for each extreme period (Table 1) to provide the necessary sea surface conditions that will lead to the related forcing fluxes:

\begin{array}{l} a t_{s} = \frac{T_{12}}{T_{0}} a t_{w} = \frac{T_{9.1}}{T_{0}} \end{array}

where T₁₂ and T_9.1 are the averages of the alkenone-reconstructed SST data from marine cores 434G, 293G (Rodrigo-Gámiz et al., 2014), 161-977A (Martrat et al., 2007), GeoB7702-3 (Castañeda et al., 2010), 9501 (Almogi-Labin et al., 2009), MNB-3 (Gogou et al., 2007), and SL-152 (Emeis unpublished data; Gogou et al., 2016) from the Mediterranean basin (Figure 2) for the periods 12,300–11,700 cal BP for the Younger Dryas and 9,400–8,800 cal BP for the S1a sapropel deposition intervals, respectively. These intervals were selected because they represent “time steps,” where the conditions of Younger Dryas and S1a sapropel formation are fully developed in all core sites of the Mediterranean basin. The abovementioned cores were selected as the sources of reconstructed SST information because they were available for both the past periods of interest and represent the temperature in both the eastern and western subregions. T₀ is the average spring value of the Mediterranean SST from the modern climatology (MODB) because the alkenone-reconstructed core data express the spring or winter SST (Emeis et al., 2000).

TABLE 1

Table 1. Values of the correction coefficients applied to the reference model to reach extreme case sea surface conditions similar to paleoceanographic conditions.

FIGURE 2

Figure 2. Marine core sites and transects used in the present study.

Similar to at_s and at_w, as_s and as_w were defined as the ratios of average near-surface salinities recorded in 12 and 9.1 ky BP, respectively (based on the reconstructions by Emeis et al., 2000 with a ratio of 0.25‰) to the mean of modern climatology (MODB). Here, it must be noted that the abovementioned coefficients are numerical indices not carrying any spatial information. These coefficients are applied as a spatially homogeneous correction to the current-era fields (Equations 7 and 8), and the spatial variability of the paleoclimatic simulations is generated by the model dynamics. Thus, the spatial patterns of the SST and salinity are not imposed by the corresponding core reproductions, which have been used to determine a horizontally homogeneous temporal change of the boundary conditions.

To reproduce the seasonality of 12 and 9.1 ky BP periods, the monthly values of the heat fluxes were transformed using the factors $r_{s} = \frac{Q_{S_{s}}}{Q_{S_{0}}}$ , $r_{w} = \frac{Q_{S_{w}}}{Q_{S_{0}}}$ , where the terms r_s and r_w (Figure 3A) are given by the ratio of monthly average insolation at 38^o latitude in 12 and 9.1 ky BP (Q_{S_s} and Q_{S_w}, respectively) to the corresponding insolation during the instrumental period (Q_S₀), based on the solar radiation calculations of Berger (1978) and Berger and Loutre (1991), exploiting the algorithm of Huybers and Eisenman (2006). Note that hereafter the subscript “s” denotes the experiment held under SBF (strong buoyancy loss and enhanced mixing conditions), whereas the subscript “w” signifies WBF (buoyancy gain and highly stratified conditions).

FIGURE 3

Figure 3. (A) The term r gives the ratio of monthly average insolation at 38° latitude in 12 and 9.1 ky BP to the corresponding insolation during the instrumental period, (B) the net heat flux with the applied diagnosed fluxes for each experiment. The green line corresponds to the reference period, the blue to strong buoyancy forcing (SBF; Younger Dryas-like conditions), and the red one to the weak buoyancy forcing (WBF; S1a sapropel deposition-like conditions).

It should be noted here that, for the extreme climatic simulations, the exchanges through the Dardanelles Straits were not taken into consideration. This may not be consistent with the current conditions but is more compatible with the past periods that have been used to obtain the heat and freshwater forcing fields, as the water exchange between the Mediterranean Sea and the Black Sea was fully established only after the S1a sapropel deposition (~8.7 ka BP, Sperling et al., 2003; 7.5 ka BP, Triantaphyllou et al., 2016) and the Black Sea was not a major freshwater contributor for the S1 deposition (e.g., Rohling et al., 2015). The momentum forcing and bathymetry remain unchanged in the context of hereby described experiments.

SBF: Model Set Up and Run

To implement the experiment held under SBF, simulating conditions similar to those of the Younger Dryas period, the monthly fields of heat flux (Q_N) and net evaporation (E_N) of the current period were transformed as:

\begin{array}{l} Q_{N_{s}} = r_{s} (Q_{S} - \frac{Q_{U}}{a}) & (5) \end{array}

\begin{array}{l} E_{N_{s}} = k r_{s} E_{N} & (6) \end{array}

where Q_{N_s} and E_{N_s} are the estimated monthly values of net heat flux and evaporation for the SBF period, respectively, r_s is the carrier of the seasonality difference between modern and Younger Dryas periods, Q_S is the solar radiation flux, and Q_U is the net heat flux from the other three heat flux components. The parameter k (Equation 6) is defined as the ratio of the net evaporation during the paleoceanographic period under investigation, over the current values. For the case of the Younger Dryas simulation, k has been empirically estimated to 1.8 after several numerical experiments. As the sea surface salinities of this period are comparable with the current values of the Red Sea, we used the net evaporation rate of the latter as a reference point to empirically estimate the parameter k. After spinning up the model, the mean sea surface salinity was compared to the mean core data reconstructions of the Younger Dryas, and the most suitable parameter value was selected. The increase in net evaporation E_N is equivalent to a lesser increase of evaporation E combined with a concurrent decrease of precipitation P or a change of E−P. The parameter α (Equation 5) was set equal to at_s and expresses the extra heat loss during Younger Dryas climatic conditions. For the extreme climatic simulations, it was assumed that the ratio of the heat fluxes is proportionate with the ratio of the SSTs and that is expressed by the parameter a.

Thus, the formulae for the surface boundary conditions (Equations 1 and 2) were used to determine the heat flux F_{H_s} and the salt flux F_{S_s} for the SBF experiment event at the air–sea interface and were modified as:

\begin{array}{l} F_{H_{s}} = Q_{N_{s}} + C_{T} (a t_{s} T^{c} - T_{1}) & (7) \end{array}

\begin{array}{l} F_{S_{s}} = a s_{s} E_{N_{s}} S_{1} + C_{S} (a s_{s} S^{c} - S_{1}) & (8) \end{array}

The modern-climatology values T^c and S^c were multiplied by the ratios at_s and as_s, respectively, to simulate mixed-period surface values, similar to those found in the Younger Dryas period. The model was initialized from the current climatology and then spun up for a total of 400 years while reaching a steady state after 160 years.

Similar to the baseline model runs, after the completion of the initial spin-up, the salt and heat monthly diagnosed fluxes (Figure 3B) were estimated, and the model was turned toward the diagnosed atmospheric fluxes and run for another 1,000 years.

WBF: Model Set Up and Run

According to the SBF simulation experiment, the circulation of the Mediterranean Sea was simulated under WBF (actually, weak buoyancy gain experiment). The final state of the Mediterranean, after the SBF experiment, was defined as the initial conditions of the WBF simulation. Given that the basin obtained the characteristics of a well-mixed basin after a 400-year-long period of spin-up run, the new forcing conditions were applied after gradually reducing the net evaporation and increasing the net heat flux over the basin until they reached the WBF experiment forcing values.

Similar to the calculations for the SBF experiment, the monthly fields of Q_{N_w} and E_{N_w} for the stratified conditions were estimated by replacing the term r_s with r_w in Equations (5) and (6).

For the surface fluxes (Equations 7 and 8), the coefficients at_s and as_s were replaced with at_w and as_w, respectively. In Equation (5) the parameter α was set equal to 1. Furthermore, in regard to the application of Equation (6), numerical experiments showed that the abrupt replacement of the parameter k = 1.8 (as used in the SBF experiment) to k = 0.8 (a realistic estimate of the reduced evaporation over the Mediterranean during the S1a deposition period (Myers, 2002; Meijer and Tuenter, 2007) caused numerical instability and blow up of the simulation. To reach the S1a period-like conditions, it was necessary to vary the evaporation ratio k as follows: For the first 80 years, the value k = 1 was used. For the next 200 running years, the coefficient k was set to 0.5, indicating half of the current evaporation rate. Finally, the model ran for 520 more years with k = 0.8. The model spin-up lasted 800 years totally. A steady state was reached after 120 years of the last k modification (the 440th year of spin-up).

Then, in accordance with previous experiments and to enable the study of inherent centennial variability, the model was switched to the diagnosed boundary conditions and ran for another 1,000 years.

The seasonal cycle of the diagnosed surface heat fluxes for each experiment is shown in Figure 3B. The mean annual net heat flux in the Gibraltar Strait for each period was calculated to be +1.0 Wm⁻² for the current period (a positive sign denoting eastward net heat flux, which corresponds to a surface heat loss), +5.5 Wm⁻² for the strong buoyancy, and −6.8 Wm² for the WBF.

Moreover, regarding the net freshwater fluxes applied for the experiments, the ratios of the diagnosed net E–P fluxes to the current value were calculated to be 2.8 and 0.63 for the SBF and WBF experiments respectively (computed on a mean annual basis).

Analysis of the Results

The basic characteristics of the mean circulation reproduced by each experiment and the interdecennial/intercentennial variability of the basin are described as follows.

Deep-Water Residence Time

The residence time of Eastern Mediterranean deep water was defined as the mean annual ratio V_M/R of the last 100 years of the model running with the diagnosed fluxes, where V_M is the volume of the Eastern Mediterranean Sea waters below 300 m (the depth below which the Aegean Sea is considered to influence the whole eastern basin) and R is the inflow of dense waters into the Eastern Mediterranean through the corresponding straits connecting them to the deep water sources of the basin, as represented by a double integral:

\begin{array}{l} R = \int_{l_{1}}^{l_{2}} \int_{- h (x)}^{0} U d z d x \end{array}

where

• For the reference climate simulation and the SBF experiment, U is the velocity of the incoming to the eastern basin dense waters through the Otranto Strait (for Adriatic-originated dense waters) and the Cretan Arc Straits (for Aegean-originated waters).

• For the WBF runs, U is the velocity of the incoming dense water to the eastern basin via the Sicily Strait since (as it is shown below) our simulation revealed that, during that period, the Western Mediterranean imports the buoyancy from the Eastern basin, in complete contrast to its current function.

Buoyancy Calculations

To investigate the thermohaline functioning of the Mediterranean Sea (concentration or dilution basin characteristics) under different extreme climatic conditions, the integrated lateral advective buoyancy fluxes through the Gibraltar and Sicily Straits were calculated for the last 500 years of simulations with the diagnosed fluxes by

\begin{array}{l} F = \int_{l_{1}}^{l_{2}} \int_{- h (x)}^{0} - g (ρ - ρ_{o}) u d z d x \end{array}

where l₂ and l₁ are the lateral limits of the transect, h(x) is the bottom depth along the transect, ρ is the density at each point of the transect, ρ_o is the mean density of the corresponding basin, and u is the velocity component normal to the transect at each point.

Variability Analysis

To investigate the inherent spatial and temporal variability of the Mediterranean Sea, the EOF analysis was applied to the mean annual density fields of the last 500 years running with the diagnosed fluxes for each experiment, at different depths (30, 100, 200, 400, 800, 1,200, and 2,000 m), breaking each field data into “modes of variability.” Following the terminology of Björnsson and Venegas (2001), the spatial patterns that are found are referred to in the text as EOFs and the time series as expansion coefficients. The importance of each EOF pattern is referred to as a contribution.

In all experiments, eight modes were found to represent 99.5% or more of the variability for all depths. For each mode, the power spectrum of the expansion coefficient was calculated to provide a description of the temporal variability of the mode. The spectra presented here calculated with seven sequential segments at a 50% overlay. When spectral shapes suggested that most of the energy lies in centennial variability, spectra were also computed using only three segments at a 50% overlay to provide the information of the lowest possible frequencies. However, for the purpose of clarity and consistency, the spectral results hereby presented are calculated using seven segments with a length of 125 years per segment.

To identify the most dominant frequencies of variability, the spectra of the simulated variability are compared with the theoretical continuous red-noise spectra of the same variance. To calculate the latter, first, the decorrelation time scale of each EOF expansion coefficient is computed, and then, the red-noise spectrum of a process exhibiting this decorrelation timescale is estimated (Hasselmann, 1976).

Finally, a cross-correlation between the expansion coefficients of modes at different depths was used to identify and cluster the layers with similar variability.

Results

Overturning Functioning of the Mediterranean Sea Under Different Climatic Conditions

The simulation experiments produced three different thermohaline circulation patterns (Figures 4, 5).

FIGURE 4

Figure 4. The Mediterranean zonal overturning transport (estimated via streamfunction computation) for (A) the reference period, (B) the SBF, and (C) the WBF, as produced by the model experiments under the diagnosed fluxes. The transport is given in Sverdrups.

FIGURE 5

Figure 5. Functioning of the Mediterranean Sea under the climate conditions of the three experiments: reference (top), SBF (middle), and WBF (lower panels). (A,C,E) reveal the overall functioning of the Mediterranean Sea as witnessed by the exchange at Gibraltar Strait. The pink arrow denotes the inflow of Atlantic water, whereas the light blue arrow corresponds to the outflow from the Mediterranean. The density distribution along the transect of Figure 2 along the Mediterranean is shown in (B,D,F). In all figures, the values are annual means, averaged over the last 100 years run with the diagnosed fluxes for each experiment.

The mean millennial overturning circulation of the Mediterranean Sea for each numerical experiment reveals that both the basins behave as concentration basins in the reference and SBF experiments, whereas the eastern Mediterranean clearly becomes a dilution basin under the WBF conditions. A comparison of the zonal—vertical transport (computed as streamfunction fields) of the reference period (Figure 4A) to the corresponding calculations from the instrumental period (Adloff et al., 2011, 2015) suggests that the western basin shows significant differences while the eastern Mediterranean exhibits a very similar overturning circulation: the present simulation produces more waters that reach deeper layers in comparison to the results of (Adloff et al., 2011, 2015). Given the necessary compromises (primarily in terms of model resolution and the temporal average of atmospheric forcing) accepted to sustain a millennial simulation and independence from external variability and the good results from the comparison with climatology (see the section “Baseline (Reference Period): Model Validation”), we consider the results acceptable for our purpose.

As expected, the SBF simulation (Figure 4B) produces more dense water in the eastern basin compared to the reference experiment (Figure 4A), whereas the results in the western basin are rather surprising, due to the production of dense water extending over a larger area of the basin; however, most of the produced water does not reach the maximum depths of the basin, possibly due to a higher stratification in the deeper layers.

In the WBF experiment, the eastern Mediterranean switches to a dilution basin, exporting surface waters to the west (Figure 4C). The western basin exhibits some weak dense-water formation cells at the longitude of the Gulf of Lions; however, at both connecting Straits (Gibraltar and Sicily) the streamfunction sign suggests an exchange typical of dilution basins.

The thermohaline functioning of the Mediterranean basin as a whole can be assessed by examining the exchange in the Strait of Gibraltar (Figure 5). The reference simulation produces (as expected) a two-layer exchange, with the surface inflow of Atlantic waters and the outflow of Mediterranean waters below the depth of 100 m. This corresponds to a striking but consistent three-layer buoyancy exchange, with buoyancy import at the top and bottom layers and a thin layer of weak buoyancy export from 80 to 120 m, at the region of mixing between the two water masses (Figure 5A). Under SBF, the exchange is intensified and the interface between the two layers deepens to 130 m depth (as the net import from the Atlantic has to balance the enhanced net evaporation). Atlantic water flows into the basin through a layer extending from the surface to a depth of about 130m, whereas below that depth dense Mediterranean water is exported to the Atlantic Ocean. Under these conditions, the buoyancy-exporting layer is diminished at the interface depth. A comparison of the corresponding density transects along a zonal section of the Mediterranean reveals overall denser waters in both the basins during both experiments, and a zonal density gradient with denser waters in the eastern basin. In contrast, for the WBF (during a warm and humid climate), light waters and strong stratification dominate in the Eastern Mediterranean, and dense water is exported from the western basin to the eastern one (Figure 5F). During this period, the thermohaline functioning of the whole basin is reversed and the Mediterranean Sea acts as an estuarine basin, exporting upper and intermediate waters to the Atlantic (Figures 5E,F). Sluggish intermediate Mediterranean water from 30 to 180 m depths moves toward the Atlantic, which feeds the Mediterranean basin through deeper layers. The buoyancy export takes place through the surface and near-bottom layers, whereas the sluggish, thick intermediate layer is a net but an insignificant buoyancy importer.

Moreover, according to the calculation described in section Deep-Water Residence Time, the residence time of the Eastern Mediterranean deep waters was calculated to be centennial (~100 years), semicentennial (~55 years), and intercentennial (semi-millennium, ~540 years), for reference climate conditions, the SBF experiment, and the WBF simulation, respectively.

Having described the millennial-scale climate of the basin for the three experiments, we can proceed to examine the inherent variability of the circulation at the interannual and intercentennial periods. Firstly, we will examine the variability of the thermohaline functioning of the whole basin, as evidenced by the Gibraltar exchanges (section Buoyancy and Its Variability in the Gibraltar and Sicily Straits). After that, the analysis is extended to the spatiotemporal variability of the interior of the basin, as revealed by the EOF analysis of the density fields (section EOF Analysis).

Buoyancy and Its Variability in the Gibraltar and Sicily Straits

Integrated buoyancy fluxes through the Gibraltar and Sicily Straits for the last 500 years of simulations using the diagnosed fluxes are presented in Figures 6, 7. To reveal the inherent variability of the exchange (and thus the corresponding variability of the functioning of the basin), variance-preserving spectra of the buoyancy fluxes were calculated for each period.

FIGURE 6

Figure 6. Transect-integrated buoyancy fluxes (A,C,E) and variance-preserving spectra (B,D,F) through the Gibraltar Strait of the last 500 years of simulations with the diagnosed fluxes for the reference (A,B), SBF (C,D) and WBF (E,F) experiments. The dashed gray lines surrounding the mean spectral estimate correspond to the 90% confidence intervals. For comparison, the corresponding normalized red noise of the same variance is overlaid (red lines). Positive values correspond to eastward buoyancy flux.

FIGURE 7

Figure 7. (A–F) Transect-integrated buoyancy flux through the Sicily Strait (as in Figure 6). Positive values correspond to a northward buoyancy flux.

According to the traditional oceanographic notation, northward and eastward directions are defined as positive to comprehend the direction of buoyancy flux through the Straits.

For the Gibraltar Strait (flow along the east–west axis), the buoyancy exchange corresponds to an inflow for the reference simulation (mean of 0.88 × 10⁷ N/s, Figure 6A) and also an even stronger inflow for the SBF experiment (mean 3.74 × 10⁷ N/s, Figure 6C); for the WBF simulation, the negative values indicate the net buoyancy export from the Mediterranean Sea (mean −1.04 × 10⁷ N/s, Figure 6E), i.e., the import of dense waters from the Atlantic. These calculations suggest that the Mediterranean Sea shows two “opposite” phases between the experiments. It behaves as a concentration basin under the reference and SBF conditions and as a dilution basin during the WBF experiment.

For the Sicily Strait, where the flow is mainly meridional, the negative buoyancy flux values observed during the reference period (mean −0.67 × 10⁷ N/s, Figure 7A) and SBF (mean −2.59 × 10⁷ N/s, Figure 7C) simulations indicate that dense water from the eastern basin is exported to the western. The positive flux during the WBF period (mean 4.84 × 10⁶ N/s, Figure 7E) shows an opposite functioning. The Eastern Mediterranean under these climatic conditions functions also as a dilution basin in agreement with the analysis mentioned in section Overturning Functioning of the Mediterranean Sea Under Different Climatic Conditions.

The long-term simulations revealed that the integrated buoyancy fluxes at the Straits oscillate on different timescales, but in all the experiments, a decennial variability dominates in both the Gibraltar and Sicily Straits (Figures 6B,D,F, 7B,D,F).

The variance-preserving spectra of the buoyancy flux through the Gibraltar Strait (Figures 6B,D,F) reveals that the variability is mainly intercentennial (peaking at a period of about 20–25 years), whereas interannual periods are responsible for a significantly smaller part of the variance. We note a relatively higher contribution of interannual variability (at periods of around 5 years) in the SBF experiment compared to the reference period and the diminishing contribution of interdecennial variability during the WBF experiment.

In the Sicily Strait, the interannual and interdecennial variabilities dominate the variance of the buoyancy exchange (Figures 7B,D,F). For all experiments, the spectral peaks are concentrated in the periods between 10 and 25 years, whereas a secondary peak appears at a period of about 40 years for the reference period. As in the case of Gibraltar, the interdecennial periods do not contribute much to the variance of the buoyancy exchange.

EOF Analysis

There is probably no better way to describe the mean spatiotemporal variability of a basin than EOF analysis. The eight lowest modes of spatiotemporal variability obtained from the EOF analysis in seven depths characterizing the shallow (30 and 100 m), intermediate (200 and 400 m), and deeper layers (800, 1,200, and 2,000 m) of the Mediterranean, account for >99.5% of the density variability. In the reference period, the variability of density is observed to be a synthesis of several modes with important contribution of each mode even at deeper layers—the first two EOF modes at 1,200-m depth express only 60% of the variability. In contrast, the density variability of the deeper layers at extreme climatic simulations can be described by only the first two EOF modes as they represent more than the 90% of the variability (Figure 8A).

FIGURE 8

Figure 8. (A) Empirical orthogonal function (EOF) mode contribution and (B) qualitative characterization of the dominant temporal variability of the first eight EOF modes of density variability at all depths for the three experiments (reference buoyancy forcing, SBF, and WBF).

Qualitative characterization of the more pronounced oscillation scales, as produced by the variance-preserving spectra of the expansion coefficients, for all depths and for the three experiments is shown in Figure 8B. Interannual variability refers to the main spectral peaks in the period between 2 and 9 years, the decennial between 10 and 20 years, the interdecennial between 21 and 80 years, the centennial between 80 and 120, and the intercentennial over 120 years. In all the experiments, the interannual and decennial variabilities dominate in upper layers. The deeper layers of the current period show also interdecennial and centennial variability. Considering that the variability at all depths is a contribution of several EOF modes during the current period, it is shown that the variability in upper layers is a synthesis of interannual and decennial modes and in the deeper layers the variability is more a synthesis of centennial and decennial/interdecennial oscillations (Figures 8A,B). In case of SBF conditions, the upper layers are dominated by the decennial and interannual variabilities, whereas in deeper layers the first two EOF modes exhibit intercentennial and interdecennial oscillations, respectively. Higher EOF modes with a low contribution to deeper layers show a decennial variability. For the WBF simulation, the modes of the upper layers show mostly an interannual variability. In deeper layers, the first two modes express intercentennial and centennial oscillations, respectively. Higher modes with a very low contribution to the variability have an interannual and a decennial variance.

Spatial and Temporal EOF Variability and Clustering

In each experiment, expansion coefficients of the modes from the different depths were subjected to a cross-correlation analysis to identify modes representing vertically coherent motions. Finding the depths that exhibit a similar temporal variability (correlation > 95%) leads to the identification of respective horizontal EOF modes representing vertically coherent motions.

In the following sections, the most dominant EOF modes exhibiting vertically coherent variability are presented.

Reference Period

In the deeper layers of the reference period, the variability is dominated by interdecennial periods of density variability with a maximum of ~63 years. This variability is recorded in the second and first mode at the layers of 800 m depth and the deeper layers, respectively (Figure 9) and corresponds to an oscillation of almost the entire Mediterranean Sea in phase agreement between the depths of 800 and 2,000 m. The contribution of this variability ranges from 24.8 to 51.8% in the upper and deeper layers, respectively. This mode of oscillation exhibits its peak variability in the Western Mediterranean, with reduced amplitudes in the Eastern basin, and probably reflects the dense-water formation processes, which are more energetic and exhibit great variability in the western basin. This variability can also be associated with the renewal time of the deep water of the Mediterranean Sea.

FIGURE 9

Figure 9. The second EOF mode of 800 m depth (A) and the first EOF modes of 1,200 and 2,000 m of the current experiment (B,C), whose expansion coefficients are well-correlated. Here the variance-preserving spectrum of the 2,000 m first EOF is presented (green line), along with its 95% confidence intervals (dashed gray lines) and the normalized red noise (red line) (D).

The EOF analysis of the reference period also reveals an interannual variability, with a moderate coherence from the upper to deeper layers, showing strong opposite phases in a part of the Eastern Mediterranean, southwest of Crete. Figure 10 presents the two clusters of intermediate (200 and 400 m depth) and deeper depths (800 and 1,200 m). Upper depths (200 and 400 m) are in phase with the mode of 1,200 m depth and out of phase with the first EOF of 800 m. This variability is related to the dense water export from the Aegean Sea, which appears to be related to a gyre-like circulation in the Ionian basin (based on the EOF amplitude distribution around the basin, observed at all depths, but more clearly in the deeper layers), a mechanism reminiscent of the BIOS mechanism.

FIGURE 10

Figure 10. As in Figure 9, for the two clusters of (A) intermediate (first EOF of 200 m and third EOF of 400 m depth) and (B) deeper layers (first EOF of 800 m and second EOF of 1,200 m) with a 2-year cycle of variability as shown from the variance-preserving spectrum of the EOF modes of 400 (A) and 1,200 m (B), respectively.

SBF Experiment

As revealed through the EOF analysis of the results of the SBF experiment, the major variability caused by extremely cold and dry atmospheric conditions over the Mediterranean is found in the first mode and is related to a vertical exchange of buoyancy between the deep intermediate (800 m) and the very deep layers of the basin (>1,200 m). Thus, this variability is related to a variation of the deep stratification of the Mediterranean and is characterized by very long temporal scales (Figure 11). In the correlated depths of 800, 1,200, and 2,000 m, this variability is expressed by the first EOF modes and represents over 83% of the variability of the basin in deep waters (83.9, 89.6, and 89.9% at 800, 1,200, and 2,000 m, respectively).

FIGURE 11

Figure 11. The first EOF modes of 800 (A), 1,200 (B), and 2,000 m (C) depth of the SBF experiment, whose expansion coefficients are well-correlated. Here, the variance-preserving spectrum of the 2,000 m first EOF is presented as (D).

The second mode of variability in the deep layers describes the buoyancy exchange between eastern and western basins (Figure 12). The variance-preserving spectra of the expansion coefficient indicate a strong decennial (~14 years) and interdecennial variability of density. The interdecennial variability may correspond to the water residence time of the Eastern Mediterranean, as calculated earlier. The 14-year peak is near to a 13-year cycle that is present at the buoyancy variability of Sicily Strait and may represent the pulses of feeding the western basin with the dense water of the Eastern Mediterranean Sea. Note that this mode exhibits a strong amplitude at 1,200 m in the Tyrrhenian Sea and minimal signature in the Algero-Provencal basin, but its signal is strong throughout the western basin at the deeper layers, a fact signifying the role of the lateral advection of Levantine waters and the subsequent formation in the Gulf of Lions and dispersion in the bottom layers.

FIGURE 12

Figure 12. The second EOF modes of 1,200 (A) and 2,000 m (B) depth layers of the SBF experiment, whose expansion coefficients are well-correlated. Here, the variance-preserving spectrum of the 2,000 m first EOF is presented as (C).

Stratified Period Experiment

In the period of extreme stratification (caused by excessive warm and humid conditions), the first EOF modes from 200 to 2,000 m depth show that the whole Mediterranean Sea is in phase and oscillates in intercentennial timescales (Figure 13). The variance contribution of the first mode ranges from 46.4% at 200 m to 86.5% at 2,000 m depth. The increased values of the Eastern basin suggest that this variability is driven by this basin, and the intercentennial variability is in accordance with the residence time of the Eastern basin waters that has been estimated to be ~540 years.

FIGURE 13

Figure 13. The first EOF modes from 200 to 2,000 m depth (A–E) of the WBF experiment, whose expansion coefficients are well-correlated. Here, the variance-preserving spectrum of the 2,000 m first EOF is presented as (F).

The second EOF modes of deeper layers (800, 1,200, and 2,000 m depth) of the stratified period experiment exhibits an interdecennial variability (~63 years) (Figure 14). All three depths are in phase with the negative values in the Western Mediterranean and close to zero, which gradually are becoming positive with depth in the Eastern. The two subbasins exhibit a phase difference, suggesting that this mode represents a lateral buoyancy exchange; furthermore, this process appears to be driven by the Western Mediterranean (in accordance with a reversal of the buoyancy exchange between the eastern and western basins during highly stratified conditions similar to the S1a sapropel formation period, as described in section Overturning Functioning of the Mediterranean Sea Under Different Climatic Conditions). Furthermore, a signal of opposite phase to that of the Western Mediterranean is concentrated south of Crete at 800 m, but gradually fills the deeper layers of the eastern basin.

FIGURE 14

Figure 14. The second EOF modes of 800, 1,200, and 2,000 m depth (A–C) of the WBF experiment, whose expansion coefficients are well-correlated. Here, the variance-preserving spectrum of the 2,000 m second EOF is presented as (D).

Discussion

Although the extreme climate simulations are not designed to reproduce the paleocirculations of S1a sapropel deposition and Younger Dryas event (our adjustment does not involve any modification of the wind-stress fields and sea level of these periods). Topper and Meijer (2015) have used numerical simulations to show that the size of the Gibraltar restriction may affect the salinity and temperature of the whole basin but that the size does not have much impact on the inter-basin exchanges and circulation, especially of the Eastern basin. Thus, it is possible to assess the validity of our results based on paleoceanographic analyses of the Younger Dryas and S1a periods, keeping in mind that our simulations cannot be treated as reproductions of the paleocirculation.

Comparison of SBF Results to the Sedimentary Records From Younger Dryas

The overturning circulation of the SBF period in the Mediterranean Sea is presented as more energetic than the current one and as more intense in the eastern basin. The buoyancy exchange in the Gibraltar Strait under these conditions is much more intense than today. In agreement, grain size records from the northern and the southern parts of the Gulf of Cadiz reveal enhanced Mediterranean outflow water velocity (Toucanne et al., 2007) and a strong advection of the Mediterranean outflow water to the Atlantic during the Younger Dryas period (El Frihmat et al., 2014), respectively. Paleoceanographic analysis also suggests an even stronger export of LIW into the western basin during the Younger Dryas event (Dubois-Dauphin et al., 2017).

The residence time of the Eastern Mediterranean deep waters during the strong buoyancy conditions calculated on a semicentennial timescale. This calculation is also reflected in the interdecennial cycles of the inherent oscillations of deep waters, from first and second EOF modes of variability of deeper layers.

Comparison of WBF Results to S1a Sapropel Deposition Period

This study suggests that the Mediterranean Sea exhibited an estuarine behavior during the WBF experiment. Circulation pattern reversal during the sapropel deposition has already been proposed by Thunell and Williams (1989). Although other model studies for the S1 sapropel deposition proposed that the Mediterranean Sea circulation remained anti-estuarine (Myers et al., 1998; Myers, 2002; Meijer and Tuenter, 2007; Adloff et al., 2011), benthic foraminiferal carbon isotope variations indicate a significant presence of North Atlantic Central Water mass in the southern part of the Gulf of Cadiz (El Frihmat et al., 2014) in agreement with our results. Moreover, based on the contourite grain size from the Gulf of Cadiz (core MD99-2341, 582 m depth), a sluggish Mediterranean water outflow is inferred during the Early Holocene, synchronous with sapropel S1 (Toucanne et al., 2007), and records from the Alboran Sea indicate relative low oxygenation (Jiménez-Espejo et al., 2015). Recent findings from neodymium isotopic composition in foraminifera and corals of the Western Mediterranean Sea (Dubois-Dauphin et al., 2017) imply that, during S1a-like conditions, the Western Mediterranean is not very sensitive to hydrological contributions from the eastern basin and that the intermediate water of the western basin is mainly locally formed.

The results of the WBF period simulation suggest also that, in the Sicily Strait, the western basin feeds the eastern with denser waters (Figure 5). This functioning could be an explanation for the total lack of sand grains and layer barren of foraminifera in Malta Graben (Ferraro et al., 2018), which indicates different processes on the seafloor.

The residence time of the Eastern Mediterranean deep waters under WBF was calculated to be semimillennial. This calculation is reflected in the inherent oscillation of deep waters, from the first EOF modes of the variability of the intermediate and deeper layers. The estimate of intercentannial variability (Figure 13) relates well with the dominant variability, 300–600 years, of the enrichment of the redox-sensitive trace metals during the early S1 in the Eastern Mediterranean (Jilbert et al., 2010).

Records from the Adriatic Sea (Tesi et al., 2017) indicated that a weakening of North Adriatic Deep Water production controlled the sapropel development in the Adriatic and Ionian Seas and that the Adriatic S1 onset was synchronous with the other sites of the Eastern Mediterranean, suggesting a wide-basin, physical-driven mechanism, which hampered the dense water formation over the entire Eastern Mediterranean Sea. Notably, during S1 sapropel deposition, the sediment core findings of the Aegean Sea reveal that sapropel was deposited in dysoxic/oxic conditions (Abu-Zied et al., 2008; Triantaphyllou et al., 2009b, 2016; Triantaphyllou, 2014) and that local deep-water formation for the North Aegean Sea during S1 is documented by Schmiedl et al. (2010). This functioning of the Aegean Sea may relate with the interdecennial internal variability to the region and with limited dense water production during the stagnant periods of the Mediterranean Sea (Figure 14), suggesting EMT-like events during the stratified period. Relevant EMT-type events have been also recognized over the last five centuries (Incarbona et al., 2016).

Comparing the Three Experiments

The model results revealed that the Mediterranean Sea exhibits high sensitivity to climatic conditions, its functioning may alternate from concentration to estuarine and that the basin presents inherent interannual, interdecennial, and intercentennial variability. Maybe the most striking result regarding the mean overturning circulation is that, under the reference period, deep water formation in the western basin appears to be stronger and deeper than that during the SBF experiment. This result requires further investigation; however, this paradox has been also pointed out by Cisneros et al. (2019), who, combining oceanographic mooring measurements and sediment cores from the western Mediterranean, argue that the strongest dense-water formation events do not coincide with the coldest winters over the region. Investigation of this result is beyond the scope of this study; however, it is a very intriguing finding.

In all experiments, the upper layers of the basin are dominated by the interannual/decennial variability. The long-term simulations revealed that the integrated buoyancy fluxes in the Straits oscillate on different timescales, but in all the experiments, a decennial variability dominates in the Gibraltar and Sicily Straits.

The second EOF modes of deep waters under both the paleoclimate scenarios showed that the western and eastern basin functions on opposite phases in decennial (for the mixed case) and centennial (for the stratified case) timescales of variability. Especially in the case of the stratified conditions, the phase difference between the western Mediterranean and the region south of Crete suggests the possible existence of a mechanism similar to BIOS (Borzelli et al., 2009; Gačić et al., 2010), which refers to the buoyancy exchange, not between the Adriatic and the Aegean Sea but, between the Western and Eastern basins (where the deep water formation response in the Eastern basin is centered in the Aegean Sea, most probably due to the lack of buoyancy input through the Dardanelles). This observation requires further investigation (especially the processes involved in this result), which however exceeds the context of this study.

More than 90% of the deep layer variability of the Mediterranean basin under extreme climatic conditions can be described by the first two EOF modes. The deeper Mediterranean basin under the current climate conditions is a synthesis of modes of variability, which is more than those during the other two experiments (the two first modes express only ~60%). The above suggests that, during both extreme climatic conditions, the deeper layers tend to exhibit larger lateral scales of circulation features. Having excluded momentum forcing and vorticity exchange between atmosphere and ocean as a potential source of this result, the sources of this behavior could be different for the two extreme cases, ranging from different internal Rossby radii associated to the different stratification profiles to the isolation of the deeper layers during the highly stratified conditions.

A comparison of the spectra of the lateral buoyancy exchanges through the Straits and of the spectra of the density variability through the EOF analysis is rather revealing: the buoyancy flux variability at the Straits is dominated by interdecennial periods, throughout all experiments, whereas the deep-water density EOF temporal variability is dominated by centennial periods, especially for the WBF experiment. This difference could partly be attributed to the fact that the Strait exchanges are affected by the variability at the surface/intermediate layers, which is of a higher frequency than the deeper layers and is partly due to the impact of barotropic motions on the Strait exchanges, which, being faster in response, would tend to shift the spectrum to higher frequencies.

A dominant oscillation mode is identified both at the reference and the WBF experiments and exhibits depth coherent, in-phase inter-decadal oscillations in the deeper layer, which are more energetic in the western basin and diminish in the eastern Mediterranean. This motion is related to the overturning mechanisms of the western Mediterranean, which still operate under the WBF conditions.

Another interesting mechanism of inherent variability is the exchange of buoyancy between the intermediate and deep layers observed at the SBF simulation. This mechanism reflects the alternation of deep water formation (which enhances deep stratification, carrying buoyancy upward) and subsequently a gradual homogenization through turbulent mixing, which gradually homogenizes the water preconditioning for the next formation event.

Finally, EOF modes exhibiting phase differences between the eastern and the western basins reveal another inherent mechanism of a lateral buoyancy exchange through the Sicily Channel. It is a mechanism known to the Mediterranean oceanographic community, which has witnessed at least a part of it through the recognition of the role of inter-basin water exchange on deep water formation processes in the western Mediterranean (e.g., Amitai et al., 2021).

This study has been only a preliminary effort to investigate the inherent dynamics of the Mediterranean Sea under extremely different stratification conditions. The heat and freshwater atmospheric forcing for the two experiments have been adjusted based on two extreme cases within the Holocene period. However, our simulation cannot be regarded as efforts to reproduce the exact circulation during the Younger Dryas and the S1a sapropel deposition periods, as the atmospheric heat and freshwater fluxes have been diagnosed from SST and salinity records, the momentum forcing has remained the same as in the current period (a dynamical inconsistency for the atmospheric circulation), the perpetual year assumption for the forcing violates the very important role of the interannual variability on the Mediterranean deep water formation processes, and the bathymetric changes due to sea-level variability (important especially at the connecting Straits) have not been considered. An effort to reproduce the circulation during these two periods of the Holocene requires all the above issues to be addressed. However, this study has shed light on the characteristic trends of the inherent variability of the Mediterranean and furthermore revealed possible processes similar to the BIOS mechanism acting between the two major Mediterranean basins.

Conclusions

The aim of the present study was to investigate the spatiotemporal characteristics of the inherent variability of the Mediterranean Sea at two extremely different climate states from the Holocene period, as well as the present climate state. To achieve that aim, three long-term simulation experiments were designed: the instrumental period of 1973–1993 was selected based on the simulations by Skliris and Lascaratos (2004); Skliris et al. (2007) as a reference. The numerical experiment simulating the Mediterranean under a strong buoyancy loss to the atmosphere used sea surface conditions similar to the Younger Dryas period, a cold and dry period for the atmosphere. The experiment simulating the weak buoyancy gain from the atmosphere used the sea surface conditions from the period of S1a sapropel deposition, a period characterized by a maximum stratification and warm and humid atmospheric conditions.

The model used in these studies was modified accordingly by removing all climatological relaxation terms, adopting a perpetual year approach (to remove external forcing) and applying monthly diagnosed heat and freshwater fluxes at the surface, as estimated from the sea surface properties given by instrumental climatology for the reference period and sedimental record for the two extreme climate experiments. To simulate as closely as possible the Younger Dryas and S1a periods, solar irradiance was modified accordingly, and the Mediterranean was isolated from the Black Sea.

The mean hydrographic conditions arising from the two extreme case experiments show that, during the SBF experiment, overturning is accelerated, the export of dense Eastern Mediterranean waters to the western basin is intensified, and the residence times become significantly shorter than present-day values. On the other hand, the WBF experiment, through which the Mediterranean gains buoyancy from the atmosphere and switches to a dilution basin, results in an estuarine Eastern Mediterranean, a western basin exhibiting local cells of limited dense-water formation and multi-centennial residence times of the deep waters.

The results of the simulations showed that the inherent variability of the Mediterranean (as revealed through the water density fields) is dominated by inter-decennial oscillations at the present climate, a result supported by a multitude of observational and simulation studies. The Mediterranean imports the buoyancy from the Atlantic, and this lateral exchange exhibits an inherent to the basin interannual variability dominated by inter-decadal periods (with a peak at about 25 years). The inherent variability of the Gibraltar exchange peaks at the same frequencies during both extreme case experiments despite a higher buoyancy import during the SBF experiment and the change of sign during the WBF experiment. The buoyancy exchange in the Sicily and Gibraltar Straits exhibits similar behavior.

Regarding the internal variability of the Mediterranean throughout the SBF and WBF experiments, the major difference from the reference simulation (representing the present climate) is an increased contribution of intercentennial variability to the total variance, especially in the case of the Mediterranean exhibiting an estuarine flow (WBF experiment). In the SBF case, a significant contribution comes from the interannual and inter-decadal periods, which is probably related to the shortening of the residence times of deep waters.

The EOF analysis has revealed the mechanisms of either a lateral buoyancy exchange between the adjacent basins, vertical exchange between the intermediate and deeper layers as well as BIOS-like features for the present climate.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

AS contributed to developing the concept of this work, designing the numerical experiments, performed all the model runs and the production of the figures, and produced a major part of the analysis and the manuscript. VZ and ET contributed to the design of the numerical experiments, the analysis of the results, and editing of the manuscript. AG and MT contributed to the design of the experiments, selecting the periods of the paleoceanographic record to be used for the diagnosis of the atmospheric forcing, and also for the manuscript. IM has contributed to the modifications of the POM model version to run under the diagnosed fluxes and its implementation. NS has provided the model code and forcing fields for the simulations of the current conditions and his expertise and guidance. All authors contributed to the article and approved the submitted version.

Funding

This work has been funded by the project MedEcos (Decadal-Scale Variability of the Mediterranean Ecosystem) and jointly funded by the EU F.P.6 ERA-network MarinERA and by the Greek General Secretariat of Research and Technology, Ministry of Development (project number 2010SE01380001).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors are grateful to Vincenzo Artale for his valuable advice and suggestions. The comments of the two reviewers contributed greatly to improving the work. AS acknowledges support through a PhD fellowship (HCMR), by the Action National Network on Climate Change and its Impacts – CLIMPACT, funded by the Public Investment Program of Greece (GSRT, Ministry of Development and Investments).

References

Abu-Zied, R. H., Rohling, E. J., Jorissen, F. J., Fontanier, Ch., Casford, J. S. L., and Cooke, S. (2008). Benthic foraminiferal response to changes in bottom-water oxygenation and organic carbon flux in the eastern Mediterranean during LGM to Recent times. Mar. Micropaleontol. 67, 46–68. doi: 10.1016/j.marmicro.2007.08.006