Enhancing Ocean Biogeochemical Models With Phytoplankton Variable Composition

Chlorophyll (Chl) is widely taken as a proxy for phytoplankton biomass, despite well-known variations in Chl:C:biomass ratios as an acclimative response to changing environmental conditions. For the sake of simplicity and computational efficiency, many large scale biogeochemical models ignore this flexibility, compromising their ability to capture phytoplankton dynamics. Here we evaluate modelling approaches of differing complexity for phytoplankton growth response: fixed stoichiometry, fixed stoichiometry with photoacclimation, classical variable-composition with photoacclimation, and Instantaneous Acclimation with optimal resource allocation. Model performance is evaluated against biogeochemical observations from time-series sites BATS and ALOHA, where phytoplankton composition varies substantially. We analyse the sensitivity of each model variant to the affinity parameters for light and nutrient, respectively. Models with fixed stoichiometry are more sensitive to parameter perturbations, but the inclusion of photoacclimation in the fixed-stoichiometry model generally captures Chl observations better than other variants when individually tuned for each site and when using similar parameter sets for both sites. Compared to the fixed stoichiometry model including photoacclimation, models with variable C:N ratio perform better in cross-validation experiments using model-specific parameter sets tuned for the other site; i.e., they are more portable. Compared to typical variable composition approaches, instantaneous acclimation, which requires fewer state variables, generally yields better performance but somewhat lower portability than the fully dynamic variant. Further assessments using objective optimisation and more contrasting stations are suggested.


INTRODUCTION
Although phytoplankton largely drive the oceanic carbon cycle, including the export of carbon (C) from the surface to depth, direct observations of their carbon biomass are rare. The most widely observed metric of phytoplankton is chlorophyll (Chl), because of its distinctive optical properties (Macintyre et al., 2000), but Chl:C:nutrient ratios vary widely in response to fluctuations in ambient light, temperature, and nutrient levels (Geider and La Roche, 2002;Mongin et al., 2006;Martiny et al., 2013;Jakobsen and Markager, 2016). Inaccurate estimates of phytoplankton C biomass limit our ability to quantify C export, and therefore our understanding of how climate change is affecting marine ecosystems (Polovina et al., 2008;Arteaga et al., 2016). Nonetheless, for the sake of computational efficiency and simplicity, many large scale ocean biogechemical (OBGC) models use a simplistic (Monod, 1949) type formulation for phytoplankton growth, assuming constant composition, especially C:nutrient ratio (Le Quere et al., 2005;Bopp et al., 2013;Totterdell, 2019). Various approaches have been developed for modelling internal variable phytoplankton composition, from relatively simple elemental quotas (Droop, 1968), to resource allocation among a few different physiological processes (Smith et al., 2011(Smith et al., , 2016, and more detailed macromolecular allocation schemes (Inomura et al., 2020). In regions such as the Mediterranean Sea or the North Atlantic, constant composition models are able to capture the biogeochemical observations (Faugeras et al., 2003;Ward et al., 2013), but lack of variable composition hinders the models' ability to consistently reproduce various biological as well as chemical observations in oligotrophic regions (Steinberg et al., 2001;Mongin et al., 2006;Ayata et al., 2013).
Compared to models assuming constant composition, those accounting for acclimation processes (individual-level physiologic response and associated variations in Chl:C:nutrient ratios) better reproduce observations at oligotrophic timeseries sites (Schartau et al., 2001;Mongin et al., 2006;Ayata et al., 2013). However, including flexible composition can be computationally expensive; the dynamics of C and nutrients bound to phytoplankton (i.e., internal stores) are typically described using the Droop "quota" model (Caperon, 1968;Droop, 1968), which requires a separate state variable for each element or nutrient resolved (Ward, 2017;Chen and Smith, 2018b). Apart from computational cost, added complexity also increases the number of uncertain parameters (Kwiatkowski et al., 2014;Ward, 2017) and can make models less portable (Friedrichs et al., 2007). Recently, a computationally efficient "Instantaneous Acclimation" (IA) model was shown to capture phytoplankton seasonality, including variable composition, at two contrasting stations in a 0-D setup (Smith et al., 2016). Further tests in a 0-D setup suggested that IA may be suitable for application in large-scale OBGC models (Ward, 2017).
When applied in 3-D setups, particularly at the global scale as with Earth system climate modelling, OBGC models pose the formidable challenge of capturing a wide range of oceanic regions with computational efficiency and minimal parameter tuning. To assess the potential for improving the biogeochemical realism of such 3-D models, here we evaluate different formulations for the flexible composition and acclimative response of phytoplankton, in terms of model performance in a spatially explicit 1-D setup, which unlike the previous 0-D studies couples the physical processes of advection and diffusion with the biogeochemical dynamics. We apply the IA-based FlexPFT model (Smith et al., 2016) and three controls: a Droop-quota model, which captures variable composition by calculating phytoplankton C and nitrogen (N) separately (and is therefore less computationally efficient), a typical Monod-type fixed stoichiometry model, and an intermediate Monod-type model with photoacclimation. In order to fully evaluate how different aspects of acclimation may affect growth response, each phytoplankton model is incorporated into an otherwise identical Nutrient Phytoplankton Detritus (NPD) model. We compare the performance of the four model variants at two oligotrophic sites, where phytoplankton composition is known to deviate from the average "Redfield ratio" (Steinberg et al., 2001;Mongin et al., 2006;Ayata et al., 2013) typically assumed in fixed stoichiometry models. Also, zooplankton grazing rate and phytoplankton growth are tightly coupled in these oligotrophic regions (Jackson, 1980;Gutiérrez-Rodríguez et al., 2011;Jiang et al., 2021), and both vary little over the seasonal cycle (Cáceres et al., 2013). Therefore, a reasonable reproduction of the observations is possible without explicitly representing zooplankton in the models, which facilitates interpretation of our results concerning phytoplankton physiology without the complications of considering various grazing-related processes at the same time. The two sites differ in physical conditions; ALOHA (A Long term Oligotrophic Habitat Assessment, 22.75 o N, 158 o W), with permanent stratification (Kavanaugh et al., 2018) and BATS (Bermuda Atlantic Time Series, 31.67 o N, 64.167 o W) with its seasonal cycle of deep mixing (Dave and Lozier, 2010). Apart from model formulation, parameter values in OBGC models can strongly affect model performance (Ward et al., 2010). To assess the response of model skill to parameter perturbations, we explore each model's response to variations of the two parameters that are most relevant to the key differences in models, photosynthesis and nutrient uptake, respectively. Previous studies have clarified how model complexity, in terms of the number of compartments (Friedrichs et al., 2007) and parameters (Ward et al., 2013), affects the applicability of models for different oceanographic sites. However, little is known about how the complexity of formulations for phytoplankton growth may affect the suitability of models for different oceanic environments. Finally, we perform cross-validation experiments to test each model's portability (i.e., the robustness of its response) between these two challenging sites.
Our goal is to evaluate and compare how different formulations affect model results and how adding phytoplankton variable composition may improve simulations of biogoechemical variables, in particular Chl, dissolved inorganic nitrogen (DIN), and primary production (PP). Specifically, we assess how light-harvesting pigments (FS vs. FSPA), variable C:N ratio (FSPA vs. DQ), or physiological acclimation (IA vs. DQ) may enhance model performance, as well as the sensitivity of different phytoplankton growth formulations to parameter perturbations, and how these affect portability.

MATERIALS AND METHODS
The IA model is based on the assumption that the physiological allocation factors, and thereby the cellular nutrient and chlorophyll contents (expressed as N:C and Chl:C ratio, respectively) are instantly optimised to achieve maximal net growth rate (Smith et al., 2016). We compare the IA model with the Monod-type Fixed Stoichiometry (hereafter FS) variant, where phytoplankton growth rate is limited by the ambient resource levels based on a prescribed set of parameters that describe a frozen acclimation state and the "Dynamic-Quota" (hereafter DQ) variant, which assumes that growth rate depends on the non-optimal but variable internal cellular stores of nutrients, similar to a classical Droop-quota model, and optimal pigment density in the chloroplast. In addition to the two wellknown phytoplankton growth models, we consider a "Fixed Stoichiometry with Photoacclimation" (FSPA), which includes photoacclimation but not variable internal N:C ratio. This formulation is included in order to isolate and understand the effects of photoacclimation alone on model performance.
All model variants are run within the Nutrient flexible Phytoplankton Detritus (NFlexPD) model (Kerimoglu et al., 2021), written in the Framework for Aquatic Biogeochemical Models (FABM, Bruggeman and Bolding, 2014). This allows switching between different formulations to describe phytoplankton growth and uptake within a common modelling framework. Therefore, despite the different assumptions, the model variants follow a similar set of differential equations for the dynamics of the state variables: dissolved inorganic nitrogen (DIN), C and N bound to phytoplankton (Phy C and Phy N ), Detritus (Det C and Det N ), Dissolved Organic Nitrogen (DON), and Dissolved Organic Carbon (DOC). Although the NFlexPD model allows multiple phytoplankton types and different nutrients, in this study we consider one phytoplankton type and one limiting nutrient (nitrogen) to avoid over-complication and facilitate the interpretation of results. The model framework is illustrated in Figure 1 and the differential equation set describing the fluxes between each pool can be found in Kerimoglu et al. (2021). In the next subsection, we provide the description of phytoplankton growth formulations for each model variant.

Model Variants
The four model variants follow a common formulation to describe phytoplankton net growth rate (µ): Where µ g , R Chl and R N are the cellular gross growth rate, the costs associated with chlorophyll maintenance and synthesis, and cost for nutrient uptake, respectively. Since C-fixation, which determines phytoplankton growth, occurs in the chloroplast, the gross growth rateμ g and the fractional allocation of resources to the chloroplast, f C (see below), determine µ g : The gross growth rate depends on the fractional day length, L D , the potential maximum growth rate,μ 0 , and the light limited growth within the photosynthetic apparatus, L I : L I is a saturating function of daytime average irradianceĪ, the chlorophyll density in the chloroplastθ , and the light affinity α: The cost of maintaining and synthesising chlorophyll, R Chl is calculated by scaling the chloroplast-specific cost with the fractional allocation of resources to the chloroplast, f C , which determines the growth rate. Hence, similar to µ g (Equation 2): The next subsection describes how specific physiological processes that contribute to phytoplankton growth are modelled across different model variants. All equations described herein are also described in Kerimoglu et al. (2021) and Smith et al. (2016).

Instant Acclimation (IA)
Phytoplankton are able to sustain high growth rates despite low nutrient concentration, by varying their elemental composition to produce cells that contain less of the nutrients that are in short supply. Both the IA and DQ represent this flexibility. The classical way to represent this flexibility is via the Droop (1968) quota model, where gross growth rate is calculated as a function of the nutrient quota, Q: where µ ∞ , and Q 0 are the growth rate at infinite cell quota, and the subsistence cell quota at which growth becomes zero, respectively. Based on the study by , assuming a fixed amount of cellular nitrogen, N, is bound in the structural material (Q s ), a fraction of cellular N in excess of Q s is allocated for nutrient uptake, as represented by f V , and the remainder is assumed to be allocated to chloroplasts for growth (see Equation 2), as represented by f C : In turn, the quota-dependent "nutrient limitation" term in brackets in Equation (6), is replaced by f C in our approach, as described by Equation (2). Increasing f V therefore decreases the gross growth rate (Equation 2), while increasing the rate of cellular nutrient uptake, V: whereV is the protoplast specific nutrient uptake. The cost associated with nutrient uptake R N (Equation 1), is the product of the prescribed cost of N uptake, ζ N and V itself: Following  the optimal cell quota, Q 0 can be calculated as: Thus, the optimal Q 0 depends on the ratio of light-to nutrientlimitation (μ net V ). By balancing growth and uptake, µQ = V, via the balanced growth equation (Burmaster, 1979), it is possible to calculate the optimal f V , f 0 V : where ζ N is the cost of assimilating DIN [mol C (mol N -1 )]. Thus, f V increases with decreasing cellular N quota, Q; i.e., more nutrient limited cells allocate more resources towards nutrient uptake,V. This process rate depends on maximum uptake ratê V max and nutrient affinityÂ: The second fractional resource allocation is towards nutrient affinity, f A . Increasing f A increases nutrient affinity at the expense of decreasing maximum uptake rate: WhereÂ 0 andV 0 are the potential maximum affinity and nutrient uptake rate, respectively. Substituting these assumptions to Equation (12),V The optimum f A maximisesV, therefore f A can be dynamically calculated: Apart from changing the Nutrient:C composition to optimise their growth, phytoplankton are also able to adjust their Chl:C ratio, i.e., photoacclimate under changing light conditions. The Chl:C ratio within the chloroplast increases under low light to enhance light-harvesting efficiency, and decreases under high light to free up resources for other uses. In the IA and DQ model, this process is based on the rate of the net carbon fixation rate within the chloroplast : WhereR Chl is the cost of light harvesting within the chloroplast, R Chl M is the cost of maintaining chlorophyll and ζ Chl is the cost of chlorophyll synthesis.
From Equation (18), reducingθ will also reduce the cost of light harvesting, but will increase the light limitation, as described in Equation (3). In the IA model,θ adjusted instantaneously to maximise theμ net : where, W 0 is the 0-branch of the Lambert-W function andĪ C is the critical daytime average irradiance level, below which it is not worthwhile to synthesise chlorophyll : The cellular Chl:C ratio is obtained by multiplyingθ by the allocation of resources towards chloroplast, f C : With the assumption of fixed f V , nutrient uptake rate for the whole cell can become unrealistically high, as DIN and hence Phy N increase, yielding unrealistically high Q. Therefore, a downregulation term is required to limit Q values. Following earlier examples (e.g., Morel, 1987;Grover, 1991) this is described as a liner function of the difference between a maximum prescribed quota, Q max and the actual Q, relative to the effective capacity, given by the difference between the Q max and subsistence quota, Q s : The nutrient uptake rate, V in Equation (8) then becomes: which also affects the respiratory cost of nutrient uptake, R N (Equation 9): Another important difference from the IA model is that the DQ variant tracks both N and C bound to phytoplankton (i.e., Phy C , Phy N ) explicitly, similar to the "DA" variant in Kerimoglu et al. (2021), therefore Q becomes: In the classical Droop approach, increasing Q to high levels yields diminishing returns in terms of growth. The DQ variant herein, with constant allocation factors as prescribed, retains the typical saturating relationship between growth rate and Q from the Droop equation, albeit somewhat flatter, as discussed in detail in the Supplementary Material. In the DQ model, Q is calculated simply as the ratio of Phy N to Phy N . Kerimoglu et al. (2021) have demonstrated that the IA approach employed in a 1D setup yields results almost identical to an otherwise identical model that explicitly resolves the dynamics of C and N biomass (their "DA" variant). Therefore, the most important difference between our IA variant and this DQ control is that the nutrient vs. carbon uptake and nutrient affinity vs. maximum nutrient uptake rate are optimised only in the former.

Fixed Stoichiometry (FS)
The FS variant resembles the typical phytoplankton growth model based on fixed Chl:C:N as commonly used in global OBGC models (i.e., models in Le Quere et al., 2005;Laufkotter et al., 2015;Totterdell, 2019), and is the same as described by Kerimoglu et al. (2021). It describes nutrient limitation with a rectangular hyperbolic nutrient uptake function of ambient nutrient concentration (Monod, 1949): Following Button (1978) and Smith et al. (2009), the half saturation value K N can be diagnosed from the solution of the IA variant, as a function ofV max (Equation 14) andÂ (Equation 13): Based on this identity, K N is calculated from the spatio-temporal biomass-weighted average of f A , and the prescribed values of parametersV 0 andÂ 0 used in the IA variant. The light limitation is described by the same saturation function used in the other variants (Equation 4), but based on prescribed chlorophyll density,θ, based on the IA run, making θ constant, as illustrated in Figure 1. In order to ensure consistency with the IA model, θ is calculated in its expanded form, based on the resource allocation towards growth (f C ) described in Equation (7): Throughout each simulation, the terms Q fixed and f V are set to constant values, which are their respective biomass-weighted means from the IA variant. The cellular nutrient uptake rate, V, is calculated based on the balanced growth equation (i.e., V = µ · Q fixed ) as in the IA variant.

Fixed Stoichiometry With Photoacclimation (FSPA)
This variant is formulated to isolate the effect of photoacclimation on model performance and response in the oligotrophic ocean. As in the FS variant, nutrient limitation is calculated using the rectangular hyperbolic function and directly determines the phytoplankton growth. However, the Chl:C ratio within the chloroplast,θ adjusts instantaneously via Equation (19)

Simulations and Model Metrics
We simulate two oligotrophic stations with extensive time-series observations: ALOHA and BATS. In order to simulate realistic conditions with a 1-D setup, the FABM-NflexPD is coupled to the General Ocean Turbulence Model (GOTM, Burchard et al., 2006). FABM acts as an interface between the biogeochemical and hydrodynamic models, which allows switching between different model formulations (i.e., versions). The hydrodynamic model solves the advection-diffusion-reaction equations, provides physical variables, such asĪ, to the biogeochemical model, and handles the model input and output. Meanwhile, FABM provides the rate of sinking/floating of biogeochemical state variables, which are solved as residual vertical advection (Bruggeman and Bolding, 2014). For the vertical advection scheme we use the third-order TVD with ULTIMATE QUICKEST limiter, and the Runge-Kutta method for the time integration scheme with no flux boundary condition. FABM also provides feedback to physics, such as the light absorption, which is applied in this study.
As the initial conditions for the hydrodynamical model, we use in situ temperature and salinity profiles (obtained from https://hahana.soest.hawaii.edu/hot/hot-dogs/cextraction.html and batsftp.bios.edu/BATS/bottle/ for ALOHA and BATS, respectively). Meteorological forcing, from the European Centre for Medium-Range Weather Forecasts (ECMWF), ERA-5 hourly reanalysis, with horizontal resolution of 0.25 o ×0.25 o (https:// cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5single-levels?tab=overview), include: wind speed, air pressure, air temperature, humidity, cloud cover, shortwave radiation, and precipitation, each calculated using the Fairall et al. (1996) method. The model domain, split into 100 levels with surface zooming, extends to 500 m for ALOHA and 450 m for BATS, where some temperature and salinity profiles are limited to 450 m. To describe the background turbidity, we assume Jerlov type IA for ALOHA, based on field measurements in the North Pacific (Paulson and Simpson, 1977), and Jerlov type I for the very clear water at BATS (Kullenberg, 1984).
All model variants are run for 8 years, from 1st January 2008 to 31st December 2016. The first 3 years, as spinup period, are forced using repeating climatology, i.e., mean hourly and daily values for the meteorology, and the monthly means of observed temperature and salinity from 1st January 2008 to 31st December 2016. For the rest of the year, none of the forcings are nudged towards the observation (analytical method in GOTM = 1). The last 5 years of model output is compared with observations of Chl from HPLC, DIN (Nitrate + Nitrite), and PP, which are the most widely observed variables that are commonly compared to model outputs (Mongin et al., 2006;Kwiatkowski et al., 2014;Chen and Smith, 2018a). Of the widely available observables, these three are also the most directly related to the biogeochemical impacts of variable phytoplankton composition. These observations can be obtained from http://batsftp.bios.edu/ BATS/bottle/ and https://hahana.soest.hawaii.edu/hot/hot-dogs/ bextraction.html, for BATS and ALOHA, respectively. Observed Chl, DIN, and PP are used to tune the model variants.
In order to thoroughly assess how well each model variant captures the observations, and its portability (i.e., applicability to different sites without re-tuning parameter values), we perform three experiments: (i) the reference simulations of each station with individually tuned parameter sets (model runs labelled IA, FS, FSPA, and DQ); (ii) cross-validation, where phytoplanktonrelated parameters tuned for station ALOHA are applied at station BATS and vice versa, (model runs labelled IA-X, FS-X, FSPA-X and DQ-X); and (iii) simultaneous runs using a common parameter set for phytoplankton-related processes at both stations, (labelled IA-S, FS-S, FSPA-S, and DQ-S). For experiments (ii) and (iii) initial conditions and some abiotic parameters, such as sinking speed and detrital degradation rate, are kept the same as experiment (i) for each station. The choice of parameters for these three experiements are described in section 2.3. We quantify model performance in terms of statistical metrics (i.e., correlation, bias, and RMSE) and a weighted cost function for overall mismatch between models and observations.
The cost function is computed as the weighted average of the squared difference between observed and simulated values, using the sample standard deviation of all observations of type m,σ m , to appropriately weight the contribution of each model-data misfit, as in previous studies (e.g., Friedrichs et al., 2006;Hemmings and Challenor, 2012;Kaufman et al., 2018). However, because the variables that we compared tend to be distributed log normally, especially at depths where their concentration are high, we use a square-root transformation to calculate the model data difference (Dadou et al., 2004;Hemmings and Challenor, 2012). Thus, the total cost function becomes: Where M is the number of variable types (M= 3; Chl, DIN, and PP), N m is the number of observations of type m, and a jm andâ jm are the modelled and observed values, respectively. The term C m in Equation (29) is included to enhance the weight of DIN compared to the cumulative weight of the other two phytoplankton-related variables that we evaluate, by setting C m = 1 for Chl and PP, and C m = 2 for DIN. All variables are evaluated at their observed time and depth. We quantify model portability in terms of the portability index (PI, Friedrichs et al., 2007) which is the ratio of total costs (Equation 29b) from simultaneous experiment (J s ) and the cross-validation (J x ) experiment: PI values approaching unity indicate increasing portability.

Parameter Fitting
As explained in the introduction, this study explores how phytoplankton growth formulation affects model performance. Based on initial model runs, we adjust three parameters separately for each stations:μ 0 , k det , and k DON . We choose a slightly lower value forμ 0 at station ALOHA, since Chl concentration at ALOHA is lower than BATS (Saba et al., 2010) and this yield better results in terms of statistical metrics compared to using the sameμ 0 . A higher k det and lower k DON are also applied in this station to prevent the depletion of DIN near the surface during the simulations. Mortality, M p , and the fraction thereof that becomes detritus, M part , as well as other abiotic parameters are kept constant for all model variants at both stations. These parameters are described on Table 1.  , Smith et al. (2016), and Robson et al. (2018), apart from w det (Yool et al., 2013). Mortality rates are higher than literature values to account implicitly for predation.
In order to ensure a fair comparison, we systematically tune the two parameters, separately for each variant, that most directly affect phytoplankton growth, α I andÂ 0 , respectively such that each model variant attains its best performance at each station. This is achieved through an analysis of the response surface of model skill for each variant, respectively. This analysis also reveals how perturbing the parameters affects the model performance, in terms of simulating the biogeochemical variables (Chl, DIN, and PP), as well as the overall model-data agreement, as quantified by the cost, for each station. The parameter ranges were made as wide as practically possible, up to the limits for which solutions can be obtained. Thus, the parameter values explored for α I space, are 2.0, 2.5, 3.0, 3.5, 4.0, 4.2, andÂ 0 space are 0.07, 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, and 1.3, resulting in 48 model simulations per variant per site. The combination of parameter values giving the lowest overall cost for each station is used in the reference simulation and cross validation experiments (running the model at BATS using the parameter set that has the lowest cost for ALOHA, and vice versa). The combination which give the lowest total cost for both stations (J BATS + J ALOHA ) is used for the simultaneous experiment (see section 2.2). The chosen parameters for these experiments are summarised on Table 2. To better understand the effect of parameter variations, we plotted the landscapes of component costs (for each observable) and the overall cost obtained over the range of parameter values.
To summarise, we perform three experiments: individual tuning, cross validation, and simultaneous using four different model variants, IA, FS, FSPA, and DQ models. Short descriptions of these experiments can be found on Table 3.

Model Skill Response Surfaces
Previous studies have shown that optimising fewer than 10 parameters tends to yield better predictive skill for models (Friedrichs et al., 2006;Ward et al., 2013). In this study, we choose to vary the two most important parameters that determine the overall phytoplankton growth response under oligotrophic conditions: the light and nutrient affinities,Â 0 and α I , respectively. Apart from identifying the most suitable parameter combinations to be applied at different stations for each model variant, exploring this parameter space allows us to investigate how different phytoplankton acclimation The values for α I andÂ 0 are chosen to give the lowest overall cost during the parameter exploration experiment (see text). Other parameters are taken from , Smith et al. (2016), and Robson et al. (2018). Parameter ζ Chl is set to a common vallue of 1, because this allowed solutions to be obtained with all photo-acclimating variants. formulations respond to parameter perturbations. Model skill response surfaces (Figure 2) suggest two distinct regions: where all model variants show relatively low costs (Â 0 > 5.0 and α I > 2.5), and where the cost increases steeply, especially at loŵ A 0 . Variants lacking variable Q especially suffer in this cost region. Despite the FSPA variant allowing photoacclimation, the cost associated with lowÂ 0 and low α I is similar to that in FS model. However, adding photoacclimation lowers the cost associated with perturbing α I , especially in the "low" cost region, to values even lower than for the variants that include both Q and photoacclimation (Figures 2A,C,D). Variants with variable N:C and Chl:C have flatter overall cost landscapes; when combining lowÂ 0 and low α I , the costs for both IA and DQ models are further reduced whenÂ 0 = 0.3 and α I < 4.0, although cost still increases for low α I andÂ 0 > 0.5. The most striking difference between the fixed stoichiometry and variable Q groups are the shapes of their cost landscapes for both non-phytoplankton (DIN) and phytoplankton-related variables (Chl and PP, Figure 3). In terms of DIN, whenÂ 0 < 0.7, models with fixed stoichiometry show a dramatic increase in DIN cost, and for the variants that include variable Q, cost increases forÂ 0 < 0.5. For Chl, these variants show an opposite pattern compared to DIN, where cost decreases with decreasingÂ 0 , with the DQ model producing higher cost and steeper landscape compared to the IA model (Figures 3E,H). Conversely, the fixed stoichiometry variants produce similar patterns for Chl and DIN costs, with a sharp increase in cost forÂ 0 < 0.1, especially in the FS variant ( Figure 3F). Inclusion of photoacclimation in the FS model reduces the Chl and PP costs, especially in the lower cost region (Â 0 > 0.5, and α I > 2,5, Figure 3G), or when bothÂ 0 and α I are low. For most of the model variants, the cost landsapes are similar for PP and Chl. However, compared to Chl, the PP cost increases less steeply with decreasing α I in the DQ model ( Figure 3L), while the IA model shows a slight increase in PP cost with increasing α I (Figure 3I). The FS and FSPA models show similar tendency in terms of the cost landscape, although the latter produces slightly lower cost. The FS and FSPA variants also show shallow declines in PP cost with decreasingÂ 0 , but sharp increases forÂ 0 > 0.1. For these variants, with nutrient affinity near the low end of its range (Â 0 = 0.07), PP cost declines sharply with increasing α I .
In summary, for the overall cost, all variants follow a similar pattern, but with less steep increases in cost as either α I and A 0 decreases in variable Q models (IA and DQ). The variable Q  models also generally produce flatter cost surfaces, both for the overall cost and its components, compared to only accounting for photoacclimation. Nevertheless, photoacclimation alone generally improves both overall and variable-specific costs.

Performance of the Three Variants
In oligotrophic regions, summertime nutrient concentrations are typically < 0.01 mmol m −3 within the euphotic zone (Steinberg et al., 2001;Anderson and Pondaven, 2003;Dave and Lozier, 2010), with vertical stratification that is destroyed by deep mixing during winter and spring (Dave and Lozier, 2010). All model variants capture these characteristics well, as seen in the DIN distributions (Figures 4F-I, 4U-X, 5I-P, 6C,D). At station BATS, winter mixing typically increases DIN from ∼ 0.05 to 0.5 mmol m −3 (Figure 4J), as captured by all model variants, but exaggerated by the IA and DQ (Figures 4F,I). At BATS all variants overestimate the average winter time DIN within the upper 50m (Figure 6C), but the IA underestimates the DIN concentration in the summer. All variants also overestimate the winter DIN concentration at ALOHA in the upper 50 m, but capture the summer-autumn concentrations ( Figure 6D). However, at 75-125 m between summer and autumn the four variants are more similar, with the DQ variant producing slightly higher DIN than other variants, especially during summertime ( Figures 5M-P, 6G,H). Although all model variants realistically capture DIN concentration and seasonality in the upper 50 m (Figures 4U-Y), none captures the sporadic spikes of DIN that occur at both stations (Figures 6G,H).
Similar to DIN, Chl in the oligotrophic region is typically scarce in the euphotic zone due to nutrient limitation. Between 50 and 200 m depth, where nutrients are more abundant, subsurface Chl maxima (SCM) occur during summer (Dave and Lozier, 2010;Mignot et al., 2014). All variants capture these characteristics qualitatively well, both seasonally and vertically at BATS (Figures 4A-E). However, at ALOHA the FS variant is unable to simulate low Chl in the upper 100 m, nor distinct summertime SCM profiles (Figures 4Q, 5R) as shown in the observation ( Figure 4T) and other variants which allow variable cell quota (Figures 4P,S). The inclusion of photoacclimation, in the FSPA model, results in better simulation of Chl profile, with low concentration in the upper 100 m and summertime SCM profiles (Figures 4C,R), despite simulating slightly deeper and lower concentration of DCM at BATS and ALOHA, respectively, compared to IA and DQ, (Figure 5C). All model variants capture typical Chl concentrations and seasonality for the SCMs that usually occur between 75 and 125 m at both stations (Figures 5A-D, 6E). However, when averaged vertically, the IA variant tends to overestimate Chl concentrations compared to other variants (Figures 5A, 6E). This result differs from Ayata et al. (2013), where all models underestimated SCM concentrations, except during blooms. Although the IA and DQ variants are qualitatively similar at BATS, generally the latter simulates slightly lower Chl and shallower SCMs at both stations, compared to the IA variant (Figures 4A,D). Comparing all model variants, from 75 to 125 m, the FS variant produces the lowest Chl because it lacks flexible quota, but the FSPA variant produces a similar range of averaged Chl concentration as IA and DQ, more similar to the former during summertime (Figures 6E,F). Despite capturing the DIN concentration in the upper 50 m, neither the DQ nor IA variant captures Chl from summer to winter at ALOHA (Figures 6B,D). The FS and FSPA perform somewhat better in this regard, but still overestimate summertime Chl at BATS (Figures 6A,B).
Unlike Chl and DIN, PP distributions in oligotrophic regions are usually skewed towards the surface and limited to the upper ∼120 m (Dave and Lozier, 2010). All variants capture the PP depth and its decline during summer at BATS, but with higher PP and later peaks than observed (Figures 4K-N, 6I). Compared to the observations and other variants, the FS variant produces briefer spikes of higher PP (>20 mgC m −3 day −1 ). The inclusion of photoacclimation prolonged the spikes of PP and deepened its summertime distribution to depths deeper than the IA and DQ models (Figures 4M,B1). When averaged horizontally, the variants with variable Q generally capture well observed PP (Figures 5Q,T,U,X) especially at BATS, but the FS and FSPA underestimates PP in the top 50 m, particularly during summer. The IA and DQ variants simulate a slight increase in PP at ∼50 m, slightly deeper than its observed depth of ∼60 m (Figures 5Q,T). Unlike BATS, PP at ALOHA peaks in the summer, and declines in spring (Dave and Lozier, 2010), with the highest rate recorded in the upper 50 m, decreasing with depth (Figures 5U-X). None of the variants completely capture these patterns, especially the FS and FSPA variants which consistently underestimate PP in the upper 50 m (Figures 6J, 5V,W), but the IA and DQ do capture winter and spring PP rates at ALOHA. The DQ variant produces generally the highest PP, which may be because it produces the highest C:N ratio (lowest Q) within the top 100 m and therefore agrees best with the observations at ALOHA (Figure 6J).
The Taylor and target diagrams (Jolliff et al., 2009) shown in Figures 7, 8 and cost values (J m , J) listed in Table 4 summarise the skills of model variants for all experiments. In general, the models perform similarly for DIN, but differ more for Chl and PP. At BATS, the statistical metrics from the latter two show some disparity in terms of standard deviations (Figures 8B,C), which are also reflected by their RMSDs in the target diagrams (Figures 8E,F). At the more stable station ALOHA, the statistical metrics differ more among the model  Table 4, the FSPA variant produces the lowest overall ( Figure 2B), Chl (Figure 3B), and DIN ( Figure 3F) costs for the individually tuned runs, as reflected in the overall cost landscape regions for highÂ 0 and α I (Figure 2C). The FSPA variant also produces the lowest Chl RMSD and bias. The absence of photoacclimation in the FS model produces the worst overall Chl, and PP costs as well as correlation and RMSD. Based on Table 4, The 'second best' model in terms of overall cost, is the IA model, which produces the lowest PP cost, RMSD, and bias, as well as the highest Chl correlation. The DQ performs slightly worse than the IA model in terms of cost (Figures 7, 8, and Table 4), however it produces the best correlation for PP, and generally performs similarly to the IA model.

Model Portability
The cross-validation experiment tests a model's ability to reproduce observations from different regions without tuning, and therefore its predictive ability (Friedrichs et al., 2006). In this experiment, all model variants differ only slightly in their statistical metrics for DIN, compared to the individual tuning. However, for Chl and PP, using the parameter set tuned for a different station generally increases cost, RMSD, and bias, while reducing correlation for all variants (Figures 7A-F, 8A-F). In the cross-validation experiment the FS and FSPA variants' statistical metrics suffer most, especially in terms of PP, where both variants greatly overestimate and underestimate the winter and seasonal values at BATS and ALOHA, respectively (Figures 6I,J). For the IA and DQ variants, the cross validation and individually tuned results differ only slightly from the observations and often capture the observations when the reference simulations do not, especially at BATS during summer within the upper 50 m (Figures 6A,B). However, at 75-125 m, the cross validations from IA often overestimate greatly the observed Chl at BATS and ALOHA (Figures 6E,F), while the DQ variant captures the observation well at BATS, but also overestimates the Chl at ALOHA. In terms of cost, the DQ variant produces the lowest discrepancy between the reference and cross validation simulation ( J= 0.0043).
The simultaneous experiment evaluates the potential applicability of each model variant at multiple stations with a common parameter set, as in a typical global or regional biogeochemical model (Friedrichs et al., 2007). The simultaneous tuning is done by choosing the parameter set that yields the  Table 4. lowest overall cost, summed over both stations (J BATS + J ALOHA ). This experiment yields relatively similar J values and statistical metrics, although the results are slightly worse in terms of J, except that for the FS and DQ variants the simultaneous runs performed better than the reference, which is most apparent for ALOHA (Figures 7B,C,E,F). As discussed in Friedrichs et al. (2007), the ideal model would have a PI ∼1 (Equation 30), indicating good performance at multiple sites without being tuned individually for each, while also having a low cost. Figures 7G,H, 8G,H reveal that in general photoacclimation greatly improves model performance, especially for Chl and DIN, and enhances portability. The inclusion of variable Q in combination with photoacclimation (IA and DQ models) yields even better performance in terms of PP and portability, at least between these oligotrophic sites. Although the fixed stoichiometry with photoacclimation (FSPA) model produces the lowest overall cost in both the reference and simultaneous experiments, its poor performance in the cross validation experiment indicates that it is less portable compared to the IA and DQ models. From Table 4 the DQ is the most portable, however in terms of the average cost over all experiments, the IA variant produces slightly lower cost (average costs for DQ and IA are 0.434 and 0.424, respectively), indicating that the DQ model is the most portable. However, compared to the DQ model, the IA model produces nominally better RMSD and bias for all observables in both cross validation and simultaneous experiments, despite having fewer state variables.

DISCUSSION
The Chl:C:nutrient ratios of phytoplankton vary as a result of both evolutionary adaptation, which has produced a variety of species and strains of different compositions, and physiological acclimation, which dynamically alters the composition of many species (Smith et al., 2011;Moreno and Martiny, 2018). However, many large-scale OBGC models assume constant elemental and pigmentary composition for phytoplankton (Bopp et al., 2013;Laufkotter et al., 2015), similar to the FS approach herein. Some models resolve diverse phytoplankton communities by explicitly accounting for as many as hundreds of differently adapted types (Follows and Dutkiewicz, 2011), and a few account for interspecific differences in composition (Dutkiewicz et al., 2019). Others represent these flexibilities in terms of physiological acclimation, typically using a Droop quota model (e.g., Moore et al., 2001), but seldom with a more consistent allocation  framework in large-scale OBGC models (Kerimoglu et al., 2017;Kwiatkowski et al., 2018;Pahlow et al., 2020). A few studies have even resolved variable composition for a diversity of phytoplankton types (e.g., Blackford et al., 2004;Ward et al., 2014;Butenschon et al., 2016). All of these approaches typically require additional state variables and calculations and therefore when applied to multiple phytoplankton functional types can be computationally expensive, as demonstrated by Kwiatkowski et al. (2014). To overcome these problems, at least for acclimative models, Smith et al. (2016) proposed the IA approach, which optimises Chl:C:nutrient ratios instantaneously for local conditions, so that variable composition can be tracked via a single state variable (biomass).
We have recently demonstrated that the IA approach behaves similarly to the fully dynamic version (DA) in a spatially explicit 1D setup (Kerimoglu et al., 2021). Here we assess four different models for phytoplankton growth against insitu observation at two oligotrophic stations with different physical conditions. Comparing the FS, FSPA, DQ, and IA models, which differ in their degrees of flexible composition and computational efficiency, has allowed us to disentangle the impact of different aspects of acclimative response on model performance and portability.

Cost Landscapes
The steep increase in cost with decreasing nutrient affinity for the FS and FSPA variants reveals their greater sensitivity, compared to the DQ and IA variants. This is because with the Monod formulation used in the FS and FSPA models, growth rate depends directly and immediately on the ambient nutrient concentration and the nutrient affinity,Â 0 , which determines nutrient uptake rate at low nutrient concentrations (Button, 1978;Healey, 1980;Smith et al., 2014). For the IA and DQ models, changes inÂ 0 do not so directly alter growth rate, because growth depends on the variable cell quota, Q. It has long been appreciated that internally stored nutrients, in excess of immediate requirements for growth, allow sustained growth under low ambient nutrient concentrations in temporally (Grover, 1991) or spatially (Kerimoglu et al., 2012) heterogeneous environments.
Previous studies have found that, compared to variable quota models, Monod-type models are more sensitive to fluctuations in nutrient supply, even in nutrient-replete regions (Cerco et al., 2004;Pauer et al., 2020). In low-nutrient regions, variations in nutrient supply (or affinity) would be expected to have even greater impact on growth rates, because of the typically saturating response for nutrient uptake (Button, 1978;Healey, 1980;Smith et al., 2014). The similar cost landscapes for DQ and IA models reveal that the IA model captures this storage effect despite its instantaneous adjustment of Q (hence C:N ratio). One reason for this is that the IA model's optimal resource allocation scheme compensates for lowÂ 0 values by increasing the fractional allocations of intracellular resources towards nutrient uptake broadly (f V ), and affinity specifically (f A ). The decoupling of growth rate from ambient nutrient concentrations may explain why the DQ and IA variants yield higher overall cost for the reference and simultaneous runs, compared to the FSPA model; i.e., regions where J DIN is high coincide with regions where J CHL and J PP are low. However, since these models have flatter overall cost landscapes, their costs differ less between the simultaneous and cross validation experiments (i.e., they are more portable), compared to the FS and FSPA models.

Model Performance With Individual Tuning
With parameters tuned individually for each station, compared to the FS variant, the photo-acclimative variants capture better common features in the oligotrophic ocean, such as the depletion of DIN and Chl in the summer, near-surface primary production, and deep SCM (Mongin et al., 2006;Mignot et al., 2014;Kavanaugh et al., 2018). Even with data-assimilation, models that have fixed stoichiometry are unable to reproduce these features (Schartau et al., 2001;Ayata et al., 2013;Ward et al., 2013). Because the FS variant cannot adjust its internal composition, i.e., decrease its Chl:C during the stratified summer, it often overestimates near-surface Chl and DIN, and underestimates SCM concentrations. Additionally, the lack of variable Q in the FS and FSPA variants limits their ability to capture continuously high PP through summer, as observed. This is because with the Monod formulation used in these models, growth and PP are directly proportional to nutrient uptake rate, and therefore decline instantaneously with ambient DIN. The four variants differ most in terms of Chl between 75 and 125 m, especially at ALOHA, where the optimal physiological acclimation in the IA variant sustains growth and hence higher biomass and Chl, with higher C:N ratio (lower Q) than the DQ during summertime nutrient depletion (Figures 6E,F). The FSPA model also captures the sustained high summer Chl concentration at the DCM, because it uses fixed Q from the biomass-averaged Q obtained from the IA model. Nevertheless, compared to the FS model, the addition of only photoacclimation in the FSPA variant indeed improves the overall model performance, yielding the lowest overall cost, especially in terms of DIN and Chl, when tuned separately for each station. This improvement in model performance may be because of the lack of variable cell quota which causes a lag in nutrient consumption, therefore consuming more nutrients and producing more Chl nearer the surface (thus capturing the observations better, as shown in Figures 6B,D). However, this hinders the model's ability to capture the seasonality of PP, which depends on Phy C as calculated from the cell quota.
A more thorough assessment of the realism of the model variants under consideration could in principle include a comparison of modelled POC:PON ratios against in situ observations (as in, e.g., Mongin et al., 2006), but that would require modelling a variety of other processes that alter the composition of detritus (e.g., digestion, excretion and detritivory by zooplankton, and differential degradation by bacteria), which would add considerable uncertainties and complicate the interpretation of model results. As a more direct alternative, it may soon be possible to compare model results against single-cell or species-level observations of cellular composition, which have only recently become available for certain types of phytoplankton in selected regions (e.g., Lopez et al., 2016;Baer et al., 2017).
In terms of PP, both the IA and DQ variants produce lower costs than the fixed quota variants. Although the FSPA variant performs substantially better than the FS in terms of Chl, it does not allow variable Q (hence, N:C), which limits its ability to capture the seasonality of PP. Additionally, inclusion of optimal resource allocation, as represented by the IA variant, indeed improves the simulation of PP, although the DQ variant captures better the maintenance of high PP during summer (Figures 4K,N,Z,C1), resulting in better PP cost at ALOHA than BATS. This may be because the DQ variant produces lower Q when µ is high due to the prescribed f V as well as its flatter µ vs. Q curve compared to the IA variant (see Supplementary Material), which smooths the transition between 'bloom' time and nonbloom times, as shown in the observed PP (Figures 4L,X, 6I,J).
Compared to the DQ model, the IA assumption constrains variations of Q more narrowly, which limits its ability to capture the storage effect. This is one reason why the IA model produces less Phy C (Kerimoglu et al., 2012(Kerimoglu et al., , 2021, and therefore PP, compared to the DQ model. In effect, the IA assumption restricts luxury uptake in order to allow faster growth immediately under nutrient-replete conditions (Kerimoglu et al., 2021), which can yield slower growth under subsequent nutrient scarcity compared to models that resolve more internal storage. This may also explain why the IA model overestimates Chl concentrations near the DCM.
Although almost all variants realistically simulate near-surface DIN ( Figure 6D) and elevated PP in spring, none captures the typical increase between May to August (Dave and Lozier, 2010) in the North Pacific Subtropical Gyre. This discrepancy may be due to the absence of nitrogen fixers in NFlexPD (Dore et al., 2008;Church et al., 2009;Böttjer et al., 2017) and nanophytoplankton, such as Prochlorococcus (Jiang et al., 2021), which despite being low in actual biomass, can contribute up to 40% in terms of PP (White et al., 2015;Kavanaugh et al., 2018). Indeed, it would be technically straightforward to extend the model with additional phytoplankton (or zooplankton) types, which may improve these results. However, as mentioned in the introduction, this would over-complicate the interpretation of the impacts of phytoplankton growth formulations, and their differing responses to parameter perturbations and physical settings (ALOHA vs. BATS). For this reason, we prefer to keep the model configuration as simple as possible in the current study. The lack of sporadic mixing, which lifts the thin layer of Chl towards the surface (Dore et al., 2008) and other processes, such as anticyclonic eddies or effects of El-Nino (Church et al., 2009;Kavanaugh et al., 2018) that are not well simulated by the 1D physical model may also be reasons for our underestimated PP values. The in situ PP measurements, based on 14 C incubations at both stations (Letelier et al., 1996;Steinberg et al., 2001), are also known to result in inaccurate PP values, including overestimates in the presence of slow-growing species Laws, 2013, 2014), and it is unclear whether they measure net or gross PP (Marra, 2009;Kavanaugh et al., 2018).

Model Portability
The physiological acclimation, which aims to optimise growth rate in the IA and DQ variants, is similar to dynamic optimisation of parameter values (e.g., Mattern et al., 2012), as is typically applied to capture complex adaptive behaviours in response to environmental changes (Arhonditsis and Brett, 2004), e.g., physiological plasticity or succession of plankton groups (Follows and Dutkiewicz, 2011). It is therefore expected that the IA variant, which continually re-allocates resources to optimise growth, would capture better DIN and Chl dynamics at both stations in the cross-validation and simultaneous experiments. Although earlier studies (Friedrichs et al., 2007;Kriest et al., 2012;Ward et al., 2013) found that increasing model complexity does not necessarily improve misfits or predictive capability, here we found that the inclusion of photoacclimation improves the model portability. That is, the more "complex" acclimative models overall perform better compared to the simpler FS variant. However, a previous study from Ward et al. (2013) showed that with objective parameter tuning, a biogeochemical model is sufficient to capture the observation without photoacclimation at the temperate North Atlantic. A future investigation including data assimilation would provide a more thorough and conclusive assessment of how the different aspects of photoacclimation impact model performance and portability at both temperate and oligotrophic stations.
As explained in subsection 4.1, the DQ and IA variants differ least between the simultaneous and cross validation experiments, and therefore are the two most portable. However, because for the DQ model the same value ofÂ 0 is applied for BATS and ALOHA, its cost varies less with α I than withÂ 0 . Therefore, in this particular experiment the DQ model is more portable than the IA variant. However, the IA variant generally produces better statistical metrics for all the observed variables than the DQ model. Thus, when applied to a 3D regional OBGC model, for oligotrophic regions, the IA variant can be expected to more realistically capture phytoplankton growth, nutrient uptake, and chlorophyll concentrations with fewer state variables.

Future Outlook
This study shows that adding photoacclimation improves model performance at two oligotrophic stations, and adding variable phytoplankton C:N makes models less sensitive to parameter perturbations, especially those related to affinity, and therefore improves portability. The sensitivity analysis could in future be expanded by perturbing other parameters, such as mortality and potential maximum growth, which are directly linked to phytoplankton biomass. The portability experiment could be extended by applying data assimilation to fit each model variant to the observations, as done previously in a 0-D setup (Smith et al., 2016), and by performing cross validation experiments between more contrasting (e.g., subpolar and subtropical) regions. Although both the IA and DQ variants capture Chl, PP, and DIN concentrations fairly well, the idealised models presented herein are too simplistic to fully capture oligotrophic ecosystem dynamics. Multiple plankton types inhabit oligotrophic regions, including zooplankton (Dave et al., 2015) and N-fixers (Dore et al., 2008;Karl et al., 2012) as well as other nutrients, such as phosphate (Steinberg et al., 2001;Karl et al., 2012); but our models only represent one type of phytoplankton and nutrient. Adding predation or another trophic level, depending on the grazing formulation, may change the depth of the DCM, and may also change the timing and magnitude of phytoplankton blooms . However, how this may change the portability remains an open question. Three dimensional modelling studies will be required to evaluate the impacts of physical processes such as mesoscale eddies and inter-annual variations driven by El-Nino/La-Nina events, which our 1-D setup cannot capture. For more comprehensive OBGC model applications, it will also be essential to trace C, O 2 , and alkalinity (Kwiatkowski et al., 2014). Behrenfeld et al. (2016) showed that across most of the ocean variations in Chl, the most widely observed metric of phytoplankton, result more from physiological acclimation than from variations in their biomass. Models that account for photoacclimation can help to disentangle the mechanisms underlying observed Chl variations. In a broader sense, although the effects of physiological acclimation and evolutionary adaptation can in many cases be modelled with the same mathematical approaches (Smith et al., 2011), they can be expected to have distinct effects on large scale chlorophyll distributions and biogeochemistry (Moreno and Martiny, 2018). Disentangling their relative contributions remains a formidable challenge, which will require resolving explicitly both intra-and inter-specific variability in models. Our results suggest that for such purposes, the IA approach may be useful to reduce the massive computational cost of resolving acclimation processes for many phytoplankton types.

CONCLUSIONS
We test an NPD-type model with four different variants differing with respect to phytoplankton flexible composition and acclimative response: A Monod-type variant with fixed stoichiometry (FS), an improved monod-type variant with photoacclimation (FSPA), a variant with both variable composition using a classical Droop-quota approach and photoacclimation (DQ), and finally the Instantaneous Acclimation (IA) variant. We assess whether adding variable composition and acclimation can enhance generality and portability at two oligotrophic stations BATS and ALOHA. Additionally, we conduct a sensitivity analysis to explore how different phytoplankton growth formulations respond to parameter perturbations, especially those related to nutrient and light affinity. When individually tuned, we find that FSPA model captures Chl and DIN observations at both stations better than the other three variants, but the DQ and IA variants with variable C:N and Chl:C ratios capture the observed PP better and are less sensitive to parameter perturbations. Therefore, flexible composition and photoacclimation can be expected to enhance portability and hence potential applicability at large scales. The addition of optimal resource allocation slightly enhances model performance as quantified by statistical metrics. However, the DQ model, which calculates the dynamics of phytoplankton C and N explicitly, is more portable than the IA approach. It should be noted that these experiments are done without any formal parameter optimisation and limited to oligotrophic sites. Further studies using objective optimisation (data assimilation) and including more contrasting regions would provide a more comprehensive and thorough assessment.

AUTHOR CONTRIBUTIONS
PA performed the model runs, analysed the outputs, and prepared the figures. All authors designed the study and contributed equally to the writing of the manuscript.