FABMNflexPD 2.0: testing an instantaneous acclimation approach for modeling the implications of phytoplankton ecophysiology for the carbon and nutrient cycles
 ^{1}Institute for Chemistry and Biology of the Marine Environment, University of Oldenburg, Oldenburg, Germany
 ^{2}GEOMAR Helmholtz Centre for Ocean Research Kiel, Kiel, Germany
 ^{3}Earth SURFACE Research Center, Research Institute for Global Change, JAMSTEC, Yokosuka, Japan
 ^{a}present address: Dept. of Earth, Ocean and Ecological Sciences, School of Environmental Sciences, University of Liverpool, Liverpool, UK
 ^{1}Institute for Chemistry and Biology of the Marine Environment, University of Oldenburg, Oldenburg, Germany
 ^{2}GEOMAR Helmholtz Centre for Ocean Research Kiel, Kiel, Germany
 ^{3}Earth SURFACE Research Center, Research Institute for Global Change, JAMSTEC, Yokosuka, Japan
 ^{a}present address: Dept. of Earth, Ocean and Ecological Sciences, School of Environmental Sciences, University of Liverpool, Liverpool, UK
Correspondence: Onur Kerimoglu (kerimoglu.o@gmail.com)
Hide author detailsCorrespondence: Onur Kerimoglu (kerimoglu.o@gmail.com)
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 physicalflux calculations to account for spatial differences in Q between model grid cells.
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 readjustment 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 FABMNflexPD 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.
FABMNflexPD 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 FABMNflexPD 2.0, which tracks both C and N. We present a detailed description of a Cbased 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.
FABMNflexPD 2.0 differs from FABMNflexPD 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
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
where Q is the phytoplankton N quota (N:C ratio). Dilution and sinking terms describe fluxes in and out of the system and are nonzero only for the test T4 (see below).
2.1 C and N content of phytoplankton
As mentioned above, in this Cbased 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):
where Q_{0} 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:
The second term in Eq. (1a), ${F}_{{\mathrm{Phy}}_{\mathrm{C}}{\mathrm{Det}}_{\mathrm{C}}}$, represents the mortality of phytoplankton:
Its N counterpart is found by multiplication with Q:
where m_{C} [${\mathrm{m}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}{\mathrm{mmolC}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$] is the Cbased specific mortality rate. The hydrolysis and remineralization fluxes are calculated as firstorder reactions:
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:
For the IA variant, the exact value of this flux is unknown due to the nonexistent Phy_{N} pool and the corresponding flux. By substituting Phy_{N} with Phy_{C}⋅Q and applying the product rule we get
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):
and plugging Eqs. (7) and (9) into Eq. (1b.3) yields
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 redefinition of ${F}_{\mathrm{DIN}{\mathrm{Phy}}_{\mathrm{N}}}$ for the IA variant:
which essentially redirects part of the fluxes associated with Phy_{N} to DIN. Plugging this into Eq. (4b) and recognizing that $\frac{\mathrm{d}\phantom{\rule{0.125em}{0ex}}Q}{\mathrm{d}\phantom{\rule{0.125em}{0ex}}t}$ consists of partial derivatives with respect to DIN, daily average irradiance, $\stackrel{\mathrm{\u203e}}{I}$, fractional day length, L_{D}, and temperature, T (see Appendix A),
and reorganizing the $\frac{\mathrm{d}\phantom{\rule{0.125em}{0ex}}\mathrm{DIN}}{\mathrm{d}\phantom{\rule{0.125em}{0ex}}t}$ on the righthand side
The partial derivatives of Q with respect to DIN, $\stackrel{\mathrm{\u203e}}{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{\u203e}}{I}$, L_{D}, and T over time, i.e., $\mathrm{d}\phantom{\rule{0.125em}{0ex}}\stackrel{\mathrm{\u203e}}{I}/\mathrm{d}\phantom{\rule{0.125em}{0ex}}t$, $\mathrm{d}\phantom{\rule{0.125em}{0ex}}{L}_{\mathrm{D}}/\mathrm{d}\phantom{\rule{0.125em}{0ex}}t$, and $\mathrm{d}\phantom{\rule{0.125em}{0ex}}T/\mathrm{d}\phantom{\rule{0.125em}{0ex}}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}\phantom{\rule{0.125em}{0ex}}E/\mathrm{d}\phantom{\rule{0.125em}{0ex}}t\approx ({E}_{i}{E}_{i\mathrm{1}})/\mathrm{\Delta}t$ for $E=\mathit{\{}\stackrel{\mathrm{\u203e}}{I},\phantom{\rule{0.33em}{0ex}}{L}_{\mathrm{D}},\phantom{\rule{0.33em}{0ex}}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:
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 circulardependency error in the current implementation, an intermediate module collects the necessary terms from the two modules and sets the righthand sides for DIN at once.
2.3 Test setups and model operation
For all tests, the model is operated in a spatially homogeneous onebox setup, using the 0D driver of FABM (Bruggeman and Bolding, 2014). With this 0D setup, numerical solutions are obtained using a fourthorder 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{\u203e}}{I}$) is provided (as implemented in K21) as a sinusoidal function of day of year (t) to represent a seasonal cycle typical of a highlatitude environment in the Northern Hemisphere:
where ${\stackrel{\mathrm{\u203e}}{I}}_{min}=\mathrm{1.6}\phantom{\rule{0.125em}{0ex}}\mathrm{mol}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$ and ${\stackrel{\mathrm{\u203e}}{I}}_{max}=\mathrm{110}\phantom{\rule{0.125em}{0ex}}\mathrm{mol}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\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{\u203e}}{I}$ is calculated numerically as a finitedifference approximation:
$$\begin{array}{}\text{(14a)}& {\displaystyle \frac{\mathrm{d}\phantom{\rule{0.125em}{0ex}}\stackrel{\mathrm{\u203e}}{I}}{\mathrm{d}\phantom{\rule{0.125em}{0ex}}t}}\approx {\displaystyle \frac{{\stackrel{\mathrm{\u203e}}{I}}_{i}{\stackrel{\mathrm{\u203e}}{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}\phantom{\rule{0.125em}{0ex}}\stackrel{\mathrm{\u203e}}{I}}{\mathrm{d}\phantom{\rule{0.125em}{0ex}}t}}=({\stackrel{\mathrm{\u203e}}{I}}_{max}{\stackrel{\mathrm{\u203e}}{I}}_{min}){\displaystyle \frac{\mathit{\pi}}{\mathrm{365}}}\mathrm{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).
2.3.2 T2 – T3
For these tests, we apply the numerical and analytical time derivatives for $\stackrel{\mathrm{\u203e}}{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_{0}, ${\widehat{V}}_{\mathrm{0}}$, and ${\widehat{A}}_{\mathrm{0}}$ vary according to allometric relationships (Table 1). Scalings of Q_{0} and ${\widehat{V}}_{\mathrm{0}}$ are based on a combination of cellspecific scalings of subsistence quotas (Edwards et al., 2012, “marine species”), maximum uptake rates (Marañón et al., 2013), and cellspecific C content (MendenDeuer 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 nonzero 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:
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:
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_{0}, ${\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.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 2fold 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.
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{\u203e}}{I}$ and $\mathrm{d}\stackrel{\mathrm{\u203e}}{I}/\mathrm{d}t$.
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 reallocation 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.
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 nearidentical. 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).
3.4 T4: comparing DA and IA variants in a nonclosed 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.
In this study, we present FABMNflexPD 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, FABMNflexPD 1.0 (K21, Kerimoglu et al., 2021), only resolves the N cycle and does not close the C cycle. FABMNflexPD 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 Reestablishing 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ándezCastro 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 reestablish 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{\u203e}}{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 nonnegligible (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.
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 righthand sides of a biogeochemical model. Therefore, reducing the number of state variables might offer computational advantages.
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 reestablished 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.
To facilitate the solutions of the $\partial Q/\partial E$ ($E=\mathit{\{}\mathrm{DIN},\stackrel{\mathrm{\u203e}}{I},\mathrm{T},{\mathrm{L}}_{\mathrm{D}}\mathit{\}}$) in Eq. (4b.3), we introduce a new variable Z and rewrite Eq. (6) in terms of Z (Smith et al., 2016, S16 in the following):
In Eq. (A2), the common term $\partial Q/\partial Z$, as in S16, is
Recalling ${\widehat{V}}_{\mathrm{N}}$ from K21, Eq. (17),
We set the potential maximum rates of N and C acquisition numerically equal to the maximumrate parameter μ_{0} (Pahlow et al., 2013):
The partial derivative of Z with respect to DIN is
To calculate the partial derivative of Z with respect to $\stackrel{\mathrm{\u203e}}{I}$, $\partial Z/\partial \stackrel{\mathrm{\u203e}}{I}$, we recall ${\widehat{\mathit{\mu}}}_{\mathrm{net}}$, L_{I} and $\widehat{\mathit{\theta}}$ from K21, Eqs. (20)–(23) and (26),
where W_{0} is the 0 branch of Lambert's W function and α and ζ_{Chl} are model parameters (initial Chlspecific slope of P–I curve and cost of Chl synthesis, respectively, Table 3 in K21). $\partial Z/\partial \stackrel{\mathrm{\u203e}}{I}$ can then be derived by canonical application of the chain rule:
The day length derivatives are
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}:
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{\u203e}}{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{\u203e}}{I}$ unlike in S16; (3) and changes in Q due to T were not accounted for by Smith et al. (2016).
To run the model and reproduce the results presented in this study, FABM (submodule version that matches GOTM v6.0.0) needs to be installed with its 0D driver as the “host” (see https://github.com/fabmmodel/fabm/wiki/Buildingandinstalling, last access: 25 December 2022). The version of the FABMNflexPD used in this paper has been stored in a Zenodo repository, accessible under https://doi.org/10.5281/zenodo.6600755 (Kerimoglu, 2022). Instructions for compiling FABMNflexPD for GOTMFABM and our 0D setup are provided in README.md. The src
folder contains the Fortran code. The model was implemented as two separate modules: the phy.F90
module that describes phytoplankton growth and the abio.F90
module that describes everything other than phytoplankton (see Fig. 1 in Kerimoglu et al., 2021). The phytoplankton module can reproduce the behavior of both the IA and DA variants considered in the paper by setting model parameters accordingly. The testcases
folder contains the configuration (yaml) file that was used to produce the results for T4 in this paper, which can be simplified or extended to conduct the other tests.
OK and SLS conceived and designed the study with contributions from MP; OK and MP extended the model with fluxes required to satisfy the mass balance; OK implemented the model in FABM; PA configured and performed the runs with multiple PFTs; OK drafted the paper; MP and SLS contributed to writing and revision of the text.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We acknowledge the developers of the opensource software used in this study, foremost FABM and GOTM.
This research has been supported by the Deutsche Forschungsgemeinschaft (grant no. KE 1970/21, PI: Onur Kerimoglu) and the Japan Society for the Promotion of Science (grant no. PI: Sherwood Lan Smith).
This paper was edited by Christopher Horvat and reviewed by two anonymous referees.
Anugerahanti, P., Kerimoglu, O., and Smith, S. L.: Enhancing Ocean Biogeochemical Models With Phytoplankton Variable Composition, Front. Mar. Sci., 8, 675428, https://doi.org/10.3389/fmars.2021.675428, 2021. a
Arrigo, K. R.: Marine microorganisms and global nutrient cycles, Nature, 437, 343–348, https://doi.org/10.1038/nature04158, 2005. a
Ayata, S. D., Lévy, M., Aumont, O., Sciandra, A., SainteMarie, J., Tagliabue, A., and Bernard, O.: Phytoplankton growth formulation in marine ecosystem models: Should we take into account photoacclimation and variable stoichiometry in oligotrophic areas?, J. Marine Syst., 125, 29–40, https://doi.org/10.1016/j.jmarsys.2012.12.010, 2013. a
Bonachela, J. A., Allison, S. D., Martiny, A. C., and Levin, S. A.: A model for variable phytoplankton stoichiometry based on cell protein regulation, Biogeosciences, 10, 4341–4356, https://doi.org/10.5194/bg1043412013, 2013. a
Bonachela, J. A., Klausmeier, C. A., Edwards, K. F., Litchman, E., and Levin, S. A.: The role of phytoplankton diversity in the emergent oceanic stoichiometry, J. Plankton Res., 38, 1021–1035, https://doi.org/10.1093/plankt/fbv087, 2016. a
Bruggeman, J. and Bolding, K.: A general framework for aquatic biogeochemical models, Environ. Modell. Softw., 61, 249–265, https://doi.org/10.1016/j.envsoft.2014.04.002, 2014. a, b, c, d
Burchard, H., Bolding, K., Kühn, W., Meister, A., Neumann, T., and Umlauf, L.: Description of a flexible and extendable physicalbiogeochemical model system for the water column, J. Marine Syst., 61, 180–211, https://doi.org/10.1016/j.jmarsys.2005.04.011, 2006. a
Burmaster, D. E.: The Continuous Culture of Phytoplankton: Mathematical Equivalence Among Three SteadyState Models, Am. Nat., 113, 123–134, https://doi.org/10.1086/283368, 1979. a
Burson, A., Stomp, M., Akil, L., Brussaard, C. P. D., and Huisman, J.: Unbalanced reduction of nutrient loads has created an offshore gradient from phosphorus to nitrogen limitation in the North Sea, Limnol. Oceanogr., 61, 869–888, https://doi.org/10.1002/lno.10257, 2016. a
Caperon, J.: Population growth response of Isochrysis Galbana to nitrate variation at limiting concentrations, Ecology, 49, 866–872, 1968. a
Chen, B. and Smith, S. L.: Optimalitybased approach for computationally efficient modeling of phytoplankton growth, chlorophylltocarbon, and nitrogentocarbon ratios, Ecol. Model., 385, 197–212, https://doi.org/10.1016/j.ecolmodel.2018.08.001, 2018. a
Droop, M.: Vitamin B12 and marine ecology. IV. The kinetics of uptake, growth and inhibition in Monochrysis lutheri, J. Mar. Biol. Assoc. UK, 48, 689–733, 1968. a
Dutkiewicz, S., Cermeno, P., Jahn, O., Follows, M. J., Hickman, A. E., Taniguchi, D. A. A., and Ward, B. A.: Dimensions of marine phytoplankton diversity, Biogeosciences, 17, 609–634, https://doi.org/10.5194/bg176092020, 2020. a
Edwards, K. F., Thomas, M. K., Klausmeier, C. A., and Litchman, E.: Allometric scaling and taxonomic variation in nutrient utilization traits and maximum growth rate of phytoplankton, Limnol. Oceanogr, 57, 554–566, https://doi.org/10.4319/lo.2012.57.2.0554, 2012. a
FernándezCastro, B., Pahlow, M., MouriñoCarballido, B., Marañón, E., and Oschlies, A.: Optimalitybased Trichodesmium Diazotrophy in the North Atlantic Subtropical Gyre, J. Plankton Res., 38, 946–963, https://doi.org/10.1093/plankt/fbw047, 2016. a
Flynn, K. J.: Do we need complex mechanistic photoacclimation models for phytoplankton?, Limnol. Oceanogr., 48, 2243–2249, 2003. a
Follows, M. J., Dutkiewicz, S., Grant, S., and Chisholm, S. W.: Emergent Biogeography of Microbial Communities in a Model Ocean, Science, 315, 1843–1846, https://doi.org/10.1126/science.1138544, 2007. a
Forsythe, W. C., Rykiel, E. J., Stahl, R. S., Wu, H.I., and Schoolfield, R. M.: A model comparison for daylength as a function of latitude and day of year, Ecol. Model., 80, 87–95, 1995. a
Fulton, E. A., Smith, A. D., and Johnson, C. R.: Effect of complexity on marine ecosystem models, Mar. Ecol. Prog. Ser., 253, 1–16, https://doi.org/10.3354/meps253001, 2003. a
Garcia, N. S., Bonachela, J. A., and Martiny, A. C.: Interactions between growthdependent changes in cell size, nutrient supply and cellular elemental stoichiometry of marine Synechococcus, ISME J., 10, 2715–2724, https://doi.org/10.1038/ismej.2016.50, 2016. a
Geider, R. and La Roche, J.: Redfield revisited: variability of C:N:P in marine microalgae and its biochemical basis, Eur. J. Phycol., 37, 1–17, https://doi.org/10.1017/S0967026201003456, 2002. a
Geider, R., Maclntyre, H., and Kana, T.: A dynamic regulatory model of phytoplanktonic acclimation to light, nutrients, and temperature, Limnol. Oceanogr., 43, 679–694, https://doi.org/10.4319/lo.1998.43.4.0679, 1998. a
Grover, J.: Resource competition in a variable environment: phytoplankton growing according to the variableinternalstores model, Am. Nat., 138, 811–835, 1991. a
Inomura, K., Omta, A. W., Talmy, D., Bragg, J., Deutsch, C., and Follows, M. J.: A Mechanistic Model of Macromolecular Allocation, Elemental Stoichiometry, and Growth Rate in Phytoplankton, Front. Microbiol., 11, 1–22, https://doi.org/10.3389/fmicb.2020.00086, 2020. a
Kerimoglu, O.: OnurKerimoglu/fabmnflexpd: FABMNflexPD Version 2.0 release candidate 0 (v2.0rc0), Zenodo [code], https://doi.org/10.5281/zenodo.6600755, 2022. a
Kerimoglu, O., Hofmeister, R., Maerz, J., Riethmüller, R., and Wirtz, K. W.: The acclimative biogeochemical model of the southern North Sea, Biogeosciences, 14, 4499–4531, https://doi.org/10.5194/bg1444992017, 2017. a
Kerimoglu, O., Große, F., Kreus, M., Lenhart, H.J., and van Beusekom, J. E.: A modelbased projection of historical state of a coastal ecosystem: relevance of phytoplankton stoichiometry, Sci. Total Environ., 639, 1311–1323, https://doi.org/10.1016/j.scitotenv.2018.05.215, 2018. a
Kerimoglu, O., Anugerahanti, P., and Smith, S. L.: FABMNflexPD 1.0: assessing an instantaneous acclimation approach for modeling phytoplankton growth, Geosci. Model Dev., 14, 6025–6047, https://doi.org/10.5194/gmd1460252021, 2021. a, b, c, d, e
Klausmeier, C., Litchman, E., Daufresne, T., and Levin, S.: Optimal nitrogentophosphorus stoichiometry of phytoplankton, Nature, 429, 171–174, https://doi.org/10.1038/nature02454, 2004. a
Kwiatkowski, L., Aumont, O., Bopp, L., and Ciais, P.: The Impact of Variable Phytoplankton Stoichiometry on Projections of Primary Production, Food Quality, and Carbon Uptake in the Global Ocean, Global Biogeochem. Cy., 32, 516–528, https://doi.org/10.1002/2017GB005799, 2018. a
Lenton, T. M. and Klausmeier, C. A.: Biotic stoichiometric controls on the deep ocean N:P ratio, Biogeosciences, 4, 353–367, https://doi.org/10.5194/bg43532007, 2007. a
Litchman, E., Klausmeier, C. A., and Yoshiyama, K.: Contrasting size evolution in marine and freshwater diatoms, P. Natl. Acad. Sci. USA, 106, 2665–2670, 2009. a
Marañón, E., Cermeño, P., LópezSandoval, D. C., RodriguezRamos, T., Sobrino, C., HueteOrtega, M., Blanco, J., and Rodriguez, J.: Unimodal size scaling of phytoplankton growth and the size dependence of nutrient uptake and use, Ecol. Lett., 16, 371–379, https://doi.org/10.1111/ele.12052, 2013. a
MendenDeuer, S. and Lessard, E.: Carbon to volume relationships for dinoflagellates, diatoms, and other protist plankton, Limnol. Oceanogr., 45, 569–579, 2000. a
Moreno, A. R. and Martiny, A. C.: Ecological Stoichiometry of Ocean Plankton, Annu. Rev. Mar. Sci., 10, 43–69, https://doi.org/10.1146/annurevmarine121916063126, 2018. a
Pahlow, M., Dietze, H., and Oschlies, A.: Optimalitybased model of phytoplankton growth and diazotrophy, Mar. Ecol. Prog. Ser., 489, 1–16, https://doi.org/10.3354/meps10449, 2013. a, b
Pahlow, M., Chien, C.T., Arteaga, L. A., and Oschlies, A.: Optimalitybased nonRedfield plankton–ecosystem model (OPEM v1.1) in UVicESCM 2.9 – Part 1: Implementation and model behaviour, Geosci. Model Dev., 13, 4663–4690, https://doi.org/10.5194/gmd1346632020, 2020. a
Redfield, A.: On the proportions of organic derivatives in sea water and their relation to the composition of plankton, in: James Johnstone Memorial Volume, edited by: Daniel, R., 177–192, University Press of Liverpool, 1934. a
Redfield, A.: The biological control of chemical factors in the environment, Am. Sci., 46, 205–221, 1958. a
Smith, S. L., Merico, A., Hohn, S., and Brandt, G.: Sizingup nutrient uptake kinetics: combining a physiological tradeoff with sizescaling of phytoplankton traits, Mar. Ecol. Prog. Ser., 511, 33–39, https://doi.org/10.3354/meps10903, 2014. a
Smith, S. L., Pahlow, M., Merico, A., AcevedoTrejos, E., Sasai, Y., Yoshikawa, C., Sasaoka, K., Fujiki, T., Matsumoto, K., and Honda, M. C.: Flexible phytoplankton functional type (FlexPFT) model: sizescaling of traits and optimal growth, J. Plankton Res., 38, 977–992, https://doi.org/10.1093/plankt/fbv038, 2016. a, b, c, d, e, f, g, h
Taherzadeh, N., Kerimoglu, O., and Wirtz, K. W.: Can we predict phytoplankton community size structure using size scalings of ecophysiological traits?, Ecol. Model., 360, 279–289, https://doi.org/10.1016/j.ecolmodel.2017.07.008, 2017. a
Ward, B. A.: Assessing an efficient “Instant Acclimation” approximation of dynamic phytoplankton stoichiometry, J. Plankton Res., 39, 803–814, https://doi.org/10.1093/plankt/fbx040, 2017. a, b, c
Wirtz, K. W. and Kerimoglu, O.: Autotrophic Stoichiometry Emerging from Optimality and Variable Colimitation, Frontiers in Ecology and Evolution, 4, 131, https://doi.org/10.3389/fevo.2016.00131, 2016. a
instantaneous acclimationapproach, in which the elemental composition of the phytoplankton is assumed to adjust to an optimal value instantaneously. Through rigorous tests, we evaluate the consistency of this scheme.