Mixing in the Black Sea detected from the temporal and spatial variability of oxygen and sulfide – Argo float observations and numerical modelling

The temporal and spatial variability of the upper ocean hydrochemistry in the Black Sea is analysed using data originating from profiling floats with oxygen sensors and carried out with a coupled three-dimensional circulationbiogeochemical model including 24 biochemical state variables. Major focus is on the dynamics of suboxic zone which is the interface separating oxygenated and anoxic waters. The scatter of oxygen data seen when plotted in density coordinates is larger than those for temperature, salinity and passive tracers. This scatter is indicative of vigorous biogeochemical reactions in the suboxic zone, which acts as a boundary layer or internal sink for oxygen. This internal sink affects the mixing patterns of oxygen compared to the ones of conservative tracers. Two different regimes of ventilation of pycnocline were clearly identified: a gyre-dominated (cyclonic) regime in winter and a coastal boundary layer (anticyclonic eddy)-dominated regime in summer. These contrasting states are characterized by very different pathways of oxygen intrusions along the isopycnals and vertical oxygen conveyor belt organized in multiple-layered cells formed in each gyre. The contribution of the three-dimensional modelling to the understanding of the Black Sea hydro-chemistry, and in particular the coast-to-open-sea mixing, is also demonstrated. Evidence is given that the formation of oxic waters and of cold intermediate waters, although triggered by the same physical process, each follow a different evolution. The difference in the depths of the temperature minimum and the oxygen maximum indicates that the variability of oxygen is not only just a response to physical forcing and changes in the surface conditions, but undergoes its own evolution.


Introduction
The Black Sea (Fig. 1) is the world's largest marine anoxic basin (Skopintsev, 1975;Murray et al., 1995Murray et al., , 1989Sorokin, 1982;Jørgensen et al., 1991). Anoxic conditions were formed about 8000 years ago because of the reconnection of Mediterranean and Black Sea and the intrusion of saltier Mediterranean water which followed the reconnection (Deuser, 1974). The micro-structure profiling measurements (Gregg and Yakushev, 2005) supported earlier estimates from analysis of tracers (Lewis and Landing, 1991) and numerical simulations (Stanev et al., 1997) revealing extremely low vertical turbulent exchange (coefficients of about 1-4 × 10 −6 m 2 s −1 ), approaching the values of molecular exchange. This explains the formation of a strong vertical stratification (Fig. 2) and the Cold Intermediate Layer (CIL), a permanent layer at 50-100 m (see the blue line in Fig. 2a). This layer with low temperature in its core acts as a near-surface thermal reservoir, similar to the North Atlantic Subtropical Mode Water. Another consequence of the strong stratification is the decoupling of the Black Sea in oxygenated surface layer and a sulfide-containing deep layer (Oguz et al., 2005). The oxic and sulfidic waters are separated by a layer between the isopycnal depths σ t = 15.40 and σ t = 16.20 known as the "suboxic zone" (Murray et al., 2005). This zone was defined by Murray et al. (1989Murray et al. ( , 1995 as the water layer, where the concentrations of dissolved oxygen and hydrogen sulfide are below the detection limit. In the study of Konovalov et al. (2003) these limits have been specified as 30 nmol L −1 for sulfide and 3 µmol L −1 for oxygen. As demonstrated by Stanev et al. (2013), Argo floats' observations in the Black Sea equipped with an optode sensor can resolve these low oxygen concentrations. For a recent analysis of the spatial and temporal variability of chemical properties in the Black Sea based on historic data, see Tugrul et al. (2014).
The suboxic zone is an important biogeochemical interface between the surface and deep waters. The processes controlling its origin and variability have been extensively discussed by Saydam et al. (1993), Basturk et al. (1994) and Yakushev et al. (1997Yakushev et al. ( , 2007. Murray et al. (1989) suggested that the suboxic zone might be a new feature caused by the reduced fresh water input from rivers and the resulting change in the ventilation, but Buesseler et al. (1994) claimed that it was most likely a permanent feature in the Black Sea. Later, Konovalov and Murray (2001) demonstrated that the depth of the upper boundary of suboxic zone was governed by the balance between the ventilation of thermocline with oxygen-rich surface water and oxygen consumption by the oxidation of organic matter.
The sensitivity of the suboxic zone upon the physical drivers was revisited by the numerical experiments of Yakushev et al. (1997Yakushev et al. ( , 2009) who demonstrated the role of the cycling of nitrogen-manganese-iron-sulfur elements. Within the oxygen-deficient part of the water column de-nitrification results in a sharp decrease of the nitrate concentration ( Fig. 2c) from the peak to trace values and zero concentrations around 110 m. Below the suboxic zone hydrogen sulfide and ammonium are found in high concentration filling the anoxic zone down to the bottom (Fig. 2b, c).
A series of complicated redox processes takes place in this transition zone -the manganese cycling is one of the most important ones. This process was first presented correctly in the classic paper by Spenser and Brewer (1971) and was studied later by Murray et al. (1995), Rozanov (1995) and Yakushev et al. (2007). The manganese cycling, in which reduced Mn can be oxidized only by oxygen, while oxidized forms of Mn can be reduced by sulfide and iron was first incorporated in a biogeochemical model by Yakushev (1998) and later implemented in the model of Debolskaya and Yakushev (2002) accounting for the enhanced sinking rate of Mn particles.
The temporal and spatial variability of the suboxic zone reflects fundamental changes associated with the sulfide production in the water column and sulfide oxidation at the base of the suboxic zone (Neretin et al., 2001), which are shaped by circulation and mixing. The circulation in the Black Sea is structured usually in two connected gyre systems encompassing the basin (the main current is known as the Rim Current) and two quasi-permanent sub-basin eddies: the Batumi and Sebastopol ones. The Rim Current, which is a relatively narrow jet current, is associated with a difference of ∼ 0.2 m between sea level in the coastal and open sea. Between the jet current and the coast a number of coastal anticyclonic eddies occur. The horizontal transport almost doubles in winter, which is the season of more intense circulation. During this season, the circulation usually occurs as a one-gyre system in the entire basin. The baroclinic eddy dynamics, and in particular the balance between vorticity in the open sea and coastal sea, control the transition between summer and winter circulation ; in the warm part of the year the anticyclonic sub-basin-scale eddies are more pronounced (Stanev and Staneva, 2000; for a short overview of Black Sea dynamics, see Stanev, 2005).
While the temporal and spatial variability of circulation is relatively well known, this is not the case for the oxygen because of insufficient amount of appropriate observations. The present paper aims to fill this gap and to demonstrate that continuous observations over large areas and coupled 3-D numerical modelling can advance our understanding of the temporal and spatial variability of oxygen in the Black Sea. The amalgamation of observations and numerical modelling is considered thus as the basic methodological concept of the present research. The reason to put the oxygen in the middle of our research focus is that we want to use in a coherent way high-resolution in time and space observations in parallel with numerical simulations. The recent observational activities using Argo floats equipped with oxygen sensors (Stanev et al., 2013) provided such data.
Although the Black Sea has been used as the major playground to test oxygen models and hypoxia dynamics (Pena et al., 2010), there are still a number of open questions, including (1) what is the affordable, with available computational power, vertical and horizontal resolution which would guarantee a reasonable model-representation in the suboxic layer, and (2) what is the consistence of the numerical simulations with continuous observations? The answers to (1) and (2) could open discussion of the needed sampling capabilities of available observational platforms and design of appropriate models. The present study takes a step in this direction. However, there is one very basic question to be answered first: which are the important features in the dynamics of oxygenated and sulfide water, which cannot be understood fully using available in the past data and simple concepts? One such 1-D concept is the assumption that tracers in the Black Sea well follow the isopycnals that is the vertical stratification in density coordinates (one profile only) adequately represents the stratification in the entire basin. While this could apply to passive tracers and reflects the fact that the vertical displacement of pycnocline would result in similar displacement of water properties (see Fig. 16b of Stanev et al., 2004), it is not fully applicable to the non-conservative tracers (oxygen is one of them). The basic hypothesis of the present study is that the diapycnic and isopycnic mixing (that is the mixing across and along isopycnic surfaces) are closely linked with different dynamical states of the Black Sea circulation and largely control the patterns of oxygen and sulfide. This hypothesis is consistent with the ideas of Gnanadesikan (1999) that the overturning circulation in the ocean is determined by the rates of isopycnal and diapycnal mixing. In the case of the Black Sea there are some specificities. The vertical diffusion coefficient is very small, which enables one to speculate that the diapycnal diffusion is very low. However, Konovalov et al. (2003) reported massive lateral injections of oxygen-enriched waters of the Bosphorus plume into the suboxic and anoxic layers. This could serve as an indication that the isopycnic diffusion could play the major role in propagating tracers' anomalies. It will be shown in the present study that this process is not only important for the region around the Bosphorus Strait, but dominates the oxygen exchange between the coastal and open sea. If the isopycnic diffusion were very strong it would smear tracers' anomalies and the along-isopycnal gradients would be small. This is known from the previous observations of temperature and  (Stanev et al., 2004) and is consistent with the dominant role of the isopycnic diffusion.
The available observations from profiling stations could identify the importance of isopycnic diffusion for oxygen, which is not a conservative tracer. Unlike the temperature, oxygen does not change the vertical stratification, however it is subject to internal sources and sinks associated with biogeochemistry reactions. One layer of active reactions is in the suboxic zone, which acts as a very narrow chemistry boundary layer from where the isopycnic diffusion propagates the oxygen anomalies. Thus unlike in the case of temperature, salinity and passive tracers, the biogeochemical reactions play an important role at these depths. They determine the shape of the chemistry boundary layer, thus ensuring a coupling between physical and biochemical processes. Finally, because the gradients of chemical tracers cannot change the stability of stratification (they do not change the physical mixing), large anomalies of oxygen can occur on isopycnic surfaces.
The demonstration of the role of diapycnic and isopycnic mixing for the dynamics of non-conservative tracers (in this case oxygen and sulfide) and the examination of the individual aspects of the hypothesis formulated above will be achieved by addressing the following specific research objectives: (1) to describe the coupling of the 1-D Redox Layer Model (ROLM, Yakushev et al., 2007) with a 3-D numerical circulation model (General Estuarine Ocean Model, GETM), and (2) to describe the temporal and spatial variability of oxygen as seen in profiling floats observations and numerical simulations. The paper is organized as follows. Section 2 presents a brief analysis of the Black Sea hydrochemistry. The used physical and biogeochemical models are described in Sect. 3 and Appendices A-D. Section 4 describes the model validation against observations followed by analyses of the major results in Sect. 5 and short conclusions.

Diapycnic and isopycnic mixing -general considerations and definitions
The T -S analysis developed by descriptive oceanographers has given one illustration of how the pathways of oceanic motion can be identified from observed tracer fields (Wüst, 1935). Later, the physical understanding of water mass formation was substantially developed following the pioneering work of Walin (1982) who demonstrated that buoyancy has to be supplied to allow water to cross isopycnal surfaces. Soon after, McDougall (1984) proposed a solid theoretical (mathematical) formulation of the problem. His "water mass transformation" equation represents a balance between downgradient advection on an isopycnal and diffusion through the along-isopycnal and vertical mixing. Recently it has been realized that the rates of isopycnal and diapycnal mixing con-trol the overturning circulation (Gnanadesikan 1999;Sijp et al., 2006). This fundamental result indicates that in order to describe the vertical overturning correctly one has to well understand both isopycnal and diapycnal mixing. As shown by McDougall (1984) and Zika et al. (2009) along-isopycnal mixing depends upon the ratio of alongisopycnal curvature to along-isopycnal gradient of tracer field and smooths out the curvature of this field. Thus gradients and curvature of tracers observed on isopycnals can identify isopycnal mixing. The vertical mixing depends upon the curvature of the tracer with respect to density. This is one of the most useful characteristics of ocean mixing displayed usually by the T -S diagrams, and this will be extensively used in the present study. As indicated by McDougall (1984), if the T -S curve obtained from a vertical CTD cast is locally straight, the diapycnal eddy diffusion would not alter temperature on this isopycnal surface. This theoretical finding is another expression of an axiom of oceanography since at least Sverdrup, Johnson and Fleming (1942) that two conservative properties will plot along a linear relationship.
The above short synthesis on the ocean isopycnal and diapycnal mixing applies to temperature and salinity, which are conservative tracers. For the Black Sea case, Murray et al. (1991) argued that all waters deeper than the CIL are formed by two end-member mixing of the CIL and Bosphorus outflow. However, the major interest in this paper is on the non-conservative tracers oxygen and sulfide. Let us assume that oxygen and salinity dynamics is described as a balance between vertical diffusion and transport (Munk, 1966) with the addition of a reaction term for oxygen: where O 2 , S, w, ν t and R are concentration of oxygen, salinity, vertical velocity, coefficient of turbulent diffusion and reaction term, respectively. We admit that the above model does not fully describe all processes of interest. It is used just to facilitate the understanding of further results and can in principle be extended by substituting the terms in the lefthand side by other physically more consistent terms. Knowing the profiles of S and O 2 from observations and using the above equations we can estimate the reaction term as The term in brackets is proportional to the curvature of the oxygen with respect to salinity d 2 O 2 ∂S 2 =  McDougall (1984) and Zika et al. (2009) and is proportional to the diffusion in S-space. For the suboxic zone the S-space is almost the same as the density space because the role of temperature anomalies there is minor. In this way the above simple theory unifies the oxygen consumption and the diapycnic (diahaline) diffusion. Therefore the question is worth asking, what do we learn about mixing and biogeochemical reactions from the distribution of S and O 2 in the Black Sea? It is worth recalling in this context that in estuarine physics (bearing in mind that the Black Sea can be considered as a big estuary), mixing plots, in which a dissolved constituent is plotted against salinity, are commonly used to interpret conservative and non-conservative processes (Liss, 1976;Loder and Reichard, 1981). A bend in the mixing curve is indicative of the reactive or nonconservative nature of tracers.

The survey
Two Argo floats were deployed on 7 May 2010 in the northern Black Sea (Stanev et al., 2013). The floats known as the Navigating European Marine Observer (NEMO) profilers were equipped with temperature, salinity and oxygen sensors. The sensor for temperature and salinity was CTD SBE 41, and that for oxygen was Aanderaa Oxygen Optode 3830. The trajectories of the two floats (NEMO-0144 and NEMO-0145) are shown in Fig. 1. Their initial positions were close to each other, but they departed one from another soon after deployment. The float NEMO-0144 entered the central basin and ceased operations on 28 December 2012; the other one, NEMO-0145 moved along the western and southern coasts and continued operations until 27 November 2012. Both floats were programmed to sample the water column from the sea surface down to 500 m and transmit the data every 5 days via satellite. This deployment enabled for the first time more than 2 yr long continuous observations in the Black Sea. The fact that the two floats operated in two dynamically very different areas -the central gyre and the continental slope area (Fig. 1) -gives the unique opportunity to compare synchronously the differences between changes of hydrochemistry in these different zones. In what follows, we will refer for brevity to these two floats as the central gyre and jet current floats, respectively, having in mind that the first has spent most of its lifetime in the ocean interior and the second has mostly sampled the Rim Current area (see Fig. 1).

Overall presentation of observations -the suboxic zone as a boundary layer
In the following the vertical stratification in the Black Sea will be illustrated from the profiles measured by NEMO-0144 during the time it operated in the central gyre. This choice was made because we want to illustrate the most general characteristics of tracers' distribution, and this is better seen in the basin interior. The observations show a pronounced three-layer structure. In the upper layer the high oxygen concentrations extend down to about σ t = 14 (Fig. 3a). This isopycnic depth, above which the temperature decreases almost linearly with the increasing density ( Fig. 3b), marks approximately the upper boundary of CIL (Fig. 3b). A rapid, overall linear, increase of oxygen is observed in the second layer between σ t = 14.2 and σ t = 15.2. The decrease of oxygen concentrations in the third layer (below σ t = 15.2) represents the transition to anoxic conditions and is the lower boundary of the Black Sea aerobic zone. The layer where a linear regression between oxygen and density occurs is actually the upper halocline. The base of this layer is where the isothermal part of temperature versus density relationship starts (Fig. 3c).
The slope of the regression line describing the approximate relationship between hydrogen sulfide and density in the upper part of the anoxic zone (below σ t = 16.5) is similar to the respective oxygen versus density regression line (between σ t = 14.2 and σ t = 15.2). Thus the oxygen and sulfide profiles can be considered as a unique O 2 -H 2 S curve (Fig. 3a). It resembles the variable D = aO 2 -H 2 S introduced by Stanev (1989), where a is a stoichiometric constant. The O 2 -H 2 S curve undergoes a jump between σ t = 15.4 and σ t = 16.5 which is the transition zone between the layers where oxygen is the dominant oxidizer, and sulfide is the dominant reducer (redox zone). Figure 3a nicely illustrates that the suboxic zone has the typical characteristics of a boundary layer and that this boundary layer is easily detectable in density space. The straight regression curves above and below the suboxic zone are also reminiscent of the classical theories (Sverdrup et al., 1942;McDougall, 1984) demonstrating that for conservative tracers (temperature and salinity) the straight T -S curve ( Fig. 3c) means that diapycnal eddy diffusion does not alter temperature locally.
It follows from the simple theory presented in Sect. 2.1 (see Eq. 3) that the "jump" in the O 2 -H 2 S curve is intimately associated with the reaction term and that the latter is proportional to the diapycnal mixing. Therefore this curvature reflects the chemical transformations of matter in the redox layer where there are different oxidizers such as nitrate, Mn(IV), Fe(III), reducers such as NH 4 , Mn(II), Fe(II) and intermediate species acting as oxidizers or reducers (NO 2 , Mn(III), S 2 O 3 , S 0 ). In this simple theoretical example (Eq. 3) the diapycnal mixing is balanced by the biogeochemical reactions. The curvature of profiles below σ t = 15.2 in the suboxic layer indicates strong changes of the balance between physical processes and consumption (or matter transformation). Recall that the diapycnic mixing in the Black Sea is small; however, there is also a mixing along isopycnals which can close the balance. Thus the interplay between physical and biogeochemical processes should create specific mixing patterns of non-conservative tracers in the suboxic zone. We will try to identify in the following these patterns from the continuous Argo observations.

Evidence of diapycnic and isopycnic mixing
The difference between the mean vertical oxygen stratification in the central gyre and in the jet current region observed by the two floats reveals sharper stratification in the central gyre and deeper penetration of oxic water in the jet current region (Fig. 4). Since during the periods of intense mixing (especially in winter), the surface density value increases to ∼ 13 (outcropping of isopycnal surfaces) whilst during the stratified period the surface density values decrease to ∼ 11, the oxygen profiles under mixed and stratified conditions are analysed separately. We will refer below to the corresponding Vertical profiles of the standard deviation of oxygen from the mean (for each profile) values. For each Argo float we have separated the profiles into two sets: the red profile corresponds to stratified conditions while the blue one corresponds to mixed situation (this explains why the blue profiles which correspond to winter conditions have a surface density of 13 instead of 11 for the red ones). We made this distinction between the stratified and mixed conditions in order to average similar types of profiles, so as to ensure working with homogeneous data sets. profiles as "mixed" and "stratified", respectively (see the legend in Fig. 4). Below the CIL the mixed and stratified parts of the profiles (red and blue lines) display approximately the same mean stratification because the seasonal variability almost does not reach these depths.
The standard deviation from the mean profiles ( Fig. 4b) gives an idea about how stable the above mentioned linear regression is. In the upper layers (down to σ t = 13) the standard deviation reaches about 20-40 µM during the stratified period that is about one-tenth of the mean value. This large signal versus noise ratio, expressed as the ratio of the mean over the standard deviation, renders a statistical significance to the mean profiles ( Fig. 4a) in the upper layer. Below σ t = 13 in winter and σ t = 14 in summer the standard deviation continuously increases for the two floats reaching a maximum of about 40 µM for NEMO-0144 at σ t = 14 and above 60 µM for NEMO-0145 at σ t = 14.6. These maxima occur during winter which is associated with the intense circulation during this season. The magnitude of the standard deviation and the signal-to-noise ratio demonstrate that (1) the variations are pronounced down to σ t = 16, and (2) a statistical significance of the mean oxygen stratification down to the upper boundary of the suboxic zone. Below this depth, where the curvature of vertical profiles tends to increase (at σ t = 15.1 and σ t = 15.5 for the central gyre and jet current area, respectively), the mean oxygen value are comparable with the magnitude of the standard deviation. From this result it follows that the mean profiles in the suboxic zone are not enough representative for the individual locations sampled by the two floats.
The curvature of oxygen profiles in the jet current area appears smaller than in the central gyre (compare the two curves in Fig. 4a). The pronounced differences between the second derivatives of the mean oxygen curves of central gyre and jet current areas in the suboxic layer manifests the possible strong coupling between the biogeochemistry processes and turbulent exchange (see Eq. 3) at the transition between two zones.
Small gradients of oxygen on isopycnic surfaces between the jet current and the central gyre would indicate a large contribution of the isopycnic diffusion, as is the case with passive tracers (Stanev et al., 2004). However, the large scatter in individual profiles ( Fig. 3a and Fig. 4b), as well the large difference between profiles in the central gyre and jet current area, reaching values of more than 50 µM in the quasi-linear part of the mean profiles, manifest the substantial difference between passive and non-conservative tracers.
As a conclusion we can say that the strong decrease of oxygen with increasing density expressed by the regression line in Fig. 3a gives only an overall idea about the oxygen stratification. The large differences between profiles of NEMO-0144 and NEMO-0145 are indicative of the rigorous biogeochemistry processes in the area between central gyre and jet current. It appears that in the suboxic zone (below σ t = 15) the contribution of biogeochemistry processes is comparable to that of seasonal and mesoscale variability, and thus the direct use of what is known for passive tracers becomes questionable. The large curvature of profiles in Fig. 4a at these depths demonstrates that the transition to anoxic conditions in the presence of low diapycnic diffusion should be accompanied with substantial isopycnal mixing. Later we will give an answer of the question where this occurs.

The numerical model
The General Estuarine Transport Model (GETM, Burchard and Bolding, 2002) solves the equations for the three velocity components u, v and w, temperature T , salinity S and sea surface height ζ . The applicability of this model to the problems addressed here is justified by the generic formulation of turbulent mixing using General Ocean Turbulence Model (GOTM; Burchard et al.,1999). GOTM solves the equations for turbulent kinetic energy (TKE) k and the dissipation rate ε: where υ ε and υ k are the turbulent diffusivities of energy dissipation and TKE, respectively. The shear stress production P and the buoyancy production B are where N 2 = −g ρ 0 ∂ρ ∂z is the Brunt-Väisälä frequency. For the turbulent diffusivities of temperature, salinity and momentum the relation of Kolmogorov and Prandtl is used. The turbulent diffusivities for k and ε are All parameters in the above formulas are the same as in He et al. (2012); see also (Burchard et al., 1999) and Appendix A for the performance of GOTM. The numerics of GETM are described by Burchard and Bolding (2002); its implementation to the Black Sea is described by Peneva and Stips (2005). Because the present study addresses the complex redox dynamic including the suboxic zone it is of utmost importance to select an appropriate model. One available and well-documented model is the ROLM (RedOx Layer Model) which is a model for the dynamics of the oxic-anoxic interface zone.
The redox processes implemented in this model were first modelled by Yakushev and Neretin (1997) with a focus on the nitrogen and sulfur cycles. In their model the sulfur cycle is driven by a continuous supply of oxygen provided by a downward diffusive oxygen flux. Debolskaya and Yakushev (2002) advanced the model by incorporating a simplified manganese cycling in which the particulate manganese is used for oxidizing the hydrogen sulfide, and the dissolved manganese is oxidized by oxygen. Yakushev et al. (2007) added in ROLM a parametrization of the processes of formation of organic matter during both photosynthesis and chemosynthesis.
ROLM considers specific important consumers of oxygen. The reduced and intermediate species of metals (Mn(II), Mn(III), Fe(II)) and sulfur (S 2 O 3 and S 0 ) are important sinks of O 2 in the suboxic zone in addition to O 2 consumption for mineralization of organic matter and nitrification. If these processes are not included, the model would produce increased values of O 2 at these depths, and therefore unrealistically deeper positions of the oxygen surfaces. Furthermore, Yakushev et al. (2007) demonstrated that under low oxygen conditions (O 2 < 30 µM) 50 % of O 2 was consumed by the processes connected with oxidation of Mn(II), Mn(III), Fe(II), S 2 O 3 and S 0 . The parametrized complex redox dynamics is also needed for the oxygen modelling above the suboxic zone. We admit that the ecosystem parametrization in ROLM does not include all details of the ecosystem functioning, but nevertheless it allows us to parametrize the main features of the dissolved oxygen fate in the upper layer.
In its 1-D version ROLM uses the following set of equations for non-conservative substances C i : where C i is the concentration of "ith" model variable, A C V is the vertical turbulent diffusion coefficient for biogeochemical tracers, which is assumed to be equal to the vertical turbulent diffusion coefficient for temperature and salinity ν t , W C is the sinking velocity of particulate matter, and W Mn is the accelerated velocity of sinking of particles with Mn hydroxides. The term R C i = j Rate B j C i is the production-destruction of C i caused by biogeochemical interactions. Rate B j C i describes the rate of biogeochemical production and consumption of C i , by B j . The parametrization of biogeochemical processes was described by Yakushev et al. (2007).
In the present paper the 1-D coupled Black Sea GOTM-ROLM (He et al., 2012) is implemented in GETM. Because the biogeochemical part used here is essentially the same as in the study above, only some details concerning the simulation of turbulence, vertical resolution and model performance are presented in Appendices A and B. From these results it follows that the simulated diffusion coefficients are in agreement with the ones based on analysis of observed data. Furthermore, a very fine vertical resolution of about 2 m is needed to realistically simulate the biogeochemical process in the upper 150 m because the vertical gradients (e.g. in the suboxic zone) are extremely sharp. Keeping a resolution of 2 m in the 3-D model that is about 1000 vertical levels in the deep part of the sea, which is around 2000 m deep, is a big computational challenge. Due to this reason a simplified approach to the modelling of the Black Sea hydrochemistry is proposed in the following: (1) focus on the upper ocean down to 200 m, which enables us to address with sufficient vertical resolution the dynamics of suboxic zone and, (2)

Upper ocean patterns
The simulated February-mean SST in 2010 (Fig. 6a) is lower than 8 • C over the western and central Black Sea reaching 10 • C around the eastern coast. This low SST enables cold intermediate water to form over large areas (for the water mass formation in the Black Sea see Ovchinnikov and Popov, 1987;Stanev et al., 2003). The lowest SST occurs in the northwestern shelf; this low-temperature area extends southward because of the transport along the western coasts. The surface oxygen pattern is almost opposite to the distribution of SST because low SST enhances the dilution of oxygen in the surface layer. This pattern is not shown here because it is very similar to that at 20 m (see below). The oxygen gradients in the upper mixed layer are much stronger than the ones of SST and more affected by the hydrodynamics, as seen along the eastern Black Sea coast by the tonguelike coast-guided pattern in Fig. 6c. The difference between the surface and saturated oxygen concentration displays the pattern of air-sea oxygen flux (Fig. 6b). This field reveals signatures of jet current and meanders, which is observed in the area west of the Crimea Peninsula.
The highest oxygen values are simulated in the northwestern shelf area (Fig. 6c), which is explained by the minimum in the air temperature in this region. Similarly, the eastern coast guides the propagation of the low-oxygen water. With the increasing depth the distribution of O 2 does not change down to 20 m, however at 40 m the patterns are rather different (Fig. 6d). At this depth, the oxygen values in the eastern Black Sea drop to about 200 µM, which is due to the vertical circulation bringing low-oxygen water into the surface layer. The above analysis demonstrates that spatial variability, which is missing in the 1-D models (He et al., 2012), could substantially impact the oxygen dynamics.

The thermohaline fields in the intermediate layers
The numerical simulations compare well with the ones reported by Stanev et al. (2003) and Stanev and Kandilarov (2012), which used the same horizontal grid, as well as with the GETM modelling of the Black Sea circulation of Peneva and Stips (2005). This is the reason why the simulated hydrodynamics will not be extensively addressed here, leaving space for more details on the modelling of hydrochemistry. The circulation is predominantly cyclonic and structured in persistent western and eastern gyre; the jet current is stronger in winter than in summer (e.g. Staneva et al., 2001). Some dominant signatures in the circulation patterns are described in more detail further in the text when the impact of circulation on the biogeochemistry is presented. Among them is the strong meander west of the Crimea Peninsula, which could be considered as the model analogue of Crimea Eddy. The second important sub-basin-scale eddy is the Batumi eddy, which is better pronounced in summer.
The CIL in the Black Sea (Fig. 7) is a permanent layer at 50-100 m (see also Figs. 2a and 3) acting as a near-surface thermal reservoir, similar to the North Atlantic Subtropical Mode Water, which is renewed due to vertical mixing during cold winters (Stanev et al., 2003;Gregg and Yakushev, 2005). Winter cooling penetrates down to about 100 m in the western Black Sea (Fig. 7) and the cold water area extends over the entire basin. Only in the easternmost regions is the SST higher than the temperature in the deeper layers. In summer the cold water is overlain by the warm surface water forming the CIL (Fig. 7b). It is noticeable until the next winter, giving a clear evidence of how persistent this water mass is. The variations in the depth of its lower boundary correlates with the oscillations of the depth of the pycnocline (compare the upper and bottom panels in Fig. 7). Its thickness in the central basin is smaller than in the continental slope area reflecting the impact of upwelling and downwelling, respectively. The temperature in the core of the CIL is in most cases higher than 8 • C, which does not support the usually accepted definition of Blatov et al. (1984) describing the CIL as a permanent layer at 50-100 m with temperatures lower than 8 • C. However, this gives a modelling support to the recent evidence (Stanev et al., 2013) that for the considered period the cold intermediate water was substantially warmer than previously known.
At the base of the CIL the main pycnocline which is mostly due to the strong vertical salinity gradients, acts as an obstacle for the propagation of temperature signal into the deeper layers. Therefore its lower boundary follows the topography of pycnocline.
The variability patterns in the continental slope area and the central gyre reveal substantial deviations reaching 2 • C in the seasonal thermocline. The pattern of the subsurface temperature anomalies demonstrates that the differences are rather a consequence of the dynamics and not of the local atmospheric forcing. Overall, the jet current zone is cooler than the central gyre (Fig. 7b).

The mesoscale dynamics as seen in the data from Argo profiling floats
The mesoscale oxygen dynamics in the Black Sea were previously analysed by Stanev et al. (2013). In order to better demonstrate some features associated with the isopycnic mixing we present here the oxygen, temperature and salinity data from the two floats in density coordinates (Fig. 9). Although the whole data series from the two floats is shown in the figure, only the data which are representative of central gyre are commented below for NEMO-0144 -that is the data which are measured after 6 July 2010 when the float en-tered the deep sea. Before entering the central gyre this float showed similar characteristics as the jet current flow. The missing data during the warm part of the years reflects the outcropping of the isopycnic surfaces. The rest of the missing data (white strips) is due to errors in the observations or reaching the bottom. Both floats reveal an increase of surface oxygen concentrations during the periods of winter convection -the latter periods are better identified in the temperature diagrams (Fig. 8c,  d). With increasing the time during the re-stratification season this oxygen maximum displaces to lower isopycnic layers. In depth coordinates it occurs at about 25 m (Stanev et al., 2013). The suboxic zone identified in Fig. 8a, b from the position of the 5 µM oxygen surface is at about σ t = 15 and σ t = 16 for NEMO-0144 and NEMO-0145, correspondingly.
The physical conditions dominating the temporal and spatial change of oxygen are better understood from the plots of temperature, salinity and depth of isopycnic surfaces (Fig. 8c-h). The CIL (dark blue coloured layer in Fig. 8c, d) is clearly observed in the records of both floats. The isolines in Fig. 8 c, d display well the formation and evolution of the CIL: isopycnic surfaces of σ t = 14 outcrop and lowtemperature water is injected into the deeper isopycnic layers. This process is observed around the southern coast where NEMO-0145 spent part of its lifetime, but also in the central gyre. The second demonstrates clearly that the CIL is formed not only in the coastal zone but also in the central gyre, which supports the observations of Ovchinnikov and Popov (1987) and numerical simulations of Stanev et al. (2003). The CIL is refilled every year in February-March. Its thickness in density coordinates is almost equal in the central gyre and coastal area (compare Fig. 8c and d during 2010). It is noteworthy that after the relatively warm winter in 2011 the core of this layer weakened, and this was recorded by the two floats. As explained by Stanev et al. (2013) who used data from all available floats in the Black Sea, the temperature of the CIL was warmer than 8 • C not only during 2011 but during most of the last decade. This suggests revisiting the accepted definition of this layer as an intermediate water mass with temperature lower than 8 • C (Blatov et al. 1984). However the extremely cold winter of 2012 resulted in a strong refill, which was recorded also in the zone around the southern coast. This situation, characterized by a thick (in density space) layer of water colder than 8 • C, persisted along the southern coast until December 2012. It is thus instructive to use the comparison between the cold water content in 2010 and 2012 (compare the left-and rightmost parts in Fig. 8d) as an indication of the approximate range of the interannual change in the thermal conditions. The temporal variability of the Black Sea stratification in density coordinates differs largely from the analysis in depth coordinates described by Stanev et al. (2003). They showed that (1) the variations of the depth of the lower boundary of CIL were very pronounced in the continental slope area and much weaker in the central gyre and (2) the CIL in the central gyre was thinner than in the continental slope area. These features reflect the impact of general upwelling and downwelling, respectively. They also demonstrate that the distribution of temperature correlates with the oscillations of the depth of the pycnocline. The isopycnic analysis does not reveal this because Fig. 8a-f represent the actual mixing patterns, which do not always correlate with the oscillations of pycnocline (Fig. 8g, h). Salinity (Fig. 8e, f) changes very little and in a smooth way below σ t = 15.5 (NEMO-0144) and σ t 16.2 (NEMO-0145). The horizontal and straight isohalines manifest that below these depths the density is almost fully dependent upon salinity (Fig. 8e and Fig. 8f). However the density surfaces undergo substantial vertical displacement, which is  (Fig. 8g, f). It follows from the above analysis that temperature and salinity tend to be homogenized along isopycnic surfaces. In the presence of a very small diapycnic mixing this becomes possible due to the isopycnic mixing, which results in very low gradients along isopycnals, something which is consistent with the earlier knowledge in this field.
The above conclusions do not fully apply to oxygen, and this is better seen in the observations from the jet current float NEMO-0145 revealing vigorous change of oxygen along the isopycnals. In comparison, the central gyre float shows lower amplitudes of oxygen at the same surfaces. The oxygen isolines 5 µM and 50 µM are found at a distance of about 0.5 σ t in the interior part of the gyre and at larger isopycnic distance in the jet current area, exceeding during some intervals 1 σ t . This difference is possibly associated with the large magnitude of mixing at the periphery of coastal anticyclones .
The results presented above manifest different pattern of isopycnic mixing of non-conservative tracers (from those of temperature), which is associated with the biochemistry processes and mesoscale dynamical oscillations. The oxygen isoline of 5 µM penetrated down to 170 m in February 2012 (Stanev et al., 2013) and reached σ t = 16.8. During this time NEMO-0145 float was between the Sakarya Canyon and Synop, an area representative for the coastal jet dynamics. The short duration of this oxygen change and the correlation with similar changes in temperature and salinity at that time made the coastal eddies the most plausible candidates to explain this event. However, there are many oscillations of oxygen ( Fig. 8b) with almost the same amplitudes as the abovementioned event, most of which do not clearly correlate with the oscillations of pycnocline. The opposite is also true -not all strong oscillations of the pycnocline (Fig. 8h) resulted in a clear change of oxygen. As a conclusion one can say that the oscillations of oxygen have relatively larger amplitudes and are much "noisier" than those of temperature. The absence of the "counterpart event" in the halocline depth reveals that the hydrochemical fields do not always follow the evolution of the physical system.
There are also some pronounced regularities in the oxygen variability. While the isoline of 50 µM was recorded at almost the same isopycnic levels in the central gyre and in the jet current area, the 5 µM isoline was deeper in the latter area. This demonstrates that the suboxic zone fills wider density space in the jet current area, however during short periods it shrinks substantially.

Comparison between the numerical simulations and observations
The observations and 3-D simulations follow the temporal variability in the vertical profiles known from the 1-D simulations (see Fig. 7 of He et al., 2012). In both data sets the mean oxygen values in the continental slope area are higher than in the central gyre. The observed standard deviation of oxygen from the mean oxygen profiles (NEMO-0144 and NEMO-0145) and the numerically simulated standard deviation displayed in Fig. 9 demonstrate that the simulated pattern of temporal variability replicate relatively well the observed one; the magnitudes of the maxima are also similar. The oxygen stratification is stronger than in the observations; the shallower penetration of oxygen signal in the continental slope area in comparison to the observations can be explained by the lower intensity of coastal circulation simulated in the model. This calls for using in the future simulations finer horizontal resolution (full eddy resolution) and not only eddypermitting resolution as used here. The temporal variability of the upper layer oxygen can be well illustrated with the help of the oxygen versus temperature relationships (Fig. 10). These profiles are characterized by almost vertical curves in the late fall and beginning of winter -that is with a relatively homogeneous oxygen distribution in the cold season. The lowest oxygen values at sea surface approach in summer 250 µM both in the model and observations. The warming of surface waters is accompanied with a decrease in the surface-oxygen concentration and the formation of oxygen subsurface maximum (see also Fig. 8a,  b). One could expect that this oxygen maximum is the analogue of the CIL (e.g. oxygen-rich water overlain in summer by less oxygenated water). However, the core of this layer is at shallower depth than the core of the CIL (Stanev et al., 2013, compare also Fig. 8a, b and Fig. 8c, d). Furthermore, a closer analysis of observations demonstrates that some summer values in the core of the oxygen maximum are higher than the ones in winter. This excludes the possibility that subsurface oxygen maximum is just a direct consequence of the oxygen-rich water created by winter convection, which remains overlain in summer by the low-oxygen (because of the high SST) surface water. The observed absolute oxygen maximum in the subsurface layers is consistent with the findings of Yakushev et al. (2007), Konovalov et al. (2006) and Gregoire and Soetart (2010)  oxygen in the Black Sea was due to phytoplankton, which shows larger production in summer. The conclusion from the above comparison between numerical simulations and observations is that the model is able to replicate almost all observed features of oxygen dynamics. This qualifies the numerical simulations as a useful tool to investigate the upper-ocean hydrochemistry of the Black Sea and its relationship with the dominant physical processes.

Seasonal variability
Although the 3-D oxygen dynamics in the Black Sea has been widely addressed by Grégoire et al. (1998), Grégoire and Lacroix (2001) and Oguz and Salihoglu (2000), the associated diapycnic and isopycnic mixing have still not been sufficiently described. This is a fundamental problem as indicated by the major novelty seen in profiling float observations, which is that the temporal and spatial variability of oxygen does not fully follow the one of physical variables identifying active biogeochemical interactions (or sinks). This result is supported by the numerical simulations, which will be presented in the following. The comparison between the 1-D model simulations and the basin mean 3-D simulations (Fig. 11) demonstrates some important differences: (1) the oxygenated upper layer in the 1-D simulations is thinner than in the 3-D simulations, (2) the vertical gradients are stronger in the 1-D simulations; (3) the suboxic zone is deeper and thicker in the 3-D simulations. The overall result of the comparison between the simulations from the 1-D and 3-D models is that the basinmean vertical column in the former shows weaker variations. Although the same boundary conditions are used in the two models, the oxic zone is deeper in the 3-D simulations, demonstrating clearly the role of circulation for establishing the vertical stratification of hydrochemical fields. Recall that the relatively smooth behaviour of the seasonal signal in Fig. 11b is rather a consequence of the basin averaging. The comparison between the fields in the central gyre and the jet current zone (see further in the text) reveals that the winter cooling reaches larger depths along the basin periphery and the anoxic layer is strongly eroded there by the sinking oxygen-rich water. In the 1-D model (without spatial dynamics) the temporal variability of anoxic waters remains low. It becomes obvious that the most unrealistic drawback of the 1-D simulations is the very thin suboxic zone. This has been substantially corrected in the 3-D simulations and resulted in a better agreement of the numerical simulations with the observations. Furthermore, the suboxic zone in the 3-D simulation is thinner in spring and thicker in summer, while its thickness almost does not change in the 1-D simulation (Fig. 11). These differences between the 1-D and 3-D simulations are interpreted as a consequence of the circulation.
The representation of the same signals in isopycnal coordinates gives a clearer idea about the major sources of the variability in the upper-ocean oxygen content. The low surface density values disappear in winter and an outcropping of density surfaces down to σ t = 14.2 are simulated on 15 February. This is accompanied by an injection of oxygenrich water into the pycnocline during winter, as seen in the basin-averaged conditions. The secondary subsurface maximum appearing above the major one (at about σ t = 14) replicates successfully a similar maximum observed by the profiling floats in the upper seasonal thermocline (see Fig. 8). Thus one important difference between the 1-D and 3-D simulations is that in the 3-D simulations the main subsurface oxygen maximum appears in a narrower σ t interval (Fig. 11). Furthermore, without the biological processes the core of the oxygen-rich surface water would constantly diffuse in time before the next convective event. This is actually the case for temperature, for which there are not internal sources. We admit that the inflow from the Bosphorus Strait tends to increase the temperature in deep layers, however this signal penetrates deeper than the processes considered here.

Analysis along vertical sections -the oxygen conveyor belt
In the upper mixed layer the simulated oxygen patterns in winter are characterized by the surface maximum and a very sharp oxygen stratification at 50 m marking the depth of the direct ventilation from the sea surface (Fig. 12a). The deepening of the oxygen isoline 5 µM, describing approximately the position of the upper boundary of suboxic zone, down to about 100 m in winter along the western and eastern coasts presents an important evidence of water ventilation between the jet current and the coastal zone. This surface remains almost horizontal at around 80 m (Fig. 12a, b) in the central gyre. In summer, the surface oxygen values drop to about 250 µM, which is consistent with the observations from profiling floats (see Figs. 8 and 10). The subsurface oxygen maximum covers in summer the entire deep part of the Black Sea with its core at about 30 m (Fig. 12b). While the area of coastal downwelling is well revealed in the simulations by the deepening of the oxygenated waters in the continental slope area of the western basin in winter (Fig. 12a), this situation is not always the case, as seen in the summer crosssection (Fig. 12b). The thickness of the suboxic zone, considered here as the space between the oxygen and sulfide isolines 5 µM, increases in summer. Its lower boundary deepens from 100 m in February to about 140 m in July. This change is associated with the decrease of density stratification in summer -the density isoline of 16 kg m 3 is at about 125 m in February and 155 m in July. The presentation of the Black Sea hydrochemical properties in density coordinates (Vinogradov and Nalbandov, 1990;Konovalov et al., 2005) reduces the data scatter seen in geopotential coordinates, and renders analyses less dependent on the specific location (Stanev et al., 2004). Analyses in density coordinate have also been extensively used to evaluate the variability of physical and biogeochemical fields in 1-D models (Oguz et al., 1996(Oguz et al., , 2000(Oguz et al., , 2001(Oguz et al., , 2002aYakushev and Neretin, 1997;Yakushev et al., 2007). However, under the assumption that simple tracer-density relationships are valid basin-wide (e.g. the regression curve in Fig. 3a), the vertical sections in density coordinates (Fig. 12c, d) would show horizontal (in density space) isolines of oxygen and sulfide. This would mean the same oxygen value on a constant density surface. Obviously, this is not the case in the simulations, giving an idea about the role of isopycnic mixing. The isolines of both variables (O 2 and H 2 S) almost mirror the density plots in Fig. 7c, d which could illustrate that in the areas of shallower depth of isopycnic surfaces the oxygen concentrations would be larger because the shallower position is closer to the sea surface where the oxygen content is larger. In contrast, the deeper position of isopycnic surface would tend to decrease the concentration of oxygen because in the deeper layers the consumption of oxygen would be larger. This presents an illustration of the possible coupling between physical and biogeochemical mechanisms in creating boundary layers for non-conservative tracers.
The agreement between oxygen distribution and density (compare Fig. 12c and Fig. 7c) is particularly well seen in winter above the suboxic zone (σ t = 14.5 − 15). In summer the along-isopycnal contrasts are smaller, which reflects the lower intensity of the summer circulation. During this season it is the lower boundary of the suboxic zone (H 2 S = 5 µM) which mirrors better the shape of density surfaces (in particular for the levels below σ t = 14.5) than the oxygen isolines in the layers above do.
The above results motivated us to analyse the differences between the stratification in particular time and along particular section and the basin mean stratification (Fig. 13). This difference would "remove" the mean situation and make clearer the contribution of the 3-D dynamics. The "anomaly" representation makes the origin of the surface oxygen-rich water (largest positive anomalies in winter) very clear. Not quite trivial is the negative anomaly originating from the sea surface, which gives an indication of localized convection. In these "chimneys" (e.g. Fig. 13a) the vertical mixing results in a consumption of oxygen by sulfide (through the other elements intermediates), leading to lower oxygen values. This serves as a modelling indication that the signatures of pronounced mixing between oxic and anoxic waters triggered by winter convection could reach the surface layer. Such indications can be found in the observations shown Fig. 8a, b, in particular at the end of 2010.
The anomaly-plots in winter reveal several thinner zones which are localized around the areas where the pycnocline is in its highest position. In these areas the upward motions change the stratification with respect to the mean one resulting in the formation of multi-layer vertical structure of oxygen anomaly. In summer, the oxygen anomaly at the same section is represented by the following: (1) a negative anomaly in the coastal zone (Fig. 13c, d), (2) oxygen anomaly following the main pycnocline represented by the thin blue strip, and (3) thin secondary maxima and minima in the upper 50 m layer. Obviously, the horizontal anomalies from the mean (Fig. 13a, c) correlate well with the shape of the pycnocline (Fig. 7c, d). The diapycnal mixing in the area of multiple layers acts as the major agent controlling the matter exchange. This is justified by the presentation of the same cross-sections in isopycnic coordinates (Fig. 13e-h). The intermediate layers of positive anomaly at about σ t = 15 (Fig. 13e, f) provide the major source for isopycnal mixing.
From the upwelling zones, the oxygen-rich waters propagate into the direction of downwelling zones -that is into the coastal zone, as well as into the basin interior the latter separating the western and eastern gyres. As a conclusion of the above discussion, which is mostly focused on the winter situation, one could summarize that the outcropping appears to be the major mechanism pumping oxygen-rich water into the pycnocline. The upwelling gives the second mechanism to increase the oxygen concentrations. Double cells are formed (one in each gyre). Propagating approximately along the domed isopycnal surfaces these waters tend to distribute basin-wide the positive oxygen anomalies.
In summer the oxygen anomalies are overall positive in the eastern basin where a two-layer structure occurs (Fig. 13g). The numerical simulations give a clear indication of lowoxygen water along the north-south section (Fig. 13h) manifesting a large consumption of oxygen over relatively large areas in the western basin. The continental slope area acts as a guide for the diapycnal mixing in summer, forming one branch of the vertical oxygen conveyor belt. Furthermore, a number of local anomalies are seen in the area of CIL and below it. They give an evidence of diapycnic mixing in the basin interior and the resulting intrusions which follow the neutral surfaces.

Analysis along horizontal sections
The horizontal distribution of oxygen at 50 m (Fig. 14a, b) reveals an important seasonal change consisting in the following: while the coastal maximum in February is observed mostly in the area of the northwestern shelf (Fig. 14a), the summer maximum is in the interior part of the western basin ( Fig. 14b) and follows along the southern and eastern costs. The pattern of the oxygen-rich waters at 50 m in winter is explained by the extreme cooling of the western basin and the resulting ventilation. In the same area, in summer, the lowoxygen waters propagate along the continental slope starting from the Crimea Peninsula and continuing along the western and part of the southern coasts (Fig. 14b). The western gyre gives the origin of the oxygen maximum in summer, which reaches the area of Kerch Strait following the coast (Fig. 14b).
The role of the sub-basin circulation features is better seen at 80 m (Fig. 14c, d). The winter pattern indicates an increased concentration to the southeast of the Crimea Peninsula (Fig. 14c), which is caused by the downward propagation of water in this area. The downward motions explain the increased oxygen concentrations in most of the coastal area (Fig. 14c), which is well seen around the western and southern coasts as very thin high-oxygen strips. In the easternmost part of the sea, the anticyclonic circulation associated with the Batumi eddy (Fig. 14d) plays the role of the oxygen source for the deeper layers in summer.
The distribution of oxygen on selected isopycanl depths (σ t = 14.5 and σ t = 15.5) in February (Fig. 14e, g) and July (Fig. 14f, h) reveals the major pathways of penetration of oxygen into the pycnocline and its consumption. The summer and winter cases manifest different mechanisms of upper-ocean hydrochemistry. In winter, the outcropping in the central gyre enhances the penetration of oxygen from the surface into the deeper layers of pycnocline (Fig. 14e, g). In this season the direction of the isopycnal transport (from high to low values) is from the ocean interior towards the rim current. This is opposite to the direction of Ekman drift which is directed from the central gyre to the coast (Stanev et al., 2014). Several small-scale exceptions are noteworthy such as the region south of the Crimea Peninsula and along the western and southern coasts where the oxygen-rich waters at σ t = 14.5 propagates from the coast in the direction of Rim Current. This feature does not penetrate much deeper (e. g. down to σ t = 15.5) as seen in Fig. 14g. The conclusion from the discussion of the winter situation is that the oxygen patterns can be easily explained based on the known Black Sea dynamics. The upwelling in the central gyre is stronger in winter because of the more intense winter circulation. This results in a shallower position of the pycnocline, which acts as a preconditioning factor. Thus the fluxes of oxygen from the well oxygenated in this season surface layers reach deeper isopycnic levels.
In summer the sea surface could not be considered as an important source of oxygen because of the high temperatures and the resulting low dilution rates. Furthermore, the Black Sea circulation in this season is less intense, which changes the dynamical control. The oxygen patterns seen in isopycnic coordinates (Fig. 14f) again reveal a boundary of aligned layers with high gradients of oxygen along the southern and eastern coasts. The situation is completely different along the western continental slope (a relatively mild slope) which, at this isopycnic depth, acts as a sink of oxygen. Obviously, the numerical simulations reveal complex mixing patterns with jet current area acting either as a source or sink for the isopycnic mixing (recall that in the case of large isopycnic mixing the distribution of O 2 would be homogeneous on isopycnal surfaces). Furthermore, the direction of the along-isopycnic gradients (coast versus open sea) changes at σ t 15.5 (compare Fig. 14g and h), demonstrating how complicated the vertical structure of isopycnal mixing can appear.
The opposing oxygen patterns in winter and summer ( Fig. 14e and f) necessitate an explanation of the reason of the simulated differences. One step in this direction is to examine the correlation between the patterns in depth and density coordinates. Overall, the comparison between Fig. 14a and e (winter) reveals a negative correlation, while the comparison between Fig. 14b and f (summer) reveals a positive correlation. This would mean that it is not the intensity of the circulation in summer and the oxygen sources from above, which dominates the distribution of oxygen but rather the mixing bringing oxygen-rich water from the core of subsurface oxygen maximum layer into the pycnocline. The strengthening of the coastal anticyclones and the weakening of the general circulation in summer (Stanev and Staneva, 2000;Staneva et al., 2001) could be considered as the physical candidate to explain the above-described shifts of the oxygen dynamics between the two different states. Summarizing, while the winter state is dominated by the intense basin-wide circulation, outcropping and large oxygen sources at sea surface, the summer one is dominated by the mesoscale circulation and mixing. Thus two contrasting regimes control the isopycnal diffusion in different seasons: a gyre-driven (cyclonic) regime in winter and a coastal boundary layer-driven one in summer. The latter is associated with the enhancement of anticyclones between the jet current and the coast. The difference between the profiles from the two Argo floats in the jet current and central gyre demonstrates the possible coastal-to-open-sea oxygen exchange and will be used in the following to validate the simulations. The contrast in the pycnocline seen in the density coordinates (Fig. 15a) is about one-fifth of that in depth coordi-nates (Fig. 15b). This supports similar estimates based on the observed profiles of passive tracers (Stanev et al., 2004). Knowing that the oxygen transport along the isopycnals is of diffusive character the distribution of areas with high (in the jet current area) and low (in the central gyre) oxygen values would give an overall idea about its direction. Above σ t = 13 the oxygen difference is presented in Fig. 15a by the summer curves because only in this season is the density lower than 13. The gradients are strongest in the layer σ t = 14 − 16 (geometrical depths of 50-100 m, see Fig. 15b) with a core at about σ t = 15 (Fig. 15a). It is noteworthy that the maximum contrast between the two areas is located just above the layer of maximum salinity difference (Fig. 15c). The position of the latter as observed in 2010-2011 was confirmed by our independent analysis of salinity profiles (not shown here) measured by earlier Argo floats (Korotaev et al., 2006). These floats were not equipped with oxygen sensors, however they were operating for longer time. This finding indicates that salinity does not only control the depth of winter cooling by providing an extremely stable vertical stratification, but also guides the intrusions from the coastal zone. Evidence for that is given by the correlation between the oxygen and salinity contrasts during different seasons. There is a negative correlation between salinity and oxygen contrasts in the two studied areas: salinity is lower in the jet current zone and higher in the central gyre -the opposite is valid for oxygen (compare Fig. 15b and c).
One could hypothesize thus that the maximum coastal-toopen-sea contrasts could be associated either with (1) the dynamics of the frontal area (Rim Current), or (2) the coastal boundary layer and sources of salinity and oxygen there. This first assumption would lead to the expectation that the intensification of the circulation in winter (steeper frontal area) would support stronger coastal-to-open-sea gradients of oxygen in this season. Just the opposite is observed on the isopycnic surfaces -that is the correlation between salinity and oxygen contrasts does not change sign, but larger contrasts appear in summer than in winter (Fig. 15a). This evidence supports the second hypothesis: the pronounced seasonal variability of oxygen and salinity difference between the coastal and open sea water is rather associated with the anticyclonic circulation occurring between the Rim Current and the coast. This type of circulation supports the downwelling, which is stronger in summer (Stanev and Staneva, 2000;Staneva et al., 2001). The downwelling supports larger concentrations of oxygen and low-salinity coastal water during this season. The evidence of this source is well manifested by the numerical simulations (Fig. 14f) and applies not only to the Bosphorus Strait area, as known from previous studies (Konovalov, 2003), but also for the larger part of the jet current area.
It is noteworthy that the smallest coastal-to-open-sea contrast is at about σ t = 14, where the transition between the seasonal thermocline and the pycnocline occurs. This is also the layer where the CIL core is observed. Above this isopycnal depth the temperature and density are almost linearly correlated in the warm part of the year (Fig. 3b). In the layers below, salinity and density are linearly correlated because the water mass below σ t = 14.5 is almost isothermal. Therefore at these depths the oxygen stratification is controlled by the salinity gradients and the oxygen consumption, which explains the linear regression in Fig. 3a. At the depth of the suboxic interface the signature of mixing is enhanced because of the rigorous biogeochemical reactions making this layer act as boundary layer, something which is specific for oxygen and not for temperature.

Conclusions
In this study we analysed the temporal and spatial variability of oxygen and sulfide in the Black Sea using continuous oxygen observations from profiling floats and numerical simulations with an appropriate 3-D physical-biogeochemical coupled numerical model. The major focus of the analyses was put on the dynamics of suboxic zone. The validation of the numerical model demonstrates that it replicates the statistics seen in the observations. Theories about the water mass formation were briefly presented and their relevance to the Black Sea conditions addressed with a focus on the behaviour of non-conservative tracers. In particular, the differential properties (vertical gradient and curvature of oxygen profiles) were analysed in isopycnic coordinates which allowed us to elucidate their relevance to the isopycnic and diapycnic mixing. It was also demonstrated that the chemistry interface separating oxygenated and anoxic waters seen in density coordinates gave a clear evidence of different mixing controls in the upper mixed layer, in the upper part of the halocline (approximately where the CIL is observed) and around the suboxic zone. In the latter zone a clear evidence of coupling between the physical and chemical mechanisms was given, which contributed to the formation of a specific boundary layer of mixed physical-biogeochemical nature (non-conservative mixing).
Numerical simulations provided a useful tool supporting the understanding of the spatial and temporal variability of the suboxic zone. However, an extremely fine vertical resolution of about 2 m was needed to address adequately its dynamics. In order to keep the computational requirements for the 3-D simulations reasonable the model was set up only for the upper 200 m part of the water column using appropriate boundary conditions. Even under these limitations the numerical model is able to reproduce the main features of water mass circulation in the Black Sea -the cyclonic currents and upwelling in the interior gyre, as well as the coastal anticyclonic circulation. This realistic physical background is very basic for the realistic simulations of horizontal and vertical transport and mixing of chemical properties, as well as their temporal-spatial variability. The temporal variabil-ity of the upper and lower boundaries of the suboxic zone showed a clear seasonal character.
Further to the numerical modelling study of Stanev et al. (2004) dealing with passive tracers, it appeared that the biogeochemical properties are not only displaced up and down by the vertical motion. The simulations revealed that pronounced isopycnic mixing dominated not only the Black Sea thermodynamics at the depth of the CIL and below it, as known from earlier studies in the area of the Bosphorus Strait, but also its biogeochemical state presented here by the oxygen gradients on isopycnic surfaces in a wider transition area between the jet current and central gyre. One of the basic conclusions from this paper is that the hypothesis of property alignment to the isopycnal surfaces depicted only very coarsely the oxygen distribution (in particular in the suboxic zone) and that using a 3-D numerical model is the better choice to enhance the realism of simulations. It was also demonstrated that neglecting the mesoscale variability or events, such as convection of cold and oxygenated waters, when addressing the Black Sea hydrochemistry largely disregarded physical processes associated with the alongisopycnal mixing. This limits the applicability of simplified theories, e.g. 1-D interpretations of data, or 1-D modelling.
A new concept has been proposed about the control of dynamics on the Black Sea hydrochemistry based on the analysis of numerical simulations, which is well supported by observational data. This concept was manifested by the existence of an oxygen conveyor belt organized in multiple cells formed in each gyre, which was demonstrated from the analysis of vertical cross-sections presented in Sect. 6.2.1 both in geopotential and density coordinates. These cells tend to distribute basin-wide the oxygen anomalies between the oxygen sources and sinks. It was also demonstrated that these layered cells were largely guided by the continental slope because of the large mixing there. Dynamical controls appeared very different for summer and winter: a gyre transport dominance in winter versus mixing dominance in summer. These contrasting regimes are characterized by very different pathways of oxygen intrusions along the isopycnals (Sect. 6.2.2).
Evidence was also given that the formation of oxic waters and cold intermediate waters, which were triggered by the same physical process, followed different evolution. The first substantial difference between the CIL and the core of oxygen maximum was that the depth of the former was larger. This was explained by the different distributions of sources and sinks of heat and oxygen. The second difference was the subsurface maximum of oxygen in summer, which was of biological nature. This and some other results presented in this study indicate that the distribution of oxygen did not always follow that of temperature, salinity or passive tracers, manifesting the difference between conservative and nonconservative mixing. Its variability is not only a response to mesoscale variability in the circulation and seasonal changes in the surface conditions, but undergoes its own natural evolution.

E. V. Stanev et al.: Mixing in the Black Sea
The differences between the central gyre and jet current areas revealed strong contrasts enhancing the isopycnic mixing, which appeared quite diverse in the different layers, along the steep continental slope or over the regions where the bottom topography is characterized by milder changes. This extreme variety called for deeper analyses and more focused study on the temporal-spatial variability of the biogeochemical fields.
Although the contribution of the 3-D modelling to the understanding of the Black Sea hydro-chemistry, and in particular the coast-to-open-sea diapycnal mixing, appeared very clear, some further developments are also needed, as for instance using more realistic riverine fluxes, increasing the horizontal resolution (moving from eddy-permitting to eddyresolving coupled models), giving deeper consideration of light regime and sedimentation, and assessing the success of using simpler and computationally more efficient biogeochemistry models. Among the computational challenges is the extension of the model grid down to the bottom (about 2000 m), which could be facilitated by the use of adaptive grids in the vertical, using better numerics, or addressing in depth the interplay between the processes dominant on the shelf and in the upper ocean (e.g. hypoxia) and those studied in the present research (dynamics of suboxic zone).

Appendix B: 1-D coupled physical-biogeochemistry model
This appendix is motivated by the need to discuss the optimal vertical resolution of the model, which has to be consistent with the complex biochemical processes, in particular the ones taking place in the very thin (about 20-30 m) suboxic zone. Three experiments with different vertical resolutions are carried out using 1-D coupled General Ocean Turbulence Model (GOTM) and ROLM. The basic model run with a vertical resolution of 2 m (experiment A) is essentially the same as the one described by He et al. (2012). The maximum water depth is 200 m -that is, 100 model levels are used. Such resolution would be very computationally intensive if used in the 3-D model, therefore results of two sensitivity experiments (B and C) with vertical resolution of 5 and 20 m, correspondingly, are first discussed in order to demonstrate the resolution limits.
All experiments used the same initial and boundary conditions and the meteorological forcing is computed using atmospheric data from the ECWMF reanalysis (http://www. ecmwf.int) as has been explained by He et al. (2012). At the sea surface oxygen flux is computed as in Yakushev et al. (2007): where O 2 sat is the oxygen saturation concentration computed as a function of temperature and salinity (UNESCO, 1986), Sc is the Schmidt number for oxygen, and k 660 is the reference gas-exchange transfer velocity computed as k 660 = 0.365W 2 + 0.46W where W is the wind speed magnitude (Schneider et al., 2002). Fluxes of PO 4 and NO 3 (see Sect. 3.1 for the description of variables) are prescribed using the concept developed by Fonselius (1974) who showed that 6700 tons of phosphorus had to be added annually to the Black Sea to balance the deposited bottom phosphorus, ensuring thus that the phosphorus concentration remains at a constant level. The corresponding flux is Q PO 4 = 0.13 mmol m −2 d −1 . Using similar considerations, the flux of nutrients due to the rivers' input and atmospheric deposition is estimated as Q NO 3 = 1.5 mmol m −2 d −1 (Yakushev et al., 2007).
At The model has been integrated in time until the solution has reached a quasi-periodic state. Four model variables are chosen to analyse the effect of resolution (Fig. B1). At first glance the major difference between the simulations with different vertical resolution is the coarse resolution of the vertical profiles, in particular in the 20 m resolution case. Not less important, however, is that in the simulations with 2 m resolution (experiment A) O 2 and H 2 S do not overlap, which is consistent with the concept that the suboxic zone decouples the oxic and anoxic waters. In experiment B, oxygen and sulfide overlap slightly, whereas the coexistence of O 2 and H 2 S in experiment C at 55-85 m totally disagrees with the observations. The maximum of nitrate (6 µM at about 55 m) is sharper in experiment A, which agrees well with the observations presented by Yakushev et al. (2007). This layer is more "diffuse" in experiment B and badly resolved in experiment C. The present result demonstrates a similarity with the sensitivity experiments of Yakushev et al. (2009). However, while in the later work the abnormally high vertical diffusion is the reason for the disappearance of suboxic zone, the same effect in the present simulation is caused by the coarse vertical resolution. Recall that with the model parameters used (see He et al., 2012 and the references therein) the simulated suboxic layer is thinner than the one seen in the observations. However, it is not aimed here to improve the performance of the 1-D model because, as shown further in the paper, this drawback is rather a consequence of the oversimplifications in the 1-D frame and disappeared in the 3-D model.