Accounting for Dissolved Organic Nutrients in an SPBEM-2 Model: Validation and Veriﬁcation

: Modern models of the Baltic Sea eutrophication describe only a bioavailable fraction of the nutrient input from land, thus introducing uncertainty into forcing. In order to alleviate this uncertainty, the coupled 3D hydrodynamical-biogeochemical St. Petersburg Eutrophication Model (SPBEM) has been expanded with variables representing dissolved organic nutrients. The model modiﬁcation involves an explicit description of the labile and refractory fractions of dissolved organic nitrogen and phosphorus, in addition to their particulate forms, represented by the detritus variables. The modiﬁed SPBEM-2 allows for a full account of the total amounts of nutrients reported in ﬁeld measurements and presented in environmental documents. Particularly, a model description of detritus, as the only bulk organic matter variable, has been replaced by more realistic parameterizations with adequate rates of settling and mineralization. The extensive validation and veriﬁcation of the model performance in the Gulf of Finland from 2009 to 2014, based on over 4000 oceanographic stations, shows that SPBEM-2 plausibly reproduces all the major large-scale features and phenomena of the ecosystem dynamics in the Gulf of Finland, especially in its surface productive layer. These demonstrated capabilities of SPBEM-2 make the model a useful tool, both in studies of biogeochemical interactions and in historical and scenario simulations.


Introduction
As has been demonstrated by the recent compilation of multiple-model computations of the Baltic Sea eutrophication [1], one of the major uncertainties in historical and scenario simulations is generated by the considerable differences in the prescribed nutrient land loads. These differences originate from the models' formulations and calibration, including their assumptions concerning the bioavailable fractions contained in the external nutrient inputs. As a result, the external nutrient load ranges between models in absolute numbers up to 50% for nitrogen and almost three-fold for phosphorus [2]. The dissolved organic matter (DOM) has for a long time been described in the box and 1D simulation models [3][4][5]. However, the limited computing resources, especially when developing three-dimensional ecosystem models, has typically constrained the total number of modeled variables. The rationale behind cutting down the list of variables by only describing the dynamics of the bioavailable fraction, which is considered to be rather easily mineralized, in these older models [2,6,7] was the presumption of the insignificant role of refractory fractions in nutrient cycling.
The comparative simulations of DOM, conducted for the surface waters in idealized marine environments [8], have shown that the photo-transformation of organic nitrogen from a stable fraction into a labile one was much higher in the estuaries (0.12 mmol N m −3 day −1 ) compared to those in the coastal areas (0.04 mmol N m −3 day −1 ). Applying these rates to corresponding areas in the Gulf of Finland and assuming a light penetration depth of 10 m, the annual "recovery" flux of bioavailable organic nitrogen can be estimated at 2000-6000 tons N year −1 , that is, according to our simulations, about 5-15% of a nitrogen input due to the nitrogen fixation [9]. Aarnos et al. [10] have experimentally estimated the nitrogen photo-transformation in the Gulf of Finland and found that it would barely exceed 38 to 50 tons N year −1 , that is, about 13% to 23% of the annual atmospheric deposition. These two-orders-of-magnitude differences in the estimated rates of nitrogen photo-transformation clearly illustrate the large uncertainties concerning the integral nitrogen fluxes, which are especially important for long-term projections.
In some form(s), DOM is included both in the regional [11][12][13][14] and global [15][16][17][18][19][20][21] simulation models. However, these representations are less commonly accepted than those which have been used for a long time, for instance, in the parameterization of the phytoplankton growth rate or zooplankton excretion [22][23][24]. Even definitions of the different fractions differ from model to model.
In a recent review [25], different fractions of DOM were categorized into five types, accordingly to a typical time scale of their existence: the labile DOM lasts for 1 day, the semi-labile lasts for 1-2 years, the semi-refractory lasts for 20 years, the refractory lasts for 16,000 years, and the ultra-refractory lasts for 30,000 years. The lifetime of different organic fractions is defined by two major mechanisms of their transformation and destruction [26][27][28]. Firstly, the biochemical destruction and eventual mineralization of the particulate and dissolved organic compounds by microorganisms is generally directed from more complex higher-molecular-weight substances to simpler lower-molecular-weight substances and inorganic molecules. The second mechanism is the photochemical transformation degrading refractory fractions into the more labile ones. However, this classification is rather far from being commonly accepted [29,30]. Therefore, we will give our reasoning and working definitions for the dissolved organic nutrient fractions considered in the present study.
In an attempt to eliminate uncertainties concerning external nutrient loads, the St. Petersburg model of eutrophication (SPBEM) has been modified with an explicit description of the total amounts of organic nutrients, including both dissolved and particulate forms [31]. In particular, besides N and P detritus variables, the dynamics of the labile and refractory fractions of dissolved organic nitrogen and phosphorus have been described, with four additional equations. The modification was developed and implemented based on numerical experiments in the Gulf of Finland, made with plausible initial and realistic boundary conditions, conducted for recent years. An explicit description of the dissolved organic nutrients made the representation of biogeochemical cycles more realistic and allowed for a full account of the nutrient input from external sources. However, splitting the total inputs originating from different sources into three fractions (particulate, dissolved labile, and dissolved refractory) still requires a prescription of fractionation coefficients, knowledge of which is, so far, very poor. At the same time, an increase in the number of model variables and the relevant processes of their interactions, while other parameterizations (including their coefficients) remained almost unchanged, did not result in a significant improvement (nor deterioration) of the model performance. Here, we present a further adjustment of those parameterizations, especially the ones describing the detritus dynamics. However, no comparison is presented in this study to the previous versions of the SPBEM family. SPBEM and SPBEM-2 are not directly comparable because of the differences in their hydrodynamic modules, sets of described variables, and structures of biogeochemical fluxes. On the other hand, a comparison to the previous simulation with SPBEM-2 [31] would merely reflect the calibration results of a rather technical nature.
Our choice of the Gulf of Finland as an object of modeling was determined on the basis of several important reasons. Firstly, this region is one of the most "loaded" and eutrophicated areas of the Baltic Sea, where the largest Baltic metropolis, St. Petersburg, is situated at the mouth of the largest Baltic river, the Neva River, both generating a very high nutrient input from the land (ca. 100,000 tons of N and 4000 tons of P annually). Therefore, it is important to have a eutrophication model of this area, as it is necessary for scenario simulations of the changing nutrient land loads. Secondly, the relatively limited size of the Gulf of Finland allows it to be covered by a high-resolution grid which remains small enough to permit multiple numerical experiments at a modest computational expense. Thirdly, the extensive database of hydrographic and hydrochemical observations allows for an extensive and confident model-data comparison.

Modification of the SPBEM-2 Model
The St. Petersburg model of eutrophication (SPBEM) implemented in this study is a further development of the biogeochemical formulations suggested in [6,32], which gradually evolved into the 3D SPBEM model [9,[33][34][35]. An updated version of the SPBEM-2 model was presented in [31], where the hydrodynamic module was based on the MIT general circulation model (MITgcm), with corresponding adjustments. In [31], the configuration of the MITgcm, used to simulate the hydro-thermodynamics of the Gulf of Finland, was described in detail, including the initial and boundary conditions. In addition, a detailed validation of the hydrodynamic module was performed. In addition to changes in the hydrodynamic module, the same study further presented a modification of the biogeochemical module of the SPBEM family, in which the labile and refractory fractions of dissolved organic nitrogen and phosphorus were introduced as new dependent variables and were described with four additional equations. Here, we further develop the dissolved organic matter sub-model and provide a full description of it.
The original formulation and initial versions of SPBEM intentionally described only the bioavailable fractions of organic nitrogen and phosphorus, both in the external inputs and cycling through the pelagic and sediment subsystems. Moreover, all the "dead" organic matter, both in its particulate and dissolved states, was assumed to be lumped together in detritus variables. Here, we retain variables representing particulate organic nutrients but expand the model with four variables representing the labile and refractory fractions of dissolved organic nitrogen (DON) and phosphorus (DOP) ( Table 1).  In the current version of SPBEM-2, the scheme of biogeochemical fluxes ( Figure 1, Table 2), described by the biogeochemical translocation functions Φ (the "right-hand sides"), in a system of advection-diffusion equations is defined as follows: The equations for the sediment variables are written as:  Table 1) and nutrient fluxes (see Table 2).

Flux Definition Units
Phototransformation of N-and P-dissolved organic refractory fractions mg N (P) m −3 day −1 Most of the parameterizations implemented in Equations (1)- (9) and (14)- (17) were drawn from our previous publications, particularly from [6]. For the sake of the readers' convenience, a full mathematical description of these formulations is repeated in Appendix A, together with new formulations for dissolved organic nutrients, where the initial forms of the formulations were taken from the carbon incarnation of BALTSEM (BALTSEM-C) [13]. Here, we discuss the assumptions behind the chosen definitions and interpret the newly introduced variables.
The new translocation functions (Equations (10)-(13)), describing biogeochemical transformations of the dissolved organic nitrogen (DON) and phosphorus (DOP), include the following sources and sinks described below. Since our original formulation [6,32,34] does not explicitly resolve diurnal processes, we do not account for dissolved organic substances whose lifetime is a day or less. We are also not interested in the geological scales of thousands of years. Therefore, only two fractions of DON and DOP were introduced into the modified SPBEM-2. The first fraction is considered to be "labile," with a characteristic lifetime of about 1 year. The sources of labile organic matter are: inputs from land, part of the excretion of heterotrophs (E N(P) ), part of mineralized detritus (W DN(P) ), and influx due to the photochemical transformation of refractory organic matter (W DRN(P) ). A sink for labile organic matter is mineralization (W DLN(P) ) into dissolved inorganic nitrogen and phosphorus.
where a mN(P) is mineralization rate and b mN(P) is a temperature constant (see Table A3).
The second fraction includes all the remaining dissolved organic compounds considered "refractory," and its transformation is determined only by the process of photochemical oxidation. In the present version of our model, we disregard the autochthonous production of refractory material considered, for instance, in the ERSEM model [36]. Thus, the only source of refractory organic nutrients here is that part of the land inputs which comes with the river runoff as humic substances and iron-humate complexes [37,38]. A small part of the refractory DON and DOP is transferred into a bioavailable labile form by the photochemical transformation (W DRN(P) ).
W DRN(P) = a uv exp(ε uv I z )RDON(P) (19) where I z is the photosynthetically active radiation on the z level (see Equation (A6)), a uv is a rate of photo-transformation, and ε uv is a factor describing faster extinction of the ultraviolet radiation with depth (see Table A3). As a result of this modification of the dissolved organic matter sub-model, the directions and interpretations of the biogeochemical fluxes have been changed, compared to [31]. In the presented version of SPBEM-2, part of the decomposed detritus is immediately mineralized to ammonium and phosphate, and part is hydrolyzed to a labile form of organic matter. Heterotroph excretion fluxes have also been revised. In the new version, some of the catabolic products are excreted in organic form and enter the pool of labile organic matter. However, we still disregard both the extracellular excretion (exudation) and possible utilization of DON and DOP by the phytoplankton.
Adding equations for the dissolved form of organic matter required a revision of the sedimentation and mineralization rates of detritus. Because in the previous model's version (SPBEM), the variables of particulate nitrogen detritus and particulate phosphorus detritus were the only form of dead organic matter, its sedimentation rate was underestimated in comparison with the observed values in order to account for the DOM implicitly included in detritus variables. The introduction of dissolved organic nutrients made it possible to adjust the detritus sinking velocity closer to the reported ranges [23,39,40], increasing it, for example, from 2.6 to 6.6 m day −1 (at 10 • C). These values may seem underestimated in comparison to the sinking velocities of diatom aggregates and fecal pellets of 50-100 m day −1 [39], which, however, are considered to be a minor component of the total downward flux of the organic matter in the Baltic Sea [41]. Such an augmentation of the detritus removal from the surface layers required a significant intensification of the detritus destruction rate, from 0.004 to 0.088 day −1 (at 10 • C). As a result, the nutrient retention capacity in the water column was maintained almost the same as in the comparably well-validated simulation in [31].
A full mathematical description of the parameterizations used in Equations (1)- (17), including all the coefficients and constants, is given in the Appendix A.

Data
The data used in this study originate mostly from the international database collected during the Gulf of Finland Year 2014 [42], which contains observations made by the Estonian, Finnish, and Russian scientific institutions and monitoring authorities for the period from 1999 to 2014. Data from occasional field surveys, conducted by other research vessels, were taken from the Baltic Environmental Database (BED, http://nest.su.se/bed). In total, data from over 4000 oceanographic stations were used. The average monthly values of river runoff and nutrient land loads were taken from the HELCOM Pollution Load Compilation database (HELCOM PLC-Water, http://nest.su.se/helcom_plc). The annual values were decomposed into monthly values using historical seasonal cycles, according to [32]. The separation into labile and refractory components was conducted in accordance with [31], assuming a 30% bioavailability of dissolved organic nitrogen and 80% bioavailability of dissolved organic phosphorus. The recent annual nitrogen atmospheric deposition was taken from [43] and decomposed into monthly values, as in [32].

Model Set-Up and Initial and Boundary Conditions
The computational grid covers the Gulf of Finland from Neva Bay to the meridian of 24.08 • E, where the open boundary is located ( Figure 2). The horizontal resolution of the spherical grid is 2 in latitude and 4 in longitude (~2 nautical miles), and the vertical resolution for the z-coordinate is 3 m from the surface to the bottom.
The initial fields of the temperature, salinity, and biogeochemical variables of the model were built from the observed vertical profiles of these oceanographic parameters, which were averaged for the winter months of 2002-2012 by interpolation into the nodes of the computational grid. The values of the sea level elevation and components of the current velocity were initially set equal to zero. The initial fields of benthic variables were taken from the final distributions, simulated from 2009 to 2014 using the previous version of SPBEM-2 and discussed in [31]. Both the initial conditions and conditions at the western boundary were obtained by processing the data from the international database [42] and BED. To resolve both the seasonal and interannual variations in the characteristics at the open boundary, the observed vertical profiles from the area between 23 • E and 24.7 • E longitude were step-wise horizontally averaged over successive 15-day time intervals. Subsequently, the averaged data were interpolated into a continuous time series. In this way, the boundary conditions were generated for the temperature, salinity, inorganic nitrogen, phosphorus, silicon, and organic matter. During the simulations, the boundary conditions for various forms of organic matter (two forms of dissolved organic nutrients, phytoplankton, and detritus) were determined as a fraction of the total organic amount, estimated from observations as the difference between the total amount of nitrogen and phosphorus and their inorganic forms. The further splitting of the organic nutrients between labile and refractory forms was made to be identical to their ratios, which were simulated in the inner grid cell nearest to the boundary.
Atmospheric forcing was set based on the ERA-Interim reanalysis fields (https://www.ecmwf.int). The model uses fields of pressure, wind speed components, air temperature, humidity, short-wave and long-wave incoming radiation, and precipitation.
Accounting for a water residence time in the Gulf of Finland of 2-3 years [32], both the pelagic and benthic fields were further adjusted during a spin-up run of three years, performed with the repeating boundary conditions for 2009. Numerical experiments were run for the period from 1 January 2009 to 31 December 2014.

Model Performance
We demonstrate the model performance by (a) validating it with a quantitative model-data comparison for several variables, and (b) verifying the internal consistency of SPBEM-2 and its capability to reproduce important feedback loops ("the vicious circle" [43]).
For validation, several methods of demonstrating the capabilities of SPBEM-2 were employed. A direct model-data comparison, with a seasonal and inter-annual resolution, was conducted for the most data-rich international monitoring stations, LL7, LL5, and LL3A ( Figure 2). Considering that the real rugged bottom topography is approximated by a two-dimensional step function on a 2 × 2 nautical miles grid, a comparison with simulations was conducted as follows. First, we considered a 3 × 3 square of grid cells, where the central node was chosen as close to the average position of the monitoring station in question as possible. Further, we selected the cell, the depth of which was as close to the depth of the monitoring station as possible, and used the simulated profiles in this cell for comparison with the observed profiles. Unfortunately, all three stations with frequent monitoring are situated in the Western Gulf of Finland. Despite a large volume of observations in the Eastern Gulf of Finland, there are no intensive monitoring stations with an amount of measurements sufficient for making a seasonal to interannual model-data comparison, as described above for the western part of the gulf. Besides, a comparison at the sparse fixed stations gives less information about the quality of the simulation of the horizontal gradients in the distribution of ecosystem components. Therefore, we made a 3D statistical inter-comparison for the surface and bottom layers. For every single observation in the database, a corresponding daily averaged simulated variable was selected for the closest grid point as a match to this observation. The resulting pair-wise datasets were used for the calculation of some simple statistics. Apparently, with a grid size of 2 × 2 nautical miles, a grid cell height of 3 m, and daily averaging, we cannot expect a reproduction of the diurnal extremes in the observed variables, neither in these datasets, nor at the monitoring stations.
The capability of the model in simulating the biotic part of the ecosystem is demonstrated through a comparison with the measurements of phytoplankton primary production (PP) performed for years by the Russian State Hydrometeorological University in the Eastern Gulf of Finland. Measurements were carried out by the light-dark (or bulk) oxygen evolution method, which has been used for almost a century [44][45][46]. The method consists of the registration of changes in oxygen concentration using the high-precision Winkler method, with a 0.1% precision in oxygen determinations [47,48], following the 24 h incubation of natural communities in clear and dark bottles. The primary production is calculated as the sum of the rate of change in the oxygen concentration in clear bottles (equal to the primary production minus respiration) and that in dark bottles, i.e., the respiration. The measurements were conducted at a network of stations ( Figure 2) from the end of July to the beginning of August. The selection of simulated values for comparison was conducted as follows. First, for every experimental estimate of PP, simulated values were compiled from a quadrate of 3 × 3 grid cells surrounding the station of measurements and for three days in row, enveloping the time of measurements. Then, the simulated phytoplankton growth rates in nitrogen units, as well as the nitrogen fixation by cyanobacteria, were summed up and converted into carbon units, with a factor of 6.625. Finally, this dataset was averaged and used in a statistical model-data comparison for all the measurements in every simulated year. At the same time, we did not consider the routine chlorophyll observations as a good measure for comparison with the simulated phytoplankton biomass, expressed in nitrogen units, when such a comparison is conducted with quota coefficients that are fixed relative to carbon and chlorophyll [49,50]. A large intrinsic seasonal and interspecies variation of these quotas [51][52][53] make this type of model-data comparison too uncertain and unconvincing [2].
In addition to presenting compiled datasets as the usual graphs (Figures 3-5), we used such simple statistics as the average, standard deviation, and coefficient of linear correlation, as well as more integral characteristics, such as the cost function (CF, Equation (18), [7]): where Mean Mod and Mean Obs are the averages of simulated and observed variables, respectively, and STD Obs is the standard deviation of observed variables.

Results and Discussion
As shown in [31], the hydrodynamic module of SPBEM-2 rather plausibly reproduces the thermohaline characteristics and water dynamics of the Gulf of Finland at the simulated seasonal and interannual scales. Therefore, this study focused solely on a validation and verification of the biogeochemical part.

Validation
As demonstrated in Figures 3-6 and Table A4, SPBEM-2 quite realistically reproduces the nutrient and oxygen dynamics in the surface layer. The lowest values of the coefficients of the linear correlation in the surface layer were found for the total nitrogen (TN) at the LL5 and LL7 stations, which were 0.38 and 0.16, respectively. At the same time, the model-data differences in the overall averages of TN did not exceed a few percent, while the difference in standard deviations was much smaller than the standard deviation of the observations (Table A4). As seen in Figures 5 and 6, the low values of the correlation coefficients can be explained by some biases between the simulated and observed timing of extremes. Such biases could have a different origin, including analytical problems with the measurements of the total nitrogen [54] and inaccuracies in the prescribed boundary conditions, both in the land loads and at the open boundary. However, these differences do not significantly affect the general conclusion concerning the plausibility of the model performance. The usefulness of combining the usual graphical presentation and simple statistical estimates for the analysis of the model's skill should also be noted. On the other hand, the highest correlation coefficients belonged to a comparison of oxidized forms of nitrogen (NO) in the surface layer at LL3A and LL7, reaching values of 0.8 and 0.9, respectively. The differences in the mean values and standard deviations varied from rather insignificant at LL3A to almost non-existent at LL5 and LL7 ( Figure 6; Table A4). The graphical comparison allows for a further refinement of such a conclusion, with respect both to the timing and concentration of the winter nutrient maxima (Figures 3-5).
In the bottom layers, the results of the model-data comparison were mostly determined by the overestimated oxygen concentrations in the model (Figure 6, see also [31]), especially during the summer minima. The discrepancy increased with the depth of the compared monitoring stations (cf. Figures 3-5), reaching, on average, 1.7 mg O 2 L −1 at LL7. Additionally, in the model, the hypoxic conditions (O 2 < 2.86 mg O 2 L −1 ) were found less frequently than in observations. One of the evident reasons for such a general deviation is that the smoothed model bottom is shallower than the depth of the monitoring stations that, once upon a time, had been intentionally chosen to represent the deepest point in the area. For instance, at LL7, we compare measurements from samples taken at 80-90 m horizons to simulated concentrations from the model bottom layer at 73.5 m, and such a depth difference is important within the marked vertical gradients of the oxygen concentration (see also integral model-data comparison in Section 3.2). Another possible reason could also be a larger retention of nutrients in the surface layer, with a lesser detritus transport to the deep sediments and, correspondingly, a reduced biochemical oxygen demand (BOD). However, for the time being, we have postponed a corresponding enquiry and refining until a further model development, particularly of its sediment module, has been achieved.
Generally, the maximal discrepancies, both graphical and statistical, were found in a comparison of the deep layer phosphates and, consequently, the total phosphorus dynamics (Figures 3-6; Table A4). The simulated concentrations in the bottom layer are systematically underestimated. This underestimation is a direct consequence of the underestimated flux from the bottom, which, in accordance with our parameterization (A61-A63), is substantially reduced because of the overestimated oxygen concentration. On the contrary, the simulated dynamics of the NO variable did measure up to observation substantially better, even despite some mismatch in the extrema, formally indicated by the practically non-existent linear correlation between the simulated variables and observations (Table A4), particularly at LL3A and LL5. Additionally, an overestimation of the mean values of NO at the deeper stations, LL5 and, especially, LL7, might indicate a somewhat underestimated rate of denitrification as a result of the overestimated availability of oxygen for the oxidation of organic matter.
The assessment of the model capability to reproduce the seasonal spatial-temporal dynamics, in general, for the entire Gulf of Finland was made by comparing six-year monthly averages and standard deviations of all the available measurements and, corresponding to them, simulated values extracted and compiled from 3D fields, as explained in Section 2.4. The results of this model-data comparison are presented in Figure 7 for the surface layer. For the bottom layer (Figure 8), we used only oceanographic stations deeper than 20 m, that is, deeper than the average depth of the seasonal thermocline in order to prevent the inclusion of the waters occurring within the seasonal mixed layer into the comparison. December is excluded from this comparison because only two oceanographic stations were available for this month during the entire period of 2009-2014.
These comparisons further confirm the main conclusion that could be made from an analysis of the monitoring stations: so far, the SPBEM-2 better reproduces the nutrient and oxygen dynamics in the surface layers (Figure 7) compared to the deep layers ( Figure 8). Most discrepancies occurring between the measured and simulated inorganic nutrient concentrations were almost quantitatively equal to those discrepancies found for the total nitrogen and phosphorus concentrations, which indicates a very good reproducibility of the organic nutrient concentrations. In the surface layer, both the winter maxima and summer minima values are quite realistic. However, the noticeable discrepancies between the observation and simulation for inorganic phosphorus and nitrogen in April and May indicate a delayed phytoplankton spring bloom in the model. In the deep layers, the overestimated oxygen concentration results in under and overestimated concentrations of phosphate and inorganic oxidized nitrogen, respectively.

Verification
To demonstrate the causal relationships between oxygen availability and nutrient redox alterations, we extracted and presented both simulation and observation results as time-depth contour plots, built on all the vertical profiles horizontally averaged over the entire modeled domain (Figures 2 and 9, cf. also Figure 2 in [55]). Among the other expected patterns of the hypothetical "vicious circle of the Baltic Sea eutrophication" [43,54,56], the expansion of the long-lasting oxygen deficit and hypoxia in 2009-2010, which extended over rather shallow and voluminous areas and were interrupted by the intensive vertical mixing in winter 2011, looks especially conspicuous. The mixing brought up phosphate reserves accumulated in the oxygen-deficient areas, thus increasing the winter surface DIP concentrations far above the average values (see also , while the DIN concentration remained at a more typical level. As a result, the winter (March) DIN:DIP ratio dropped to rather low values, with corresponding consequences for the expansion and dynamics of both the nitrogen limitation and nitrogen fixation by cyanobacteria ( Figure 10) and, consequently, for the phytoplankton primary production (Figure 11, Table 3).    Table 3. Comparison of the simulated primary production (g C m −2 day −1 ) with estimates obtained using the bulk light-dark oxygen method at the end of July until the beginning of August in the Eastern Gulf of Finland (positions of the stations with measurements are indicated in Figure 2). Primary Production, g C m −2 day −1

Observations Model
Year As is well known, from an analysis of both observations [57][58][59][60][61][62] and simulations [33,34,63], and as demonstrated in Figures 10A and 12, the phosphorus limitation in Neva Bay and some other estuaries is gradually replaced by the nitrogen limitation offshore. This knowledge and these modelling results make questionable the conclusion concerning the present-day nitrogen limitation of these areas and the emergence of a phosphorus limitation only in the distant future (cf. Figure 6 in [64]). However, such an average distribution can be substantially altered by the nutrient exchange with the Northern Baltic Proper [55] and internal biogeochemical interactions. Similar to the succession of events in 1995-1997 [65,66], the extensive hypoxia of 2009-2010 resulted in an expansion of the nitrogen-limited surface waters, with very elevated stocks of the winter Redfield excess of inorganic phosphorus (calculated in weight units as eDIP = DIP -DIN/7; see e.g., [54] and references therein) in 2011 ( Figure 10B, right). Consequently, the simulated rate of nitrogen fixation by cyanobacteria in 2011 increased many-fold, while the area of elevated rates substantially expanded ( Figure 10C). It should also be noted that the larger-than-usual expansion of cyanobacteria accumulations over the Gulf of Finland in 2011, compared to 1979-2014 summer distributions, is unequivocally found also in the multi-decadal time-series of the satellite-detected accumulations of cyanobacteria [67]. This massive contribution of cyanobacteria to the phytoplankton primary production in 2011 is also quite clearly seen in comparison with the field estimates ( Figure 11, Table 3). Even in the Eastern Gulf of Finland, where the nitrogen fixation rates were lower than in the rest of the Gulf of Finland ( Figure 10C), the simulated primary production was almost two-fold higher than the one estimated by the light-dark oxygen evolution method. Most likely, the measured primary production is underestimated because the filamentous cyanobacteria could hardly be caught in the incubation bottles adequately. On the other hand, in other years, the simulated phytoplankton primary production was expected to be somewhat smaller than the measured production because the phytoplankton growth rate, simulated in nitrogen units, is by definition equal to the net primary production, while all the measurements represent some values intermediate between the gross and net primary production [68].
Considering the results of the validation and verification in general, it already appears that the presented version of SPBEM-2 is capable of plausibly reproducing all the major large-scale features and phenomena of the nutrient dynamics in the Gulf of Finland, especially in its surface productive layer that, for hypsographic reasons, contain and transform the major pools of nutrients. The expansion of SPBEM-2 with dissolved organic nutrients allows for a full account of the land loads in both historical and scenario simulations, thus reducing the uncertainty in forcing.

1.
A modified version of the coupled 3D hydrodynamical-biogeochemical St. Petersburg Eutrophication Model (SPBEM-2) has been presented, including a full description of the biogeochemical module, with all the equations and coefficients.

2.
The SPBEM-2 has been expanded with four new variables (the labile and refractory fractions of the dissolved organic nitrogen and phosphorus) and a redefinition of the detritus variables, allowing for a full account of the total amounts of nutrients as they are reported in field measurements and presented in environmental documents, thus alleviating one of the important sources of uncertainty in the boundary conditions. 3.
An extensive model-data comparison of the model performance in the Gulf of Finland from 2009 to 2014, based on over 4000 oceanographic stations, was conducted, showing that SPBEM-2 is capable of plausibly reproducing all the major large-scale features and phenomena of the seasonal and interannual nutrient dynamics in the Gulf of Finland, especially in its surface productive layer.

4.
Validation and verification of the capabilities of SPBEM-2 were conducted, showing the model to be a useful instrument in studies of the biogeochemical interactions determining the state of eutrophication and a necessary tool in both historical and scenario simulations. 5.
From the results and experience obtained from the performed model-data comparison, it was shown that it is very important to combine both formal (statistical) and expert (graphical) means of comparison in order to properly reveal and explain both the realistic model behavior and its deficiencies. The net primary production PP i (mg N m −3 day −1 ) by the i-th (i = 1,2,3) autotroph group (Table 1) is: where T ( • C) is the water temperature and I z (W m −2 ) is the mean daily integral photosynthetically active radiation on the z level: which is calculated from its value at the upper surface I sur and the light extinction ε due to the water, salinity, heterotrophs, autotrophs, and detritus: The simple parameterization of the reversible ammonium-induced inhibition of the nitrate uptake, shown in Equation (A8), accounts for the fast (hours, i.e., "instant" at the model time scale) synthesis of nitrate reductase and produces a sharper transition between the inhibited and uninhibited nitrate uptake than many other parameterizations do.
The nitrogen fixation F NF (mg N m −3 day −1 ) by cyanobacteria (i = 1) is defined by the water temperature, ambient N/P ratio, and phosphate concentration: The temperature-dependent mortality M i and sinking S i (mg N m −3 day −1 ) of cyanobacteria (i = 1): and other autotroph groups (i = 2, 3): are increased by γ + 1 times from non-limitation (Λ i = 1) to the highest limitation (Λ i = 0). The consumption G k (mg N (P) m −3 day −1 ): is defined by the filtration rate: where f is the concentration of available food: The consumed food is split into assimilated K k : and unassimilated U k : The excretion of ammonium E N (mg N m −3 day −1 ) and phosphate E P (mg P m −3 day −1 ) by heterotrophs is related to the assimilation and accounts for the difference between the heterotrophs' stoichiometry from that of the food: The mortality of heterotrophs M Z (mg N m −3 day −1 ), the closure term in the model, is density-dependent: The destruction of the detritus W DN(DP) variables (mg N (P) m −3 day −1 ) and dissolution of detritus silica W DS (mg Si m −3 day −1 ) are described as temperature-dependent, with first-order kinetics: The mineralization of labile dissolved organic matter, i.e., LDON and LDOP, is also given as: Similarly, the mineralization of sediment variables B k (mg N (P)(Si) m −3 day −1 ) is described as: The process of photo-transformation (mg N (P) m −3 day −1 ), which turns a refractory fraction of dissolved organic matter into a labile fraction, is described as: The temperature-dependent sinking of the detritus variables S k (mg N (P) (Si) m −3 day −1 ): is modulated by the contribution of their major sources into the total flux c D = M i + M Z + U N , whereas the numerical values of the algal sinking velocities a si are used as a dimensionless weight factor. The ammonium V N , nitrate V O , phosphate V P , and silica V S uptakes by autotrophs (mg N (P) (Si) m −3 day −1 ) are defined by the rates of primary production and nitrogen fixation: The nitrification W O is both photo-inhibited at high light intensities and quickly decreases to zero at very low oxygen concentrations: The pelagic denitrification W R starts sharply at a "threshold" oxygen concentration O crd , ends in an anoxic condition, and its rate depends on the nitrate concentration: The sediment mineralization fluxes (Equations (A35)-(A37)) are split into several pathways, and the proportions depend on the water oxygen concentrations: The sediment nutrients B k are buried B Bk (mg N (P) (Si) m −3 day −1 ) at a constant rate b ur : The oxygen production and consumption fluxes are proportional to the corresponding nitrogen fluxes.