The relative importance of selected factors controlling the oxygen dynamics in the water column of the Baltic Sea

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 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 Turbulent 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 longterm 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.


Introduction
The Baltic Sea is a semi-enclosed and brackish sea, which together with other 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. Pulses of water streaming in at the bottom through the Danish straits transport salty and oxygen rich water from the North Sea into the Baltic Sea (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. Only in 1993 and 2003 such major inflows took place (Jakobsen, 1995;. Inflows from the North Sea are currently the principle 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 2 ml/l at around a depth of 100 m, or even more shallow than that (HELCOM, 2003). At the same time, the area covered by hydrogen sulphide extends from the main eastern Basin of the Gotland 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 nearbottom water, although the level is lower compared to the long-term measurements. The oxygen conditions in the bottom waters of the Baltic Proper continues to be bad during -2006as well (HELCOM, 2007. The zones on the seabed with anoxic areas where hydrogen sulphide forms increase both in size and volume. More phosphorous consequently diffuses out of the sediments and into the deep waters of the Baltic. Additional to the above mentioned horizontal advection of oxygen the principal 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 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.
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 more susceptible to oxygen deficiency. These models should 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 present study aims to assess the relative importance of different factors controlling the oxygen cycle in the water column of the Baltic Sea by the use of a 1-D water column model. Thus, the relative importance of following factors is investigated in detail: -the significance of the principal hydrographic situation is studied by comparing several stations with very different hydrographic characteristics; -the importance of the accuracy of hydrographic characteristics (temperature/salinity structure) -by comparing simulations relaxed with measured profiles and 3-D model results; -the effect of the vertical turbulent exchange -by varying the parameters of the turbulence model; -the influence of the atmospheric forcing -by multiplying the wind speed by a factor from the interval [0.5;1.5]; -the importance of the parameterisation of the air-sea oxygen exchange -by analysing the impact of different available parameterisations; -the relative importance of limiting nutrients.
The study is organized as follows. In Sect. 2 we describe briefly the 1-D model and characterise the type of the method used to model the system, while in Sect. 3 we provide the model setup and forcing. Section 4 shows the effect of the air-sea oxygen parameterisations on the surface oxygen dynamics. In Sect. 5 are presented model results at different stations 4 and comparisons between observations and model results. The model sensitivity analysis is presented in Sect. 6. The last section includes a discussion and some conclusions.

Model description
We use the coupled 1-D ecosystem model of Burchard et al. (2006) to simulate the oxygen and nitrogen cycles in some selected stations of the Baltic Sea. As a physical part of the 1-D ecosystem model the GOTM (General Ocean Turbulence Model, www.gotm.net) is applied.
The 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 ( const k = min ). We have found that from the large number of well-tested turbulence models implemented in GOTM, the -ε κ model is 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 (ten state variables) is used in this study (Neumann, 2000;Neumann et al., 2002). This model is of Eulerian-type, 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). For example, the ERSEM (European Regional Seas Ecosystem Model, Baretta et al., 1995) is an Eulerian-type model of higher complexity. In the model, the oxygen utilisation and production is connected with nitrogen conversation. The oxygen concentration controls processes as denitrification and nitrification. . It is illustrated that the parameter min k 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 have not been reported.
The validity of a 1-D approximation in the Baltic Proper is confirmed also by some other model results (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, 2001), it is a good tool for basic studies, improving the model parameterisation and investigation of some system properties.

Model forcing and setup
The model is run for a six year period, from January 1, 1998 to December 31, 2003 and the initial profiles are approximated from available oceanographic measurements. 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 . However, several inflows of less strength occurred during the period  has been thoroughly estimated in the report of Laamanen et al. (2004) and it is assumed to be 5 m in our calculations.
Nutrient fluxes at the air-sea surface have been adjusted in order to parameterise lateral nutrient fluxes which are neglected in the 1-D model. Thus, much higher values than the real ones are used in calculations. In order to highlight the differences between the physical conditions at the studied stations, we fix the surface fluxes and initial concentrations of ammonium, nitrate, and phosphate for all numerical simulations. The estimation of the nutrient values is done on the base of sensitivity analysis. Statistical and graphical techniques are applied to compare quantitatively the multiple executions of the model (Sect. 6.4).
The computed temperature and salinity profiles have been relaxed towards observed profiles (BED data) or profiles calculated with GETM model (www.getm.eu, Stips et al. 2005). The optimal relaxation time is about 5 days. The model is 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.

Improvement of the model
In this section we discuss the effect of parameterisation of the air-sea exchange on oxygen dynamics. The oxygen exchange with the atmosphere is usually described by where s T is the surface temperature and 1 a , 2 a are constants (Neumann et al., 2002;Burchard et al., 2006). First, we have implemented the model for the station BY15 in the central Gotland Sea. In Fig. 2 a, b are shown the surface temperature and oxygen time series, respectively. The model is in a good accordance with the data over the full six year period, especially in describing the seasonal variability. However, it captures the variation only with lower amplitudes of surface oxygen concentrations during summer (Fig. 2b). The difference between simulated and observed surface concentrations of oxygen is more pronounced during summer of 1999, 2001-2003 when the surface temperature reaches about 23 C˚ (Fig. 2a). This discrepancy is due to an overly simplified computation of the oxygen surface fluxes. So, we have modified the parameterisation of the surface oxygen flux in the BIO_IOW module.
In this paper, the piston velocity is calculated by 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 : The Schmidt number Sc is defined as ratio between the kinematic viscosity and the molecular diffusivity of oxygen. We have applied the following expression for Sc (Stigebrandt, 1991 ) of simulated and measured surface oxygen concentrations is also given in Table 1. It measures the size of discrepancies between simulated and observed values. A complete description of the above mentioned statistics can be found in Taylor (2001) among other sources. The statistics are calculated on the basis of the available measurements of the 8 surface oxygen content during the studied six year period at the station BY15 and the corresponding simulated values. As shown in Table 1, 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. We have found the best agreement with the observation data in case IV (new piston velocity and nonlinear oxygen saturation). The improvement is caused approximately to the same amount by both new piston velocity and new nonlinear oxygen saturation, as can be seen from Table 1 (please compare case II (only nonlinear oxygen saturation) and case III (only new piston velocity)). It is worth to note, that even without any primary production (case V) the improved model predicts reasonable surface oxygen concentrations that are better than in the case with a linear dependence of sat O on temperature and constant piston velocity (case I). This suggests that the parameterisation of the air-sea oxygen exchange has a major effect on the surface oxygen dynamics.

Model results and validation
The 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.

Water column structure
The annual temperature variation in the surface water of the Baltic Sea is great, having differences of up to 20˚C. For example, in Fig (Lass et al., 2003;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 the 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 BY5 is shown in Fig. 4 Table 2. The statistics are calculated on the basis of the available measurements of the full water column during the year 11 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, BY31, 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 3). It should be noted that the number of observations at each principle 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 BY5, BY15 and BY31 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, this information supports 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 5 shows the modelled time series of surface oxygen for stations BY0, BY1 and BY5 compared with the BED and FIMR data (see Fig. 2c for station BY15).
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  Fig. 7a) and BY1 (20 m depth in Fig. 8a). The discrepancy between calculated and observed concentrations of oxygen is higher at station BY0 (Fig. 9a).
This is expected because the influence of horizontal advection is more pronounced at BY0 than at the other selected stations.

Bottom layer and inflow dynamics
The . This seems to indicate that the bottom oxygen concentration at BY0 is 13 mainly determined by the inverse temperature dependent oxygen saturation and less by the sediment oxygen demand.
Second, the picture at BY1 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 (Fig. 8). A similar episode 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 BY1 is 57 Third, in the Baltic Proper, at BY5 the increase of oxygen usually corresponds to a sudden increase in salinity (Fig. 7). This relation is true also for bottom temperature and oxygen except for the winter inflow in 2003. 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  bring less oxygenated water, however they ventilate to some extend the bottom water at Bornholm station (Fig. 7). 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 BY5 and BY15, the bottom oxygen flux is of higher importance for the oxygen concentration. Therefore, a better parameterisation of the oxygen depletion in the bottom layer might be essential.
In summary

Summary statistics for full water column
Summary statistics of the interannual model performance (  (Figs. 8 a and 9 a). 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 2 confirms the information obtained by the time-series plots (Figs. 6-9). It should be noted here that 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 , which would require the consideration of horizontal oxygen transport.

Chlorophyll a simulation
Biological activity is another factor controlling oxygen concentrations. The interannual variability of simulated and observed average phytoplankton concentrations, shown as average chlorophyll a (Chla) is given in Fig. 10 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 BY5 (Fig. 10 a) and in the beginning of April for BY15 (Fig. 10 b). To some extent this result coincides with HELCOM (1996) report stating that the spring bloom of phytoplankton develops earlier at the western part of the Baltic Sea then in its eastern and northern parts. In these 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 a weak evidence of a summer bloom in the model results at BY5 (Fig. 10 a), however, it is not simulated for BY15 (Fig. 10 b) by the model. Typically, the autumn bloom is predicted to develop in September/October. The autumn peak is well phased and corresponds to all presented observation data. There is a reasonable agreement between the modelled and observed average Chla in 2003 at BY5, however, in all other years the model predicts lower bloom peaks than the observed ones at both stations BY5 and BY15. A part of the discrepancy between calculated concentrations of chlorophyll a and observed values could be explained by the simplified parameterisation used for chlorophyll in the model, which is a simple linear function of the N-content (Janssen et al., 2004). Still one has to keep in mind that comparing in-situ and model data involves many uncertainties, as 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 Chla 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. 10, 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 cycles of oxygen typical for the Baltic Sea.

Sensitivity analysis
Statistics, such as correlation coefficient, R , normalised standard deviation, σ , and the normalised "unbiased" root mean squared difference, S (normalised by r σ ) are used to compare the multiple model runs with the reference (observation) data. The difference between normalised RMSD and potential bias is denoted with S . The RMSD is a measure of the average magnitude of the difference, while S 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 S are referred as "pattern statistics". The three pattern statistics are related to one another by (Taylor, 2001) R S σ σ2 The normalised standard deviation and the correlation coefficient from the model to reference field comparisons may be displayed on a single Taylor diagram (for example, see Fig. 11).
The Taylor diagram is a polar coordinate diagram with polar angle proportional to ) ( arccos R and radial distance from the origin proportional to σ . Therefore the reference field point has the polar coordinates (1.0, 0). The model to reference comparison points are then assessed by how close they fall to the reference point. This distance is equal to S . The relationship (5) 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. The same as for statistics presented in Table 2

Effect of vertical turbulent exchange
The results of 11 separate model runs with different values of min k are shown in Fig. 11. It is a  Fig. 11 demonstrates that min k is an important parameter for predicting all the presented state variables. 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. Figure 11

Effect of relaxation to temperature and salinity profiles
As it has been mentioned in Sect. 3, the model is forced by prescribed depth profiles of temperature (T ) and salinity ( S ) among the other forcing. The relaxation to the T and S profiles is necessary for 1-D simulations in an environment where lateral processes cannot be neglected (Reissmann et al., 2009). It is found that the model performance depends on the salinity relaxation time scale rather than that of temperature. All model results presented above have been calculated by applying the observation data of BED for relaxation. The best fit for oxygen is found for a relaxation time of 5 days. In order to study how the variability of T and S tracer concentrations used for relaxation will affect the oxygen dynamics in the different stations, we have applied also profiles from 3-D model simulations. In Fig. 12  From the Taylor diagram shown in Fig. 12 it can be seen that the forcing with BED data gives slightly better results. The close coincidence of the oxygen comparison symbols for BY0 and BY1 (yellow and green) points to the low sensitivity of the oxygen dynamics at these stations to the prescribed salinity field. The influence of the S T / forcing data is more pronounced for the other two stations and in particular for BY5 where 8 . 0 = σ and 95 . 0 = R in the case of the 3-D model profiles, however the agreement with observation data is rather better in the case of using the BED profiles (see Table 3). Despite the underestimation of salinity, the good results for oxygen demonstrate that it is possible to utilise 3-D model data for S T / relaxation in all cases when observation data is scarce or absent.

Effect of atmospheric forcing
In order to investigate the model sensitivity to variations of the atmospheric forcing, we present results from five different cases and compare them with the observation 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.13). The value of min k is fixed to its best fit value which is different for each particular station (see the values of min k already reported in Sect. 6.1). The close grouping of the comparison points for BY15 (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 ones 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 BY0 and BY1, even though the results for σ and S are worse for BY5. Another possible inference drawn from Fig. 13 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 surface have to be adjusted in order to parameterise lateral nutrient fluxes.
A Taylor diagram is drawn in Fig. 14  to note, 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 3.

Summary and Conclusions
In As it has been shown, these physical factors also clearly dominate the biological production and respiration of oxygen at the surface layer.
The largest mismatch with observations is found in 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 21 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 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        Table 4.