Sensitivity of oxygen dynamics in the water column of the Baltic Sea to external forcing

A 1-D biogeochemical/physical model of marine systems has been applied to study the oxygen cycle in four stations of different sub-basins of the Baltic Sea, namely, in the Gotland Deep, Bornholm, Arkona and Fladen. The model consists of the biogeochemical model of Neumann et al. (2002) coupled with the 1-D General Ocean Turbulence Model (GOTM). The model has been forced with meteorological data from the ECMWF reanalysis project for the period 1998–2003, producing a six year hindcast which is validated with datasets from the Baltic Environmental Database (BED) for the same period. The vertical profiles of temperature and salinity are relaxed towards both profiles provided by 3-D simulations of General Estuarine Transport Model (GETM) and observed profiles from BED. Modifications in the parameterisation of the air-sea oxygen fluxes have led to a significant improvement of the model results in the surface and intermediate water layers. The largest mismatch with observations is found in simulating the oxygen dynamics in the Baltic Sea bottom waters. The model results demonstrate the good capability of the model to predict the time-evolution of the physical and biogeochemical variables at all different stations. Comparative analysis of the modelled oxygen concentrations with respect to observation data is performed to distinguish the relative importance of several factors on the seasonal, interannual and long-term variations of oxygen. It is found that natural physical factors, like the magnitude of the vertical turbulent mixing, wind speed and the variation of temperature and salinity fields are the major factors controlling the oxygen dynamics in the Baltic Sea. The influence of limiting nutrients is less pronounced, at least under the nutrient flux parameterisation assumed in the model. Correspondence to: S. Miladinova (svetla.miladinova@jrc.ec.europa.eu)


Introduction
The Baltic Sea is a semi-enclosed and brackish sea, which combined with specific physical as well as socio-economic characteristics makes it very sensitive to anthropogenic pressures (Bonsdorff et al., 2001).Eutrophication remains the most pressing problem in the region, as nitrogen and phosphorous inputs are still high, despite considerable efforts to reduce discharges.The inflow of salty and oxygen rich water from the North Sea through the Danish straits is due to episodic pulses (Omstedt et al., 2004).The strong pulses are driven by special atmospheric forcing conditions, which cause large and long-lasting sea level differences between the Kattegat and the Western Baltic.Since the early 1980s, the Baltic Sea has experienced long-lasting stagnation periods with absence of strong pulses (Matthäus and Nausch, 2003).Only in 1993 and 2003 such major inflows took place (Jakobsen, 1995;Feistel et al., 2003b).Inflows from the North Sea are currently the main source of oxygen in the deep water.The deepwater basins in the Baltic Proper suffer severely from long-term oxygen depletion.Oxygen deficiency has prevailed over very large areas.In the central Baltic Proper the oxygen concentrations are less than 90 [mmol O 2 /m 3 ] at around a depth of 100 m, or even shallower than that (HEL-COM, 2003).At the same time, the area covered by hydrogen sulphide extends from the main eastern Basin of the Gotland Sea towards the Northern Central Basin (Fig. 1).Typically in August, oxygen is depleted in the bottom water of the Bornholm Basin and the western Gotland Basin.In the Arkona Basin the oxygen situation is good in the near-bottom water, although the level is lower when compared to the longterm measurements.
Additional to the above mentioned horizontal advection of oxygen the most important natural physical factors affecting the concentrations of oxygen in the marine environment are temperature and salinity.Oxygen concentrations decrease with increasing temperature and salinity (Quinlan, Published by Copernicus Publications on behalf of the European Geosciences Union.1980).The other major factor controlling oxygen concentrations is the biological activity in the water and at the seafloor: photosynthesis producing oxygen and respiration and nitrification consuming oxygen.The present study is motivated by the need to explore the processes that are controlling the oxygen dynamics in the Baltic Sea.For this purpose, the 1-D ecosystem model of Burchard et al. (2006) is used.Marine ecosystem models which involve the interaction of physical and biogeochemical processes are useful tools for assessing and predicting the trends in oxygen variation and for identifying the areas susceptible to oxygen deficiency.These models must take into account the most important biogeochemical processes and the physical control of the ecosystem driven by advection and diffusion.Efficient models of marine systems can simulate the seasonal evolution, inter-annual variability and spatial heterogeneity across the range of coastal and eutrophic situations with little or without re-parameterisation.Although the usual way to develop such models is to couple circulation models with biological models, simplified model systems based on 1-D water column models (e.g.those of Burchard et al., 2006;Kühn and Radach, 1997;Blackford et al., 2004) are very helpful tools for model development.
Depending on the scientific question, they can be also reliable in studying marine ecosystem dynamics of coastal marine areas.The validity of a 1-D approximation in the Baltic Proper is confirmed also by other model studies (Vichi et al., 2004;Omstedt and Axell, 1998;Stigebrandt, 1987).They are mainly related to periods when advection is negligible (so-called stagnant periods).Despite that a 1-D model exhibits limitations in simulating seasonal and interannual variability of the deep water mixing and the formation of density currents (Axell, 1998), it is a good tool for basic studies, aimed at improving the model parameterisation and for the investigation of some system properties.We have identified a number of questions that need explicit attention: -How well are the air-sea fluxes and vertical mixing parameterised?
-What alternatives can be ruled out?
-How sensible are the model results to variation in hydrographic characteristics, atmospheric forcing and nutrient loads?
-How accurate is the model at different locations?
There are two main reasons for applying the 1-D approach of Burchard et al. (2006).First, the model has to remain as simple as possible, both conceptually and computationally for a future incorporation into a 3-D circulation model.Second, the model has to be able to describe main physical and biogeochemical processes within the water column and simulate the dynamics of oxygen on time scales of several years.Before coupling the 1-D ecosystem model with a 3-D model we have to assess the relative importance of different factors controlling the oxygen cycle in the water column of the Baltic Sea.Once the ecosystem model is used in the 3-D modelling framework it will not be possible anymore to clearly separate the contribution of horizontal and vertical transport processes, which will make it very hard to understand the reasons for any disagreements.

Study area
The strong density stratification in the Baltic Sea suppresses vertical mixing of the water and the transport of oxygen from the surface to the bottom.During very exceptional conditions when the inflow lasts long enough (over two weeks) saline water from the North Sea can reach far enough into the Baltic Sea.The saline water is only very slowly mixed with Baltic Sea water and it flows through the Arkona and Bornholm basins in about six months, then further to the central basin of the Baltic Sea, the Gotland Deep.There it is replacing the old Baltic Sea water, which often contains little or no oxygen but some hydrogen sulphide instead (Feistel et al. 2003b).The medium-strength inflows are important as well because they have the potential to renew intermediate layers of the Baltic Proper halocline (Feistel et al., 2006).Since one of our purposes is to explore the influence of the principal hydrographic situation on the oxygen cycle in the water column of the Baltic Sea, the oxygen and nitrogen cycles have been simulated at several stations with very different hydrographic characteristics.For a detailed presentation, we selected four stations with a quite different location in the Baltic Sea, namely: -Gotland (GS, 249 m depth), a very deep central station of the Baltic Proper (20 Each of the first three stations might be considered as a representative station for the corresponding basin (Reissmann, 2006).The regional characteristics of the salinity, potential temperature and oxygen content are represented well by the hydrographic measurements in the corresponding central stations.

Model description and methods
We used the coupled 1-D ecosystem model of Burchard et al. (2006) to study the sensitivity of oxygen dynamics at above selected stations of the Baltic Sea to external forcing.As the physical part of the 1-D ecosystem model GOTM (General Ocean Turbulence Model; www.gotm.net)was applied.Turbulence is modelled with a two-equation turbulence model; one equation for the turbulent kinetic energy and one equation for the dissipation rate of the turbulent kinetic energy.The model includes a simple parameterisation of deepwater mixing.In order to parameterise unresolved turbulence production by internal wave shear, internal wave breaking or Kelvin-Helmholtz instability under stably stratified conditions, a lower limit to the turbulent kinetic energy is set (k min =const).From the large number of well-tested turbulence models implemented in GOTM, the often used κ−ε model is considered to be a very appropriate tool to model the dynamical vertical structure and the actual turbulent diffusive vertical transport in some Baltic Sea stations.
A biogeochemical model of medium complexity (The Baltic Sea Research Institutes Ecosystem Model (ERGOM) with ten state variables) was used in this study (Neumann, 2000;Neumann et al., 2002).This model is of Euleriantype, so all state variables are expressed as concentrations, no matter whether they are dissolved chemicals (e.g.nutrients, oxygen) or particles (e.g.phytoplankton cells).In the model, the oxygen utilisation and production is connected with nitrogen conversion.The oxygen concentration controls processes as denitrification and nitrification.If oxygen is depleted, than nitrate is used to oxidize detritus, and if nitrate vanishes sulphate is reduced to hydrogen sulphide.Hydrogen sulphide is accounted for as negative oxygen concentrations (2H 2 S=−O 2 ).Reduction of nitrate (denitrification) is counted as a loss of nitrogen in the model.In detail, the state variables are: ammonium, nitrate, phosphate, flagellates, diatoms, blue-green algae, detritus, zooplankton, oxygen and sediment detritus.
The model of Neumann et al. (2002) is coupled to the physical model as BIO IOW module of the GOTM package.The GOTM-BIO IOW model was tested by Burchard et al. (2006) for the Gotland station (GS) with water depth of about 250 m.The comparisons between model results and observation data from COMBINE program (under the umbrella of HELCOM) for the period 1983-1991 showed that the hindcasting of interannual variability of nutrients nitrate and phosphate, and phytoplankton was not satisfactory.Burchard et al. (2006) found that the κ−ε model predicts too shallow mixed layers in the Baltic Sea when applied without limitation of turbulent kinetic energy, k min [m 2 /s 2 ].It was illustrated that the parameter k min can act as a tuning parameter of the model (Burchard et al., 1998;Burchard et al., 2006).However, more complete and accurate studies of model sensitivity analysis and/or model skill assessment were not reported.
We run the model for a six year period, from 1 January 1998 to 31 December 2003 whereby the initial values are approximated from available oceanographic measurements.The selected stations (Fig. 1), which might be representative for the corresponding region, are from North to South: Landsort Deep (LS), Gotland Sea (GS), Bornholm Sea (BS), Arkona Sea (AS) and Kattegat (FS).The simulation period includes stagnant (1998)(1999)(2000)(2001)(2002) and fluctuant (2003) periods.The only major inflow to the Baltic Sea during the investigated period was in 2003 (Feistel et al., 2003b).However, several inflows of less strength occurred during the period (Matthäus and Nausch, 2003).
Depth profiles of temperature and salinity along with surface meteorological data and nutrient components were used to force the model.The meteorological forcing data were taken from the ECMWF (European Centre for Medium Range Weather Forecast, www.ecmwf.int)data server (ERA-40 re-analysis data).The frequency of the meteorological data is six hours.Data sets of temperature, salinity, concentrations of oxygen and chlorophyll a were extracted from the Baltic Environmental Database (BED) via the internet based software NEST (http://nest.su.se/bed).The initialization of some initial parameters of the BIO IOW module was done by the use of BED data, as well.Finish Institute of Marine Research (FIMR) Baltic Sea monitoring data (http://www.fimr.fi/en/tietoa/helcomseuranta/en GB/bmp/ data) was also used for model verification.The water transparency of the Baltic Sea, measured as Secchi depth, had been thoroughly estimated in the report of Laamanen et al. (2004) and it was assumed to be 5 m in our calculations.The computed temperature and salinity profiles were www.ocean-sci.net/6/461/2010/Ocean Sci., 6, 461-474, 2010 relaxed towards observational profiles (BED data) or profiles calculated with the GETM model (www.getm.eu;Stips et al., 2005).The optimal relaxation time is about 5 days.The model was run using a two year repeating cycle of forcing data for 1998 as a "spin-up" period in order to achieve a quasi-equilibrium state and obtain reasonable initial conditions.Nutrient fluxes at the air-sea surface were adjusted in order to parameterise lateral nutrient fluxes which are neglected in the 1-D model.Thus, much higher values than the real ones were used in the calculations.In order to highlight the differences between the physical conditions at the studied stations, we fixed the surface fluxes and initial concentrations of ammonium, nitrate, and phosphate for all numerical simulations.The estimation of the nutrient values was done on the base of the sensitivity analysis (Sect.4.5).
Statistics, such as correlation coefficient, R, normalised standard deviation, σ = σ m σ r (σ r and σ m are the standard deviations of the reference and the model field, respectively), and normalised "unbiased" root mean squared difference, Q (normalised by σ r ) are used in the following sections in order to compare the multiple model runs with the reference (observational data).The difference between normalised RMSD and potential bias is denoted with Q.The RMSD is a measure of the average magnitude of the difference, while Q may be conceptualized as an overall measure of the agreement between the amplitude ( σ ) and phase (R) of two temporal patterns.For this reason, R, σ and Q are referred as "pattern statistics".The three pattern statistics are related to one another by (Taylor, 2001) (1) The normalised standard deviation and the correlation coefficient from the model to reference field comparisons may be displayed on a single Taylor diagram (see, for example Fig. 2).On it the distance from the origin is the normalised standard deviation, σ , while the azimuth angle is proportional to arccos (R).Therefore the reference field point has the polar coordinates (1.0, 0).Model to reference comparison points are assessed by how close they fall to the reference point.This distance is equal to Q.The relationship (1) makes the Taylor diagram useful because the individual contribution of misfits of amplitude may be compared to misfits in phase to distinguish how they contribute to the normalised unbiased RMSD.All calculations have been done on the basis of all the available measurements of a selected water column (see, stations in Fig. 1) during the period 1998-2003 and the corresponding model results.It is important to note that the model and reference fields were not log-transformed or averaged in all the presented comparisons.4 Sensitivity analysis

Effect of air-sea exchange
The oxygen exchange with the atmosphere is usually described by where F [mmolO 2 /m 2 s] is the air-sea oxygen flux, V [m/s] is the transfer (piston) velocity, O and O sat [mmolO 2 /m 3 ] are surface and saturation oxygen concentrations, respectively.In the BIO IOW module the piston velocity is assumed as a constant and the saturation oxygen concentration is calculated by where T s is the surface temperature and a 1 , a 2 are constants (Neumann et al., 2002;Burchard et al., 2006).
Another option to estimate the piston velocity is to apply the model of Liss and Merlivat (1986), which includes three regimes (smooth surface, rough surface and breaking waves) depending on the magnitude of wind speed, w: V = 5.9(2.85w− 9.65)/Sc 0.5 at 13 < w[m/s] : V = 5.9(5.9w− 49.3)/ScThe Schmidt number Sc is defined as ratio between the kinematic viscosity and the molecular diffusivity of oxygen.We have applied the following parameterisation for Sc (Stigebrandt, 1991) Equation 5 is valid in the interval 0 < T s < 40 • C and thus, it is applicable in the case of a non-freezing sea surface.Instead of the linear dependence of O sat on temperature like in the BIO IOW module, we have used the correct formula of oxygen solubility of Weiss (1970).
In order to investigate in detail the effect of parameterisation of the air-sea exchange on the surface oxygen dynamics, we considered four cases with different parameterisations of the air-sea exchange and additionally a case without phytoplankton growth and grazing (Table 1).The mean absolute error (mean value of the absolute differences between simulated and measured values) represents the magnitude of the difference between the BED observation data and our numerical simulations, while the linear correlation coefficient, R, measures the strength and the direction of a possible linear relationship between them.A root mean square difference (RMSD [mmol O 2 /m 3 ]) of simulated and measured surface oxygen concentrations is also given in Table 1.It is a measure of discrepancies between simulated and observed surface concentration of oxygen.A complete description of the above mentioned statistics can be found in Taylor (2001).The statistics are calculated on the basis of the available measurements of the surface oxygen content during the studied six year period at GS and the corresponding simulated values.In case I (constant piston velocity and linear oxygen saturation) both mean absolute error and RMSD reach the highest values, while the correlation coefficient has the lowest value (Table 1).We have found the best agreement with the observation data in case IV (new piston velocity and nonlinear oxygen saturation).The comparison of statistics in case II and III with case I shows that the model improvement is caused approximately to the same extend by both variable piston velocity (Eqs.4 and 5) and nonlinear oxygen saturation (Weiss, 1970).
An estimation of the influence of biological activity on the surface oxygen is given as case V in Table 1.The lack of biological activity does not affect considerably the correlation coefficient but it does affect the mean absolute error and RMSD.It is worth noting that even without any primary production (case V) the improved model (case IV) predicts reasonable surface oxygen concentrations.Moreover, statistics in case V are much better than these in case I.The statistics presented in Table 1 clearly indicate that the parameterisation of the air-sea oxygen exchange has a major effect on the surface oxygen dynamics.

Effect of vertical turbulent exchange
The results of 10 separate model runs with different values of k min are shown in Fig. 2. It is a Taylor diagram of the sensitivity of the model to the vertical turbulent exchange.
The diagram shows the model to reference statistics for the oxygen, phosphorus, ammonium, nitrate and chlorophyll a fields during the period 1998-2003 at Bornholm station (BS).The parameter investigated here is the minimum turbulent kinetic energy, k min , which is used in the turbulence model as a parameterisation to account for unresolved mixing processes such as internal waves (Burchard et al., 2006).The colour bar represents 10 different values of k min × 10 7 in the interval [5;30].The selection of this interval is based on the values of k min used in the forth numerical experiment of Burchard et al. (2006).Generally, the best model performance is achieved for oxygen (the highest R values and the smallest Q values).Limiting nutrients have an intermediate goodness of fit (R values ranging from 0.4 to 0.8 and Q values from 0.65 to 1) and chlorophyll a has the highest misfit with the observed values.The spread of comparison points in Fig. 2 demonstrates that k min is an important parameter for predicting all the presented state variables.We found the best fits for nutrients in the interval 6 × 10 −7 ≤ k min ≤ 8 × 10 −7 [m 2 /s 2 ] (blue colour in Fig. 2).It appears, that chlorophyll a is not strongly influenced by the variation of k min .Since our interest is mainly related to the oxygen dynamics, we will discuss in detail the sensitivity of oxygen to changes in the vertical turbulent mixing.[m 2 /s 2 ]) and underestimates it at high k min (> 15.10 −7  [m 2 /s 2 ]).The value of σ changes rapidly with increasing k min , while the value of R does not.In other words, the vertical turbulent mixing has a higher influence on the amplitude rather than on the phase of the simulated oxygen field.Both minimum of the total RMSD (indicated by "♦") and minimum of the unbiased RMSD (dotted lines on the Taylor diagram are the isolines of Q) are found for k min = 1.10 −6 [m 2 /s 2 ].Thus, the bias between modelled and reference fields has also a minimum at this point.
The best fit between the model and reference oxygen fields was found for k min ×10 7 [m 2 /s 2 ]: 8 at GS; 25 at AS; 80 at FS; 5 at LS.These particular values of k min are used in all simulations presented below.It is obvious that k min is an important model parameter and one must decide carefully how to parameterise it when one couples the GOTM-BIO IOW model with a 3-D circulation model of the Baltic Sea.
There is a decreasing trend of the optimal k min for oxygen (80; 25; 10; 8; 5)×10 −7 [m 2 /s 2 ] with distance from the entrance of the Baltic Sea, which might reflect the decrease in the effective vertical exchange in the Baltic.The strength of the density stratification expressed as the observed mean vertical density difference (bottom to surface), ρ t [kg/m 3 ] for the period 1998-2003, shows a similar spatial pattern: 11.56 at FS; 8.17 at AS; 8.63 at BS; 6.4 at GS and 6.41 at LS.

Effect of relaxation to temperature and salinity profiles
As it has been mentioned in Sect.3, the model is relaxed to prescribed depth profiles of temperature (T ) and salinity (S).The relaxation to T /S profiles is necessary for 1-D simulations in an environment where lateral processes cannot be neglected (Reissmann et al., 2009).The model performance depends on the salinity relaxation time scale rather than on that of temperature.
We have performed two separate runs with two different types of T /S relaxation profiles: (i) experimental profiles from BED; (ii) calculated profiles from a 3-D simulation.A Taylor diagram is drawn in Fig. 3 which is comparing T , S and O (oxygen) calculated by the use of BED data for T /S relaxation (comparison points are denoted with small letters) and 3-D model values (denoted with capital letters).The pattern statistics of T , S and O are normalised by the standard deviation of the corresponding observation field.The different colour of the letters in the diagram refers to the different stations.The overall impression from Fig. 3 is that T /S fields in the run (ii) partially disagree with the experimental ones (reference).This is especially valid for salinity with σ ≈ 0.45 and Q ≈ 0.57 at BS and GS, however the model salinity field is well phased for all stations (R ≥ 0.88).It is not a surprise that all comparison points of T /S in the run (i) are very close to the reference.
It can be seen in Fig. 3 that the forcing with BED data gives slightly better results for oxygen.The close coincidence of the oxygen comparison symbols for FS and AS (yellow and green letters) points to the low sensitivity of the oxygen dynamics at these stations to the prescribed salinity field.The influence of the T /S forcing data is more pronounced for the other two stations and in particular for BS where σ = 0.8 and R = 0.95 in the run (ii).However the agreement of the simulated oxygen with the observation data is better in the run (i).Despite the underestimation of salinity, the good results for oxygen demonstrate that it is possible to utilise 3-D model data for T /S relaxation in all cases when observation data is scarce or absent.

Effect of atmospheric forcing
In order to investigate the model sensitivity to variations in atmospheric forcing, we present results from five different cases and compare them with the observational data.The normalised pattern statistics of oxygen have been calculated for the period 1998-2003 after varying the wind speed values in the ERA-40 re-analysis data.Namely, the wind speed has been rescaled by a factor of 0.5, 0.8, 1.0, 1.2, and 1.5 (plotted with different colours in Fig. 4).The close grouping of the comparison points for GS (circles) indicates that the oxygen dynamics at this deep station is not sensitive to a possible uncertainty in the forcing data.We get significant changes in the modelled oxygen for all other stations.Particularly, when the wind speed is scaled down the comparison points are farther away from the reference than when it is scaled up.In summary, one can conclude that an increase of wind speed by a factor of 1.2 has led to a general improvement in the model performance.For the scaling factor of 1.5 the correlation is slightly improved for FS and AS, even though the results for σ and Q are worse for BS.Another possible inference drawn from Fig. 4 could be that the wind speed magnitude of the ERA-40 reanalysis could be underestimated.

Effect of limiting nutrients
In the model, the nutrient load is taken into account via initial concentrations and surface fluxes of nitrate, phosphate and ammonium.For the 1-D model considered here, the nutrient fluxes at the air-sea interface have to be adjusted in order to parameterise lateral nutrient fluxes.A Taylor diagram is drawn in Fig. 5 for testing the model sensitivity to limiting nutrients.It shows the model to reference statistics for oxygen (red) and chlorophyll a (green) at BS.The results of 150 separate model runs are shown on the diagram and the corresponding intervals from which the initial concentrations and the surface fluxes of nutrients are randomly chosen are given in Table 2.The surface fluxes of nutrients are assumed as constants during a single model run.The average values (for the upper 20 m) of chlorophyll a are used for comparison.It appears that both oxygen and chlorophyll a are weakly sensitive to the variation in the concentrations of nutrients.Moreover, only the amplitude of the model oxygen field is sensitive, while the phase remains approximately unchanged (R ∼ = 0.95).The low sensitivity of the oxygen and chlorophyll a fields to a relatively big variation in the values of the nutrient surface fluxes could be explained with the simple flux parameterisation used (as a constant).Typically, the surface water concentrations of nutrients in the Baltic Sea are very low in summer and high in winter.It worth noting, that at these points the unbiased RMSDs have also a minimum.
The initial concentrations and surface fluxes of nutrients for which we have found the best fits for oxygen and chlorophyll a are given in Table 2. 5 Model results and validation

Water column structure
The annual temperature variation in the surface water of the Baltic Sea is great, with differences of up to 20 • C. For illustration, in Fig. 6a is shown the surface temperature at station BS.The surface temperature at GS (also at FS and AS) behaves in the same way like that at GS, while the bottom one is approximately constant (7 • C) at both stations (it decreased to 3 • C only after the cold inflow of 2003).At BS the surface salinity is about 7.5 (7 at GS) and the bottom salinity varies slightly between 15 and 17.5 (12 and 13 at GS) and reaches a peak of 19.2 after the inflow in 2003.A halocline separates the surface waters with lower salinity (6-9) from the deep waters with higher salinity (15-20), (for all stations except for FS, where the surface salinity varies between 16 and 30 and the bottom one between 33 and 35) and excludes the deep water from vertical mixing.The halocline begins at a depth of about 10-20 m in the Fladen station, 30-40 m in the Arkona basin, 35-50 m in the Bornholm basin, and 60-70 m in the Gotland basin (IOW, 2003;Wasmund et al., 1998).
For illustrating the seasonal variability of the density stratification, in Fig. 6b is shown the comparison between the simulated and observed density difference, ] (where ρ b and ρ s are the bottom and surface density, respectively) at BS.It particularly indicates the less stratified winter period and the presence of more stable conditions in summer (Lass and Mohrholz, 2003;Mohrholz et al., 2006;Sellschopp et al., 2006).The variability of ρ t is simulated quite well, because of the applied salinity relaxation.
In summer, additionally a thermocline forms at about 15-20 m depth and the temperature of the intermediate water between thermocline and halocline usually remains the same as during winter (4-10 • C).The thermocline exists until October, then in autumn the surface water starts cooling and sinking until it reaches the temperature of maximum density.The thermocline and the related density differences in the upper layer disappear and finally wave and wind actions mix the whole layer above the halocline.
The vertical oxygen distribution at BS is shown in Fig. 7 for selected representative days during the year 2001.It is nearly constant in the layer above the halocline except for the summer months.Moreover, the concentrations of oxygen are higher in the layer below the thermocline (cold intermediate layer) than in the other water layers.In the halocline oxygen decreases rapidly, so the halocline acts as a barrier for oxygen transport into the deeper waters.Thus, one can distinguish three main layers of the water column at BS (as well as at the other three stations): -surface (mixed) layer, where the temperature, the salinity and the oxygen concentrations are more or less vertically constant; -intermediate layer (the depths below thermocline till the end of halocline), where the temperature, the salinity and the oxygen concentrations change significantly; -bottom layer, where the temperature, the salinity and the oxygen concentrations become approximately constant.
In the surface layer, the calculated oxygen concentrations are in perfect agreement with the measurements.Then, in the intermediate layer the model well predicts the trends in the vertical distribution of oxygen.In the bottom layer, the calculated concentrations of oxygen are reasonable however not match well the observations.Negative oxygen is the amount of the oxygen equivalent to the amount of hydrogen sulphide produced through reduction of sulphate.Generally, the vertical structure of oxygen is highly correlated with the measurements in each period of the year.Correlation coefficient, R, normalised standard deviation, σ and RMSD of oxygen concentrations are given in Table 3.
The statistics are calculated on the basis of the available measurements for the full water column during the year 1998 at five stations and the corresponding simulated values.In addition to the statistics for the four studied stations, the statistics for the Landsort station (LS), 440 m depth (see Fig. 1), is also presented to support the model validation.The measured oxygen concentrations of each observation have been interpolated on the computational grid of the water column and then R, σ , and RMSD are calculated (the same procedure has been done for the statistics presented in Table 4).
It should be noted that the number of observations at each station is about 15 per year and the number of observation points in the water column related to the station depth is also similar for all stations.Therefore, we can consider the statistics of these stations as equally reliable.The model-data agreement is perfect for BS, GS and LS and nearly perfect for the other two stations.The relatively low values of the RMSD in comparison to the variability of the data indicate a close match between predicted and observed concentrations.In summary, the statistical quantities support our conclusions that the model successfully reproduces the vertical water column variability of the oxygen.

Surface and intermediate layer
The model results are analysed at the previously identified three main water column layers for the period 1998-2003.Figure 8 shows the modelled time series of surface oxygen for all stations compared with the BED and FIMR data.The time interval between two subsequent major ticks in all time series plots is 2 months.At the surface, the modelled oxygen is in a near-perfect agreement with the observations.The observed increase of surface oxygen concentrations with decreasing surface temperature is well captured by the model.Some peaks of the oxygen concentrations are underestimated for FS (Fig. 8d) for the years 1999, 2002 and 2003.This is partially caused by the inflow of oxygenated surface water from the North Sea at this station.FS is placed in the Kat- tegat, close to the North Sea, where the surface water salinity is affected by the irregular inflows and outflows of salty or brackish water, respectively.Despite this little misfit, the adequate accordance between simulations and data indicates that the parameterisation (2), ( 4) and ( 5) used to compute the surface oxygen flux is appropriate for the Baltic Sea and the time evolution of surface oxygen is mostly determined by the gas exchange at the surface.The oxygen time series in the intermediate layer are shown in Fig. 9.In the intermediate layer the model matches very well the data, see the stations GS (50 m depth in Fig. 9a), BS (40 m depth in Fig. 9b) and AS (20 m depth in Fig. 9c).The discrepancy between calculated and observed concentrations of oxygen is higher at station FS (Fig. 9d).This is expected because the influence of horizontal advection is more pronounced at FS than at the other selected stations.

Bottom layer and inflow dynamics
The model performance in deep water layers, at the bottom, is not really satisfactory (Fig. 9).The sediment oxygen demand is only partially taken into account in the model and therefore the simulated bottom oxygen is approximately constant in time at the deep stations GS and BS.A modified model of Neumann et al. (2002) including a non-Redfield stoichiometry has not led to a significant improvement of the simulated near bottom oxygen at a central Gotland Sea station (Kuznetsov et al., 2008).The introduction of a prognostic sediment layer is still an ongoing development for this model.Contrary to the surface layer, the horizontal advection of oxygenated water is a very important component of the oxygen dynamics in the bottom layer.This can be clearly seen by the sudden increases in bottom oxygen in Figs.9a,b, which are linked to the major inflow in January 2003.Notwithstanding that the main hydrographic conditions of the Baltic Sea are characterised by permanent salinity stratification, these conditions are not the same for the different regions of the Baltic Sea.The bottom temperature is about constant except during inflow events for the stations of the Baltic Proper and varies seasonally at AS and FS.The variation of bottom salinity is shown in Fig. 10, exhibiting very distinct features in the different regions.In the Kattegat area (Fig. 10a) frequent small variations occur (∼1), in the Arkona basin (Fig. 10b) we find large quasi seasonal fluctuations (∼10) and in the Bornholm Sea (Fig. 10c) there is a steady decrease interrupted by irregular inflow events.
When one compares the evolution of bottom temperature and salinity with bottom oxygen at FS, AS and BS, several features of these time series can be underlined.First, at the Ocean Sci., 6, 461-474, 2010 entrance area (small bottom salinity fluctuations) where station FS is placed, we find a strong negative correlation between observed temperature and oxygen fields of R(T ,O 2 ) = −0.82.This seems to indicate that the bottom oxygen concentration at FS is mainly determined by the inverse temperature dependent oxygen saturation and less by the sediment oxygen demand.
Second, the picture at AS is more complicated, as both temperature and salinity variations are influencing the bottom oxygen dynamics.For example, at the end of August 1999 there is a summer inflow of warm and salty water from the Belt Sea, at the same time oxygen is depleted and therefore it reaches a minimum.A similar event happens in July/August 2002.Contrary to these low oxygen summer inflows, the normal winter inflow of high salinity and low temperature is oxygen saturated (see the inflow in January 2003).As a consequence of these different inflow types, the correlation between observed temperature and oxygen fields in the bottom layer at AS is R(T ,O 2 ) = −0.57.
Third, in the Baltic Proper, at BS the increase of oxygen usually corresponds to a sudden increase in salinity (Fig. 10c).The correlation between salinity and oxygen in the bottom layer seems to indicate that the increase of the near-bottom oxygen is due to infrequent pulses of North Sea inflow.The warmer inflows in October 1999 and December 2001 (Feistel et al., 2003a) bring less oxygenated water, however they ventilate to some extend the bottom water at Bornholm station.The decrease of oxygen during stagnation periods which is caused by the sediment oxygen demand is practically not captured by the model.It appears that at BS and GS, the bottom oxygen flux is of higher importance for the oxygen concentration.Therefore, a better parameterisation of the oxygen consumption in the bottom layer is essential.
In summary the seasonal variability of bottom oxygen at stations FS and AS is only partially matched by the model; it is only capable to capture the variation to some extent, with a reduced range of amplitudes and with a phase shift of 1-2 months.The oxygen dynamics in the more central Baltic stations (BS and GS) cannot be reproduced by a 1-D model.However, the discrepancy between model and observational data is not only due to the omitted horizontal advection as the 3-D circulation model used by Neumann et al. (2002) predicts also too high values of the near bottom oxygen at BS and GS (stations 213 and 271 in Fig. 13 of Neumann et al., 2002) for the period from 1983 until 1990.Unfortunately, it seems that the simulation of the horizontal transport of their 3-D model is too diffusive (see the near bottom salinity in Fig. 5 of Neumann et al., 2002), so that likely in the simulations no inflowing oxygen rich water has arrived at the bottom of the Gotland Sea.In Neumann and Schernewski (2005) the vertical resolution has been increased, which led to some improvement of the near bottom oxygen concentrations.The calculated time series of near bottom oxygen intersects with the observational data without show- ing any inflow dynamics.Evidently, the near bottom oxygen dynamics and near bed consumption are not well considered and further adjustments to the model are necessary.However, even the correct accounting of the sediment oxygen demand will not lead to improved 1-D simulations, as we have to consider advection by applying a 3-D model or at least parameterise the effect of the inflow events on the oxygen concentrations for the 1-D runs.

Summary statistics for a 6 year period
Summary statistics of the interannual model performance (Table 4) shows a high correlation between the observed and modelled values; R and σ are close to one, the RMSD are relatively small, although they are higher than those for the year 1998 (Table 3).Especially, the RMSD at GS is approximately doubled due to the existence of permanent mismatch between model results and observational data in deep water layers (Fig. 9a).The summary statistics are generally less favourable for FS and AS than for BS and GS.The modelled values of oxygen underestimate the measured ones ( σ = 0.71at AS and σ = 0.65 at FS).The low values of the correlation coefficient at FS and AS are expected because of the phase shift in the bottom oxygen time series (Fig. 9c, d).This discrepancy is mainly due to the effect of inflow dynamics and oxygen sediment demand which are fairly considered in the model.Unfortunately, there is not enough observation data to check this assumption.
Thus, the statistics presented in Table 4 confirms the information shown in the time series plots (Figs. 8 and 9).The agreement between modelled and observed oxygen concentrations will be a little better if we would exclude the year 2003 from the comparisons.This exclusion could be justified for our 1-D simulations because of the occurrence of the major inflow event in January 2003 (Feistel et al., 2003b), which would require the consideration of horizontal oxygen transport.

Chlorophyll a simulation
Biological activity is another major factor controlling oxygen concentrations.The interannual variability of simulated and observed average phytoplankton concentrations, shown as average chlorophyll a (Chl a) is given in Fig. 11.The time series of calculated Chl a concentrations and in-situ data of BED and FIMR correspond to the water column average values (from the surface to 20 m depth).Also presented in the figure are the monthly mean values taken from satellite images (Environmental Marine Information System (EMIS) database, http://emis.jrc.ec.europa.eu/).The model predicts a spring bloom mostly composed of diatoms and flagellates in the beginning of March for BS (Fig. 11a) and in the beginning of April for GS (Fig. 11b).To some extent this result coincides with the HELCOM (1996) report stating that the spring bloom of phytoplankton develops earlier at the western part of the Baltic Sea than in its eastern and northern parts.In the latter areas, a strong spring bloom develops in April/May, followed by a small summer bloom in July/August and an autumn bloom in October/November.After mild winters, the spring bloom could appear earlier.Also, the regional differences in timing and strength of the spring bloom are related to the mixing depth (Wasmund et al., 1998) and the strength of the winter deep mixing (Janssen et al., 2004).There is weak evidence of a summer bloom in the model results at BS (Fig. 11a), however, it is not simulated for GS (Fig. 11b) by the model.Typically, the autumn bloom is predicted to develop in September/October.The autumn peak is well phased and the timing corresponds to all presented observation data.There is a reasonable agreement between the modelled and observed average Chl a in 2003 at BS, however, in all other years the model predicts lower bloom peaks than the observed ones at both stations BS and GS.Still one has to keep in mind that comparing insitu and model data involves not only model but also data uncertainties.It should be obvious that the typical random pull of a bucket of water out of a patchy plankton bloom might lead to a drastic over-or underestimation of the real mean Chl a concentrations in the measurement area.This could be overcome only by rather expensive measurement methods as for example taking about 100 random samples within the comparison region in order to establish statistical means and confidence intervals for the measurements.Additionally, as depicted in Fig. 11, there is not a good agreement between both measured data types (in-situ and satellite data).The satellite data are often missing the spring bloom peak, which might be related to cloud cover during that time.An interesting finding is that the model shows better succession in the phytoplankton content for the years when in-situ and satellite data match better.This might be an indication that in such a case of agreement the observed data are more representative of the real situation in the field.Despite the above mentioned limitations of the model, we can conclude that under the influence of atmospheric forcing and at different hydrographic characteristics the model reproduces the annual and interannual variability of oxygen typical for the Baltic Sea.

Summary and Conclusions
In the present work we have examined the influence of some important physical and geochemical factors on the oxygen concentrations at several regions of the Baltic Sea.For this purpose we have used the GOTM-BIO IOW model.The model has been forced with meteorological data for a six year period.Modifications in the parameterisation of the air-sea oxygen fluxes have led to a significant improvement of the model results in the surface and intermediate water levels.Sensitivity analysis has been performed in order to examine the influence of turbulent mixing, hydrographic forcing (salinity and temperature profiles used for relaxation), atmospheric forcing (wind speed), and nutrient loads.The normalised standard deviation, the correlation coefficient and the normalised unbiased RMSD from each model to reference field comparison are displayed as Taylor diagrams.
There can be no doubt that the 1-D model used here allows to distinguish the most important driving mechanisms of the oxygen dynamics in the Baltic Sea.An additional strong argument supporting 1-D modelling is the fact that the used ecosystem model ERGOM has been applied in many 3-D studies of the Baltic Sea (Neumann, 2000;Neumann et al., 2002;Neumann andSchernewski, 2005, 2008), but the apparent poor parameterisation and simulation of the surface oxygen dynamics in the model passed unnoticed for more than 10 years.We would therefore argue that before any ecosystem model is used in 3-D simulations, it should be vigorously tested in a 1-D framework.We could confirm that the major factors controlling the oxygen dynamics in the Baltic Sea are natural physical factors, like the magnitude of the vertical turbulent mixing, wind speed and the variation in temperature and salinity.The influence of limiting nutrients is less pronounced, at least under the nutrient flux parameterisation assumed in the model.
The interesting fact that the minimum kinetic energy used in the turbulence model giving the best fit of simulations to observations is decreasing with the distance from the entrance of the Baltic Sea, namely, k min = (80; 25; 10; 8; 5)×10 −7 [m 2 /s 2 ], could be a hint to unresolved mixing due to e.g.breaking internal waves as the strength of the density stratification is decreasing in a similar way.Further this clearly underlines the fact that the use of a spatial and temporal constant k min in 3-D applications is inappropriate, an improved parameterisation is urgently needed.
A model validation has been done by evaluating the agreement between predicted values of oxygen and observation data from the BED and FIMR data bases.The correlation with observation data is good and consistent for all stations and with low values of the RMSD (Tables 3 and 4).Specifically, the oxygen dynamics of the surface mixed layer is simulated in close agreement with the observations.The fact that the oxygen dynamics at the surface can be accurately simulated by a 1-D model has been already shown by Vichy et al. (2004) for the BS during the stagnation period 1979-1990and by Kuznetsov et al. (2008) ) for the Central Baltic Deep during 1978-1993.However, it comes certainly as a surprise that even the very dynamic transitional stations FS and AS are very well simulated by the 1-D model, which is ignoring completely the advection of oxygen.And this remains true even in the case when a major inflow event appears like this in 2003.Therefore, it can be concluded that in the surface layer the dynamics of the mixed layer and the oxygen exchange with the atmosphere are the controlling factors of near surface oxygen development.As it has been shown, these physical factors also clearly dominate the biological production and respiration of oxygen in the surface layer.
The largest mismatch with observations occurs when simulating the bottom water oxygen dynamics.This is of course expected, as the bottom oxygen concentrations in the Baltic Sea are not only determined by the local sediment oxygen demand, but largely influenced by inflowing oxygenated water from the North Sea.As we have not taken into account the horizontal advection of oxygen in the 1-D model, we could not simulate the increase of bottom oxygen during inflow events.Nevertheless, it is obvious that the oxygen consumption at the sediment interface demands for an improved parameterisation.However, one has to keep in mind that when incorporating a better sediment oxygen demand parameterisation in a 1-D model, the results of the simulation could even become worse because of the fact that an eventual higher consumption will not be counterbalanced by oxygen transport.The statistical properties of the modelled nutrient and phytoplankton concentrations are also reasonable.This demonstrates the good capability of the model to predict the oxygen dynamics at all selected stations.

Fig. 2 .
Fig. 2. Taylor diagram of the model sensitivity to the vertical turbulent exchange parameterisation (different values of k min are used).Model to reference statistics are presented for oxygen (denoted with circles "•"), phosphorus (denoted with triangles " "), ammonium (asterisks " * "), nitrate (diamonds " ") and chlorophyll a (squares " ") fields in the period 1998-2003 at BS.The colour bar represents 10 different values of k min × 10 7 in the interval [5;30].The minimum value of the RMSD for oxygen is indicated by black diamond ("♦").The black circle at a correlation of 1.0 and normalised standard deviation of 1.0 represents the reference.
Figure 2 clearly indicates that the model overestimates the seasonal variability at low k min (< 7.10 −7

Fig. 3 .
Fig. 3. Taylor diagram comparing T , S and O calculated by the use of BED for T /S relaxation (comparison points are denoted with small letters) and 3-D model values for T /S relaxation (denoted with capital letters).The pattern statistics of T , S and O are normalised by the standard deviation of the corresponding observation field.Different letter colour in the diagram refers to a particular station: BS -blue; GS -red, AS -green; FS -yellow.

Fig. 4 .
Fig. 4. Taylor diagram comparing oxygen fields calculated using different meteorological forcing at the principal stations for the period 1998-2003.Symbols representing the different stations are labelled in the figure.Different symbol colours categorise the variation of wind speed scaling: 0.5 -blue; 0.8 -red; 1.0 -green; 1.2 -yellow; 1.5 -magenta.

Fig. 5 .Fig. 6 .
Fig. 5. Taylor diagram for the model sensitivity to limiting nutrients showing model to reference statistics for the oxygen (red) and the chlorophyll a (green) field for the period 1998-2003 at BS.The comparison points with the minimum RMSD values are indicated by a black diamond ("♦").The ranges of the intervals in which vary the initial concentrations and surface fluxes of nutrients are given inTable 2.

Fig. 7 .
Fig. 7. Vertical oxygen profiles at BS in some selected days of 2001.Calculated results are presented with a solid line, while circles connected with a dashed line show the observation data of BED.

Fig. 8 .
Fig. 8. Simulated six year time series of surface oxygen.Calculated results are presented with a thick solid line, FIMR data with asterisks, and BED data with diamonds: (a) -at GS; (b) -at BS; (c) -at AS; (d) -at FS.

Fig. 9 .
Fig. 9. Oxygen time series in the intermediate and bottom layer for a six year period.Calculated results are presented with a thick solid line: magenta colour (intermediate); red colour (bottom).BED data is denoted by black squares (intermediate) connected with a thin line and blue circles (bottom) connected with a thin line.FIMR data is presented with green asterisks.(a) -at GS; (b) -at BS; (c)at AS; (d) -at FS.

Fig. 10 .
Fig. 10.Time series of bottom salinity for a six year period.Calculated results are presented with a thick solid blue line, observation data of BED with magenta squares and FIMR data with green asterisks.(a) -at FS; (b) -at AS; (c) -at BS.

Fig. 11 .
Fig. 11.Modelled (thick solid line) and in-situ data (denoted with blank magenta diamonds and green asterisks) of average Chl a [mg/m 3 ] at: (a) -BS, (b) -GS.Data from satellite images (EMIS database) is presented with filled black circles connected with a dash line.

Table 1 .
The impact of the air-sea oxygen exchange parameterisation and primary production on surface oxygen concentrations.

Table 2 .
Ranges of initial concentrations and surface fluxes of limiting nutrients used in the sensitive analysis and the corresponding values, for which the minimum of the RMSD has been found.

Table 3 .
Correlation coefficient, R, normalised standard deviation, σ , and root mean square difference, RMSD [mmolO 2 /m 3], of the simulated and measured oxygen concentrations in the full water column for n days in the year 1998.

Table 4 .
The same statistics like in Table3but for a six year period.