Drivers , mechanisms and long-term variability of seasonal hypoxia on the Black Sea northwestern shelf – is there any recov ry after eutrophication ?

The Black Sea northwestern shelf (NWS) is a shallow eutrophic area in which the seasonal stratification of the water column isolates the bottom waters from the atmosphere. This prevents ventilation from counterbalancing the large consumption of oxygen due to respiration in the bottom waters and in the sediments, and sets the stage for the development of seasonal hypoxia. A three-dimensional (3-D) coupled physical– biogeochemical model is used to investigate the dynamics of bottom hypoxia in the Black Sea NWS, first at seasonal and then at interannual scales (1981–2009), and to differentiate its driving factors (climatic versus eutrophication). Model skills are evaluated by a quantitative comparison of the model results to 14 123 in situ oxygen measurements available in the NOAA World Ocean and the Black Sea Commission databases, using different error metrics. This validation exercise shows that the model is able to represent the seasonal and interannual variability of the oxygen concentration and of the occurrence of hypoxia, as well as the spatial distribution of oxygen-depleted waters. During the period 1981–2009, each year exhibits seasonal bottom hypoxia at the end of summer. This phenomenon essentially covers the northern part of the NWS – which receives large inputs of nutrients from the Danube, Dniester and Dnieper rivers – and extends, during the years of severe hypoxia, towards the Romanian bay of Constanta. An index H which merges the aspects of the spatial and temporal extension of the hypoxic event is proposed to quantify, for each year, the intensity of hypoxia as an environmental stressor. In order to explain the int rannual variability of H and to disentangle its drivers, we analyze the long time series of model results by means of a stepwise multiple linear regression. This statistical model gives a general relationship that links the intensity of hypoxia to eutrophication and climater lated variables. A total of 82 % of the interannual variability of H is explained by the combination of four predictors: the annual riverine nitrate load ( N ), the sea surface temperature in the month preceding stratification ( Ts), the amount of semi-labile organic matter accumulated in the sediments ( C) and the sea surface temperature during late summer ( Tf). Partial regression indicates that the climatic impact on hypoxia is almost as important as that of eutrophication. Accumulation of organic matter in the sediments introduces an important inertia in the recovery rocess after eutrophication, with a typical timescale of 9.3 yr. Seasonal fluctuations and the heterogeneous spatial distribution complicate the monitoring of bottom hypoxia, leading to contradictory conclusions when the interpretation is done from different sets of data. In particular, it appears that the recovery reported in the literature after 1995 was overestimated due to the use of observations concentrated in areas and months not typically affected by hypoxia. This stresses the urgent need for a dedicated monitoring effort in the Black Sea NWS focused on the areas and months concerned by recurrent hypoxic events. Published by Copernicus Publications on behalf of the European Geosciences Union. 3944 A. Capet et al.: Seasonal hypoxia in the Black Sea northwestern shelf

In order to explain the interannual variability of H and to disentangle its drivers, we analyze the long time series of model results by means of a stepwise multiple linear regression.This statistical model gives a general relationship that links the intensity of hypoxia to eutrophication and climaterelated variables.
A total of 82 % of the interannual variability of H is explained by the combination of four predictors: the annual riverine nitrate load (N ), the sea surface temperature in the month preceding stratification (T s ), the amount of semi-labile organic matter accumulated in the sediments (C) and the sea surface temperature during late summer (T f ).Partial regression indicates that the climatic impact on hypoxia is almost as important as that of eutrophication.
Accumulation of organic matter in the sediments introduces an important inertia in the recovery process after eutrophication, with a typical timescale of 9.3 yr.
Seasonal fluctuations and the heterogeneous spatial distribution complicate the monitoring of bottom hypoxia, leading to contradictory conclusions when the interpretation is done from different sets of data.In particular, it appears that the recovery reported in the literature after 1995 was overestimated due to the use of observations concentrated in areas and months not typically affected by hypoxia.This stresses the urgent need for a dedicated monitoring effort in the Black Sea NWS focused on the areas and months concerned by recurrent hypoxic events.

Introduction
Hypoxia is a major global environmental problem affecting marine waters and is among the most widespread deleterious anthropogenic influences.The phenomenon of hypoxia is drawing increased attention as dead zones in the coastal ocean have spread exponentially since the 1960s and have been reported from more than 400 sites (Diaz and Rosenberg, 2008) (e.g., Baltic Sea, Black Sea, Kattegat Sea, Gulf of Mexico, East China Sea), 39 % of which are located in Europe (Zhang et al., 2010).
Hypoxia is considered to occur when the concentration of dissolved oxygen decreases below a threshold of 62 mmol m −3 (= 2 mg L −1 ).However, it has been shown that low oxygen concentrations already affect living organisms above that threshold (Vaquer-Sunyer and Duarte, 2008).
Hypoxic conditions may be encountered in various parts of the ocean (e.g., the coastal area, the shelf and the deep ocean), where water ventilation is not able to renew the oxygen consumed for the degradation of organic matter.The coastal area which receives the input of organic matter from rivers and adjacent land is particularly exposed to this phenomenon.In this shallow area, where the bottom is occupied by a benthic ecosystem, the occurrence of hypoxic conditions may be lethal to these organisms and catastrophic for the benthic biodiversity (Levin et al., 2009).
Hypoxia has important impacts on biogeochemical cycling and may significantly disturb ecosystem functionality as, for instance, it -directly affects living organisms (e.g., benthic organisms, fishes) (e.g., Vaquer-Sunyer and Duarte, 2008), -may, in extreme cases, lead to anoxia with the production of greenhouse gases (e.g., CH 4 , N 2 O), or toxic substances (e.g., sulfide) (e.g., Naqvi et al., 2010), -alters the cycling of chemical elements such as nitrogen or phosphate (e.g., Conley et al., 2009), -modifies the sedimentary geochemical cycling through the removal of bioturbating infauna (e.g., Levin et al., 2009), -and alters the food web structure by changing the balance of chemical elements (e.g., N, C, P) and by killing some components and hence reducing the transfer of energy towards the higher trophic levels (e.g., Ekau et al., 2010).
The phenomenology of hypoxia events may be classified according to the temporal scales of hypoxia (e.g., permanent, seasonal, episodic, diel) and geomorphological characteristics of the water body (Kemp et al., 2009).The Black Sea northwestern shelf (NWS) is stratified, non-tidal, and affected by recurrent seasonal hypoxia.According to the classification of Kemp et al. (2009) the situation of the Black Sea NWS system is comparable with, for instance, Chesapeake Bay (e.g., Boesch et al., 2001), Pensacola Bay (Hagy III and Murrell, 2007), the Changjiang plume (Chen et al., 2007), and the northern Gulf of Mexico (Turner et al., 2008(Turner et al., , 2012)).
While it is generally understood that hypoxia in those systems is permitted by stratification and fueled by humaninduced eutrophication, non-linearities, indirect and feedback effects involving a number of external and internal factors complicate the situation and prevent direct links between a driver and the intensity of the resulting hypoxia (Kemp et al., 2009;Steckbauer et al., 2011).For instance, we may expect that a temperature rise will intensify hypoxia by strengthening stratification, reducing oxygen solubility and enhancing respiration rates.However a temperature change may also have indirect effects on hypoxia: an increased stratification may for instance reduce the availability of nutrients in the surface layer (and hence limit primary production) and/or intensify the respiration in the upper layer, which would actually reduce the oxygen consumption in bottom waters (Conley et al., 2009).
The organic matter accumulated in the sediments during eutrophication periods introduces an inertia in the response of hypoxia to reduced nutrient load, whose timescale is difficult to estimate although essential to the study of the system recovery.Moreover, hypoxic conditions affect the sediment biogeochemistry by reducing denitrification (Kemp et al., 2009) and causing the release of dissolved inorganic phosphorus from ferric complexes (Slomp and Van Cappellen, 2007).These two processes change the nutrient content of the upper water and hence the level of primary production that can be sustained.
Finally, the modified food web structure of ecosystems affected by hypoxia may also alter the response of the ecosystem to nutrient load (e.g., filter feeders, bioturbators, fish/jellyfish trophic pathways).
For all the reasons mentioned above, which are not exhaustive, hypoxia is a complex process whose dynamics (spatial and temporal scales), understanding (e.g., drivers, impact) and mitigation (e.g., recovery) are challenging to assess (Cardoso et al., 2010).Sound management allowing the preservation and restoration of hypoxia-damaged ecosystems to a healthy state therefore necessitates a good understanding of these interconnected processes (Boesch et al., 2001).
In situ data collection is a key component for this assessment but shows some limitations (Cardoso et al., 2010) that can be alleviated when combined with numerical modeling.Indeed, coupled physical-biogeochemical models are valuable tools for investigating some mechanisms such as the effects of the physics (governed by atmospheric conditions)and biogeochemistry (including the sediments) on the occurrence of hypoxia (Kemp et al., 2009;Peña et al., 2010).
In this paper, a three-dimensional (3-D) physicalbiogeochemical model is used to investigate the space-time dynamics of hypoxia and its drivers (atmospheric versus eutrophication) on the Black Sea NWS.
The Black Sea is an almost enclosed sea with restricted exchanges with the Mediterranean Sea through the Bosphorus Strait.It contains a wide shelf (depth < 100 m) on its northwestern part, connected to the deep sea (depth 2000 m) (Fig. 1).The deep sea is permanently stratified by a halocline, and hence water ventilation in this area does not extend deeper than ∼ 100-150 m, leading to anoxia of waters below that depth.Moreover, hypoxia on the NWS was first observed in the early 1970s at 30 m depth in a region located between the Danube and Dniester discharge area (Zaitsev, 1997).
The losses of bottom macro-fauna attributed to hypoxia for the 1972-1990 period have been estimated at around 60 million tons, including 5 million tons of fish.The severe decrease in that period of the population of the mussel Mytilus Galloprovincialis is also attributed to hypoxia (Zaitsev, 1997).In particular hypoxia affected a unique congregation of Phyllophora named "Zernov's Phyllophora field" after its first establishment in 1908 by S. A. Zernov (Zernov, 1909).It is not clear whether hypoxia affected directly the Phyllophora field or rather indirectly through the decline of filter feeders that may have decreased water transparency.In return, the decrease of Zernov's field constituted an additional positive feedback on hypoxia, as it has been identified as an important generator of oxygen at depths of 20-60 m (Zaitsev, 1997).Zernov's field is sustaining an ecosystem known as the Phyllophora biocenosis, offering habitat and spawning grounds to 47 species of fish, 110 species of invertebrates and other living organisms, among which 9 species are listed in the Red Book of Ukraine (Kostylev et al., 2010).Since recruits from this area disperse outward from the seaweed, the ecological importance of this field, and thus the damaging effect of hypoxia, extends beyond its actual boundaries.
The main aims of this study are (1) to understand and identify the drivers of hypoxia at seasonal and interannual scales, (2) to propose a statistical model to predict the interannual variability of hypoxia and to quantify the contributions of the identified drivers to this variability, (3) to assess the timescale of recovery of the system submitted to eutrophication due to the sediment inertia, and (4) to give some recommendations for the monitoring of hypoxia on the Black Sea NWS.
Since the numerical model presented here has the ambition to be used for the prediction of hypoxia in answer to management strategies, several error metrics, issued from the comparison of model results and in situ observations, are combined to illustrate the model performances regarding different aspects of the oxygen dynamics.
The paper is structured as follows.The 3-D coupled model is presented in Sect.2, along with the data and metrics used for the validation.Section 3 reviews the key processes driving the dynamics of oxygen in the Black Sea NWS and details the mechanisms which lead to the seasonal occurrence of hypoxia.Section 4 evaluates the model performance by comparing the variabilities of the oxygen distribution resolved by the model to those pictured from in situ data.Based on the previous description of the oxygen dynamics, Sect. 5 evaluates the interannual variability of hypoxia in response to eutrophication and meteorological external pressures.Section 6 discusses the practical implications of this study in relation notably to the implementation of the European Marine Strategy Framework Directive (Cardoso et al., 2010).Finally, a summary of the new findings provided by this study and some recommendations are given in Sect.7 2 Materials and method

The coupled model
The physical model is the GeoHydrodynamic and Environmental Research (GHER) 3-D hydrodynamical model described in, e.g., Beckers (1991) and Delhez (1996) and used in the particular case of the Black Sea in Grégoire et al. (1998Grégoire et al. ( , 2004)), Stanev andBeckers (1999), andBeckers et al. (2002).It resolves the currents and the physical state of the sea, which are temperature, salinity, internal mixing and surface elevation.The implementation used in the frame of these simulations is described in Capet et al. (2012), and validated in detail with satellite sea surface temperature (SST) and sea level elevation data.
The biogeochemical model simulates oxygen, nitrogen, phosphorus, silicate and carbon cycling and explicitly represents processes in anoxic and suboxic conditions.It is described in detail in Grégoire et al. (2008).Compared to Grégoire et al. (2008), the model has been coupled (online) in 3-D with the GHER 3-D hydrodynamic model and extended with a dynamic representation of the benthic compartment (in Grégoire et al. (2008), a reflective boundary condition was used for describing benthic degradation).This representation is based on the comparative analysis performed by Soetaert et al. (2000) on the representation of benthic processes that is both reliable and tractable for long-term simulations in a 3-D ocean model.Hence, the benthic compartment is described by a semi-empirical model whose state variables are 2-D variables (vertically integrated C and N content with two degrees of liability).From these 2-D stock variables, the fluxes of solutes are estimated using bulk parameters (i.e., the fraction of organic matter degraded by denitrification and anaerobic oxidation, oxygen consumption in nitrification).As in Soetaert et al. (2000), these bulk parameters are estimated from Monte Carlo simulations performed with an existing 1-D diagenetic model of the Black Sea sediments developed and validated by Wijsman et al. (1999).The organic matter accumulated in the sediment stocks is subjected to resuspension from the action of wave-and current-induced bottom stress, which is modeled using the setting presented in Stanev and Kandilarov (2012).

A. Capet et al.: Seasonal hypoxia in the Black Sea northwestern shelf
The coupling between the three components of the model (hydrodynamic, biogeochemical, benthic) is fully online.For instance, the transport of biogeochemical variables is determined at each time step from the advection and diffusion fields computed by the hydrodynamical model.Additionally, the rate of each biogeochemical process depends on the temperature field, and the oxygen solubility is also dependent on the salinity.On the other hand, the penetration of solar radiation, which governs the heating of the upper layer, is influenced by the vertical distribution of phytoplankton and organic matter computed by the biogeochemical model.Similarly the benthic component receives from the overlaying waters the export of detritus and diatoms and provides the flux of solutes issued from sediment degradation.The dynamics of benthic detritus depends on environmental conditions of bottom waters (temperature, currents including waves, oxygen, nitrate and ammonium).The coupling between the benthic and pelagic compartment is described in details in a companion paper in which the modeled benthic fluxes of oxygen, nitrate, ammonium, silicate and phosphate are compared with available observations.
The model domain covers the whole Black Sea, on a 15 km Arakawa-C grid, with 31 vertical levels using double sigma coordinates.Riverine water flows and nutrients' (NO x , PO 4 , SiO 2 ) monthly loads are issued from Ludwig et al. (2009) and the Black Sea commission data (after 2000).The total Black Sea riverine discharge is delivered through the 6 main rivers: Danube, Dniester, Dnieper, Rioni, Kizilirmak and Sakarya (the first three flowing directly on the NWS, Fig. 1).Averaged and constant atmospheric deposits of nitrate and phosphate are imposed all over the basin (NO x : 0.78 g m −2 yr −1 ; PO 4 : 0.7 mg m −2 yr −1 ; Kanakidou et al. (2012)), which represents 3 Gmol N yr −1 for the NWS (Fig. 1).Average concentrations of nitrogen organic forms are estimated for the three shelf rivers (Cauwet et al., 2002;Reschke et al., 2002;Walling and Fang, 2003).When multiplied by the considered discharges, this leads to a total load of 10.1 Gmol N yr −1 (7, 0.8 and 2.3 for Danube, Dnieper and Dniester, respectively).
Atmospheric forcings are issued from the ERA-Interim reanalysis provided by the ECMWF data center 1 with a 6 h and 0.75 • resolution.

Data sets
A total of 14 123 oxygen measurements, available in the NOAA World Ocean Database 2 (WOD), are retrieved for the area of investigation (Fig. 1) for 1980 to 2000.Additional data reported to the Black Sea Commission (BSC) are included in the data set, making it possible to cover the years 2000-2010 (636 points).Since hypoxia on the NWS usually occurs on the bottom of stratified areas at the end of summer, the assessment of the occurrence of hypoxia requires data collected below the mixed layer, during the summer months (August-October).Additionally, to quantify the intensity of the hypoxic event in term of spatial extension, a sufficient spatial coverage over the bathymetric gradient is required.While a satisfactory part of the WOD data meets these criteria (Table 1), a major proportion of BSC data are located in very coastal areas, and are usually not available for the period of hypoxia development.This means that basically no observation is available to assess the potential recovery of the Black Sea NWS from hypoxia after 2000.The data located in very coastal areas (i.e., bathymetry, Z, < 17 m; see Table 1) are not considered for the validation procedure because of the higher relative importance of processes unresolved by the model (e.g., small-scale hydrodynamics, benthic flora and fauna)

Model skill metrics
The data described in Sect.2.2 (Table 1) are used to assess the model ability to resolve the dynamics of oxygen using a point-to-point approach: each observation is compared with the model predictions at the data point (model results which are weekly averages are interpolated in time and space to obtain a value corresponding to the position and time of sampling).This level of validation is penalizing, as it requires perfect matches in space and time, but avoids the problem of representativity occurring when deriving volume averages from two samples with contrasted distribution, i.e., the sparse and strongly irregular distribution of observations, and the dense and regular distribution of model values (see also the discussion in Sect.5.6).
The following metrics are combined to address distinct aspects of the comparison between model and observations, according to the recommendations of Allen et al. (2007) and Stow et al. (2009).Besides, these can be compared with wellaccepted reference values, allowing the performances of the model to be judged in a more objective way.The mean and standard deviation of a variable X are noted X and σ X , respectively.
-The percentage bias compares the mean value of observations (O) and predictions (P ), and expresses the average bias as a percentage value of the observations average, P B = P −O O • 100.It indicates a global overestimation (underestimation) when positive (negative) and the model is considered as excellent (in that regard) when |P B| < 10 % (Allen et al., 2007).
-The standard deviation ratio, r σ = σ O σ P , compares the distribution of observations and predictions around their respective average.It indicates a larger variability in the observation (the predictions) when r σ > 1 (r σ < 1).The water column in this region is permanently stratified.The potential energy anomaly, i.e., the volume-specific difference in potential energy between the stratified and mixed state of the water column, is defined by where g is the gravity constant, Z the total water depth and ρ the density computed for the given value of temperature, T , and salinity, S. The thermal (φ T ) and haline (φ S ) contributions to the stratification can be computed by Haline stratification is present during the whole year and peaks in May, when river discharges are maximum, while the warming of surface waters enhances the stratification from the beginning of spring (early April) until mid-fall (mid-October) (Fig. 3).
The duration of this enhanced stratification period varies from year to year and is a key factor in determining the severity of bottom hypoxia.For further use, we define the stratification and mixing periods (Fig. 3) by using a criterion on the maximum value of the Brunt-Vaïasala frequency along the vertical: max Z z=0 −g δρ(z) ρ(z) ≥ and < 0.05 s −1 for stratification and mixing period, respectively.The mixed layer depth, z ρ , is defined as the depth where a density difference, ρ, compared to the 3m density reaches a threshold of 0.0125 kg m −3 as proposed for the Black Sea by Kara et al. (2009).z ρ reaches ∼ 10 m on average during the stratification period.
During the stratification period, bottom waters are virtually separated from surface waters.This clearly appears when looking at the differences between surface and bottom temperature (Fig. 4).The low-rising curve of bottom temperature indicates a limited warming of the bottom layer by the surface layer, which has the consequence that the average summer bottom temperature is determined by the winter/spring SST rather than by the summer SST.This is important since summer bottom temperature acts on the metabolic and respiration rates below the mixed layer depth and on the benthic diagenesis, hence directly on the intensity of the consumption of the oxygen pool locked below the thermocline.Obviously, the absence of tidal currents over the bottom also leads to less mixing.

Oxygen solubility and surface fluxes
The surface oxygen concentration is constantly regulated by the air-sea fluxes (Fig. 4), which tend to bring the surface concentration towards the saturation value.This saturation value depends on the oxygen solubility, which is mainly determined by the SST.Therefore, the seasonal SST signal explains the larger part of the surface concentration variability: the seasonal variations of the spatially averaged difference, [O 2 ] surf − [O 2 ] surf,sat , is less than half of that of [O 2 ] surf (σ = 13 and 30 mmol m −3 , respectively).
In spring and summer, surface waters are oversaturated in oxygen due to enhanced photosynthesis and to the decreased solubility of oxygen as SST increases.This leads to a release of oxygen towards the atmosphere (Fig. 4).In fall, the oxygen saturation increases as SST decreases, and the oxygen content of surface waters decreases as a result of the vertical mixing with the underlying oxygen-depleted waters.Surface waters are thus undersaturated and oxygen is taken up from the atmosphere.This ventilation of surface waters lasts until the next spring.
The exchanges of oxygen with the atmosphere of the investigated area are almost balanced over the year: the total annual uptake of 372 Gmol yr −1 is compensated by a release of 407 Gmol yr −1 , resulting in an annual net release of 35 Gmol yr −1 (Fig. 5).This annual net surface flux is variable and averages to an uptake of 2 Gmol yr −1 when the average is computed over the 1981-2009 period.While atmospheric fluxes have a marginal contribution to the annual oxygen budget, these play a crucial role in the seasonal mechanism of hypoxia as the intensity of the winter/spring ventilation fixes the initial concentration of oxygen in the bottom waters before their occlusion by the thermocline.

Biogeochemical processes and annual oxygen budget
Figure 5 presents the contribution of biogeochemical processes (budgets integrated over the whole shelf) to the simulated oxygen budget for the mixed and stratified period (distinguishing in this last case the layers above and below z ρ for the vertical integration).
Of the amount of oxygen produced by the (net) photosynthesis of autotrophs, 85 % is consumed for the remineralization of organic matter (in the water column and sediments), 8 % is consumed for the oxidation of ammonium and of the reduced substances released by the benthic layer, 4 % is exported laterally and 3 % is delivered to the atmosphere.About one-third of the organic matter degradation occurs in the sediments, but this number may vary according to the equilibrium status of the sediments compared to the nutrient load, as detailed in Sect.5.3.
Both processes are thus intensified in summer, but more markedly for the respiration, which shifts the summer balance towards oxygen consumption.
More importantly, at that time, we simulate different vertical distribution of autotrophs and heterotrophs: photosynthesis mainly occurs in the surface mixed layer (74 % of summer production), where light is available, while respiration of the sinking organic matter occurs principally below the mixed layer depth (66 % including benthic diagenesis).Additionally, 80 % of the summer chemical consumption occurs below the mixed layer depth.
These different distributions explain the summer depletion of oxygen below the mixed layer.
Finally, the concentration of oxygen in the bottom waters at the end of the stratification period is determined by the content of oxygen before stratification, which is fixed by winter ventilation (and hence the SST during this period), by the amount of oxygen consumed during the stratified period essentially through the heterotrophs respiration in the water column and the sediments, and by the duration of the stratified period.The importance of these factors is confirmed by the statistical analysis of the interannual variability described in Sect. 5.

Model validation
Since the present study is centered on the understanding and predictions of hypoxic events, the validation exercise is concentrated on the assessment of model performances in simulating the oxygen dynamics.The validation of the simulated physical variables is described in Capet et al. (2012) where the interannual variability  of the model- simulated SST and sea surface elevation is compared with satellite observations.This validation exercise demonstrates the model ability to correctly reproduce the interannual variability of SST and associated hydrodynamical fields such as the dynamics of the mixed layer.This means that the model has the ability to simulate the dynamics of water ventilation and resulting oxygenation, which is an important factor governing the oxygen dynamics (see Sect. 3).The simulated benthic fluxes of oxygen (as well as nitrate, ammonium, silicate and phosphate) have recently been compared (not shown) with all the data sets available on the Black Sea NWS (Friedl et al., 1998;Friedrich et al., 2002) and unpublished data kindly provided by I. Akoumianaki (EU Sesame project data set) and A. Lichtschlag (EU HYPOX project data set).This means that a good representation of the oxygen distribution and its variations in time and space by the model suggests a good representation of the pelagic sources and sinks of oxygen (e.g., photosynthesis, respiration, chemical reactions).
In the following, we first present the global model skills in simulating oxygen.Then, the validation procedure is focused on answering the following questions: 1. Is the model able to represent the seasonal cycle of oxygen?
2. Is the model able to represent the interannual variability of oxygen?
3. Is the model able to represent the location of oxygendepleted waters?
4. Is the model able to resolve the threshold of hypoxia?
The skill metrics are given in Table 2.

Global assessment of model performances
The global (i.e., considering all the data) percentage bias is virtually nul (0.03 %), indicating that the model resolves the oxygen dynamics without a persistent bias.The standard deviation ratio is higher than one, indicating a higher variability in the observations.This is expected since model results correspond to weekly averages, estimated on a 15 km × 15 km grid, while observed values are precisely located in space and time, and are therefore subject to small-scale variabilities which are either unresolved by the model (e.g., smallscale physical structure) or smoothed by the weekly averaging (e.g., the diurnal cycle of surface photosynthesis).The global N-S statistic is penalized by the large errors (Fig. 7c) associated with point-to-point validation, but still qualifies the model as "good" (Sect.2.3).

Is the model able to represent the seasonal cycle of oxygen?
Figure 6 shows the monthly distribution (median, percentiles) of the oxygen concentration deduced from observations and model results as well as the associated distribution of bias and percentage of hypoxic records.Model skill statistics are then computed on the medians obtained for each month (Table 2).These indicate excellent skill of the model regarding the timing and the range of seasonal variability.However, a seasonal signal can be seen in the monthly biases (Fig. 6c), which is very close to that of the oxygen saturation in surface waters (Fig. 4).The overestimation of the surface oxygen concentration in spring and the underestimation in fall occur when the surface waters are oversaturated and undersaturated, respectively.This suggests an underestimation of the air-sea fluxes intensity, possibly due to the use of averaged wind speed (hence of lower magnitude) to compute the exchange coefficients appearing in the formulation of the air-sea flux (Ho et al., 2006).This also explains the slightly larger seasonal variability in the model reflected by the standard deviation ratio.

Is the model able to represent the interannual variability of oxygen content?
Following the same approach, observations and predictions are gathered by years (Fig. 7), discarding the years with less than 10 available observations.The error metrics reveal an excellent ability of the model to represent the interannual variability present in the WOD data, but this ability is considered as good when including the BSC data in the computation.It should be noted, however, that the very sparse distribution of data after 2000 prevents a reliable estimate of the model skill during this period.
The standard deviation ratio indicates a slightly wider range of interannual variability in the observations, reflecting that some of the complexity of the oxygen dynamics (Sect.1) is not represented by the model.

Is the model able to represent the location of oxygen depleted waters?
The distribution of observed and simulated vertical oxygen profiles along the bathymetric axis (Figs. 8 a and b respectively) indicates a zone of low oxygen content along the bottom between approximately the 15 m and the 45 m isobaths in both cases.In waters deeper than 45 m, hypoxic conditions still appear at some sites in the data but not anymore in the model, while the model predicts hypoxia in the very coastal area (total water depth, Z, < 20 m) contrary to the data.The disagreement between model and data along the coast may be explained by several factors, such as the misrepresentation of coastal small-scale hydrodynamic features due to the insufficient horizontal resolution of the model and of the atmospheric forcings, and also the ignorance in the model of the role of the benthic flora (e.g., the seagrass Zostera) as oxygenator of the water.
A first mapping of the area of hypoxia over the NWS was proposed by Zaitsev (1997) (redrawn in Fig. 9a), who located bottom hypoxic waters near the rivers mouths, mainly north of the St. George Danube entrance, extending until the Ukrainian coast in some cases and, to a smaller extent, along www.biogeosciences.net/10/3943/2013/Biogeosciences, 10, 3943-3962, 2013 the Romanian coast.Unfortunately Zaitsev (1997) did not mention how this mapping was derived (e.g., oxygen concentration, Phyllophora extent), but we can note that the horizontal extent of hypoxia predicted by the model based on an oxygen threshold agrees quite well with that estimated from Zaitsev (1997) as shown by the comparison of Fig. 9a and c.The reduction of the southeastern extension of the area of hypoxia in the model compared to the mapping proposed by Zaitsev (1997) can be explained by the absence of distinction of the different branches of the Danube in the model.In the model, all the Danube discharges occur through the northernmost Chilia branch, which is the most important one (Fig. 1).It is noteworthy, however, that the model predicts hypoxia in the bay of Constanta in agreement with Zaitsev (1997) without the direct inputs of the southernmost St. George branch.

Is the model able to represent the threshold of hypoxia?
To address this question the number of hypoxic records ([O 2 ] < 62 mmol m −3 ) are compared along their seasonal (Fig. 6a) and interannual (Fig. 7a) distributions.It should be noted that only a few observations are below that threshold (174 measurements in total, none in the BSC database).
Meanwhile, the error statistics reported in Table 2 indicate that the distribution of hypoxic records among the different months and years is very well reproduced by the model.

An index to quantify the annual intensity of hypoxia
The repercussion of hypoxia in a given region depends not only on the spatial extension of the hypoxic event but also on its duration in time, which is a critical factor for the mortality of living organisms (Vaquer-Sunyer and Duarte, 2008).An index of the severity of hypoxia should thus merge these two characteristics of spatial extension and time duration.This information is rarely available from local observations since they usually provide either a snapshot view of the status of hypoxia (e.g., field cruise) or a local assessment (e.g., time series obtained at a fixed station).In some cases, spatial maps of areas affected by low oxygen concentrations are obtained by merging non-synoptic data sets.
In the present case, the model resolves both the spatial and temporal aspects of hypoxia, which allows appraising the temporal variability of the simulated benthic hypoxia area.Figure 11a presents examples of this seasonal variations for selected years, illustrating the limitation of defining hypoxia by a snapshot map.
In order to take into account temporal dynamics, we quantify the annual intensity of hypoxia with the index H defined by in which A(t) is the area affected by bottom hypoxia at time t, D is the timescale of the hypoxic event, and D its average computed over the simulation period (102 days).D corresponds roughly to the time during which hypoxia occurs over a surface larger than half of the maximal area (see Fig. 11a).
The normalization of the integral by D gives to H the units of an area.Eqs. ( 4) and ( 5) imply that H is equivalent to the maximal hypoxic area, max A(t), if D = D for this year, and higher (lower) if D is larger (smaller) than D, i.e., if the timescale of the hypoxic event is larger (shorter) than its average value.These formulations provide an easy algebraic relationship between H , D and max A(t), graphically illustrated in Fig. 11b.This figure highlights that H integrates both temporal and spatial aspects in its definition since high values of H are obtained for either large max A(t) or for large D.

Empiric linear models
Several studies use empirical models (e.g., linear or multiple linear regressions) to link the interannual variability of hypoxia to a selection of external drivers (Conley et al., 2007(Conley et al., , 2009;;Turner et al., 2012).These studies used observations and selected the potential drivers among available data, according to the current understanding of hypoxia.
In our case the process of hypoxia is completely and explicitly resolved by the model, and all potential drivers and processes may be extracted from the model results and forcings.Using this ideal "data set", we propose deriving an empirical relationship linking the intensity of hypoxia, represented by the H index (Sect.5.1), to its drivers.This empirical relationship is not based on a mechanistic description of the processes but has the following advantages: 1.It requires a reduced number of data for the prediction of hypoxia.
2. It is more easily comparable to studies addressing other areas, for instance by comparing the regression parameters of similar drivers among sites.This would allow deriving a relationship between these parameters and the characteristics (e.g., geomorphology) of the sites, extending the scope of hypoxia understanding across the global coastal areas.
3. It provides an easy way to identify the drivers and to assess their respective weight in governing the variability of hypoxia, differentiating for instance eutrophication from atmospheric impact.
However, even if the drivers involved in this regression analysis are assumed to have a causal link with hypoxia, we should be careful when using this formulation for making predictions in a range of conditions for which the model has not been fitted.This is an inherent limitation of statistical models.
A multiple linear regression procedure (e.g., Legendre and Legendre, 1998) consists in deriving a regression equation expressing the response variable Y as a linear combination of a set of predictors (or explanatory variables) X 1 , . . ., X N :  (1981-1991: light gray; 1992-2000: mid-gray; and 2001-2009: dark gray).Their distribution shows the variability between the years and the lack of relationship between max A(t) and D. In order to compare the case of the Black Sea NWS with other sites affected by hypoxia, the percentages given below the colored scale correspond to a scaling of H by the whole NWS area, depicted in Fig. 1.
where Y is the estimated value of Y and the coefficients a 1 , . . ., a N are computed automatically in order to minimize the residuals (Y − Y ) 2 .The predictors may be selected arbitrarily based on a priori assumptions, or selected automatically using statistic metrics.Here, we use a stepwise regression procedure (backward elimination based on p = 0.01 significance threshold) to identify among a preselected set of potential predictors issued from the model (i.e., winter SST; mean summer SST; late summer SST (Aug-Sept-Oct); March SST; mean winter organic carbon stock in the sediments; average summer winds: magnitude, zonal and meridional components; minimum summer average Brunt-Vaïasala frequency; water and nutrients riverine discharge) an optimal subset of predictors such that no additional predictor would increase significantly the coefficient of determination of the regression, R (computed as the ratio of the variance of the response explained by the predictors and the true variance of Y ) (Legendre and Legendre, 1998).
The stepwise procedure indicates that H can be retrieved from a linear combination of the annual nitrate load (N), the winter sediment organic content (C), March SST (T s ), and late summer SST (T f , averaged for Aug-Sept-Oct) (Fig. 12).In its standardized forms, the regression equation is writes as follows: with X * = X−µ X σ X .This regression equation explains 82 % (R 2 = 0.82) of the variability of H .The use of standardized predictors allows comparing the standard regression coefficients, a i , to measure the relative impacts of the different predictors on the variability of H .However, the confidence intervals associated with these coefficients (shown in Fig. 12) must be considered to ascertain within satisfying probabilities the preponderant role of a variable over another.
Climatic and eutrophication impacts on hypoxia may be distinguished by performing partial regressions, i.e., estimating how much of the variation of hypoxia can be attributed exclusively to climate (eutrophic) drivers once the effect of eutrophication (climate) has been taken into account.This analysis indicates that both factors have a similar impact, with eutrophication-related variables (N, C) explaining slightly more (R 2 = 0.36) than the climatic related predictors (T s , T f , R 2 = 0.27).A higher part of the H variance is attributed to climatic drivers when reducing the number of years considered for the partial regression.This indicates that climatic drivers mostly drive a shorter timescale variability while eutrophication drivers act on a longer term.Interestingly, the role of temperature introduces an oscillation in the interannual variations of hypoxia, which is related to the East Asia/West Russia atmospheric teleconnection pattern (Capet et al., 2012).

Interpretation of the selected predictors
A proper interpretation of regression coefficients requires that the selected predictors are not significantly correlated (Legendre and Legendre, 1998).Here, the predictors selected by the stepwise regression are each related to different processes, and their relation to hypoxia may be interpreted as follows.
Nitrate load is the main fuel of eutrophication.Model results show that N is significantly correlated to the maximum surface oxygen concentration simulated during the spring bloom (ρ = 0.7) and to the range of bottom oxygen decrease from spring to autumn (ρ = 0.86).
March SST, T s , is a better predictor of bottom hypoxia than the summer or winter SST.This important role of the SST in March has already been discussed in Sect. 3 and is due to two main reasons: 1.The March SST fixes the temperature of bottom waters for the duration of the stratification period (correlation between March SST and bottom summer temperature is ρ = 0.82) and hence affects the summer rates of organic matter respiration.
2. The March SST controls oxygen solubility and hence the intensity of oxygenation before the stratification period.Therefore, it fixes the oxygen content of bottom waters in spring (ρ = 0.85), before their occlusion by the thermocline.
The late summer SST, T f , plays an important role in determining the duration of the stratified period.As the lowest oxygen concentrations are reached at the end of summer, an anomalous extension of the stratification period causes increased damages to the ecosystem.In particular, T f is related to the annual release of reduced compounds from the sediments (ρ = 0.67), including H 2 S (UkrSCES, 2002).
While a part of the sediment organic matter is labile and remineralized within a timescale of months, a part of the sediment organic stock is more refractory (semi-labile) and requires more time for its remineralization.During years of high nutrient loads, the semi-labile stock accumulates year after year and continues to cause oxygen consumption several years after high-riverine-load events.The predictor C, the amount of organic matter in the sediments computed for the winter season, reflects this semi-labile component, as the labile fraction of organic matter is almost entirely remineralized in summer or resuspended at the begin of winter.The mean summer benthic oxygen demand is significantly correlated to C (ρ = 0.82), indicating the important role of this accumulated stock in benthic diagenesis.

Impact of sediment inertia on the timescale of recovery
The fact that the accumulation of organic matter in the sediments keeps record of the nutrient load history and may significantly intensify hypoxia is an important mechanism which has to be considered when addressing the long-term variability of hypoxia in response to a given regime of nutrient load.In the case of decreasing load (which has been the case in study area since the early 1990s) the sediment introduces an inertia in the process of recovery.In the case of a maintained or increasing level of nutrient load, this mechanism causes an increase of the response ratio "hypoxic area / annual nutrient load", as illustrated in the Gulf of Mexico by Turner et al. (2008) The evolution of C as a function of the annual nitrate load history can be reproduced by a parametric model using, in addition to N, the March SST (as a proxy of bottom temperature, Sect.3) and the history of the silicate-to-nitrate ratio (i.e., Si : N) of riverine load: in which y represents any particular year.β 0 can be seen as a bulk parameter representing the natural decrease rate of C due to resuspension, diagenesis and burial (i.e., for N = 0 and T * =0 we have C(y) = C(0) exp −β 0 y ); β has the same interpretation as β 0 but additionally considers the effect of Biogeosciences, 10, 3943-3962, 2013 www.biogeosciences.net/10/3943/2013/temperature through a Q 10 formulation; and α represents the part of riverine load that feeds C, and depends on the ratio Si : N.This ratio is a proxy of the proportion of diatoms in phytoplankton which conditions the amount of organic matter sinking to the sediment (this effect has been found important to represent the evolution of C after 2002 when the Si : N ratio of riverine nutrients increases substantially and modifies the composition of phytoplankton communities).
Using the interannual evolution of C simulated by the 3-D coupled model, parameters α 0 , β 0 , Q 10 and α Si : N are estimated using a non-linear fitting procedure.The reconstruction of C history from the initial condition, C 0 , using Eq. ( 8) and the fitted parameters provides a coefficient of determination R 2 = 0.95 (Fig. 13).The timescale of sediment memory after a high-load period can be estimated from the parameter β 0 which reaches a value of 0.16(±0.02,p > 0.1), corresponding to a characteristic timescale of 1/β = 9.3 yr.
Important information that can be derived from Eq. ( 8) is the value of the stock of semi-labile benthic carbon reached at equilibrium (noted C eq ) for a given value of N (noted N eq ).To derive an estimation for the values C eq (N eq ), Eq. ( 8) is integrated for 10 yr using a value of N that increases linearly from its value reached in 2009 to the selected N eq (imposing typical temperature and Si conditions).Then, the integration in continued for 10 additional years while maintaining a constant load level (N = N eq ).Fig. 13 shows the levels of C eq reached for different values of N eq .It can be seen that in some cases 10 yr of integration is not long enough to reach steady state in the sediments.

Assimilation capacity
The main motivation for studying hypoxia is to derive management strategies (e.g., river loading restriction policies ) to avoid it, or at least to keep hypoxia intensity below a threshold, H 0 , that ensures healthy functioning of the coastal and marine ecosystem, as defined by the Marine Strategy Framework Directive (MSFD) established by the European Union (Cardoso et al., 2010).
As N is the only driver identified in Sect.5.2 that can eventually be regulated, we propose focusing on the dependence of H as a function of N, by taking a closer look at the distribution H (N ).
Obviously, the original time series H does not only depend on N, but is also affected by the specific effects of unbalanced sediment and climatic anomalies.Using the time series of these drivers and the coefficients obtained in Sect.5.2, we can rewrite Eq. 7 to obtain H eq , the intensity of hypoxia that correspond to the particular situation of average atmospheric conditions (i.e., climatic drivers equal to their average for 1981-2009) and sediment stock in equilibrium with the river nitrogen load.
in which the values C eq corresponding to various values of N are estimated from Sect.5.3 (Fig. 13).H eq is now only depending on N and may thus by used to assess the value of H reached for a given nitrogen discharge (N eq ) under climate conditions typical of the last 30 yr and with a sediment in equilibrium with N eq .The resulting distribution, H eq (N eq ), can be described by a power law (H eq = a • (N eq − N O ) b ; R 2 = 0.90; Fig. 14).Unlike the linear relationship, implicitly assumed within Sect.5.2, the use of a power law, in addition to resulting in slightly enhanced R 2 , also presents the conceptual advantage of allowing H eq (N 0 ) = 0 for N 0 > 0, whether a linear relationship would result in H eq (N = 0) > 0, which makes no sense.
Atmospheric drivers have an impact on H which compares to that of eutrophication (Sect.5.2).This is shown in Fig. 14 by indicating the effect on H eq of one standard deviation of the atmospheric drivers around their average value (i.e., H eq ± [a T s + a T f ] • σ H ).This important effect of the climatic drivers underlines that, even in the framework of pertinent nutrient management, the system is not free of risk of hypoxia occurrence due to atmospheric variability.

Practical relationship
We have shown and explained the role of the different predictors appearing in the formulation of H , and in particular the importance of the accumulation of organic matter in the sediments.However, as the predictor C can not be easily estimated without a model, we propose for practical purposes a relationship that directly links the value of H , corresponding to equilibrated organic sediment content, to accessible predictors: , i.e., H corrected from the effect of unbalanced sediment contents.The confidence intervals are computed for p < 0.1.

The recovery process: representativity of monitoring coverage
There are different and contradictory ideas about the recovery process from hypoxia on the Black Sea NWS after eutrophication in the 1990s and 2000s.For instance, the study of Mee (2006), which is cited in number of reviews (e.g., Diaz, 2001;Diaz and Rosenberg, 2008;Kemp et al., 2009;Rabalais et al., 2010), finds a clear recovery during the 1992-2000 period, closely following the reduction of nitrate load.
A quite different picture of the trend of hypoxia during the last decades is presented in the Ukrainian national report (UkrSCES, 2002).More specifically, Fig. 9 8).After 2009, the parametric model is integrated for 10 yr using a value of N that linearly increases from its value reached in 2009 to different levels of N eq .Then, the integration is continued for 10 additional years, maintaining constant the level of nitrate load (N = N eq ).
Fig. 14.Level of hypoxia, H eq , reached at "equilibrium" for a range of nitrate riverine load, N eq .By equilibrium we mean that the atmospheric variables remain at averaged conditions  and that the sediment is in equilibrium with the nitrate load (Sect.5.3).
The colored dots indicate the distribution of H eq obtained from the 3-D model results (Eq. 11).This distribution can be described by a non-linear relationship between H eq and N eq (R 2 = 0.90).The effect of climatic variability is illustrated by adding (subtracting) the effect of one standard deviation for the two climatic predictors to the central line of H eq .the interannual variability of the extent of hypoxic areas provided by these two studies (and additionally by the model presented in this paper) and evidences two main important differences: 1.The maximum surface is 2 times smaller in UkrSCES (2002) (20 000 km 2 , in 1983).
We can note that the model estimations are closer to the Ukrainian report concerning the recovery of the system from hypoxia during the 1990s.These different interpretations of the same phenomenon may come from the use of different data sets by the authors.Unfortunately, we were not able to find details about the data and the methods used to estimate the extent of hypoxia in the study of Mee (2006).We only find that the presented values before 1993 match those presented by Zaitsev (1993), in which no more information is given concerning data coverage and computation methods.
To illustrate how sensitive the conclusions concerning hypoxia and its recovery are to the distribution of data, Fig. 10 compares the interannual evolution of the summer (July-August) oxygen concentration as obtained from the available data, from the model outputs restricted to data locations and from the model outputs averaged over the whole area (Fig. 1).
The months July and August were chosen because 1. those are the months with the highest number of observations, with data for almost every year (Fig. 6), 2. the model bias is minimal for these months (Fig. 6c), 3. and bottom hypoxic events are frequent.
Observations and model results obtained at observation locations both indicate a rise of the average summer oxygen concentrations during the 1990s.However, this trend is not visible anymore when the model results are considered over the whole area, which indicates that a part of this apparent recovery may be non-factual and rather results from the inhomogeneity of monitoring after the mid-1990s.
Similarly, Fig. 7 shows that there is a significant correlation (ρ = 0.63, p < 0.01) between the total number of observations performed during a given year and the percentage of hypoxic events revealed by these data.As an example, no hypoxic event is recorded between 1995 and 2000, but a careful look at the data distribution during this period reveals that no measurement has been made north of 45 • 30 where hypoxia usually occurs.This suggests that the observation coverage does not always allow appraising the "real" interannual variability of hypoxia and biased the interpretation of its potential recovery.
Rather than the consensual picture of a recovery from hypoxia closely following the decrease in nutrient load in the Black Sea NWS, the present study indicates a delayed recovery, attributed to the sediment's memory of eutrophication events (Sect.5.2).This slower recovery may explain why biological signs of recovery only started to be visible in the 2000 decade (Mee et al., 2005;Langmead et al., 2009;Kostylev et al., 2010).Furthermore, the present-day occurrence of hypoxia is supported directly by in situ observations (Friedrich et al. (2002); J. Friedrich, personal communication, 2010;M. E. Mihailov, personnal communication, 2010) and indirectly by the occurrences of massive fish mortality (D. Banaru, personnal communication, 2012).

Discussion
The present study shows that over the whole period  the impact of atmospheric drivers on hypoxia is commensurate with that of eutrophication (Sect.5.2).It is interesting to remark that, while the eutrophication predictors generally have lower values in the last decade, the climatic predictors have larger values (not shown) and are thus mainly responsible for the high values of H occurring in the last decade (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009).This stresses that a possible climate change will complicate the management of hypoxia through riverine load restriction policies.It must be noted that in Sect. 5 the effect of temperature is interpreted as a natural variability, which, while indicative, may not reflect the long-term effect of a firm warming trend on the system.It can be seen in Fig. 11 that the high H indexes in the last decade are mainly due to longer temporal extension of hypoxia, while high indexes in the eutrophication period were related to larger spatial extension.Because of these different temporal dynamics, the eutrophication-and climate-driven hypoxia may bear a distinct impact on the ecosystem.
Lethal and non-lethal effects of low oxygen concentrations on living organisms strongly depend on the tolerance of the considered species, of the low level of oxygen concentration and of the duration of the hypoxic event (Vaquer-Sunyer and Duarte, 2008).In particular, the damaging effect of hypoxia on a given ecosystem is not expected to vary linearly with the index H defined in Sect.5.1.H is rather defined as an index to quantify the importance of hypoxia as an environmental stressor.The choice of a sustainable level of H that meets the requirements of good environmental status (GES), as defined by the European Marine Strategy Directive (Cardoso et al., 2010), is therefore a very delicate issue which requires the combination of appropriate tools, as well as dedicated and site-specific studies.
The response of H to eutrophication and climatic variations given in Sect. 5 may not be valid when using values of the drivers (predictors) far out of the range covered by the present hindcast model experiment.For instance, we suspect that for low nutrient load conditions -not encountered during the period 1981-2009 and, hence, for which the model could not have been validated -the effect of benthic macro-biota could play an important role in regulating hypoxia.Nevertheless, Sect. 5 gives useful indications of the sensitivity of hypoxia to the riverine nitrate load.These can be used to construct scenarios of river load restrictions to be tested with further modeling experiments in order to give an estimation of the maximum nitrate load which would allow the Black Sea NWS to fully recover from hypoxia.To be realistic, these scenario model experiments should be forced with realistic climatic scenario.

Conclusions
The Black Sea NWS has been affected by seasonal bottom hypoxia for decades.In the worst cases, the hypoxic zone covers all of the shelf north of 45.5 • N and extends along the Romanian coast in the bay of Constanta.The oxygendepleted waters are located at the bottom between the coast and the 45 m isobaths.
These hypoxic events have a dramatic impact on the benthic biota and have been linked to the almost disappearance of the mussel Mytilus Galloprovincialis in the 1980s, the reduction of the surface of Zernov's Phyllophora field and the loss of important commercial fish stocks.Bottom hypoxia alters the health of the ecosystem and is a serious barrier for reaching GES in the Black Sea.However, its understanding, mapping and mitigation still require tools able to address its temporal and spatial scales of variability and to incorporate its various drivers.
In this paper, we use a combination of mechanistic modeling and statistical tools to address the complexity of hypoxia.A 3-D coupled hydrodynamical-biogeochemical model is used to represent the oxygen dynamics during the 1981-2009 period in response to the physical and biogeochemical processes submitted to the historical nutrient load and atmospheric variability.
The ability of the coupled model to simulate the dynamics of oxygen, at seasonal and interannual scales as well as its spatial variations, is asserted from the computation of various error metrics between model results and in situ measurements gathered from the WOD (18 735 points, 1980WOD (18 735 points, -2000) ) and from the BSC (733 points, 2000-2009).In particular, the comparison of the observed and simulated number of occurrence of hypoxic points ([O 2 ] < 62 mmol m −3 ) for each month and year indicates that the model is able to render this particular threshold.
The mapping of the surface of hypoxia from punctual measurements is a delicate issue that may lead to discrepancies between sources (e.g., Mee, 2006;UkrSCES, 2002)  on the status of hypoxia.For instance, a concentration of observations south of the typical hypoxia area leads in some cases to the conclusion that a recovery occurred after the early 1990s, shortly after the reduction of nutrient load.This quick recovery is not reflected by the model results, which rather depict a recovery delayed by the inertia of the sediment layer and mitigated by the important role of atmospheric conditions (see below).There is thus an urgent need for dedicated oxygen-monitoring efforts concentrated on the sites and periods of hypoxia occurrence as well as an integration of the collected data sets in an internationally accessible database.For instance, almost no observations are available during the periods and in the regions of seasonal hypoxia occurrence after 2000.
The coupled model has been used to understand the dynamics of hypoxia at seasonal scales and in particular its driving factors.Seasonal hypoxia occurs because stratification, resulting first from river discharge (April) and then from the warming of surface waters (May), isolates bottom waters from the surface, with the consequence that oxygen depletion caused by the respiration of sinking organic matter can not be balanced by turbulent re-oxygenation.Just before the onset of stratification, March SST significantly impacts the intensity of hypoxia by 1. controlling the oxygen solubility in the last period of vertical mixing, and hence the oxygen concentration of bottom waters before their occlusion, 2. and setting the temperature of bottom waters for the stratification period, and hence affecting the summer rates of organic matter respiration.
Pelagic respiration and benthic diagenesis cause the decrease of oxygen concentrations in the bottom waters until the end of the stratified period (October).While pelagic respiration is mostly related to nitrate load, benthic diagenesis of accumulated organic sediments introduces a memory of eutrophication with a typical timescale estimated at 9.3 yr, which contributes to maintain hypoxia after high-nutrientload events.
The lowest oxygen concentrations are reached at the end of the stratification period.An anomalously long duration of the stratification, due to anomalously high SST in late summer, decreases the bottom oxygen level to very low values and sustains the release of reduced substances to the overlaying waters such as hydrogen sulfide with potentially dramatic consequences for the benthic ecosystem.
The intensity of hypoxia, as an environmental stressor, is quantified for each year by a index H , defined to account both for the spatial and temporal extension of the hypoxic event, as both these aspects matter regarding the repercussions of hypoxia on the living organisms.
Using the results of the 3-D model and the external forcings as an ideal "data set", long time series  of physical and biogeochemical variables serve as predictors in a multiple regression model representing the interannual variations of H .
Within this large set, a stepwise regression procedure isolates four driving predictors: the annual load of nitrate, the semi-labile carbon content in the sediments, the SST in March and the SST in the late summer.These four predictors allow reconstructing 82 % of the hypoxia intensity interannual variability.Partial regressions show that the drivers linked to eutrophication (i.e., nitrate load, sediment organic content) explain a slightly larger part (36 %) of the interannual variability in hypoxia intensity than climate predictors (March SST, late summer SST, 27 %).
From this regression we derive a non-linear relationship that links H to the riverine nitrogen load for the particular conditions of sediments in equilibrium with the river nitrogen load and for a typical present climate.This relationship has the advantage of simplicity but should be used with care outside its conditions of derivation (e.g., a different climate, low levels of nitrate load not encountered during the period 1981-2009).
Finally, this relationship could be used to identify the level of the riverine nitrate load that is critical in terms of hypoxia generation and to link these with management strategies of river discharges that integrate socio-economic aspects.This is an important step towards reaching the GES of marine waters as defined by the European Marine Strategy Framework Directive (Cardoso et al., 2010).

Table 1 .-
Figure2illustrates the seasonal evolution of a vertical oxygen profile simulated at an arbitrary site affected by seasonal hypoxia (point A in Fig.1).The figure shows how stratification allows the settlement of bottom hypoxic conditions by preventing the bottom oxygen-depleted waters from mixing with the surface ventilated waters.

Fig. 1 .
Fig. 1.This map shows the area of investigation over which spatial averages are computed, the area in which the model predicts hypoxia and anoxia for water depth > 120-150 m (averaged situation for the period 1981-2009), and the entrances of the rivers considered in the model.The three distributaries of the Danube are also indicated (from north to south: Chilia, Sulina and Saint George).The dot (A) locates the profile presented in Fig. 2.

Fig. 2 .
Fig. 2. Typical seasonal evolution of the oxygen vertical profile in the northern part of the Black Sea NWS (point A in Fig. 1).Hypoxic conditions (i.e., oxygen concentrations below the threshold of 62 mmol m −3 ) are underlined by the color scale.The mixed layer depth is represented by the position of the black line and the intensity of the stratification, as estimated by the Brunt-Vaïasala frequency, by the line thickness.

Fig. 3 .
Fig. 3. Seasonal evolution of the vertical stratification as appraised by the potential energy anomaly (Eq.1), i.e., the volume-specific difference in potential energy between stratified and mixed state of the water column.The red and blue lines depict the thermal and haline contributions to the potential energy anomaly (Eqs. 2 and 3).Averages are computed over the investigation area presented in Fig. 1.Vertical lines indicate the beginning and the end of the stratification period defined by a criterion on the Brunt-Vaïasala frequency (Sect.4.1).

Fig. 4 .
Fig. 4. Upper panel: seasonal cycle of the surface and bottom oxygen concentration and of the surface saturation value determined by the sea surface temperature.During periods of oversaturation (i.e., the sea surface oxygen concentration is higher than the saturation value), the air-sea oxygen flux is directed to the atmosphere, while in the case of undersaturation this flux is directed to the sea.Lower panel: seasonal cycle of the surface and bottom temperature.

Fig. 5 .
Fig.5.Contribution of biogeochemical processes to the oxygen budget for the mixing and stratification period (left and right respectively).For the mixing period (a), processes are integrated over the whole water column, while during the stratified period the integration is performed separately above (b) and below (c) the mixed layer depth.The following processes are considered: net photosynthesis (Auto., including oxygen production associated with nitrate reduction); heterotroph respiration (Hetero.); the oxygen air-sea fluxes (Surf., counted positive when directed towards the sea); benthic diagenesis (Benth.);and oxygen consumption in chemical reactions (Chem., including nitrification and the oxidation of reduced substances).

Fig. 6 .
Fig. 6.The ability of the model to resolve the seasonal cycle of oxygen concentration is addressed by merging observations and predictions by month.For this comparison, model results are interpolated at the location and time of data sampling.(a) Proportion of records below the hypoxic threshold ([O 2 ]< 62 mmol m −3 ).(b) Distribution (median, percentile) by month of the oxygen concentrations: the central bars indicate the medians, the filled areas extend to the 0.25 and 0.75 percentile (respectively, Q 25 and Q 75 ), the dotted bars indicate the 1.5 • (Q 75 − Q 25 ) range, and the extreme values are marked by an "x".(c) Distributions of biases: positive (negative) bias correspond to overestimation (underestimation) by the model.

Fig. 7 .
Fig. 7.The ability of the model to resolve the interannual variability of oxygen concentration is addressed by merging observations and predictions by year.For this comparison, model results are interpolated at the location and time of data sampling.(a) Proportion of records below the hypoxic threshold ([O 2 ] < mmol m −3 ).(b) Distribution (median, percentile) by year of the oxygen concentrations: the central bars indicate the medians, the filled areas extend to the 0.25 and 0.75 percentile (respectively, Q 25 and Q 75 ), the dotted bars indicate the 1.5 • (Q 75 − Q 25 ) range, and the extreme values are marked by an "x".(c) Distributions of biases: positive (negative) bias correspond to overestimation (underestimation) by the model.

Fig. 8 .
Fig. 8. (a) Observed and (b) predicted vertical distribution of oxygen concentrations along the bathymetric axis.

Fig. 9 .Fig. 10 .
Fig. 9. (a) Area affected by hypoxia, redrawn from Zaitsev (1997), and locations of hypoxic records from the WOD database.(b) Spatial extension of the area over which the model predicts hypoxia lasting more than one month.Variability due to meteorological effects is filtered out by averaging the results over 6 yr periods.(c) Extension of the surface affected by bottom hypoxia, as reported in the literature(Mee, 2006; UkrSCES, 2002)  and simulated by the 3-D model.

Fig. 11 .
Fig. 11.(a)Seasonal dynamics of the hypoxic area, A(t), for selected years (colored solid lines) and climatological average (black solid line).The dotted lines indicate the time during which hypoxia covers more than half of the maximal area simulated during the considered year (i.e., the timescale of hypoxia, D, defined by Eq. 4).(b) The H index (colored lines) represents the annual intensity of hypoxia as an environmental stressor and combines the aspects of spatial and temporal extension of the hypoxic event.Its definition (Eq.4) allows this graphical representation of the relationship between H , max A(t), and D. The circles locate the hypoxic events simulated for each year(1981-1991: light gray; 1992-2000: mid-gray; and 2001-2009: dark gray).Their distribution shows the variability between the years and the lack of relationship between max A(t) and D. In order to compare the case of the Black Sea NWS with other sites affected by hypoxia, the percentages given below the colored scale correspond to a scaling of H by the whole NWS area, depicted in Fig.1.

Fig. 12 .
Fig. 12. Interannual variations of the annual hypoxic index, H (Eq. 4), diagnosed from the results of the 3-D coupled model and computed from the regression Eq. (7).The standardized regression coefficients on the right panel indicate the relative role of each predictor in explaining the variability of H , with error bars indicating the p = 0.5 and 0.9 confidence intervals.Partial regressions assign a comparable part of the variance of H to eutrophication (N, C) and climate (T , D) drivers.

Fig. 13 .
Fig. 13.Interannual variations of the sediment organic content in winter, C, simulated by the 3-D coupled model and by the parametric model described by Eq. (8).After 2009, the parametric model is integrated for 10 yr using a value of N that linearly increases from its value reached in 2009 to different levels of N eq .Then, the integration is continued for 10 additional years, maintaining constant the level of nitrate load (N = N eq ).
due to different spatial coverage of the data.The spatial heterogeneity of the data among years could lead to biased conclusions www

Table 2 .
Model skill statistics obtained when considering all the observation/model prediction pairs (Global), annual medians (Years, m) annual number of hypoxic events (Years, n H ), monthly medians (Months, m), monthly number of hypoxic events (Months, n H ). Error metrics are defined in Sect.2.3: ρ: Pearson correlation; N-S: Nash-Sutcliffe efficiency; and r σ : standard deviation ratio.The numbers in parentheses are obtained when only considering the WOD database (Sect.2.2).