Articles | Volume 16, issue 1
Development and technical paper
04 Jan 2023
Development and technical paper |  | 04 Jan 2023

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

Onur Kerimoglu, Markus Pahlow, Prima Anugerahanti, and Sherwood Lan Smith

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 Martiny2018). 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 Roche2002), 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; Arrigo2005; Burson et al.2016).

The potential relevance of such variability in the cellular composition of phytoplankton for biogeochemical cycles was recognized decades ago (Redfield1934, 1958), and evidence has been building ever since (Lenton and Klausmeier2007; 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 Smith2018). 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; Flynn2003) but possibly even more (Bonachela et al.2013; Wirtz and Kerimoglu2016; 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 (Caperon1968; Droop1968) 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 Bolding2014), 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:


– 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;


– 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;


– simulation with multiple phytoplankton groups;


– 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





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

(5) Phy N = Phy C Q , { IA }

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 (PhyC) is dynamically tracked via Eq. (1a), whereas PhyN 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, V^N, and net photosynthesis in the chloroplast, μ^net (Eq. 10 in K21):

(6) Q = Q 0 2 1 + 1 + 2 Q 0 ( μ ^ net / V ^ N + ζ N ) , { IA }

where Q0 and ζN are the subsistence N quota and cost of N uptake, respectively (Table 1). The first term in Eq. (1a), FDIC-PhyC, represents net phytoplankton growth:

(7) F DIC - Phy C = μ Phy C .

The second term in Eq. (1a), FPhyC-DetC, represents the mortality of phytoplankton:

(8) F Phy C - Det C = m C Phy C 2 .

Its N counterpart is found by multiplication with Q:

(9) F Phy N - Det N = F Phy C - Det C Q ,

where mC [m3mmolC-1d-1] is the C-based specific mortality rate. The hydrolysis and remineralization fluxes are calculated as first-order reactions:


2.2 N fluxes between DIN and phytoplankton

For tracking the PhyN for the DA variant (Eq. 1b) and the dissolved organic nitrogen, DIN, for both variants (Eq. 4b), the flux from DIN to PhyN 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:

(12) F DIN - Phy N = V N Phy C . { DA }

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

(1b.2) d Phy N d t = d ( Phy C Q ) d t = d Phy C d t Q + Phy C d Q d t ;

here, the first term reflects the N equivalent of the change in PhyC, 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):

(1b.3) d Phy N d t = F DIC - Phy C Q - F Phy C - Det C Q - D Phy C Q + Phy C d Q d t .

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

(1b.4) d Phy N d t = μ Phy C Q - F Phy N - Det N - D Phy C Q + Phy C d Q d t .

In Eq. (1b.4), μPhyCQ can be identified with FDIN-PhyN in the IA variant because VN=μQ for balanced growth. Following Smith et al. (2016), we assign the last term to FDIN-PhyN too, yielding a re-definition of FDIN-PhyN for the IA variant:

(12.2) F DIN - Phy N = Phy C μ Q + d Q d t , { IA }

which essentially redirects part of the fluxes associated with PhyN to DIN. Plugging this into Eq. (4b) and recognizing that dQdt consists of partial derivatives with respect to DIN, daily average irradiance, I, fractional day length, LD, and temperature, T (see Appendix A),

(4b.2) d DIN d t = F DON - DIN - Phy C μ Q + d Q d t - D ( DIN - DIN in ) = F DON - DIN - Phy C μ Q + Q DIN d DIN d t + Q I d I d t + Q L D d L D d t + Q T d T d t - D ( DIN - DIN in ) , { IA }

and reorganizing the dDINdt on the right-hand side

(4b.3) d DIN d t = F DON - DIN - Phy C μ Q + Q I d I d t + Q L D d L D d t + Q T d T d t - D ( DIN - DIN in ) 1 + Phy C Q DIN . { IA }

The partial derivatives of Q with respect to DIN, I, LD, 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 I, LD, and T over time, i.e., dI/dt, dLD/dt, and dT/dt. 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-1) values, divided by the integration time step, i.e., dE/dt(Ei-Ei-1)/Δt for E={I,LD,T} (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:

(4b.4) d DIN d t = F DON - DIN - j Phy C j μ j Q j + Q j I d I d t + Q j L D d L D d t + Q j T d T d t - D ( DIN - DIN in ) 1 + j Phy C j Q j DIN . { IA }

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 1Descriptions, 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=CESDSC, where SC is the allometric scaling coefficient for this parameter. Values for C and SC are as in Smith et al. (2016), and other parameters as in Kerimoglu et al. (2021).

Download Print Version | Download XLSX

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 Bolding2014). 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 (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:

(13) I ( t ) = I min + ( I max - I min ) 0.5 1 + sin 2 π t , t = t 365 - 0.25 ,

where Imin=1.6molm-2d-1 and Imax=110molm-2d-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, LD, is unity, and we ignore light attenuation.


– in the first case, the time derivative of I is calculated numerically as a finite-difference approximation:

(14a) d I d t I i - I i - 1 Δ t ,

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


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

(14b) d I d t = ( I max - I min ) π 365 cos 2 π t .

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

Figure 1Daily 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 I (Eqs. 14a and 14b) and also for variable day length, LD, and temperature, T. LD, 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 Tmin=2 and Tmax=20C.

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), Q0, V^0, and A^0 vary according to allometric relationships (Table 1). Scalings of Q0 and V^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 Lessard2000, “protist plankton excluding diatoms”). Scaling of A^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., DIN=TotalN-DON-DetN-QPhyC). 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 wDet>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 Dmax=1.0 during winter and Dmin=0.001 during summer:

(15) D = D min + 0.25 ( D max - D min ) 1 - sin ( 2 π t ) 2 , t = t 365 - 0.15 ,

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, ExtC and ExtN, which trace the amounts of N and C exported from and imported into the system:

(16) d Ext X d t = D Phy X + Det X + DO X + ( DI X - DI X in ) + w Det H SML Det X , X { C , N } ,

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 Q0, V^0, and A^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 2T1: 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 LD (Fig. 3). Note that variations in LD also affect the seasonal cycles of I and dI/dt.

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


The IA and DA variants produce almost identical results (Figs. 45). This similarity is expected (Ward2017). On closer inspection, some differences can be detected, such as slightly higher PhyN at the peak of the spring bloom and slightly higher DetN shortly after the spring bloom in the DA variant. The differences are due to the re-allocation of part of the fluxes between PhyN 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 4T2: 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 5T2: 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 Grover1991; Litchman et al.2009).

Figure 6T3: (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 7T4: annual variations in global C and N (totalC+ExtC and totalN+ExtN, see Eq. 16), total C and N (PhyCj+DIC+DOC+DetC and PhyNj+DIN+DON+DetN), 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 (Burmaster1979), that is, assuming that cells reach the equilibrium state instantaneously.

Owing to the lack of a state variable PhyN, 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 PhyN. It is impossible to maintain mass balance mechanistically without adding a state variable PhyN. 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 PhyN to DIN, the resulting differences compared to explicitly resolving PhyN are relatively small, as was also shown previously (Ward2017).

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 Bolding2014).

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., dQ/dt, which in turn requires the calculation of individual components of this change driven by different environmental factors, namely, DIN, I, T, and LD. 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 8T2: total dQ/dt and its components as contributed by the changes in DIN, I, LD, 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 Bolding2014). 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 Q/E (E={DIN,I,T,LD}) 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):


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

(A3) Q Z = - Q 0 4 Z Z ( 1 + Z ) .

Recalling V^N from K21, Eq. (17),

(A4) V ^ N = ( 1 - f A ) V ^ 0 f A A ^ 0 DIN ( 1 - f A ) V ^ 0 + f A A ^ 0 DIN = V ^ 0 A ^ 0 DIN V ^ 0 + A ^ 0 DIN 2 , f A = 1 1 + A ^ 0 DIN V ^ 0 .

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

(A5) V ^ 0 = μ ^ 0 = μ 0 f ( T ) , f ( T ) = exp - E a R 1 T / K - 1 T ref / K .

The partial derivative of Z with respect to DIN is

(A6) Z DIN = Z V ^ N d V ^ N d DIN = - μ ^ net Q 0 2 A ^ 0 DIN 2 1 + A ^ 0 DIN V ^ 0 .

To calculate the partial derivative of Z with respect to I, Z/I, we recall μ^net, LI and θ^ from K21, Eqs. (20)–(23) and (26),


where W0 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). Z/I can then be derived by canonical application of the chain rule:

(A9)ZI=Zμ^netμ^netI+μ^netθ^dθ^dI=Zμ^netμ^netI(becauseμ^netθ^=0by definition),(A10)Zμ^net=Q02V^N,(A11)μ^netI=LD(1-θ^ζChl)αθ^(1-LI).

The day length derivatives are


The temperature derivative of Z is obtained via the derivatives with respect to μ^0, V^0, and RChl:


We would like to clarify that (1) replacement of μ^g (in S16, μ^I) with μ^net in our model (see K21) results in the appearance of (1-θ^ζChl) when computing μ^net/I; (2) LD used to be implicit in S16 but not it is explicit (K21 Eq. 21) and therefore appears for μ^net/I unlike in S16; (3) and changes in Q due to T were not accounted for by Smith et al. (2016).

Code availability

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, last access: 25 December 2022). The version of the FABM-NflexPD used in this paper has been stored in a Zenodo repository, accessible under (Kerimoglu2022). Instructions for compiling FABM-NflexPD for GOTM-FABM and our 0D setup are provided in 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.

Author contributions

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.

Competing interests

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 open-source software used in this study, foremost FABM and GOTM.

Financial support

This research has been supported by the Deutsche Forschungsgemeinschaft (grant no. KE 1970/2-1, PI: Onur Kerimoglu) and the Japan Society for the Promotion of Science (grant no. PI: Sherwood Lan Smith).

Review statement

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,, 2021. a

Arrigo, K. R.: Marine microorganisms and global nutrient cycles, Nature, 437, 343–348,, 2005. a

Ayata, S. D., Lévy, M., Aumont, O., Sciandra, A., Sainte-Marie, J., Tagliabue, A., and Bernard, O.: Phytoplankton growth formulation in marine ecosystem models: Should we take into account photo-acclimation and variable stoichiometry in oligotrophic areas?, J. Marine Syst., 125, 29–40,, 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,, 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,, 2016. a

Bruggeman, J. and Bolding, K.: A general framework for aquatic biogeochemical models, Environ. Modell. Softw., 61, 249–265,, 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 physical-biogeochemical model system for the water column, J. Marine Syst., 61, 180–211,, 2006. a

Burmaster, D. E.: The Continuous Culture of Phytoplankton: Mathematical Equivalence Among Three Steady-State Models, Am. Nat., 113, 123–134,, 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,, 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.: Optimality-based approach for computationally efficient modeling of phytoplankton growth, chlorophyll-to-carbon, and nitrogen-to-carbon ratios, Ecol. Model., 385, 197–212,, 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,, 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,, 2012. a

Fernández-Castro, B., Pahlow, M., Mouriño-Carballido, B., Marañón, E., and Oschlies, A.: Optimality-based Trichodesmium Diazotrophy in the North Atlantic Subtropical Gyre, J. Plankton Res., 38, 946–963,, 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,, 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,, 2003. a

Garcia, N. S., Bonachela, J. A., and Martiny, A. C.: Interactions between growth-dependent changes in cell size, nutrient supply and cellular elemental stoichiometry of marine Synechococcus, ISME J., 10, 2715–2724,, 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,, 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,, 1998. a

Grover, J.: Resource competition in a variable environment: phytoplankton growing according to the variable-internal-stores 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,, 2020. a

Kerimoglu, O.: OnurKerimoglu/fabm-nflexpd: FABM-NflexPD Version 2.0 release candidate 0 (v2.0-rc0), Zenodo [code],, 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,, 2017. a

Kerimoglu, O., Große, F., Kreus, M., Lenhart, H.-J., and van Beusekom, J. E.: A model-based projection of historical state of a coastal ecosystem: relevance of phytoplankton stoichiometry, Sci. Total Environ., 639, 1311–1323,, 2018. a

Kerimoglu, O., Anugerahanti, P., and Smith, S. L.: FABM-NflexPD 1.0: assessing an instantaneous acclimation approach for modeling phytoplankton growth, Geosci. Model Dev., 14, 6025–6047,, 2021. a, b, c, d, e

Klausmeier, C., Litchman, E., Daufresne, T., and Levin, S.: Optimal nitrogen-to-phosphorus stoichiometry of phytoplankton, Nature, 429, 171–174,, 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,, 2018. a

Lenton, T. M. and Klausmeier, C. A.: Biotic stoichiometric controls on the deep ocean N:P ratio, Biogeosciences, 4, 353–367,, 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ópez-Sandoval, D. C., Rodriguez-Ramos, T., Sobrino, C., Huete-Ortega, 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,, 2013. a

Menden-Deuer, 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,, 2018. a

Pahlow, M., Dietze, H., and Oschlies, A.: Optimality-based model of phytoplankton growth and diazotrophy, Mar. Ecol. Prog. Ser., 489, 1–16,, 2013. a, b

Pahlow, M., Chien, C.-T., Arteaga, L. A., and Oschlies, A.: Optimality-based non-Redfield plankton–ecosystem model (OPEM v1.1) in UVic-ESCM 2.9 – Part 1: Implementation and model behaviour, Geosci. Model Dev., 13, 4663–4690,, 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.: Sizing-up nutrient uptake kinetics: combining a physiological trade-off with size-scaling of phytoplankton traits, Mar. Ecol. Prog. Ser., 511, 33–39,, 2014. a

Smith, S. L., Pahlow, M., Merico, A., Acevedo-Trejos, E., Sasai, Y., Yoshikawa, C., Sasaoka, K., Fujiki, T., Matsumoto, K., and Honda, M. C.: Flexible phytoplankton functional type (FlexPFT) model: size-scaling of traits and optimal growth, J. Plankton Res., 38, 977–992,, 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 eco-physiological traits?, Ecol. Model., 360, 279–289,, 2017. a

Ward, B. A.: Assessing an efficient “Instant Acclimation” approximation of dynamic phytoplankton stoichiometry, J. Plankton Res., 39, 803–814,, 2017. a, b, c

Wirtz, K. W. and Kerimoglu, O.: Autotrophic Stoichiometry Emerging from Optimality and Variable Co-limitation, Frontiers in Ecology and Evolution, 4, 131,, 2016. a

Short summary
In classical models that track the changes in the elemental composition of phytoplankton, additional state variables are required for each element resolved. In this study, we show how the behavior of such an explicit model can be approximated using an instantaneous acclimation approach, 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.