References

Biogeosciences Discussions This discussion paper is/has been under review for the journal Biogeosciences (BG). Please refer to the corresponding final paper in BG if available. Abstract The Black Sea northwestern shelf (NWS) is a shallow eutrophic area in which seasonal stratification of the water column isolates bottom waters from the atmosphere and prevents ventilation to compensate for the large consumption of oxygen, due to respiration in the bottom waters and in the sediments. 5 A 3-D coupled physical biogeochemical model is used to investigate the dynamics of bottom hypoxia in the Black Sea NWS at different temporal scales from seasonal to interannual (1981–2009) and to differentiate the driving factors (climatic versus eutrophication) of hypoxic conditions in bottom waters. Model skills are evaluated by comparison with 14 500 in-situ oxygen measurements 10 available in the NOAA World Ocean Database and the Black Sea Commission data. The choice of skill metrics and data subselections orientate the validation procedure towards specific aspects of the oxygen dynamics, and prove the model's ability to resolve the seasonal cycle and interannual variability of oxygen concentration as well as the spatial location of the oxygen depleted waters and the specific threshold of hypoxia. 15 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, receiving large inputs of nutrients from the Danube, Dniestr and Dniepr rivers, and extends, during the years of severe hypoxia, towards the Romanian Bay of Constanta. In order to explain the interannual variability of bottom hypoxia and to disentangle 20 its drivers, a statistical model (multiple linear regression) is proposed using the long time series of model results as input variables. This statistical model gives a general relationship that links the intensity of hypoxia to eutrophication and climate related variables. The use of four predictors allows to reproduce 78 % of hypoxia interannual variability: 25 the annual nitrate discharge (N), the sea surface temperature in the month preceding stratification (T), the amount of semi-labile organic matter in the sediments (C) and the duration of the stratification (D). Eutrophication (N, C) and climate (T , D) predictors 18398 explain a similar amount of variability (∼ 35 %) when considered separately. A typical timescale of 9.3 yr is found to describe the inertia of sediments in the recovering process after eutrophication. From this analysis, we find that under standard conditions …


Conclusions References
Tables Figures

Back Close
Full explain a similar amount of variability (∼ 35 %) when considered separately.A typical timescale of 9.3 yr is found to describe the inertia of sediments in the recovering process after eutrophication.
From this analysis, we find that under standard conditions (i.e.average atmospheric conditions, sediments in equilibrium with river discharges), the intensity of hypoxia can be linked to the level of nitrate discharge through a non-linear equation (power law).
Bottom hypoxia does not affect the whole Black Sea NWS but rather exhibits an important spatial variability.This heterogeneous distribution, in addition to the seasonal fluctuations, complicates the monitoring of bottom hypoxia leading to contradictory conclusions when the interpretation is done from different sets of data.We find that it was the case after 1995 when the recovery process was overestimated due to the use of observations concentrated in areas and months not typically affected by hypoxia.This stresses the urging need of a dedicated monitoring effort in the NWS of the Black Sea focused on the areas and the period of the year 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 since 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, Black, Kattegat Sea, Gulf of Mexico, East China Sea), 39 % of which are located in Europe (Zhang et al., 2010).
Hypoxia is usually defined as an oxygen concentration below the 65 mmol m −3 (= 2 mg l −1 ) threshold.However, it has been shown that low oxygen concentrations already affect living organism 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 Introduction

Conclusions References
Tables Figures

Back Close
Full renew the oxygen consumed for the degradation of the organic matter.The coastal area which receives the input of organic matter from river discharges 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 it: directly affects living organisms (e.g.benthic organisms, fishes) (e.g.Vaquer-Sunyer and Duarte, 2008), may, in extreme cases, leads 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 pore water chemistry by entraining the death of bioturbators (e.g.Levin et al., 2009), 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 and geomorphological characteristics of the water body (Kemp et al., 2009).Among this diversity of hypoxic sites, the Black Sea north-western shelf (NWS) is to be classified with the permanent seasonal hypoxia in a stratified body, which according to Kemp et al. (2009) makes it comparable, among others, to the Chesapeake Bay (e.g.Boesch et al., 2001), the Pensacola Bay (Hagy III and Murrell, 2007), the Introduction

Conclusions References
Tables Figures

Back Close
Full While on a whole it is recognized that hypoxia in those systems is permitted by stratification and fueled by human-induced 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 discharge whose time scale is difficult to estimate, although essential to the study of the system recovery.Moreover, hypoxic conditions affect the sediments 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 discharge (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.recovering) are challenging to assess (Cardoso et al., 2010).
Sound management allowing the preservation and restoration of hypoxia-damaged ecosystems to a healthy state therefore needs a good understanding of these interconnected processes (Boesch et al., 2001).Introduction

Conclusions References
Tables Figures

Back Close
Full 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 to investigate 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 3-D physical-biogeochemical model is used to investigate the spacetime 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 north-western 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 that ∼ 100-150 m leading to anoxia of waters below that depth.Moreover, hypoxia on the NWS was first observed in the early 70's at 30 m depth in a region located between the Danube and Dniestr 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 the 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).The 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 Introduction

Conclusions References
Tables Figures

Back Close
Full  , 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 mains 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 time scale of recovery of the system submitted to eutrophication due to the sediment inertia, (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 aspect 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 gives the results of the validation, in answer to specific questions related to the oxygen dynamics.Section 4 reviews the key processes interacting in the settlement of seasonal hypoxia.Based on this understanding Sect. 5 evaluates the interannual variability of hypoxia in response to eutrophication and meteorological external pressures.A summary of the new findings provided by this study and some recommendations are given in Sect.7 in relation notably to the implementation of the European Marine Strategy Framework Directive (Cardoso et al., 2010).Introduction

Conclusions References
Tables Figures

Back Close
Full (GeoHydrodynamics and Environment Research laboratory) three-dimensional, nonlinear and baroclinic model which has been successfully applied to explore the general circulation of the Bering Sea (e.g.Deleersnijder and Nihoul, 1988), the North Sea (e.g.Delhez, 1996), the Mediterranean Sea (e.g.Beckers, 1991) and the Black Sea (Gr égoire et al., 1998;Beckers et al., 2002;Gr égoire and Friedrich, 2004;Vandenbulcke et al., 2009;Capet et al., 2012).It uses a refined turbulent closure scheme and solves the motions of the free surface.The biogeochemical model simulates oxygen, nitrogen, phosphorus, silicate and carbon cycling and explicitly represents processes in anoxic and suboxic conditions.It is described in details and validated in Gr égoire et al. (2008).In this paper, we will focus on the oxygen cycle.
The sediment layer is represented by horizontal variables.The model simulates resuspension processes as proposed by Stanev and Kandilarov (2012) and the variable benthic diagenesis is represented using a vertically integrated dynamic sediment model (based on Wijsman et al., 1999) where the bottom fluxes of dissolved constituents are parameterized based on mass budget considerations (Soetaert et al., 2000).
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 discharges are issued from Ludwig et al. (2009) and the Black Sea commission data.The total Black Sea riverine discharge is delivered trough the 6 main rivers locations: Danube, Dniestr, Dniepr, Rioni, Kizilmark and Sakarya (the three first 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 yr −1 for the NWS (Fig. 1).Nitrogen discharge in organic form (both dissolved and particulate) are estimated to 10.1 Gmol yr −1 (7, 0.8 and 2.3 for Danube, Dniepr and Dniestr, respectively) in average and imposed as fixed concentrations (Cauwet et al., 2002;Reschke et al., 2002;Walling and Fang, 2003).Introduction

Conclusions References
Tables Figures

Back Close
Full Atmospheric forcings are issued from the ERA-interim reanalysis provided by the ECMWF data center1 with a 6-h and 0.75 • resolution.

Datasets
14 123 oxygen measurements, available in the NOAA World Ocean Database2 (WOD) are taken in the area of investigation between 1980 and 2000.Additional data reported to the Black Sea Commission (BSC) are included in the data set, allowing 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 status of hypoxia requires data located at some distance from the coast, below the mixed layer and measured during the summer months (August-October).While a satisfactory part of the WOD data meets these criteria (Table 1), a major proportion of BSC data are located in very coastal area, and are usually not available for the period of hypoxia development.This means that basically no observations are available to assess the potential recovery of the Black Sea NWS from hypoxia after 2000.

Model skill metrics
The data described in Sect.2.2 (except those located in very coastal areas (i.e.depth < 17 m, see Table 1) are used to assess the model ability to resolve the oxygen dynamic using a point-to-point validation procedure which consists in comparing the observations 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 Introduction

Conclusions References
Tables Figures

Back Close
Full and strongly irregular repartition of observations, and the dense and regular repartition of model values (see also the discussion in Sect.5.5).
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); Stow et al. (2009).Besides, these can be compared with well accepted reference values allowing to judge in a more objective way the performances of the model.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 observa- -The standard deviation ratio, , compares the distribution of observations and predictions around their respective average.It indicates a larger variability in the observation (resp.the predictions) when r σ > 1 (resp.r σ < 1).

Model validation
We first present the global model skills in simulating oxygen.Then, the validation procedure is focused to answer 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 oxygen depleted 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 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 × 15 km grid, while observed values are precisely located in space and time, and are therefore subjects to small scale variabilities which are either unresolved by the model (e.g.small scale 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. 3c) associated to 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 2 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 skills 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. 2c), which is very similar to the Introduction

Conclusions References
Tables Figures

Back Close
Full seasonal cycle of the oxygen saturation state of the sea surface (Fig. 8).This signal in the bias is due to the underestimation by the model of the air-sea flux kinetic as explained in Sect.4.2.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. 3), 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 repartition 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) are not represented by the model.

Is the model able to represent the location of oxygen depleted waters?
The repartition of observed and simulated vertical oxygen profiles along the bathymetric axis (Fig. 4a A first mapping of the area of hypoxia over the NWS was proposed by Zaitsev (1997) (redrawn on Fig. 5a).They 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, on a smaller extent, along the Romanian coast.Unfortunately Zaitsev (1997) did not mention how they derived this mapping (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. 5a, c.
The reduction of the south-eastern 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 rivers in the model.In the model, all the discharges occur through the Chilia branch which is the most important one (Fig. 1).It is noteworthy, however, that the model predicts hypoxia in the Constanta Bay in agreement with Zaitsev (1997) without the direct inputs of the 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 in the seasonal (Fig. 2a) and interannual (Fig. 3a) repartitions.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  allows the settlement of bottom hypoxic conditions by preventing bottom oxygen depleted waters to mix with surface ventilated waters.In this region, the water column is permanently stratified, due to the presence of the river plume all along the year, with an enhanced stratification from the beginning of spring (early April) until mid-fall (mid-October), due to the warming of surface waters.
The contribution of thermal and haline effects to the stratification of the water column in the northern part of the BS-NWS is illustrated in Fig. 7. Haline stratification is present during the whole year and peaks in May, when river discharges are maximum while thermal stratification is only effective in summer, and maintains a stratification maximum until August.
The duration of this enhanced stratification period varies interannually and is a key factor in determining the severity of bottom hypoxia (see Sect. 5).A criteria on the Brunt-Vaïasala frequency −g δρ ρ > 0.05 s −1 is used to define the period of enhanced stratification, referred to as the stratification (resp.mixing) period in the following (Fig. 7).The mixed layer depth z ρ , defined by a density difference of ∆ρ| 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. 8).The low rising slope of bottom temperature indicates a limited warming of the bottom layer by the surface layer, which has as consequence that the average summer bottom temperature is determined by the winter/spring Sea Surface Temperature (SST) rather than by the summer SST.This is important since bottom temperature in summer acts on metabolic and respiration rates below the mixed layer depth and on 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.Introduction

Conclusions References
Tables Figures

Back Close
Full

Oxygen solubility and surface fluxes
The surface oxygen concentration is constantly regulated by the air-sea fluxes (Fig. 8) which tend to bring the surface concentration towards the saturation value.This saturation value depends on the oxygen solubility and is therefore mainly determined by the surface temperature.Therefore, the seasonal SST signal explains the larger part of the surface concentration variability: the seasonal variations of spatial averaged difference, [DOX] surf − [DOX] surf,sat , is less than half of that of [DOX] surf (σ = 13 and 30 mmol m −3 , resp.).
In spring and summer, surface waters are over-saturated in oxygen due to the decreased solubility as SST increases and to photosynthesis by phytoplankton bloom.This leads to a release of oxygen towards the atmosphere (Fig. 8).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 under-saturated and oxygen is uptaken from the atmosphere.This ventilation process lasts until the next spring.
The overestimation of the surface oxygen concentration in spring and the underestimation in fall (as described by the seasonal bias shown in Fig. 2) occur when the surface water are over-saturated and under-saturated, respectively.
This suggests an underestimation of the air-sea fluxes intensity due to the use of averaged wind (hence of lower magnitude) to compute the exchange coefficients appearing in the formulation of the air-sea flux (Ho et al., 2006).
The net oxygen transfer with the atmosphere is balanced when integrated over the year (net 2 Gmol influx from the atmosphere, Fig. 9) and its net contribution to the annual oxygen budget may be seen marginal.The winter/spring ventilation of the water column, however, is a critical process for the occurrence of bottom hypoxia in summer as it fixes the concentration of oxygen in the bottom waters before their occlusion by the thermocline.

Conclusions References
Tables Figures

Back Close
Full

Biochemical processes and annual oxygen budget
Figure 9 gives 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).
Over 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 river discharges, as detailed in Sect.5.3.While the stratification period accounts for 55 % of the autotroph oxygen production, it gathers 62 % of the oxygen consumption by respiration and benthic diagenesis.
The summer depletion of oxygen content below the mixed layer depth is mainly due to the different vertical repartition of photosynthesis and respiration during this period: 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).80 % of the summer chemical consumption occurs below the mixed layer depth.
Finally, the concentration of oxygen in 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.Introduction

Conclusions References
Tables Figures

Back Close
Full

Diagnostic variables of hypoxia
The status of hypoxia in a system is usually represented by referring to the average oxygen concentration measured for a particular month (Fig. 10) or to the extension of the surface area characterized by low level of oxygen (Fig. 5).The severity of hypoxia for the ecosystem will additionally depend on the duration of the hypoxic event (Vaquer-Sunyer and Duarte, 2008), a variable difficult to assess from observations (except when using time series obtained at a fixed station but in this case the spatial dimension of the problem can not be appraised).For the following investigations, we define an annual hypoxic index which merges the spatio-temporal dimension of the problem: where A(t) is the hypoxic area at time t and the integral is computed over the year.

Empiric linear models
Several studies used empirical models (e.g.linear or multiple linear regressions) to link the interannual variability of a diagnostic of hypoxia to a selection of external drivers (Conley et al., 2007(Conley et al., , 2009;;Turner et al., 2011).These studies used observations and select 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.Using this ideal "data set", we propose to derive an empirical relationship linking the intensity of hypoxia as defined in Eq. ( 1) to its drivers.Compared to the model presented here, this empirical relationship has not the complex mechanistic description of processes but has the following advantages: (1) it requires a reduced number of data for the prediction of hypoxia, (2) it is more easily

Conclusions References
Tables Figures

Back Close
Full comparable to studies addressing others areas by for instance 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 site 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 be fitted.This is a inherent limitation of a statistical model.
A multiple linear regression procedure (e.g.Legendre and Legendre, 1998) consists in deriving a regression equation expressing the response variable Y (i.e. the intensity of hypoxia) as a linear function of a set of predictors (or explanatory variables) X 1 , . . ., X N : Where Y is the estimated value of Y and the coefficients a 1 , . . ., a N are computed automatically in order to minimize the misfits (Y − Y ) 2 .
Among a preselected set of potential predictors issued from the model (i.e.winter SST, 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, river discharges in water, nitrate and phosphate) a stepwise regression procedure (backward elimination based on p = 0.05 significance threshold) is used to identify a subset of predictors for which no additional predictor will increase significantly the coefficient of determination of the regression (R 2 , computed as the ratio of the variance of the response explained by the predictors Y and the true variance of Y ) (Legendre and Legendre, 1998).Introduction

Conclusions References
Tables Figures

Back Close
Full It has been found (Fig. 11) that the hypoxia intensity index, H, can be retrieved from a linear combinations of the annual nitrate discharge (N), March SST(T ), the duration of the stratification period (D) and the winter sediments stock of semi-labile detritus (C).After standardization, this regression equation writes: This regression equation explains 79 % (R 2 = 0.79) of the variability of H.Moreover, when using standardized predictors, the standard regression coefficients a i allow to compare the importance of the different predictors in determining the variability of H.
However, the confidence intervals associated to these coefficients must be considered to ascertain within satisfying probabilities the preponderant role of a variable over another.For instance, Fig. 11 shows that the higher value of a T compared to a N is not significant at the p > 0.9 confidence and does not allow to conclude that the winter ventilation explains more of the variability of H compared to the annual nitrate discharge.
In order to disentangle the impact on hypoxia of climatic factors and eutrophication, we perform a partial regression and estimate how much of the variation of hypoxia can be attributed exclusively to climate drivers once the effect of eutrophication has been taken into account.This indicates that eutrophication and climatic drivers have a commensurate impact on hypoxia, each accounting for ∼ 36 % of H. variance (R 2 = 0.35, 0.37, resp.).Interestingly, the impact from temperature introduce a periodicity 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
Using the value of the regression coefficients to interpret the relative role of the predictors in governing the variability of the response requires that the selected predictors are not significantly correlated (Legendre and Legendre, 1998).This is the case here Introduction

Conclusions References
Tables Figures

Back Close
Full since each of them is related to different processes.The interpretation of the predictors selected by the stepwise regression procedure can be done as follows.
Nitrate discharge 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 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. 4 and is due to two main reasons: (1) the March SST fixes the temperature of bottom waters from spring to the end of summer (i.e.period of stratification) (correlation between March SST and bottom summer temperature is ρ = 0.82) and hence affects the rates of organic matter respiration, (2) the March SST controls oxygen solubility in the last days of the mixed period and hence fixes the oxygen content of bottom waters in spring (ρ = 0.85), before their occlusion by the thermocline.
The duration of the stratified period, D, responds to mean summer SST (ρ = 0.64), and to the freshwater budget (i.e.river discharge, precipitation and evaporation).D determines the time scale of isolation of bottom waters.As the lowest oxygen concentrations are reached at the end of summer, anomalous extensions of the stratification period cause increased damages on the ecosystem.In particular, D is related to the annual release of reduced compounds from the sediments (ρ = 0.65).
While a part of sediment organic matter is labile and respired at the time scale of months, a part of it is semi-labile and is respired more slowly.During years of eutrophication, the semi-labile stock accumulates year after year and continues to cause oxygen consumption several years after an eutrophication event (even if management strategies reduce nutrient discharge).The predictor C reflects this semi-labile component and the impact of the sediments memory.The mean summer benthic oxygen demand is significantly correlated to C (ρ = 0.82) indicating the important part of this accumulated stock in benthic diagenesis.Introduction

Conclusions References
Tables Figures

Back Close
Full

Impact of sediment memory on the time scale of recovery
The evolution of C as a function of the annual nitrate discharge history can be reproduced by a parametric model using, in addition to N, the March SST (as a proxy of bottom temperature, Sect.4) and the history of the silicate to nitrate ratio (i.e.Si : N) of river discharged nutrients: where 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 init exp −β 0 y , β has the same interpretation as β 0 but additionally considers the effect of temperature through a Q 10 formulation, α represents the part of riverine discharge that feeds C, and depends on the ratio Si : N of riverine discharge.This ratio is a proxy of the proportion of diatoms and flagellates 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 using Eq. ( 4) and the fitted parameters provides a coefficient of determination r 2 = 0.95 (Fig. 12).
The time scale of sediment memory after a deposition event can be estimated from the parameter β 0 which reaches a value of 0.16(±0.02,p > 0.1) corresponding Introduction

Conclusions References
Tables Figures

Back Close
Full N (noted N eq ).To derive an estimation for the values C eq (N eq ), Eq. ( 4) 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, we continue the integration for 10 additional years maintaining a constant level (N = N eq ). Figure 12 shows the levels of C eq reached for different values of N eq .We can see 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 how to avoid it, or at least to restrain trough riverine discharge restriction policies the phenomenon of hypoxia below a threshold H 0 ensuring an healthy functioning of the coastal and marine ecosystem, as defined by the Marine Strategy Framework Directive (MSFD) defined by the European Union (Cardoso et al., 2010).
We propose to use the coefficients of the regression Eq. ( 3) to predict the level of hypoxia corresponding to different levels of nitrate discharge.The aim is to derive from this regression a general relationship between the level of hypoxia reached at equilibrium (noted H eq ), i.e. when the system is submitted yearly to a repeated level of nitrate discharge, N eq , with a benthic stock in equilibrium with this discharge (as defined in Sect.5.3) and when the average atmospheric conditions (March SST, duration of stratification) remain unchanged compared to the actual situation.From Eq. (3), we first derive H eutroph. as the variability of the hypoxic index imparted to eutrophication (due to C, N), and thus obtained by removing the variability due to meteorological conditions: Since H eq represents the level of hypoxia reached with a sediment response in equilibrium with the repeated level of nitrate discharge, H eq can be derived from Eq. ( 7) by

BGD Introduction Conclusions References
Tables Figures

Back Close
Full removing the effect of a non-equilibrated sediments: Equation ( 4) gives a value C eq corresponding to each value N eq (Fig. 12, which can be used within Eq. ( 8) to provide a relationship between H eq and N eq .Figure 13 illustrates this relationship when using in Eq. ( 8) the time series of H and N given by the 3-D model.The resulting distribution of H eq can be described by a power law (H = H 0 +a•N b , R 2 = 0.85, Fig. 13) and provides the average level of hypoxia corresponding to a long lasting discharge level.Atmospheric variability, however, has a commensurate effect on hypoxia (Sect.5.2), indicating that even in the framework of a pertinent nutrient management, the system is not free of risk of hypoxia occurrence due to atmospheric variability.The effect on H eq of one standard deviation of the atmospheric drivers from their average value is indicated on Fig. 13 by adding the terms (a

The recovery process: representativity of monitoring coverage
There are different and contradictory ideas about the recovering process from hypoxia on the Black Sea NWS after eutrophication in the 90's and 2000's.For instance, 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 discharge while a quite different picture is presented in the Ukrainian national report UkrSCES (2002).
Unfortunately, this publication does not specify neither the data, nor the methods used to estimate the extent of hypoxia.The presented values before 1993 match those presented by Zaitsev (1993), in which no more information are given concerning data coverage and computation methods.
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. 5  studies (and additionally by the model presented in this paper) and evidences two main important differences: (1) the maximum surface is 2 time smaller (20 000 km 2 , in 1983) in UkrSCES ( 2002), (2) the surface in 1999 and 2000 is similar to that of 1989 (14 000 km 2 ) suggesting no recovery contrary to Mee (2006).We can note that the model estimations are closer to the Ukrainian report concerning the recovery of the system from hypoxia during the 1990's period.These different interpretations of a same phenomenon may come from the use of different data sets to base their conclusions.Unfortunately, we were not able to find some 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 are given concerning data coverage and computation methods.
To illustrate the sensitivity to the distribution of data of the conclusions concerning hypoxia and its recovery, Fig. 10 compares the interannual evolution of the summer (July-August) oxygen concentration obtained from available data, model outputs restricted to data locations and 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. 2), (2) the model bias is minimal for these months (Fig. 2c), (3) bottom hypoxic events are frequent.
Observation and model results obtained at data points both indicate a rise of the average summer oxygen concentrations during the 90's.However, this trend is not visible anymore when model results are averaged over the whole area, indicating that this apparent recovery may be nonfactual and rather results from the inhomogeneity of monitoring after the mid-nineties.
Similarly, Fig. 3 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 during 1995-2000, but a careful look on 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 to Introduction

Conclusions References
Tables Figures

Back Close
Full appraise the "real" interannual variability of oxygen and biased the interpretation of hypoxia and its potential recovery.Rather than the consensual picture of a recovery from hypoxia closely following the decrease in nutrients discharge in the Black Sea NWS, the present study indicates a delayed recovery, attributed to the sediments memory of eutrophication events (Sect.5.2.1).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).
In addition, the present study attribute a atmospheric impact on hypoxia which is commensurate to that of eutrophication (Sect.5.2).This can explain why hypoxia is still present nowadays on the NWS, as supported directly by in-situ measurements (Friedrich, J., 2010, unpublished data collected during the EU HYPOX project and Mihailova, M. E., 2010, unpublished data) and indirectly by the occurrences of massive fish mortality (D. Banaru, personal communication, 2012).

Discussion and perspectives
Lethal and non-lethal effects of low oxygen concentrations on living organisms strongly depend on the tolerance of the considered species, on the level of oxygen concentration and the duration of the hypoxic event (Vaquer-Sunyer and Duarte, 2008).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 a dedicated and site-specific study.
The linear responses of hypoxia intensity to eutrophication and climatic variations (Sect.5.2) are not expected to hold far out of the range covered by the present hindcast model experiment.We suspect that for low nutrients discharges conditions, not encountered during the period 1981-2009 and for which the model could not have been validated, the effect of benthic macro-biota could play an important role in regulating Introduction

Conclusions References
Tables Figures

Back Close
Full hypoxia.Nevertheless, Fig. 13 gives an indication of the sensitivity of hypoxia to nitrate discharge that can be used to construct scenarios model experiments.
To be realistic, these scenarios of nutrients emissions should include estimations of population growth, socio-economics and climate change scenario.It must be noted, in particular, that the effect of temperature is interpreted here as a natural variability, which, while indicative, may not reflect the long term effect of a firm warming trend on the system.

Conclusions
The Black Sea NWS is affected by annual seasonal bottom hypoxia since decades.In the worse case, the hypoxic zone covers all the shelf North of 45.5 • N and extends along the Romanian coast in the Bay of Constanta.The oxygen-depleted 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 80's, the reduction of the surface of the 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 incorporating 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 physical and biogeochemical processes submitted to the historical nutrient discharge 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 World ocean data base (18735 points, 1980-2000) and from the Black Sea Introduction

Conclusions References
Tables Figures

Back Close
Full commission (733 points, 2000-2009).In particular, the comparison of the observed and simulated number of occurrence of hypoxic points ([O 2 ]< 65 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) due to different spatial coverage of the data.The spatial heterogeneity of the data among years could also biased the conclusions on the status of hypoxia.For instance, a concentration of observations south of the typical hypoxia area lead in some cases to conclude that a recovery occurred after the early 90's, shortly after the reduction of nutrient discharge.This quick recovery is not reflected by the model results which rather depict a recovery delayed by the inertia of the sediments layer and mitigated by the important role of atmospheric conditions (see below).There is thus an urgent need of dedicated oxygen monitoring effort concentrated on the spatial and temporal cover of hypoxia including the management of an international accessible database.For instance, almost no observations are available in the period and region of seasonal hypoxia occurrence after 2000.
The coupled model has been used to understand the dynamic 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 affects hypoxia intensity by (1) controlling the oxygen solubility in the last period of vertical mixing, hence the oxygen concentration of bottom waters before their occlusion, (2) settling the temperature of bottom waters during the whole stratified period which affects respiration rates and consumption of oxygen below the mixed layer depth.
Pelagic respiration and benthic diagenesis decrease oxygen concentration in the bottom waters until the end of the stratified period (October).While pelagic respiration Introduction

Conclusions References
Tables Figures

Back Close
Full is mostly related to nitrate discharge, benthic diagenesis of accumulated organic sediments introduces a memory of eutrophication, and is estimated to keep on sustaining hypoxia with a typical time scale of 9.3 yr after high nutrient discharge events.The lowest oxygen concentrations are reached at the end of the stratification period.An over-lasting stratification period therefore significantly enhances the severity of hypoxia and its damages on the benthic biota directly and also indirectly through the release of reduced substances such as hydrogen sulfide.
For each year, we define a diagnostic of hypoxia as the surface extent of hypoxia integrated over its duration in time.
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 the hypoxia diagnostic.Among this large set, a stepwise regression procedure isolates as driving predictors the annual discharge of nitrate, the semi-labile carbon content in the sediments, the SST in March and the duration of the stratification period.These four predictors allow to reconstruct 78 % of the hypoxia intensity interannual variability.Partial regressions show that the drivers linked to eutrophication (i.e.nitrate discharge, sediment organic content) and to climate (March SST, stratification duration) have commensurate effects on hypoxia, explaining respectively 35 and 37 % of the variability of hypoxia intensity.
From this regression we derive a non-linear relationship that gives the level of hypoxia reached for a range of nitrate discharge, assuming that average atmospheric conditions remain at the level recorded during the 1981-2009 period and that the sediment is in equilibrium with the nitrate discharge.
This relationship has the advantage of simplicity but should be used with care outside its conditions of derivation.This is the case for low levels of nitrate discharge not encountered during the period 1981-2009.This relationships could be used to design further experiments testing more precisely nutrient discharge restriction policies along with realistic socio-economics and Introduction

Conclusions References
Tables Figures

Back Close
Full  Full  Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | et al.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | tions average, P B = P −O O • 100.It indicates a global overestimation (resp.underestimation) when positive (resp.negative) and the model is considered as excellent (in that regard) when |P B| < 10 % (Allen et al., 2007).
O) 2 , relates the model errors to the variability in the observations.The squaring increases the weight of big errors.Allen et al. (2007) classifies the performance levels as > 0.65 excellent; [0.65, 0.5] very good; [0.5, 0.2] good , and < 0.2 poor.-The widely used Pearson correlation coefficient, ρ, reflects whether P and O varies accordingly around their respective average.Significant correlations (p < 0.01, two-tailed distributions) can be concluded for ρ > 0.71 and 0.47 for datasets of 12 and 29 elements, respectively, which corresponds to the number of months and years considered in Sect.3Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | , b, respectively) indicates in both cases a zone of low oxygen content along the bottom between approximately the 15 m and the 45 m isobaths.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 (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.Discussion Paper | Discussion Paper | Discussion Paper | Figure6shows the seasonal evolution of a vertical oxygen profile simulated at a site affected by seasonal hypoxia (point A in Fig.1).This figure illustrates how stratification kg m −3(Kara et al., 2009) reaches ∼ 10 m in average during the stratification period.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | compares the interannual variability of the extent of hypoxic areas provided by these two Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | climate change scenario, with the aim of reaching a sustainable level of seasonal hypoxia meeting European marine strategy requirements.Zhang, J., Gilbert, D., Gooday, A. J., Levin, L., Naqvi, S. W. A., Middelburg, J. J., Scranton, M., Ekau, W., Pe ña, A., Dewitte, B., Oguz, T., Monteiro, P. M. S., Urban, E., Rabalais, N. N., Ittekkot, V., Kemp, W. M., Ulloa, O., Elmgren, R., Escobar-Briones, E., and Van der Plas, A. K.: Natural and human-induced hypoxia and consequences for coastal areas: synthesis and future development, Biogeosciences, 7, 1443-1467, doi:10.5194/bg-DiscussionPaper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .Fig. 7 .Fig. 9 .Fig. 11 .
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, and the river discharges considered in this study.The dot (A) locates the profile presented in Fig. 6.

Table 1 .
Number of observations available from the World Ocean Database (WOD) and Black Sea Commission (BSC) datasets.The second column considers all the data, while the next columns only consider subsets selected following hypoxia-related criterion (see the text for details): data collected at a water depth Z > 17 m (3rd column) and below the summer mixed layer (observation depth, z > 15 m, 4th column) and during the months of strongest hypoxia (August, September, October) (second line of each data set).Bold number indicates the data used for validation in Sect.3.

Table 2 .
Model skill statistics obtained when considering all the observation/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, r σ : standard deviation ratio.The numbers in parentheses are obtained when only considering the WOD database (Sect.2.2).