Recent decline of the Black Sea oxygen inventory

Recent decline of the Black Sea oxygen inventory A. Capet, E. V. Stanev, J.-M. Beckers, J. W. Murray, and M. Grégoire OGS, National Institute of Oceanography and Experimental Geophysics, Trieste, Italy Laboratory of Oceanology, University of Liège, Liège, Belgium HZG, Helmholtz-Zentrum Geesthacht, Hamburg, Germany GHER, GeoHydrodynamics and Environment Research, University of Liège, Liège, Belgium School of Oceanography, University of Washington, Seattle, WA, USA


Introduction
The Black Sea deep waters constitutes the world's largest reservoir of toxic hydrogen sulfide. 100 m of ventilated surface waters are all that separate this reservoir from the atmosphere. This situation results from the permanent halocline (Öszoy and Ünlüata, 1997) that separates the surface layer (of low salinity due to river inflow) from the deeper layer (of high salinity due to inflowing Mediterranean seawater), restraining ventilation to the upper layer ( Fig. 1).
In the lower part of the halocline, a permanent suboxic layer separates the Black Sea surface oxygenated waters ([O 2 ] > 20 µM) from the deep sulfidic waters ([H 2 S] > 20 µM) (Murray et al., 1989(Murray et al., , 1995Tugrul et al., 1992). More precisely, Murray et al. (1989) considered a threshold of 10 µM of oxygen because they analyzed highquality oxygen data. The threshold of 20 µM of oxygen was applied later to analyze historical oxygen data of lower quality. The upper (O 2 disappearance) and lower (H 2 S onset) interfaces of this suboxic layer are controlled by different biogeochemical and physical processes Stanev et al., 2014), and undergo uncorrelated vertical migrations (Konovalov and Murray, 2001). Sinking organic matter is mainly respirated aerobically within the oxycline: the lower part of the oxygenated layer where oxygen concentration decreases downwards to 20 µM. Increasing flux of organic matter, induced by a period of high nutrient load from the 1970s to the late 1880s, resulted in higher oxygen consumption above the suboxic layer and a shoaling of the upper suboxic interface (Codispoti et al., 1991;Konovalov and Murray, 2001;Tugrul et al., 2014).
After reduction of nutrient inputs around 1990 (Kroiss et al., 2006), the Black Sea was described as a recovering ecosystem (Mee et al., 2005;Oguz et al., 2006). This perspective was supported by improved eutrophication indices in the open sea (Kideys, 2002) as well as the stabilization of the upper suboxic interface in the 1990s (Konovalov and Murray, 2001). However, the timescale of the expected recovery (i.e., the timescale associated with the chain of biogeochemical mechanisms relating oxycline penetration depth to riverine nutrient loads) is not quantitatively understood. Several processes cause the oxycline depth to respond with a time lag to the reduction of riverine nutrient inputs. First, nutrients are mainly delivered to the northwestern shelf, where the accumulation of organic matter in the sediments buffers the riverine inputs, with slow diagenetic processes controlling and delaying the nutrient outflow across the seaward boundary (Capet et al., 2013). Second, the intermediate oxidation-reduction cycling of nitrogen, sulfur, manganese, iron and phosphorus that separates oxygen from hydrogen sulfide (Shaffer, 1986;Codispoti et al., 1991;Konovalov et al., 2006;Yakushev et al., 2007) can delay the response of the lower suboxic interface to changing nutrient fluxes by several years .
In addition to these biogeochemical factors, the dynamics of the upper and lower interfaces of the suboxic layer are controlled by physical processes Stanev et al., 2014). In the Black Sea, dense waters formed by winter cooling and mixing (Staneva and Stanev, 2002) do not sink to the deepest layer, as in the Mediterranean Sea, but accumulate on top of the permanent halocline (Fig. 1). The resulting cold intermediate layer (CIL) is a major feature of the Black Sea vertical structure. Cold intermediate water formation and advection by the cyclonic basin-wide Rim Current (Öszoy and Ünlüata, 1997;Capet et al., 2012) ventilate the oxycline and thereby influence variability in the depth of the upper suboxic interface . Recently, atmospheric warming (Oguz et al., 2006) was shown to reduce the ventilation of the lower oxic layer (Tugrul et al., 2014;Pakhomova et al., 2014). At deeper levels, the dense sinking plume formed by the Mediterranean inflow through the Bosporus, which entrains water from the overlying CIL, injects fingers of oxygenated water directly into the deeper part of the suboxic layer and upper sulfidic layer and thus acts to control the depth of the lower suboxic interface (Konovalov and Murray, 2001;Konovalov et al., 2003Konovalov et al., , 2006Glazer et al., 2006).
Previous long-term analyses of the vertical migration of the suboxic interfaces either ended (1955-1995Konovalov and Murray, 2001) or started (1985-2015Pakhomova et al., 2014) with the eutrophication period, excluding the largescale overview required to grasp the interactions of eutrophication and climate factors. Those analyses lacked a comprehensive consideration of the natural spatial and seasonal variability of the vertical distribution of oxygen. In the presence of large gradients, uneven data distribution may induce artificial signals when interannual trends are assessed from direct annual averages. In the stratified Black Sea, properties expressed in terms of depth coordinates (m) present a high spatial variability due to mesoscale features (Kempe et al., 1990) and to the general curvature of Black Sea isopycnals (Öszoy and Ünlüata, 1997;Stanev et al., 2014). As an alternative, using density (isopycnal levels, σ θ ) as vertical coordinate is generally considered a stable solution to assess the vertical migration of the chemocline on a decadal scale (Tugrul et al., 1992;Saydam et al., 1993;Murray et al., 1995). However, the spatial confinement of the lateral oxygen injections associated with the Bosporus plume, as well as the spatial variability of diapycnal ventilating processes (Zatsepin et al., 2007), imposes a horizontal structure to the oxygen penetration depth when expressed in terms of density (Stanev et al., 2004;Glazer et al., 2006). As this spatial gradient might scale with the temporal variations (a range of 0.17 kg m −3 was observed during the Knorr 2003 campaign, Glazer et al., 2006), it has to be considered when deriving interannual trends.
The present study describes the application of the DIVA (Data-Interpolating Variational Analysis) detrending procedure (Troupin et al., 2012;Capet et al., 2014) to untangle the temporal and spatial variability of three indices related to the Black Sea oxygenation status: the depth and density level of oxygen penetration and the oxygen inventory. These values were diagnosed from a composite historical data set of oxygen vertical profiles. We review the evolution of those indices through the past 60 years and discuss the respective controls of eutrophication and climate factors.

Data
We gathered a composite set of 4385 ship-based vertical profiles (oxygen, temperature and salinity) obtained between 1955 and 2005 in the Black Sea using CTD rosette bottles, continuous pumping profilers (Codispoti et al., 1991) and in situ analyzers (Glazer et al., 2006)  To complement the analysis of ship-based casts, we considered profiles originating from 10 Argo autonomous profilers (May 2010-December 2015. Only good quality-checked real-time data were considered (Carval et al., 2014). Two of these floats (Argo ID 7900465 and 7900466) have been presented and discussed by Stanev et al. (2013), where the consistence and comparability of Argo and historical profiles is asserted within a 1 µM error range.
Several studies address the error of Argo real-time oxygen data (e.g., Bittig and Körtzinger, 2015;Takeshita et al., 2013;Johnson et al., 2015). Demonstrating that the Black Sea real-time Argo data are precisely (i.e., at fine scales) comparable with historical Winkler data, or identifying the relevant correction, is beyond the scope of the present study which addresses monthly to decadal timescales. Evenly distributed small-scale error (e.g., difference between ascending and descending profiles due to sensor time response) were thus filtered by the temporal smoothing. However, a systematic error is not strictly excluded which could reach an under- estimation of 10 µM (Virginie Thierry, IFREMER, personal communication, January 2016). Therefore, we evaluated a "worst-case" scenario in the analysis of Argo data by considering a systematic underestimation of oxygen concentration by 10 µM.
Although most of the floats drifted along the basin periphery, some were also advected in the central part (Fig. 3). These trajectories highlight the range of spatial variability for the diagnostics described in Sect. 2.2.

Profile analysis
From each profile we derived (1) the depth and (2) the potential density anomaly σ θ where oxygen concentration went below 20 µM and (3) the oxygen inventory, integrated above this limit (Fig. 1). The threshold value of 20 µM used to define the upper interface of the suboxic layer was suggested to compare oxygen observations issued from sensors with different detection limits (Konovalov and Murray, 2001). To evaluate how a 10 µM underestimation by Argo profilers would affect the main conclusions, oxygen penetration depths an density levels for Argo were also computed using a threshold of 10 µM.
The CIL cold content ( Fig. 1) was diagnosed from corresponding salinity and temperature profiles following Capet et al. (2014). It indicates on the intensity of CIL formation smoothed over 4-5 years, i.e,. the residence time of cold intermediate waters (Staneva and Stanev, 2002;Piotukh et al., 2011;Capet et al., 2014): where ρ is the density and c the heat capacity and T CIL = 8.35 • C (Stanev et al., 2013).

DIVA analysis
Climatologies for the whole period and interannual trends were identified for the three oxygen diagnostics by applying the DIVA detrending algorithm on the ship-based data set (see details in Appendix A). In short, the DIVA interpolation software (http://modb. oce.ulg.ac.be/mediawiki/index.php/DIVA; Troupin et al., 2012) computes a gridded climatology obtained by minimizing a cost function which penalizes gradients and misfits with observations. The DIVA detrending algorithm (Capet et al., 2014) computes trends for each year, i.e., the average difference between data pertaining to this year and the spatial analysis at these data locations. This procedure allows one to account for the sampling error associated with spatial/temporal variability.  , accounting for the temporal variability of these diagnostics and the uneven distribution of data (see Sect. 2.3).

Spatial variability
The spatial distribution of the oxygen penetration depth (Fig. 4a) reflects the general curvature of the Black Sea vertical structure. A range of approximately 70 m was observed between oxygen penetration depth in the periphery (150 m) and in the central part (80 m).
A significant spatial variability remains when expressing oxygen penetration in terms of potential density anomaly σ θ (Fig. 4b) The spatial distribution of the oxygen inventory (Fig. 4c) follows that of the oxygen penetration depth. The range of spatial variability reaches 12 mol O m −2 , i.e., between 17 mol O m −2 in the central part and 29 mol O m −2 in the periphery.
The ranges of spatial variability derived from these spatial analyses agreed with those depicted by the Argo profilers (Fig. 5), bearing in mind the different timescales under consideration.

Temporal variability
Between 1955 and 2005, the oxygen penetration depth rose by an average rate of 7.9 m per decade (Fig. 6a) The oxygen inventory, integrated from the surface down to the suboxic upper interface, decreased by 44 % during the last 60 years (Fig. 6c), considering the ship-based estimate for 1955 (27 mol O m −2 ) and the Argo estimate for 2015 (15 mol O m −2 ). The few ship-based profiles available after the mid-1990s revealed the lowest oxygen inventories recorded during the time frame covered by the present study.
The temporal signals departed from these linear trends between 1988 and 1996, during which deeper oxygen penetration (both in terms of depth and density) and higher oxygen content were observed.

Oxygen inventory and CIL cold content
Positive relationships between oxygen inventory and CIL cold content were obtained for all periods (Fig. 7). Considering a given level of CIL cold content, the corresponding oxygen inventory decreased significantly from period 1955-1975 to period 1986-1998 (Fig. 7b).
The relationship between oxygen inventory and CIL cold content for the period 1999-2015 does not differ significantly from that obtained for the period 1986-1998 (Fig. 7b). This comparison should be considered with caution, however, as oxygen profiles for the period 1999-2015 originate mainly from Argo floats whose sampling rate is much higher than ship-based casts.
High CIL cold content is much more frequent during the period 1986-1998, while low CIL cold content is more frequent during 1999-2013.  Fig. 6).

Discussion
The spatial analysis of oxygen penetration depth showed that the use of density coordinates does not eliminate the sampling error associated with uneven spatial coverage (Fig. 4). Deeper oxygen penetration (on a density scale) in the Bosporus area were expected, in relation with the intermediate lateral injections associated with the Bosporus plume. In addition, deeper oxygen penetration (on a density scale) in the southern and eastern periphery suggests the occurrence of diapycnal ventilation along the steep bathymetry (Zatsepin et al., 2007). The aggregation of the most recent ship-based profiles in the Bosporus area and in the southeastern region (Fig. 3), might have led to an overestimation of the basinaverage oxygen penetration depth in the last decade, hence to an underestimation of the shoaling trend of the Black Sea oxic layer.
Considering spatial variability revealed a clear shoaling trend for oxygen penetration depth. This shoaling can be seen on both depth and density scales (Fig. 6a, b). This confirms the hypothesis that the shoaling of oxygen penetration depth is not due to a general shoaling of the main halocline, but is associated with a shifted biogeochemical balance in the oxygen budget (Codispoti et al., 1991;Konovalov and Murray, 2001;Tugrul et al., 2014).
Using σ θ coordinates depicts clearer temporal variations (Figs. 5 and 6). The shoaling rate varies in time and was more intense during 1970-1985 and from 1996 onwards. Argo diagnostics using different oxygen threshold show a larger discrepancy in the case of pycnal coordinates. The cooccurrence of density and oxygen gradients (Fig. 1) results in a higher sensitivity to the sensor accuracy for the σ θ diagnostic for oxygen penetration. However, even a systematic underestimation by 10 µM of oxygen concentration by Argo profilers does not invalidate our results.
The positive correlations between CIL cold content and oxygen inventory observed for all the periods illustrate the ventilation of intermediate layers by CIL formation and advection (Fig. 7b). In the early 1990s, the transient recovery of the three oxygenation diagnostics (Figs. 6a, b, c, 7a) provided arguments supporting the stability of the oxic interface (Tugrul et al., 1992;Buesseler et al., 1994). This stabilization matched the convenient perception of a general recovery of the Black Sea ecosystem after the reduction of nutrient load around 1990 (Kroiss et al., 2006). However, Fig. 7 indicates that the oxygenation diagnostics obtained for the period 1986-1998 were associated with much higher ventilation rates (i.e., higher CIL cold content) than during the previous periods. If, in response to nutrient reduction, the biogeochemical oxygen consumption terms had been lower during the period 1986-1998 than previously, the increased ventilation during that period would have resulted in higher oxygen inventories. Instead, oxygen inventories observed during 1986-1998 are lower than those observed in the previous decade for similar levels of CIL cold content. We conclude that high CIL formation rates during this period (Piotukh et al., 2011;Capet et al., 2014) provided enough ventilation to mask ongoing high oxygen consumption.
The fact that the relationship between oxygen inventories and CIL content for the last period 1999-2015 is similar to that of 1986-1998 indicates a stabilization in the biogeochemical oxygen consumption terms. Higher air temperature in this last period (Oguz and Cokacar, 2003;Oguz et al., 2006;Pakhomova et al., 2014), by limiting winter convective ventilation events (Capet et al., 2014) led to the lowest oxygen inventories ever recorded for the Black Sea (Fig. 6c).
Forecasted global warming, without excluding transient high ventilation periods, will limit CIL water formation (Capet et al., 2014) and reduce the oxygenation of the Black Sea intermediate layers. At the same time, uncertainties remain regarding the capacity of re-flourishing economies of the lower Danube watershed to recover their productivity in a more sustainable, less polluting form. Economic development in the Danube Basin could reverse the improving situation of eutrophication if nutrients are not managed properly (Kroiss et al., 2006). Under these conditions, there is no reason to expect that the oxycline shoaling observed over the past 60 years will stabilize.
There are reasons to worry about a rising oxycline in the Black Sea. First, biological activity is distributed vertically on the whole oxygenated layer, as indicated by zooplankton dial migration (Ostrovskii and Zatsepin, 2011). The reduction of the oxygenated volume described in this study could therefore have impacted on Black Sea living stocks by reducing carrying capacity and increasing predation encounter rates. It would be timely to estimate now the impact that a further shoaling of the oxic interfaces would bear on the Black Sea resources for the fishing industry.
Second, under present conditions, a massive atmospheric release of hydrogen sulfide caused by a sudden outcropping of anoxic waters remains unlikely, due to the stability of the Black Sea pycnal structure. Such outcropping event of sulfidic waters would have dramatic ecological and economical consequences (Mee, 1992). On 27 September 2005, an anomalous quasi-tropical cyclone was observed over the western Black Sea that led, in a few days, to the outcropping of waters initially located at 30 m depth (Efimov et al., 2008). Two years earlier, sulfide was measured in the same area (western central gyre) at around 80 m (Glazer et al., 2006). Because global warming is expected to increase the occurrence of extreme meteorological events (Beniston et al., 2007), every meter of oxycline shoaling would bring the Black Sea chemocline excursion events closer to the realm of possibility.  Cleveland et al., 1992) between oxygen inventory and CIL cold content for the different periods (confidence interval α = 0.99). The positive relationships observed during each period illustrate the ventilating action of CIL formation as a source of oxygen to the intermediate levels. The shift of these relationships towards lower oxygen inventories indicates shift in the oxygen budgets (higher consumption) that are independent of the intensity of CIL formation.

Conclusions
The present study evidenced the decline of the Black Sea oxygen inventory during the second half of the 20th century and first decade of the 21st, and highlighted the threat that further atmospheric warming casts upon the vertical stability of the Black Sea oxygenated layer.
Further works are urgently required to assess how actual nutrient emission policies adequately prevent, in the context of forecasted warming, the ecological and economical damages that would arise from a further shoaling of the oxic interface.
Spatially resolved biogeochemical models are needed to integrate explicitly the interacting processes affecting the Black Sea oxycline.
It is also essential (1) to determine to which extent the shoaling of the oxygen penetration depth entrains a shoaling of the sulfidic onset depth, (2) to set up a continuous monitoring of the Black Sea oxygen inventory and the intensity of winter convective ventilation (through CIL cold content), and (3) to clarify and quantify the interplay of diapycnal and isopycnal ventilation mechanisms and, in particular, the role played by the peripheral permanent/semipermanent mesoscale structures and how this relates to the intensity of the Rim Current (Stanev et al., 2014;Kubryakov and Stanichny, 2015). We propose that these objectives might be answered by maintaining in the Black Sea a minimum population of both moored and drifting autonomous profilers equipped with oxygen and sulfidic sensors.
Appendix A: The DIVA detrending algorithm DIVA (Data-Interpolating Variational Analysis) is a method for spatial interpolation. Its principle is to construct an analyzed field ϕ that satisfies a set of constraints expressed in the form of a cost function over a spatial domain . The cost function is made up of (1) an observation constraint, which penalizes the misfit between data and analysis, and (2) a smoothness constraint, which penalizes the irregularity of the analyzed field (gradients, Laplacian, etc.).
Let us assume that we work with data anomalies, i.e., a reference (or background) field is subtracted from the data points prior the analysis. For N data anomalies d i at locations (x i , y i ), the cost function reads, in Cartesian coordinates, where µ i , α 0 and α 1 are coefficients related to characteristics of the data set. ∇ is the horizontal gradient operator and ∇∇ϕ : ∇∇ϕ = i j (∂ 2 ϕ/∂x i ∂x j ) (∂ 2 ϕ/∂x i ∂x j ), the generalization of the scalar product of two vectors. The first term of Eq. (A1) measures the spatial variability (curvature, gradient and value) of the analyzed field and is identified as the smoothness constraint. The second term is a weighted sum of data-analysis misfits and is identified as the observation constraint: it tends to pull the analyzed field towards the observations. The coefficients of Eq. (A1) can be determined from (1) the relative weights w i attributed to each observation d i , (2) the correlation length L and (3) the signalto-noise ratio λ (Troupin et al., 2012). The analyses presented in this study were achieved with equal weights w i = 1, L = 0.8 • and λ = 0.5. The minimization of Eq. (A1) is solved over with a finite-element technique (Brasseur et al., 1996) which excludes data influence across land points (Troupin et al., 2010).
The detrending algorithm, presented in Capet et al. (2014) with synthetic and real case studies, proceeds as follows. Input data can be classified amongst the different classes C j (e.g., 1990, 1991, . . . ) of a given group C (e.g., the year). The observation constraint of the functional Eq. (A1) can then be rewritten by including an unknown trend value for each class (d C 1 , d C 2 , . . . ): If the function ϕ(x, y) were known, minimization with respect to each of the unknowns d C j would yield and similarly for the other classes: the trend for each class is the weighted misfit of the class with respect to the overall analysis.
Using an analysis without detrending as a first guess for ϕ, trends are computed for each classes in each group and subtracted from the original data. Following this, a new analysis is performed, the trends are recalculated, and the iterations continue until a specified convergence criterion is fulfilled. The procedure can be generalized with several groups of classes. The present study considered years and months.
The DIVA software and up-to-date related information can be found on http://modb.oce.ulg.ac.be/mediawiki/index.php/ DIVA.