Climate change induced a new intermittent regime of convective ventilation that threatens the Black Sea oxygenation status

Abstract. The Black Sea is entirely anoxic, except for a thin (∼ 100 m) ventilated surface layer. Since 1955, the oxygen content of this upper layer has decreased by 44 %. The reasons hypothesized for this decrease are, first, a period of eutrophication from mid 70's to early 90's and, second, a reduction in the ventilation processes, suspected for the recent years. Here we show that the Black Sea convective ventilation regime has been drastically altered by atmospheric warming during the last decade. Since 2008, the prevailing regime is below the range of variability recorded since 1955, and is characterized by consecutive years during which the renewal of intermediate waters does not occur. Oxygen records from the last decade indicate a clear relationship between cold water formation events and oxygenation status at different pycnal levels, suggesting a leading role of convective ventilation in the oxygen budget of the upper intermediate layers. We thus suggest that this regime shift has a significant impact on the oxygenation structure of the Black Sea and on its biogeochemical balance.


Here, we combine different data sources to analyze the variability in the Black Sea intermediate layer ventilation over the last 65 years and in particular, investigate the existence of a statistically significant shift in the CIL formation regime, in regards to the variability observed over this extended period. The hypothesis of a significant regime shift is tested against the more traditional linear and periodic interpretation of the observed trends (e.g. Belokopytov, 2011), as the consequence for Black Sea ventilation and the future of the Black Sea oxygenation status in particular, are drastically different. Our analysis thus aims to 70 isolate annual convective ventilation as a particular component of the complex Black Sea deoxygenation dynamics (Konovalov and Murray, 2001;Capet et al., 2016) that typically involves intricate biogeochemical and physical processes.
Section 2 details the datasets considered to characterize the Black Sea CIL and oxygenation dynamics and the method of regime shifts analysis. In section 3, we analyse the long-term dynamics of the CIL through the lense of regime shift analysis.
In section 4, we use the most detailled datasets to relate CIL formation to changes in the Black Sea oxygenation conditions. 75 Section 5, then concludes with a short discussion on 2 Material and methods

The CIL cold content
While annual CIL formation rates are difficult to assess directly from observations, the status of the CIL can be quantified locally on the basis of vertical profiles of temperature and salinity. This simple indicator, based on routinely monitored vari-80 ables, provide a suitable metric to combine various sources of data while summarizing an essential aspect of the thermo-haline conditions. The CIL cold content C is defined as the heat deficit within the CIL, integrated along the vertical: with ρ, the in-situ density, c p , the heat capacity of sea water and T CIL = 8.35 • C, the temperature threshold which together with a density criterion ρ > 1014.5 kg m −3 , defines the CIL layer over which the integration is performed (Stanev et al., 2003, 85 2014; Capet et al., 2014). Although the use of a given temperature threshold to define the occurence of convective mixing is subject to discussion, the existence of a fixed threshold to define the CIL as a distinct watermass, and as a signature of its formation process, is evident given the fixed value of ∼ 9 • C that characterize the underlying deep waters (Stanev et al., 2019).
The above definition has been chosen for correspondence with past litterature.
the temperature at a fixed depth or the depth of a given iso-thermal surface. Although C is a deficit, we inverted the sign of C in comparison with previous literature (Stanev et al., 2003;Piotukh et al., 2011;Capet et al., 2014) for the convenience of working with a positive quantity. Large C values thus correspond to large heat deficit in the CIL, ie. to low temperature in a well-formed CIL layer, which is characteristic of cold years. A decrease in C corresponds to a weakening of cold water formation (typically for warm years), an increase in the intermediate water temperature and/or a decrease in the vertical extent 95 of the CIL.
C has been estimated for each year over the 1955-2019 period using four data sources summarized in Table 1. These sources include in-situ historical (ship-casts) and modern (Argo) observations, as well as empirical and mechanistic modelling ( Fig.   1). Annual and spatial average values for the deep sea (depth > 50 m) were derived from each data set, while considering the artifacts induced by uneven sampling in the context of pronounced seasonal and spatial variability.Each data source has 100 particular assets and drawbacks, and involves specific data processing to reach estimates of annual and spatial C averages as described below. All processed annual time-series are made available in netcdf format on a public repository (see 'Data availability').
Argo profilers: The assets of autonomous Argo profilers are a high temporal resolution and the continuous coverage of recent years, which offers unprecedented means to explore the CIL dynamics at fine spatial and temporal scales (Akpinar et al., 2017;105 Stanev et al., 2019). The drawbacks are the mingled spatial and temporal variability inherent to Argo data, the incomplete spatial coverage, and the lack of data prior to 2005. These data were collected and made freely available by the Coriolis project and programmes that contribute to it (http://www.coriolis.eu.org Argo profiles (Fig. 1). All available profiles in a given year were averaged to produce the annual Argo time series C Argo .
Although homogeneous seasonal sampling can be assumed, we note that the uneven spatial coverage of Argo profiles might induce a bias in the inferred trends. This potential bias stems from the horizontal gradient in C, that is structured radially from the central (lower C) to the peripheral (higher C) regions of the Black Sea (Stanev et al., 2003;Capet et al., 2014).
Argo samplings are more abundant in the peripheral regions, which suggests that C Argo might be slightly biased towards high 115 values.
In-situ ship-cast profiles: The asset of ship-based profiles is their extended temporal coverage. The drawbacks are the difficulty to untangle spatial and temporal variability (as for any non-synoptic data source), the uneven sampling effort and the low data availability posterior to 2000. The C Ships time series was provided by the application of the DIVA detrending methodology on ship-cast profiles extracted from the World Ocean Database (Boyer et al., 2009) in the box 40°-47°30' N, 120 27°-42°E for the period 1955-2011. DIVA is a sophisticated data interpolation software (Troupin et al., 2012) based on a variational approach. The detrending methodology (Capet et al., 2014) provided inter-annual trends, here representative for the central basin, cleared from the artifacts induced by the combination of uneven sampling and pronounced variability along  the seasonal and spatial dimensions. We redirect the reader to the reference study for further details on data sources, data distribution and methodology.

125
Atmospheric predictors: The statistical model considered here consists in a linear combination of time-lagged winter air temperature anomalies. This model was obtained by a stepwise selection amongst potential descriptor variables (incl. summer and winter air temperature, winds and fresh water discharge), in order to reproduce the inter-annual variability of C Ships (Capet et al., 2014) and proposed as an alternative to the winter severity index defined earlier (Simonov and Altman, 1991). C Atmos is thus naturally representative of the same quantity, i.e. annually and spatially averaged C. The asset of this approach is the 130 opportunity to fill the gaps between observations in past years, using atmospheric reanalysis of 2m air temperature provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) for the period 1980-2013. Its drawbacks lie in its empirical nature and indirect relationship to observable sea conditions. C Atmos was only extracted for the years covered in the reference study, considering that the potential non-linearity in the air temperature-C relationship may be exacerbated for the low C values typical of recent years.

135
3D hydrodynamic model: The Black Sea implementation of the 3D hydrodynamic model GHER has been used in several studies (Grégoire et al., 1998;Stanev and Beckers, 1999;Vandenbulcke et al., 2010;Capet et al., 2012). The present simulation uses the set-up presented in Capet et al. (2012) and has been produced without any form of data assimilation, on the basis of the ERA-Interim set of atmospheric forcing provided by the ECMWF data center (Dee et al., 2011). Aggregated weekly outputs of the GHER3D model are made available on a public repository (see 'Data availability'). C 3dM odel was derived from synoptic More precisely, the standard deviations estimated from the different series are similar (∼ 100 M J m −2 ), despite their distinct temporal coverage. The root mean square errors that characterize the disagreement between the different data sources remain 155 below this temporal standard deviation (in all but one case, see appendix A for details). This justifies to merge the different sources in a unique composite time series, enabling a robust long term analysis of the variablity in the Black Sea CIL formation.

Regime shift analysis and descriptive model selection
The inter annual variability of the Black Sea CIL formation is analyzed in the framework of regime shift analysis (Zeileis et al., 2003). The rationale behind this approach is to identify periods over which the times series depicts statistically distinct regimes, which point the addition of a new change point (ie. a new distinct period) in the regime shift model remains informative, we considered the Akaike Information Criterion (AIC), that assesses the RSS of each model but gives a penalty for the introduction of additional parameters (Akaike, 1974). The model with the smallest AIC should be favored for interpretation. The same 170 approach was used to compare the regime shifts models to the linear and periodic models.

Oxygen
BGC-Argo oxygen observations were obtained from the Coriolis data center for a period extending from 01/01/2010 to 01/01/2020. Only descending Argo profiles were considered, to minimize discrepancies associated with oxygen sensors response time (Bittig et al., 2014). To minimize the impact of spatial variability, oxygen concentrations were considered using

Descriptive models
The composite time series C i is depicted on Fig. 3 with calibrated linear and periodic models (with i as a year index).
The poor statistics associated with a linear model, C i = l 1 .i + l 2 (adjusted R 2 = 0.05, AIC=794), dismiss the perception of a linear trend extending over the entire period. Using the periodic model, gives a better representation of the cold content inter-annual variability (AIC=763), and provides broad characteristics of C i : the mean value, p 1 = 222 ± 190 12 M J m −2 , the amplitude of inter-annual variability, p 2 = 114 ± 16 M J m −2 , and the periodicity of pseudo-oscillations, p 3 = 43.04 ± 0.02 years.
However, both the linear and periodic models overestimates C in the recent years, as the composite time series C i shows a departure from its usual range of variability during the last decade. This is evidenced by ranking the 65 years of C i on the basis of their cold content. It is remarkable that, among the ten years with the least cold content, eight occurred after 2010.   Figure 4 gives the AIC obtained for regime shifts models considering from 1 to 5 change points (ie. partitioning the entire period into 2 to 6 segments). All of those regime shift descriptions appear as statistically more informative, sensu AIC, than a linear or periodic interpretation of the time series (and than a combination of linear and periodic model, AIC=758, not shown).
In particular the 4-segments model (ie. 3 change points) should be favored for interpretation.

Regimes shifts in the CIL cold content 200
The evolution of C i over 1955-2019 is thus best described by discriminating four periods (P 1 -P 4 , Fig. 5), objectively identified through regime shift analysis.  A "routine" regime is identified that is consistent for periods P 1  and P

CIL formation as a basin-wide ventilation process
The intra-annual resolution provided by the 3D model and Argo time series (Fig. 1) suggests that the CIL renewal, that was taking place systematically each year in precedent regimes, has now become occasional. Here we focus on the latest period, better detailed in our datasets, to characterize the CIL formation as a basin wide ventilation process and its relationship with 220 changes in oxygen concentration at different σ levels.
Basin wide CIL formation and destruction rates were computed from the synoptic 3D model outputs, as differences between weekly C values (Fig. 6). The seasonal sequence depicts CIL formation peaks from late December to March, typically reaching Following our attempt to summarize large datasets and to characterize a basin-scale annual oxygenation rate, we computed 235 the difference between the median oxygen concentration in November between successive years. The rationale behind this approach is that CIL formation typically extends from December to March (Fig. 6).
In order to obtain a general indication on the (pycnal) depth of penetration of the convective ventilation associated with CIL formation, we assessed the Pearson correlation coefficient between this annual oxygenation index, and a first order assessment of CIL formation, obtained as the annual difference in the C composite time series. The correlation between oxygenation 240 and CIL formation is high near the surface and decrease regularly as deeper pycnal levels are considered. Those correlations remains significant (p < 0.1) until σ = 15.4kg m −3 (Fig. 8).

Discussion
The regime shift paradigm describes an abrupt and significant change in the observable outcome of a non linear-system, as could result from a threshold in this system response to an external forcing. A periodic model, on the other hand, supposes either a 245 linear response to periodically varying external forcing, or an oscillation resulting from internal dynamics. It is our hypothesis, supported by the quantitative consideration presented above, that the first should be favored for interpreting the recent evolution of the Black Sea CIL dynamics. This difference in interpretation is fundamental in what regards the expected consequences on the Black Sea oxygenation status and in particular the threat on Black Sea marine populations, whose ecological adaptation  1981 1983 1985 1987 1989 1991 1993 1995 1997 1999 2001 2003 2005 2007 2009  due to the presence of reduced iron and manganese species, ammonium and finally hydrogen sulfide (Pakhomova et al., 2014).
Beyond the changes in convective ventilation highlighted above, it thus appears as a lead priority to assess the biogeochemical consequences of this new thermo-haline dynamics of the Black Sea. In particular, the influence of CIL formation on the biogeochemical components of the oxygen budget should be addressed in more details, asking for instance how the presence or absence of CIL formation influences on planktonic growth, trophic interactions, and organic matter respiration rates. We 13 https://doi.org/10.5194/bg-2020-76 Preprint. Discussion started: 30 March 2020 c Author(s) 2020. CC BY 4.0 License. voluntarily adopted here a wide integrative point of view, so as to highlight the large scale relevance of the depicted regime shift on Black Sea oxygenation. However, we still consider that a detailed assessment of all components of the oxygen budget, ie. ventilating processes and biogeochemical production/consumption terms, is required in order to infer the future evolution 265 of the Black Sea oxygenation status (e.g. Grégoire and Lacroix, 2001).
Although, there are clear indications of a long-term warming trend in the Black Sea (Belkin, 2009), it remains a delicate task to strictly dissociate the contributions of global warming from that of regional atmospheric oscillations (Kazmin and Zatsepin, 2007;Capet et al., 2012). One such assessment in the neighboring Mediterranean sea (Macias et al., 2013) concluded that global warming trend and regional oscillation contributed to the recent regional sea surface temperature trend , for 270 respectively 42% and 58%. While corresponding assessment will have to be routinely reevaluated for the Black Sea as time series expand, it may conservatively be considered that global warming had a significant contribution to warming winters in the Black Sea. This contribution is expected to increase in the next decades (Kirtman et al., 2014).

Conclusions
We have analyzed the variability of the Black Sea ventilation over the last 65 years and investigated the existence of regime 275 shifts in this dynamics. To this aim, we have produced a composite time series of the CIL cold content (C), that is considered as a proxy for the intensity of the convective ventilation resulting from the formation of dense oxygenated waters. This composite time series is built from four different data products issued from observations and modelling, so as to optimize its temporal extent. The consistency between those products, and in particular the close correspondence between observational and mechanistic predictive time series supports the reliability of the composite series and its adequacy to describe the evolution of the 280 Black Sea subsurface convective ventilation during the last 65 years.
The composite time series was analyzed to detect different regimes, corresponding to periods characterized by significantly distinct averages. Over the last 65 years, we identified 3 main regimes: 1) a routine regime prevailing during 1955-1984 and 1996-2008 that is consistent with the full-period average C, 2) a cold regime (high C, 1984-1996)  Given that cold winter air temperature is the leading driver of CIL formation (Oguz and Besiktepe, 1999;Ivanov et al., 2000;Capet et al., 2014), given that CIL formation constitutes 290 a predominant ventilation mechanism for the Black Sea intermediate layer, and assuming that oxygen conditions constitutes an environmental structuring factor affecting the ecosystem organization, its vigor and its resilience, this shift in the Black Sea ventilation regime may be associated with global warming and is expected to affect its biogeochemical balance and to threaten marine populations adapted to previously prevailing ventilation regimes.
fast and clear response that stems from the specific Black Sea geomorphology makes it a natural laboratory to study this dependency and related phenomena. Here, we showed that non-linear dynamics and feed-backs in ventilation mechanisms resulted in a significant shift of the average ventilation regime, in response to rising air temperature. Since the temporal extent of low oxygen conditions is critical for ecosystems, we stress the importance to assess the potential for similar ventilation regime shifts in other oxygen deficient basins.

300
Data availability. The data used are listed in Table 1 Argo data were were collected and made freely available by the Coriolis project 1 and programs that contribute to it. Era-Interim atmospheric conditions were obtained from ECMWF interface 2 . Aggregated weekly outputs of the GHER 3D model, as well as processed annual time series from the different sources are publicly available on a ZENODO repository : https://doi.org/10.5281/zenodo.3691960.
Appendix A: Comparison of the C) time series issued from different data sources

305
The C time series are noted C m i for source m (i is the year index). Each pair of time series (C m i , C n i ) are compared over the years i ∈ I m,n for which C m i and C n i are both defined. The following statistics are given for each pair of data sources in Fig.  A1 : -N m,n the number of elements in I m,n -Pearson correlation coefficient : The RMS difference between time series : -The percentage bias : For a better appreciation of variation scales, the temporal standard deviation is also shown for each data source.
The last value of the atmospheric predictor time series (C Atmos 2013 , Fig. 1. of the main manuscript) was not considered in the composite time series, as it was based on the two, rather than three, predictor values available at the time of publication (hence 320 the larger associated uncertainty). It is remarkable, however, that the published prognostic values for 2012 and 2013, match with independent Argo estimates.

Appendix B: Regime shift analysis
The identification of change points in a time series is the natural first step towards identification of a regime shift (Andersen et al., 2009). The basic change point problem can be expressed as follows: to identify the change point i = k in a sequence x i 325 of independent random variable with constant variance, such that the expectation of x i is µ for i < k and µ + ∆µ k otherwise.
Obviously, this problem can be generalized to several change points. The procedure for change point detection is stepwise and has been achieved following the methodology described in the documentation of the R package strucchange (Zeileis et al., 2003).
First, the existence of at least one significant change point had to be tested. strucchange provides seven statistical tests 330 to compute the p-value at which the null hypothesis of no change points can be rejected. The presence of change points can be tested on the basis of F-Statistic tests or generalized fluctuation tests (Zeileis et al., 2003, and references therein). Table A1 provides the significance level at which the null hypothesis can be rejected for each test implemented in the strucchange R package. Among the seven tests considered to assess the presence of at least one significant change point in C i , 6 rejects the null hypothesis with a significance level p < 0.05.

335
Second, the locations of the N most likely change points were identified. In this study, we considered from one to five change points. The change point locations can be estimated by finding the index values k n = [k 1 , .., k N ] that maximize a likelihood ratio, defined as the ratio of the residual sum of squares for the alternative hypothesis (i.e. change points at k n with ∆µ kn = 0) to that of the null hypothesis (i.e. no change point, ∆µ kn = 0).
Finally, the choice of a given model (ie. number of change points) had to be considered using the AIC criterion (Sect. 2.2).

Appendix B: Statistics
Formally, the methodology to identify and date structural change is designed for normal random variables, two conditions which can not be granted for environmental time series such as considered here. The following sections detail why 1) the departure of C distribution from a gaussian distribution, 2) the autocorrelation in C i , and 3) the biases between source-specific components of the composite time series C i , do not affect the conclusions drawn above.

B1 Normality
Skewness in the distribution of C values and its departure from normality is visible at low C values (not shown), as expected for physical reasons: C is a vertically integrated property, naturally bounded by zero. However, the Shapiro-Wilk test that measures the correlation between the quantiles of C i and those of a normal distribution, indicates no significant departure from normality: W = 0.975, p = 0.23. For completeness, a Box-Cox transformation (λ = 0.7) of the original data has been tested 350 which slightly enhances the Shapiro-Wilk test (W = 0.98, p = 0.39), but brings no sensible alteration in the conclusions of the structural change analysis.

B2 Autocorrelation
Similarly, C i can not be considered as a random variable. In particular, we introduced in Sect. 1 the CIL preconditioning and partial renewal mechanisms, both physical reasons for which autocorrelation may be expected in C i . Indeed, correlations 355 between the original and k-lagged time series are, at first glance, significant up to k = 5 (Tab. C1, the confidence interval above which autocorrelation can be considered to be significant is given by 1.96/ √ N = 0.25). However, it should be considered that the regime shift evidenced in the main manuscript may itself induce apparent autocorrelation statistics. To evidence that this is indeed the case encountered here, the correlation between the original and lagged time series of the residuals stemming from the 4-segments model are indicated in Tab. S C1. The fact that no significant autocorrelation persists when change points 360 are considered indicates that the non-randomness of C i does not jeopardize conclusions drawn from the application of the structural change methodology.

C1 Biases between components of the composite time series
Given that biases exist between different data sources, it might be argued that the composite time series is skewed by the uneven temporal cover of the different sources. For instance if a strongly biased source would solely cover a given period, the 365 composite series would be biased over that period. To ensure that this issue does not affect the presented conclusions, C unbiased i was constructed as C i , but removing from each component C m i the bias identified identified with the longest C