Hotter and Weaker Mediterranean Outflow as a Response to Basin-Wide Alterations

Time series collected from 2004 to 2020 at an oceanographic station located at the westernmost sill of the Strait of Gibraltar to monitor the Mediterranean outflow into the North Atlantic have been used to give some insights on changes that have been taking place in the Mediterranean basin. Velocity data indicate that the exchange through the Strait is submaximal (that is, greater values of the exchanged flows are possible) with a mean value of −0.847 ± 0.129 Sv and a slight trend to decrease in magnitude (+0.017 ± 0.003 Sv decade−1). Submaximal exchange promotes footprints in the Mediterranean outflow with little or no-time delay with regards to changes occurring in the basin. An astonishing warming trend of 0.339 ± 0.008°C decade−1 in the deepest layer of the outflow from 2013 onwards stands out among these changes, a trend that is an order of magnitude greater than any other reported so far in the water masses of the Mediterranean Sea. Biogeochemical (pH) data display a negative trend indicating a gradual acidification of the outflow in the monitoring station. Data analysis suggests that these trends are compatible with a progressively larger participation of Levantine Intermediate Water (slightly warmer and characterized by a pH lower than that of Western Mediterranean Deep Water) in the outflow. Such interpretation is supported by climatic data analysis that indicate diminished buoyancy fluxes to the atmosphere during the seven last years of the analyzed series, which in turn would have reduced the rate of formation of Western Mediterranean Deep Water. The flow through the Strait has echoed this fact in a situation of submaximal exchange and, ultimately, reflects it in the shocking temperature trend recorded at the monitoring station.

Time series collected from 2004 to 2020 at an oceanographic station located at the westernmost sill of the Strait of Gibraltar to monitor the Mediterranean outflow into the North Atlantic have been used to give some insights on changes that have been taking place in the Mediterranean basin. Velocity data indicate that the exchange through the Strait is submaximal (that is, greater values of the exchanged flows are possible) with a mean value of −0.847 ± 0.129 Sv and a slight trend to decrease in magnitude (+0.017 ± 0.003 Sv decade −1 ). Submaximal exchange promotes footprints in the Mediterranean outflow with little or no-time delay with regards to changes occurring in the basin. An astonishing warming trend of 0.339 ± 0.008 • C decade −1 in the deepest layer of the outflow from 2013 onwards stands out among these changes, a trend that is an order of magnitude greater than any other reported so far in the water masses of the Mediterranean Sea. Biogeochemical (pH) data display a negative trend indicating a gradual acidification of the outflow in the monitoring station. Data analysis suggests that these trends are compatible with a progressively larger participation of Levantine Intermediate Water (slightly warmer and characterized by a pH lower than that of Western Mediterranean Deep Water) in the outflow. Such interpretation is supported by climatic data analysis that indicate diminished buoyancy fluxes to the atmosphere during the seven last years of the analyzed series, which in turn would have reduced the rate of formation of Western Mediterranean Deep Water. The flow through the Strait has echoed this fact in a situation of submaximal exchange and, ultimately, reflects it in the shocking temperature trend recorded at the monitoring station.

INTRODUCTION
The Mediterranean Sea (MedSea hereinafter) is the only basin away from polar regions where open-ocean deep convection reaching the ocean bottom happens, despite its location in temperate latitudes. The convection is the result of the net buoyancy flux to the atmosphere, mainly determined by the fresh water deficit (evaporation minus precipitation and river runoff, E − P) of the basin. It drives an open thermohaline circulation cell starting and ending in the Strait of Gibraltar (SoG, hereinafter), the only relevant connection of the MedSea with the open ocean. The cell starts with a surface inflow of Atlantic water that travels the whole Mediterranean Sea to be transformed into the very salty Levantine Intermediate Water (LIW) in the easternmost basin of the MedSea (Figure 1). There, it sinks and initiates the way back to the SoG and to the open ocean as a dense underflow, thus completing the thermohaline cell (see Hopkins (1985) or Tsimplis et al. (2006), for instance, for a review on this general topic).
In addition to intermediate waters, deep waters are formed in different places of the MedSea, mainly in the Adriatic and Aegean Seas in the eastern basin (Eastern Mediterranean Deep Water, EMDW; Schlitzer et al., 1991;Roether et al., 1996), and in the Gulf of Lions in the western basin (Western Mediterranean Deep Water, WMDW; Stommel, 1972). The basic thermohaline cell becomes more complex since both basins have their own closed thermohaline circulations (Tsimplis et al., 2006) and bottom water renewal processes. The EMDW is mainly renewed by upwelling inside the eastern basin itself Roether and Schlitzer, 1991) and only a very small fraction seems to be drained out to the western basin through the Sicily Strait (Astraldi et al., 1999). On the contrary, a significant fraction of WMDW finds its way out to the Atlantic Ocean through the SoG (Stommel et al., 1973;Kinder and Parrilla, 1987), even though the depth of the SoG's main sill (Camarinal sill, 290m, Figure 1) is less than Sicily's west and east sills (360m and 430m, respectively, Astraldi et al., 1999). Bernoulli aspiration of WMDW in the Alboran Sea by the energetic lower layer flow at the SoG (Bryden and Stommel, 1982;Naranjo et al., 2012) is responsible for the differentiated pattern of both basins, as the Strait of Sicily does not hold any comparable bottom current.
Levantine Intermediate Water and WMDW make up the main bulk of the Mediterranean outflow in the SoG (Q 2 hereinafter, assumed negative as it flows westwards), to which small fractions of other intermediate and deep waters may eventually contribute . The Mediterranean outflow has to flow out beneath a slightly greater Atlantic inflow (Q 1 , positive), moving over a very irregular seafloor with steep sills and marked contractions (Figure 1). This strongly constraining topography leads to hydraulic control of the exchanged flows, an idea first suggested by Bryden and Stommel (1984) and further addressed by Armi and Farmer (1985). A relevant result of this theory is that the exchanged flows attain a maximum value that cannot be exceeded whenever the exchange is hydraulically controlled (maximal exchange, Armi andFarmer, 1985, 1987). In other words, |Q 2 | has an upper limit. From this point of view, the twolayer hydraulic theory provides a useful frame to interpret the observations collected at Espartel sill (Figure 1), the westernmost sill of the SoG, which constitute the dataset analyzed and exploited in the present study. Next Section revises some relevant aspects of this theory applied to the SoG.

The Hydraulically-Controlled Exchange Model
The objective of this work is to infer changes in the MedSea from single-point observations at the site where it joins with the global ocean (Figure 1). Within the frame of two-layer theory, the usual approach is to consider the MedSea as a homogeneous basin connected to the Atlantic Ocean, also homogeneous, through a channel (the SoG). Actually, the outflow consists of different water masses of very similar density, which justifies homogeneity for dynamic purposes as a first order approach. The approach, however, cannot be applied to the water masses analysis carried out later on.
The simplest model for the SoG is a non-rotating channel of rectangular section with no mixing between layers (Armi, 1986;Farmer and Armi, 1986). Critical section is a key concept: it is a section where the composite Froude number, G, defined as is G = 1. Then, the section exerts hydraulic control on the flow. In equation [1], F i is the internal Froude number of layer i, g = ρ ρ 0 g is the reduced gravity with ρ = ρ 2 − ρ 1 the density difference between layers (also between basins), and ρ 0 a reference density. The meaning of the remaining variables is shown in Figure 2.
If the flows are hydraulically controlled at two critical sections connected by a region of sub-critical flow, then the exchange is maximal (Armi and Farmer, 1987). The condition G = 1 is usually achieved in notable topographic sections, generally the narrowest and the shallowest. For steady exchange, the SoG fulfills these requirements at Tarifa Narrows and at the main sill of Camarinal (narrowest and shallowest sections, respectively, Figure 1) and, therefore, it is candidate to sustain maximal exchange (Farmer and Armi, 1988).
The constrain imposed by the maximal exchange together with the climatology over the basin, determine the global characteristics of the Mediterranean and link them to the flow properties at the SoG, a result of obvious interest to our study. The issue was first addressed by Bryden and Kinder (1991), who assumed a channel of triangular cross-section and a density difference brought about by salinity exclusively ( ρ = β (S 2 − S 1 ) , β the haline contraction coefficient of seawater). The salinity of the MedSea, S 2 , was left as an output of the model. For a net evaporation in the basin of 0.6 m/year (equivalent to 0.04Sv), they obtained steady-state inflow Q 1 = 0.92Sv and outflow Q 2 = 0.88Sv, and S = S 2 − S 1 ≈2, close to the actual values.
At short time-scales, however, the exchange is markedly timedependent. Tides provide a net barotropic flow strong enough to override the hydraulic control at Camarinal sill, particularly during the intensified flows of spring tides (Farmer and Armi, 1988;Sannino et al., 2004;Sánchez-Garrido et al., 2011). The loss of the control leads to the release and further radiation into the MedSea of packets of large amplitude internal waves, a well-known feature of this region (Farmer and Armi, 1988;Brandt et al., 1996;Vlasenko et al., 2009;Sánchez-Garrido et al., 2011). When this happens, the hydraulic control of Camarinal is transferred to Espartel sill, where it is in force most of the time Sannino et al., 2009). In addition, the loss of control boosts the contribution of tidally induced eddy-fluxes (positive correlations of vertical interface  Sammartino et al., 2015).
fluctuations and tidal currents) to the exchange (Bryden et al., 1994;Sannino et al., 2004;Vargas et al., 2006) at either critical section. Therefore, large eddy-fluxes at a control section lead to hydraulic control intermittency. Alternately, very weak eddyfluxes (or no eddy-fluxes at all) would suggest lastingness of the control. The frequent loss of hydraulic control at Camarinal sill (Farmer and Armi, 1988;García Lafuente et al., 2000Vargas et al., 2006;Sánchez-Garrido et al., 2008 and the large associated eddy-fluxes there (up to 0.5Sv or 40% to 60% of the total exchange; Bryden et al., 1994;Sannino et al., 2004;Vargas et al., 2006) contrast with those at Espartel (∼0.04Sv or 5%; Sánchez-Román et al., 2009;Sammartino et al., 2015). This particular behavior stresses the potential of Espartel sill replacing Camarinal sill as a long-lasting hydraulic control section.
Should the critical value G = 1 be achieved in Espartel section, it would basically be done through F 2 (equation [1]), as the large width and thickness of the surface layer (Figure 1) imply very low velocity u 1 . Similar reasoning applies to F 1 in Tarifa Narrows due to the large thickness of layer 2. These approximations, already suggested by Farmer and Armi (1986), may be interpreted as if Tarifa Narrows controls the inflow Q 1 , which is the active layer there, while Espartel sill does the same with the outflow Q 2 . It underlines the experimental benefits of monitoring the outflow in this section in order to discriminate the state of the exchange: F 2 = 1 for possible maximal exchange (F 1 in Tarifa Narrows must be inspected), and F 2 < 1 for submaximal exchange. It also emphasizes the importance of the collected data in relation to the global properties of the MedSea, particularly in case of submaximal exchange, as it is outlined next.

Buoyancy Fluxes, Water Formation and Outflow
In our simple MedSea-SoG system the other key variable is the amount of Mediterranean water (Q M ) formed due to buoyancy fluxes (B 0 ) to the atmosphere, which should be equal to |Q 2 | in an ideal steady-state on yearly basis. In maximal exchange (G = 1), |Q 2 | has attained its maximum. If, eventually, Q M > |Q 2 | during FIGURE 2 | Sketch of the two-layer exchange in the Strait of Gibraltar. Subscripts 1 and 2 stand for variables in the upper and lower layers, whose densities are ρ 1 and ρ 2 , respectively. In a steady state, Q i 's are constant but u i and h i are not. The three notable topographic sections, Espartel and Camarinal sills and Tarifa narrows, have been indicated in both, the sketch and the inset map. Camarinal sill and Tarifa narrows are the candidate critical sections for steady exchange. In case of time-dependent flow (tides) Espartel sill and Tarifa narrows are much more likely candidates for critical sections (see text). a given period, the excess of formed water cannot flow out unless either the interface goes up or g increases (equation [1]), which requires larger density difference ρ. The interface cannot rise because it would block partially the inflow, preventing the water balance of the basin Q 1 + Q 2 = (E − P) The excess of water has to remain in the basin contributing to increase its density (hence, ρ) in successive years until reaching a new equilibrium. This is how the hydraulics of the SoG determines the MedSea properties. If, on the contrary, Q M < |Q 2 | the equilibrium is easily maintained as |Q 2 | has no constrain for diminishing. But the exchange would change to submaximal. If the exchange was already submaximal, there is no impediment for |Q 2 | to increase if Q M > |Q 2 | or to decrease if Q M < |Q 2 |, making G (F 2 in the previous approximation) increase or decrease accordingly. As pointed out by Garrett et al. (1990), it is only if the exchange is submaximal that changes occurring in the MedSea can bring about a non-lagged footprint in the outflow at the SoG (i.e., at Espartel section). If it were maximal, such changes would be only reflected in Q 2 after several years (Garrett et al., 1990).
The formed water Q M is proportional to the buoyancy flux B 0 , given by Gill (1982) where g is gravity acceleration, α and β the thermal expansion and haline contraction coefficients of seawater, c p its heat capacity, and ρ 0 and S a reference density and surface salinity, respectively (g = 9.81ms −2 , α = 2.6·10 −4 K −1 , β = 7.44·10 −4 , c p = 3984JKg −1 K −1 , ρ 0 = 1028kgm −3 and S = 37 in this study). Q net is the net heat flux into the sea (Wm −2 , positive downwards), and (E − P) is the net evaporation which includes river run-off (ms −1 , positive upwards). They give the buoyancy flux B 0 in kgm −1 s −3 (Nm −2 s −1 ). Negative heat flux (cooling) and positive net evaporation cause positive buoyancy flux to the atmosphere.
A crude estimation of the amount of water of density ρ 2 that could be formed from water of density ρ 1 under a constant and homogeneous buoyancy flux B 0 acting on a basin of area A MED (say, the MedSea) would be (for Q net ≈−3Wm −2 and (E − P)≈0.7 m/year (=2.22·10 −8 ms −1 ), representative values of the Mediterranean basin, B 0 ≈8·10 −6 Nm −2 s −1 from equation [2], which gives the realistic value of Q M ≈1.03Sv in equation [3] with A MED = 2.5 · 10 12 m 2 and ρ≈2kgm −3 ). Obviously, the actual processes of deep water formation are much more complex (pioneer paper by MEDOC Group (1970), or Houpert et al. (2016 for a recent comprehensive study). The simplicity of equation [3], however, provides a linear relationship for the expected dependence of Q M on B 0 . The simple two-layer model summarized above links Q M and Q 2 , thus extending this relationship to the outflow and the buoyancy flux, and is a well-suited first approximation for the discussion of the results presented in this study. The model, however, is a steady-state model that could be possibly extended to very slowly (longterm) varying conditions such as trends, but not to short-term scales. Therefore, the present study focuses on these long-term changes or trends.

Currentmeter and Temperature/Salinity Data
The monitoring station at Espartel sill is located at 35 • 51.71 N and 5 • 58.22 W (Figure 1) at ∼360 m depth. The instrumented mooring line is equipped with an up-looking Acoustic Doppler Current Profiler (ADCP RDI Workhorse Long Ranger 75 kHz), embedded in a 1.5 m wide buoy deployed at ∼17 m above the seafloor, a Conductivity-Temperature sensor (CT, Seabird SBE37-SMP) and a single point current meter (Nortek Aquadopp DW), both clamped approximately 3 m below. Sampling interval for all instruments has been set to 30 min. The station was first deployed in October 2004, and has been working since, with a large gap in 2011-12 due to a line loss, and few shorter gaps due to different technical issues. The dataset used in this study covers the period from October 2004 to February 2020.
The accurate calibration of CT probes is critical because the instrumental drift of the sensors (O(10 −3 ) • Cyr −1 for temperature; O(10 −3 )yr −1 for salinity) is comparable to the trends of the Mediterranean water masses characteristics reported in the literature (see Vargas-Yañez et al., 2017, for instance). CT data in this study come from a bunch of 6 different instruments, all which have been regularly calibrated throughout the probe life every one or two years. Details about ADCP datasets are presented in Supplementary Appendix A, which includes a short history of the station and other relevant issues.

Biogeochemical Data
After the gap, in year 2012, autonomous sensors for continuous recording of pH (at total scale and reference temperature of 25 • C) and CO 2 partial pressure (pCO 2 , µatm) (SAMI-pH and SAMI-pCO 2 , Sunburst Sensors, LLC.) were incorporated to the monitoring line and placed at the same depth as the CT probe. Unfortunately, the instruments were lost in year 2017 after an accident underwent by the mooring line and they have not been replaced until very recently. It reduces the pH data availability from August 2012 to March 2017, and to a shorter period in the case of pCO 2 data. The sampling interval was initially set to 60 min, but a battery run off happened in summer 2013 (which caused a six-month data gap) advised changing the interval to 120min to extend the battery life.
According to manufacturer, precision and accuracy of measurements were ∼0.001 and ± 0.003 units for pH and ∼ 1 µatm and ± 3 µatm for pCO 2 , approximately. Repeated measurements conducted in the laboratory with the probes submersed in aged seawater and high CO 2 seawater to validate these specifications provided values of ∼0.003 and ± 0.005 units for pH and ∼3 µatm and ± 5 µatm for pCO 2 . SAMI-devices were nonetheless, regularly checked upon periodic servicing of the line on board and pH data obtained by the spectrophotometric technique (Clayton and Byrne, 1993) in samples collected in oceanographic cruises carried out in the area were used to validate the series and identify possible drifts and drops of the autonomous sensors (see section "Other Biogeochemical Properties"). Thus, samples at the mooring site were taken within the Mediterranean outflow in August 2012 (4), May 2013 (4), November 2014 (2), June 2015 (5), and March 2017 (9) under the frame of the Gibraltar Fixed Time Series monitoring program (GIFT, Flecha et al., 2019). Samples were taken directly from Niskin bottles in 10 cm path-length optical glass cells and submitted to a Shimadzu UV-2401PC spectrophotometer containing a 25 • C-thermostated cells holder after addition of m-cresol purple as indicator. Precision and accuracy of these discrete and independent pH T25 measurements were determined from measurements of certified reference material (CRM batches #97 and #136 provided by Prof. Andrew Dickson, Scripps Institution of Oceanography, La Jolla, CA, United States) and were equivalent to ± 0.0048 and ± 0.0047 units for pH T25 .

Other Data
Data of mass and heat fluxes to the atmosphere from the fifth generation of ECMWF atmospheric reanalysis of the global climate (ERA5, Hersbach et al., 2020), distributed by the Copernicus Climate Change Service (C3S), have been retrieved in order to estimate buoyancy fluxes in the Mediterranean basin.

WATER PROPERTIES OF THE OUTFLOW AT THE MONITORING STATION
Temperature, Salinity, Density distance implies that the instrument is sampling the densest (hence, coldest) water that leaves the MedSea. Temperature fluctuates noticeably at tidal timescale due to mixing (lightgray dots in Figure 3). It is not local mixing (the instrument is far enough from the interface for local mixing to reach its depth) but advection of waters mixed upstream in the Tangier basin (Figure 1), where intense mixing happens (Wesson and Gregg, 1994). Conventional numerical filtering is not advisable to remove these fluctuations because the low-passed series smooths out the footprint of the coldest samples, which are particularly relevant since they inform about the deepest waters that are leaving the MedSea. A possibility of removing tidal fluctuations without losing this information is to extract the coldest sample in each semi-diurnal tidal cycle and decimate the series to a sample per cycle (García Lafuente et al., 2007;Naranjo et al., 2012Naranjo et al., , 2017. Light-blue dots in Figure 3 show the resulting series where most semidiurnal and diurnal tidal variability has been removed. Dashed-blue line in Figure 3 shows a clear trend of 0.144 • C/decade in θ (see also Table 1). This value, however, is misleading as the trend during the first half of the record is much less (0.029 • C/decade, dashed-green line), whereas it is more than double in the second half (0.339 • C/decade, dashed-red line). It seems to be around year 2013 that the trend changed. And it did in an order of magnitude with regards to the before-2012 value! Although this pattern was anticipated in Naranjo et al. (2017), they could not give definitive ratification as they used a correspondingly shorter length of this time series (up to year 2015). The extended series available nowadays does confirm the existence of that trend, which is so large and unexpected that has motivated the present study. Similar analysis has also been carried out for salinity, S. Only salinities corresponding to the samples of minimum θ have been selected in order to use the same subset of data for both variables. The analysis yields a positive trend of 0.023 decade −1 for the whole series (Table 1), which increases to 0.040 decade −1 in the second part of the series after the year 2012 gap. Before 2012, the trend does not differ from 0 within the 95% confidence level ( Table 1). Because of the opposite effects of temperature and salinity, the density trend depends on the relative weight of θ and FIGURE 4 | θ − S diagram of the data recorded by the station. Top-left diagram shows the data sampled every 30 min during the first (blue dots) and last (red dots) three years, respectively (see text for definition of "year"). Black rectangle indicates the portion of the diagram enlarged in the right panel. Dots in this diagram correspond to the samples of minimum θ extracted from the whole series. They have been colored by years, according to the code in the legend. Large dots are the mean value of a given year, labelled next to them, and follow the same color code. The percentage of data available each year is indicated in the bottom-left plot: 100% correspond to 17520 30 min-data or 705 minimum-θ-data per year.
S trends. In this case, the former prevails and the trend of density anomaly σ θ is negative: −0.009 kgm −3 /decade for the whole series (Table 1). It changes markedly to −0.034 kgm −3 /decade in the second part of the series as a consequence of the sharp increase of θ during the same period.
These results are illustrated in the θ-S diagrams of Figure 4. The available data have been grouped together by years. As the monitoring station was first deployed in October 2004, years have been defined from October 1 to September 30 of the following year: year 2005 goes from October 1, 2004 to September 30, 2005 and so on. Light-blue dots in the top-left panel show the data sampled at t = 30min during the first three years (2005, 2006, 2007, see legend), and red dots do the same for the last three years (2018,2019,2020). The increase of potential temperature of the water flowing past the station is palpable. The rise of salinity and the diminution of density are less obvious but still detectable.
Right panel shows the samples of minimum θ and their yearly mean (large dots). The lack of data in some years can bias the mean, as it could be the case of year 2012 with a low percentage of data (bottom-left plot), although, curiously, it remains close to the neighbor years. The distribution of large dots clearly suggests two different periods, a pattern that would still stand even if some means are biased: a first one from 2005 to 2014 when years group close together around θ = 13.15 • C, S = 38.4, and a second period of continuous and marked rise of mean θ and also of S to a lesser extent, with the exception of year 2016. The partition agrees with the previous one found in the analysis of trends.
Other Biogeochemical Properties Figure 5A shows the pH time series in total scale at a reference temperature of T = 25 • C (pH T25 ) measured at the monitoring station. No data during the first part of the series are available, as the biogeochemical sensors were only deployed in summer 2012. Therefore, they illustrate the temporal evolution of this variable in the Mediterranean outflow until 2017.
These series have been contrasted against the discrete spectrophotometric pH values deduced from the seawater samples collected in the outflow layer during the five cruises that overlapped the monitoring period. The resulting pH T25 values, 1 | Summary of estimated trends (with 95% confidence level), mean and std (in italics) and number of data used in each calculation (in bold in brackets) for different variables of the study: Q 2 is the total outflow, Q 2s−ν is the slowly part of it and Q 2E−F is the tidally driven contribution of the eddy-fluxes; h is the interface height, θ is the potential temperature, S the salinity, σ θ the density anomaly, F 2 the lower layer Froude number, pH T25 the seawater pH and pCO 2 the CO 2 partial pressure.

Variable
Trend ± 95% confidence level mean ± std (Number of data) 2004-2020 2004-2012 2013-2020 All trends are for decade −1 (Sv for the flows, m for h, • C for θ, kgm −3 for σ θ , and µatm for pCO 2 ; S, F 2 and pH T25 are dimensionless). Shaded cells indicate trends no different from 0 at the 95% confidence level. Three values are given for each variable, corresponding to the whole series, and the first and second parts, respectively (as indicated in the Table heading), except for pH T25 and pCO 2 , whose data only covers a portion of the second part. which come from independent samples, are shown in Figure 5A in red symbols, and agree with the overall trend of the series quite satisfactorily.
During the four-and-half-year period, pH T25 data ranged between a maximum of 7.9193 measured in October 2013 to a minimum of 7.8520 in January 2017 and exhibits the expected trend to diminish in the present scenario of ocean acidification, a trend also confirmed by the values obtained through spectrophotometry (red diamonds in Figure 5A). The decrease is not steady but step-like, the most obvious jump occurring in early 2016. Overall, the series shows a remarkable negative trend of −0.0462 ± 0.0006 decade −1 , significant at the 95% confidence level. Figure 5B shows complementary in situ pCO 2 data from the SAMI sensor, which exhibits a trend inversely and significantly correlated to that of pH T25 (not shown). During the period of available data, the pCO 2 in the outflow at the depth of the monitoring station showed a remarkable rise of 186.8 ± 1.3 µatm decade −1 . Nevertheless, and in concordance with the pH T25 series, pCO 2 evolution did not exhibit a smooth increase but a more evident rise from early 2016.

The Interface Between Atlantic and Mediterranean Waters
The first step to estimate the Mediterranean outflow is the computation of the height of the interface between Atlantic and Mediterranean waters, which implicitly requires a definition of interface. Much has been written about this issue and different solutions have been proposed (Bryden et al., 1994;García Lafuente et al., 2000Tsimplis and Bryden, 2000;Sánchez-Román et al., 2009Naranjo et al., 2014;Sammartino et al., 2015). The depth of null velocity, the obvious definition in a two-way steady exchange, is not applicable due to tidally driven reversal of the flows, which makes both layers to flow in the same direction during a portion of the tidal cycle and prevents zero-velocity from being achieved. Selecting a material surface is another possibility. Salinity is the best candidate and specific isohalines have been proposed and used to this aim (Bryden et al., 1994;García Lafuente et al., 2000, 2013. It requires time series of S profiles throughout the water column. If they are not available, which is the case in this study, the procedure is not affordable. A third possibility is to link the interface position with the depth of maximum vertical shear of horizontal velocity (Tsimplis and Bryden, 2000;Sánchez-Román et al., 2009Naranjo et al., 2014;Sammartino et al., 2015). The technique is well suited to the available series of ADCP velocity profiles and has been adopted here. The detailed description of the procedure can be seen in Sammartino et al. (2015) and has been followed to the letter in this work. Figure 6B shows the time evolution of the interface height h 2 at the monitoring station, whose mean value is 172.8 m from the seafloor ( Table 1). A marked subinertial variability stands out (the substantially larger tidal oscillations have been filtered out) overlapping an unclear, but still recognizable, seasonal cycle with the interface higher in spring and deeper in late summer-autumn. A comprehensive description of these issues can be seen in Sammartino et al. (2015) and are not addressed here. For the purpose of this work, the result of interest is the trend of the interface to become shallower at a rate of 1.96 m/decade. Following the splitting of the whole series in two parts suggested by θ (Figure 3), the same differentiation has been accomplished with this series. Table 1 indicates that most of the interface elevation took place in the first part before the data gap of year 2012, whereas the trend is not significantly different from 0 at the 95% confidence level during the second one.

The Estimated Outflow
The outflow Q 2 has been also calculated following the procedure in Sammartino et al. (2015) after revising the lower part of the vertical profiles of horizontal velocity. The revision has been motivated by the availability of near-bottom velocity data collected by a new instrument incorporated into the monitoring station since September 2016. The data offer an ameliorated description of the bottom boundary layer and compel to modify the vertical profile of velocities within a few tens of meters above the seafloor with regards to the profile used in Sammartino et al. (2015). A detailed explanation of how it has been recalculated in light of these new data is presented in Supplementary Appendix A.
The low-passed (subinertial) outflow Q 2 computed with the new profiles ( Figure 6A, blue line) differ only slightly from the one computed using Sammartino et al. (2015)'s velocity profiles. Its mean value (2004 to 2020) is -0.85 ± 0.13 Sv (standard deviation), whereas it is −0.82 ± 0.16 Sv if Sammartino et al. (2015)'s profiles were used instead. The ameliorated estimation gives a 3.7% greater and, also, less fluctuating outflow. Of interest is the positive trend identified in the whole series of Q 2 (Figure 6A, dashed-blue line, Table 1), estimated in +0.016 Sv/decade ( Table 1). Since Q 2 is negative, this trend diminishes the outflow. Care must be taken, however, because the trend was negative before the gap of year 2012 (Table 1), meaning that Q 2 was increasing during the first part, to reverse the tendency in the second one.
Red line in Figure 6A displays the slowly varying part of the outflow (Q 2s−ν ), which results from filtering out tidal fluctuations from the horizontal velocity and interface motions, and computing the transport with the low-passed series. The result differs from the total outflow (blue line) because of the already mentioned eddy-fluxes (Q 2E−F ; see Vargas et al., 2006, for details about Q 2s−ν and Q 2E−F computations), which hardly contribute by ∼5% to the outflow (sienna line in Figure 6A, see Table 1) in agreement with Sánchez-Román et al. (2012) and Sammartino et al. (2015). Both terms exhibit positive trends for the whole series (Table 1), the trend of Q 2s−ν being comparable to the one of Q 2 , whereas the trend of Q 2E−F is much smaller. During the second part the sign of both trends is different (Table 1), a fact that will be discussed later.

The Froude Number at Espartel Section
The Froude number at Espartel has been computed in order to assess the hydraulic state of the exchange. The lack of information in the upper layer only allows for the estimation of F 2 in equation [1]. Notwithstanding and according to the arguments given in Section "The Hydraulically-Controlled Exchange Model", F 2 ≈ G is a reasonable approximation and, therefore, F 2 is a good indicator of that state. It can be evaluated in terms of Q 2 (instead of u 2 ) and the interface height h 2 as (see Supplementary Appendix B and Supplementary Figure B1 for details) for a realistic trapezoidal cross section of bases W b (=2100m) at the sea bottom and W 0 (=7000m) at the mean height of the interface H 0 (=173m). Dimensionless constant c is defined as c = (W 0 − W B ) H 0 = 28.32. Figure 7 displays the time series of F 2 computed from the series of Q 2 and h 2 presented in Figure 6.
Except for an isolated event by the end of year 2009, F 2 does not reach the critical value of 1. Including F 1 in the estimation of G in equation [1] does not change the result because F 1 hardly reaches 0.1 due to the low vertically averaged velocity u 1 in this section. Despite being close, Espartel does not meet the condition of critical section and, therefore, the exchange is not maximal, but submaximal. Greater exchanged flows are possible. Figure 7 indicates a decreasing trend of F 2 , which could have been anticipated from the tendencies of Q 2 to diminish in magnitude and h 2 to increase (Figure 6), as both them act in the same direction to diminish F 2 (equation [4]). The simultaneity of such trends in Q 2 and h 2 are counterintuitive: the higher the interface, the larger the cross-area for the Mediterranean water to flow out. The expected result would be an increased outflow and not a diminution. The explanation of the paradox is in Table 1. The rise of the interface took place during the first part of the series previous to the data-gap, when Q 2 had negative trend that increased the magnitude of the outflow. It is the intuitive physical behavior, which leads to a positive trend of F 2 in this first part (Table 1). Thus, F 2 increased driven by the augmented outflow, which approached maximal exchange conditions that, in the end, did not attain. The interface had no significant trend during FIGURE 7 | Time evolution of the lower layer Froude number (F 2 ) computed according to equation [4]. Dashed blue line is the linear trend for the whole series and dashed green and red lines show the trends before and after the gap of year 2012, respectively. All numerical values, including 95% confidence levels, are indicated inside the text-boxes with the same color code. Mean and standard deviation for the whole series, and for the first and second parts are also indicated. the second part whereas the one of Q 2 changed sign ( Table 1), reducing progressively the outflow and, therefore, F 2 .  as 2008 or 2011, when no WMDW was reported to have been formed (Houpert et al., 2016).

Buoyancy Fluxes and Water Formation
In the Gulf of Lions, buoyancy flux of three out of the seven last years, namely 2014, 2016 and 2020, had the largest negative anomaly of the whole series ( Figure 8B). Except for year 2020, the same applies to the whole basin (green bars). They showed minima well below the mean, which makes the second half of the period be negatively anomalous in both, the Gulf of Lions and the whole basin. These results agree with Margirier et al. (2020), who indicate the absence of intense deep convection in the Gulf of Lions between 2014 and 2017 (their study ends in year 2018). Interestingly, this period coincides with the second part of the oceanographic observations, when the strong trend of θ was detected in the monitoring station. The coincidence encourages a similar partition of the B 0 series and estimating the trends before and after year 2013 separately. Table 2 shows the results of such analysis. In addition to the three-winter-month averages presented in Figure 8 (DJF, Winter3 block in Table 2), two other averages have been considered: five-winter-month (NDJFM, Winter5 block) and all year round (Annual block). Trends are given in Nm −2 s −1 decade −1 and the cells in Table 2 have been tagged according to the sign of the estimated trend. They are nosignificant at the 95% confidence level, which is not surprising considering the reduced number of data-points in the fits. There is, however, a remarkable pattern in the Table: all trends during the period 2013-2020 are negative except for the positive one found in the Levantine basin in block Winter3 (although it is the less positive of the three trends in Winter3 block). Moreover, most trends in the period 2005-2012 are positive (Gulf of Lions and MedSea in block Winter5 are the exceptions). All in all, the conclusion that during the first part of the series the trend of the buoyancy was positive and changed to negative during the second part is reasonably supported. Equation [3] would indicate the same signs in the trends of Q M . Within the framework of the two-layer SoG-MedSea model with balance between water production and export, during the first part Q 2 should increase (negative trend, as it is a negative quantity) and decrease during the second part (positive trend), which is what Table 1 reflects.

Fraction Composition of the Outflow
The time evolution of the properties of the water flowing pass the monitoring station may be the result of changes in the proportion of the Mediterranean waters in the outflow, or a concomitant evolution of the properties of each of the Mediterranean waters that form the outflow, or both. We focus on the WMDW and LIW, which are the main bulk of the outflow at the monitoring station, the rest being a very small fraction of Atlantic water, mainly North Atlantic Central Water (NACW) entrained by the Mediterranean flow a short distance upstream in the Tangier Basin (García Lafuente et al., 2007. An estimation of the fraction of the different water masses in the outflow has been carried out following García Lafuente et al. (2007). The triplet of points used by these authors was [12.8 • C, 38.45], [13.22 • C, 38.56] and [15 • C, 36.2], representing WMDW, LIW and NACW, respectively. The water masses characteristics change from place to place and with time, and there are no unique values to represent them. A requirement for the fraction analysis to be feasible is that the θ-S points of the water samples to be decomposed lay inside the triangle defined by the triplet of reference to avoid fractions greater than 1 or less than 0, which have no meaning. The selection of the triplet is a tradeoff between fulfilling this requirement and maintaining values as representative of the implied water masses as possible. The points used in García Lafuente et al. (2007) do not fulfill the first condition with the present series (left panel of Figure 9) because data of the available series at that time were grouped closer to WMDW-type (cooler, see Figure 4) (Figure 9).
The first hypothesis that the time evolution of water characteristics stems from changes in the fractions of each contributor implicitly assumes that the θ-S characteristics of the triplet remain unchanged. Right panel of Figure 9 shows a sudden increase of LIW-type fraction from year 2015 onwards, along with a corresponding diminution of the WMDW. The fraction of AW was kept essentially unaffected. LIW would become the main TABLE 2 | Estimated trends of averaged buoyancy flux (Nm −2 s −1 /decade) during the period of available data in the monitoring station with the 95% confidence interval.
There are no numerical values in the literature for such a recent period to check with, but these trends in WMDW are clearly unrealistic. For instance, Vargas-Yañez et al. (2017) report 0.03 • C/decade and 0.02 decade −1 in the 600m-to-bottom layer in the Alboran Sea during 1943-2015; Houpert et al. (2016) indicate annual changes of 0.0032 ± 0.0005 • C and 0.0033 ± 0.0002 for θ and S during the period 2009-2013 in the 600-2300 m layer in the Gulf of Lions (equivalent to trends of 0.032 ± 0.005 • C/decade and 0.033 ± 0.002decade −1 ). They are one order of magnitude less than the rate necessary to validate the second hypothesis, which therefore should be discarded. Certainly, the trend could be partially due to the warming and salinification of LIW instead of WMDW. To this regard, a noticeable warming of LIW in the Gulf of Lions from 2007 to 2018 has been recently reported (Margirier et al., 2020). Even so, a consistent diminution of the proportion of WMDW in the Mediterranean outflow at the depth of the monitoring station is still necessary for that warming to manifest itself in the form of the observed trend in the Mediterranean outflow. Thus, the approach followed in this work is that the recent θ-S changes in the monitoring station (Figure 3) have been mainly caused by the progressive increase of LIW fraction in the outflow and the corresponding diminution of WMDW from around year 2013 onwards.
Estimates of PH and PCO 2 in the Outflow Water Masses Figure 5 illustrates a remarkable acidification of the Mediterranean outflow at the monitoring station, particularly from 2016 onwards. As in the case of θ discussed above, the observed pH T25 and pCO 2 trends could be the result of similar trends in the LIW and WMDW, the main components of the outflow, the change with time of the fractions of each water mass in the outflow, or both. In order to gain insights on the origin of this issue, the fraction analysis performed above (Figure 9, dashed-green box) was combined with the multi-linear regression model developed by Flecha et al. (2015) to obtain the individual pH T25 and pCO 2 values characterizing the LIW and the WMDW separately. This procedure has been applied successfully to discriminate between the components of the outflow and to estimate the values of oceanographic properties in each water body separately (Flecha et al., 2012(Flecha et al., , 2015(Flecha et al., , 2019. Figure 10 displays the time evolution of the pH T25 and pCO 2 assigned to each water mass in the outflow from 2012 to 2017. The analysis provides time-averages of 7.8718 ± 0.0004 and 7.9030 ± 0.0003 for pH T25 of the LIW and WDMW, respectively (Figure 10), values that are slightly lower than those reported by Flecha et al. (2015) who used a shorter time series spanning from 2012 to 2015. The difference can be interpreted as a gradual acidification attributed to the increase in CO 2 content in the Mediterranean waters. In fact, during the last portion of the available pH data (June 2015 to March 2017), the time-averages of the pH T25 series dropped to 7.8654 ± 0.0005 in the case of LIW and 7.8978 ± 0.0004 in WMDW, which were accompanied by a rise of pCO 2 (Figures 10C,D). The four-and-half-year time series in this study also revealed a striking and sudden decrease of pH T25 by the beginning of year 2016 in the outflow. It occurred concomitantly with a progressive rise in pCO 2 (Figure 5), the trends in both properties remaining relatively stable from then onwards. Signatures of the pH T25 decrease and pCO 2 rise can be indeed detected in the LIW and WMDW separately ( Figure 10) and therefore, it is plausible to assume that both water masses are contributing to the overall patterns detected in the outflow. However, water fraction analysis reveals that the warmer LIW prevailed in the outflow from January 2016 with the corresponding decline of WMDW fraction (Figure 9, green box), a situation more evident in the spring of this year. Fractions kept increasing (LIW) and decreasing (WMDW) rather steadily from September 2016 onwards.
According to former studies that used shorter series of this SAMI pH dataset (up to 2015(up to , Flecha et al., 2015 or discrete pH T25 measurements (Flecha et al., 2019), LIW was invariably characterized by pH T25 values lower than those in WMDW and relatively constant. These results agree with previous studies that report pH minima in the LIW in several Mediterranean subbasins in relation to the rest of Mediterranean waters (Alvarez et al., 2014;Hassoun et al., 2015). Such minima have been attributed to the remineralization of organic matter during aging of the LIW. As this water flows at intermediate depths isolated from the atmosphere in its way back to the Atlantic ocean, anthropogenic carbon absorption is prevented, making its pH levels remain stable (Alvarez et al., 2014). On the contrary, deep convection and/or shelf water cascade that drive the winter formation of WMDW from the very surface in the Gulf of Lions accumulates anthropogenic carbon in the newly formed water, which favors a WMDW acidification trend (Hassoun et al., 2015;Flecha et al., 2019). Patterns of pH T25 and pCO 2 up to 2015 in Figure 10 corroborate these findings: lower but relatively constant pH T25 in the LIW versus higher but decreasing-over time pH T25 in the WMDW. But Figure 10 also suggests a novel acidification trend in the LIW since 2016, which, along with the suggested prevalence of LIW fraction over that of WMDW in the observed outflow (Figure 9), would result in the pH T25 and pCO 2 trends visible in Figure 5 during the last one year and a half of data.

DISCUSSION AND CONCLUSION
Despite a few vicissitudes inherent to field work in harsh environments, the monitoring station deployed at Espartel sill in the SoG in late 2004 is providing worthy information about the outflow of the Mediterranean water and its properties. The outflow in turn reflects changes in the MedSea and provides an integrated overview of what is going on in the basin. Investigating likely links between changes in the MedSea properties and the outflow is the objective of this study.
A sequence of misfortunes in the experimental work caused an important data gap around year 2012. A detailed inspection of the series of collected and derived variables suggests two distinct periods in the series, with differentiated patterns of time evolution: a first period before the gap (2004-2012) and a second one after the gap (2013-2020). The most striking result is the trend of θ, which increased markedly starting by year 2013 and kept on up to the date of writing this work. In fact, the change was so outstanding (Figure 3) that not only motivated the inspection of the other variables searching for similar patterns, but also the present study.
The total outflow Q 2 showed an overall small positive trend of +0.016Sv/decade (Table 1), which reduces its size. During the first part, however, the trend was negative, −0.068 Sv/decade (the outflow increased), and then it changed sign after the gap in the second part ( Table 1). The contribution of the eddy-fluxes Q 2E−F to the total outflow deserves a little attention. In Espartel section Q 2E−F hardly reaches −0.04Sv or ∼5% (Table 1), in agreement with previous findings Sammartino et al., 2015). The percentage is very similar to the one in the eastern section of the SoG (García Lafuente et al., 2000;Baschek et al., 2001;Sannino et al., 2004). The smallness of Q 2E−F in these two sections points at a two-layer exchange hydraulically controlled at Espartel sill (instead of Camarinal) and Tarifa Narrows and the possibility of maximal exchange. Notwithstanding, the latter was not achieved (F 2 ≈0.74 < 1, Figure 7).
An interesting fact in this regards is the behavior of the two components of the total outflow, Q 2s−ν and Q 2E−F . In the first part of the series, the trend of Q 2 originated in Q 2s−ν , Q 2E−F having no contribution ( Table 1). In the second part, however, the positive trend of Q 2s−ν was noticeably higher than the trend of Q 2 and a negative trend of Q 2E−F was necessary to conciliate the mismatch. As Q 2E−F is negative, eddy-fluxes increased in magnitude, whereas the slowly-varying contribution diminished during this part. This is the expected behavior of each component of the flow if it separates away from the critical condition. Such a separation agrees with the decreasing trend of F 2 during the second part (Table 1 and Figure 7). Although F 2 showed overall tendency to decrease, that is, to accentuate the submaximal exchange (remind that F 2 ∼G due to the smallness of F 1 ), it had positive trend that pushed the flow towards critical during the first part. In the end, it did not reach it because of the sign change of the trend in the second part. In view of the results in Table 1, the conclusion is that F 2 mimics the behavior of Q 2 : a first part with Q 2 increasing in magnitude that drove F 2 to higher values until Q 2 changed trend in the second part, reversing F 2 trend as well.
In the simple steady-state model of the exchange that includes the MedSea, the rate of production of Mediterranean water, Q M , must equal the outflow |Q 2 |. Since the former depends on the buoyancy losses of the MedSea, buoyancy fluxes B 0 were computed to get some insight on the variability of Q M production ( Figure 8A). The estimated B 0 agrees with the available data of annual deep water formation in the Gulf of Lions (Houpert et al., 2016;Margirier et al., 2020), which supports the use of B 0 as the proxy of Q M production postulated by equation [3]. It must be recalled that the steady-state nature of the model limits its application to long-term changes. In other words, trends of B 0 that could be inferred from data in Figure 8 are expected to be echoed by similar trends in Q 2 according to the model, but year-to-year fluctuations do not have to fulfill this expectation. Therefore, the discussion is restricted to these long-term changes.
Even though no definitive pattern can be confirmed, there are reasonable indications that the two-part differentiation carried out in the oceanographic variables is still applicable to B 0 ( Table 2): a first part in which buoyancy fluxes tended to increase followed by a second one with the opposite trend. This second part sustained the three years of larger B 0 negative anomaly out of the sixteen year examined in the Gulf of Lions ( Figure 8B), with the expectable result of no (or weak) WMDW formation those years, as stated in Margirier et al. (2020).
All in all, the scenario that emerges is as follows. At the time of the beginning of the series, the outflow was submaximal (in terms of two-layer model) as a consequence of the weakness of winter buoyancy loss since 1988 in the northwestern MedSea (Herrmann et al., 2010;Margirier et al., 2020). It prevented strong convection and the formation of significant volumes of WMDW during the 1990s and early 2000s and enabled the increase of heat and salt content in the region. Within this scenario of preconditioned water, further favored by the arrival of the socalled Eastern Mediterranean Transient to the Gulf of Lions , the large buoyancy loss in 2004-2005 winter ( Figure 8A) propitiated the formation of an extraordinary volume of WMDW (López-Jurado et al., 2005;Schroeder et al., 2008), which left a clear footprint in the θ series of the monitoring station in early spring, 2005(García Lafuente et al., 2007. It was not the only consequence. The Q M production exceeded the outflow Q 2 at that time, forcing it to increase (a feasible result, since it was submaximal) and likely triggered the period of negative trend of Q 2 (greater outflow) found in our analysis ( Table 1). Successive winters of moderate or large volume of deep water formation (2012, for instance, as suggested by Figure 8A, see also Houpert et al., 2016) kept this trend until the production relaxed or, even, stopped. According to Houpert et al. (2016), there was no formation of deep water in year 2013, nor was there any appreciable formation between years 2014 and 2017 (Margirier et al., 2020). This situation of very few Q M formation caused |Q 2 | to diminish (particularly its Q 2s−ν contribution, Table 1) and reversed the trend, which changed to positive (reduced outflow) up to the end of the series. Froude number F 2 diminished accordingly, thus moving away from the critical value it had been approaching during the first part. The increased contribution of eddy-fluxes Q 2E−F to Q 2 ( Table 1) would agree with the augmented departure from the critical value, in the same manner as the relatively fast response of Q 2 to buoyancy changes point at submaximal exchange through the SoG (Garrett et al., 1990). This last aspect, which is applicable to the whole period, is in full agreement with the subcritical value exhibited by F 2 in Figure 7 throughout the series.
The time evolution of θ in Figure 3, whose sharp increase in the second part of the analyzed period motivated the present study, suggests a double origin for θ trends in Table 1. During the first part of the series, the estimated trend was +0.029 ± 0.009 • C/decade, which is very similar to trends reported in the literature for LIW or for WMDW (Vargas-Yañez et al., 2017;Houpert et al., 2016). But the one-order-magnitude greater trend in the second part (+0.339 ± 0.008 • C/decade) is well above any reported value in the Mediterranean Sea for deep waters. A more appropriate explanation for this trend is the progressive replacement of WMDW by the slightly warmer LIW in the outflow Q 2 , a situation that seems to persist nowadays. Biogeochemical data also agree with this hypothesis: pH T25 data indicate a noticeable acidification of the outflow (−0.0462 ± 0.0006 decade −1 ), which is particularly evident from year 2016. The decline of pH T25 with time and the concomitant rise in pCO 2 of LIW and WMDW detected when these variables were discriminated for both water masses (Figure 10), the striking CO 2 rise in the LIW from that particular year and its greater participation in the outflow, all of them contribute to shape the acidification trend observed in the monitoring station. The reasons for the continuous enrichment in carbon content of the LIW are out of the scope of this study, but some basin processes can be invoked. Although properties of deep waters suggest enhanced ventilation in the Gulf of Lions from early 2000 to 2016 (Li and Tanhua, 2020), deep water ventilation is highly variable in time and space. For instance, in the very severe winter of 2011-12 (Figure 8), the uplift of the LIW layer in this region caused by the massive formation of WMDW reported in Schroeder et al. (2016) would have brought this water mass closer to the surface and exposed it partially to atmospheric influence. The likely increment of anthropogenic carbon fraction associated with this rising, enhanced by the expected increase in respiration that occurs in heavily ventilated waters, would have resulted in a noticeable pH T25 decrease. On the other hand, Schroeder et al. (2016) suggest that WMDW formed in the Gulf of Lions may take around 33 months to reach the vicinity of the SoG, a figure that could be applicable to that LIW exposed to atmospheric effects. This time delay, which is also supported by age calculations carried out for both water masses in the Alboran Sea (Rhein and Hinrichsen, 1993;Flecha et al., 2019), would relate the pH T25 decrease observed in the time series in year 2016 with that cold event of winter 2012. Unrevealing the origin of changes in the properties of the carbon system detected in the LIW from 2016 onwards requires a thorough assessment of the whole bulk of biogeochemical properties recorded in the area over the last four years, particularly dissolved oxygen and organic matter content, a study that is presently under way. The analysis is particularly relevant due to the significant acidification rates recently measured in waters of the Levantine Basin (Hassoun et al., 2019), which will eventually reach the SoG and exit the basin as part of the outflow, with implications in the biogeochemistry of the North Atlantic.
The prevailing negative anomaly of buoyancy flux B 0 in the whole basin during the second part (green bars, Figure 8B) suggests a reduced formation of Mediterranean water Q M , which in turn gave rise to the concomitant reduction of the outflow Q 2 (Table 1), as discussed earlier. However, the B 0 anomaly was not homogeneously distributed over the basin. It was more accentuated in the Gulf of Lions (blue bars) than in the Levantine basin (red bars) where, in fact, the timeaveraged anomaly in the second part was nearly null. Therefore, the eventual reduction of Q M would have partially had its origin in the scarce production of WMDW in this period, confirmed by Margirier et al. (2020), with the collateral effect of leaving room for LIW to replace progressively the WMDW fraction in Q 2 , as suggested by Figure 9. In this Figure, the evolution of LIW fraction resembles the evolution of θ in Figure 3 because the variability of temperature prevails over that of salinity.
A remark must be made regarding this conclusion. Fractions in Figure 9 refer to the water outflowing at the depth of the monitoring station, close to the bottom, and should not be extrapolated to the rest of the Mediterranean layer above the station. If LIW contributes to the outflow more than WMDW, as it is generally accepted, the fractions above the depth of the station (but yet in the Mediterranean layer) will change to give LIW the dominant role. In the deepest portion of the outflow layer, though, the proportion of the WMDW used to prevail, as suggested by its fraction during the first part of the series (Figure 9; see also García Lafuente et al., 2007). But the recent data collected by the station indicated that, even at great depth, LIW would be becoming the prevailing water mass in the mixing. It is an expectable result if the diminished buoyancy flux to the atmosphere during the latest years were causing a WMDW production shortage that reduced its participation in the outflow. Harsh winters to come would halt this situation by increasing the formation rate of WMDW and, likely, reverse the recent outstanding θ trend observed in the monitoring station at Espartel section. Further research is needed to assess and confirm these changes as well as the mentioned acidification trend in the LIW from 2016 (and, hence, in the outflow through the SoG), which highlights the importance of maintaining the observation program whose data have propitiated this study.

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

AUTHOR CONTRIBUTIONS
JG-L designed the study and wrote the bulk of the manuscript. SS, IEH, and SF also contributed to the writing. SS, JG-L, CN, and IN processed the physical time series. IEH and SF acquired and processed biogeochemical data. RS-L and MJB organized most of the oceanographic surveys for servicing the monitoring station. All authors participated in the data acquisition and in the field campaigns and contributed to the interpretation, result discussion and revision of the manuscript. JG-L and IN carried out the manuscript edition. (CTM2006_02326/MAR), and INGRES3 (CTM2010_21229-C02-01/MAR), Special Action CTM2009-05810-E/MAR, funded by the Spanish Government. Since year 2016, the monitoring station is part of the Spanish Oceanographic Institute (IEO) internal project STOCA in which frame the station is presently serviced with the participation of Málaga University and researchers of the Interdisciplinary Thematic Platform of the CSIC WATER:iOS, with funding provided by the Ministry of Science and Innovation (EQC2018-004285-P) and the European Commission through the project COMFORT (H2020-820989). The mooring line is included in the Mediterranean Sea monitoring network of HYDROCHANGES project sponsored by CIESM, which is part of the CMEMS Marine Copernicus In Situ TAC Program, registered as WMO ID 6202100 in the JCOMM in situ Observations Program Support.