Black Sea biogeochemistry: Response to decadal atmospheric variability during 1960–2000 inferred from numerical modeling

Marine Environmental Research 77 (2012) 90e102

Contents lists available at SciVerse ScienceDirect

Marine Environmental Research journal homepage: www.elsevier.com/locate/marenvrev

Black Sea biogeochemistry: Response to decadal atmospheric variability during 1960e2000 inferred from numerical modeling Yunchang He a, *, Emil V. Stanev a, Evgeniy Yakushev b, Joanna Staneva a a b

HZG, Geesthacht 21502, Germany NIVA, Gaustadalleen 21, NO-0349 Olso, Norway

a r t i c l e i n f o

a b s t r a c t

Article history: Received 9 June 2011 Received in revised form 15 February 2012 Accepted 23 February 2012

The long-term variability of the physical and biochemical structure of oxic and suboxic layers in the Black Sea was studied using a one-dimensional coupled hydrophysical and biogeochemical model. The focus was on the correlation between atmospheric forcing (2 m air temperature and dew point temperature, surface level pressure, surface wind) affected by the North Atlantic Oscillation in and the regional responses. The quality of model performance was demonstrated using observed vertical and temporal distribution of biogeochemical variables. It was shown that during 1960e2000, the long-term variability of simulated winter-mean SST in the Black Sea correlated reasonably well with the variability of 2 m air temperature. Furthermore, the thermal state of the upper ocean impacted largely on the variability of biogeochemical variables, such as oxygen, nitrate and phytoplankton concentration. The tele-connection between North Atlantic Oscillation and Black Sea biogeochemistry was manifested in a different way for the specific time-interval 1960e2000; the corresponding regime shifts were thus associated with the large scale forcing. One such extreme event occurred in 1976 leading to a pronounced shift in the oxygen and hydrogen sulfide state. Crown Copyright Ó 2012 Published by Elsevier Ltd. All rights reserved.

Keywords: Black Sea Biogeochemical model Atmospheric forcing Decadal variability Oxygen

1. Introduction Permanent anoxic conditions have been observed in lakes, fjords, and ocean basins like Black Sea, Baltic Sea, Cariaco Trench, characterized by extremely stable stratification. In the Black Sea, the biogeochemical system adjusted to the physical stratification in a unique way, building thus the largest anoxic environment in the world. Major changes of biogeochemical variables were observed across levels of constant density, while gradients along density surfaces were extremely small (Spencer and Brewer, 1971). This specificity of the vertical stratification justified the approach used here, which was to address dominating processes using onedimensional numerical model. Oxygen in the Black Sea is almost at atmospheric saturation in the upper 30e40 m; it drops sharply below the detection limit of 3 mM at about the depth of 60 m (Fig. 1). The first appearance of

Abbreviations: NAO, North Atlantic Oscillations; OM, Organic Matter; AT, Air Temperature; DPT, Dew Point Temperature; CIL, Cold Intermediate Layer; GOTM, General Ocean Turbulent Model; ROLM, RedOx Layer Model; SST, Sea Surface Temperature. * Corresponding author. Tel.: þ49 0 4152 87 1523; fax: þ49 0 4152 87 1565. E-mail addresses: [email protected], [email protected] (Y. He).

hydrogen sulfide is at about 80 m and then it increases continuously to 42 mM at about 150 m. Vertical distribution of nitrate is dominated by a pronounced maximum at about 60 m; ammonia reveals a continuous increase starting approximately at the same level where sulfide occurs. The same profiles as in Fig. 1a, c are shown in Fig. 1b, d in density coordinates because such representation is typical for many studies on the Black Sea biogeochemistry. It is noteworthy that the biogeochemical structure of the water column is characterized by no-overlap of dissolved oxygen and hydrogen sulfide layers. The corresponding layer between oxic and anoxic waters is known as the suboxic layer (the zone without oxygen and hydrogen sulfide, Murray et al., 2007). Yakushev et al. (2007) showed that the manganese cycle in the Black Sea could explain this layering. The present-day vertical structure is very stable and has been established as a consequence of the re-connection between the Black Sea and Mediterranean. Some analyses of long-term variations of chemical and ecological systems in the Black Sea (Konovalov and Murray, 2002; Oguz, 2005) reveal regime shifts, which could be due to either anthropogenic, or climatic forcing. Gregoire and Soetaert (2010) quantified the chemical flows including oxygen and sulfide in the entire Black Sea through a coupled physical-biogeochemical model.

0141-1136/$ e see front matter Crown Copyright Ó 2012 Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.marenvres.2012.02.007

Y. He et al. / Marine Environmental Research 77 (2012) 90e102

91

Fig. 1. Vertical profiles of chemical and ecological variables in the Black Sea plotted in depth (left) and density (right) coordinates. Data used to plot this figure originate from (Leg 8, Stn 5) of R/V Knorr Cruise in the Black Sea during 2003 (http://www.ocean.washington.edu/cruises/Knorr2003/data/Vol_Oxygen_Sulfide/08-5.html).

So far, to our knowledge, numerical modeling of long-term changes has not been used to identify different responses of the Black Sea biogeochemical system to external and internal forcing. Identification of major events dominating interannual variability of physical state is addressed in a number of publications (e.g. Staneva and Stanev, 2002). Yet, the response of biogeochemistry to interannual variability has not been fully considered in numerical modeling. This motivates us to address here the possible regime changes of biogeochemical structure resulting from the changes in the atmospheric state. This issue is relevant to the response of the Black Sea to global warming, which could lead to increasing stability of vertical stratification, reduced ventilation of deep layers and finally changes in the oxygen balance. Identification of possible regime changes constitutes the first goal of present study. Our second goal is to test the capability of a complex biogeochemical model to reproduce interannual variability triggered by atmospheric forcing. In the present paper we use a 1-D model. This simplification of the coupled physical-biogeochemical system presents a first step, which is needed in order to develop initial understanding of the possible reaction of natural anoxic-ocean

systems to atmospheric forcing. Testing the performance of the 1-D biogeochemical model and validation against observations will also be illustrated here with a main focus on the upper layer variability. The paper is structured as follows. We first describe the numerical model and experiments carried out, followed by a discussion of results and short conclusions. 2. Model and data 2.1. Description of the model The biogeochemical model used in this study is the Redox-Layer Model (ROLM, Yakushev et al., 2007). This model includes production and decay of organic matter and reduction and oxidation of nitrogen, sulfur, manganese, and iron. The individual state variables and the corresponding abbreviations used for them are given in Table 1 (see also Fig. 2). Unlike other existing models for the Black Sea (Oguz et al., 1998; Gregoire et al., 2008), organic mater production due to both photosynthesis and chemosynthesis is

92

Y. He et al. / Marine Environmental Research 77 (2012) 90e102

Table 1 State variables and their abbreviations used in the text. Variable

Abbreviation

Dissolved Oxygen Hydrogen Sulfide Elemental Sulfur Thiosulfate Sulfate Ammonia Nitrite Nitrate Particulate Organic Nitrogen Dissolved Organic Nitrogen Phosphate Particulate Organic Phosphorus Dissolved Organic Phosphorus Bivalent Manganese Trivalent Manganese Quadrivalent Manganese Bivalent Iron Trivalent Iron Phytoplankton Zooplankton Aerobic Heterotrophic Bacteria Aerobic Autotrophic Bacteria Anaerobic Heterotrophic Bacteria Anaerobic Autotrophic Bacteria

O2 H2S S0 S2O3 SO4 NH4 NO2 NO3 PON DON PO4 POP DOP MnII MnIII MnIV FeII FeIII Phy Zoo Bhe Bae Bha Baa

included, thus the feedbacks between the upward fluxes of nutrients and organic mater production are accounted for. A detailed description of the model is presented by Yakushev et al. (2007). ROLM was online coupled with the General Ocean Turbulent Model (GOTM; Burchard et al., 1999). The latter is a onedimensional water column model for the most important hydrodynamic and thermodynamic processes related to vertical mixing. The basic set of Eqs. (1e7) for temperature T, salinity S, mean horizontal velocity components u and v, pressure p, turbulent kinetic energy (TKE) k and the TKE dissipation rate 3 are:


vT  gT ðT  T0 Þ þy vz  
vS vS v y0t þ y0  gS ðS  S0 Þ ¼ vt vz vz vT v ¼ vt vz





y0t

0

(1)

(4)

vP ¼ g r vz

(5)

vk v2 k ¼ yk 2 þ P þ B  3 vt vz

(6)

v3 v2 3 3 ¼ y3 2 þ ðc3 1 P þ c3 3 B  c3 2 3 Þ vt k vz

(7)

where y0t and y0 are molecular and turbulent diffusivities of temperature and salinity, y and yt are molecular and turbulent viscosity (momentum), f is the Coriolis parameter, and y3 and yk are the turbulent diffusivities of energy dissipation and TKE, respectively. Relaxation with time scale 1/g(T, S) towards prescribed profiles (T, S) was enabled in the deep layers to compensate for the unresolved processes in the one-dimensional model such as gravity currents originating from the Bosphorus Straits, coastal-deep ocean exchange across the Black Sea Rim Current, deep ocean upwelling, etc. The shear stress production P and the buoyancy production B are:

P ¼ yt



  2  vu 2 vv þ vz vz

(2)

(3)

(8)

B ¼ y0t N2

(9)

Where N2 ¼ ðg=r0 Þðvr=vzÞ is the Brunt-Väisälä frequency. The turbulent viscosity and diffusivity are computed using the relation of Kolmogorov and Prandtl

k2

k2

3

3

yt ¼ cm ; y0t ¼ c0m

The turbulent diffusivities for k and



vu v vu ðyt þ yÞ ¼ þ fv vt vz vz

  vv v vv ðyt þ yÞ  fu ¼ vt vz vz

yk ¼

cm k2

sk

3

; y3 ¼

cm k2

s3

3

(10) 3

are

(11)

All parameters in the above formulas are listed in Table 2, see also (Burchard et al., 1999).

Fig. 2. Flow-chart of the biogeochemical processes represented in ROLM (from Yakushev et al., 2007).

Y. He et al. / Marine Environmental Research 77 (2012) 90e102 Table 2 Parameters of the K-3 Model (Rodi, 1987).

Table 3 Parameter names, notations and used values in the model (for more details, see Yakushev et al., 2006).

Model constant

Value

Definition

c3 1 C3 2 c3 3 (B < 0)

1.44 1.92 0.4

c3 3 (B > 0)

1

s3 sk

1.3 1

Empirical coefficient of dissipation equation Empirical coefficient of dissipation equation Empirical coefficient of dissipation equation for stable stratification Empirical coefficient of dissipation equation for unstable stratification Schmidt number for TKE diffusivity Schmidt number for TKE diffusivity

The surface fluxes were calculated by means of bulk aerodynamic formulas described by Roussenov et al. (1995) using atmospheric data at sea surface and current simulated SST. ROLM solves one-dimensional vertical diffusion equations for 24 non-conservative substances of the type:

vCi v ¼ vt vz



y0t þ y


vCi vz

 

v ððWC þ WM Þ*Ci Þ þ RCi ; vz

93

(12)

where Ci is the concentration of i-th model variable, WC is the sinking velocity of particulate matter; WM accounts for the contribution of settling of Mn hydroxides, and RCi describes biogeochemical production or consumption. Parameter values used in ROLM are given in Yakushev et al. (2007). Recently, some parameters concerning phytoplankton growth, respiration, mortality and excretion rates have been updated (Yakushev, personal communication), which enabled a better consistence of simulations with the observations. These parameters and their values used in the present study are given in Table 3. Coupling between physical and biogeochemical model uses standard numerical routines described in GOTM (Burchard et al., 1999), and more details about that are given in Yakushev et al. (2007). The numerical integration of biogeochemical state variables uses an Euler scheme with process splitting. The problem considered in the present paper is similar to the one addressed by Burchard et al. (2005), who used a much simpler ecosystem model (NPZD-type of model). Therefore we will not present below a detailed description of the numerics but will only refer to the major result of this study, which states that the second-order modified PatankareRungeeKutta scheme gives good results even for long integration time steps. Additional model runs to check the dependence of the solution originating from the coupled ROLMeGOTM model upon the numerics have also been carried out using available GOTM algorithms. The comparison of different numerics did not reveal major differences. Because the present study can be considered as an initial research step to address longterm variability in the Black Sea, we preferred the simpler and slightly smoother solution provided by the first order method. As far as the computational efficiency of different schemes is concerned this is not an issue for 1-D models, but could become important consideration in more computationally demanding 3Dmodels. The vertical resolution was 1 m, which enables sufficient resolution for the thin suboxic layer, the physical and biochemical time steps were 10 and 20 min, correspondingly. The maximum depth considered in the model was 200 m. 2.2. Forcing and initial data 2.2.1. Forcing data Forcing data for the physical model consisted of radiation (including clouds), 10 m wind speed, 2 m air temperature and dew point temperature, as well as sea level pressure. They were

Parameter name Phytoplankton Maximum specific growth rate Extinction coefficient Incident light Optical light Specific respiration rate Specific rate of mortality Specific rate of excretion Zooplankton Maximum specific rate of grazing of Zoo on Phy Half-saturation constant for the grazing of Zoo on Phy for Phy/Zoo ratio Maximum specific rate of grazing of Zoo on POP Half-saturation constant for the grazing of Zoo on POP in dependence to ratio POP/Zoo Specific respiration rate Maximum specific rate of mortality of Zoo Food absorbency for Zoo Ratio between dissolved and particulate excretes of zooplankton OM Specific rate of decomposition of POM to DOM P Half-saturation constant for uptake of PO4 N Strength of ammonium inhibition of nitrate uptake constant Half-saturation constant for uptake of NO3 þ NO2 Half-saturation constant for uptake of NH4 Specific rate of decomposition of DON Specific rate of decomposition of PON Specific rate of the 1st stage of nitrification Specific rate of the 2nd stage of nitrification Specific rate of the 1st stage of denitrification Specific rate of the 2nd stage of denitrification Specific rate of thiodenitrification Specific rate of anammox S Specific rate of oxidation of H2S with O2 Specific rate of oxidation of S0 with O2 Specific rate of oxidation of S2O3 with O2 Specific rate of OM sulfate reduction with sulfate Specific rate of OM sulfate reduction with thiosulfate Mn MnII oxidation with O2 constant MnIII oxidation with O2 constant MnIV reduction with sulfide constant MnIII reduction with sulfide constant Fe Fe oxidation with O2 constant Fe oxidation with NO3 constant Fe oxidation with MnIV constant FeIII reduction by sulfide

Notation

Value

KNF k_Erolv I0 Iopt KFN KFP KFD

2.3 d1 0.1 80 25 0.01 d1 0.02 d1 0.02 d1

KFZ KF

0.5 d1 200

KPZ KPP

0.6 d1 200

KZN KZP UZ HZ

0.05 d1 0.001 d1 0.7 0.6

KPD

0.01 d1

KPO4

0.01 mM

Kpsi

1.46

KNO3

0.20 mM

KNH4 KND4 KNP4 KN42 KN23 KN32 KN24 KT K_annamox

0.02 mM 0.1 d1 0.04 d1 1.70 d1 8.0 d1 0.02 d1 0.01 d1 0.9 mM1 0.04 d1

K_hs_ox K_s0_ox K_s23_ox K_s4_rd

0.2 d1 4.0 d1 1.5 d1 107 d1

K_s23_rd

1.0 d1

K_mn_ox K_mn_ox2 K_mn_rd K_mn_rd2

2 d1 18 d1 22 d1 2 d1

K_fe_ox K_fe_nox K_fe_mnox K_fe_rd

4 d1 1 d1 5 d1 0.05 d1

produced by the European Center for Medium-Range Weather Forecasts (ECMWF) and were kindly made available for the period 1958e2002 in the frame of SESAME (Southern European Seas: Assessing and Modeling Ecosystem changes) project by the INGV, Italy. The temporal sampling rate was 6 h. The data had spatial resolution of 0.25 . From these data we extracted model forcing corresponding to an open ocean location (32.625 E, 43.177 N). The comparison with the atmospheric data used by Stanev et al. (2003) demonstrated that the ones used in the present study wellreplicated the regional climate and its variability. Atmospheric variability was characterized by a pronounced seasonal cycle exemplified in Fig. 3-1 as the daily means over the

94

Y. He et al. / Marine Environmental Research 77 (2012) 90e102

Fig. 3-1. 45-year mean meteorological forcing data of: wind speed at 10 m (a), sea level pressure (SLP) and air temperature (AT) at 2 m (b), dew point temperature (DPT) at 2 m and total cloud cover (TCC) (c). Because of problems with the quality of cloud-data in the atmospheric analysis, monthly climatological cloud cover is used here, which explains the low temporal resolution of the corresponding curve. Observed meteorological data is from ECWMF.

45-year model integration-period. The curves in Fig. 3-1 represent a typical signal for temperate latitudes (maximum in temperature and minimum in atmospheric pressure, winds and clouds in summer). However, as seen in Fig. 3-2, individual years differed considerably from the climatic year (1976 is shown as an illustration). In particular, wind magnitude was two times higher than in the climatic data, air temperature was too warm and, importantly, no extreme events which are responsible for the winter mixing were noticeable. The Black Sea regime shifts appear to be sporadic events forced by the strong transient decadal perturbations, and therefore differ from the multi-decadal scale cyclic events observed in pelagic ocean ecosystems under low-frequency climatic forcing (Oguz and Gilbert, 2007). The North Atlantic oscillation (NAO) index is known to be the most prominent mode of the low-frequency variability controlling atmospheric circulation and climate over the North Atlantic and Eurasia (Marshall et al., 1997). The positive NAO index is associated with stronger north-south pressure gradient, more moisture and heat transported in Scandinavian area and cold and dry air conditions in the Mediterranean and Black Sea region. Conversely, the negative NAO index supports warmer air temperatures in the southern Europe (Fig. 4-1a). The NAO index for the period 1960e2000, defined as the normalized sea level pressure difference measured in meteorological stations located at Gibraltar and Iceland (Jones et al., 1997), was obtained from the Climate Prediction Center of National Weather Service in NOAA (http://www.cpc.noaa.gov/ products/precip/CWlink/pna/nao_index.html). Monthly values (Fig. 4-1) were averaged for winter season (DecembereFebruary). In the present study, all winter-mean atmospheric forcing data were presented similarly for DecembereFebruary. As seen in Fig. 4-1 cold and dry air intrusions dominate the Black Sea region during times of positive NAO index, and conversely, warmer air temperatures dominate the Black Sea region during times of negative NAO index. All data in Fig. 4-1 were filtered with a five-year moving average. It is clearly seen that there were three peaks of the NAO index: 1976, 1984 and 1995, correspondingly.

From the 1960s to 1976, the negative NAO index dominated and the air temperature had a decreasing trend (as well as the dew point temperature, not shown here). For the following seven years the NAO index changed from positive to negative. At the same time, the air temperature increased a little (with about 1  C), and the same trends appeared for the dew point temperature. For the third period (1985e1995), a positive NAO index was very well pronounced and air and dew point temperature started to decrease. After 1995, the situations started to inverse again. The overall similarity between the NAO index and sea level pressure is clear because of the definition of the index. Nevertheless, a correlation between two curves in Fig. 4-1b was too weak around 1985. More interesting was the correlation between regional wind and global atmospheric variability (Fig. 4-1c, d), keeping in mind that the zonal wind component was about two times stronger than the meridional one. Both wind components revealed an overall negative correlation with the NAO index. However, this correlation was positive for the meridional component during 1960e1976. Afterwards correlation was markedly negative. The response of the regional zonal wind to the NAO demonstrated a negative correlation during 1960e1976. After that point the NAO did not show a similar response until 1984. During the last period the correlation tended to be negative, but the trend was not clear. The impact of these changes on the Black Sea physical and biogeochemical system will be addressed below. 2.2.2. Initial conditions Initial conditions corresponded to the annual mean temperature and salinity profiles for the interior part of the basin. We relaxed temperature and salinity below the depth of CIL to climatology. At every time step temperature and salinity profiles used for relaxation were interpolated onto the actual model grid. This procedure ensured that the simulations remained consistent with the dominating stratification in deep layers, which cannot be reproduced with the used one-dimensional model. At the same time upper layer characteristics were free to adjust to surface forcing.

Y. He et al. / Marine Environmental Research 77 (2012) 90e102

95

Fig. 3-2. Meteorological forcing data (green line) in 1976: wind speed at 10 m (a), sea level pressure (SLP) (b) and air temperature (AT) at 2 m (c). Blue dash line shows the 45-year mean value. The shaded area shows the times when extreme events were responsible for the shift of biogeochemical process in 1976. Data source is same with Fig. 3-1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

At the sea surface all chemical and ecological constituents were set to zero except for O2, NO3 and PO4. At the lower boundary, which is at the depth of 200 m constant values were specified: NH4 ¼ 20 mM, H2S ¼ 60 mM, MnII ¼ 8 mM, FeII ¼ 0.4 mM, PO4 ¼ 4.5 mM, which corresponds to the typical Black Sea conditions.

2.3. Numerical experiments and model validation 2.3.1. Experiments We defined two experiments, which we integrated for 45 years, each forced with the data described above. The difference between them was that in the first one called experiment we used the full

Fig. 4-1. Interannual variability of air temperature (a), SLP (b) and wind components (u10-c and v10-d) in winter (DecembereFebruary) compared with the winter NAO index. Data were five-year moving averaged. The shaded area shows in 1973e1982.

96

Y. He et al. / Marine Environmental Research 77 (2012) 90e102

atmospheric signal, which contained interannual variability. In the second experiment we substituted the atmospheric forcing by the mean year forcing (45 times repetition of the mean year). Therefore we refer to this experiment as to Quasi-steady one. 2.3.2. Validation of simulations The performance of the model is illustrated below by validating the simulations against observations from the Knorr Black Sea Cruise in 2003. Because the processes in the Black Sea align to isopycnals the simulated results in Fig. 5 are presented in depth and density coordinates. Simulated vertical stratification reveals the well-known vertical distribution of temperature, salinity, oxygen and hydrogen sulfide. Overall, the comparison with observations was quite good, however some differences need to be mentioned: the cold intermediate layer (CIL) in the model was a little deeper and more diffuse than observed and the oxygenated water penetrated slightly deeper in the upper layer. However, considering that the spatial-temporal variability in the real basin was not fully included in the model, one important missing element was the upwelling in the interior displacing water masses upwards, the agreement between observations and simulations is considered very good overall. The model also reproduces the isopycnic depth and thickness of the suboxic layer. This result justifies a deeper analysis of the simulations with the aim to identify and quantify processes which are not easily explained using the available observations only. 2.4. Upper ocean response of physical variables to interannual atmospheric variability Decadal and interannual variability of the atmosphere, in particular in the European region, is markedly pronounced in winter months, however the deep layers of the Black Sea did not show pronounced sensitivity to the interannual variations of atmospheric forcing (Stanev and Peneva, 2002; Stanev et al., 2004; Tsimplis et al., 2004). This was due to the strong stratification decoupling of surface and deep layers. On the contrary, the interannual variability in the upper layers and down to the depth of cold intermediate layer (CIL) was well noticeable. As demonstrated by Staneva and Stanev (2002) the time of replenishment of CIL ranged

between 5 and 12 years. Furthermore, the magnitude of the interannual and seasonal signal of the temperature in the core of CIL was comparable. Recently, the response of a number of biogeochemical parameters to long-term atmospheric variability has been addressed by Oguz et al. (2006). In the following we describe the response of simulated variables to atmospheric forcing (Fig. 4-2). The comparison between the NAO index, from one side and SST and CIL temperature anomaly, from another, in winter during 1962e2000 clearly revealed the sensitivity of the physical system to changes in the NAO. SST anomaly tracked closely the air temperature (compare with Fig. 4-1a). Correlation with the NAO index was negative (a positive NAO index was accompanied with a negative temperature anomaly). However, this response did not follow the same path during different periods, which we identify by the three peaks of the NAO index in 1976, 1984 and 1995. From 1960s to 1973, the negative NAO index dominated; the SST and CIL temperature had a decreasing trend. During the following seven years (1974e1982) the NAO index changed from positive to negative, the trend of SST anomaly increased by 0.5  C, while the CIL temperature anomaly increased by about 0.1  C which, for the intermediate layer, is a substantial change. During 1983e1995, a positive NAO index showed a clearer negative correlation both with SST and CIL temperature. The above results can be better understood bearing in mind that the Black Sea intermediate water shows a delayed response to the atmospheric variability (Stanev et al., 2003). The response of the CIL to atmospheric variability described above can thus be explained by the fact that both the NAO index and the CIL heat content are very sensitive to the interannual variability during winter (see also Stanev and Peneva, 2002; Staneva and Stanev, 2002). Because the temporal variability of the biogeochemical system is strongly dependent upon the heat content and mixing in the upper and intermediate layer, we will focus our further analyses on: (1) the response of the biochemical characteristics to interannual variability in the severity of winters, and (2) the different correlations between forcing and response, which could identify shifts in the state of physical and biogeochemical variables of Black Sea. In summary, the interannual variability in the Black Sea diagnosed from the numerical simulations did not reveal a very clear correlation between analyzed variables and the NAO index, in

Fig. 4-2. Interannual variability of the simulated anomalies in winter (mean for DecembereFebruary) of SST (a) and CIL temperature (b) compared with the winter NAO index. Data were five-year moving averaged. The shaded area shows in 1973e1982.

Y. He et al. / Marine Environmental Research 77 (2012) 90e102

97

Fig. 5. Comparisons of vertical distribution of temperature (a), salinity (b), oxygen (c) and hydrogen sulfide (d) between simulation and observation (dash line: observation; solid line: simulation). Observations used for inter-comparison are taken from Knorr Cruise in April 2003.

particular when shorter periods are concerned. As explained by Oguz et al. (2006), the NAO is not the only climatic index relevant to the region of Black Sea, therefore a perfect correlation between this index and regional climate is not to be expected.

and that this reaction could be balanced by the flux of H2S from below. Thus, the main observed biogeochemical features of the redox layer in the central Black Sea known are well reproduced by the 1-D model.

3. Analysis of simulated biogeochemistry responses

3.2. Interannual variability: analysis on dominating processes

3.1. Vertical stratification

The vertical distribution of oxygen in the Black Sea (Fig. 1) consists of two layers. In the upper layer oxygen is 90e110% saturated, depending upon the temperature and salinity. This dependency explains the reason why seasonal variations of oxygen in this layer are well pronounced (Konovalov et al., 2001). The lower layer is the oxic/anoxic zone where oxygen decreases rapidly to zero. Long-term variations in the structure of the seasonal thermocline have been addressed in a number of studies. Due to this important characteristic of vertical stratification, which severely restricts the vertical mixing, oxygen was consumed faster than it could be replaced by vertical and lateral fluxes at a depth of around 70 m. In the real basin, vertical fluxes are usually associated with the winter mixing, which changes very strongly from year-to-year (Stanev and Peneva, 2002; Staneva and Stanev, 2002). According to Belokopitov (1998), extreme cold winters occurred in 1964, 1972, 1976, 1985 and 1993, which motivates consideration of the impact of the meteorological conditions on the biogeochemical state of the Black Sea. The seasonal signal, as seen in the simulated oxygen concentration (Fig. 7), reached 70 m depth on average. The penetration was deeper than 70 m in the middle of 1970se1986, and also in 1994. The isoline of 5 mM reached a maximum depth of 90 m in 1976. This time coincided with the turning point associated with the trends of air temperature (Fig. 4-1a), SST and CIL temperature (Fig. 4-2). Sulfide demonstrated an immediate response to the ventilation event in 1976 with 5 mM isolines sinking deeper than 100 m. A similar, but weaker response was observed in 1981. Analysis of simulations revealed a dramatic increase of TKE during these years, which made the deeper penetration of the surface signal possible.

Here we describe the model performance at the depths above and around the depth interval where suboxic layer is usually observed (Fig. 6). In the numerical simulations, the position of the suboxic zone was at the density layer of 15.2e15.7 kg/m3. The vertical distribution of H2S was similar in summer and winter. Large variations were observed in the distributions of oxygen in surface layers mostly due to temperature changes. A maximum oxygen concentration of about 380 mM occurred at the sea surface, similar to observations (Yakushev et al., 2007). Oxygen reached zero concentration (less than 0.5 mM) at st 15.2 kg/m3 that were a little shallower than in observations (15.3e16.0 kg/m3, Murray et al., 2005). Vertical distribution of nitrogen compounds (NH4, NO2 and NO3) showed that the onset of NH4 appeared at about st 15.4 kg/m3. The vertical gradient of Mn2 increased greatly where the onset of H2S occurred. The maximum of NO3 appeared approximately at the depth where the O2 dropped below 10 mM. This might serve as an indication that NO3 was an important oxidant in redox layer, which was below the depth where O2 was depleted. This result agrees with similar distributions described by Yakushev et al. (2006). The other oxidants were NO2. The depth of redox layer was between the depth of 65 m and 90 m. In the model, the oxidation of Mn2 resulted in the formation of alternative electron acceptors, and had a sinking rate to accelerate the downward transport. O2 dropped to 0 at the same depth where Mn2 increased that replicates well the observations. Furthermore, Yakushev and Debolskaya (2000) suggested that the reduction of Mn4 by H2S was very intense

98

Y. He et al. / Marine Environmental Research 77 (2012) 90e102

Fig. 6. Winter-mean (DecembereFebruary) in 1960 of the vertical profiles of chemical variables in depth coordinates (a, b) and density coordinates (c, d).

It is noteworthy that the winter of 1976 was not the coldest winter during the 45 years of analysis, and the winter-mean minimum of air temperature occurred after 1990 (Fig. 4-1a). However, penetration of O2 deeper than 90 m in 1976 was a consequence of the low temperature and extreme wind speed at the same time (see Fig. 3-2). More detailed analysis of forcing data revealed that the minimum of 2 m AT was w5  C lower than climatic mean values; maximum wind speed was about 5 m/s higher that the climatic norm. These events were not easily to detect in the NAO index and monthly mean climatic values. However, the ocean response reflected this climate variability originating from the atmosphere, revealing large variations in temperature and salinity at CIL with colder and fresher water penetrating down to 95 m in 1976 (Fig. 8). An unexpected feature accompanying the two events (in 1976 and 1981) was associated with a short-term rapid increase of H2S

below 110e120 m. Analysis of simulations revealed that this “perturbation” of the system was due to sinking particulate mater and enhanced bacterial respiration in anoxic waters. These changes resulted in longer-term evolution, with different trends in surface and deep layers. To explain these situations, three stages of oxidation of H2S with O2 (where S0 and S2O3 are intermediate chemicals) between 80 and 100 m were identified by Volkov (1974).

2H2 S þ O2 /2S0 þ 2H2 O

(13)

þ 2S0 þ O2 þ 2H2 O/S2 O2 3 þ 2H

(14)

 2 S2 O2 3 þ 2O2 þ 2OH /2SO4 þ H2 O

(15)

Y. He et al. / Marine Environmental Research 77 (2012) 90e102

99

It has been proven by Yakushev et al. (2006) that the dominant role in the oxidation of sulfide is due to the reduction of MnO2 and Mn(III), 50.2% and 42.6%, respectively. A considerable number of intermediate chemicals such as S0 and S2O3 also play an important role. Due to the large sinking rate of 0.5 md1 for S0 (Yakushev et al., 2006), they were rapidly delivered to deeper layers where: þ S0 þ H2 O þ H2 S/S2 O2 3 þ 2H :

(18)

The OM ((CH2O)106(NH3)16H3PO4) was divided in ROLM into DOM and POM with different rates of mineralization with different electron acceptors in oxic and anoxic conditions. In oxic conditions (e.g. above 90 m), OM was oxidized by O2 with the release of phosphate and ammonia. Below 110 m (e.g. anoxic conditions), sulfate reduction occurred with the production of phosphate and sulfide:

OM þ 106O2 /106CO2 þ 16NH3 þ H3 PO4 þ 106H2 O

Fig. 7. Simulated yearly distribution of the concentration of oxygen and hydrogen sulfide between 55 m and 150 m during of 1960e2000. Color scale shows the concentration of oxygen; the black lines give the concentration of hydrogen sulfide. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Furthermore, Mn species at the redox interface appeared to be of utmost importance. In particular, MnO2 and Mn(III) reached the maximum at suboxic layer and played an active role in the biogeochemical balance at about 90 to 105 m. At these levels sulfide oxidation proceeds as follows:

MnO2 þ HS þ 3Hþ /S0 þ Mn2þ þ 2H2 O;

(16)

2Mn3þ þ HS /S0 þ 2Mn2þ þ Hþ :

(17)

(19)

OM þ 53H2 SO4 /106CO2 þ H3 PO4 þ 16NH3 þ 53H2 S þ 106H2 O (20) This explains the H2S peak below 110 m in 1976. To now there has been little evidence of such a phenomenon, and further observations with accompanying 3-D model analyses would be useful. Some encouraging preliminary results from analyses of Argo floats in the Black Sea with oxygen sensors give promising results in this direction (Stanev et al, paper in preparation). 3.3. Long-term variability of the Black Sea hydrochemical conditions and the NAO index Numerical simulations gave a clear indication of long-term changes in the biogeochemical state of the Black Sea. Because the only external forcing in the model comes from the atmosphere and because there was a clear response of the physical state to the NAO,

Fig. 8. Simulated yearly mean variations of temperature (a) and salinity (b) at 95 m in 1960e2000.

100

Y. He et al. / Marine Environmental Research 77 (2012) 90e102

Fig. 9. Simulated interannual variability in the anomalies of oxygen, phytoplankton and nitrate concentrations in winter simulated at the sea surface (0e5 m) (left) and the CIL (60e90 m) (right) vs. the winter NAO index (DecembereFebruary). a. Oxygen at the surface; b. Oxygen at the CIL; c. Phytoplankton at the surface; d. Phytoplankton at the CIL; e. Nitrate at the surface; f. Nitrate at the CIL. The shaded is the interval 1973e1982.

we compare below yearly variability of oxygen, phytoplankton and nitrate concentrations at both the surface and the CIL with the NAO index. As shown in Fig. 9, anomalies of the major biogeochemical variables both at the sea surface and at the depth of the CIL were consistent with the change in the NAO index. The main changes in SST and CIL temperature had an opposite (as compared to NAO index) phase from 1960 to 2000 (see Fig. 4-2) with correlation coefficients of 0.84 and 0.59, correspondingly (Fig. 10a and b). The solubility of oxygen at the sea surface depended mainly on SST, such that the surface oxygen concentration and SST were

significantly correlated (R2 ¼ 0.99, Fig. 10c). At deeper levels it was not directly the solubility that dominated O2 concentration, but oxygen flux from above due to ventilation. Thus, the correlation coefficient between oxygen anomaly and temperature at CIL was only 0.34 (Fig. 10d). From 1973 to 1982, the oxygen anomaly at the surface remained relatively stable when phytoplankton and nitrate showed maxima around 1980. Furthermore, the correlation between the general trend of the NAO index and variability of both phytoplankton and nitrate anomalies at the surface were negative (Fig. 9c and e). The phytoplankton concentration at the surface showed a clear correlation with surface nitrate anomaly (R2 ¼ 0.70, Fig. 10e), that was much better than the relationship between phytoplankton and oxygen anomaly at the surface. The correlation between phytoplankton and oxygen anomaly at the CIL (R2 ¼ 0.88, Fig. 10f), was better than with the nitrate anomaly. This may indicate that nitrate was the limiting factor for phytoplankton at the surface, and oxygen played more important role for phytoplankton growth than nitrate at the depths of the CIL. It is noteworthy that the anomalies of oxygen, phytoplankton and nitrate at the CIL reached their maxima during 1973e1982. A similar trend was revealed for surface phytoplankton and nitrate. Consistent with the above analysis, we conclude that the NAO index, especially during times of intensity shifts, could be regarded as a potential long-term forcing explaining the decal variability of biogeochemical processes. 4. Discussion and conclusions

Fig. 10. Scatter plots revealing the correlation of annual anomalies of SST and NAO index (a), CIL temperature and the NAO index (b), oxygen at surface and SST (c), oxygen and temperature at the CIL (d), phytoplankton and nitrate at the surface (e), phytoplankton and oxygen at the CIL (f). Plots correspond to winters from 1962 to 2000. Correlation coefficients are also shown. NAO index is from NOAA.

We demonstrated in this paper that the coupled numerical model simulations revealed a pronounced interannual variability of state variables in the Black Sea, which showed that changes in the NAO index during the past four decades could be considered as an important driver for the variations of biogeochemistry in the Black Sea. Frequent abrupt shifts of the atmospheric state during 1975e1985 were of utmost importance for this period. The oxygen

Y. He et al. / Marine Environmental Research 77 (2012) 90e102

content at the sea surface revealed a stable increasing trend until 1996, demonstrating clearly the negative correlation of surface oxygen with SST. A maximum oxygen anomaly at the CIL was simulated in 1985, which also revealed a negative correlation with the CIL temperature. Our explanation for this finding is that biogeochemical responses to physical forcing are not always simple (linear). We demonstrated that the phytoplankton and nitrate anomalies at the surface had a maximum at around 1979, which was caused by the strong vertical mixing resulting in a larger supply of nutrients in 1976. Furthermore, phytoplankton revealed a good correlation with nitrate at the surface, while at the depth of the CIL it showed a better correlation with oxygen. This gives an indication that the limitation of phytoplankton growth changed in different layers. Our findings can be summarized as follows: (1) The NAO index had an overall increasing trend, except for the abrupt transition during the decade from 1973 to 1982. Consistent with that, the AT showed a decreasing trend for the whole period. However, a temporal increase was observed between 1973 and 1982 when wind speed had maximum, correspondingly minimum magnitude. Associated with these changes (shifts) in the atmospheric forcing, a pronounced minimum of SST during winter 1976 occurred. The trend of SST and temperature at the CIL were similar to that of the AT. The turning points of the trends for both variables occurred during 1973e1982. (2) Oxygen fluxes from above ventilated the redox layer periodically with the seasonal signal reaching w70 m. However, extreme events are more important for the biogeochemistry of the redox zone. The largest event was identified in 1976, resulting in a deepening of redox layer by more than 20 m. This event was a consequence of extremely low AT and high wind magnitude. At the same time the sulfide zone was displaced deeper to about the lower boundary of redox layer. (3) Below 110 m, a depth not reached by convection, a peak of H2S was identified which was almost synchronous (slightly delayed) with the timing of the surface signal. However, this cannot be explained as a direct consequence of the colder and fresher water intrusions into deeper water because the ventilation from above brings oxygen-rich water which oxidizes H2S. Assuming a direct response to the atmospheric forcing would lead to the opposite result, that is, to a decrease of the H2S concentration. Analysis on the biogeochemical variables simulated in the model showed that direct ventilation due to the extreme climate perturbation did not dominate the temporal variability of the entire biogeochemical system (including layers below the suboxic zone). Instead, it was the enhancement of the fluxes of particulate matter associated with the upper layer changes which enabled large amounts of Mn in its different forms to sink rapidly below the oxic layer, resulting in the simulated increase of the H2S concentration there. In the present study the model simulations were carried out in a strongly idealized model environment forced only by atmosphere, without considering river runoff. We also did not consider anthropogenic impacts. This could explain why simulated variability patterns did not always follow observations (Konovalov and Murray, 2002; Oguz, 2005). In this regard, very pronounced responses to anthropogenic forcing were studied by Lancelot et al. (2002) and Gregoire et al. (2008) who demonstrated that eutrophication increased gradually beginning in the 1960s. A strong decreasing trend started in the 1990s due to economic changes in the region. Consideration of anthropogenic and natural variations in isolation as well as their collective effects is an important area for

101

future study. However, such an approach would need a more realistic model representation. Another direction of future development in this field has been proposed in the framework of long-term 3-D numerical modeling (Gregoire and Beckers, 2004). This is mandatory because horizontal transport processes are expected to also play an important role. Furthermore, considering that most processes in the Black Sea are controlled by wind (Stanev et al., 2003) one could expect that the spatial variability of the atmospheric forcing will also have a pronounced impact in a 3-D modeling framework. Work is in progress on 3-d simulations using ROLM coupled with the General Estuarine Transport Model (GETM), one part of which is GOTM (He et al., paper in preparation). Preliminary results demonstrated that proper resolution of spatial dynamics was critical for Black Sea biogeochemistry. One additional perspective is model validation against data from Argo floats with oxygen sensors (Stanev et al., paper in preparation) and using Argo floats with new sensors measuring other biogeochemical variables. Using new data will enable better calibration of model parameters; a problem not yet fully resolved. There are preliminary indications that one of the major novelties of the present research e the increase of H2S synchronous with strong ventilation events e seems to be replicated by the recent observations with Argo floats. If this effect is proven by future observations, more focus in the future studies should be placed on the control of Black Sea biogeochemistry by sinking particulate matter. Acknowledgments We thank O. Podymov for his useful suggestions on the model coupling, J. Su for his comments on the discussion part and R. Kandilarov for his proframing help. Careful reading and useful comments of A. Dale are highly appreciated. Thanks are also due to two anonymous reviewers who for their comments. Atmospheric model data were produced by the European Center for MediumRange Weather Forecasts (ECMWF) and were kindly made available for the period 1958e2002 in the framework of the SESAME project by the INGV, Italy. The study was funded by project of HYPOX (No. EC Grant 226213 (gs1)). References Belokopitov, V., 1998. Long-term variability of cold intermediate layer renewal conditions in the Black Sea. NATO TU-Black Sea Project Ecosystem Modelling as a Management Tool for the Black Sea, Symposium on Scientific Results. NATO ASI Series 2/47 (2), 47e52. Burchard, H., Bolding, K., Villareal, M.R., 1999. GOTM, A General Ocean Turbulence Model. Theory, Applications and Test Cases, European Commission Report EUR 18745 EN: 103 pp. Burchard, H., Deleersnijder, E., Meister, A., 2005. Application of modified Patankar schemes to stiff biogeochemical models for the water column. Ocean Dynamics 55, 326e337. Gregoire, M., Beckers, J.M., 2004. Modeling the nitrogen cycle in an enclosed environment (the Black Sea): transport versus biogeochemical processes and exchanges across the shelf break. Biogeosciences 1 (1), 30e61. Gregoire, M., Raick, C., Soetaert, K., 2008. Numerical modeling of the deep Black Sea ecosystem functioning during the eutrophication phase. Progress in Oceanography 76 (3), 286e333. Gregoire, M., Soetaert, K., 2010. Carbon, nitrogen, oxygen and sulfide budgets in the Black Sea: a biogeochemical model of the whole water column coupling the oxic and anoxic parts. Ecological Modelling 221, 2287e2301. Jones, P.D., Jonsson, T., Wheeler, D., 1997. Extension to the North Atlantic oscillation using early instrumental pressure observations from Gibraltar and South-West Iceland. International Journal of Climatology 17, 1433e1450. Konovalov, S.K., Murray, J.W., 2002. Variations in the chemistry of the Black Sea on a time scale of decades (1960e1995). Journal of Marine Systems 31, 217e243. Konovalov, S.K., Ivanov, L.I., Samodurov, A.S., 2001. Fluxes and budget of sulphide and ammonia in the Black Sea anoxic layer. Journal of Marine Systems 31 (1e3), 203e216. Lancelot, C.L., Staneva, J.V., Van Eeckhout, D., Beckers, J.-M., Stanev, E.V., 2002. Modelling the Danube-influenced North-western continental shelf of the Black

102

Y. He et al. / Marine Environmental Research 77 (2012) 90e102

Sea. II: Ecosystem response to changes in nutrient delivery by the Danube River after its damming in 1972. Estuarine, Coastal and Shelf Science 54, 473e499. Marshall, J., Kushnir, Y., Battisti, D., Chang, P., Hurrell, J., McCartney, M., Visbeck, M., 1997. A North Atlantic climate variability: phenomena, impacts and mechanisms. International Journal of Climatology 21 (15), 1863e1898. Murray, J.W., Stewart, K., Kassakian, S., Krynytzky, M., Dijulio, D., 2005. Oxic, Suboxic and Anoxic Conditions in the Black Sea. Black Sea Overview. Murray, J.W., Stewart, K., Kassakian, S., Krynytzky, M., DiJulio, D., 2007. Oxic, suboxic, and anoxic conditions in the Black Sea. In: Yanko-Hombach, V., Gilbert, A.S., Panin, N., Dolukhanov, P.M. (Eds.), The Black Sea Flood Question: Changes in Coastline, Climate and Human Settlement. Springer, Dordrecht, pp. 1e22. Oguz, T., Gilbert, D., 2007. Abrupt transitions of the top-down controlled Black Sea pelagic ecosystem during 1960e2000: evidence for regime-shifts under strong fishery and nutrient enrichment modulated by climate-induced variations. Deep-Sea Research I. Oguz, T., Ducklow, H., Malanotte-Rizzoli, P., Murray, J.W., 1998. Simulations of the Black Sea pelagic ecosystem by one-dimensional, vertically resolved, physicalbiochemical models. Fisheries Oceanography 7 (3/4), 300e304. Oguz, T., Dippner, J.W., Kaymaz, Z., 2006. Climatic regulation of the Black Sea hydrometeorological and ecological properties at interannual-to-decadal time scales. Journal of Marine Systems 60, 235e254. Oguz, T., 2005. Black Sea ecosystem response to climatic teleconnections. Oceanography 18 (2), 122e133. Rodi, W., 1987. Examples of calculation methods for flow and mixing in stratified fluids. Journal of Geophysical Research 92, 5305e5328. Roussenov, V., Stanev, E.V., Artale, V., Pinardi, N., 1995. A seasonal model of the Mediterranean Sea general circulation. Journal of Geophysical Research 100 (C7), 13515e13538.

Spencer, D.W., Brewer, P.G., 1971. Vertical advection, diffusion and redox potentials as controls on the distribution of manganese and other trace metals dissolved in waters of the Black Sea. Journal of Geophysical Research 76, 5877e5892. Stanev, E.V., Peneva, E.L., 2002. Regional sea level response to global climatic change: Black Sea examples. Global and Planetary Change 32, 33e47. Stanev, E.V., Bowman, M.J., Peneva, E.L., Staneva, J.V., 2003. Control of Black Sea intermediate water mass formation by dynamics and topography: comparisons of numerical simulations, survey and satellite data. Journal of Marine Research 61, 59e99. Stanev, E.V., Staneva, J.V., Bullister, J.L., Murray, J.W., 2004. Ventilation of the Black Sea pycnocline. Parameterization of convection, numerical simulations and validations against observed chlorofluorocarbon data. Deep-Sea Research 51/12, 2137e2169. Staneva, J.V., Stanev, E.V., 2002. Water mass formation in the Black Sea during 1991e1995. Journal of Marine System 32, 199e218. Tsimplis, M.N., Josey, S.A., Rixen, M., Stanev, E.V., 2004. On the forcing of sea level in the Black Sea. Journal of Geophysical Research 109, C08015. Volkov, I.I., 1974. Geokhimiya sery v osadkah okeana. (Geochemistry of Sulfur in the Ocean Sediments). Nauka, Moscow, 272 pp. Yakushev, E.V., Debolskaya, E.I., 2000. Particulate Manganese as a Main Factor of Oxidation of Hydrogen Sulfide in Redox Zone of the Black Sea. IOC Workshop Report No. 159. Kluwer Academic Publishers, 592e597 pp. Yakushev, E.V., Chasovnikov, V.K., Debolskaya, E.I., Egorov, A.V., Makkaveev, P.N., Pakhomova, S.V., Podymov, O.I., Yakubenko, V.G., 2006. The northeastern Black Sea redox zone: hydrochemical structure and its temporal variability. Deep Sea Research Part II: Topical Studies in Oceanography 53, 1769e1786. Yakushev, E.V., Pollehne, F., Jost, G., Kuznetso, I., Schneider, B., Urnlauf, L., 2007. Analysis of the water column oxic/anoxic interface in the Black and Baltic seas with a numerical model. Marine Chemistry 107, 388e410.