A novel acclimative biogeochemical model and its implementation to the southern North Sea

Ecosystem models often rely on heuristic descriptions of autotroph growth that fail to reproduce various stationary and dynamic states of phytoplankton cellular composition observed in laboratory experiments. Here, we present the integration of an advanced phytoplankton growth model within a coupled 3-dimensional physical-biogeochemical model, and the implementation of the model system to the Southern North Sea (SNS) defined on a relatively high resolution (∼ 1.5-4.5 km) curvilinear grid. The autotrophic growth model, recently introduced by Wirtz and Kerimoglu (2016), is built up on a set of 5 novel concepts for the allocation of internal resources and operation of cellular metabolism. The coupled model system consists of the general estuarine transport model (GETM) as the hydrodynamical driver, a lower trophic level model and a simple sediment diagenesis model. We force the model system with realistic atmospheric and riverine fluxes, background turbidity caused by suspended particulate matter and open ocean boundary conditions. For a simulation for the period 2000-2010, we show that the model system satisfactorily reproduces the physical and biogeochemical states of the system, as inferred from 10 comparisons against data from long-term monitoring stations, sparse measurements, continuous transects, and remote sensing data. In particular, the model shows high skill both in coastal and off shore waters, and captures the steep gradients in nutrient and chlorophyll concentrations observed prevalently across the coastal transition zone. We show that the cellular chlorophyll to carbon ratio show significant seasonal and lateral variability, the latter amplifying the steepness of the transitional chlorophyll gradient, thus, pointing to the relevance of resolving the physiological acclimation processes for an accurate description of 15 biogeochemical fluxes. Wirtz, K.W. and Kerimoglu, O.: Optimality and variable co-limitation controls autotrophic stoichiometry, Frontiers in Ecology and Evolution, doi:10.3389/fevo.2016.00131, 2016.


Introduction
Modelling the biogeochemistry of coastal and shelf systems requires the representation of a multitude of interacting processes, not only within the water but also at the adjacent earth system components such as the atmosphere (e.g., nitrogen deposition), land (e.g., rivers), sediment (e.g., diagenetic processes), and biochemical processes in water (see., e.g., Cloern et al., 2014;Emeis et al., 2015).For being able to reproduce the large scale spatial and temporal distribution of biogeochemical variables in coastal systems, a realistic representation of hydrodynamical processes is often critically important, at least those relevant to coarse grids were employed (≥ 7 km).While showing good skill in reproducing off-shore dynamics, these models seemed to have a relatively limited performance at the shallow, near-coast regions (when reported).The BLOOM-Delft3D (Los et al., 2008) on the other hand, is one of the rare examples with a finer grid (down to 1 km at the Dutch coasts) at the cost of a relatively smaller domain, similar to ours.Although this model system performs decently at both coastal and off-shore areas, its performance within the German Bight has not been fully assessed.Moreover, none of these models provide elaborate descriptions of the stoichiometric regulation of autotrophs, as mentioned above.Therefore, our new model system is expected to fill two important gaps by; 1) exemplifying for the first time to the best of our knowledge, implementation of a highly complex phytoplankton growth model at an ecosystem scale, coupled to a hydrodynamic model and other biogeochemical compartments; 2) establishing the capacity to reproduce the biogeochemistry of the German Bight both at coastal and offshore regions with a single parameterization and model setup.
For a 11 year hindcast simulation of the period 2000-2010, we show that the model can adequately capture the spatiotemporal variability of the physical and biogeochemical features of the SNS based on comparisons against various data sources.
Importantly, the model can reproduce the steep chlorophyll and nutrient gradients prevalently observed across the Waddensea-German Bight continuum.We show that the chlorophyll gradients are linked with nutrient, hence, productivity gradients, but also further amplified by the high chlorophyll to carbon ratios at the shallower regions owed to the high turbidity.

Data
Data from monitoring stations all reflect surface measurements.Extensive analyses of the data from Helgoland Roads have been provided by Wiltshire et al. (2008) and from Sylt by Loebl et al. (2007).Temperature, salinity, dissolved inorganic nitrogen Biogeosciences Discuss., doi:10.5194/bg-2017Discuss., doi:10.5194/bg- -104, 2017 Manuscript under review for journal Biogeosciences Discussion started: 30 March 2017 c Author(s) 2017.CC-BY 3.0 License.and phosphorus data obtained from the International Council for the Exploration of the Sea (ICES, www.ices.dk)were used for the validation of the physical and biogeochemical model, by means of point-wise comparisons within the surface and bottom layers, i.e., upper and lower 5 meters.
Vertically resolved Scanfish data and continuous Ferrybox measurements were all gathered within the Coastal Observing System for Northern and Arctic Seas (COSYNA, Baschek et al., 2016, and references therein).Satellite data employed here are provided by the European Space Agency (ESA), Ocean Color-Climate Change Initiative version 2.0, where the NASA OC4.V6 algorithm had been applied to the MERIS, MODIS and SeaWiFS products for the estimation of chlorophyll concentrations (Grant et al., 2015).

Model
Major processes taken into account by the model are the lower trophic food web dynamics, phytoplankton ecophysiology and basic biogeochemical transformations in water, and the transformation of N-and P-species in benthos (Fig. 2 and Section 2.2.1).Physical processes are resolved by the coupled 3-D hydrodynamical model, GETM (Section 2.2.2).Turbidity caused by suspended particulate matter (SPM), nutrient loading by rivers and atmospheric nitrogen deposition were considered as model forcing (Section 2.2.3).The model grid and riverine fluxes considered in this study are shown in Fig. 1.

Biogeochemical model
The pelagic module, the Model for Adaptive Ecosystems in Coastal Seas (MAECS), is a lower trophic level model that resolves cycling of carbon, nitrogen (N) and phosphorus (P), and importantly, acclimation processes of phytoplankton, a detailed description of which is provided by Wirtz and Kerimoglu (2016) and in the supplementary material (A).The acclimation module features a whole set of physiological traits (x), which control accessory and assimilation of multiple resources in autotrophs.
These comprise, affinity for DIN and DIP, protein pools invested to nutrient uptake and light harvesting, and specific P-and N-uptake activities.Their instantaneous or transitory acclimation follows an extended optimality principle, where δ x corresponds to a flexibility constant, and the two terms in brackets describe the direct effects of trait changes on the specific phytoplankton growth rate V C (in units of cellular C) and the indirect effects through changes in the Chl:C:N:P stoichiometry, expressed by the quotas q, respectively (Wirtz and Kerimoglu, 2016).As a result of trait variations formulated in Eq.1, Chl:C:N:P stoichiometry is continuously varied depending on ambient light and nutrient conditions and on the metabolic demands of autotrophic cells.
Other components of the pelagic module are similar to standard descriptions in state-of-the-art ecosystem models.Phytoplankton take up nutrients in the form of dissolved inorganic material (DIM).Losses of phytoplankton (B) and zooplankton (Z) due to mortality are added to the particulate organic matter (POM) pool, which degrades into dissolved organic material (DOM), before becoming again DIM and closing the cycle.The benthic module describes only the dynamics of macronutrients N and P. Degradation of OM to DIM is described as a one step, first order reaction.Denitrification is described as a proportion of POM degradation, limited by DIN and dissolved oxygen (DO) availability in benthos.As DO is not directly modeled, it is estimated from temperature in order to mimic the seasonality of the hypoxia-driven denitrification.The model accounts for the sorption-desorption dynamics of phosphorus as an instantaneous process also as a function of temperature based on the correlation observed in the field (Jensen et al., 1995).
Further details are provided in Appendix A.

Hydrodynamic model and model coupling
The General Estuarine Transport Model (GETM) was used to calculate various hydrodynamic processes, as well as the transport of the biogeochemical variables.A detailed description of GETM is provided by Burchard and Bolding (2002); Stips et al. (2004).GETM utilizes the turbulence library of the General Ocean Turbulence Model (GOTM) to resolve vertical mixing of 10 density and momentum profiles with a k-ε two equation model (Burchard et al., 2006).GETM was run in baroclinic mode, drying of cells at the Wadden Sea.We used 20 terrain-following layers and a curvilinear grid of 144x98 horizontal cells with a horizontal resolution of approximately 1.5 km at the south-east corner and 4.5 km at the north-west corner (Fig. 1).The curvilinear grid focuses on the German Bight, and roughly follows the coastline (Fig. 1) for an optimal representation of alongand across-shore processes.Similar gridding strategies were applied successfully in other coastal setups with the GETM model (Hofmeister et al., 2013;Hetzel et al., 2015).We employed integration time steps of 5 and 360 seconds for the 2-D and 3-D processes, respectively.
Integration of model forcing was realized through the Modular System for Shelves and Coasts (MOSSCO, http://www.mossco.de),which among others, provides standardized data representations.Meteorological forcing originated from an hourlyresolution hindcast by COSMO-CLM (Geyer, 2014).Boundary conditions for surface elevations and currents are extracted from an hourly resolution hindcast by TRIM-NP (Weisse et al., 2015).For temperature and salinity, daily climatologies from HAMSOM (Meyer et al., 2011) are used, all of which are available through coastDat http://www.coastdat.de.
Two-way coupling of the biological model with GETM was achieved via the Framework for Aquatic Biogeochemical Models (FABM, Bruggeman and Bolding, 2014) as one of the coupling standards adopted in MOSSCO.The pelagic module is defined in the 3-D grid of the hydrodynamic model, whereas the benthic module is defined in 0-D boxes for each water column across the lateral grid of the model domain (Fig. 1).Each benthic box interacts with the bottom-most pelagic box of the corresponding water column in terms of a uni-directional flux of POM from the pelagic to the benthic states, and a bi-directional flux of DIM depending on the concentration gradients.
For the integration of the source terms, a fourth order explicit Runge-Kutta scheme was used with an integration time step of 360 seconds, as for the 3-D fields in GETM.Exchange between pelagic and benthic variables was integrated with a first order explicit scheme at a time step identical to that of the biological model.

Model forcing and boundary conditions
Light extinction is described according to: where, I 0 is the photosynthetically available radiation at the water surface, and the first and second terms describe the attenuation at the red and blue-green portions of the spectrum.We assume that the partitioning of the two (a) and the attenuation length scale of the red light (η 1 ) are constant over space and time as in Burchard et al. (2006), and that the attenuation of blue-green light is due to SPM (as described by η 2 ) and organic matter (sum term).We chose a = 0.58 and η 1 = 0.35, which correspond to Jerlov class-I type water, thus clear water conditions (Paulson and Simpson, 1977), given that the attenuation by SPM and organic matter is explicitly taken into account.For calculating attenuation due to SPM, a daily climatology of SPM concentrations defined over the model domain was utilized, like in ECOHAM (Große et al., 2016).The SPM field was constructed by multiple linear regression of salinity, tidal current speed and depth for each Julian day (Heath et al., 2002).
Then, η 2 , or the inverse of SPM caused attenuation coefficient was calculated according to: where, the attenuation for background turbidity, K w = 0.16 m −1 and specific attenuation coefficient for SPM, SP M = 0.02 m 2 g −1 according to Tian et al. (2009).For calculating the attenuation due to organic matter in Eq.( 2), phytoplankton, POC and DOC were considered (Table A3).
Freshwater and nutrient influxes were resolved for eleven major rivers along the German, Dutch, Belgian and British coasts (Fig. 1).For eight of these rivers, Radach and Pätsch (2007) and Pätsch and Lenhart (2011) present a detailed quantitative analysis of nutrient fluxes.Besides the fluxes in inorganic form based on direct measurements, fluxes in organic form have been accounted for, first by calculating the total organic material concentration by subtracting dissolved nutrient concentrations from total nitrogen and total phosphorus, then by assuming 30 % of the organic material to be in particulate form (i.e., POM; Amann et al., 2012).Further, 20 % of POM is assumed to describe phytoplankton biomass (Brockmann, 1994), C:N:P ratio of which was assumed to be in Redfield proportions.Finally, no estuarine retention/enrichment was assumed, following Dähnke et al. (2008).All river data except for the river Eider were available in daily resolution, however with gaps.Short gaps (<28 assumed to be at zero-gradient.For DIM, DOM and POM, monthly values of ECOHAM (Große et al., 2016), interpolated to 5m depth intervals are used as clamped boundary conditions.

Quantification of Model Performance
For the comparisons with the station data for DIN, DIP and Chl concentrations, Pearson correlation coefficients were calculated for all temporal matches.For the evaluation of model performance against the ICES data for temperature, salinity, DIN and DIP, correlation scores and model standard deviations normalized to measured standard deviations are displayed as Taylor diagrams, where the correlation score and the normalized standard deviation correspond to the angle and distance to the center (Jolliff et al., 2009).For this purpose, temporal matching was identified at daily resolution, vertical matching were obtained by comparing the measurements within the upper 5 meters from the sea surface and within the 5 meters above the sea floor with the model estimates at the top-most and bottom-most layers, and finally lateral matching by calculating the average of the values from four cells surrounding the measurement location, inversely weighted with respect to the Cartesian distance.A comparison of simulated salinities with the FerryBox measurements along the cruise between Cuxhaven (at the mouth of river Elbe) and Immingham (at the mouth of river Humber) (Petersen, 2014), demonstrates that the model captures the general salinity patterns (Fig. 3).However, the freshwater plume in the German Bight as simulated by GETM seems to extend further from the coast than observed, which suggests moderate over-estimation of horizontal mixing by the model.given the proximity of these stations to major rivers Elbe, North Sea Canal and Rhine (Fig. 1), these mismatches are likely to be related with our assumption that the non-dissolved fractions of the total nitrogen and total phosphorus to be all in labile form (Sect. 2.2.3).iii) Timing of spring blooms, nutrient draw-down and regeneration are mostly well reproduced, as also reflected by high correlation coefficients in general.Earlier replenishment of phosphorus relative to nitrogen is often reproduced, although with delays in some coastal stations like Norderney, which probably reflects the oversimplification of the benthic processes with respect to the description of oxygen-driven iron-phosphorus complexation kinetics (Sect.A), which has been suggested to be the main driver for the phenomena in the coastal areas (Jensen et al., 1995;van Beusekom et al., 1999;Grunwald et al., 2010).Finally, measurements from S. Amrum suggest that the classical summer-low, winter-high phosphorus pattern, as also predicted by our model in general, is entirely reversed, calling for a more detailed investigation.
There are a number of caveats when using coastal time-series data for model validation, which start from the sampling problem: short term fluctuations common for near-shore waters can only be tentatively represented by measurement frequencies of several weeks or months.For the special case of Sylt, for example, data reflect the average monthly concentrations, which naturally smooths out the short-lived blooms.However, we acknowledge the potentially inadequate description of certain processes that might have led to the overestimation of chlorophyll concentrations at Sylt, such as the grazing formulation of zooplankton and the representation of the light climate.In reality, effects of temperature on mesozooplankton occurs through phenological shifts (e.g., Greve et al., 2004) that might have a determining role on the maximum chlorophyll concentrations (van Beusekom et al., 2009), which can probably be only partially reflected by the simple Q10 rule we applied for grazing     are in some cases in the form of higher concentrations below the pycnocline but in some others, appear as thin layers around the pycnocline.While deep chlorophyll maxima are visible in stratified waters that occur in deeper regions, well-mixed shallower regions mostly show vertically homogeneously distributed high chlorophyll concentrations.The MAECS simulation agrees qualitatively very well with these patterns and captures the spatial variability of the observed vertical chlorophyll distribution (Fig. 9).Former 3-D modeling studies, such as that of van Leeuwen et al. (2013), apart from capturing the presence of a deep chlorophyll maximum, were not able to reproduce the rich variability revealed by the observations.Our model-based analysis indicates that the formation and maintenance of such structures are critically dependent on the parametrization of the underwater light climate and sinking rate of phytoplankton.Sinking speed of algae in the MAECS is inversely related to the nutrient quota of the cells, which mimics the internal buoyancy regulation ability of algae depending on internal nutrient reserves (see section A1) but also indirectly emulates chemotactic migration as typical for dinoflagellates (Durham and Stocker, 2012).The critical dependence of the formation and maintenance of vertical chlorophyll structures on the functional representation of sinking underlines the relevance of an accurate description of the intracellular regulation of nutrient storages and pigmentory material.

Coastal Gradients
Temperature stratification is one of the key drivers of biogeochemical processes through its determining role on the resource environment, i.e., light and nutrient availability experienced by the primary producers.The comparison against Scanfish transect (Fig. 9) indicated that our model can capture density stratification quite accurately.Using the temperature difference between surface and bottom layers as an indicator of temperature stratification (Schrum et al., 2003;Holt and Umlauf, 2008;van Leeuwen et al., 2015), and using monthly averages across all simulated years (2000-2010), we illustrate the areal extent and seasonality of stratification within the SNS in Fig. 10.This analysis suggests that a large portion of the model domain deeper than ∼30 m becomes stratified from April to September, with maximum intensity and areal coverage in July.The areal extent and seasonality of stratification is in agreement with those reported by earlier studies (Schrum et al., 2003;van Leeuwen et al., 2015).
Average winter concentrations of surface DIN and DIP for the entire simulation period display steep gradients along the coasts of the German Bight (Fig. 11), in line with observations (e.g., Brockmann et al., 1999).Given the riverine nutrient fluxes and uninterrupted nutrient supply from the bottom layers owed to the lack of stratification (Fig. 10), it is not surprising that the surface nutrient concentrations are higher at the coasts, but the persistence of the phenomenon in the Wadden Sea regions that are not in close proximity to the riverine fluxes and the steepness of these gradients are not intuitively predictable.Steep cross-shore nutrient gradients were partially explained by an interplay between density gradients and tidal mixing: during tidal flooding, high-salinity, therefore denser off-shore water sinks below the low-salinity, therefore lighter coastal water, thereby pushing the nutrient-rich bottom waters towards the coast (Ebenhöh, 2004;Burchard et al., 2008;Flöser et al., 2011;Hofmeister et al., 2016).The relevant physical processes, i.e. tidal assymmetries in currents and mixing under coastal density gradients, and the accumulation of nutrients at the bottom layers in off-shore waters are represented in the current model framework, so that this residual density-driven circulation mechanism is likely to be responsible for the emergence of steep nutrient gradients  shown in Fig. 11.Productivity-enhanced aggregation of particulate organic matter with suspended sediments, and hence, higher sedimentation velocities within the coastal transition zone has recently been proposed to be another contributing factor for the maintenance of the coastal nutrient gradients (Maerz et al., 2016).This process has not been accounted for by our model, and representation of such a mechanism would necessitate explicit descriptions of mineral-POM interactions and variable sinking rate as a function of particle composition, structure and size (Maerz et al., 2011).An exhaustive elaboration of the mechanisms for the maintenance of such gradients, their regional variability and steepness would be out of the scope of the current work, but constitutes a potential research goal for the future.Both the satellite (ESA-CCI) images and our model estimates, averaged over the years 2008-2010 for the two halves of the growing season, suggest much higher concentrations within a thin coastal stripe relative to the off-shore concentrations.The large scale agreement in coastal gradients result in high correlation coefficients, especially during summer (Fig. 13).For the first half of the growing season, the higher range of the modelled chlorophyll values exceed those of the ESA-CCI (Fig. 13), which seems to be caused by underestimation by the satellite product, suggested by the fact that in-situ concentrations frequently reach well over 50 mg/m 3 at the coastal stations in the German Bight (Fig. 5).These lateral gradients in chlorophyll concentrations overlap with the nutrient (Fig. 11), hence, productivity gradients (not shown).
According to our simulation results of the year 2010, average Chl:C ratio displays considerable spatio-temporal variability, even when the seasonal averages are considered, i.e., omitting short-term variability (Fig. 14).Chl:C ratio is in general higher at the coasts than at off-shore.This pattern has been previously identified based on monitoring data by Alvarez-Fernandez and Riegman (2014), and reflects the photoacclimative response to stronger light limitation at the coasts, manifested by both higher organic matter (not shown, but see Fig.  wave radiation at the water surface and increased turbidity due to stronger vertical mixing near the coast.A similar seasonal amplitude in CHL:C has been found by Llewellyn et al. (2005) for the English Channel.Slightly higher Chl:C ratios during the first half of the growing season compared to the second half are likely due to lower nutrient concentrations during the second half, which results in larger investments in nutrient harvesting in the expense of light-harvesting machinery (see, e.g., Geider et al., 1997;Pahlow, 2005;Wirtz and Kerimoglu, 2016).The modeled spatio-temporal differences in Chl:C ratios reach 5 to about three fold between different seasons of the year and between off-shore and coastal areas.The latter indicates that the differential acclimative state of phytoplankton cells amplify the steepness of the chlorophyll gradients across the coastal transition shown in Fig. 13, which, as mentioned above, seem to be driven mainly by nutrient gradients.This is significant, and calls for further attention, especially given that many modelling schemes applied at ecosystem-scales do not consider photoacclimation processes.As Behrenfeld et al. (2015) recently pointed out from a global perspective, satellite-based primary

Conclusions
In this study, we described the implementation of a coupled physical-biogeochemical model to the Southern North Sea (SNS) and analyzed the model results in comparison to a large collection of in-situ and remote sensing data.The model system accounts for key coastal processes, such as the forcing by local atmospheric conditions, riverine loadings of inorganic and organic material, atmospheric nitrogen deposition, spatio-temporal variations in the underwater light climate, major benthic processes and nutrient concentrations at open boundaries, and importantly, it hosts a novel model of phytoplankton growth, which replaces otherwise heuristic formulations of photosynthesis and nutrient uptake with mechanistically sound ones (Wirtz and Kerimoglu, 2016).Based on comparisons with a number of data sources, we conclude that the model system can produce a realistic decadal hindcast of the SNS for the period 2000-2010, in terms of both the temporal and spatial distribution of key ecosystem variables, as well as a large area of validity, i.e., both in coastal and off-shore regions of the German Bight.
We emphasize that even the phytoplankton concentrations are generally well captured by our model, considering the systematic difficulties in reproducing chlorophyll concentrations by ecosystem models in general (see, e.g., Radach and Moll, 2006).This is noteworthy, given that phytoplankton is represented by a single species in our model, whereas in reality the phytoplankton composition displays systematic shifts throughout the season: the early spring composition is usually dominated by diatoms while other species or groups usually become more abundant later in the year, such as Phaeocystis at coastal regions (eg., van Beusekom et al., 2009) and dinoflagellates off-shore (Freund et al., 2012;Wollschläger et al., 2015).We argue that the ability of our model to capture both the spring and summer is, to a great extent, owed to the fact that photoacclimation and optimality in nutrient uptake processes were accounted for.In reality, environmental change, e.g., improvement of light conditions and depletion of nutrient concentrations from spring to summer, promotes the species which have more suitable traits, e.g., regarding light and nutrient utilization.In 3-D model applications so far, photoacclimation of phytoplankton has been either ignored altogether, or it was accounted for in a heuristic sense, where the change in Chl:C ratio is described based on an empirical relationships (Blackford et al., 2004;Fennel et al., 2006).In our model, adaptation of the phytoplankton community to the light and nutrient environment is represented by dynamically changing and instantaneously optimized trait values as described extensively by (Wirtz and Kerimoglu, 2016).Other potentially relevant selection factors, such as the changes in the community structure of heterotrophic grazers (Alvarez-Fernandez et al., 2012;Löder et al., 2012;Beaugrand et al., 2014) or limitation of diatoms by silicate (Loebl et al., 2009)  Our findings suggest that the steep chlorophyll gradients across the coastal transition zone is mainly driven by the nutrient gradients, but amplified by the higher Chl:C ratios at the coastal waters.The large variations in simulated Chl:C ratios within the SNS, both in a space and time, indicate that ignorance of photoacclimation can lead to potentially flawed estimates for primary production or phytoplankton biomass as was recently pointed out by Arteaga et al. (2014) and Behrenfeld et al. (2015), who used photoacclimation schemes to derive Chl:C ratios at global scales.Here we show that such considerations apply also at coastal environments, which may be critical, given the increasing recognition of the role of coastal-shelf systems in the global carbon and nutrient cycling (Fennel, 2010;Bauer et al., 2013).
Appendix A: Detailed Model Description

A1 Pelagic Module
Local source-sink terms for all dynamic variables, functional description of processes and relationships between quantities and parameters used for the pelagic module are provided in Tab.A1-A3.
Importantly, the biogeochemical model resolves photoacclimation of phytoplankton, described by dynamical partitioning of resources to light harvesting pigments (Eq.A1), enzymes involved in carboxylation reactions (Eq.A1) and nutrient uptake sites (i.e., f LH + f C + f V = 1) as in Wirtz and Pahlow (2010).Uptake of each nutrient is optimally regulated (as expressed by a i in Eq.A2-A2), and following Pahlow (2005); Smith et al. (2009), optimality along the affinity-intracellular transport trade-off A3 for the definition of parameters).As a second novelty, the growth model uniquely describes the interdependence between limiting nutrients to be variable between full inter-dependence (as in product rule) and no-interdependence (as in Liebig's law of minimum) as a function of nitrogen quota (See Eq.A2).
For a detailed explanation of the phytoplankton growth model and solution of differential expressions in Eq.A1, A1 and A2 refer to Wirtz and Kerimoglu (2016).For enabling the spatial transport of the 'property variables' of phytoplankton such as Q i , f LH and f LH , they have been transformed to bulk variables by multiplying with the phytoplankton carbon biomass, i.e., B C .Parameterization of the phytoplankton model, except θ C , fall within the range of values used for the species considered by Wirtz and Kerimoglu (2016).The exact values of the parameters were established by manual tuning, given that important phytoplankton species such as various diatom and dinoflagellate species, and Phaeocystis sp. that dominate the phytoplankton composition in the SNS (eg., Wiltshire et al., 2010) have not been studied formerly within the presented model framework.
Phytoplankton losses are due to aggregation and zooplankton grazing (see below).Specific aggregation loss rate (Eq.A2) is described as a function of DOC that mimics transparent exopolymer particles (Schartau et al., 2007) to account for particle stickiness, multiplied by the sum of phytoplankton biomass and of POM reflecting density dependent interaction, which is equivalent to a quadratic loss term.Zooplankton dynamics are described only in terms of their carbon content, assuming stoichiometric homeostasis (Sterner and Elser, 2002).Grazing is described by a Holling Type-3 function of prey concentration (Eq.A2).A lumped loss term accounts for the respiratory losses and exudation of N and P in dissolved inorganic form (Eq.A2), Biogeosciences Discuss., doi:10.5194/bg-2017Discuss., doi:10.5194/bg- -104, 2017 Manuscript under review for journal Biogeosciences Discussion started: 30 March 2017 c Author(s) 2017.CC-BY 3.0 License.which are adjusted depending on the balance between the stoichiometry of zooplankton and that of the ingested food for maintaining the homeostasis (Eq.A2).Effect of the organisms at higher trophic levels, mainly by fish and gelatinous zooplankton are mimicked by a density-dependent mortality of zooplankton, modified by a function of total attenuation of Photosynthetically Available Radiation (PAR) (Eq.A2) to account for higher predation pressure exerted by fish at the off-shore regions of the North Sea, which amounts to about two times that in the coastal regions according to the estimates based on trawl surveys (Maar et al., 2014).
Settling velocity of POM, expressed by w P OM is prescribed as a constant value, whereas that of phytoplankton, w B is assumed to be modified by their nutrient (quota) status.As decreased internal nutrient quotas likely affect the cells ability to regulate buoyancy and lead to faster migration towards deeper, potentially nutrient rich waters (Boyd and Gradmann, 2002), we assume that maximum sinking rates realized at fully depleted quotas approach to a small background value with increasing quotas as has been observed especially for, but not limited to, diatoms (Smayda and Boleyn, 1965;Bienfang and Harrison, 1984).
Finally, all kinetic rates were modified for ambient water temperature, T (K) using the Q10 rule parameterized specifically for autotrophs and small heterotrophs (=bacteria for hydrolysis and remineralization) and for zooplankton.
Table A1.Source-sink terms of the dynamic variables of the pelagic module.The index i represents the elements C, N, P. By definition, QC = Q Z C =1, Qi=Bi/BC.Dynamics of Dissolved Inorganic Carbon (DIC) is not resolved, thus (Eq.A1) not integrated for i=C.Description of processes or functional relationships (capital letters) and of parameters (small letters) are provided in Tables A2 and A3, respectively.

A2 Benthic Module
The benthic module provides simplistic descriptions of the degradation of N and P from POM to DIM, their fluxes across the benthic-pelagic interface, removal of N due to denitrification and accounts for the sorption dynamics of P.
POM degrades into DIM in one step, described as a first order reaction, the rate of which is modified for temperature using the Q10 rule.POM flux into the sediments by settling of material from the water fuels the benthic POM (bPOM) (Eq.A4).
On the other hand, diffusive flux of DIM is possibly bi-directional, depending on the concentration gradient between water

Figure 1 .
Figure 1.Bathymetry of the model domain and the location of rivers considered in this study.Gray lines display the model grid.

Figure 2 .
Figure 2. Structure of the biogeochemical model.Model components (rectangles) comprise; B: phytoplankton, Z: zooplankton, POM and DOM: particulate and dissolved organic matter, DIM(N,P): dissolved inorganic matter (nitrogen,phosphorus), Fe-P: P adsorbed in ironphosphorus complexes (See Section 2.2.1 and supplementary material A for further details).C, N, P in small circles refer to carbon, nitrogen and phosphorus bound to each component, respectively, whereas fLH and fC are the allocation coefficients for light harvesting and carboxylation (Section 2.2.1).Boxes in dashed lines indicate model forcing.
resolving the 3-D dynamics of temperature, salinity and currents and 2-D dynamics of sea surface elevation and flooding-Biogeosciences Discuss., doi:10.5194/bg-2017-104,2017 Manuscript under review for journal Biogeosciences Discussion started: 30 March 2017 c Author(s) 2017.CC-BY 3.0 License.
days) were filled by linear interpolation.Loadings from the river Eider were calculated first by merging the data measured at the stations on two upstream branches, Eider and Treene, then by filling the short gaps (<28 days) by linear interpolation, replacing the larger gaps with daily climatology, and extending for 2000-2003 by using the climatology as well.To describe DIN deposition at the water surface, sum of annual average atmospheric deposition rates of N O x and N H 3 provided by EMEP (European Monitoring and Evaluation Programme, http://www.emep.int)were used.At the open boundaries in the north and west of the model domain (Figure 1), all state variables belonging to the phytoplankton and zooplankton compartments are

Finally, comparison of
the spatial structure of the model estimates to that of the satellite (ESA-CCI) data was achieved also through Taylor Diagrams.For this purpose, temporal matching was obtained by averaging the data from both sources for the period 2008-2010 for particular seasons of the year, and lateral matching by performing a 2-dimensional linear interpolation of the satellite data to the model grid.For the comparisons against the ICES and ESA-CCI data, only the middle-99 percentile of model and measurements were considered (i.e., leaving out the first and last 0.05th percentiles).Biogeosciences Discuss., doi:10.5194/bg-2017-104,2017 Manuscript under review for journal Biogeosciences Discussion started: 30 March 2017 c Author(s) 2017.CC-BY 3.0 License. 3 Results and Discussion 3.1 Evaluation of Model Performance by in-situ Data

Figure 3 .
Figure 3. Salinity [PSU] measured by ferry-box (a) and estimated by the model (b) along the route shown in the inset

Figure 4 .
Figure 4. Comparison of modeled and measured (ICES) temperature (abbreviated T in panels a,d,b,e) and salinity (abbreviated S in a,d,c,f) at the surface (a-c) and bottom (d-f) layers for the period 2006-2010.2-D histograms show the number of occurrence of simulation-measurement pairs.

Figure 5 .
Figure 5. Observations and model estimates of surface chlorophyll, DIN and DIP concentrations at the stations located along the coasts of the German Bight, operated by Alfred Wegener Institute (Helgoland and Sylt), Landesamt für Landwirtschaft, Umwelt und ländliche Räume des Landes Schleswig-Holstein (S.Amrum, Norderelbe) and Niedersächsichser Landesbetrieb für Wasserwirtschaft, Küsten-und Naturschutz (Norderney).Pearson correlation coefficients, and corresponding number of data points are shown at the top-right corner.

Figure 6 .Figure 7 .
Figure 6.As in Fig. 5, but for the stations located along the coasts of the Netherlands, operated by Rijkswaterstaat.

Figure 8 .
Figure 8.Comparison of simulated and measured (ICES) DIN (a,d,b,e) and DIP (a,d,c,f) at the surface (a-c) and bottom (d-f) layers for the period 2000-2010.2-D histograms show the number of occurrence of simulation-measurement pairs.

Figure 9 .
Figure 9. σT (a) and fluorescence (c) measured by Scanfish recorded during 13-19 July 2010 compared to σT (b) and chlorophyll concentration (d) estimated by the model averaged throughout the same period.Black line indicates the sea floor.Cruise track shown in (e).

Figure 10 .
Figure 10.Average temperature difference [K] between the surface and bottom layers, averaged throughout 2000-2010 for each month

Figure 12 .
Figure 12.DIN:DIP ratio in water binned over 1 m bottom depth intervals according to ICES data (a,d) and matching MAECS results (b,e) at the surface (a-c) and bottom (d-f) layers, within the eastern portion of the model domain (> 5 • E, exact locations shown in c and f) 13 for chlorophyll concentrations) and SPM concentrations.Higher Chl:C ratios during the non-growing (months 10-12 and 1-3) season similarly reflects light limitation due to low amounts of incoming short Biogeosciences Discuss., doi:10.5194/bg-2017-104,2017 Manuscript under review for journal Biogeosciences Discussion started: 30 March 2017 c Author(s) 2017.CC-BY 3.0 License.

Figure 13 .
Figure 13.Comparison of satellite (ESA-CCI, a,b) and MAECS (c,d) estimates of surface chlorophyll concentrations averaged over 2008-2010 and for different seasonal intervals of the year.2-D histograms (e,f) show the number of occurrence of simulation-satellite data pairs.

Figure 14 .
Figure 14.Chlorophyll:C ratio in phytoplankton, averaged over non-growing (a) and two halves of the growing season (b,c) of 2010.
are omitted in this study, and remain to be future research goals, along with Biogeosciences Discuss., doi:10.5194/bg-2017-104,2017 Manuscript under review for journal Biogeosciences Discussion started: 30 March 2017 c Author(s) 2017.CC-BY 3.0 License.other model refinements like improving the descriptions of the benthic processes and benthic-pelagic exchange, light climate as a function of SPM dynamics and composition of riverine nutrient fluxes.

Figure A1 .sa
Figure A1.Fraction of sorbed fraction of benthic phosphorus as functions of DO (a) and T (b), regulation of benthic denitrification rate as functions of DO (c) and T (d), and DO as a function of T (b,d).