FABM-NflexPD 1.0: assessing an instantaneous acclimation approach for modeling phytoplankton growth

Coupled physical–biogeochemical models can generally reproduce large-scale 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 1-D 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 1-D 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, non-acclimative and fixed-stoichiometry (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 cost-effective improvement over a fixed-stoichiometry 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 Introduction

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 so-called “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 satellite-based 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 so-called “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 photo-acclimation, i.e., variation of Chl:C, using one more state variable for Chl bound to phytoplankton (Moore et al., 2002; Schourup-Kristensen 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 photo-acclimation 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 optimality-based 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 trade-offs (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⁻¹]) of phytoplankton.

f_V:

fractional allocation (dimensionless) to the nutrient uptake compartment (protoplast) to optimize the trade-off between photosynthesis (μ) and nutrient uptake (V), as described by Pahlow and Oschlies (2013).

f_A:

fractional allocation (dimensionless) to affinity to optimize the trade-off 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⁻¹]) to optimize the trade-off 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 above-mentioned 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 short-term (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 short-term experiments, where the response of phytoplankton may differ substantially in terms of nutrient uptake vs. carbon-based 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 short-term 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 0-D (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 fixed-stoichiometry model, in a 0-D 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 fixed-stoichiometry 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 nutrient-replete 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, non-acclimative fixed-stoichiometry (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.

Figure 1 Diagram of the FABM-NflexPD model. Abiotic components, dissolved organic matter (DOM), dissolved inorganic matter (DIM) and Det are calculated by the module abio.F90, which are then coupled to the phytoplankton simulated by the module phy.F90 that simulates the dynamics of Phy_N, Phy_C and Phy_Chl by the DA, IA and FS variants (see Sect. 2.2.2). Solid circles in the phytoplankton module represent state variables, dashed circles/ellipsoids represent diagnostically calculated variables, and solid squares (for FS) represent prescribed values. The DA variant estimates the N, C and Chl content of phytoplankton based on a resource allocation scheme, whereas the FS variant estimates only N prognostically, while C and Chl are based on prescribed values of nitrogen quota (Q) and cellular Chl:C ratio (θ) (see the text).

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 Model description

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.

$$\begin{array}{}\text{(1a)}& {\displaystyle }s\left({\mathrm{Phy}}_{\mathrm{N}}\right)& {\displaystyle }=F_{\mathrm{DIN}-{\mathrm{Phy}}_{\mathrm{N}}}-F_{{\mathrm{Phy}}_{\mathrm{N}}-{\mathrm{Det}}_{\mathrm{N}}}\text{(1b)}& {\displaystyle }s\left({\mathrm{Phy}}_{\mathrm{C}}\right)& {\displaystyle }=\underset{\text{Uptake}}{\underbrace{F_{\mathrm{DIC}-{\mathrm{Phy}}_{\mathrm{C}}}}}-\underset{\text{Mortality}}{\underbrace{F_{{\mathrm{Phy}}_{\mathrm{C}}-{\mathrm{Det}}_{\mathrm{C}}}}}\mathit{\left\{}\mathrm{DA}\mathit{\right\}}\end{array}$$

$$\begin{array}{}\text{(2a)}& {\displaystyle }s\left({\mathrm{Det}}_{\mathrm{N}}\right)& {\displaystyle }=F_{{\mathrm{Phy}}_{\mathrm{N}}-{\mathrm{Det}}_{\mathrm{N}}}-F_{{\mathrm{Det}}_{\mathrm{N}}-\mathrm{DON}}\text{(2b)}& {\displaystyle }s\left({\mathrm{Det}}_{\mathrm{C}}\right)& {\displaystyle }=\underset{\mathrm{Mortality}}{\underbrace{F_{{\mathrm{Phy}}_{\mathrm{C}}-{\mathrm{Det}}_{\mathrm{C}}}}}-\underset{\mathrm{Hydrolysis}}{\underbrace{F_{{\mathrm{Det}}_{\mathrm{C}}-\mathrm{DOC}}}}\end{array}$$

$$\begin{array}{}\text{(3a)}& {\displaystyle }s\left(\mathrm{DON}\right)& {\displaystyle }=F_{{\mathrm{Det}}_{\mathrm{N}}-\mathrm{DON}}-F_{\mathrm{DON}-\mathrm{DIN}}\text{(3b)}& {\displaystyle }s\left(\mathrm{DOC}\right)& {\displaystyle }=\underset{\mathrm{Hydrolysis}}{\underbrace{F_{{\mathrm{Det}}_{\mathrm{C}}-\mathrm{DOC}}}}-\underset{\mathrm{Remineralization}}{\underbrace{F_{\mathrm{DOC}-\mathrm{DIC}}}}\end{array}$$

$$\begin{array}{}\text{(4)}& s\left(\mathrm{DIN}\right)=\underset{\mathrm{Remineralization}}{\underbrace{F_{\mathrm{DON}-\mathrm{DIN}}}}-\underset{\mathrm{Nuptake}}{\underbrace{F_{\mathrm{DIN}-{\mathrm{Phy}}_{\mathrm{N}}}}}\end{array}$$

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:

$$\begin{array}{}\text{(5)}& F_{\mathrm{DIN}-{\mathrm{Phy}}_{\mathrm{N}}}=V\cdot {\mathrm{Phy}}_{\mathrm{C}}.\end{array}$$

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:

$$\begin{array}{}\text{(6)}& V=\mathit{\mu }\cdot Q\mathit{\{}\mathrm{FS},\mathrm{IA}\mathit{\}},\end{array}$$

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):

$$\begin{array}{}\text{(7)}& \mathit{\mu }={\mathit{\mu }}_{\mathrm{g}}-R_{\mathrm{Chl}}-R_{\mathrm{N}}={\mathit{\mu }}_{\mathrm{net}}-{\mathit{\zeta }}_{\mathrm{N}}\cdot V,\end{array}$$

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):

$$\begin{array}{}\text{(8)}& F_{\mathrm{DIC}-{\mathrm{Phy}}_{\mathrm{C}}}=\mathit{\mu }\cdot {\mathrm{Phy}}_{\mathrm{C}}\mathit{\left\{}\mathrm{DA}\mathit{\right\}}.\end{array}$$

Table 1 Definitions, expansions/values and units of terms/symbols regarding the fluxes between model compartments.

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 balanced-growth 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 chloroplast-specific net growth and protoplast-specific N uptake rates (Table 2), and Q₀ 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 fixed-stoichiometry variant, Q is a prescribed parameter (Table 2).

Table 2 Summary of differences between model variants. n/a denotes “not applicable”. ^* Prescribed; see Table 3.

2.2.2 Flexibility II: growth vs. nutrient uptake

Given the high nitrogen content in the enzymes responsible for both CO₂ fixation and nutrient uptake and assimilation (Geider and La Roche, 2002), we consider a trade-off in the allocation of nitrogen between carbon fixation and nutrient uptake for the acclimative variants, whereas this trade-off is ignored for the FS variant.

• For the acclimative variants (DA and IA), following Pahlow and Oschlies (2013), the trade-off 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 protoplast-specific 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 fixed-stoichiometry variant, the gross growth rate, μ_g is obtained by the multiplication of ${\widehat{\mathit{\mu }}}_{\mathrm{g}}$, for FS, interpreted as the light-limited 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 fixed-stoichiometry 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 protoplast-specific 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 trade-off between affinity vs. maximum uptake rate. This trade-off 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 fixed-stoichiometry variant ignores this trade-off 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 biomass-weighted spatiotemporal average K_N value so obtained was prescribed for the FS variant (Table 3).

2.2.4 Flexibility IV: photo-acclimation

Photo-acclimation 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 chloroplast-specific synthesis and maintenance costs of chlorophyll, from the gross growth rate within the chloroplast, i.e.,

$$\begin{array}{}\text{(20)}& {\displaystyle }{\widehat{\mathit{\mu }}}_{\mathrm{net}}{\displaystyle }={\widehat{\mathit{\mu }}}_{\mathrm{g}}-{\widehat{R}}_{\mathrm{Chl}},\end{array}$$

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 light-saturation of the photosynthetic apparatus, L_I:

$$\begin{array}{}\text{(21)}& {\displaystyle }{\widehat{\mathit{\mu }}}_{\mathrm{g}}{\displaystyle }=L_{\mathrm{D}}{\widehat{\mathit{\mu }}}_{\mathrm{0}}{L}_I.\end{array}$$

L_I is a saturating function of daytime average light, $\overline{I}$, and chlorophyll density in chloroplasts, $\widehat{\mathit{\theta }}$:

$$\begin{array}{}\text{(22)}& {\displaystyle }{L}_I{\displaystyle }=\mathrm{1}-\exp \left({\displaystyle \frac{-\mathit{\alpha }\widehat{\mathit{\theta }}\overline{I}}{{\widehat{\mathit{\mu }}}_{\mathrm{0}}}}\right),\end{array}$$

where α is light affinity. Returning to Eq. (20), ${\widehat{R}}_{\mathrm{Chl}}$ is given by

$$\begin{array}{}\text{(23)}& {\displaystyle }{\widehat{R}}_{\mathrm{Chl}}{\displaystyle }=\left({\widehat{\mathit{\mu }}}_{\mathrm{g}}+R_{\mathrm{M}}^{\mathrm{Chl}}\right){\mathit{\zeta }}_{\mathrm{Chl}}\widehat{\mathit{\theta }},\end{array}$$

where $R_{\mathrm{M}}^{\mathrm{Chl}}$ and ζ_Chl are the costs of chlorophyll maintenance and synthesis, respectively (Table 3).

Photo-acclimation 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:

$$\begin{array}{}\text{(24)}& \mathit{\theta }=f_{\mathrm{C}}\cdot \widehat{\mathit{\theta }}\mathit{\{}\mathrm{IA},\mathrm{DA}\mathit{\}}.\end{array}$$

Similarly, the overall respiratory cost of maintaining cellular chlorophyll is obtained by multiplying the chloroplast-specific cost by the size of the chloroplast:

$$\begin{array}{}\text{(25)}& R_{\mathrm{Chl}}=f_{\mathrm{C}}\cdot {\widehat{R}}_{\mathrm{Chl}}\mathit{\{}\mathrm{IA},\mathrm{DA}\mathit{\}}.\end{array}$$

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 chloroplast-specific 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)\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₀ 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 fixed-stoichiometry variant, $\widehat{\mathit{\theta }}$ is prescribed as the biomass-weighted average value calculated by the IA variant. Considering that θ is typically a constant “conversion factor” in classical, fixed-stoichiometry and fixed-Chl: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 biomass-weighted 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.

Table 3 Descriptions, values and units of model parameters regarding phytoplankton growth. Prescribed values for Q, K_N, f_V, and $\widehat{\mathit{\theta }}$ are based on the biomass-weighted averages estimated by the IA variant. All other parameter values are taken from within the published range (Pahlow et al., 2013; Smith et al., 2016), without particular reference to species.

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:

$$\begin{array}{}\text{(28)}& f\left(T\right)=\exp \left({\displaystyle \frac{-E_{\mathrm{a}}}{R}}\left[{\displaystyle \frac{\mathrm{1}}{T}}-{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{r}}}}\right]\right),\end{array}$$

where the gas constant R=8.3145 $\mathrm{J}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, and the activation energy, $E_{\mathrm{a}}=\mathrm{4.82}\times {\mathrm{10}}^{\mathrm{4}}$ J mol⁻¹, 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

$$\begin{array}{}\text{(29)}& \begin{array}{rl}I\left(z\right)=& I_{\mathrm{0}}\left[A\exp \left({\displaystyle \frac{-z}{{\mathit{\eta }}_{\mathrm{1}}}}\right)+(\mathrm{1}-A)\right.\\ & \left.\exp \left({\displaystyle \frac{-z}{{\mathit{\eta }}_{\mathrm{2}}}}-\underset{z}{\overset{\mathrm{0}}{\int }}\sum _i{k}_i{c}_i\left(z^{\prime }\right)\mathrm{d}{z}^{\prime }\right)\right],\end{array}\end{array}$$

where A, η₁ and η₂ 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² mmolN⁻¹ 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.0 m, η₂=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):

$$\begin{array}{}\text{(30)}& {\displaystyle \frac{\partial c_i}{\partial t}}+{\displaystyle \frac{\partial }{\partial z}}\left(w_i{c}_i-K_z{\displaystyle \frac{\partial c_i}{\partial z}}\right)=s\left(c_i\right),\end{array}$$

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⁻¹ 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).

Figure 2 Atmospheric variables. (a) Astronomically estimated instantaneous irradiance at the water surface and (b) prescribed air temperature.

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 third-year 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.

3 Results

Daytime-averaged 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⁻³) 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.

Figure 3 Abiotic environment. (a–c) daytime-averaged photosynthetically active radiation, $\overline{I}$ [E m⁻² d⁻¹], (d–f) water temperature T [^∘C] and (g–i) DIN [mmolN m⁻³], as simulated by the FS (a, d, g); IA (b, e, h) and DA (c, f, i) 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).

Figure 4 Phytoplankton C, N and Chl concentrations: (a–c) Phy_C [mmolC m⁻³], (d–f) Phy_N [mmolN m⁻³], (j–l) Phy_Chl [mgC m⁻³]; and phytoplankton N:C (Q) and Chl:C (Θ) ratios: (g–i) Q [molN molC⁻¹] and (m–o) Θ [gChl gC⁻¹], as simulated by the FS (left); IA (center) and DA (right) variants.

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 chloroplast-specific respiratory costs $\widehat{\mathit{\theta }}$ (Eq. 23). This can be seen in the flattening of the light-saturation 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).

Figure 5 Phytoplankton physiological variables. (a–c) Chlorophyll density in chloroplasts, $\widehat{\mathrm{\Theta }}$ [gChl gC⁻¹]; (d–f) fractional allocation to affinity, f_A [–]; (g–i) fractional allocation to nutrient uptake, f_V [–]; nutrient limitation term of the FS variant, L_N [–] (j) and fractional allocation to carbon fixation of the IA and DA variants, f_C [–] (k–l); and (m–o) light saturation of photosynthesis, L_I [–] as simulated by the FS (left); IA (center) and DA (right) variants.

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.

Figure 6 Phytoplankton growth, uptake and respiration rates. (a–c) Net growth rate, μ [d⁻¹], (d–f): specific uptake rate, V [mmolN mmolC⁻¹ d⁻¹] and respiration costs of (g–i) N uptake, R_N [d⁻¹] and (j–l) chlorophyll maintenance and synthesis, R_Chl [d⁻¹] as simulated by the FS (left); IA (center) and DA (right) variants.

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⁻³ near the surface. In the FS variant, DIN remains higher (minimum concentration of ∼0.7 mmolN m⁻³ 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 non-existent 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₀, 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.

Figure 7 Upper 50 m averages of critical variables. (a) DIN [mmolN m⁻³], (b) phytoplankton Q [mmolN mmolC⁻¹], (c) f_C [–] (in addition L_N [–] for FS, shown with pale broken line), and (d) Phy_C [mmolC m⁻³], Phy_N [mmolN m⁻³] (e), and μ [d⁻¹] (f), as simulated by the FS (dashed green line), IA (fine-dashed dark blue line) and DA (continuous orange line) variants.

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 chloroplast-specific growth rate, $\widehat{\mathit{\mu }}$, which is maximized for the acclimative variants through photo-acclimative flexibility (Sect. 2.2.4), is always higher than that calculated by the FS variant, as expected (not shown). As the chloroplast-specific 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.

Figure 8 Vertical profiles on 1 February (indicated as F in panel label), March (M), April (A), July (J), and 15 October (O) and November (N) for DIN [mmolN m⁻³] (a), phytoplankton Q [mmolN mmolC⁻¹] (b), f_C [–] (c), Phy_C [mmolC m⁻³] (d), Phy_N [mmolN m⁻³] (e), μ [d⁻¹] (f), and V [molN molN⁻¹ d⁻¹] (g) as simulated by the IA (fine-dashed dark blue line), DA (continuous orange line) and the FS (dashed green line) variants, when the prescribed Θ (Table 3) is scaled with f_C, according to Eq. (24).

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 nutrient-replenished 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 re-shuffling 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 time-lag effect: in the DA variant, the nutrient-starved 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 biomass-weighted 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).

Figure 9 Detrital C:N [molC molN⁻¹] (a–c) in the entire water column and (d–f) in the bottom layer (d–f), as simulated by the FS (a, d); IA (b, e) and DA (c, f) variants.

Simulated process rates determining ecosystem functioning, such as the water-column-integrated 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}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) and IA (45.66 $\mathrm{mmolC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) variants are, respectively, 8.1 % and 13.9 % smaller than that of the DA variant (53.06 $\mathrm{mmolC}{\mathrm{m}}^{-\mathrm{2}}{\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}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) is the highest, followed by the 8.2% lower IA (4.39 $\mathrm{mmolN}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) and 14.3 % lower FS (4.1 $\mathrm{mmolN}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) variants.

Figure 10 (a) Water-column-integrated NPP rate [mmolC m⁻² d⁻¹] and (b) water-column-integrated NDD rate [mmolN m⁻² d⁻¹], as simulated by the IA (fine-dashed dark blue line), DA (continuous orange line) and the FS (dashed green line) variants.

4 Discussion

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 down-regulation of nutrient uptake, which is needed to avoid unrealistically high nutrient quotas in a Droop scheme. Often, down-regulation is formulated as some function (linear, e.g., Grover, 1991 or non-linear, 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 down-regulation term nor any prescribed maximum quota value. This is because the optimization of growth, subject to the growth vs. nutrient uptake trade-off (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 trade-offs, 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, 1-D 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 3-D 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 3-D setups. Recent applications of these models in realistic 3-D 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 steady-state, this is a natural consequence of any “Droop-like” 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 0-D (box) setup. Here, for the first time, we have tested this approach in an idealized 1-D 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 1-D 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 3-D 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 nutrient-replete phase provides a competitive advantage during subsequent periods of nutrient scarcity (Grover, 1991; Litchman et al., 2009). Similarly, diffusion or active movement of nutrient-rich cells from the nutrient-replete to nutrient-rich 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 re-shuffling 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 non-acclimative 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 light-limited 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, trade-offs exist between multiple objectives. Such trade-offs 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 longer-term simulation), (ii) selection among existing species that had been pre-adapted to these conditions and (iii) individual-level acclimation. Optimality-based acclimative models can thus capture some key community-level 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 system-level 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}{\mathrm{m}}^{-\mathrm{2}}{\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}{\mathrm{m}}^{-\mathrm{2}}{\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 well-known 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 eco-evolutionary 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 well-adapted 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 near-zero 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 plankton-related 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. Fixed-stoichiometry models erroneously predict constant elemental composition of detrital matter production, as demonstrated by our FS variant in this study. The so-called 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 fixed-stoichiometry models.

4.4 Present implementation, challenges and perspectives

Moving a coupled hydrodynamic–biogeochemical models from a 0-D setup to a spatially explicit setup can be error prone and time consuming. FABM provides an easy-to-use 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 0-D to n-D 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 top-down 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 N-only 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 co-limitation 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 cost-effective 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 carbon-based 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 mass-balance 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 non-acclimative variant. Our acclimation scheme consists of four acclimative flexibilities: variability of internal nutrient quota, optimization of uptake vs. growth trade-off, optimization of maximum uptake vs. affinity trade-off 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.

5 Conclusions

In this study, we present a FABM implementation of the “NflexPD” model and the behavior of three variants it can emulate: a fixed-stoichiometry (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, 1-D 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 non-acclimative, fixed-stoichiometry 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 ecosystem-scale 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.

Appendix A: Details of derivations

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

$$\begin{array}{}\text{(A1)}& {\displaystyle }\mathit{\mu }{\displaystyle }={\widehat{\mathit{\mu }}}_{\mathrm{net}}\cdot L_{\mathrm{N}}-{\mathit{\zeta }}_{\mathrm{N}}\cdot \mathit{\mu }\cdot Q.\end{array}$$

Reorganizing yields

$$\begin{array}{}\text{(A2)}& {\displaystyle }\mathit{\mu }{\displaystyle }={\displaystyle \frac{{\widehat{\mathit{\mu }}}_{\mathrm{net}}\cdot L_{\mathrm{N}}}{\mathrm{1}+{\mathit{\zeta }}_{\mathrm{N}}\cdot Q}}.\end{array}$$

Substituting this with μ in

$$\begin{array}{}\text{(A3)}& {\displaystyle }{R}_{\mathrm{N}}{\displaystyle }={\mathit{\zeta }}_{\mathrm{N}}\cdot V={\mathit{\zeta }}_{\mathrm{N}}\cdot \mathit{\mu }\cdot Q,\end{array}$$

we obtain a V-independent expression for R_N:

$$\begin{array}{}\text{(A4)}& {\displaystyle }{R}_{\mathrm{N}}{\displaystyle }={\mathit{\zeta }}_{\mathrm{N}}\cdot {\displaystyle \frac{{\widehat{\mathit{\mu }}}_{\mathrm{net}}\cdot L_{\mathrm{N}}}{\mathrm{1}+{\mathit{\zeta }}_{\mathrm{N}}\cdot Q}}\cdot Q.\end{array}$$

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

$$\begin{array}{}\text{(A5)}& \mathit{\mu }=f_{\mathrm{C}}{\widehat{\mathit{\mu }}}_g-{\mathit{\zeta }}_N{f}_V\widehat{V}-({\widehat{\mathit{\mu }}}_g+R_M^{\mathrm{Chl}}){\mathit{\zeta }}_{\mathrm{Chl}}\widehat{\mathit{\theta }}{f}_{\mathrm{C}}.\end{array}$$

Reorganizing yields

$$\begin{array}{}\text{(A6)}& \mathit{\mu }=f_{\mathrm{C}}\left[{\widehat{\mathit{\mu }}}_g(\mathrm{1}-{\mathit{\zeta }}_{\mathrm{Chl}}\widehat{\mathit{\theta }})-{\mathit{\zeta }}_{\mathrm{Chl}}\widehat{\mathit{\theta }}{R}_M^{\mathrm{Chl}}\right]-{\mathit{\zeta }}_N{f}_V\widehat{V}.\end{array}$$

Substituting the term in square brackets with ${\widehat{\mathit{\mu }}}_{\mathrm{net}}$ based on Eq. (7) and expanding f_C using Eq. (11) yields

$$\begin{array}{}\text{(A7)}& \mathit{\mu }=\left(\mathrm{1}-{\displaystyle \frac{Q_{\mathrm{0}}}{\mathrm{2}Q}}-f_V\right){\widehat{\mathit{\mu }}}_{\mathrm{net}}-{\mathit{\zeta }}_N{f}_V\widehat{V}.\end{array}$$

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.

Appendix B: FS variant with a variable chloroplast size

Figure B1 Phytoplankton (a) Chl concentration, Phy_Chl [mgC m⁻³]; (b) net growth rate, μ [d⁻¹]; (c) Chl:C, Θ [gChl gC⁻¹]; (d) respiration cost of chlorophyll maintenance and synthesis, R_Chl [d⁻¹] as simulated by the FS variant, when the prescribed Θ (Table 3) is scaled with f_C, according to Eq. (24).

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 chloroplast-specific chlorophyll density and respiration costs in order to represent spatiotemporal variations of the cellular Chl:C ratio and proportional respiration costs.

Figure B2 Like Fig. 7 but when, for the FS variant, prescribed $\widehat{\mathit{\theta }}$ (Table 3) is scaled with L_N, i.e., replacing f_C with L_N in Eq. (24).

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 nutrient-rich 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₀ (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 nutrient-rich 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 nutrient-scarce 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.