the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
FABMNflexPD 1.0: assessing an instantaneous acclimation approach for modeling phytoplankton growth
Prima Anugerahanti
Sherwood Lan Smith
Coupled physical–biogeochemical models can generally reproduce largescale patterns of primary production and biogeochemistry, but they often underestimate observed variability and gradients. This is partially caused by insufficient representation of systematic variations in the elemental composition and pigment density of phytoplankton. Although progress has been made through approaches accounting for the dynamics of phytoplankton composition with additional state variables, formidable computational challenges arise when these are applied in spatially explicit setups. The instantaneous acclimation (IA) approach addresses these challenges by assuming that $\mathrm{Chl}:\mathrm{C}:\mathrm{nutrient}$ ratios are instantly optimized locally (within each modeled grid cell, at each time step), such that they can be resolved as diagnostic variables. Here, we present the first tests of IA in an idealized 1D setup: we implemented the IA in the Framework for Aquatic Biogeochemical Models (FABM) and coupled it with the General Ocean Turbulence Model (GOTM) to simulate the spatiotemporal dynamics in a 1D water column. We compare the IA model against a fully dynamic, otherwise equivalently acclimative (dynamic acclimation; DA) variant with an additional state variable and a third, nonacclimative and fixedstoichiometry (FS) variant. We find that the IA and DA variants, which require the same parameter set, behave similarly in many respects, although some differences do emerge especially during the winter–spring and autumn–winter transitions. These differences however are relatively small in comparison to the differences between the DA and FS variants, suggesting that the IA approach can be used as a costeffective improvement over a fixedstoichiometry approach. Our analysis provides insights into the roles of acclimative flexibilities in simulated primary production and nutrient drawdown rates, seasonal and vertical distribution of phytoplankton biomass, formation of thin chlorophyll layers and stoichiometry of detrital material.
1.1 Modeling phytoplankton and their cellular composition
In early ecosystem models, the elemental composition, i.e., proportion of carbon (C), nitrogen (N) and phosphorus (P) content of phytoplankton, was generally assumed constant, and at least since the work of Dugdale (1967) their growth was typically described by the socalled “Monod” model (Monod, 1949), which assumes a saturating response of the rate of carbon assimilation (and hence of nutrient uptake) to the ambient nutrient concentration, described by a rectangular hyperbolic function. Similarly, specific chlorophyll (Chl) content, i.e., the Chl:C ratio, was assumed to be constant when comparing the simulated phytoplankton biomass against the in situ or satellitebased chlorophyll measurements. In many primary production modules coupled to general circulation models that are actively being used for various purposes to this date, phytoplankton C:Chl and/or $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratios are assumed to be constant (see, e.g., the models in Laufkötter et al., 2015).
The inadequacy of these simplifying assumptions was made clear decades ago by the discovery that phytoplankton elemental composition (e.g., Gerloff and Skoog, 1954) and chlorophyll content (e.g., Platt and Jassby, 1976) are variable. $\mathrm{Chl}:\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratios of phytoplankton have since been found to vary widely in many laboratory experiments (e.g., Kruskopf and Flynn, 2006) and field observations (e.g., Martiny et al., 2013; Burson et al., 2016). Since the work of Caperon (1968) and Droop (1968), the socalled “quota” (or variable internal stores or “Droop”) model has been widely employed to describe the dynamics of carbon and nutrients bound to phytoplankton, using a separate state variable for each element or nutrient resolved. For describing variable Chl:C ratios, acclimation models, most commonly that of Geider et al. (1998), but also others (e.g., Pahlow and Oschlies, 2009; Wirtz and Kerimoglu, 2016) are being increasingly employed in biogeochemical model frameworks. Such models typically couple a description of variable N:C (or other nutrient:C) with photoacclimation, i.e., variation of Chl:C, using one more state variable for Chl bound to phytoplankton (Moore et al., 2002; SchourupKristensen et al., 2014; Kerimoglu et al., 2017; Kwiatkowski et al., 2018). Some models assume a constant N:C ratio, while describing the variations in Chl:C, e.g., using only the photoacclimation portion (e.g., Moore et al., 2004) of the model by Geider et al. (1998) or using an empirical function (e.g., Oschlies and Schartau, 2005) that was earlier proposed by Cloern et al. (1995).
Models that account for variations in cellular composition are in principle more likely to provide more realistic estimates of phytoplankton biomass and biogeochemical fluxes: when the variabilities in Chl:C and C:nutrient ratios are realistically represented by the models, their calibration on the basis of in situ and satellite Chl observations becomes more accurate (Behrenfeld et al., 2009; Ayata et al., 2013; Kerimoglu et al., 2017), and their estimates of biosynthesis rates of C and nutrients, consequently the drawdown of nutrients, and elemental composition of the export flux can be better reproduced (Anderson and Pondaven, 2003; Mongin et al., 2003), respectively. However, the mechanistic basis of some of the models remains questionable, given their parameterization of certain processes using heuristic or empirically inspired functions (Flynn et al., 2015). Moreover, schemes that require additional state variables, due to the need to calculate their transport as tracers, impose substantial computational costs. Especially for models that contain many phytoplankton functional types or clones (e.g., 350 in Dutkiewicz et al., 2020), such additional computational costs may severely limit the kinds of simulations and in silica experiments that can be conducted.
1.2 An optimalitybased resource allocation model
For the prediction of growth, nutrient uptake and acclimative variations of pigment and nutrient content of phytoplankton in response to changes in resource environment, such as the availability of mineral nutrients and light, “resource allocation models” (RAMs) have been used (Shuter, 1979; Laws and Chalup, 1990; Armstrong, 1999; Klausmeier et al., 2004; Pahlow, 2005; Wirtz and Kerimoglu, 2016). This approach is based on the expectation that evolution has produced organisms that strive to maximize their net growth rate by optimally allocating their resources to cellular functions. The dependence of all such functions on common resources therefore implies ecophysiological tradeoffs (Smith et al., 2011). In this study, we specifically consider four physiological variables for describing the acclimative flexibilities involved in phytoplankton growth, as described by Pahlow et al. (2013):
 Q:

N quota (i.e., N:C ratio [molN molC^{−1}]) of phytoplankton.
 f_{V}:

fractional allocation (dimensionless) to the nutrient uptake compartment (protoplast) to optimize the tradeoff between photosynthesis (μ) and nutrient uptake (V), as described by Pahlow and Oschlies (2013).
 f_{A}:

fractional allocation (dimensionless) to affinity to optimize the tradeoff between nutrient affinity (A) and maximum uptake rate within the nutrient uptake compartment (V_{max}), as described by Pahlow (2005) and Smith et al. (2009).
 $\widehat{\mathit{\theta}}$:

Chl:C ratio in chloroplasts ($\widehat{\mathit{\theta}}$, [gChl molC^{−1}]) to optimize the tradeoff between energy gained by light harvesting and energetic costs of chlorophyll synthesis and maintenance, as described by Pahlow et al. (2013).
1.3 Instantaneous acclimation approach
As in most previous models of flexible phytoplankton composition, the abovementioned model by Pahlow et al. (2013) explicitly resolved the dynamics of the carbon, nitrogen, and chlorophyll within phytoplankton biomass. This approach is well suited for simulating the shortterm (i.e., hours to days) dynamics of growth and hence for testing model assumptions against the results of batch culture experiments (e.g., Pahlow, 2005; Pahlow et al., 2013). Resolving the transient dynamics is important for such shortterm experiments, where the response of phytoplankton may differ substantially in terms of nutrient uptake vs. carbonbased growth and chlorophyll synthesis.
By contrast, oceanic (or even freshwater) observations are rarely available at such fine temporal resolution. The lack of observations at sufficient temporal resolution to test shortterm model dynamics motivated the development of the instantaneous acclimation (IA) approach (Smith et al., 2016) as a way to potentially capture growth response at longer timescales while requiring substantially fewer calculations. IA is based on the balanced growth assumption, which Burmaster (1979) showed was able to reconcile the ability of the Droop, Monod and Michaelis–Menten models to capture phytoplankton growth response at steady state, as measured by continuous culture experiments. The key assumption is that growth and nutrient uptake are at all times strictly balanced with respect to the internal C:N stoichiometry of the cell (see Sect. 2.2 below for details). Based on this assumption, IA calculates only one specific rate for both growth and nutrient uptake. Smith et al. (2016) applied this assumption in a 0D (box) model, adequate for reproducing sparse oceanic observations, but did not evaluate its performance compared to fully dynamic models of flexible composition.
Ward (2017) compared the results of a phytoplankton model with instantaneously adjusting quota against a fully dynamic model with explicit state variables for each element resolved, and a fixedstoichiometry model, in a 0D setup. He found that for a wide range of realistic forcing dynamics, the instantaneous approach yielded results practically indistinguishable from the fully dynamic model, whereas these results differed considerably from those of the fixedstoichiometry model. To our knowledge, the IA approach has yet to be tested in a spatially explicit model, where the inclusion of transport terms may lead to additional complications. In a spatially structured environment, transport of cells with a certain internal state to a zone where the typical (average) cellular composition differs, can result in a spatial storage advantage (Grover, 2009). A typical example of this is nutrientreplete cells (as represented by high N:C) at the deeper layers diffusing towards the surface mixed layer (SML) across the thermocline where the cells are typically nutrient starved (e.g., Kerimoglu et al., 2012). In principal, this effect can be resolved only by explicitly tracing the constituents of the cell dynamically.
1.4 Objectives of this study
This study presents a novel implementation of the IA approach in the Framework for Aquatic Biogeochemical Models (FABM; Bruggeman and Bolding, 2014) , and an assessment of its behavior compared to two other established variants (Fig. 1): the first is the widely used, nonacclimative fixedstoichiometry (FS) variant, which resolves only the N bound to phytoplankton explicitly. The second variant is the dynamic acclimation (DA) variant, which resolves the C and N bound to phytoplankton fully dynamically, with two state variables. The comparisons of the three model variants were conducted to answer the following two specific questions: (i) how do the simulations performed with the IA variant differ from those of the fully dynamic DA variant, and (ii) compared to the FS variant, do the results of the IA variant differ substantially? While answering these questions, we aimed to gain mechanistic understanding of the dynamics driving the difference between the model results.
In the following sections, we describe the general structure of the model, the details of the physiological flexibilities mentioned above for each model variant, and the setup to simulate the model. Then we show the results of the simulated patterns of phytoplankton in terms of carbon, nitrogen, chlorophyll, cell quota (Q) and Chl:C ratio, as well as the fractional allocations. We finally discuss the advantages, as well as the challenges and limitations of implementing the IA approach.
2.1 General structure
For describing the cycling of N, we consider a simple model structure (Fig. 1) with four compartments: C and N bound to phytoplankton (Phy_{C}, Phy_{N}), detritus (Det_{C}, Det_{N}), dissolved organic matter (DOC, DON) and dissolved inorganic nitrogen (DIN). Note that our model does not resolve the dynamics of dissolved inorganic carbon (DIC) per se.
The coupled set of differential equations (s(x) short for $\frac{\mathrm{d}x}{\mathrm{d}t}$) that describe the dynamics of state variables are provided in Eqs. (1)–(4), where each of the coupled C or N terms are annotated with the processes they represent. The formal definition and exact formulation of the flux terms (F_{FROM−TO}) in Eqs. (1)–(4) that are trivial (i.e., all except ${F}_{\mathrm{DIN}{\mathrm{Phy}}_{\mathrm{N}}}$ and ${F}_{\mathrm{DIC}{\mathrm{Phy}}_{\mathrm{C}}}$) are provided in Table 1. For equations applying only to a subset of our model variants, the variants are indicated near the equation number in curly brackets ({}). In addition, Table 2 provides an overview of how the model variants differ.
It should be noted that the Phy_{C} is resolved as a state variable only by the DA variant (Eq. 1b). The terms ${F}_{\mathrm{DIN}{\mathrm{Phy}}_{\mathrm{N}}}$ and ${F}_{\mathrm{DIC}{\mathrm{Phy}}_{\mathrm{C}}}$ have central importance to this study and deserve explanation. ${F}_{\mathrm{DIN}{\mathrm{Phy}}_{\mathrm{N}}}$ represents the net N flux from the DIN to phytoplankton and is given by the product of the phytoplankton carbon biomass, Phy_{C} and the specific nutrient uptake rate, V:
For the FS and IA variants, balanced growth (Burmaster, 1979) is assumed, such that V is directly linked to net growth rate, μ, via the nutrient quota, Q:
whereas for the DA variant, V is calculated explicitly (Eq. 12). Net growth rate, μ, is obtained by subtracting the respiration costs associated with chlorophyll maintenance and synthesis, R_{Chl}, and nutrient uptake, R_{N}, from the cellular gross growth rate, μ_{g} (Eq. 13):
where ζ_{N} is the cost of N assimilation (Table 3) and R_{Chl} is the cost of chlorophyll synthesis and maintenance (Sect. 2.2.4).
${F}_{\mathrm{DIC}{\mathrm{Phy}}_{\mathrm{C}}}$ is required only by the DA variant that explicitly resolves the dynamics of Phy_{C} (Eq. 1b). It is given by the product of net growth rate, μ with Phy_{C}, as is typical in quota models (Caperon, 1968; Droop, 1968):
2.2 Flexibilities represented by the model variants
We compare the behavior of three model variants that differ in their representation of the physiological flexibilities. These variants are as follows:

DA explicitly describes the dynamics of nitrogen and carbon bound to phytoplankton, and the acclimation mechanisms introduced in Sect. 1.2, here as represented by flexibilities in growth vs. nutrient uptake; nutrient affinity vs. maximum uptake; and chlorophyll density in chloroplasts; each of which are explained in detail in the following sections. A full description of this variant (including diazotrophy) can be found in Pahlow et al. (2013).

IA assumes that the nitrogen quota (molar N:C ratio) adjusts instantaneously to its optimal value locally (i.e., at any point in time and space) but is otherwise identical to the DA variant with respect to the acclimation mechanisms. A full description of this variant can be found in Smith et al. (2016).

FS assumes no physiological acclimation or quota variability whatsoever.
In the following, representations of the acclimative flexibilities by each model variant are explained in detail.
2.2.1 Flexibility I: nutrient quota
Flexibility in the elemental composition of phytoplankton (Q) is a result of acclimation processes, such as synthesis of enzymes or pigments, which differ in elemental composition (e.g., Geider and La Roche, 2002), in response to changes in resource (light and nutrient) availability.

For the dynamic acclimation variant, Q, is the ratio of the phytoplankton N and C state variables:
$$\begin{array}{}\text{(9)}& Q={\displaystyle \frac{{\mathrm{Phy}}_{\mathrm{N}}}{{\mathrm{Phy}}_{\mathrm{C}}}}\mathit{\left\{}\mathrm{DA}\mathit{\right\}}.\end{array}$$ 
For the instantaneous acclimation variant, Q is assumed to adjust instantaneously to its balancedgrowth optimum (Q^{o}) according to Pahlow and Oschlies (2013):
$$\begin{array}{}\text{(10)}& {Q}^{\mathrm{o}}={\displaystyle \frac{{Q}_{\mathrm{0}}}{\mathrm{2}}}\left[\mathrm{1}+\sqrt{\mathrm{1}+{\displaystyle \frac{\mathrm{2}}{{Q}_{\mathrm{0}}({\widehat{\mathit{\mu}}}_{\mathrm{net}}/\widehat{V}+{\mathit{\zeta}}_{\mathrm{N}})}}}\right]\mathit{\left\{}\mathrm{IA}\mathit{\right\}},\end{array}$$where ${\widehat{\mathit{\mu}}}_{\mathrm{net}}$ and $\widehat{V}$ are the chloroplastspecific net growth and protoplastspecific N uptake rates (Table 2), and Q_{0} and ζ_{N} and are the subsistence quota, and cost of N uptake (Table 3), respectively. Note that this solution differs slightly from the solution proposed by Smith et al. (2016), where the cost of chlorophyll maintenance and synthesis was ignored (see Appendix A for details).

In the fixedstoichiometry variant, Q is a prescribed parameter (Table 2).
2.2.2 Flexibility II: growth vs. nutrient uptake
Given the high nitrogen content in the enzymes responsible for both CO_{2} fixation and nutrient uptake and assimilation (Geider and La Roche, 2002), we consider a tradeoff in the allocation of nitrogen between carbon fixation and nutrient uptake for the acclimative variants, whereas this tradeoff is ignored for the FS variant.

For the acclimative variants (DA and IA), following Pahlow and Oschlies (2013), the tradeoff is specified in terms of the fraction of cellular nitrogen reserves allocated to nitrogen uptake (f_{V}), which linearly increases V and decreases μ_{g}, through decreasing the resources available for carbon fixation, f_{C}, which is interpreted as the relative size of the chloroplast (Pahlow and Oschlies, 2013).
$$\begin{array}{}\text{(11)}& {f}_{\mathrm{C}}=\left(\mathrm{1}{\displaystyle \frac{{Q}_{\mathrm{0}}}{\mathrm{2}Q}}{f}_{V}\right)\mathit{\{}\mathrm{IA},\mathrm{DA}\mathit{\}},\end{array}$$where f_{V} is the fractional allocation towards nutrient uptake for the DA variant (see Eq. 6 for IA variant):
$$\begin{array}{}\text{(12)}& V={f}_{V}\cdot \widehat{V}\mathit{\left\{}\mathrm{DA}\mathit{\right\}},\end{array}$$where $\widehat{V}$ is the protoplastspecific N uptake rate (see below). The cellular gross growth rate is then determined by scaling the gross growth rate within the chloroplast ${\widehat{\mathit{\mu}}}_{\mathrm{g}}$ (see Sect. 2.2.4) by the relative size of the chloroplast, f_{C}:
$$\begin{array}{}\text{(13)}& {\mathit{\mu}}_{\mathrm{g}}={f}_{\mathrm{C}}\cdot {\widehat{\mathit{\mu}}}_{\mathrm{g}}.\end{array}$$Note that, for calculating the effective flux from DIN to Phy_{N} (Eq. 5), only the DA variant uses V as calculated by Eq. (12), while the IA variant calculates the uptake rate from the growth rate, based on the balanced growth assumption (Eq. 6). However, the IA variant still needs the V as calculated by Eq. (12) for calculating the costs of nutrient uptake (Eq. 7).
Both acclimative variants assume that f_{V} maximizes the net specific growth rate under balanced growth conditions. Following Pahlow and Oschlies (2013), this optimal value is found as (see Appendix A)
$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\mathrm{d}\mathit{\mu}}{\mathrm{d}{f}_{V}}}=\mathrm{0}\Rightarrow {f}_{V}=\left({\displaystyle \frac{{Q}_{\mathrm{0}}}{\mathrm{2}Q}}\right){\mathit{\zeta}}_{\mathrm{N}}(Q{Q}_{\mathrm{0}})\mathit{\{}\mathrm{IA},\mathrm{DA}\mathit{\}}.\end{array}$$ 
For the fixedstoichiometry variant, the gross growth rate, μ_{g} is obtained by the multiplication of ${\widehat{\mathit{\mu}}}_{\mathrm{g}}$, for FS, interpreted as the lightlimited potential growth rate, with a nutrient limitation term L_{N}, formulated as a hyperbolic function of ambient DIN concentration, following the Michaelis–Menten–Monod model (Johnson and Goody, 2011; Monod, 1949):
$$\begin{array}{}\text{(15)}& {\mathit{\mu}}_{\mathrm{g}}={\widehat{\mathit{\mu}}}_{\mathrm{g}}\cdot {L}_{\mathrm{N}}={\widehat{\mathit{\mu}}}_{\mathrm{g}}\cdot {\displaystyle \frac{\mathrm{DIN}}{{K}_{\mathrm{N}}+\mathrm{DIN}}}\mathit{\left\{}\mathrm{FS}\mathit{\right\}}.\end{array}$$Thus, for the FS variant, μ (Eq. 7), and hence, through the balanced growth assumption, V (Eq. 6), are directly linked to the external nutrient concentration (Eq. 15) as in typical fixedstoichiometry models. Given the fact that both L_{N} (Eq. 15) for the FS variant and f_{C} (Eq. 11) for the acclimative variants have an equivalent role (in scaling ${\widehat{\mathit{\mu}}}_{\mathrm{g}}$ to μ_{g}), and they both represent nutrient limitation, we consider them to be comparable, i.e., L_{N}∼f_{C}.
2.2.3 Flexibility III: nutrient affinity vs. maximum uptake rate

Originally introduced for describing the substrate uptake by bacteria, “affinity” of a microorganism “can be viewed as a measure of effective collusion between substrate and transport site” (Button, 1978), which can be practically found from the initial slope (i.e., before saturation) of the uptake rate with respect to the substrate concentration (Button, 1978). The term has been used for describing the nutrient uptake by phytoplankton (Aksnes and Egge, 1991) and recognized to be a measure of competitive ability under low concentrations. The maximum nutrient uptake rate, on the contrary, can taken to be a measure of competitiveness under high nutrient concentrations. The protoplastspecific N uptake rate, $\widehat{V}$, can be described by a function of maximal uptake rate, ${\widehat{V}}_{max}$, and nutrient affinity, $\widehat{A}$:
$$\begin{array}{}\text{(16)}& \widehat{V}={\displaystyle \frac{{\widehat{V}}_{max}\cdot \widehat{A}\cdot \mathrm{DIN}}{{\widehat{V}}_{max}+\widehat{A}\cdot \mathrm{DIN}}}\mathit{\{}\mathrm{IA},\mathrm{DA}\mathit{\}}.\end{array}$$The acclimation variants introduce a tradeoff between affinity vs. maximum uptake rate. This tradeoff is captured by the fractional allocation of resources to affinity (f_{A}), which increases affinity, $\widehat{A}={f}_{A}{\widehat{A}}_{\mathrm{0}}$, while decreasing maximum uptake rate, ${\widehat{V}}_{max}=(\mathrm{1}{f}_{A}){\widehat{V}}_{\mathrm{0}}$, so that Eq. (16) becomes
$$\begin{array}{}\text{(17)}& \widehat{V}={\displaystyle \frac{(\mathrm{1}{f}_{A}){\widehat{V}}_{\mathrm{0}}\cdot {f}_{A}{\widehat{A}}_{\mathrm{0}}\cdot \mathrm{DIN}}{(\mathrm{1}{f}_{A}){\widehat{V}}_{\mathrm{0}}+{f}_{A}{\widehat{A}}_{\mathrm{0}}\cdot \mathrm{DIN}}}\mathit{\{}\mathrm{IA},\mathrm{DA}\mathit{\}}.\end{array}$$f_{A} is set to its optimum value, which maximizes $\widehat{V}$, and hence also V (Pahlow, 2005):
$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\mathrm{d}\widehat{V}}{\mathrm{d}{f}_{A}}}=\mathrm{0}\Rightarrow {f}_{\mathrm{A}}={\left[\mathrm{1}+\sqrt{{\displaystyle \frac{{\widehat{A}}_{\mathrm{0}}\cdot \mathrm{DIN}}{{\widehat{V}}_{\mathrm{0}}}}}\right]}^{\mathrm{1}}\mathit{\{}\mathrm{IA},\mathrm{DA}\mathit{\}}.\end{array}$$ 
The fixedstoichiometry variant ignores this tradeoff entirely, by describing the nutrient limitation with the Michaelis–Menten–Monod function (Eq. 15). Following Button (1978) and Smith et al. (2009), the K_{N} parameter in Eq. (15) can be expressed as a function of V_{max} and $\widehat{A}$, according to
$$\begin{array}{}\text{(19)}& {K}_{\mathrm{N}}={\displaystyle \frac{{\widehat{V}}_{max}}{\widehat{A}}}={\displaystyle \frac{(\mathrm{1}{f}_{A})\cdot {\widehat{V}}_{\mathrm{0}}}{{f}_{A}\cdot {\widehat{A}}_{\mathrm{0}}}}\mathit{\left\{}\mathrm{FS}\mathit{\right\}}.\end{array}$$Based on Eq. (19), corresponding K_{N} values were diagnosed from the solution of the IA variant (i.e., using the locally optimized f_{A} values as calculated with Eq. (18), and ${\widehat{A}}_{\mathrm{0}}$ and ${\widehat{V}}_{\mathrm{0}}$ parameters specified for the IA and DA variants). The biomassweighted spatiotemporal average K_{N} value so obtained was prescribed for the FS variant (Table 3).
2.2.4 Flexibility IV: photoacclimation
Photoacclimation is based on the net carbon fixation rate within the chloroplast, ${\widehat{\mathit{\mu}}}_{\mathrm{net}}$ (equivalent to 𝒜 in Pahlow and Oschlies, 2013), which is obtained by subtracting the chloroplastspecific synthesis and maintenance costs of chlorophyll, from the gross growth rate within the chloroplast, i.e.,
where ${\widehat{\mathit{\mu}}}_{\mathrm{g}}$ is given by the product of day length as a fraction of 24 h, L_{D}, potential turnover rate, ${\widehat{\mathit{\mu}}}_{\mathrm{0}}$, and the lightsaturation of the photosynthetic apparatus, L_{I}:
L_{I} is a saturating function of daytime average light, $\overline{I}$, and chlorophyll density in chloroplasts, $\widehat{\mathit{\theta}}$:
where α is light affinity. Returning to Eq. (20), ${\widehat{R}}_{\mathrm{Chl}}$ is given by
where ${R}_{\mathrm{M}}^{\mathrm{Chl}}$ and ζ_{Chl} are the costs of chlorophyll maintenance and synthesis, respectively (Table 3).
Photoacclimation is mainly represented in terms of the chlorophyll density in chloroplasts, $\widehat{\mathit{\theta}}$. Increasing $\widehat{\mathit{\theta}}$ reduces light limitation (Eq. 22) but at the expense of greater respiration costs (Eq. 23). In turn, for obtaining the cellular Chl:C ratio, θ is calculated by multiplying $\widehat{\mathit{\theta}}$ times f_{C}, i.e., size of the chloroplast:
Similarly, the overall respiratory cost of maintaining cellular chlorophyll is obtained by multiplying the chloroplastspecific cost by the size of the chloroplast:
Although θ (Eq. 24) is only a diagnostic quantity, R_{Chl} (Eq. 25) directly determines the net growth rate through Eq. (7). Therefore, scaling of the chloroplastspecific respiration rate, ${\widehat{R}}_{\mathrm{Chl}}$, by f_{C} can be considered to be an acclimative quality implied by variable f_{V} and Q, which, in combination (Eq. 11), determine the chlorophyll maintenance cost through Eq. (25).

In the acclimation variants (IA and DA), $\widehat{\mathit{\theta}}$ is assumed to adjust instantaneously to its optimal value, which maximizes ${\widehat{\mathit{\mu}}}_{\mathrm{net}}$. Following Pahlow et al. (2013), this optimal value is
$$\begin{array}{}\text{(26)}& \begin{array}{rl}& \widehat{\mathit{\theta}}=\\ & \left\{\begin{array}{ll}{\displaystyle \frac{\mathrm{1}}{{\mathit{\zeta}}_{\mathrm{Chl}}}}+{\displaystyle \frac{{\widehat{\mathit{\mu}}}_{\mathrm{0}}}{\mathit{\alpha}\overline{I}}}\left(\mathrm{1}{W}_{\mathrm{0}}\left[\left(\mathrm{1}+{\displaystyle \frac{{R}_{\mathrm{M}}^{\mathrm{Chl}}}{{L}_{\mathrm{D}}{\widehat{\mathit{\mu}}}_{\mathrm{0}}}}\right)\mathrm{exp}\left(\mathrm{1}+{\displaystyle \frac{\mathit{\alpha}\overline{I}}{{\widehat{\mathit{\mu}}}_{\mathrm{0}}{\mathit{\zeta}}_{\mathrm{Chl}}}}\right)\right]\right),& \overline{I}{\overline{I}}_{\mathrm{C}}\\ \mathrm{0},& \overline{I}\le {\overline{I}}_{\mathrm{C}}\end{array}\right)\\ & \mathit{\{}\mathrm{IA},\mathrm{DA}\mathit{\}},\end{array}\end{array}$$where W_{0} is the 0 branch of the Lambert W function, $\overline{I}$ is the daytime average irradiance (i.e., $\widehat{I}={\overline{I}}_{\mathrm{24}\mathrm{h}}/{L}_{\mathrm{D}}$), and ${\overline{I}}_{\mathrm{C}}$ is the critical daytime average irradiance level, above which chlorophyll synthesis is worthwhile (Pahlow et al., 2013):
$$\begin{array}{}\text{(27)}& {\overline{I}}_{\mathrm{C}}={\displaystyle \frac{{\mathit{\zeta}}_{\mathrm{Chl}}{R}_{\mathrm{M}}^{\mathrm{Chl}}}{\mathit{\alpha}{L}_{\mathrm{D}}}}.\end{array}$$ 
For the fixedstoichiometry variant, $\widehat{\mathit{\theta}}$ is prescribed as the biomassweighted average value calculated by the IA variant. Considering that θ is typically a constant “conversion factor” in classical, fixedstoichiometry and fixedChl:C models, in Eqs. (24) and (25), we assume that the size of the chloroplast, f_{C}, is constant too. For the sake of consistency with the IA variant, f_{C} for FS is diagnosed from its expanded form, i.e., $\mathrm{1}\frac{{Q}_{\mathrm{0}}}{\mathrm{2}Q}{f}_{V}$ (Eq. 11). Hence, in addition to the prescribed value of Q (see Sect. 2.2.1), the biomassweighted mean of f_{V}, as calculated by the IA variant is prescribed (Table 3). Given the comparability of the terms (Sect. 2.2.2), diagnosing f_{C} from L_{N} comes into question, which is elaborated on in Appendix B.
2.2.5 Temperature scaling
Kinetic rate constants (m, r_{hyd}, r_{rem} in Table 1, and ${\widehat{V}}_{\mathrm{0}}$, $\widehat{A}$, ${\widehat{A}}_{\mathrm{0}}$ and ${R}_{\mathrm{M}}^{\mathrm{Chl}}$ in Table 3) are prescribed for a reference temperature of T_{r}=20 ^{∘}C =293.15 K and scaled to the ambient temperature in water, T (in K), according to the Arrhenius function:
where the gas constant R=8.3145 $\mathrm{J}\phantom{\rule{0.125em}{0ex}}{\mathrm{mol}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}$, and the activation energy, ${E}_{\mathrm{a}}=\mathrm{4.82}\times {\mathrm{10}}^{\mathrm{4}}$ J mol^{−1}, such that every 10^{∘} increase or decrease in T approximately doubles or halves the reference rates.
2.3 Coupling with the hydrodynamical host
The model is implemented in FABM (Bruggeman and Bolding, 2014), so that it can be used, without modification, in combination with various hydrodynamical hosts. In this study, we performed simulations of an idealized water column, using the General Ocean Turbulence Model (i.e., GOTM Burchard et al., 2006). GOTM calculates and provides the relevant physical quantities, such as I (needed in Eq. 22) and T (needed in Eq. 28). I is attenuated with depth (z) by various substances in water, according to
where A, η_{1} and η_{2} represent the differential attenuation length scales of red and blue light (Burchard et al., 2006), and k_{i} is the specific attenuation coefficient of the biological quantities, which we set as 0.03 m^{2} mmolN^{−1} for Phy_{N} and Det_{N}. In order to account for background attenuation, we set the “light extinction method” to “Jerlov type IB”, corresponding to A=0.67 η_{1}=1.0 m, η_{2}=17 m, characterizing water of medium clarity (Paulson and Simpson, 1977). Our results are qualitatively insensitive to these parameter settings. Besides providing necessary environmental variables, GOTM calculates the transport rates of the biological quantities, according to the general equation (Burchard et al., 2006):
where K_{z} is the eddy diffusivity calculated by GOTM, the source terms, s(c_{i}), correspond to Eqs. (1)–(4), and advection rates, w_{i}, are all set to 0.0, except that of detritus for which a sinking rate of −2.0 m d^{−1} was specified. Note that the latter value was arbitrarily chosen to induce a downward flux in this idealized setup, and that in reality, it depends on the average size and density of detritus particles being modeled and displays a vast range (Guidi et al., 2008).
2.4 Idealized setup and simulations
We consider an idealized water column of 100 m depth. In order to mimic an environment that is characterized by strong seasonality, with deep mixed layers in spring and summer stratification, we force the model with astronomically calculated shortwave radiation at 60^{∘} N latitude and a repeating annual cycle of air temperature that ranges between 4–20 ^{∘}C as described by a scaled sinusoidal function (Fig. 2).
All other meteorological variables (wind speed, air pressure, humidity and cloud cover) are assumed to be constant, and the model ignores precipitation and evaporation losses, as well as tidal variations. Starting from initial conditions and annually repeating meteorological forcing as described above, each model variant was run for 3 years. The thirdyear results were nearly identical to those for the second year, indicating that an equilibrium annual cycle was reached. In the following, we elaborate the seasonal dynamics during the third year.
Daytimeaveraged irradiance, $\overline{I}$ and water temperature T simulated by different model variants are very similar with subtle differences (Fig. 3a, d vs. b, e, vs. c, f), because each variant calculates slightly different phytoplankton biomass (see below), resulting in differences in attenuation of light and associated heating. Seasonal and vertical distributions of DIN as estimated by the model variants are similar (Fig. 3g–i). DIN depletion (<1 mmolN m^{−3}) during summer is confined to the upper 25 m as estimated by the FS variant, whereas it extends 5–10 m deeper, as estimated by the IA and DA variants.
With all three model variants, phytoplankton growth patterns are characterized by an intense surface bloom in spring, followed by a gradual deepening of the biomass maxima (Fig. 4a–c). Biomass concentration as estimated by the IA and DA variants during summer is greater than with the FS variant (Fig. 4a, b, c). Compared to the FS variant, the acclimation response in the other two variants tends to produce steeper gradients over both depth and time, because of combined dependencies on the three dynamically optimized allocation factors (f_{A}, f_{V}, and $\widehat{\mathit{\theta}}$). This effect is most pronounced for Phy_{Chl}, which differs the most between the FS and the other two variants. With the FS variant, given the constant N:C (Q) and cellular Chl:C (θ) (Fig. 4g, m), C, N and Chl bound to phytoplankton clearly display identical patterns (Fig. 4a, d, j; note that apparent differences in contour plots are due to contour limits not matching these ratios). IA and DA on the other hand simulate slightly different patterns for C, N and Chl bound to phytoplankton (Fig. 4b, e, k and c, f, l), because of the seasonally and vertically variable $\mathrm{Chl}:\mathrm{C}:\mathrm{N}$. Decoupling of Phy_{N} from Phy_{C} is mainly monotonic and is driven by increasing Q with depth (Fig. 4h–i). On the other hand, decoupling of Phy_{Chl} from Phy_{C} follows a more complex pattern, because of the unimodal distribution of θ across the water column (Fig. 4n–o). As a result of this unimodality, Chl simulated by the IA and DA variants forms a distinct, thin layer below the thermocline (Fig. 4k–l).
During summer, $\widehat{\mathit{\theta}}$ follows a complex but roughly unimodal distribution across depth (Fig. 5b–c): intermediate values at the surface first increase with depth to reach a maximum and then sharply decrease with increasing depth. The low values of $\widehat{\mathit{\theta}}$ towards the surface reflect the optimization, which reduces pigment density when light is abundantly available because of the chloroplastspecific respiratory costs $\widehat{\mathit{\theta}}$ (Eq. 23). This can be seen in the flattening of the lightsaturation function L_{I} (Eq. 22). In the deep layers, as $\overline{I}$ approaches ${\overline{I}}_{\mathrm{C}}$, irradiance becomes insufficient to support the synthesis and maintenance of chlorophyll, and $\widehat{\mathit{\theta}}$ rapidly converges to 0. f_{A} and f_{V} simulated by the IA and DA variants (Fig. 5e–f, h–i) increase with nutrient limitation (Fig. 5j–l) as expected (Smith et al., 2016). The fraction of resources available for carbon fixation, f_{C}, displays a similar pattern in all model variants and is roughly the inverse of f_{V}: high during winter throughout the water column and in the deeper layers during summer and low in the upper layers during summer (Fig. 5j–l). For the FS variant, the pattern of the nutrient limitation term, L_{N}, is similar to the patterns of f_{C} for IA and DA variants (Fig. 5), although its magnitude in the summer is higher than other variants, as can be explained by the incomplete DIN depletion (Fig. 3g; see below). Light saturation of photosynthesis, L_{I}, displays a similar pattern in all variants (Fig. 5m–o) and mainly reflects irradiance levels (Fig. 3a–c). However, compared to the FS variant, the intermediate L_{I} values in the IA and DA variants penetrate deeper (Fig. 5n, o vs. m), because the optimization of $\widehat{\mathit{\theta}}$ enhances light harvesting ability at these intermediate depths (Fig. 5b, c).
During winter and spring blooms, the net cellular growth rate, μ, as estimated by the FS variant, temporarily exceeds those estimated by the acclimative variants (Fig. 6a–c; see below for the explanation). The IA and DA variants estimate higher nutrient uptake rates, V, in surface layers during the spring bloom, and in deeper layers during summer (Fig. 6d–f). Negative V in the bottom layers as estimated by the FS and IA variants is a direct result of the balanced growth assumption (Eq. 6) and can be interpreted as exudation. Respiratory costs of nutrient uptake, R_{N}, (Fig. 6h–i) are much lower than R_{Chl} (Fig. 6j–l). For the FS variant, R_{N} drops below 0 in the deeper (>50 m) waters, implying negative respiration, which is a model artifact, as a result of ${\widehat{\mathit{\mu}}}_{\mathrm{net}}$ becoming negative (see Eq. A4 in Sect. A1) due to the fixed $\widehat{\mathit{\theta}}$. However, these negative values are small and therefore do not have a significant effect on the model results, as evidenced by a sensitivity experiment, where ${\widehat{\mathit{\mu}}}_{\mathrm{net}}$ was constrained to positive values for the FS variant (results not shown). In comparison to the acclimative variants, R_{Chl} of the FS variant is smaller during the spring bloom but larger during summer, the reasons for which are explained below.
For the most part, primary production and relevant dynamics take place within roughly the upper 50 m in the simulated system (Figs. 4–6). A comparison of average quantities in this zone (Fig. 7), in combination with vertical profiles throughout the water column during different times of the year (Fig. 8), as estimated by the three model variants, reveals differences between model variants that are not resolved by the contour plots (Figs. 4–6). In both the IA and DA variants, DIN concentrations are almost entirely depleted before the onset of winter mixing in early November, with minimum concentrations of ∼ 0.005 mmolN m^{−3} near the surface. In the FS variant, DIN remains higher (minimum concentration of ∼0.7 mmolN m^{−3} near the surface (Figs. 7a, 8Ja, Na). Q and f_{C}, as estimated by the IA and DA variants, are nearly identical throughout the season (Fig. 7b, c), but differences arise during winter. For DA, Phy_{C} and Phy_{N}, and hence Q, become vertically homogeneous due to rapid turbulent mixing (Figs. 4c, f, i, 8Fb). However, under the instantaneous acclimation assumption in the IA variant, no matter how well mixed the water column may be, vertical gradients persist for the optimal Q values between the surface and deeper layers during winter (Fig. 8Fb, Nb).
During winter and the spring bloom in March–April, nutrient limitation is almost nonexistent for the acclimative variants, as indicated by f_{C} approaching unity (Fig. 7c), whereas for the FS variant, a degree of nutrient limitation persists (as indicated by L_{N}<f_{C}), due to the saturating behavior of the Monod function to the nutrient concentrations. During late summer (July to October), nutrient limitation becomes less severe for the FS variant than for the acclimative variants in the surface layers (i.e., L_{N}>f_{C}, Fig. 8Jc, Oc). The relatively high L_{N} (minimum: 0.12) of the FS variant results from the incomplete DIN depletion as simulated by the FS variant as mentioned above and the linear response of the Monod function to substrate concentrations at low levels (Eq. 15). In contrast, for the IA and DA variants, Q approaches Q_{0}, and f_{V} approaches its maximum value of 0.5, causing (through Eq. 11) severe nutrient limitation, as f_{C} approaches zero (minimum: 0.005) near the surface.
The cellular net growth rate, μ, as estimated by the FS variant is slightly faster than those of the acclimative variants during winter–spring near the surface (e.g., Fig. 8Ff, Mf) but becomes slower right after the spring bloom (e.g., Fig. 8Af) and stays low throughout the summer (Figs. 7f, 8Jf). It should be noted that the chloroplastspecific growth rate, $\widehat{\mathit{\mu}}$, which is maximized for the acclimative variants through photoacclimative flexibility (Sect. 2.2.4), is always higher than that calculated by the FS variant, as expected (not shown). As the chloroplastspecific chlorophyll maintenance and synthesis costs, ${\widehat{R}}_{\mathrm{Chl}}$ is scaled to the cellular level (through multiplication with f_{C}, Eq. 25), the resulting R_{Chl} for the FS becomes lower than those of the acclimative variants, given that the prescribed f_{C} of the FS variant during this time period is smaller than the dynamically calculated values by the acclimative variants (Figs. 7c; 8Fc, Mc). The lower R_{Chl} of the FS variant, in turn, explains the higher μ during the spring bloom (Fig. 7f). When the chloroplast size of the FS variant is assumed to be proportional to L_{N}, as explained in the Appendix B, estimated growth rate becomes similar to those of the acclimative variants (Fig. B2f). During summer, this effect becomes reversed: high R_{Chl} as estimated by the FS variant in the surface layers (Fig. 6j vs. k–l) contributes to the relatively low μ estimated by this variant (Figs. 7; 8Jf): in addition to the higher $\widehat{\mathit{\mu}}$, the IA and DA variants achieve lower R_{Chl} (Fig. 6j–l) through lower $\widehat{\mathit{\theta}}$ (Fig. 5a–c) and f_{C} (Figs. 7c, 8Jc) at the surface.
During the spring bloom, C bound to phytoplankton, Phy_{C}, simulated by the FS variant exceeds those of the IA and DA variants (Fig. 7d), whereas the differences between the N bound to phytoplankton, Phy_{N}, as simulated by different variants are much smaller (Fig. 7e). This discrepancy between C and N content of phytoplankton is due to the decoupling in the acclimative variants: due to the lower value of the prescribed Q of the FS variant (based on the spatiotemporal average of the values simulated by IA) during winter–spring season (Fig. 7b), a larger amount of C biomass can be synthesized per N taken up in comparison to the acclimative variants, explaining therefore the higher Phy_{C} simulated by the FS. The sensitivity of Phy_{C} of the FS variant is evidenced also by a strong reduction of Phy_{C} (in contrast to relatively unaltered Phy_{N}) during the spring bloom in response to a doubling of the prescribed Q (not shown). During summer, the FS variant estimates considerably lower values of Phy_{C} compared to the IA and DA variants (Fig. 7d), whereas the simulated Phy_{N} concentrations remain similar (Fig. 7e). Therefore, the higher Phy_{C} concentrations simulated by the acclimative variants during this period are promoted by lower Q (Figs. 7b, 8Jb) in the surface layers.
Differences between the IA and DA variants emerge especially right after the spring bloom and autumn destratification. After the spring bloom, growth rate simulated by the IA variant (until May) becomes lower than that by the DA variant (Fig. 7f) near the surface (Fig. 8Af). The main reason for this difference is the slightly lower f_{C} of the IA variant during the winter–spring period (i.e., from December to May) near the surface (Fig. 8Fc, Mc, Ac) except for a short period at the peak of the bloom (Fig. 7c). The lower f_{C} of IA during this period is, in turn, driven by slightly lower Q (Figs. 7b, 8Mb, Ab), which also leads to slightly higher f_{V} (see Eq. 14). As pointed out above, the higher Q simulated by the DA variant before the spring bloom near the surface is maintained by the homogenizing effect of vertical transport (which does not occur with the IA variant), and after the spring bloom following the onset of stratification, the persistently higher Q of the DA variant near the surface reflects the lagged response captured by dynamically tracing C and N content of phytoplankton.
Following the weakening of stratification in early November (Fig. 3), a new phytoplankton bloom develops, especially as reflected by Phy_{N} in all variants but also by Phy_{C} as simulated by the DA variant (Fig. 7d, e). This bloom is driven by the entrainment of DIN and phytoplankton biomass below the thermocline into the SML (compare Fig. 8Oa, d, e vs. Fig. 8Na, d, e). Under these nutrientreplenished conditions, μ is predominantly limited by light, as in winter (Fig. 8Ff), and therefore monotonically increases towards the surface (Fig. 8Nf), as simulated by all variants. On the other hand, vertical distribution of Q as simulated by the IA and DA variants, becomes qualitatively different: due to the rapid turbulent mixing of Phy_{C} and Phy_{N} as simulated by the DA variant, Q is homogeneously distributed within the SML (Fig. 8Nb), but such homogenization does not occur in the IA variant, and Q is determined by the locally optimized f_{V}. Therefore, in the DA variant, a high nutrient uptake at the bottom of the SML (Fig. 8Ng), in combination with mixing within the SML, can support growth near the surface (through Q, Fig. 8Nb), whereas in IA, growth and uptake dynamics are always coupled by definition and determined by local physiological states only, as in the FS variant. The decoupling of (growth and uptake) rates and reshuffling of Q as simulated by the DA variant appear to allow faster uptake of nutrients in comparison to the IA variant within the SML (Fig. 8Ng). A related mechanism potentially contributing to the higher nutrient uptake rates is again a timelag effect: in the DA variant, the nutrientstarved phytoplankton (i.e., the low Q; see Fig. 8Ob) in the SML correspond to a higher nutrient demand.
C:N of detritus, as estimated by the FS variant, approaches a constant equilibrium value throughout the water column by the end of the first year and remains there during the third year (Fig. 9a, d). This is as expected, and this value is simply equal to the reciprocal of the prescribed constant (N:C) quota of phytoplankton, calculated as the biomassweighted average of the Q estimated by the IA variant (Table 3). The C:N ratio of detritus, as estimated by the IA and DA variants, increases during summer (Fig. 9b, c and e, f), driven by the lower phytoplankton quotas during summer (Fig. 4).
Simulated process rates determining ecosystem functioning, such as the watercolumnintegrated net primary production (NPP) and nutrient drawdown (NDD) rates, also differ between the model variants. FS estimates higher NPP rates during winter and the spring bloom (Fig. 10a), consistent with the higher Phy_{C} it estimates during this period (Fig. 7d). While the NPP estimates of IA and DA are very close between the late summer (starting from September) to the spring bloom (in early March), right after the spring bloom, IA estimates suddenly decrease, as a consequence of reduced net specific growth rate, μ (Fig. 7f), as was pointed out above. Interestingly, this difference between the IA and DA is larger than the differences in μ and contrasts with the differences in Phy_{C} averaged over the top 50 m (Fig. 7d, f) but can be explained by the higher vertical covariance between the Phy_{C} and μ in DA than in IA (Fig. 8Ad, Af). Annual average NPP rates as estimated by the FS (48.77 $\mathrm{mmolC}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$) and IA (45.66 $\mathrm{mmolC}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$) variants are, respectively, 8.1 % and 13.9 % smaller than that of the DA variant (53.06 $\mathrm{mmolC}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$). NDD rates (Fig. 10b) are similar during the spring bloom, but the acclimative variants become higher during summer. After the autumn mixing, NDD as simulated by the DA variant shows a spike not well reproduced by the IA and FS variants, which is driven by the fast uptake rates simulated by the DA variant throughout the SML, contrasting with those simulated by the IA variant constrained to the surface layers (Fig. 8Ng). Annual average NDD rate simulated by the DA variant (4.78 $\mathrm{mmolN}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$) is the highest, followed by the 8.2% lower IA (4.39 $\mathrm{mmolN}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$) and 14.3 % lower FS (4.1 $\mathrm{mmolN}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$) variants.
4.1 Modeling variable phytoplankton composition
Elemental composition and pigment density of phytoplankton are known to vary, at both the organismal and community levels (Halsey and Jones, 2015), as demonstrated in the laboratory and under in situ conditions (Moreno and Martiny, 2018). Such variations in phytoplankton and hence detrital C:nutrient ratios have implications for C and nutrient export fluxes, including the functioning of the biological carbon pump in the ocean. Notwithstanding, in many biogeochemical models coupled to general circulation models, primary producers are still unrealistically represented with a constant “Redfield” $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratio and/or constant Chl:C ratio (see, e.g., the models used in Laufkötter et al., 2015). More detailed “quota” models exist; however, these approaches are often challenged by two major limitations: (i) dependence on formulations that lack a clear mechanistic basis, and (ii) their requirement for additional state variables, which increase computational costs.
A concrete example of the first problem, i.e., dependence on heuristic formulations, is the downregulation of nutrient uptake, which is needed to avoid unrealistically high nutrient quotas in a Droop scheme. Often, downregulation is formulated as some function (linear, e.g., Grover, 1991 or nonlinear, e.g., Geider et al., 1998) of “relative quota”, with reference to a prescribed maximum value. The acclimation scheme used in this study (IA and DA variants), requires no such explicit downregulation term nor any prescribed maximum quota value. This is because the optimization of growth, subject to the growth vs. nutrient uptake tradeoff (Sect. 2.2.2), accomplishes this regulation by balancing the marginal benefits of investing into nutrient uptake vs. photosynthesis. This RAM approach, which links various cellular functions via tradeoffs, has proven successful at reproducing various $\mathrm{Chl}:\mathrm{C}:\mathrm{N}:\mathrm{P}$ measurements obtained in laboratory experiments (e.g., Klausmeier et al., 2004; Pahlow et al., 2013; Wirtz and Kerimoglu, 2016). Furthermore, given its mechanistic basis, this approach can be expected to reproduce biological feedbacks more realistically (Flynn et al., 2015).
Earlier studies had pointed out that representation of variables in $\mathrm{Chl}:\mathrm{C}:\mathrm{N}$ ratios of phytoplankton in models resulted in better reproduction of field observations (e.g., Doney et al., 1996; Christian, 2005; Ayata et al., 2013; Chen and Smith, 2018). Consistent with those studies, implementation of the model introduced here for simulating two oligotrophic ocean sites suggested that the portability of phytoplankton growth models is enhanced by the variable cellular composition (Anugerahanti et al., 2021). As demonstrated by these studies, 1D setups, as we also used here, are ideal computational environments for examining the behavior of phytoplankton growth models: while resolving the essential features of aquatic environments, foremost the seasonally variable vertical structuring of resources and transport rates, they increase the computational costs minimally in comparison to the 3D models. On the other hand, realistic representation of the horizontal gradients or investigation of the effects of phytoplankton on the biogeochemical functioning at larger scales does require 3D setups. Recent applications of these models in realistic 3D setups (Kerimoglu et al., 2017; Pahlow et al., 2020) have indicated that accounting for acclimation enhances the ability of models to reproduce field observations and large scale patterns. Moreover, a consistent representation of phytoplankton composition allows identification of potential alterations in trophic transfer efficiencies as mediated by changes in food quality of prey in response to environmental change (Kerimoglu et al., 2018; Kwiatkowski et al., 2018).
Regarding the second problem, i.e., the computational costs of resolving additional state variables, Smith et al. (2016) proposed the instantaneous acclimation approach, according to which the elemental composition of phytoplankton varies but instantaneously, such that tracking these variations do not require additional state variables. As in Smith et al. (2016), we considered the same specific acclimation mechanisms of Pahlow et al. (2013) but under the assumption that the N quota adjusts to an optimal value locally, under strictly “balanced growth” (Burmaster, 1979, see Sect. 2.2.1). While at steadystate, this is a natural consequence of any “Drooplike” model (Burmaster, 1979), assuming this behavior to hold under transient conditions is merely an approximation. Ward (2017), using a classical Droop approach, showed that this approximation holds well under a wide range of conditions in a 0D (box) setup. Here, for the first time, we have tested this approach in an idealized 1D setup and shown that in many respects the IA model and the fully explicit DA variant behave similarly. Our preliminary experiments demonstrated that, even in an environment characterized by periodic perturbations of stratification during summer, the behavior of the two variants remains similar (results not shown). This is significant considering that IA requires only one state variable, whereas DA requires two state variables. Thus, it can be concluded that IA provides improved realism over a computationally equivalent FS approach, which ignores variations in cellular composition. For simulating a few years of the dynamics of the single phytoplankton group in a 1D setup as we did here, differences in computational costs relative to the fully dynamic variant are nearly negligible, but for simulating decades/centuries or millennia in a 3D setup (e.g., as in Pahlow et al., 2020), and/or when multiple phytoplankton clones/types (e.g., with different sizes) are considered (e.g., 350 in Dutkiewicz et al., 2020), differences in computational costs can indeed be substantial.
4.2 Qualitative vs. quantitative differences between model variants
The capacity to store nutrients is known to be an advantageous trait for phytoplankton in temporally fluctuating environments, where greater nutrient storage capacity, e.g., by larger cells, during the nutrientreplete phase provides a competitive advantage during subsequent periods of nutrient scarcity (Grover, 1991; Litchman et al., 2009). Similarly, diffusion or active movement of nutrientrich cells from the nutrientreplete to nutrientrich environments, e.g., from bottom towards surface layers, has been shown to favor species with greater storage capacities (Grover, 2009; Kerimoglu et al., 2012). The IA model presented in this study cannot capture this effect, since according to this approach (i) growth and nutrient uptake rates are always proportional (by definition of the “balanced growth” assumption); thus, differential benefits along a space or time continuum cannot be combined through reshuffling of physiological states. (ii) Nutrient quotas do not have inertia, and hence lagged response, as they are instantaneously adjusted to the local (light and nutrient) resource conditions, unlike in the DA variant where Q is dynamically traced (by the virtue of dynamically tracing both Phy_{C} and Phy_{N}). In fact, the DA variant we considered here presumably has a weaker storage capacity compared to a classical Droop model, because in our acclimative approach, allocation of resources to maximize growth can be expected to suppress “luxury consumption” (Droop, 1968) of nutrients. Finally, it should be noted that because of the differences in the formulation of the uptake in the IA (Eq. 6) and DA (Eq. 12) variants, and the complex interdependencies between the physiological states and process rates involved (Q, f_{V}, f_{C} μ, V, $\widehat{V}$, DIN), comparison of the response of the two variants is not straightforward during such transient phases, and a fuller understanding will require further analysis and experimentation.
Some of the differences in phytoplankton growth dynamics, as simulated by the acclimative IA and DA variants and the nonacclimative FS variant, could be reconciled by tuning the parameters. For instance, the amount of phytoplankton biomass, or the extent of nutrient depletion as simulated by the FS variant, can be increased by specifying higher resource affinities (e.g., lower K_{N} or higher α) to make up for the deficiency in the formulation of lightlimited growth (Oschlies and Schartau, 2005). However, improvements in these specific aspects typically result in greater discrepancies in other aspects, such as the timing and magnitude of the spring bloom, or winter concentrations of nutrients and phytoplankton. In other words, in terms of model performance, tradeoffs exist between multiple objectives. Such tradeoffs become more obvious when attempting to simulate multiple environments characterized by different resource conditions (e.g., multiple sites or the same site in two different time periods) with a single parameter set (Anugerahanti et al., 2021). How acclimative flexibilities impact the sensitivity of models to parameter perturbations remains an open question.
The RAM approach used here, as in “adaptive dynamics” approaches (Follows and Dutkiewicz, 2011), ambiguously reflects processes at multiple organismal scales. For instance, higher f_{A} and f_{V} and lower $\widehat{\mathit{\theta}}$ at the surface layers during summer (Fig. 5), which agrees with lower light harvesting and higher nutrient harvesting investment as found by Bruggeman and Kooijman (2007), can be attributed to (i) evolutionary adaptation of new species (which would be more relevant in a longerterm simulation), (ii) selection among existing species that had been preadapted to these conditions and (iii) individuallevel acclimation. Optimalitybased acclimative models can thus capture some key communitylevel effects of evolutionary and ecological dynamics, without explicitly resolving competing species or groups (Smith et al., 2011). The same idea underlies the recent work of Chakraborty et al. (2020), where they described the changes in community composition by assuming that the trophic strategy of the entire plankton community is optimized instantaneously.
Some features, such as the dense and thin chlorophyll layers at the thermocline as captured by the acclimative variants (Fig. 4), seem qualitatively irreproducible by the FS variant even for a single site and time period. This is because multiple dependencies are necessary for capturing this feature, namely the unimodal distribution of chlorophyll density over depth (Fig. 4) and the steep increase in chloroplast size with depth near the thermocline (Fig. 5), as well as the thermocline being the best compromise between light and nutrient limitation (Fig. 5). The FS variant includes only the last dependency, because it lacks acclimation, and is therefore unable to produce such thin chlorophyll layers. When the chloroplast size is assumed to vary and diagnosed by the nutrient limitation term, such that the vertical Chl:C increases monotonically with depth, the vertical distribution of Chl can be partially captured (Appendix B).
We also found differences in systemlevel metrics such as NPP and NDD (e.g., Bergeron and Tremblay, 2014; Johnson et al., 2017) rates as simulated by different variants. For both metrics, DA estimates were about 10 % higher than the FS and IA variants, with FS estimates systematically skewed towards earlier in the season. It should be noted that, for the FS variant, prescribed Q, which, in this study was based on the results of the IA variant but normally is effectively a free parameter (although the common approach is to set it to the Redfield proportions), largely determines the estimated Phy_{C} and related quantities, such as NPP rates. For instance, doubling the Q of FS results in only a few percent further underestimation (relative to DA) of the NDD rate (annual average: 3.89, instead of the original 4.1 $\mathrm{mmolN}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$, which corresponds to 18.6 % lower than the DA estimate, instead of the original 14.3 %), whereas it leads to more than 50 % lower estimates of NPP rate (23.26 $\mathrm{mmolC}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$) in comparison to that of the DA variant. Some FS variants are based on C and not N as in this study (i.e., the explicit state variable is the C bound to phytoplankton). For those models, instead of the NPP, NDD rates may be more sensitive to prescribed Q. In contrast to the FS variant, with the IA variant, the total C and N content, and growth and nutrient uptake rates of the phytoplankton (thus, the NPP and NDD rates) are determined by the same set of parameters governing the fully explicit DA variant.
4.3 Physiological flexibility and environmental feedbacks
The wellknown links between the composition of phytoplankton and the biogeochemistry of their ambient environments imply feedbacks, which are important in ecology, environmental science and water quality studies. These feedbacks can be mediated by both physiological acclimation and evolutionary adaptation (Moreno and Martiny, 2018), with the latter typically understood to operate on much longer timescales. However, acclimation and adaptation do interact in ecoevolutionary dynamics, and for plankton they may even occur on similar timescales (Smith et al., 2011; Edelaar and Bolnick, 2019). Disentangling their effects is challenging, and debate continues as to the relative roles of acclimation and evolutionary adaptation in determining the observed patterns of variation. For example, although Sharoni and Halevy (2020) attribute observed seasonal variations in the elemental composition of detritus to seasonal sorting among various welladapted species, that conclusion was based on the assumption that acclimation implies a lack of nutrient limitation, which is not the assumption underlying most acclimative models, including ours. For example, the nearzero values of f_{C} in the upper 25 m during summer months (Fig. 5k, l) indicate extreme nutrient limitation, which prevents growth in the surface layers (Fig. 6b, c). In any case, only models that account for the relevant flexibilities and variations in the composition of phytoplankton can be expected to capture such feedbacks in a general yet realistic sense, which is necessary to correctly assess the relative roles of planktonrelated processes in biogeochemical cycles.
An important link between flexibility and environmental feedbacks is the role of phytoplankton in determining the elemental composition of particulate matter (Redfield, 1934). Key mechanisms involve the activities of nitrogen fixers and denitrifiers (Redfield, 1958). However, given the differences in stoichiometry of macromolecules involved in various cellular functions (Geider and La Roche, 2002), a consistent description of the acclimation of phytoplankton is necessary to represent realistically the variabilities in elemental composition of particulate matter and hence export fluxes. Fixedstoichiometry models erroneously predict constant elemental composition of detrital matter production, as demonstrated by our FS variant in this study. The socalled Droop models have been shown to capture the observed seasonal increase in detrital C:N ratios during summer, reflecting nutrient limitation of phytoplankton (e.g., Mongin et al., 2003). Representing the growth and uptake terms consistently using the RAM framework, the DA variant resolves the seasonal and vertical variations in the elemental composition of particulate matter (Fig. 9). With some exceptions, estimates of the IA variant are nearly identical to those of the DA variant, thereby implying that a more realistic representation of these can be achieved at no additional computational cost compared to a fixedstoichiometry models.
4.4 Present implementation, challenges and perspectives
Moving a coupled hydrodynamic–biogeochemical models from a 0D setup to a spatially explicit setup can be error prone and time consuming. FABM provides an easytouse coupling layer that connects a hydrodynamic model with multiple biogeochemical submodels. FABM specifies how the these models communicate by separating the hydrodynamics and biogeochemical models, with FABM acting as a glue layer in between. The biogeochemical model in this framework operates locally in space where the local source and sink terms are computed based on the local state and environment, making it feasible to scale up from 0D to nD and swap different hydrodynamic models. FABM also provides mechanisms to pass other environmental data, such as temperature, salinity and pH, from different submodules, as long as the biogeochemical models register any dependencies during initialization. Therefore, complex description of the biogeochemical models can be partitioned into several submodules. The modular implementation of our model in FABM, specifically the isolation of the phytoplankton module (Fig. 1), is expected to facilitate studies with multiple phytoplankton types. For example, without changing the model code or recompiling but just through changing a configuration file, it is possible to include further types (see Bruggeman and Bolding, 2014), which can be parameterized, e.g., according to cell size (as in, e.g., Smith et al., 2016; Dutkiewicz et al., 2020). Moreover, the isolated phytoplankton module can be relatively easily coupled with or incorporated into existing models, especially those implemented in FABM.
Currently, the model simplistically accounts for the grazing losses to higher trophic levels with a quadratic mortality term (Table 1), without describing explicitly the dynamics of predators. This limitation may prohibit realistic applications to highly productive ecosystems, where the strength of topdown control exhibits strong seasonality (e.g., Maar et al., 2014; Sailley et al., 2015). However, this problem can be easily resolved by adapting an existing zooplankton module available for FABM, such as the Nonly resolving module in the nutrient–phytoplankton–zooplankton–detritus (NPZD) example provided in the standard FABM library (Bruggeman and Bolding, 2014). An explicit consideration of zooplankton can expected to introduce additional complexities: depending on how zooplankton C and N colimitation is described, variabilities in phytoplankton stoichiometry may affect zooplankton growth (e.g., Mitra et al., 2007; Branco et al., 2018; Kerimoglu et al., 2018), and in turn, depending on the parameterization of zooplankton excretion and remineralization processes, subsequent phytoplankton blooms may occur. While it was our explicit aim to avoid such complicated indirect effects and focus on the direct effects of acclimation mechanisms on phytoplankton growth in this study, coupling the presented model to a larger ecosystem model including herbivores and their predators would allow investigating the propagation of these effects throughout the food web in a costeffective manner.
For simplicity, we have traced only N here fully (e.g., no explicit DIC but only DIN; see Eq. 4) and the model is therefore conservative with respect to N but not with respect to C. When multiple nutrient elements in the dissolved inorganic material pool (e.g., C, N, and P) are resolved, maintaining mass balance becomes more complicated under the IA assumption (see Smith et al., 2016; Ward, 2017). A FABM implementation of a carbonbased version of the model that resolves the C and N cycles is being currently developed, which we are planning to present in a separate study. The extended model will be able to resolve C, N, P and micronutrient cycles based on a common massbalance formalism and therefore allow us to investigate the validity of assuming instantaneous optimization of $\mathrm{C}:\mathrm{N}:\mathrm{P}:\mathrm{micronutrient}$ ratios under various environmental conditions (relevantly, see Bonachela et al., 2013). However, for various ecological applications, especially those resolving multiple phytoplankton types, tracing only one nutrient element, as in the current study, may be sufficient and more convenient.
In the current study, we focused on the differences between the fully acclimative IA and DA variants, and an entirely nonacclimative variant. Our acclimation scheme consists of four acclimative flexibilities: variability of internal nutrient quota, optimization of uptake vs. growth tradeoff, optimization of maximum uptake vs. affinity tradeoff and optimization of chlorophyll density in chloroplast density (and as an additional half step, size of the chloroplast; see Appendix B). In a future study, we are planning to investigate the relative importance of each of these flexibilities for the organismal fitness under various environmental conditions: such an assessment would not only help the model developers to prioritize the research needs but may also provide insights into the evolution of these acclimative flexibilities.
In this study, we present a FABM implementation of the “NflexPD” model and the behavior of three variants it can emulate: a fixedstoichiometry (FS) variant that lacks any acclimative flexibility and explicitly tracks only N bound to phytoplankton; a dynamic acclimation (DA) variant that resolves various acclimative flexibilities by explicitly tracking the C and N in phytoplankton; and the instantaneous acclimation (IA) variant that resolves the same flexibilities as the DA variant but by tracking the N in phytoplankton as in the FS variant.
By applying the NflexPD model coupled to an idealized, 1D water column model, we aimed to understand (i) whether and how the behavior of the IA and DA variants differ; and (ii) whether and how the behavior of the acclimative variants differ from the nonacclimative, fixedstoichiometry variant. With regard to the first of our objectives, we found that behavior of IA is stable and in many respects very similar to that of DA, although differences arise during the spring and autumn transitions, due to the lagged response and vertical transport of nutrient quotas in the DA variant. With this, our study provides a proof of concept that the IA approach is applicable in spatially explicit setups and hints at conditions under which deviations from the fully explicit variant can be expected. With regard to the second objective, we found substantial differences between the behavior of the FS and acclimative variants: with the particular parameterization we showcased here, the acclimative variants estimated smaller spring blooms but sustained growth during summer and stronger nitrogen depletion in the surface layers, as well as steeper chlorophyll layers at the thermocline; and unlike the FS variant, they can reproduce the variabilities in C:N of particulate matter. Moreover, a subset of quantities estimated by the FS variant, such as the phytoplankton biomass and NPP rates, was found to be strongly sensitive to the prescribed parameters such as Q, which in this study was derived as a spatiotemporal average from the IA variant but is typically an adjustable parameter, thus implying a higher degree of freedom. These qualitative differences provide insights on the impact of acclimative flexibilities on model response and their ecosystemscale implications. The model implementation presented here tracks only N as dissolved nutrient, which restricts its utility in biogeochemical studies that require a complete representation of the cycling of multiple elements, but it can be readily used in various ecological contexts.
A1 R_{N} for FS variant
According to Eq. (7), ${R}_{\mathrm{N}}={\mathit{\zeta}}_{\mathrm{N}}\cdot V$. For the DA and IA variants, V can be calculated externally (Eq. 12); hence, it can be R_{N}. For the FS variant on the other hand, there is no explicit solution for V, but it can only be calculated as a function of μ, ($V=\mathit{\mu}\cdot Q$, Eq. 6), and since μ in turn depends on R_{N} ($\mathit{\mu}={\mathit{\mu}}_{\mathrm{net}}{R}_{\mathrm{N}}$, Eq. 7), R_{N} cannot be directly calculated. Expanding the terms in Eq. (7) according to Eqs. (6), (13) and (20) yields
Reorganizing yields
Substituting this with μ in
we obtain a Vindependent expression for R_{N}:
It can be verified that, when this term is substituted in $\mathit{\mu}={\mathit{\mu}}_{\mathrm{net}}{R}_{\mathrm{N}}$, it yields $\mathit{\mu}={\mathit{\mu}}_{\mathrm{net}}{\mathit{\zeta}}_{\mathrm{N}}\cdot \mathit{\mu}\cdot Q={\mathit{\mu}}_{\mathrm{net}}{\mathit{\zeta}}_{\mathrm{N}}\cdot V$, i.e., Eq. (7), implying that using R_{N} in Eq. (A4) for the FS variant makes Eq. (7) valid for the FS variant as well.
A2 Optimal Q and f_{V}
In Eq. (7), substituting μ_{g}, R_{N} and R_{Chl} with the expanded forms in Eqs. (13), (20) and (25), respectively, and subsequently expanding θ, using Eq. (24) yields
Reorganizing yields
Substituting the term in square brackets with ${\widehat{\mathit{\mu}}}_{\mathrm{net}}$ based on Eq. (7) and expanding f_{C} using Eq. (11) yields
At this point, it can be readily recognized that Eq. (A7) is equivalent to Eq. (5) in Pahlow and Oschlies (2013), with the only difference being their denotation of ${\widehat{\mathit{\mu}}}_{\mathrm{net}}$ as ${\widehat{\mathit{\mu}}}^{I}$. Note that their formulation of respiration losses within the chloroplast as a fraction of gross growth with respect to chloroplast (i.e., ${\widehat{\mathit{\mu}}}^{I}={\widehat{\mathit{\mu}}}_{g}^{I}(\mathrm{1}{\mathit{\zeta}}^{C})$ in their notation) differs from the more precise formulation we used here, which considers a base loss rate independent of gross growth. However, considering that ${\widehat{\mathit{\mu}}}_{\mathrm{net}}$ (just like their ${\widehat{\mathit{\mu}}}^{I}$) is independent of Q and f_{V}) the solutions provided by Pahlow and Oschlies (2013) for ${f}_{V}^{o}$ (i.e., their Eq. 9, our Eq. 14) and Q (their Eq. 10, our Eq. 10) can be directly used only after replacing ${\widehat{\mathit{\mu}}}^{I}$ in the original solutions with ${\widehat{\mathit{\mu}}}_{\mathrm{net}}$ for the latter.
Given the similar roles of f_{C} in the IA and DA variants and the nutrient limitation term, L_{N}, in the FS variant for calculating μ_{g} (see Sect. 2.2.2), L_{N} can be considered as a proxy for the relative size of the chloroplast. Therefore, f_{C} in Eqs. (24) and (25) can be replaced by L_{N} for scaling the chloroplastspecific chlorophyll density and respiration costs in order to represent spatiotemporal variations of the cellular Chl:C ratio and proportional respiration costs.
When this is done, unlike the original results shown in the main text (Fig. 4m), a spatiotemporally variable Chl:C ratio (Fig. B1c) is obtained. Monotonically increasing L_{N} with depth during summer (Fig. 5j) reduces Chl at the surface and enhances it at the deeper layers relative to the Chl pattern obtained with constant Chl:C (compare Fig. 4m vs. Fig. B1a). However, due to the missing unimodal signal through $\widehat{\mathit{\theta}}$ as accounted for by the IA and DA variants (see Fig. 5b, c), the resulting Chl pattern is still qualitatively different from those estimated by the truly acclimative variants (compare Fig. B1a vs. Fig. 4k, l). Furthermore, the relatively higher value of L_{N} during the spring bloom under nutrientrich conditions (Fig. 5j) relative to the prescribed, constant value of f_{C}=0.44 used for the case with constant chloroplast size (hence, constant Chl:C) shown in the main text as yielded by the prescribed values of f_{V}, Q and Q_{0} (Table 3 and Eq. 11), results in greater R_{Chl} (compare Fig. B1d vs. Fig. 6j). Hence, net cellular growth rate, μ, becomes slightly lower than in the constant chloroplast case during the spring bloom (compare Fig. B1b vs. Fig. 6a). On the other hand, during summer, relatively lower values of L_{N} make R_{Chl} lower and μ greater compared to the constant chloroplast case.
The dynamics of the Phy_{C} within the top 50 m as simulated with this flavor of the FS variant with variable chloroplast size are almost identical to those simulated by the standard, “vanilla” version with constant chloroplast size (compare Fig. B2d with 7d). Relatively higher R_{Chl} at nutrientrich conditions during winter and early spring makes the winter Phy_{C} concentrations (Fig. B2d) lower in comparison to the standard case (Fig. 7d). On the other hand, relatively lower R_{Chl} at nutrientscarce summer conditions makes the Phy_{C} concentrations (Fig. B2d) slightly higher than the standard case (Fig. 7d). As a result, the average DIN concentrations in the surface layers at 50 m become slightly lower than the standard case (Fig. B2a vs. Fig. 7a), which is better observed in lower L_{N} (Fig. B2c vs. Fig. 7c) due to the strong response of the function at low concentrations.
Despite the differences in details explained above, especially based on the preserved qualitative differences in simulated Phy_{C} concentrations between the FS and acclimative variants, it can be concluded that the overall conclusions are insensitive to the assumption regarding the size of the chloroplast of the FS variant.
For running the model and reproducing the results presented in this study, FABM and GOTM need to be downloaded and installed. See https://github.com/fabmmodel/fabm/wiki/GOTM (last access: 4 October 2021). for the instructions. The version of the FABMNflexPD used in this paper has been stored in a Zenodo repository at https://doi.org/10.5281/zenodo.4761937 (Kerimoglu, 2021). Instructions for compiling FABMNflexPD for GOTMFABM and a 0D setup are provided in README.md. The “src” folder contains the Fortran codes. The model was implemented as two separate modules: the “phy.F90” module that describes phytoplankton growth and the “abio.F90” module that describes everything other than phytoplankton (Fig. 1). The phytoplankton module can reproduce the behavior of all three different variants considered in the paper through optional parameters. The “testcases” folder contains the configuration (yaml) file that was used to produce the results presented in this paper, thereby providing examples of how each variant can be initiated.
The underlying research data can be accessed from https://doi.org/10.17882/83699 (Kerimoglu et al., 2021).
OK and SLS designed the study; OK implemented the model code, performed the model runs and prepared the figures; PA tested the model; all authors contributed to the writing of the manuscript.
The authors declare that they have no conflict of interest.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We are grateful to Markus Pahlow for critically reviewing the model description, helping to correct the implementation of IA and DA variants and providing the code for the model diagram. We acknowledge the developers of the opensource software used in this study, foremost FABM and GOTM. We are grateful for the comments of two anonymous reviewers, which helped significantly to improve the manuscript.
This research has been supported by the Deutsche Forschungsgemeinschaft (grant no. KE 1970/21, PI: OK), the Japan Society for the Promotion of Science, JSPS (PI: SLS), and the Umweltbundesamt (grant no. 3718252110, PI: OK).
This paper was edited by Robert Marsh and reviewed by two anonymous referees.
Aksnes, D. L. and Egge, J.: A theoretical model for nutrient uptake in phytoplankton, Mar. Ecol. Prog. Ser., 70, 65–72, 1991. a
Anderson, T. R. and Pondaven, P.: Nonredfield carbon and nitrogen cycling in the Sargasso Sea: Pelagic imbalances and export flux, DeepSea Res. Pt. I, 50, 573–591, https://doi.org/10.1016/S09670637(03)000347, 2003. a
Anugerahanti, P., Kerimoglu, O., and Smith, S. L.: Enhancing Ocean Biogeochemical Models With Phytoplankton Variable Composition, Front. Mar. Sci., 8, 675428, https://doi.org/10.3389/fmars.2021.675428, 2021. a, b
Armstrong, R. A.: An optimizationbased model of iron–light–ammonium colimitation of nitrate uptake and phytoplankton growth, Limnol. Oceanogr., 44, 1436–1446, https://doi.org/10.4319/lo.1999.44.6.1436, 1999. a
Ayata, S. D., Lévy, M., Aumont, O., Sciandra, A., SainteMarie, J., Tagliabue, A., and Bernard, O.: Phytoplankton growth formulation in marine ecosystem models: Should we take into account photoacclimation and variable stoichiometry in oligotrophic areas?, J. Marine Syst., 125, 29–40, https://doi.org/10.1016/j.jmarsys.2012.12.010, 2013. a, b
Behrenfeld, M. J., Westberry, T. K., Boss, E. S., O'Malley, R. T., Siegel, D. A., Wiggert, J. D., Franz, B. A., McClain, C. R., Feldman, G. C., Doney, S. C., Moore, J. K., Dall'Olmo, G., Milligan, A. J., Lima, I., and Mahowald, N.: Satellitedetected fluorescence reveals global physiology of ocean phytoplankton, Biogeosciences, 6, 779–794, https://doi.org/10.5194/bg67792009, 2009. a
Bergeron, M. and Tremblay, J.E.: Shifts in biological productivity inferred from nutrient drawdown in the southern Beaufort Sea (2003–2011) and northern BaffinBay (1997–2011), Canadian Arctic, Geophys. Res. Lett., 41, 3979–3987, https://doi.org/10.1002/2014GL059649.Received, 2014. a
Bonachela, J. A., Allison, S. D., Martiny, A. C., and Levin, S. A.: A model for variable phytoplankton stoichiometry based on cell protein regulation, Biogeosciences, 10, 4341–4356, https://doi.org/10.5194/bg1043412013, 2013. a
Branco, P., Egas, M., Elser, J. J., and Huisman, J.: EcoEvolutionary Dynamics of Ecological Stoichiometry in Plankton Communities, Am. Nat., 192, E000–E000, https://doi.org/10.1086/697472, 2018. a
Bruggeman, J. and Bolding, K.: A general framework for aquatic biogeochemical models, Environ. Modell. Softw., 61, 249–265, https://doi.org/10.1016/j.envsoft.2014.04.002, 2014. a, b, c, d
Bruggeman, J. and Kooijman, S. A. L. M.: A biodiversityinspired approach to aquatic ecosystem modeling, Limonol. Oceanogr., 52, 1533–1544, https://doi.org/10.4319/lo.2007.52.4.1533, 2007. a
Burchard, H., Bolding, K., Kühn, W., Meister, A., Neumann, T., and Umlauf, L.: Description of a flexible and extendable physicalbiogeochemical model system for the water column, J. Marine Syst., 61, 180–211, https://doi.org/10.1016/j.jmarsys.2005.04.011, 2006. a, b, c
Burmaster, D. E.: The Continuous Culture of Phytoplankton: Mathematical Equivalence Among Three SteadyState Models, Am. Nat., 113, 123–134, https://doi.org/10.1086/283368, 1979. a, b, c, d
Burson, A., Stomp, M., Akil, L., Brussaard, C. P. D., and Huisman, J.: Unbalanced reduction of nutrient loads has created an offshore gradient from phosphorus to nitrogen limitation in the North Sea, Limnol. Oceanogr., 61, 869–888, https://doi.org/10.1002/lno.10257, 2016. a
Button, D. K.: On the theory of control of microbial growth kinetics by limiting nutrient concentrations, DeepSea Res., 25, 1163–1177, https://doi.org/10.1016/01466291(78)900115, 1978. a, b, c
Caperon, J.: Population growth response of Isochrysis Galbana to nitrate variation at limiting concentrations, Ecology, 49, 866–872, 1968. a, b
Chakraborty, S., Cadier, M., Visser, A. W., Bruggeman, J., and Andersen, K. H.: Latitudinal Variation in Plankton Traits and Ecosystem Function, Global Biogeochem. Cy., 32, e2020GB006564, https://doi.org/10.1029/2020GB006564, 2020. a
Chen, B. and Smith, S. L.: Optimalitybased approach for computationally efficient modeling of phytoplankton growth, chlorophylltocarbon, and nitrogentocarbon ratios, Ecol. Model., 385, 197–212, https://doi.org/10.1016/j.ecolmodel.2018.08.001, 2018. a
Christian, J. R.: Biogeochemical cycling in the oligotrophic ocean: Redfield and nonRedfield models, Limnol. Oceanogr., 50, 646–657, https://doi.org/10.4319/lo.2005.50.2.0646, 2005. a
Cloern, J., Grenz, C., and VidergarLucas, L.: An empirical model of the phytoplankton chlorophyll : carbon ratio – the conversion factor between productivity and growth rate, Limnol. Oceanogr., 40, 1313–1321, https://doi.org/10.4319/lo.1995.40.7.1313, 1995. a
Doney, S. C., Glover, D. M., and Najjar, R. G.: A new coupled, onedimensional biologicalphysical model for the upper ocean: Applications to the JGOFS Bermuda Atlantic Timeseries study (BATS) site, DeepSea Res. Pt. II, 43, 591–624, https://doi.org/10.1016/09670645(95)001042, 1996. a
Droop, M.: Vitamin B12 and marine ecology. IV. The kinetics of uptake, growth and inhibition in Monochrysis lutheri, Journal of the Marine Biological Association of the United Kingdom, 48, 689–733, 1968. a, b, c
Dugdale, R.: Nutrient Limitation in the Sea: Dynamics, Identification and Significance, Limnol. Oceanogr., 12, 685–695, https://doi.org/10.4319/lo.1967.12.4.0685, 1967. a
Dutkiewicz, S., Cermeno, P., Jahn, O., Follows, M. J., Hickman, A. E., Taniguchi, D. A. A., and Ward, B. A.: Dimensions of marine phytoplankton diversity, Biogeosciences, 17, 609–634, https://doi.org/10.5194/bg176092020, 2020. a, b, c
Edelaar, P. and Bolnick, D. I.: Appreciating the Multiple Processes Increasing Individual or Population Fitness, Trends Ecol. Evol., 34, 435–446, https://doi.org/10.1016/j.tree.2019.02.001, 2019. a
Flynn, K. J., St John, M., Raven, J. A., Skibinski, D. O., Allen, J. I., Mitra, A., and Hofmann, E. E.: Acclimation, adaptation, traits and tradeoffs in plankton functional type models: Reconciling terminology for biology and modelling, J. Plankton Res., 37, 683–691, https://doi.org/10.1093/plankt/fbv036, 2015. a, b
Follows, M. J. and Dutkiewicz, S.: Modeling Diverse Communities of Marine Microbes, Annu. Rev. Mar. Sci., 3, 427–451, https://doi.org/10.1146/annurevmarine120709142848, 2011. a
Geider, R. and La Roche, J.: Redfield revisited: variability of C:N:P in marine microalgae and its biochemical basis, Eur. J. Phycol., 37, 1–17, https://doi.org/10.1017/S0967026201003456, 2002. a, b, c
Geider, R., Maclntyre, H., and Kana, T.: A dynamic regulatory model of phytoplanktonic acclimation to light, nutrients, and temperature, Limnol. Oceanogr., 43, 679–694, https://doi.org/10.4319/lo.1998.43.4.0679, 1998. a, b, c
Gerloff, G. C. and Skoog, F.: Cell Contents of Nitrogen and Phosphorous as a Measure of Their Availability for Growth of Microcystis Aeruginosa, Ecology, 35, 348–353, 1954. a
Grover, J.: Resource competition in a variable environment: phytoplankton growing according to the variableinternalstores model, Am. Nat., 138, 811–835, 1991. a, b
Grover, J. P.: Is storage an adaptation to spatial variation in resource availability?, Am. Nat., 173, E44–E61, https://doi.org/10.1086/595751, 2009. a, b
Guidi, L., Jackson, G. A., Stemmann, L., Miquel, J. C., Picheral, M., and Gorsky, G.: Relationship between particle size distribution and flux in the mesopelagic zone, DeepSea Res. Pt. I, 55, 1364–1374, https://doi.org/10.1016/j.dsr.2008.05.014, 2008. a
Halsey, K. H. and Jones, B. M.: Phytoplankton strategies for photosynthetic energy allocation, Annu. Rev. Mar. Sci., 7, 265–297, https://doi.org/10.1146/annurevmarine010814015813, 2015. a
Johnson, K. and Goody, R.: The Original Michaelis Constant: Translation of the 1913 Michaelis–Menten Paper, Biochemistry, 50, 8264–8269, https://doi.org/10.1021/bi201284u, 2011. a
Johnson, K. S., Plant, J. N., Dunne, J. P., Talley, L. D., and Sarmiento, J. L.: Annual nitrate drawdown observed by SOCCOM profiling floats and the relationship to annual net community production, J. Geophys. Res.Oceans, 122, 6668–6683, https://doi.org/10.1002/2017JC012839, 2017. a
Kerimoglu, O.: OnurKerimoglu/fabmnflexpd: Version 1.0 release candidate 0 (v1.0rc0), Zenodo [code], https://doi.org/10.5281/zenodo.4761937, 2021. a
Kerimoglu, O., Straile, D., and Peeters, F.: Role of phytoplankton cell size on the competition for nutrients and light in incompletely mixed systems., J. Theor. Biol., 300, 330–343, https://doi.org/10.1016/j.jtbi.2012.01.044, 2012. a, b
Kerimoglu, O., Hofmeister, R., Maerz, J., Riethmüller, R., and Wirtz, K. W.: The acclimative biogeochemical model of the southern North Sea, Biogeosciences, 14, 4499–4531, https://doi.org/10.5194/bg1444992017, 2017. a, b, c
Kerimoglu, O., Große, F., Kreus, M., Lenhart, H.J., and van Beusekom, J. E.: A modelbased projection of historical state of a coastal ecosystem: relevance of phytoplankton stoichiometry, Sci. Total Environ., 639, 1311–1323, https://doi.org/10.1016/j.scitotenv.2018.05.215, 2018. a, b
Kerimoglu, O. Anugerahanti, P., and Smith, S. L.: FABMNflexPD 1.0: cyclicequilibria simulated by three model variants in an idealized highlatitude 1D water column, SEANOE [data set], https://doi.org/10.17882/83699, 2021. a
Klausmeier, C., Litchman, E., Daufresne, T., and Levin, S.: Optimal nitrogentophosphorus stoichiometry of phytoplankton, Nature, 429, 171–174, https://doi.org/10.1029/2001GL014649, 2004. a, b
Kruskopf, M. and Flynn, K. J.: Chlorophyll content and fluorescence responses cannot be used to gauge reliably phytoplankton biomass, nutrient status or growth rate, New Phytol., 169, 525–536, https://doi.org/10.1111/j.14698137.2005.01601.x, 2006. a
Kwiatkowski, L., Aumont, O., Bopp, L., and Ciais, P.: The Impact of Variable Phytoplankton Stoichiometry on Projections of Primary Production, Food Quality, and Carbon Uptake in the Global Ocean, Global Biogeochem. Cy., 32, 516–528, https://doi.org/10.1002/2017GB005799, 2018. a, b
Laufkötter, C., Vogt, M., Gruber, N., AitaNoguchi, M., Aumont, O., Bopp, L., Buitenhuis, E., Doney, S. C., Dunne, J., Hashioka, T., Hauck, J., Hirata, T., John, J., Le Quéré, C., Lima, I. D., Nakano, H., Seferian, R., Totterdell, I., Vichi, M., and Völker, C.: Drivers and uncertainties of future global marine primary production in marine ecosystem models, Biogeosciences, 12, 6955–6984, https://doi.org/10.5194/bg1269552015, 2015. a, b
Laws, E. A. and Chalup, M. S.: A microalgal growth model, Limnol. Oceanogr., 35, 597–608, https://doi.org/10.4319/lo.1990.35.3.0597, 1990. a
Litchman, E., Klausmeier, C. A., and Yoshiyama, K.: Contrasting size evolution in marine and freshwater diatoms, P. Natl. Acad. Sci. USA, 106, 2665–2670, 2009. a
Maar, M., Rindorf, A., Møller, E. F., Christensen, A., Madsen, K. S., and van Deurs, M.: Zooplankton mortality in 3D ecosystem modelling considering variable spatial–temporal fish consumptions in the North Sea, Prog. Oceanogr., 124, 78–91, https://doi.org/10.1016/j.pocean.2014.03.002, 2014. a
Martiny, A. A. C., Pham, C. C. T. A., Primeau, F. F. W., Vrugt, J. A., Moore, J. K., Levin, S. A., and Lomas, M. W.: Strong latitudinal patterns in the elemental ratios of marine plankton and organic matter, Nat. Geosci., 6, 279–283, https://doi.org/10.1038/ngeo1757, 2013. a
Mitra, A., Flynn, K. J., and Fasham, M. J. R.: Accounting for grazing dynamics in nitrogenphytoplanktonzooplankton models, Limnol. Oceanogr., 52, 649–661, https://doi.org/10.4319/lo.2007.52.2.0649, 2007. a
Mongin, M., Nelson, D. M., Pondaven, P., Brzezinski, M. A., and Tréguer, P.: Simulation of upperocean biogeochemistry with a flexiblecomposition phytoplankton model: C, N and Si cycling in the western Sargasso Sea, DeepSea Res. Pt. I, 50, 1445–1480, https://doi.org/10.1016/j.dsr.2003.08.003, 2003. a, b
Monod, J.: The growth of bacterial cultures, Annu. Rev. Microbiol., 3, 371–394, 1949. a, b
Moore, J. K., Doney, S. C., Kleypas, J. A., Glover, D. M., and Fung, I. Y.: An intermediate complexity marine ecosystem model for the global domain, DeepSea Res. Pt. II, 49, 403–462, https://doi.org/10.1016/S09670645(01)001084, 2002. a
Moore, J. K., Doney, S. C., and Lindsay, K.: Upper ocean ecosystem dynamics and iron cycling in a global threedimensional model, Global Biogeochem. Cy., 18, GB4028, https://doi.org/10.1029/2004GB002220, 2004. a
Moreno, A. R. and Martiny, A. C.: Ecological Stoichiometry of Ocean Plankton, Annu. Rev. Mar. Sci., 10, 43–69, https://doi.org/10.1146/annurevmarine121916063126, 2018. a, b
Oschlies, A. and Schartau, A.: Basinscale performance of a locally optimized marine ecosystem model, Research, J. Marine, 63, 335–358, https://doi.org/10.1016/b9780408707008.500185, 2005. a, b
Pahlow, M.: Linking chlorophyll – nutrient dynamics to the Redfield N : C ratio with a model of optimal phytoplankton growth, Mar. Ecol. Prog. Ser., 287, 33–43, https://doi.org/10.3354/meps287033, 2005. a, b, c, d
Pahlow, M. and Oschlies, A.: Chain model of phytoplankton P, N and light colimitation, Mar. Ecol. Prog. Ser., 376, 69–83, https://doi.org/10.3354/meps07748, 2009. a
Pahlow, M. and Oschlies, A.: Optimal allocation backs Droop's cellquota model, Mar. Ecol. Prog. Ser., 473, 1–5, https://doi.org/10.3354/meps10181, 2013. a, b, c, d, e, f, g, h
Pahlow, M., Dietze, H., and Oschlies, A.: Optimalitybased model of phytoplankton growth and diazotrophy, Mar. Ecol. Prog. Ser., 489, 1–16, https://doi.org/10.3354/meps10449, 2013. a, b, c, d, e, f, g, h, i, j
Pahlow, M., Chien, C.T., Arteaga, L. A., and Oschlies, A.: Optimalitybased nonRedfield plankton–ecosystem model (OPEM v1.1) in UVicESCM 2.9 – Part 1: Implementation and model behaviour, Geosci. Model Dev., 13, 4663–4690, https://doi.org/10.5194/gmd1346632020, 2020. a, b
Paulson, C. and Simpson, J.: Irradience measurements in the upper ocean, J. Phys. Oceanogr., 7, 952–956, 1977. a
Platt, T. and Jassby, A. D.: the Relationship Between Photosynthesis and Light for Natural Assemblages of Coastal Marine Phytoplankton, J. Phycol., 12, 421–430, https://doi.org/10.1111/j.15298817.1976.tb02866.x, 1976. a
Redfield, A.: On the proportions of organic derivatives in sea water and their relation to the composition of plankton, in: James Johnstone Memorial Volume, edited by: Daniel, R., University Press of Liverpool, 177–192, 1934. a
Redfield, A.: The biological control of chemical factors in the environment, Am. Sci., 46, 205–221, 1958. a
Sailley, S. F., Polimene, L., Mitra, A., Atkinson, A., and Allen, J. I.: Impact of zooplankton food selectivity on plankton dynamics and nutrient cycling, J. Plankton Res., 37, 519–529, https://doi.org/10.1093/plankt/fbv020, 2015. a
SchourupKristensen, V., Sidorenko, D., WolfGladrow, D. A., and Völker, C.: A skill assessment of the biogeochemical model REcoM2 coupled to the Finite Element Sea Ice–Ocean Model (FESOM 1.3), Geosci. Model Dev., 7, 2769–2802, https://doi.org/10.5194/gmd727692014, 2014. a
Sharoni, S. and Halevy, I.: Nutrient ratios in marine particulate organic matter, Sci. Adv., 6, 1–10, https://doi.org/10.1126/sciadv.aaw9371, 2020. a
Shuter, B.: A model of physiological adaptation in unicellular algae, J. Theor. Biol., 78, 519–552, 1979. a
Smith, S. L., Yamanaka, Y., Pahlow, M., and Oschlies, A.: Optimal uptake kinetics: physiological acclimation explains the pattern of nitrate uptake by phytoplankton in the ocean, Mar. Ecol. Prog. Ser., 384, 1–12, https://doi.org/10.3354/meps08022, 2009. a, b
Smith, S. L., Pahlow, M., Merico, A., and Wirtz, K. W.: Optimalitybased modeling of planktonic organisms, Limnol. Oceanogr., 56, 2080–2094, https://doi.org/10.4319/lo.2011.56.6.2080, 2011. a, b, c
Smith, S. L., Pahlow, M., Merico, A., AcevedoTrejos, E., Sasai, Y., Yoshikawa, C., Sasaoka, K., Fujiki, T., Matsumoto, K., and Honda, M. C.: Flexible phytoplankton functional type (FlexPFT) model: sizescaling of traits and optimal growth, J. Plankton Res., 38, 977–992, https://doi.org/10.1093/plankt/fbv038, 2016. a, b, c, d, e, f, g, h, i, j
Ward, B. A.: Assessing an efficient “Instant Acclimation” approximation of dynamic phytoplankton stoichiometry, J. Plankton Res., 39, 803–814, https://doi.org/10.1093/plankt/fbx040, 2017. a, b, c
Wirtz, K. W. and Kerimoglu, O.: Autotrophic Stoichiometry Emerging from Optimality and Variable Colimitation, Front. Ecol. Evol., 4, 131, https://doi.org/10.3389/fevo.2016.00131, 2016. a, b, c
 Abstract
 Introduction
 Model description
 Results
 Discussion
 Conclusions
 Appendix A: Details of derivations
 Appendix B: FS variant with a variable chloroplast size
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Model description
 Results
 Discussion
 Conclusions
 Appendix A: Details of derivations
 Appendix B: FS variant with a variable chloroplast size
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References