FABM-NflexPD 2.0: testing an instantaneous acclimation approach for modeling the implications of phytoplankton eco-physiology for the carbon and nutrient cycles

The acclimative response of phytoplankton, which adjusts their nutrient and pigment content in response to changes in ambient light, nutrient levels, and temperature, is an important determinant of observed chlorophyll distributions and biogeochemistry. Acclimative models typically capture this response and its impact on the C : nutrient : Chl ratios of phytoplankton by explicitly resolving the dynamics of these constituents of phytoplankton biomass. The instantaneous acclimation (IA) approach only requires resolving the dynamics of a single tracer and calculates the elemental composition assuming instantaneous local equilibrium. IA can capture the acclimative response without substantial loss of accuracy in both 0D box models and spatially explicit 1D models. A major drawback of IA so far has been its inability to maintain mass balance for the elements with unresolved dynamics. Here we extend the IA model to capture both C and N cycles in a 0D setup, which requires analytical derivation of additional flux terms to account for the temporal changes in cellular N quota, Q. We present extensive tests of this model, with regard to the conservation of total C an N and its behavior in comparison to an otherwise equivalent, fully explicit dynamic acclimation (DA) variant under idealized conditions with variable light and temperature. We also demonstrate a modular implementation of this model in the Framework for Aquatic Biogeochemical Modelling (FABM), which facilitates modeling competition between an arbitrary number of different acclimative phytoplankton types. In a 0D setup, we did not find evidence for computational advantages of the IA approach over the DA variant. In a spatially explicit setup, performance gains may be possible but would require modifying the physical-flux calculations to account for spatial differences in Q between model grid cells.

1 Introduction

Elemental stoichiometry and pigment density of phytoplankton vary widely across environmental conditions, at both the physiological (e.g., Garcia et al., 2016) and the community level (e.g., Moreno and Martiny, 2018). The physiological flexibility is driven by an acclimative re-adjustment of cellular machinery to changes in the availability in nutrients and light and the fact that the various cellular functions have competing requirements for resources, an example being enzymes, rich in N (Geider and La Roche, 2002), needed for nutrient uptake and photosynthesis. The systematic differences in cellular composition between species may be explained by the typical composition of some species being better or worse suited than that of others for a given resource regime (Klausmeier et al., 2004; Arrigo, 2005; Burson et al., 2016).

The potential relevance of such variability in the cellular composition of phytoplankton for biogeochemical cycles was recognized decades ago (Redfield, 1934, 1958), and evidence has been building ever since (Lenton and Klausmeier, 2007; Bonachela et al., 2016; Pahlow et al., 2020). Accounting for the acclimative capacity of phytoplankton in models is relevant for predicting the response of ecosystems to environmental change (Kwiatkowski et al., 2018; Kerimoglu et al., 2018) and for model performance (Ayata et al., 2013; Kerimoglu et al., 2017; Chen and Smith, 2018). It can also endow models with desirable properties, such as improved model portability (Anugerahanti et al., 2021). However, mechanistic acclimative models typically require additional state variables: usually one for chlorophyll and one for each of the resolved nutrients (e.g., Geider et al., 1998; Flynn, 2003) but possibly even more (Bonachela et al., 2013; Wirtz and Kerimoglu, 2016; Inomura et al., 2020).

However, additional state variables can increase the computational costs significantly in spatially explicit setups (Fulton et al., 2003), especially for models with 100s of phytoplankton groups (e.g., Follows et al., 2007; Dutkiewicz et al., 2020). As a potential remedy to this problem, the instantaneous acclimation (IA) approach can be used, where the changes in cellular composition are not dynamically tracked but adjust instantaneously to the resource environment. For instance, the FlexPFT model (Smith et al., 2016) follows an IA approach, where the Chl:C and C:N ratios instantaneously assume the optimal ratios for balanced growth (Pahlow et al., 2013). Ward (2017) compared a fully explicit, classical Caperon–Droop model (Caperon, 1968; Droop, 1968) to its IA counterpart and found that across a range of environmental settings, the predictions of the two approaches matched closely in a 0D setup. Kerimoglu et al. (2021) introduced FABM-NflexPD 1.0, a 1D setup of the FlexPFT model in the Framework of Aquatic Biogeochemical Models (FABM; Bruggeman and Bolding, 2014), and showed that the predictions of the IA variant of the model largely matched those of the fully explicit dynamic acclimation (DA) counterpart, except for minor differences during the transitions from winter to spring and from autumn to winter.

FABM-NflexPD 1.0 only tracks N, which may be sufficient for some ecological applications but not for applications that require mass balance for (multiple) nutrients and carbon. Here we introduce FABM-NflexPD 2.0, which tracks both C and N. We present a detailed description of a C-based version, which is extended to account also for N fluxes resulting from instantaneous changes in cell quotas, such that mass balance is maintained for both C and N. We evaluate the consistency and robustness of the model by means of the following formal tests:

T1

– assessment of the conservation of C and N in a simplified version of the model in a 0D setup, where temperature and day length are held constant and light is provided as a sinusoidal function of time of year;

T2

– test of the IA approach in a more realistic setup, where temperature and day length also vary over time, while also accounting for light attenuation by phytoplankton, and comparison against a fully explicit DA variant;

T3

– simulation with multiple phytoplankton groups;

T4

– simulation in an open system where N and C are not conserved.

2 Model description

FABM-NflexPD 2.0 differs from FABM-NflexPD 1.0 (K21) mainly by tracing the dissolved inorganic carbon pool, DIC, such that now the model is ideally able to conserve both the total C and N in the system and not only N as in K21. Moreover, we consider dilution/mixing and sinking terms for simulating an open system, such as a chemostat or a surface mixed layer (SML). As in K21, we consider here a DA variant to compare to the IA variant but not a fixed stoichiometry variant unlike in K21. As a final difference, for the IA, we trace the C content in phytoplankton with a state variable instead of the N content as in K21. The rates of change in the state variables are

$$\begin{array}{}\text{(1a)}& {\displaystyle \frac{\mathrm{d}{\mathrm{Phy}}_{\mathrm{C}}}{\mathrm{d}t}}=F_{\mathrm{DIC}-{\mathrm{Phy}}_{\mathrm{C}}}-F_{{\mathrm{Phy}}_{\mathrm{C}}-{\mathrm{Det}}_{\mathrm{C}}}-D\cdot {\mathrm{Phy}}_{\mathrm{C}},\text{(1b)}& {\displaystyle \frac{\mathrm{d}{\mathrm{Phy}}_{\mathrm{N}}}{\mathrm{d}t}}=\underset{\mathrm{Netuptake}}{\underbrace{F_{\mathrm{DIN}-{\mathrm{Phy}}_{\mathrm{N}}}}}-\underset{\mathrm{Mortality}}{\underbrace{F_{{\mathrm{Phy}}_{\mathrm{N}}-{\mathrm{Det}}_{\mathrm{N}}}}}-\underset{\mathrm{Dilution}}{\underbrace{D\cdot {\mathrm{Phy}}_{\mathrm{N}}}},\mathit{\left\{}\mathrm{DA}\mathit{\right\}}\end{array}$$

$$\begin{array}{}\text{(2a)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\mathrm{Det}}_{\mathrm{C}}}{\mathrm{d}t}}& =F_{{\mathrm{Phy}}_{\mathrm{C}}-{\mathrm{Det}}_{\mathrm{C}}}-F_{{\mathrm{Det}}_{\mathrm{C}}-\mathrm{DOC}}-D\cdot {\mathrm{Det}}_{\mathrm{C}}\\ & -{\displaystyle \frac{w_{\mathrm{Det}}}{H_{\mathrm{SML}}}}\cdot {\mathrm{Det}}_{\mathrm{C}},\end{array}\text{(2b)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\mathrm{Det}}_{\mathrm{N}}}{\mathrm{d}t}}& =F_{{\mathrm{Phy}}_{\mathrm{N}}-{\mathrm{Det}}_{\mathrm{N}}}-\underset{\mathrm{Hydrolysis}}{\underbrace{F_{{\mathrm{Det}}_{\mathrm{N}}-\mathrm{DON}}}}-D\cdot {\mathrm{Det}}_{\mathrm{N}}\\ & -\underset{\mathrm{Sinking}}{\underbrace{{\displaystyle \frac{w_{\mathrm{Det}}}{H_{\mathrm{SML}}}}\cdot {\mathrm{Det}}_{\mathrm{N}}}},\end{array}\end{array}$$

$$\begin{array}{}\text{(3a)}& {\displaystyle \frac{\mathrm{d}\mathrm{DOC}}{\mathrm{d}t}}& {\displaystyle }=F_{{\mathrm{Det}}_{\mathrm{C}}\text{-DOC}}-F_{\text{DOC-DIC}}-D\cdot \mathrm{DOC},\text{(3b)}& {\displaystyle \frac{\mathrm{d}\mathrm{DON}}{\mathrm{d}t}}& {\displaystyle }=F_{{\mathrm{Det}}_{\mathrm{N}}\text{-DON}}-\underset{\mathrm{Remineralization}}{\underbrace{F_{\text{DON-DIN}}}}-D\cdot \mathrm{DON},\end{array}$$

$$\begin{array}{}\text{(4a)}& {\displaystyle \frac{\mathrm{d}\mathrm{DIC}}{\mathrm{d}t}}=F_{\mathrm{DOC}-\mathrm{DIC}}-F_{\mathrm{DIC}-{\mathrm{Phy}}_{\mathrm{C}}}-D\cdot (\mathrm{DIC}-{\mathrm{DIC}}_{\mathrm{in}}),\text{(4b)}& {\displaystyle \frac{\mathrm{d}\mathrm{DIN}}{\mathrm{d}t}}=F_{\mathrm{DON}-\mathrm{DIN}}-F_{\mathrm{DIN}-{\mathrm{Phy}}_{\mathrm{N}}}-D\cdot (\mathrm{DIN}-{\mathrm{DIN}}_{\mathrm{in}}),\end{array}$$

where F_x−y is the flux from x to y, where x and y are state variables, except for Phy_N in the IA variant, which is defined as

$$\begin{array}{}\text{(5)}& {\mathrm{Phy}}_{\mathrm{N}}={\mathrm{Phy}}_{\mathrm{C}}\cdot Q,\mathit{\left\{}\mathrm{IA}\mathit{\right\}}\end{array}$$

where Q is the phytoplankton N quota (N:C ratio). Dilution and sinking terms describe fluxes in and out of the system and are non-zero only for the test T4 (see below).

2.1 C and N content of phytoplankton

As mentioned above, in this C-based version of the IA model variant, only the C content of phytoplankton (Phy_C) is dynamically tracked via Eq. (1a), whereas Phy_N is defined as a function of Q in Eq. (5). Q adjusts instantaneously to its optimal value for balanced growth, as determined by nutrient uptake in the protoplast, ${\widehat{V}}_{\mathrm{N}}$, and net photosynthesis in the chloroplast, ${\widehat{\mathit{\mu }}}_{\mathrm{net}}$ (Eq. 10 in K21):

$$\begin{array}{}\text{(6)}& Q={\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}}_{\mathrm{N}}+{\mathit{\zeta }}_{\mathrm{N}})}}}\right],\mathit{\left\{}\mathrm{IA}\mathit{\right\}}\end{array}$$

where Q₀ and ζ_N are the subsistence N quota and cost of N uptake, respectively (Table 1). The first term in Eq. (1a), $F_{\mathrm{DIC}-{\mathrm{Phy}}_{\mathrm{C}}}$, represents net phytoplankton growth:

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

The second term in Eq. (1a), $F_{{\mathrm{Phy}}_{\mathrm{C}}-{\mathrm{Det}}_{\mathrm{C}}}$, represents the mortality of phytoplankton:

$$\begin{array}{}\text{(8)}& F_{{\mathrm{Phy}}_{\mathrm{C}}-{\mathrm{Det}}_{\mathrm{C}}}=m_{\mathrm{C}}\cdot {\mathrm{Phy}}_{\mathrm{C}}^{\mathrm{2}}.\end{array}$$

Its N counterpart is found by multiplication with Q:

$$\begin{array}{}\text{(9)}& F_{{\mathrm{Phy}}_{\mathrm{N}}-{\mathrm{Det}}_{\mathrm{N}}}=F_{{\mathrm{Phy}}_{\mathrm{C}}-{\mathrm{Det}}_{\mathrm{C}}}\cdot Q,\end{array}$$

where m_C [${\mathrm{m}}^{\mathrm{3}}{\mathrm{mmolC}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$] is the C-based specific mortality rate. The hydrolysis and remineralization fluxes are calculated as first-order reactions:

$$\begin{array}{}\text{(10)}& {\displaystyle }{F}_{{\mathrm{Det}}_X-\mathrm{DO}X}& {\displaystyle }=r_{\mathrm{hyd}}\cdot {\mathrm{Det}}_X,X\in \mathit{\{}\mathrm{C},\mathrm{N}\mathit{\}},\text{(11)}& {\displaystyle }{F}_{\mathrm{DO}X-\mathrm{DI}X}& {\displaystyle }=r_{\mathrm{rem}}\cdot \mathrm{DO}X,X\in \mathit{\{}\mathrm{C},\mathrm{N}\mathit{\}}.\end{array}$$

2.2 N fluxes between DIN and phytoplankton

For tracking the Phy_N for the DA variant (Eq. 1b) and the dissolved organic nitrogen, DIN, for both variants (Eq. 4b), the flux from DIN to Phy_N needs to be known. As in K21, for the DA variant, it is simply the product of a specific uptake rate and the phytoplankton C biomass:

$$\begin{array}{}\text{(12)}& F_{\mathrm{DIN}-{\mathrm{Phy}}_{\mathrm{N}}}=V_{\mathrm{N}}\cdot {\mathrm{Phy}}_{\mathrm{C}}.\mathit{\left\{}\mathrm{DA}\mathit{\right\}}\end{array}$$

For the IA variant, the exact value of this flux is unknown due to the non-existent Phy_N pool and the corresponding flux. By substituting Phy_N with Phy_C⋅Q and applying the product rule we get

$$\begin{array}{}\text{(1b.2)}& {\displaystyle \frac{\mathrm{d}{\mathrm{Phy}}_{\mathrm{N}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}({\mathrm{Phy}}_{\mathrm{C}}\cdot Q)}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}{\mathrm{Phy}}_{\mathrm{C}}}{\mathrm{d}t}}\cdot Q+{\mathrm{Phy}}_{\mathrm{C}}\cdot {\displaystyle \frac{\mathrm{d}Q}{\mathrm{d}t}};\end{array}$$

here, the first term reflects the N equivalent of the change in Phy_C, i.e., Eq. (1a), and the second term describes the effect of the change in quota over time due to imbalances between C and N uptake. Substituting Eq. (1a) in Eq. (1b.2):

$$\begin{array}{}\text{(1b.3)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\mathrm{Phy}}_{\mathrm{N}}}{\mathrm{d}t}}& =\left[F_{\mathrm{DIC}-{\mathrm{Phy}}_{\mathrm{C}}}\cdot Q-F_{{\mathrm{Phy}}_{\mathrm{C}}-{\mathrm{Det}}_{\mathrm{C}}}\cdot Q-D\cdot {\mathrm{Phy}}_{\mathrm{C}}\cdot Q\right]\\ & +{\mathrm{Phy}}_{\mathrm{C}}\cdot {\displaystyle \frac{\mathrm{d}Q}{\mathrm{d}t}}.\end{array}\end{array}$$

and plugging Eqs. (7) and (9) into Eq. (1b.3) yields

$$\begin{array}{}\text{(1b.4)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\mathrm{Phy}}_{\mathrm{N}}}{\mathrm{d}t}}& =\left[\mathit{\mu }\cdot {\mathrm{Phy}}_{\mathrm{C}}\cdot Q-F_{{\mathrm{Phy}}_{\mathrm{N}}-{\mathrm{Det}}_{\mathrm{N}}}-D\cdot {\mathrm{Phy}}_{\mathrm{C}}\cdot Q\right]\\ & +{\mathrm{Phy}}_{\mathrm{C}}\cdot {\displaystyle \frac{\mathrm{d}Q}{\mathrm{d}t}}.\end{array}\end{array}$$

In Eq. (1b.4), $\mathit{\mu }\cdot {\mathrm{Phy}}_{\mathrm{C}}\cdot Q$ can be identified with $F_{\mathrm{DIN}-{\mathrm{Phy}}_{\mathrm{N}}}$ in the IA variant because $V_{\mathrm{N}}=\mathit{\mu }\cdot Q$ for balanced growth. Following Smith et al. (2016), we assign the last term to $F_{\mathrm{DIN}-{\mathrm{Phy}}_{\mathrm{N}}}$ too, yielding a re-definition of $F_{\mathrm{DIN}-{\mathrm{Phy}}_{\mathrm{N}}}$ for the IA variant:

$$\begin{array}{}\text{(12.2)}& F_{\mathrm{DIN}-{\mathrm{Phy}}_{\mathrm{N}}}={\mathrm{Phy}}_{\mathrm{C}}\left[\mathit{\mu }Q+{\displaystyle \frac{\mathrm{d}Q}{\mathrm{d}t}}\right],\mathit{\left\{}\mathrm{IA}\mathit{\right\}}\end{array}$$

which essentially redirects part of the fluxes associated with Phy_N to DIN. Plugging this into Eq. (4b) and recognizing that $\frac{\mathrm{d}Q}{\mathrm{d}t}$ consists of partial derivatives with respect to DIN, daily average irradiance, $\stackrel{\mathrm{‾}}{I}$, fractional day length, L_D, and temperature, T (see Appendix A),

$$\begin{array}{}\text{(4b.2)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}\mathrm{DIN}}{\mathrm{d}t}}& =F_{\mathrm{DON}-\mathrm{DIN}}-{\mathrm{Phy}}_{\mathrm{C}}\left[\mathit{\mu }Q+{\displaystyle \frac{\mathrm{d}Q}{\mathrm{d}t}}\right]\\ & -D\cdot (\mathrm{DIN}-{\mathrm{DIN}}_{\mathrm{in}})\\ & =F_{\mathrm{DON}-\mathrm{DIN}}-{\mathrm{Phy}}_{\mathrm{C}}\left[\mathit{\mu }Q+{\displaystyle \frac{\partial Q}{\partial \mathrm{DIN}}}{\displaystyle \frac{\mathrm{d}\mathrm{DIN}}{\mathrm{d}t}}\right.\\ & \left.+{\displaystyle \frac{\partial Q}{\partial \stackrel{\mathrm{‾}}{I}}}{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{I}}{\mathrm{d}t}}+{\displaystyle \frac{\partial Q}{\partial L_{\mathrm{D}}}}{\displaystyle \frac{\mathrm{d}{L}_{\mathrm{D}}}{\mathrm{d}t}}+{\displaystyle \frac{\partial Q}{\partial T}}{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}\right]\\ & -D\cdot (\mathrm{DIN}-{\mathrm{DIN}}_{\mathrm{in}}),\mathit{\left\{}\mathrm{IA}\mathit{\right\}}\end{array}\end{array}$$

and reorganizing the $\frac{\mathrm{d}\mathrm{DIN}}{\mathrm{d}t}$ on the right-hand side

$$\begin{array}{}\text{(4b.3)}& {\displaystyle \frac{\mathrm{d}\mathrm{DIN}}{\mathrm{d}t}}={\displaystyle \frac{\begin{array}{c}{F}_{\mathrm{DON}-\mathrm{DIN}}-{\mathrm{Phy}}_{\mathrm{C}}\left[\mathit{\mu }Q+\frac{\partial Q}{\partial \stackrel{\mathrm{‾}}{I}}\frac{\mathrm{d}\stackrel{\mathrm{‾}}{I}}{\mathrm{d}t}+\frac{\partial Q}{\partial L_{\mathrm{D}}}\frac{\mathrm{d}{L}_{\mathrm{D}}}{\mathrm{d}t}\right.\\ \left.+\frac{\partial Q}{\partial T}\frac{\mathrm{d}T}{\mathrm{d}t}\right]-D\cdot (\mathrm{DIN}-{\mathrm{DIN}}_{\mathrm{in}})\end{array}}{\mathrm{1}+{\mathrm{Phy}}_{\mathrm{C}}\frac{\partial Q}{\partial \mathrm{DIN}}}}.\mathit{\left\{}\mathrm{IA}\mathit{\right\}}\end{array}$$

The partial derivatives of Q with respect to DIN, $\stackrel{\mathrm{‾}}{I}$, L_D, and T are obtained by canonical application of the chain rule, as detailed in Appendix A. The final terms required in Eq. (4b.3) are the changes in $\stackrel{\mathrm{‾}}{I}$, L_D, and T over time, i.e., $\mathrm{d}\stackrel{\mathrm{‾}}{I}/\mathrm{d}t$, $\mathrm{d}{L}_{\mathrm{D}}/\mathrm{d}t$, and $\mathrm{d}T/\mathrm{d}t$. When the irradiance and temperature are supplied externally, as is typically the case in coupled physical–biological models, it is not possible to obtain their temporal derivatives analytically. Hence they are numerically approximated as the discrete backward difference between their current (t=i) and previous ($t=i-\mathrm{1}$) values, divided by the integration time step, i.e., $\mathrm{d}E/\mathrm{d}t\approx (E_i-E_{i-\mathrm{1}})/\mathrm{\Delta }t$ for $E=\mathit{\{}\stackrel{\mathrm{‾}}{I},L_{\mathrm{D}},T\mathit{\}}$ (see also Sect. 2.3.1). Finally, for the case of multiple phytoplankton functional types (PFTs) indexed by j, Eq. (4b.3) can be generalized as follows:

$$\begin{array}{}\text{(4b.4)}& {\displaystyle \frac{\mathrm{d}\mathrm{DIN}}{\mathrm{d}t}}={\displaystyle \frac{\begin{array}{c}{F}_{\mathrm{DON}-\mathrm{DIN}}-{\sum }_j{\mathrm{Phy}}_{\mathrm{C}}^j\left[{\mathit{\mu }}^j{Q}^j+\frac{\partial Q^j}{\partial \stackrel{\mathrm{‾}}{I}}\frac{\mathrm{d}\stackrel{\mathrm{‾}}{I}}{\mathrm{d}t}\right.\\ \left.+\frac{\partial Q^j}{\partial L_{\mathrm{D}}}\frac{\mathrm{d}{L}_{\mathrm{D}}}{\mathrm{d}t}+\frac{\partial Q^j}{\partial T}\frac{\mathrm{d}T}{\mathrm{d}t}\right]-D\cdot (\mathrm{DIN}-{\mathrm{DIN}}_{\mathrm{in}})\end{array}}{\mathrm{1}+{\sum }_j\left[{\mathrm{Phy}}_{\mathrm{C}}^j\frac{\partial Q^j}{\partial \mathrm{DIN}}\right]}}.\mathit{\left\{}\mathrm{IA}\mathit{\right\}}\end{array}$$

As a technical remark regarding the FABM implementation of the model, Eqs. (4b.3)–(4b.4) require the combination of terms which are calculated by separate abiotic and phytoplankton modules in K21. Therefore, in order to avoid a circular-dependency error in the current implementation, an intermediate module collects the necessary terms from the two modules and sets the right-hand sides for DIN at once.

Table 1 Descriptions, values, and units of model parameters regarding phytoplankton growth. Parameters with prime (C^′) are for a cell with an equivalent spherical diameter (ESD) of 1 µm, which is the size assumed for experiments T1 and T2. For T3, where different size classes are simulated, the respective values are obtained according to $C=C^{\prime }\cdot {\mathrm{ESD}}^{S_{\mathrm{C}}}$, where S_C is the allometric scaling coefficient for this parameter. Values for C^′ and S_C are as in Smith et al. (2016), and other parameters as in Kerimoglu et al. (2021).

2.3 Test setups and model operation

For all tests, the model is operated in a spatially homogeneous one-box setup, using the 0D driver of FABM (Bruggeman and Bolding, 2014). With this 0D setup, numerical solutions are obtained using a fourth-order Runge–Kutta method with a time step of 60 s. Model forcing applied in our 0D setup varies among different tests, as explained below.

2.3.1 T1

This test is designed to assess the effect of inaccuracies incurred by the numerical approximation of the time derivatives of external forcing variables. We consider two cases with regard to photosynthetically active irradiance (PAR): for PAR:N (numerical) the time derivative is approximated numerically and for PAR:A (analytical) it is obtained analytically. In both cases, irradiance ($\stackrel{\mathrm{‾}}{I}$) is provided (as implemented in K21) as a sinusoidal function of day of year (t) to represent a seasonal cycle typical of a high-latitude environment in the Northern Hemisphere:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}\stackrel{\mathrm{‾}}{I}\left(t\right)& ={\stackrel{\mathrm{‾}}{I}}_{min}+({\stackrel{\mathrm{‾}}{I}}_{max}-{\stackrel{\mathrm{‾}}{I}}_{min})\mathrm{0.5}\left[\mathrm{1}+\sin \left[\mathrm{2}\mathit{\pi }{t}^{\prime }\right]\right],\\ & t^{\prime }={\displaystyle \frac{t}{\mathrm{365}}}-\mathrm{0.25},\end{array}\end{array}$$

where ${\stackrel{\mathrm{‾}}{I}}_{min}=\mathrm{1.6}\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ and ${\stackrel{\mathrm{‾}}{I}}_{max}=\mathrm{110}\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ define the minimum and maximum values throughout the year and t^′ represents the relative day of the year delayed by a quarter cycle to obtain the peak value in the middle of the year (see Fig. 1 for the behavior of the function with these parameters). For simplicity, we assume that temperature, T, is fixed at 10 ^∘C, and fractional day length, L_D, is unity, and we ignore light attenuation.

PAR:N

– in the first case, the time derivative of $\stackrel{\mathrm{‾}}{I}$ is calculated numerically as a finite-difference approximation:

$$\begin{array}{}\text{(14a)}& {\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{I}}{\mathrm{d}t}}\approx {\displaystyle \frac{{\stackrel{\mathrm{‾}}{I}}_i-{\stackrel{\mathrm{‾}}{I}}_{i-\mathrm{1}}}{\mathrm{\Delta }t}},\end{array}$$

where i is the time step index and Δt the time step of the numerical integration.

PAR:A

– in the second case, the temporal derivative of shortwave radiation is calculated analytically:

$$\begin{array}{}\text{(14b)}& {\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{I}}{\mathrm{d}t}}=({\stackrel{\mathrm{‾}}{I}}_{max}-{\stackrel{\mathrm{‾}}{I}}_{min}){\displaystyle \frac{\mathit{\pi }}{\mathrm{365}}}\cos \left[\mathrm{2}\mathit{\pi }{t}^{\prime }\right].\end{array}$$

The temporal derivatives found by the numerical and analytical approaches are almost identical (Fig. 1).

Figure 1 Daily average irradiance and its temporal derivative, as extracted from the simulation outputs generated for T1. PAR:N (solid blue line) is the model version where the temporal derivative of irradiance is approximated numerically; PAR:A (dashed orange line): both irradiance and its temporal derivative are calculated analytically.

2.3.2 T2 – T3

For these tests, we apply the numerical and analytical time derivatives for $\stackrel{\mathrm{‾}}{I}$ (Eqs. 14a and 14b) and also for variable day length, L_D, and temperature, T. L_D, as described by Forsythe et al. (1995), is used to calculate the instantaneous irradiance, based on the same irradiance function as in T1 (Eq. 13). The seasonal variability in T is represented by a sinusoidal function analogous to Eq. (13), with $T_{min}=\mathrm{2}$ and $T_{max}=\mathrm{20}$ ^∘C.

For T3, we compare simulations for 10 phytoplankton groups with the IA and DA variants. Phytoplankton groups represent different size classes across the range of 0.2–100 µm equivalent spherical diameter (ESD), uniformly spaced on a logarithmic scale. As in Smith et al. (2016), Q₀, ${\widehat{V}}_{\mathrm{0}}$, and ${\widehat{A}}_{\mathrm{0}}$ vary according to allometric relationships (Table 1). Scalings of Q₀ and ${\widehat{V}}_{\mathrm{0}}$ are based on a combination of cell-specific scalings of subsistence quotas (Edwards et al., 2012, “marine species”), maximum uptake rates (Marañón et al., 2013), and cell-specific C content (Menden-Deuer and Lessard, 2000, “protist plankton excluding diatoms”). Scaling of ${\widehat{A}}_{\mathrm{0}}$ is based on heuristics (Smith et al., 2014).

As a technical note regarding the implementation, using the phy_Cbased.F90 and abio_Cbased.F90 modules, the number of phytoplankton types can be modified without changing or recompiling the code, by adjusting the configuration file (see the fabm.yaml examples in the testcases folder that were employed to produce the results presented in this study), which is a feature of the modularity of FABM.

2.3.3 T4

In a closed system, where all mass is conserved, DIN can be calculated directly as the difference between the initial total mass and the sum of all other pools (e.g., $\mathrm{DIN}=\mathrm{TotalN}-\mathrm{DON}-{\mathrm{Det}}_{\mathrm{N}}-Q\cdot {\mathrm{Phy}}_{\mathrm{C}}$). This would eliminate the necessity of deriving the additional differentials in Eq. (4b.3) (solutions of which are provided in Appendix A), and the resulting code could be significantly faster, owing to one less state variable (i.e., DIN) and lower amounts of logic and calculations. However, this would work only for closed systems.

The aim of T4 is to evaluate the behavior of the model in open systems, such as chemostats, using a non-zero D in Eqs. (1a)–(4a) to represent dilution or the dynamics within a surface mixed layer (SML), using D>0 and w_Det>0, to represent mixing with the layer below the SML and sedimentation of detritus out of the SML, respectively. In order to characterize an aquatic environment in a temperate climate zone that undergoes thermal stratification in summer, we consider a cyclic seasonal pattern in D, with values approaching $D_{max}=\mathrm{1.0}$ during winter and $D_{min}=\mathrm{0.001}$ during summer:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}D& =D_{min}+\mathrm{0.25}(D_{max}-D_{min}){\left[\mathrm{1}-\sin \left(\mathrm{2}\mathit{\pi }{t}^{\prime }\right)\right]}^{\mathrm{2}},\\ & t^{\prime }={\displaystyle \frac{t}{\mathrm{365}}}-\mathrm{0.15},\end{array}\end{array}$$

where t^′ is the relative day of the year, adjusted to mimic initiation of stratification by the beginning of April.

In order to examine the mass balance of the model in such a setup, we introduce two new state variables, Ext_C and Ext_N, which trace the amounts of N and C exported from and imported into the system:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\mathrm{Ext}}_X}{\mathrm{d}t}}& =D\cdot \left({\mathrm{Phy}}_X+{\mathrm{Det}}_X+\mathrm{DO}X+(\mathrm{DI}X-\mathrm{DI}{X}_{\mathrm{in}})\right)\\ & +{\displaystyle \frac{w_{\mathrm{Det}}}{H_{\mathrm{SML}}}}\cdot {\mathrm{Det}}_X,X\in \mathit{\{}\mathrm{C},\mathrm{N}\mathit{\}},\end{array}\end{array}$$

such that the global amounts of N and C, i.e., the sums of these variables and the corresponding C and N variables in Eqs. (1a)–(4), should be conserved. See Table 1 for the values of the additional parameters that describe the dilution and sinking fluxes.

For this test, we consider two PFTs with ESDs of 1 and 10 µm, with Q₀, ${\widehat{V}}_{\mathrm{0}}$, and ${\widehat{A}}_{\mathrm{0}}$ scaled as explained for T3 above. The configuration files for this test are provided with the code (see the “Code availability” section).

3 Results

3.1 Accuracy of numerical approximation of the temporal derivative of light (T1)

The model is conservative with respect to C and mostly also N (Fig. 2), for both numerically approximated (PAR:N) and analytically calculated (PAR:A) temporal derivatives of irradiance (see the details in Sect. 2.3.1). The range of deviation (difference between maximum and minimum values obtained throughout the run) of total N is about 2-fold higher in the PAR:N run than in the PAR:A run (see Fig. 2) and corresponds to about 0.0003 % of the total N in the system. This suggests that the numerical approximation of the derivative of light does introduce some additional error.

Figure 2 T1: carbon (a) and nitrogen (b) pools for the PAR:N (solid blue line) and PAR:A (dashed orange line) simulations.

3.2 Testing the IA approach in a more realistic setup, in comparison to the fully explicit DA approach (T2)

For T2 we consider seasonal variations also in T and L_D (Fig. 3). Note that variations in L_D also affect the seasonal cycles of $\stackrel{\mathrm{‾}}{I}$ and $\mathrm{d}\stackrel{\mathrm{‾}}{I}/\mathrm{d}t$.

Figure 3 Daily average irradiance and temperature and their numerically approximated temporal derivatives used in T2.

The IA and DA variants produce almost identical results (Figs. 4, 5). This similarity is expected (Ward, 2017). On closer inspection, some differences can be detected, such as slightly higher Phy_N at the peak of the spring bloom and slightly higher Det_N shortly after the spring bloom in the DA variant. The differences are due to the re-allocation of part of the fluxes between Phy_N and DIN according to Eq. (4b.3). They remain relatively small because (1) the timescale of the optimal regulation of N uptake in the DA variant is short relative to the timescale of phytoplankton growth and the DIN changes in our simulations and (2) the strong interaction between phytoplankton and DIN leads to a negative feedback between the deviations between the IA and DA variants and the extra DIN fluxes caused by variations in Q in the IA variant.

Figure 4 T2: carbon (left) and nitrogen (right) pools for the IA (solid blue line) and DA (dashed orange line) variants with variable day length and temperature.

Figure 5 T2: differences between the IA and DA variants for the quantities shown in Fig. 4.

3.3 T3: comparing DA and IA variants in simulating multiple PFTs

Figure 6 shows results of experiment T3 with 10 phytoplankton size classes for our IA and DA variants. Annually averaged concentrations decrease with cell size, and the larger classes exhibit stronger seasonal relative variations. Under other environmental conditions, e.g., different initial conditions or temporal variability, different outcomes can emerge (see, e.g., Taherzadeh et al., 2017), but this is outside of the scope of our current study. C biomass of phytoplankton by the two variants is near-identical. However, differences do exist during the spring bloom, with concentrations of smaller groups being higher in the IA than in the DA variant and vice versa for the larger groups. Concentrations of the larger groups are higher in the DA variant, likely because the extra DIN derived in the IA variant from the increasing cell quotas during the spring bloom benefits the larger cells less than the smaller cells as imposed by the allometric relationships (as in Grover, 1991; Litchman et al., 2009).

Figure 6 T3: (a, b) C biomass of 10 phytoplankton size classes in the IA (solid line) and DA (dotted line) variants (a, c) and the difference between IA and DA (b, d). (c, d) Sums of phytoplankton C biomass simulated by the IA and DA variants (a, c) and their differences (b, d).

3.4 T4: comparing DA and IA variants in a non-closed system

In T4, we consider competition between two species in an open system forced by fluxes to and from an external environment. We simulate a surface mixed layer (SML), where mixing with the deeper layer can introduce new nutrients (DIN and DIC) and dilute all other variables, and sedimentation exports detrital C and N out of the system. As typically observed in aquatic environments located in temperate climate zones, we prescribe a seasonally varying mixing coefficient, with lower values during summer due to thermal stratification.

Under such a regime, the system captures the characteristic features of a temperate aquatic environment, with low phytoplankton biomass during winter, a strong spring bloom, depletion of DIN within the SML during summer, and an autumn bloom. Accordingly, the total N in the system shows a strong seasonal pattern (Fig. 7). However, when the external N is taken into account, the global amount of N is consistently conserved by the model also for an open system.

Figure 7 T4: annual variations in global C and N (totalC+Ext_C and totalN+Ext_N, see Eq. 16), total C and N ($\sum {\mathrm{Phy}}_{\mathrm{C}}^j+\mathrm{DIC}+\mathrm{DOC}+{\mathrm{Det}}_{\mathrm{C}}$ and $\sum {\mathrm{Phy}}_{\mathrm{N}}^j+\mathrm{DIN}+\mathrm{DON}+{\mathrm{Det}}_{\mathrm{N}}$), and other state variables that trace individual C and N pools for a seasonally varying mixing regime in an open system.

4 Discussion

In this study, we present FABM-NflexPD 2.0, a FABM implementation of the FlexPFT model introduced by Smith et al. (2016) with a few minor corrections (see the notes at the end of Appendix A). The precursor, FABM-NflexPD 1.0 (K21, Kerimoglu et al., 2021), only resolves the N cycle and does not close the C cycle. FABM-NflexPD 2.0, which we present here, can resolve both N and C cycles in a 0D setup, owing to an additional flux term to maintain the mass balance of N (Sects. 2.2 and Appendix A).

4.1 Re-establishing the mass balance

Two variants of the model are elaborated here. The dynamic acclimation (DA) variant is fully explicit in its treatment of the C and N content of phytoplankton, as in the model by Fernández-Castro et al. (2016). The instantaneous acclimation (IA) variant aims to track the phytoplankton dynamics with a single state variable, based on a balanced growth approximation (Burmaster, 1979), that is, assuming that cells reach the equilibrium state instantaneously.

Owing to the lack of a state variable Phy_N, Eqs. (1a)–(4) do not preserve mass (total N) because they ignore the contribution of the rate of change of Q to the rate of change of Phy_N. It is impossible to maintain mass balance mechanistically without adding a state variable Phy_N. However, we can re-establish mass balance to a large extent by assigning the missing N flux to the DIN compartment in Eq. (4b.3). While this may be a rather arbitrary measure for achieving N mass balance and also violates the assumptions behind the model by assigning part of Phy_N to DIN, the resulting differences compared to explicitly resolving Phy_N are relatively small, as was also shown previously (Ward, 2017).

Through detailed tests, we show that N is conserved to a very large degree (for all tests we conducted, maximum error was 0.0063 %, which was for T1). We also show that the predictions of the IA and DA variants are mostly indistinguishable. Finally, with a simulation of 10 phytoplankton size classes with our two model variants, we demonstrate that the model is well aligned with the modular coupling philosophy of FABM (Bruggeman and Bolding, 2014).

As explained in Sect. 2.2, reducing the errors in mass balance for N requires explicitly calculating the changes in N quota over time, i.e., $\mathrm{d}Q/\mathrm{d}t$, which in turn requires the calculation of individual components of this change driven by different environmental factors, namely, DIN, $\stackrel{\mathrm{‾}}{I}$, T, and L_D. Under the idealized setup of T2, changes in DIN are clearly the dominant source of variation in Q. However, contributions by other factors are non-negligible (Fig. 8). In other setups, the relative importance of various environmental factors may be different. The relevance of these secondary factors to the elemental stoichiometry of phytoplankton is an often neglected aspect. We expect our mathematically explicit treatment of this issue to inspire and contribute to future endeavors to establish an analytical framework for investigating the mechanistic underpinnings of plankton physiology.

Figure 8 T2: total $\mathrm{d}Q/\mathrm{d}t$ and its components as contributed by the changes in DIN, $\stackrel{\mathrm{‾}}{I}$, L_D, and T.

Moreover, although the model, in its current state, is not ready to be used in a spatial setup (see below), we believe that this study can provide the basis for extension of the model for spatially explicit frameworks.

4.2 Computational efficiency and application potential

Our results demonstrate that a state variable that tracks the elemental content of plankton can be effectively removed, without leading to major issues in mass balance, in a 0D setup. In comparison to a fully explicit dynamic variant, removing a state variable does not seem to result in clear advantages in computational efficiency in such a spatially truncated setup that does not require the calculation of spatial transport. Any potential reduction in computational costs owing to one less state variable in the IA variant is apparently compensated for by the additional calculations required for the derivatives in Eq. (4b.3) (see Appendix A).

The modeling framework FABM, in which the model is implemented, allows seamless coupling of models with various hydrodynamical hosts (Bruggeman and Bolding, 2014). We have also attempted an application in a 1D setup using GOTM (Burchard et al., 2006) as the hydrodynamical host and found out that N is not conserved. This is to be expected because the spatial transport calculations do not account for spatial gradients in Q, which introduces errors analogous to the difference between Eqs. (4b) and (4b.3). It may be possible to develop a mass conservative IA approach for spatially explicit models by accounting for spatial variations in Q in addition to its temporal variations. For FABM implementation, this would require additional spatial flux terms to be communicated with the hydrodynamical driver. It is not clear whether the resulting model would be faster than the fully explicit variant. On one hand, in a spatially explicit setup, reducing the number of state variables would also reduce the size of the necessary transport matrices. On the other hand, fluxes associated with spatial changes in Q would require additional logic and calculations. As mentioned above, we have not found evidence for performance advantages of the IA approach in a 0D setup. However, in spatially explicit setups, transport calculations can become computationally more demanding than calculating the right-hand sides of a biogeochemical model. Therefore, reducing the number of state variables might offer computational advantages.

5 Conclusion

Accounting for the variability in phytoplankton cellular composition is required for a realistic representation of nutrient cycling. Variable cellular composition is usually described by a dynamic acclimation (DA) approach, which requires additional state variables for the cellular constituents, thereby increasing computational costs. Smith et al. (2016) proposed the instantaneous acclimation (IA) approach, which approximates the variability in cellular composition without the need for additional state variables. As long as only one of carbon (C) or nitrogen (N) is resolved (i.e., the mass balance is closed globally), the IA approach is fully conservative and can be applied, e.g., for ecologically oriented questions (e.g., Kerimoglu et al., 2021). Here we provide a formally consistent and complete explanation of how the mass balance can be nearly re-established for when both N and C are resolved by an IA model. Through several tests in 0D setups, we demonstrate that under stable environmental conditions, the fully explicit model can be closely reproduced but that transient differences between the IA and DA variants can emerge and mass balance can be slightly compromised. A generalization of the IA approach to account also for spatial variability will require extending our (0D) IA framework towards spatially explicit setups. In our 0D setup, we did not find evidence for improved computational efficiency. However, gains in spatially explicit setups may be possible, given that the number of state variables to be transported is known to significantly affect the computational costs.

Appendix A: Analytical solutions

To facilitate the solutions of the $\partial Q/\partial E$ ($E=\mathit{\{}\mathrm{DIN},\stackrel{\mathrm{‾}}{I},\mathrm{T},{\mathrm{L}}_{\mathrm{D}}\mathit{\}}$) in Eq. (4b.3), we introduce a new variable Z and re-write Eq. (6) in terms of Z (Smith et al., 2016, S16 in the following):

$$\begin{array}{}\text{(A1)}& {\displaystyle }Q={\displaystyle \frac{Q_{\mathrm{0}}}{\mathrm{2}}}\left(\mathrm{1}+\sqrt{\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{Z}}}\right),Z={\displaystyle \frac{Q_{\mathrm{0}}}{\mathrm{2}}}\left({\displaystyle \frac{{\widehat{\mathit{\mu }}}_{\mathrm{net}}}{{\widehat{V}}_{\mathrm{N}}}}+{\mathit{\zeta }}_{\mathrm{N}}\right),\text{(A2)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial Q}{\partial \mathrm{DIN}}}& ={\displaystyle \frac{\partial Q}{\partial Z}}{\displaystyle \frac{\partial Z}{\partial \mathrm{DIN}}},{\displaystyle \frac{\partial Q}{\partial \stackrel{\mathrm{‾}}{I}}}={\displaystyle \frac{\partial Q}{\partial Z}}{\displaystyle \frac{\partial Z}{\partial \stackrel{\mathrm{‾}}{I}}},\\ & {\displaystyle \frac{\partial Q}{\partial L_{\mathrm{D}}}}={\displaystyle \frac{\partial Q}{\partial Z}}{\displaystyle \frac{\partial Z}{\partial L_{\mathrm{D}}}},{\displaystyle \frac{\partial Q}{\partial T}}={\displaystyle \frac{\partial Q}{\partial Z}}{\displaystyle \frac{\partial Z}{\partial T}}.\end{array}\end{array}$$

In Eq. (A2), the common term $\partial Q/\partial Z$, as in S16, is

$$\begin{array}{}\text{(A3)}& {\displaystyle \frac{\partial Q}{\partial Z}}={\displaystyle \frac{-Q_{\mathrm{0}}}{\mathrm{4}\cdot Z\cdot \sqrt{Z\cdot (\mathrm{1}+Z)}}}.\end{array}$$

Recalling ${\widehat{V}}_{\mathrm{N}}$ from K21, Eq. (17),

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}{\widehat{V}}_{\mathrm{N}}& ={\displaystyle \frac{(\mathrm{1}-f_{\mathrm{A}}){\widehat{V}}_{\mathrm{0}}{f}_{\mathrm{A}}{\widehat{A}}_{\mathrm{0}}\mathrm{DIN}}{(\mathrm{1}-f_{\mathrm{A}}){\widehat{V}}_{\mathrm{0}}+f_{\mathrm{A}}{\widehat{A}}_{\mathrm{0}}\mathrm{DIN}}}={\displaystyle \frac{{\widehat{V}}_{\mathrm{0}}\cdot {\widehat{A}}_{\mathrm{0}}\cdot \mathrm{DIN}}{{\left(\sqrt{{\widehat{V}}_{\mathrm{0}}}+\sqrt{{\widehat{A}}_{\mathrm{0}}\cdot \mathrm{DIN}}\right)}^{\mathrm{2}}}},\\ & f_{\mathrm{A}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\sqrt{\frac{{\widehat{A}}_{\mathrm{0}}\cdot \mathrm{DIN}}{{\widehat{V}}_{\mathrm{0}}}}}}.\end{array}\end{array}$$

We set the potential maximum rates of N and C acquisition numerically equal to the maximum-rate parameter μ₀ (Pahlow et al., 2013):

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}{\widehat{V}}_{\mathrm{0}}& ={\widehat{\mathit{\mu }}}_{\mathrm{0}}={\mathit{\mu }}_{\mathrm{0}}\cdot f\left(T\right),\\ & f\left(T\right)=\exp \left[-{\displaystyle \frac{E_{\mathrm{a}}}{R}}\left({\displaystyle \frac{\mathrm{1}}{T/\mathrm{K}}}-{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{ref}}/\mathrm{K}}}\right)\right].\end{array}\end{array}$$

The partial derivative of Z with respect to DIN is

$$\begin{array}{}\text{(A6)}& {\displaystyle \frac{\partial Z}{\partial \mathrm{DIN}}}={\displaystyle \frac{\partial Z}{\partial {\widehat{V}}_{\mathrm{N}}}}{\displaystyle \frac{\mathrm{d}{\widehat{V}}_{\mathrm{N}}}{\mathrm{d}\mathrm{DIN}}}=-{\displaystyle \frac{{\widehat{\mathit{\mu }}}_{\mathrm{net}}\cdot Q_{\mathrm{0}}}{\mathrm{2}\cdot {\widehat{A}}_{\mathrm{0}}\cdot {\mathrm{DIN}}^{\mathrm{2}}}}\left(\mathrm{1}+\sqrt{{\displaystyle \frac{{\widehat{A}}_{\mathrm{0}}\cdot \mathrm{DIN}}{{\widehat{V}}_{\mathrm{0}}}}}\right).\end{array}$$

To calculate the partial derivative of Z with respect to $\stackrel{\mathrm{‾}}{I}$, $\partial Z/\partial \stackrel{\mathrm{‾}}{I}$, we recall ${\widehat{\mathit{\mu }}}_{\mathrm{net}}$, L_I and $\widehat{\mathit{\theta }}$ from K21, Eqs. (20)–(23) and (26),

$$\begin{array}{}\text{(A7)}& {\displaystyle }\begin{array}{rl}{\widehat{\mathit{\mu }}}_{\mathrm{net}}& =L_{\mathrm{D}}{\widehat{\mathit{\mu }}}_{\mathrm{0}}{L}_{\mathrm{I}}(\mathrm{1}-{\mathit{\zeta }}_{\mathrm{Chl}}\widehat{\mathit{\theta }})-R^{\mathrm{Chl}},\\ & R^{\mathrm{Chl}}=f\left(T\right)\cdot R_{\mathrm{M}}^{\mathrm{Chl}}{\mathit{\zeta }}_{\mathrm{Chl}}\widehat{\mathit{\theta }}\end{array}\text{(A8)}& {\displaystyle }\begin{array}{rl}{L}_{\mathrm{I}}& =\mathrm{1}-\exp \left({\displaystyle \frac{-\mathit{\alpha }\widehat{\mathit{\theta }}\stackrel{\mathrm{‾}}{I}}{{\widehat{\mathit{\mu }}}_{\mathrm{0}}}}\right),\widehat{\mathit{\theta }}={\displaystyle \frac{\mathrm{1}}{{\mathit{\zeta }}_{\mathrm{Chl}}}}+{\displaystyle \frac{{\widehat{\mathit{\mu }}}_{\mathrm{0}}}{\mathit{\alpha }\cdot \stackrel{\mathrm{‾}}{I}}}\cdot (\mathrm{1}-W),\\ & W={\mathrm{W}}_{\mathrm{0}}\left[\left(\mathrm{1}+{\displaystyle \frac{f\left(T\right)\cdot R_{\mathrm{M}}^{\mathrm{Chl}}}{L_{\mathrm{D}}\cdot {\widehat{\mathit{\mu }}}_{\mathrm{0}}}}\right)\cdot \exp \left(\mathrm{1}+{\displaystyle \frac{\mathit{\alpha }\cdot \stackrel{\mathrm{‾}}{I}}{{\widehat{\mathit{\mu }}}_{\mathrm{0}}\cdot {\mathit{\zeta }}_{\mathrm{Chl}}}}\right)\right],\end{array}\end{array}$$

where W₀ is the 0 branch of Lambert's W function and α and ζ_Chl are model parameters (initial Chl-specific slope of P–I curve and cost of Chl synthesis, respectively, Table 3 in K21). $\partial Z/\partial \stackrel{\mathrm{‾}}{I}$ can then be derived by canonical application of the chain rule:

$$\begin{array}{}\text{(A9)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial Z}{\partial \stackrel{\mathrm{‾}}{I}}}& ={\displaystyle \frac{\partial Z}{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}}\left({\displaystyle \frac{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}{\partial \stackrel{\mathrm{‾}}{I}}}+{\displaystyle \frac{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}{\partial \widehat{\mathit{\theta }}}}{\displaystyle \frac{\mathrm{d}\widehat{\mathit{\theta }}}{\mathrm{d}\stackrel{\mathrm{‾}}{I}}}\right)={\displaystyle \frac{\partial Z}{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}}{\displaystyle \frac{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}{\partial \stackrel{\mathrm{‾}}{I}}}\\ & (\mathrm{because}{\displaystyle \frac{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}{\partial \widehat{\mathit{\theta }}}}=\mathrm{0}\text{by definition}),\end{array}\text{(A10)}& {\displaystyle \frac{\partial Z}{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}}={\displaystyle \frac{Q_{\mathrm{0}}}{\mathrm{2}\cdot {\widehat{V}}_{\mathrm{N}}}},\text{(A11)}& {\displaystyle \frac{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}{\partial \stackrel{\mathrm{‾}}{I}}}=L_{\mathrm{D}}\cdot (\mathrm{1}-\widehat{\mathit{\theta }}\cdot {\mathit{\zeta }}_{\mathrm{Chl}})\cdot \mathit{\alpha }\cdot \widehat{\mathit{\theta }}\cdot (\mathrm{1}-L_{\mathrm{I}}).\end{array}$$

The day length derivatives are

$$\begin{array}{}\text{(A12)}& {\displaystyle \frac{\partial Z}{\partial L_{\mathrm{D}}}}={\displaystyle \frac{\partial Z}{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}}\left({\displaystyle \frac{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}{\partial L_{\mathrm{D}}}}+{\displaystyle \frac{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}{\partial \widehat{\mathit{\theta }}}}{\displaystyle \frac{\mathrm{d}\widehat{\mathit{\theta }}}{\mathrm{d}\stackrel{\mathrm{‾}}{I}}}\right)={\displaystyle \frac{\partial Z}{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}}{\displaystyle \frac{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}{\partial L_{\mathrm{D}}}},\text{(A13)}& {\displaystyle \frac{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}{\partial L_{\mathrm{D}}}}={\widehat{\mathit{\mu }}}_{\mathrm{0}}\cdot L_{\mathrm{I}}\cdot (\mathrm{1}-{\mathit{\zeta }}_{\mathrm{Chl}}\widehat{\mathit{\theta }}).\end{array}$$

The temperature derivative of Z is obtained via the derivatives with respect to ${\widehat{\mathit{\mu }}}_{\mathrm{0}}$, ${\widehat{V}}_{\mathrm{0}}$, and R^Chl:

$$\begin{array}{}\text{(A14)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial Z}{\partial T}}& ={\displaystyle \frac{\partial Z}{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}}\left({\displaystyle \frac{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}{\partial {\widehat{\mathit{\mu }}}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\widehat{\mathit{\mu }}}_{\mathrm{0}}}{\partial T}}+{\displaystyle \frac{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}{\partial R^{\mathrm{Chl}}}}{\displaystyle \frac{\partial R^{\mathrm{Chl}}}{\partial T}}\right)\\ & +{\displaystyle \frac{\partial Z}{\partial {\widehat{V}}_{\mathrm{N}}}}{\displaystyle \frac{\partial {\widehat{V}}_{\mathrm{N}}}{\partial {\widehat{V}}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\widehat{V}}_{\mathrm{0}}}{\partial T}}\\ & =\left[{\displaystyle \frac{\partial Z}{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}}\left({\displaystyle \frac{\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}}{\partial {\widehat{\mathit{\mu }}}_{\mathrm{0}}}}\cdot {\widehat{\mathit{\mu }}}_{\mathrm{0}}-R^{\mathrm{Chl}}\right)+{\displaystyle \frac{\partial Z}{\partial {\widehat{V}}_{\mathrm{0}}}}\cdot {\widehat{V}}_{\mathrm{0}}\right]\\ & {\displaystyle \frac{\mathrm{1}}{f\left(T\right)}}{\displaystyle \frac{\mathrm{d}f\left(T\right)}{\mathrm{d}T}},\end{array}\text{(A15)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial {\widehat{\mathit{\mu }}}_{net}}{\partial {\widehat{\mathit{\mu }}}_{\mathrm{0}}}}& =L_{\mathrm{D}}\cdot (\mathrm{1}-{\mathit{\zeta }}_{\mathrm{Chl}}\widehat{\mathit{\theta }})\left[L_{\mathrm{I}}-(\mathrm{1}-L_{\mathrm{I}}){\displaystyle \frac{\mathit{\alpha }\cdot \stackrel{\mathrm{‾}}{I}}{{\widehat{\mathit{\mu }}}_{\mathrm{0}}}}\widehat{\mathit{\theta }}\right],\\ & {\displaystyle \frac{\partial Z}{\partial {\widehat{V}}_{\mathrm{0}}}}=-{\widehat{\mathit{\mu }}}_{\mathrm{net}}{\displaystyle \frac{Q_{\mathrm{0}}}{\mathrm{2}{\widehat{\mathit{\mu }}}_{\mathrm{0}}\sqrt{{\widehat{\mathit{\mu }}}_{\mathrm{0}}\cdot {\widehat{V}}_{\mathrm{N}}}}},\end{array}\text{(A16)}& {\displaystyle \frac{\mathrm{1}}{f\left(T\right)}}{\displaystyle \frac{\mathrm{d}f\left(T\right)}{\mathrm{d}T}}={\displaystyle \frac{E_{\mathrm{a}}}{R\cdot (T/\mathrm{K}{)}^{\mathrm{2}}}}.\end{array}$$

We would like to clarify that (1) replacement of ${\widehat{\mathit{\mu }}}_g$ (in S16, ${\widehat{\mathit{\mu }}}^I$) with ${\widehat{\mathit{\mu }}}_{\mathrm{net}}$ in our model (see K21) results in the appearance of $(\mathrm{1}-\widehat{\mathit{\theta }}\cdot {\mathit{\zeta }}_{\mathrm{Chl}})$ when computing $\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}/\partial \stackrel{\mathrm{‾}}{I}$; (2) L_D used to be implicit in S16 but not it is explicit (K21 Eq. 21) and therefore appears for $\partial {\widehat{\mathit{\mu }}}_{\mathrm{net}}/\partial \stackrel{\mathrm{‾}}{I}$ unlike in S16; (3) and changes in Q due to T were not accounted for by Smith et al. (2016).