The acclimative biogeochemical model of the southern North Sea

Ecosystem models often rely on heuristic descriptions of autotrophic 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 application 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 based on a set of novel concepts 5 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 within the German Bight characterized 10 by steep salinity, nutrient and chlorophyll gradients, as inferred from comparisons against observation data from long-term monitoring stations, sparse in-situ measurements, continuous transects, and satellites. The model displays skill also in capturing the formation of thin chlorophyll layers at the pycnocline, frequently observed within the stratified regions during summer. A sensitivity analysis reveals that the vertical distributions of phytoplankton concentrations estimated by the model can be qualitatively sensitive to the description of the light climate and dependence of sinking rates on the internal nutrient reserves. A 15 non-acclimative (fixed-physiology) version of the model predicted entirely different vertical profiles, suggesting that accounting for physiological flexibility might be relevant for a consistent representation of the vertical distribution of phytoplankton biomass. Our results point to significant variability in cellular chlorophyll to carbon ratio (Chl:C) across seasons and the coastal to offshore transition. Up to 3 fold higher Chl:C at the coastal areas in comparison to those at the offshore areas contribute to the steepness of the chlorophyll gradient. The model predicts also much higher phytoplankton concentrations at the coastal 20 areas in comparison to its non-acclimative equivalent. Hence, findings of this study provide evidence for the relevance of the physiological flexibility, here reflected by spatial and seasonal variations in Chl:C, for a realistic description of biogeochemical fluxes, particularly in the environments displaying strong resource gradients. 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. 25


Introduction
Modeling the biogeochemistry of coastal and shelf systems requires the representation of a multitude of interacting processes, not only within the water but also in the adjacent earth system components such as the atmosphere (e.g., nitrogen (N) deposition), land (e.g., riverine inputs), 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 the circulation patterns and stratification dynamics: the former is needed to describe the spread of nutrient-rich river plumes and exchange at the open ocean boundaries, and the latter for being able to capture the vertical gradients in the light and Published by Copernicus Publications on behalf of the European Geosciences Union.
nutrient conditions for primary productivity.The representation of biological processes and the two-way interactions between biological, chemical and benthic compartments in models is particularly challenging, given the complexity of physiological processes displayed by individual organisms, e.g., regarding the regulation of their internal stoichiometries (e.g., see Bonachela et al., 2016) and the differences in functional traits of species constituting communities (e.g., see Litchman et al., 2010).
Three-dimensional ecosystem models often describe the processes relevant to primary production, e.g., the nutrient and light limitation of phytoplankton (B), using heuristic formulations that have been shown to be inadequate in reproducing patterns obtained in laboratory experiments.For instance, light limitation is determined not only by the instantaneously available irradiance, but also by the amount of light-harvesting apparatus, i.e., chlorophyll (Chl) pigments maintained by the phytoplankton cells, which can change considerably through a process referred to as photoacclimation.However, photoacclimation is often completely ignored in 3-D model applications, or its effects are mimicked heuristically, for instance, by describing the chlorophyll-tocarbon ratio as a function of irradiance (Blackford et al., 2004;Fennel et al., 2006), which cannot capture the dependence of chlorophyll synthesis on nutrient availability (e.g., Pahlow and Oschlies, 2009;Smith et al., 2011;Wirtz and Kerimoglu, 2016).Similarly, the interaction of limitation by different nutrient elements is described by heuristic formulations, dichotomously either by a product rule or a threshold function, which, again, cannot reproduce complex patterns observed in laboratory conditions, such as the asymmetric cellular N : C and P : C ratios emerging under N-and P-limited conditions (Bonachela et al., 2016;Wirtz and Kerimoglu, 2016).Such simplifications in the description of primary production processes, in turn, potentially lead to flawed representations of nutrient cycling.Despite the recently revived theoretical work on stoichiometric regulation and photoacclimation (e.g., Klausmeier et al., 2004;Pahlow and Oschlies, 2009;Wirtz and Pahlow, 2010;Bonachela et al., 2013;Daines et al., 2014), an implementation of a model with a mechanistic description of the regulation of phytoplankton composition at a full ecosystem scale in a coupled physical-biological modeling framework remains lacking.In this study, we therefore present a 3-D application of the Model for Adaptive Ecosystems for Coastal Seas (hereafter MAECS), to the southern North Sea (SNS), for a decadal hindcast simulation.MAECS features a photoacclimative autotrophic growth model that has been recently introduced by Wirtz and Kerimoglu (2016), which resolves the regulation of the stoichiometry and composition of autotrophs employing an innovative suit of adaptive and optimality based approaches.
The SNS is part of a shallow shelf system (Fig. 1).The southeastern portion of the SNS, known as the German Bight surrounded by the intertidal Wadden Sea, is espe-cially characterized by steep gradients with respect to both nutrients (Hydes et al., 1999;Ebenhöh, 2004) and turbidity.The latter is largely determined by suspended particulate matter (SPM) concentrations (Tian et al., 2009;Su et al., 2015).These gradients are driven by a complex interplay of riverine and atmospheric fluxes, complex topography, residual tidal currents, density gradients, biological processing of organic matter (OM), benthic-pelagic coupling and sedimentation-resuspension dynamics (Postma, 1961;Puls et al., 1997;van Beusekom and de Jonge, 2002;Burchard et al., 2008;Hofmeister et al., 2016;Maerz et al., 2016).A number of modeling studies previously addressed the biogeochemistry of the North Sea, including the German Bight.In a majority of these studies, such as ECOHAM-HAMSOM (Pätsch and Kühn, 2008), NORWECOM (Skogen and Mathisen, 2009), ECOSMO-HAMSOM (Daewel and Schrum, 2013), HAMOCC-MPIOM (Gröger et al., 2013), ERSEM-NEMO (Edwards et al., 2012;Ford et al., 2017), ERSEM-POLCOMS (de Mora et al., 2013;Ciavatta et al., 2016) and ERSEM-BFM-GETM (van Leeuwen et al., 2015;van der Molen et al., 2016), large domains and relatively coarse grids were employed (≥ 7 km).While showing good skill in reproducing offshore 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), however, 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 offshore 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, the implementation of a highly complex phytoplankton growth model at an ecosystem scale, coupled to a hydrodynamic model and other biogeochemical compartments, and gaining some first insight into the relevance of acclimation to the modeling of coastal biogeochemistry; and 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.acclimation capacity of phytoplankton, and particularly by the high chlorophyll-to-carbon ratios at the coastal regions.

Observations
Observation data from Helgoland Roads, Sylt and 17 other monitoring stations reflect surface measurements.Extensive analyses of the data from Helgoland Roads have been previously performed by Wiltshire et al. (2008) and from Sylt by Loebl et al. (2007).Sparse measurements of temperature, salinity, dissolved inorganic nitrogen (DIN), dissolved inorganic phosphorus (DIP) and chlorophyll were obtained from the online database of the International Council for the Exploration of the Sea (ICES, www.ices.dk).Continuous Scanfish and FerryBox measurements were performed within the operation of the Coastal Observing System for Northern and Arctic Seas (COSYNA, Baschek et al., 2016).Data collection, processing and quality control of the Scanfish data are described by Maerz et al. (2016) and of the FerryBox data by Petersen (2014).The satellite dataset used here is the Ocean Colour Climate Change Initiative (CCI), version 3.1, European Space Agency (ESA), available online at http://www.esa-oceancolour-cci.org/.Chlorophyll estimates of the satellite product were bias-corrected according to the product user guide (Grant et al., 2017): C bc = 10 log 10 (C)+δ , where C bc , C and δ are, respectively, the bias-corrected, raw and log 10 bias estimates for chlorophyll concentrations.B is phytoplankton, Z is zooplankton, POM and DOM are particulate and dissolved organic matter, DIM (N, P) is dissolved inorganic matter (nitrogen and phosphorus), and bAP is the P adsorbed in iron-phosphorus complexes (see Sect. 2.2.1 and Appendix A for further details).C, N and P in small circles refer to carbon, nitrogen and phosphorus bound to each component, respectively, whereas f LH and f C are the allocation coefficients for light harvesting and carboxylation (Sect.2.2.1).Boxes in dashed lines indicate model forcing.

Model
The major processes taken into account by the model are the lower trophic food web dynamics, phytoplankton ecophysiology and basic biogeochemical transformations in the water, and the transformation of N and P species in the benthos (Fig. 2 and Sect.2.2.1).Physical processes are resolved by the coupled 3-D hydrodynamical model, GETM (General Estuarine Transport Model;Sect. 2.2.2).Turbidity caused by suspended particulate matter, nutrient loading by rivers and atmospheric nitrogen deposition were considered as model forcing (Sect.2.2.3).The model grid and rivers 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 and phosphorus, and, importantly, acclimation processes involved in phytoplankton growth.In MAECS, the acclimation of phytoplankton is resolved by a scheme recently introduced by Wirtz and Kerimoglu (2016), which describes the instantaneous or transient optimization of physiological traits, x, by the extended optimality principle: www.biogeosciences.net/14/4499/2017/Biogeosciences, 14, 4499-4531, 2017 where δ x corresponds to the flexibility of traits (Eq.A18), i expands to N and P, 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.Specifically, threelevels of acclimative regulations are considered (see Fig. 2 in Wirtz and Kerimoglu, 2016): 1. Machinery allocation: we describe the changes in allocations to light-harvesting, carbon-fixation and nutrient acquisition machineries, as also in Wirtz and Pahlow (2010).These allocations correspond to the synthesis of cellular structures such as chloroplasts for absorbing light, Rubisco enzyme involved in carboxylation process and proteins for gathering nutrient molecules; therefore, we track these fractional allocations with two dynamic state variables, f LH and f C , that describe the allocations for light harvesting and carboxylation, while the allocation for nutrient uptake, f V , is assumed to be the rest, 1 Here, the flexibility term, set to , regulates the speed of optimization as determined by the differential terms in Eq. (1).
2. Nutrient affinity-processing optimality: we assume that there is a tradeoff between nutrient affinity and processing, and the optimal affinity fractions for each nutrient, f A i , are instantaneously optimized, such that d x /d t = 0 and f A i are algebraically found by setting ∂V C ∂x = 0 (Pahlow, 2005;Smith et al., 2009).
3. Nutrient uptake activity: (down-) regulation of the uptake rate of nutrients, which is often formulated as a linear function of nutrient quotas (Morel, 1987) in traditional models, is in our approach described by the instantaneously optimized uptake activity trait, a i .Assuming that energy expenditure for taking up each nutrient depends on the metabolic needs, values of a i are found by scaling their marginal growth benefits (Eq.A17).
Driven by the variations of these physiological traits, Chl : C : N : P stoichiometry varies continuously depending on ambient light and nutrient conditions and on the metabolic demands of autotrophic cells.As a further novel aspect of the acclimation model, multiple limitation is described as a queuing function, which allows formulating the co-limitation strength as a function of internal nitrogen reserve, q N , instead of prescribing it to be either high as by a product rule or low as by a threshold (Liebig) function (Wirtz and Kerimoglu, 2016).A detailed description of the phytoplankton growth module can be found in Wirtz and Kerimoglu (2016).
Equations and parameters of the model are provided in Appendix A1.
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 matter (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 matter (DOM), before becoming again DIM and closing the cycle (Appendix A1).As a relevant aspect of the model, while the sedimentation speed of POM (w POM ) is prescribed as a constant value, that of phytoplankton, w B , is assumed to be modified by its 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 converge 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).Although the phytoplankton sinking is often parameterized as a constant rate in 3-D modeling applications, similar formulations of increasing sinking rates under nutrient stress have also been used (e.g., Vichi et al., 2007).
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 and also as a function of temperature based on the correlation observed in the field (Jensen et al., 1995).Further details are provided in Appendix A2.

Hydrodynamic model and model coupling
The General Estuarine Transport Model 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) and Stips et al. (2004).GETM utilizes the turbulence library of the General Ocean Turbulence Model (GOTM) to resolve vertical mixing of density and momentum profiles with a kε two equation model (Burchard et al., 2006).GETM was run in baroclinic mode, resolving the 3-D dynamics of temperature, salinity and currents and 2-D dynamics of sea surface elevation and flooding-drying of cells at the Wadden Sea.Following Gräwe et al. (2016), we assumed the bottom roughness length to be constant throughout the domain, and z 0 = 10 −3 m.We used 20 terrain-following layers and a curvilinear grid of 144 × 98 horizontal cells, providing a hor-izontal resolution of approximately 1.5 km at the southeast corner and 4.5 km at the northwest 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 s 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 (Lemmen et al., 2017).Meteorological forcing originated from an hourly resolution hindcast by the limited area model COSMO-CLM (Geyer, 2014).Boundary conditions for surface elevations 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).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 bottommost pelagic box of the corresponding water column in terms of a unidirectional flux of POM from the pelagic to the benthic states and a bidirectional 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 s, as for the 3-D fields in GETM.The 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 (PAR) 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's type I 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, such as 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 the 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 SPM = 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 11 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) presented 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, and 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), the 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 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 from 2000 to 2003 by using the climatology as well.To describe DIN deposition at the water surface, the sum of annual average atmospheric deposition rates of oxidized and reduced nitrogen 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 (Fig. 1), all state variables belonging to the phytoplankton and zooplankton compartments are assumed to be at zero gradient.For DIM, DOM and POM, monthly values of ECO-HAM (Große et al., 2016), interpolated to 5 m depth intervals, are used as clamped boundary conditions.

Quantification of model performance
For the comparisons with the data at monitoring stations, sparse in situ measurements from the ICES database, and with the satellite dataset; Pearson correlation coefficients,  (Jolliff et al., 2009).For the comparisons against the ICES and satellite data, only the middle 99 % of simulated and measured values were considered (i.e., leaving out the first and last 0.5 %).
For the ICES data, temporal matching was identified at daily resolution, vertical matching was obtained by comparing the measurements within the upper 5 m from the sea surface and within the 5 m above the sea floor with the model estimates at the topmost and bottommost layers, and finally horizontal matching was obtained by calculating the average of the values from four nearest cells surrounding the measurement location, inversely weighted by their Cartesian distance.For the satellite data, the temporal matching was obtained by averaging the data from both sources for the period 2008-2010 for particular seasons of the year and horizontal matching by performing a two-dimensional linear interpolation of the satellite data to the model grid.The extraction of the hourly model temporally matching to the Scanfish data was based on the hourly binned average time for each cast (defined as a full downward and upward undulation cycle), and 3-D spatial matching was obtained by constructing an average vertical profile from the four closest cells to the average coordinate of each cast.For facilitating the qualitative comparison of the simulated chlorophyll and the Scanfish measurements of fluorescence, which have different units and signal strengths, normalized anomalies were used, according to pi = (p i − p )/σ p , where pi and p i are the normalized anomaly and raw value of a given data point, and p and σ p are the mean and standard deviation of all data points.

Evaluation of model performance by in situ data
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) demonstrates that the model captures the horizontal salinity distribution (Fig. 3).In particular, the contrast between the northwestern model domain characterized by the rapid flushing of the coastal freshwater input and the southeastern model domain (i.e., German Bight) characterized by a strong and permanent salinity gradient is well captured.Confinement of the salinity front during winter towards the coast and its seaward intrusion, especially during early spring, and the smaller-scale modulations that appear to be controlled by the spring-neap cycle are both reproduced by the model.
Comparison of simulated surface and bottom temperatures with those extracted from the ICES dataset for the pe-riod 2006-2010 are provided in Fig. 4. The high correlation scores and low bias attained for water temperature and salinity suggest that the model can reproduce the seasonal warming, spread of freshwater discharges and thermohaline stratification dynamics.However, in a relatively small number of instances, surface temperatures are underestimated and bottom temperatures are overestimated, which indicates that not all stratification events were captured.Almost all of these instances are found to be located either at the northeastern margin (> 4 • E and > 55 • N) or at the northwestern corner (< 4 • E and > 54 • N) of the model domain, i.e., close to the open ocean boundary (Fig. 1).
Comparisons of surface chlorophyll, DIN and DIP concentrations estimated by the model with the measurements in 19 stations scattered across the southern North Sea are shown in Figs.5-7, and the corresponding skill scores are listed in Table 1.Estimates of average nutrient concentrations and the timing of their depletion and regeneration in a majority of stations agree well with the observations, as indicated by the frequency of high and moderately high scores (Table 1).Notably, at several stations (e.g., Sylt, T8, T36, T26, T22, T11 and T12) the difference between the relative bias for DIP and DIN (i.e., B * DIP −B * DIN ) was relatively large (with 55 % being the highest at T22), suggesting a tendency for underestimating the DIN : DIP ratio, although this was not the case for the comparison against ICES measurements (see below).Relative to the nutrients, the performance of the model in estimating chlorophyll is lower, especially at the stations located along the Dutch coast (Fig. 6, Table 1).However, for about half of the 10 stations where data are available, the model performance is at moderate levels.
The comparison of model results with the DIN, DIP and chlorophyll measurements available at the ICES database at the surface and bottom layers for the entire simulation period indicates a negligible normalized mean bias (≤ 12 %) and correlation coefficients at around 0.6-0.7 for nutrients and about 60 % overestimation and correlation coefficients of about 0.3 for chlorophyll (Fig. 8, Table 1).The modeled variability for all three biogeochemical state variables is within an approximately 50 % envelope of the observed variability (Fig. 8).
For an assessment of the accuracy of the simulated vertical distributions, water density (expressed as σ T ) and fluorescence captured by a Scanfish cruise (Heincke Cruise HE331) during 13-19 July 2010 were compared to those estimated by the model (chlorophyll for fluorescence) averaged over the same time period (Fig. 9).This period was characterized by significant thermal summer stratification reaching deep into the near coastal regions of the German Bight.Thus, σ T reflects two major mechanisms that control the distribution of phytoplankton: the first is the characterization of the vertical gradients by denser water at the bottom layers, which is mainly driven by thermal stratification as suggested by temperature profiles (not shown).The second is the characterization of the horizontal gradients by lighter water at the coasts, driven by low salinity due to the freshwater flux from the rivers.The model can accurately reproduce both vertical and horizontal density gradients, although some discrepancies exist, such as the slightly underestimated depth of the pycnocline and steepness of lateral gradients at around the coastal section.Fluorescence measurements along the Scanfish track in July 2010 indicate frequent occurrences of subsurface chlorophyll maxima (Fig. 9).These are in some cases in the form of higher concentrations below the pycnocline but in some others appear as thin layers at around the pycnocline.While the deep chlorophyll maxima are prevalently found in stratified offshore regions, the well-mixed shallower regions mostly show homogeneously distributed high chlorophyll concentrations throughout the water column due to higher dissipation rates (Maerz et al., 2016).The MAECS simulation agrees qualitatively well with these patterns and captures the spatial variability of the observed vertical chlorophyll distribution (Fig. 9).

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 the Scanfish transect (Fig. 9) showed that the physical model has the potential to realistically capture the density stratification.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)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010), the areal extent and seasonality of stratification within the SNS is shown 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 a maximum areal coverage and intensity (slightly above 8 K) in July.
Simulated climatological concentrations of DIN and DIP display steep coastal gradients along the coasts of the German Bight (Fig. 11), both during the non-growing season (months 1-3 and 10-12) and the growing season (months 4-9).Within the ROFI (region of freshwater influence, Simpson et al., 1993) of the Rhine, nitrogen concentrations decrease about 5 fold (from ≥ 48 mmolN m −3 to 8-16 mmol N −3 ) within a few grid cells, corresponding to about 10-15 km distance.In the German Bight, the non-growing season is similarly characterized by a thin stripe of high nutrient concentrations along the coast, whereas during the growing season especially phosphorus becomes depleted outside a confined zone of the Elbe plume.At the offshore areas, nutrient concentrations during the growing season are considerably lower than those during the non-growing season, driven by the phytoplankton growth both directly by nutrient uptake and, for the case of nitrogen, also indirectly by fueling denitrification in the sediment.The DIN : DIP ratio in the offshore re-  gions is close to the Redfield molar ratio of 16 : 1 throughout the year, reflecting oceanic conditions, while much higher at the coastal areas, particularly during the non-growing season, reflecting the high N : P content of the continental rivers (Radach and Pätsch, 2007).This transition from high coastal to low offshore N : P ratios is qualitatively consistent with observations (e.g., Burson et al., 2016).
Both the satellite (ESA-CCI) images and our model estimates, averaged again for the non-growing and growing season, suggest steep coastal gradients in chlorophyll concentrations (Fig. 12) similar to the nutrient gradients shown above (Fig. 11).The large-scale agreement in coastal gradients results in high correlation coefficients (Fig. 12, Table 1).
The normalized mean bias is small for the non-growing season but relatively high and positive (i.e., overestimation) for the growing season.Higher model estimates at the lower range (0-10 mgChl m −3 ) are responsible for this positive bias, which is particularly the case during the first half of the growing season, where the bias is highest (Table 1).
Our simulation results indicate significant spatiotemporal variability in the Chl : C ratio, even when the seasonal averages are considered, i.e., omitting short-term variability (Fig. 13).The Chl : C ratio is generally higher at the coasts than offshore.Higher Chl : C ratios during the non-growing (months 10-12 and 1-3) season similarly reflect light limitation due to low amounts of incoming shortwave radiation at  the water surface.The simulated spatiotemporal differences in Chl : C ratios reach to about 3 fold between different seasons of the year and between offshore and coastal areas.The latter suggests that the differential acclimative state of phy-toplankton cells amplifies the steepness of the chlorophyll gradients across the coastal transition shown in Fig. 12.
For gaining a better understanding of the relevance of acclimation in capturing the coastal gradients, we considered a simplified, non-acclimative version of the model in which  after we will refer to them as the "fixed" model in short, without specifying the particular parameterization.
The annual average coastal phytoplankton concentrations estimated by the acclimative model are much higher than those estimated by the fixed model, with no significant difference between the surface and water-column-averaged values (Fig. 14a, b).In the offshore areas, the estimates of the acclimative model are higher than those of the fixed model at the surface (Fig. 14a), but slightly lower when water column averages are considered (Fig. 14b), indicating that the phytoplankton growth occurs mostly at the bottom layers in the fixed-trait model, which is consistent with the daily vertical profiles in the fixed-trait model (Fig. B4c).Importantly, these results suggest that a coastal gradient in phytoplankton concentrations is predicted by a non-acclimative model, which is presumably driven by the nutrient gradients (Fig. 11), but a much stronger gradient emerges when the acclimation processes are resolved.Specific to this example, towards the coast, phytoplankton adapt to the deteriorating light climate (Fig. B1) and increasing nutrient availability by investing more in the light-harvesting machinery, as indicated by the increasing Chl : C ratios (Fig. 13), and thereby achieving higher coastal production rates than in the case where their physiology is fixed.As a result of increasing Chl : C ratios towards the coast, the chlorophyll concentrations display an even stronger gradient than that of the biomass: at the surface layer, the increase of biomass concentrations towards the coast is about 3.5 fold (from about 10 to 35 mmolC m −3 ), while that of chlorophyll is about 7 fold (from about 2 to 14 mg m −3 ) along the transect shown in Fig. 14.

Discussion
In order to assess the performance of the new model system presented for the first time in this study, we employed several independent observation sources and types: Ferry-Box measurements to assess the horizontal distribution of salinity (Fig. 3); sparse in situ measurements from the ICES dataset for an overall evaluation of the physical and biogeochemical model (Figs. 4,8); measurements from 19 mon-  itoring stations for evaluating the estimates for DIN, DIP and chlorophyll at specific locations (Figs.5-7); Scanfish measurements for evaluating the vertical density and chlorophyll profiles (Fig. 9); and finally the satellite observations for evaluating the model skill regarding the horizontal distribution of climatological chlorophyll concentrations (Fig. 12) and attenuation of light (Fig. B1).
The physical model can provide a realistic description of the hydrodynamical processes foremost relevant for modeling the biogeochemistry of the system.Horizontal circulation patterns are captured as evidenced by the salinity distribution being in agreement with the observations (Figs. 3, 4).The density structure of the system during summer, driven by a complex interplay between the salinity gradients, heat fluxes at the surface and tidal stirring, is realistically captured, although the pycnocline depth seems to be underestimated (Fig. 9).Accordingly, temperature estimations match well with the observations, although there are cases where the stratification events are not reproduced by the model (Fig. 4), most of which are found to be within the western portion of the model domain.The areal extent and seasonality of stratification (Fig. 10) are in agreement with those reported by ear-lier studies (Schrum et al., 2003;van Leeuwen et al., 2015).For nutrient concentrations, a relative bias of ≤ 12 % and correlation coefficients between 0.58 and 0.72 correspond to a high and moderately high model skill, respectively (Table 1).For the pointwise comparisons of chlorophyll, model skill was moderate for the sparse measurements included in the ICES database, and for the stations in the German Bight (Table 1), but mostly low for the stations within the western portion of the model domain.The comparison of climatological averages of the simulated chlorophyll with those of the satellite observations resulted in high correlations for all seasons and a low to moderately low bias, except during the early growing season (Table 1).
The model captures the subsurface chlorophyll maxima occurring in the deeper parts of the model domain (Fig. 9).This phenomenon has been previously documented in the southern North Sea (Weston et al., 2005;Fernand et al., 2013).Former 3-D modeling studies, such as that of van Leeuwen et al. (2013), apart from capturing the presence of a deep chlorophyll maximum, did not reproduce the rich variability revealed by the observations.Our comparative analysis shows that the formation and maintenance of such structures are critically dependent on the parameterization of the sinking rate of phytoplankton (Fig. B4a) and underwater light climate (Fig. B4b).The sinking speed of phytoplankton in 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 Appendix A1) but also indirectly emulates chemotactic migration as is typical for dinoflagellates (Durham and Stocker, 2012).This quota dependency results in considerable spatial and seasonal variability in sinking rates (Fig. B3).The critical dependence of the formation and maintenance of verti-cal chlorophyll structures on the functional representation of sinking underlines the relevance of a consistent description of the intracellular regulation of nutrient storages.The latter, in turn, is determined by the metabolic needs, such as the intensity of light limitation, and hence, investments in the synthesis of pigmentary material (Wirtz and Kerimoglu, 2016).Indeed, the non-acclimative (fixed-trait) version of the model (Appendix B3) predicts qualitatively different vertical profiles of phytoplankton biomass (Fig. B4c), although the sinking parameterization in that simplified version is identical to that in the fully acclimative version.The non-acclimative www.biogeosciences.net/14/4499/2017/Biogeosciences, 14, 4499-4531, 2017 model version might be tuned to match the observed vertical distributions of phytoplankton; however, this would probably be at the expense of compromised performance in some other respects, such as the horizontal gradients, or timing and amplitude of chlorophyll blooms.
We conclude from the extensive model performance assessment that the model reproduces the main physical and biogeochemical characteristics of the southern North Sea, especially within the German Bight, where the model resolution is finest (Fig. 1), and the influence of fluxes at the open boundaries is relatively small, given the predominantly counterclockwise circulation pattern (Becker et al., 1992).The process of performance assessment also helped in identifying the possibilities for further model refinement.For instance, a comparison with the satellite observations revealed that the light attenuation in the offshore areas is overestimated by the model, primarily because of the contributions by the climatological SPM forcing (Fig. B1).A likely con-sequence of the overestimated attenuation is an underestimation of the depth of primary production (e.g., Fig. B4b), and this may, in turn, explain the overestimated chlorophyll concentrations in the offshore areas during the growing season (Fig. 12b, d).Another source of error regarding the SPMcaused turbidity is the fact that, at specific coastal sites, such as at the Noordwijk-10 station, the measured SPM concentrations show considerable interannual variations that can obviously not be represented by the climatological SPM forcing (Fig. B2), which may explain the particularly low correlation coefficient (0.14) obtained at this station for chlorophyll.A better representation of the SPM-caused turbidity might be achieved by an explicit description of the SPM dynamics (e.g., as in van der Molen et al., 2016).Coupling the biogeochemical model with such an SPM model would then also allow the description of the two-way interactions, i.e., not only light limitation (Tian et al., 2009), but also the acceleration of sinking of SPM by the production of transparent exopolymer particles (Schartau et al., 2007;Maerz et al., 2016).At the stations within the ROFIs of major rivers, such as the Norderelbe and S. Amrum (Fig. 5) and Noordwijk-2 and Noordwijk-10 (Fig. 6), the skill scores are relatively low (Table 1).These stations, especially Noordwijk-2 and Noordwijk-10, are located where the concentrations change dramatically within 10-15 km (e.g., Fig. 11).Accordingly, a slight error by the physical model in predicting the salinity front, e.g., because of an inadequate representation of the tidal dynamics, might result in considerable deviation of the estimated concentration of biogeochemical variables from the measurements.The relatively coarser model resolution around the Dutch coast might therefore explain the consistently lower skill scores obtained at the Dutch stations.Identifying such potential inadequacies of the physical model requires further investigation, such as an assessment of the tidal constituents at the tidal gauges (e.g., Gräwe et al., 2016).Another potential source of error for the mismatches within the ROFIs is the potential flaws in the description of riverine loadings, such as assuming that the nondissolved fractions of the total nitrogen and total phosphorus are entirely in labile form (Sect. 2.2.3).Although the earlier replenishment of phosphorus relative to nitrogen in the coastal sites is often reproduced (e.g., Sylt, Noordwijk-2, Noordwijk-10, Terschelling-4), some delays occur in stations such as Norderney, which probably reflects the oversimplification of the benthic processes with respect to the description of oxygen-driven iron-phosphorus complexation kinetics (Appendix A2), which have 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).
The model predicts steep coastal gradients in nutrient concentrations (Fig. 11), which is in line with prior observations (e.g., Brockmann et al., 1999;Hydes et al., 2004).The maintenance of these gradients during winter is explained by the limited horizontal mixing due to the density gradients caused by the freshwater influx from the land (Simpson et al.,  1993; Hydes et al., 2004) and the trapping of this nutrientrich freshwater at the coast due to the alongshore currents in the study system driven by predominantly westerly winds and the coriolis forcing (Becker et al., 1992;Simpson et al., 1993).During the warmer seasons when the offshore waters are stratified, owed to the presence of horizontal salinity gradients, a mechanism similar to the estuarine circulation (Simpson et al., 1990) was suggested to further promote these gradients along the Wadden Sea, as well as in regions far from river inputs (Burchard et al., 2008;Flöser et al., 2011;Hofmeister et al., 2016).Coastal waters remaining nutrient-replete during the growing season lead to high phytoplankton concentrations (Fig. 12), despite the higher turbidities at the coastal waters (Fig. B1).
The comparison of the present model with earlier attempts is neither in the scope of this study nor possible without a dedicated benchmarking effort, using standardized forcing data and skill performance assessment datasets and methodology (e.g., as in Friedrichs et al., 2007).Even a qualitative comparison is difficult, given that spatial and temporal binning of the data, frequently employed in model validation (such as in our Fig.12), can dramatically impact the skill scores and that pointwise comparisons with sparse observation datasets (such as in our Figs.4 and 8) are rarely performed (de Mora et al., 2013).However, the skill of the presented model in estimating the chlorophyll concentrations in the SNS can argued to be at least comparable to that of the recent modeling applications for a relevant region (e.g., Edwards et al., 2012;de Mora et al., 2013;Ciavatta et al., 2016;Ford et al., 2017, noting that all these studies had larger model domains and were evaluated for different time intervals).This is noteworthy, given that phytoplankton is represented by a single species in our model, whereas in other modeling approaches several species or groups are resolved.The inclusion of multiple functional types is motivated by the spatial and seasonal variability in the phytoplankton composition observed in the field: coastal areas of the SNS are dominated by diatoms throughout the year in some sites (e.g., Alvarez-Fernandez and Riegman, 2014) and during spring in some others, later replaced by Phaeocystis during summer (e.g., van Beusekom et al., 2009), whereas the offshore areas are often dominated by dinoflagellates especially during summer (Freund et al., 2012;Wollschläger et al., 2015).These phytoplankton groups differ from each other in a number of traits, including the physiological traits that determine their ability to access the (mineral and light) resources and build biomass.For instance, in an experimental work, two diatom species were shown to have on avwww.biogeosciences.net/14/4499/2017/Biogeosciences, 14, 4499-4531, 2017 erage more than 3-fold-higher Chl : C ratios than those of two dinoflagellate species (Chan, 1980), therefore making them more tolerant to the light-limited conditions of the turbid, coastal waters.In the presented approach, the cellular composition of the single, but acclimative, phytoplankton group dynamically approaches towards (for some traits, instantaneously adopts) the physiological state of the ideal resource competitor in a given environment, which, in nature, happens through various processes -from the plastic response of the individual cells to the species sorting at the community level.In a traditional, plankton functional type model, on the other hand, the species with the most suitable traits would become the most dominant among others, while the proximity of the physiological traits to the theoretical optima, and thus the overall productivity, would be determined by the resolution of physiological traits as represented by the defined clones.The worst-case scenario is when there is only one non-acclimative group, as illustrated in our experiment: at the turbid but nutrient-rich coastal areas, prioritization of the light harvesting over the nutrient acquisition machinery, as evidenced by the higher Chl : C ratios predicted by the acclimative model (Fig. 13), leads to better fitness and thus higher phytoplankton concentrations in comparison to their non-acclimative equivalents (Fig. 14).Moreover, because of the high Chl : C ratios at the coastal areas (Fig. 13), chlorophyll concentrations display even steeper gradients than the phytoplankton concentrations (Fig. 14).The transitional Chl : C pattern suggested by our model has been previously identified based on monitoring data by Alvarez-Fernandez and Riegman (2014).The Chl : C  ratios ranging between 0.01-0.1 gChl gC −1 at the coastal stations and 0.002-0.02gChl gC −1 at the offshore stations reported by Alvarez-Fernandez and Riegman (2014) envelop our estimated seasonal average values of 0.045 and 0.015 within the respective regions.According to the simulation results, Chl : C ratios also differ considerably between the nongrowing and growing season, with higher values during winter, due to low light availability.A similar seasonal amplitude in Chl : C has been found by Llewellyn et al. (2005) for the English Channel, with higher ratios during winter.
As mentioned above, the physiological composition is not the only relevant trait for determining the community composition in the study system.Diatoms are fast growers and are defended against the efficient microzooplankton grazers, but this comes at the cost of silicate requirement for their growth (Loebl et al., 2009) and higher sedimentation losses (Riegman et al., 1993).Phaeocystis spp.are slow growers, but, by forming large colonies, they are well defended against zooplankton (Peperzak et al., 1998).Finally, the dinoflagellates, also despite being slow growers, are mobile (Durham and Stocker, 2012) and mostly have access to alternative nutrient sources through their phagotrophic abilities (Löder et al., 2012).The representation of zooplankton with a single group may also be an oversimplification, as the microzooplankton and mesozooplankton have considerably different growth rates (Hansen et al., 1997) and functional responses to prey availabilities (Kiørboe, 2011).Moreover, effects of temperature on mesozooplankton occur 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 Q 10 rule we applied for grazing rates (Appendix A1).None of these ecophysiological aspects were taken into account in our model, and this may explain some of the discrepancies between the simulated and observed chlorophyll concentrations.In future work, inclusion of few other phytoplankton groups, each being acclimative, and one additional zooplankton group is foreseen.While the consideration of other phytoplankton traits should be straightforward, the inclusion of phagotrophy as an additional physiological allocation trait represented by a state variable is possible (e.g., as in Chakraborty et al., 2017) but would require the re-derivation of the model equations.
In this study, we described the implementation of a coupled physical-biogeochemical model in the southern North Sea 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, spatiotemporal 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 German Bight 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 offshore regions of the German Bight.
In 3-D model applications so far, photoacclimation of phytoplankton has been either ignored altogether or it has been accounted for in a heuristic sense, where the change in the 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).Our findings suggest that the steep chlorophyll gradients across the coastal transition zone are mainly driven by the nutrient gradients, but they are first amplified by the acclimative capacity, and then further by higher Chl : C ratios at the coastal waters.The large variations in simulated Chl : C ratios within the SNS, both in space and time, indicate that ignoring 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), based on the variability of Chl : C ratios at global scales.
Here we show that this warning applies especially in the coastal environments characterized by steep resource gradients, 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).
Data availability.Atmospheric forcing (COSMO-CLM), water level and currents (TRIM), temperature and salinity (HAMSOM) estimations and bathymetry of the North Sea used to construct model forcing can be accessed from www.coastdat.de.Ferrybox data are available from COSYNA data portal: http://codm.hzg.de/codm/.ICES dataset used for model validation is available at http: //ices.dk.Satellite data (CCI OC v3.1) used for model validation are available at http://www.esa-oceancolour-cci.org.For all other validation data used in this study, individual inquiries need to be made (see the Acknowledgements).Model data presented in this study can be provided by Onur Kerimoglu (kerimoglu.o@gmail.com)upon request.
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 Tables A1-A3.
Importantly, the biogeochemical model resolves the photoacclimation of phytoplankton, described by the dynamical partitioning of resources to light-harvesting pigments (Eq.A7), enzymes involved in carboxylation reactions (Eq.A8) 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, firstly in terms of the priority of each nutrient, among others (as expressed by a i in Eqs.A16-A17), and secondly following Pahlow (2005); Smith et al. (2009), along the affinity and intracellular transport rate 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 interdependence (as in product rule) and no interdependence (as in Liebig's law of the minimum) as a function of nitrogen quota (See Eq.A13).For a detailed explanation of the phytoplankton growth model and solution of differential expressions in Eqs.(A7), (A8) and (A17) 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 into bulk variables by multiplying with the phytoplankton carbon biomass, i.e., B C .Parameterization of the phytoplankton model, except θ C , falls within the range of parameter values used 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 spp. that dominate the phytoplankton composition in the SNS (e.g., Wiltshire et al., 2010), have not been studied previously within the presented model framework.
Phytoplankton losses are due to aggregation and zooplankton grazing (see below).The specific aggregation loss rate (Eq.A19) 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 densitydependent 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 Holling's type III function of prey concentration (Eq.A20).A lumped loss term accounts for the respiratory losses and exudation of N and P in dissolved inorganic form (Eq. 21), which are adjusted depending on the balance between the stoichiometry of zooplankton and that of the ingested food for maintaining homeostasis (Eq.A22).The effects of organisms at higher trophic levels, mainly of 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.A23) to account for higher predation pressure exerted by fish at the offshore 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).All kinetic rates were modified for ambient water temperature, T (K), using the Q 10 rule parameterized specifically for autotrophs and small heterotrophs (i.e., bacteria for hydrolysis and remineralization) and for zooplankton.

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, the removal of N due to denitrification and accounts for the sorption dynamics of P.
POM degrades into DIM in one step, described as a firstorder reaction (Eq.A30), the rate of which is modified for temperature using the Q 10 rule (Eq.A37).POM flux into the sediments by settling of material from the water fuels the benthic POM (bPOM; Eq.A31).The diffusive flux of DIM is possibly bidirectional, depending on the concentration gradient between water and soil (Eq.A32).Inorganic phosphorus (Eq.A28) is assumed to exist in two states: sorbed and dissolved state.The fraction of the sorbed state is given by a function of dissolved oxygen to account for the production and adsorption of Fe-P complexes in oxic conditions and their desorption in anoxic conditions (Eq.A33).Given the observed inverse relationship between temperature and oxygen concentrations in sediments (e.g., Jensen et al., 1995), DO is heuristically estimated as a function of T (Eq.A36) to capture the seasonal hypoxia events.The resulting functional relationships between the sorbed fraction of total inorganic phosphorus, T and DO are shown in Fig. A1a-b.Following the simplistic approach used for the ECOHAM model (Pätsch and Kühn, 2008), the denitrification rate is estimated from the degradation rate (Eq.A34) using empirically derived ratios and stoichiometric conversions, considering, in addition, the limitation imposed by the available DIN and the inhibition by DO (Soetaert et al., 1996).The resulting functional relationships between denitrification, T and DO are shown in Fig. A1c-d.
Table A1.Source-sink terms of the dynamic variables of the pelagic module.The index i represents the elements C, N and P. By definition, The dynamics of dissolved inorganic carbon (DIC) are not resolved; thus Eq.A4 is not integrated for i = C. Descriptions of processes or functional relationships (capital letters) and of parameters (lowercase letters) are provided in Tables A2 and A3, respectively.gies (as required by GETM), which was extracted for 5 m depth from an original 3-D climatological (daily) dataset constructed by a statistical regression approach (Heath et al., 2002) and used in ECOHAM (e.g., Große et al., 2016).Here we aim to gain insight into the realism of the light climate represented by the model, in particular with respect to the turbidity caused by SPM, as the rest of the turbidity is mainly caused by the simulated variables, in particular phytoplankton, the realism of which is discussed extensively in the main text (e.g., Figs. 8,9,12).Table A4.Source-sink terms of the dynamic variables and functional relationships of the benthic module.For POM i and DIM i , i = N, P. Description of parameters are provided in Table A5.

Dynamics
POM exchange with water DIM exchange with water Fraction of inorganic P in adsorbed phase bAP bDIP+bAP  Seitzinger and Giblin (1996).data and, in turn, lead to mismatches in the simulated phytoplankton, for instance, in the form of timing errors (Fig. 6).

B2 Phytoplankton sinking rates and vertical chlorophyll profiles
Average phytoplankton sinking rates estimated by the model at the surface layers vary considerably across seasons, with higher sinking rates at the offshore areas during summer than in the rest of the year (Fig. B3).Sinking rates at the bottom layer are lower than at the surface, both during spring and summer (Fig. B3).These patterns are driven by the nutrient quota dependence of sinking rates (Eq.A25) and low nutrient quotas at the surface layers in the stratified offshore areas during summer (Fig. 10), caused by the depletion of nutrients at the surface layers (Fig. 11).Because of the continuous supply of nutrients from the sediments, phytoplankton at the bottom layer maintain high nutrient quotas throughout the year and therefore have lower sinking rates.The range of observed mean sinking rates (−0.92-1.14m d −1 ) in the Rhine ROFI (Peperzak et al., 2003) is roughly consistent with that estimated by the model.The higher estimated sinking rates during summer than in spring, and at the surface rather than at the bottom layer, seem to be also roughly consistent with the observations made at the Yangtze River estuary (Guo et al., 2016, Fig. 6).However, both in the Yangtze River estuary and in the Rhine ROFI, some observations indicate higher sinking rates at the deeper layers during spring (Peperzak et al., 2003;Guo et al., 2016).This is explained by the species-specific differences in sinking rates (e.g., with sinking rates of diatoms being higher than those of dinoflagellates and Phaeocystis); therefore, the differences in community composition at the surface (dominated by dinoflagellates or Phaeocystis) and bottom layers (dominated by di- atoms) cannot be captured by the presented single-species model.
The vertical distribution of chlorophyll qualitatively depends on the formulation of phytoplankton light availability, sinking and the resource utilization traits of phytoplankton.To demonstrate this, we considered alternative parameterizations regarding the sinking of phytoplankton and light climate and a non-acclimative model version where the physiology of phytoplankton is fixed (see sect.B3 for a detailed description), and compare the resulting chlorophyll profiles with the original (reference) model run for 3 example days and at a deep spot inside the German Bight (Fig. B4).The vicinity of this time interval and location (55 • N, 5 • E) was characterized by the occurrence of a thin chlorophyll layer according to the Scanfish data, as captured quite realistically by the reference run (Fig. 9).
The model run with constant and low sinking rate also resulted in deep chlorophyll maxima for the first 2 days (Fig. B4a), which is, however, not concentrated at the thermocline as in the reference run (Fig. 9), but rather close to the surface with a wider vertical distribution, and on the 3rd day a monotonic profile with the higher values homogeneously distributed within the upper mixed layer.The model with constant and high sinking rate, on the other hand, resulted in profiles monotonically increasing towards the bottom, with overall low concentrations (Fig. B4a).The model run that assumed a spatially constant, low-background attenuation values characteristic of clear ocean waters resulted in a sharp increase in chlorophyll concentrations at around 15 m depth as in the case of the reference run (Fig. B4b), but then the concentrations did not decrease towards the bottom significantly, unlike in the case of the reference run.Higher specific attenuation coefficients resulted in distinct deep chlorophyll maxima in the first 2 days, although about 5 m closer to the surface, and on the 3rd day in homogeneous distribution within the upper mixed layer, again unlike in the case of the reference run (Fig. B4b).Finally in the simplified model with fixed physiologies for two different parameterizations regarding the allocation to the light-harvesting, nutrient acquisition and carboxylation machineries, phytoplankton ends up always being concentrated at the bottom layers (Fig. B4c).2. Instead of being optimized, f A i was fixed to a constant value of 0.5, implying equal investments into affinity and intracellular transport.

B3 Non-acclimative model
3. The uptake activity function, a i , was replaced with a classical linear function of individual cellular quotas In this simplified, fixed-trait model, f LH and f C , which are dynamic state variables in the full model, become parameters.We considered two conceptual assumptions for assigning their values: (1) balanced allocations to each cellular machinery (referred to as "F-bal" in Figs. 14 and B4c), achieved by setting f LH = f C = 0.333; and (2) assigning the domain-wide, volume-weighted averages obtained with the acclimative model for a specific time period (referred to as "F-avg(R)" in Figs. 14 and B4c), which were, for the year 2010 (for which the results are compared), f LH = 0.38 and f C = 0.24.The further two parameters, namely the upper bounds of nitrogen and phosphorus quotas, were set to Q max N = 0.35 and Q max P = 0.04 such that the resulting ranges of N : C and P : C ratios are similar to those obtained with the full model for the year 2010, for which the results are compared.

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 .
Figure2.Structure of the biogeochemical model.Model components (rectangles) comprise the following.B is phytoplankton, Z is zooplankton, POM and DOM are particulate and dissolved organic matter, DIM (N, P) is dissolved inorganic matter (nitrogen and phosphorus), and bAP is the P adsorbed in iron-phosphorus complexes (see Sect. 2.2.1 and Appendix A for further details).C, N and P in small circles refer to carbon, nitrogen and phosphorus bound to each component, respectively, whereas f LH and f C are the allocation coefficients for light harvesting and carboxylation (Sect.2.2.1).Boxes in dashed lines indicate model forcing.

Figure 3 .
Figure 3. Salinity (PSU) measured by FerryBox (a) and estimated by the model (b) along the route shown in the inset.Note that the lower range of salinity was truncated.

Figure 4 .
Figure 4. Comparison of modeled and measured (ICES) temperature (abbreviated T in panels a, b, c, d) and salinity (S in a, b, e, f) at the surface (left) and bottom (right) layers for the period 2006-2010.The 2-D histograms show the number of occurrences of simulationmeasurement pairs.The normalized bias (B * ), Pearson correlation coefficients (ρ) and corresponding number of data points (n) are shown on top of the scatter plots.

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

Figure 7 .
Figure 7.As in Fig. 5, but for the offshore monitoring stations operated by the Bundesamt für Seeschifffahrt und Hydrographie.

Figure 8 .
Figure 8.Comparison of simulated and measured (ICES) DIN (a, b, c, d), DIP (a, b, e, f) and chlorophyll (a, b, g, h) at the surface (left) and bottom (right) layers for the period 2000-2010.The 2-D histograms show the number of occurrences of simulation-measurement pairs.The normalized bias (B * ), Pearson correlation coefficients (ρ) and corresponding number of data points (n) are shown on top of the scatter plots.

Figure 9 .
Figure 9. (a, b) σ T measured by Scanfish and estimated by the model; (c, d) normalized anomalies of fluorescence measured by Scanfish and chlorophyll concentrations estimated by the model.The track of the cruise, which took place between 13 and 19 July 2010, is shown in (e).

Figure 10 .
Figure 10.Average temperature difference (K) between the surface and bottom layers, averaged throughout 2000-2010 for each month.Gray lines show the isobaths.

Figure 11 .
Figure 11.DIN (a, b) and DIP (c, d) concentrations at the surface layer, averaged over the non-growing (months 1-3 and 10-12, left) and growing seasons (months 4-9, right) for the entire simulation period (2000-2010).Concentrations at the bottom layer are almost identical for months 1-3 and 10-12 and similar for months 4-9.Gray lines show the isobaths.Note the different color scales used for each panel and that the scale used for DIN is 16 times that of DIP, such that the identical coloring for DIN and DIP for the same season indicates a Redfield ratio of 16 : 1.

Figure 12 .
Figure 12.Comparison of satellite (ESA-CCI; a, b) and MAECS (c,d) estimates of surface chlorophyll concentrations averaged over 2008-2010 and for the non-growing (months 1-3 and 10-12, left) and growing seasons (months 4-9, right).The 2-D histograms (e, f) show the number of occurrences of simulation-satellite data pairs.Gray lines in (a)-(d) show the isobaths.The normalized bias (B * ), Pearson correlation coefficients (ρ) and corresponding number of data points (n) are shown on top of the scatter plots.

Figure 13 .
Figure 13.Chlorophyll : C ratio in phytoplankton, averaged over the non-growing (a) and growing season (b) of 2010.Gray lines in (a)-(d) show the isobaths.

Figure 14 .
Figure 14.Annual average phytoplankton carbon (and for R, chlorophyll) concentrations in 2010 (a) at the surface layer; (b) averaged over the water column; and obtained with the following models: R is the reference (acclimative) model, F-bal is the fixedphysiology model with balanced investments and F-avg(R) is fixedphysiology model with allocation parameters as average trait values produced by R.

Figure A1 .
Figure A1.Fraction of sorbed fraction of benthic phosphorus as a function 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).
Attenuation coefficient for phytoplankton 0.015 m 2 mmolC −1 3 k POC Attenuation coefficient for POC 0.01 m 2 mmolC −1 3 k DOC Attenuation coefficient for DOC 0.0025 m 2 mmolC −1 4 Figure B1.Total light attenuation along the transect shown in the inset as estimated by the satellite product (ESA-CCI K d , solid gray line) and by the model at the surface layer (Sim K d , solid black line).Light attenuation caused only by SPM used as model forcing is shown separately (Sim K SPM , solid dotted line).All values represent averages between 2008 and 2010.The satellite estimates are from 490 nm wavelength and the model estimates represent the average within the PAR (400-700 nm) spectrum.

Figure B2 .
Figure B2.SPM concentrations measured (gray dots) and the monthly climatologies used as model forcing (black line) at the Noordwijk-10 station (location shown in Fig. 6).

For
gaining insight into the relevance of acclimation aspects of the model, we considered a simplified version of the model in which the adaptive and optimality based features of the model were excluded.For transforming the full model to an otherwise equivalent non-acclimative version the following is needed: 1. Dynamic equations Eqs.(A7)-(A8) that describe the allocation of resources to light-harvesting (f LH ), carboxylation (f C ) and nutrient acquisition (1 − f LH − f C ) traits were excluded.

Table 1 .
Skill scores obtained at each station (B * is the normalized bias, ρ is the Pearson correlation coefficients and n is the number of matching data points) against ICES and ESA-CCI data shown in Figs.5-8 and Fig.12partially (for the averages of months1-3, 10-12 and  4-9).Colors indicate skill level, with red being low, yellow being moderately low, green being moderately high and blue being high (see Sect. 2.3).