Optimality-based non-Redfield plankton–ecosystem model (OPEM v1.1) in UVic-ESCM 2.9 – Part 1: Implementation and model behaviour

Uncertainties in projections of marine biogeochemistry from Earth system models (ESMs) are associated to a large degree with the imperfect representation of the marine plankton ecosystem, in particular the physiology of primary and secondary producers. Here, we describe the implementation of an optimality-based plankton– ecosystem model (OPEM) version 1.1 with variable carbon : nitrogen : phosphorus (C : N : P) stoichiometry in the University of Victoria ESM (UVic; Eby et al., 2009; Weaver et al., 2001) and the behaviour of two calibrated reference configurations, which differ in the assumed temperature dependence of diazotrophs. Predicted tracer distributions of oxygen and dissolved inorganic nutrients are similar to those of an earlier fixedstoichiometry formulation in UVic (Nickelsen et al., 2015). Compared to the classic fixed-stoichiometry UVic model, OPEM is closer to recent satellite-based estimates of net community production (NCP), despite overestimating net primary production (NPP), can better reproduce deep-ocean gradients in the NO−3 : PO 3− 4 ratio and partially explains observed patterns of particulate C : N : P in the surface ocean. Allowing diazotrophs to grow (but not necessarily fix N2) at similar temperatures as other phytoplankton results in a better representation of surface Chl and NPP in the Arctic and Antarctic oceans. Deficiencies of our calibrated OPEM configurations may serve as a magnifying glass for shortcomings in global biogeochemical models and hence guide future model development. The overestimation of NPP at low latitudes indicates the need for improved representations of temperature effects on biotic processes, as well as phytoplankton community composition, which may be represented by locally varying parameters based on suitable trade-offs. The similarity in the overestimation of NPP and surface autotrophic particulate organic carbon (POC) could indicate deficiencies in the representation of top-down control or nutrient supply to the surface ocean. Discrepancies between observed and predicted vertical gradients in particulate C : N : P ratios suggest the need to include preferential P remineralisation, which could also benefit the representation of N2 fixation. While OPEM yields a much improved distribution of surface N* (NO−3 −16 ·PO 3− 4 +2.9 mmolm −3), it still fails to reproduce observed N* in the Arctic, possibly related to a misrepresentation of the phytoplankton community there and the lack of benthic denitrification in the model. Coexisting ordinary and diazotrophic phytoplankton can exert strong control on N* in our simulations, which questions the interpretation of N* as reflecting the balance of N2 fixation and denitrification.

cates the need for improved representations of temperature effects on biotic processes, as well as phytoplankton community composition, which may be represented by locally varying parameters based on suitable trade-offs. The similarity in the overestimation of NPP and surface autotrophic particulate organic carbon (POC) could indicate deficiencies in the representation of top-down control or nutrient supply to the surface ocean. Discrepancies between observed and predicted vertical gradients in particulate C : N : P ratios suggest the need to include preferential P remineralisation, which could also benefit the representation of N 2 fixation. While OPEM yields a much improved distribution of surface N* (NO − 3 − 16 · PO 3− 4 + 2.9 mmol m −3 ), it still fails to reproduce observed N* in the Arctic, possibly related to a misrepresentation of the phytoplankton community there and the lack of benthic denitrification in the model. Coexisting ordinary and diazotrophic phytoplankton can exert strong control on N* in our simulations, which questions the interpretation of N* as reflecting the balance of N 2 fixation and denitrification. a rather wide range, e.g. in the predicted temperature increase until the end of the current century (IPCC, 2013). Some predictions do not even agree in the sign of the projected changes, e.g. of marine net primary production, particularly in low latitudes, varying between −25 % and +40 % across current models (Laufkötter et al., 2015; see also Taucher and Oschlies, 2011). But even where many ESMs agree, their predictions are sometimes counter to observations, e.g. in the case of oceanic O 2 patterns and trends . These problems are likely rooted in uncertainties in parameter estimates (Löptien and Dietze, 2017) but also inherent model deficiencies, such as limited spatiotemporal resolution or inaccurate representation of physical and biotic processes (Keller et al., 2012;Getzlaff and Dietze, 2013).
In our view, a major limitation of the biogeochemical modules of current ESMs is that the formulations used to describe the plankton compartments are at odds with organism behaviour as observed in the laboratory. While the variability of the chlorophyll : carbon (Chl : C) ratio is considered in recent ESMs (e.g. Park et al., 2019), the carbon : nitrogen : phosphorus (C : N : P) stoichiometry of phytoplankton is still often represented by static (Redfield) ratios, entirely ignoring its highly variable nature (Klausmeier et al., 2008), which can affect model sensitivity to climate change (Kwiatkowski et al., 2018). The only model with variable C : N : P in phytoplankton in CMIP5 (Bopp et al., 2013) and CMIP6 (Arora et al., 2020) is PELAGOS (Vichi et al., 2007), which has no diazotrophs. Other models consider only variable N : P (TOPAZ2; Dunne et al., 2012) or C : P (MARBL (CESM2); Danabasoglu et al., 2020). The problem extends also to the representation of fundamental biotic processes, such as nutrient uptake or zooplankton foraging. For example, Smith et al. (2009) showed that the halfsaturation concentration of nitrate use varies systematically with nitrate concentration and suggested that optimal uptake kinetics (Pahlow, 2005) may be more appropriate than the commonly used Michaelis-Menten kinetics for simulating phytoplankton nutrient uptake. Zooplankton foraging behaviour can be characterised by a significant feeding threshold followed by a steep increase in ingestion (e.g. Kiørboe et al., 1985;Strom, 1991;Gismervik, 2005), which has also been demonstrated for a natural plankton community in the Sargasso Sea (Lessard and Murrell, 1998). This kind of feeding behaviour may be important for capturing the distribution of primary production in large ocean areas (Strom et al., 2000), but it is not represented by the Holling type II and III models (Holling and Buckingham, 1976) used in current biogeochemical models.
We have recently developed optimality-based formulations for phytoplankton and zooplankton (Pahlow and Prowe, 2010;Pahlow et al., 2013), which can describe observed plasticity of plankton organisms, yet are sufficiently simple for implementation in global biogeochemical models. Plasticity here refers to the variability of elemental composition and allocation of resources among competing require-ments for light harvesting and nutrient acquisition in phytoplankton and for foraging and digestion in zooplankton, implying variable Chl : C : N : P stoichiometry, half-saturation concentrations for nutrient uptake, and ability to fix nitrogen in phytoplankton, and zooplankton feeding thresholds and variable assimilation efficiency. The optimality concept is based on the "assumption that natural selection should tend to produce organisms optimally adapted to their environments" (Smith et al., 2011), which is particularly applicable to marine plankton, where intense mixing and the absence of physical boundaries ensure strong competition, and short generation times allow for rapid evolution. These formulations have shown their ability to describe ecosystem behaviour in 0-D and 1-D modelling studies (e.g. Fernández-Castro et al., 2016;Su et al., 2018), and to predict patterns of phytoplankton nutrient and light co-limitation based on satellite and in situ observations (Arteaga et al., 2014). In this contribution, we describe the implementation of our new optimality-based plankton-ecosystem model (OPEM) into a global 3-D ocean model component of an ESM of intermediate complexity. All of the new assumptions in OPEM are based on published experimental observations used to validate the optimality-based formulations. We view the implementation of OPEM as one step towards the ultimate goal of reconciling plankton-organism behaviour as observed in the laboratory with global marine biogeochemistry. Therefore, the variable stoichiometry of primary producers should be considered but one, albeit central, aspect of the mechanistic foundation of OPEM. The ESM employed is the University of Victoria Earth System Climate Model (UVic-ESCM) (UVic in the following; Eby et al., 2009;Weaver et al., 2001). Due to its coarse spatiotemporal resolution, UVic is a practical choice when working on long timescales (e.g. Niemeyer et al., 2017) and/or when many simulations are needed. Computational efficiency is also one of the main impediments to introducing more mechanistic formulations of biotic processes (Chen and Smith, 2018), as, e.g. the representation of variable C : N : P stoichiometry requires additional tracers, which must be mixed and advected as well. UVic has been used extensively with typical state-of-the-art fixed-stoichiometry NPZD (nutrientsphytoplankton-zooplankton-detritus)-type marine ecosystem and biogeochemistry models (e.g. Keller et al., 2012;Niemeyer et al., 2017;Oschlies et al., 2017). Here, we compare the behaviour of the OPEM with that of a previous UVic configuration, described in Nickelsen et al. (2015), modified with several improvements and bug fixes as described below. An empirically founded temperature dependence of diazotrophy is introduced in a second configuration, OPEM-H, in order to distinguish between effects of the optimalitybased physiological regulation and the temperature formulation. Since the calibration of OPEM and OPEM-H embedded in UVic presents a major challenge, it is dealt with in Part 2 (Chien et al., 2020).

Optimality-based plankton in the UVic model
UVic model version 2.9 (Weaver et al., 2001;Eby et al., 2013) in the configuration of Nickelsen et al. (2015), with the isopycnal diffusivity modifications by Getzlaff and Dietze (2013), vertically increasing sinking velocity of detritus (Kriest, 2017), and several bug fixes (some of which were already introduced by Kvale et al., 2017, see Appendix A for the new bug fixes applied here), is referred to as the original UVic in the following. We base our new configurations on this original UVic, except that we use constant halfsaturation iron concentrations and omit the upper temperature limit in the zooplankton temperature dependence. For OPEM, we replace the formulations for phytoplankton, diazotrophs, and zooplankton in the original UVic model with an optimality-based model (Pahlow et al., 2013) for phytoplankton and diazotrophs, and the optimal current-feeding model (Pahlow and Prowe, 2010) for zooplankton ( Fig. 1). Negative concentrations have always occurred in the UVic model, but they have usually been confined to small negative numbers in a few places. However, negative concentrations turned out to be a major problem for OPEM, which had to be dealt with in order to stabilise our optimality-based variablestoichiometry implementation (see Appendix B).

Phytoplankton and diazotrophs
Ordinary and diazotrophic phytoplankton are described by the optimal-growth model (OGM) of Pahlow et al. (2013), modified to account for the coarse spatiotemporal resolution of UVic and augmented with temperature and iron effects (see equations provided below). Due to the relatively long time step, the model does not resolve the dynamics of photoacclimation, and we therefore describe the Chl : C ratio of the chloroplast by its balanced-growth optimum (see Eq. C5 in Appendix C1.1). Hence, we do not need state variables for Chl. Simulating variable Chl : C : N : P stoichiometry in phytoplankton then requires three state variables, representing particulate organic carbon, nitrogen, and phosphorous (POC, PON, and POP) for each phytoplankton group and for detritus.
The OGM is a cell-quota model comprising several levels of physiological regulation. At the whole-cell level, resources are optimally allocated between nutrient acquisition and CO 2 fixation, Chl synthesis is optimised within the chloroplast, and optimal uptake kinetics (Pahlow, 2005;Smith et al., 2009) drives nutrient uptake and assimilation inside the protoplast. For all trade-offs, we define optimal as yielding maximum balanced growth of the cell. For facultative diazotrophs, N 2 fixation is switched on whenever this enhances growth. The biological model parameters of the OGM are different from the original UVic configuration. In spite of its ability to describe two additional tracers (phytoplankton C and P) and the Chl : C ratio, the OGM has only eight parameters (maximum rate V 0 , nutrient affinity A 0 , costs of N as-similation ζ N and Chl synthesis ζ Chl and maintenance R Chl M , subsistence quotas Q N 0 and Q P 0 , and the light-absorption coefficient α), i.e. the same as the phytoplankton parameters of the original UVic configuration (Nickelsen et al., 2015). In addition, two of these (V 0 and ζ Chl ) can be considered constant (Pahlow et al., 2013), leaving six parameters to be calibrated.
None of the measures against negative concentrations (Appendix B) are effective if the minimum required concentration of a tracer is greater than zero, which is the case for our phytoplankton PON and POP tracers, whose minimum (subsistence) concentrations are given by the product of POC and the N and P subsistence quotas Q N 0 and Q P 0 , respectively, which can be thought of as the subsistence PON and POP of phytoplankton. In order to circumvent this problem and also be able to benefit from the flux-corrected transport (FCT) technique (see Appendix B), we define δ tracers as the differences between actual and subsistence phytoplankton PON and POP concentrations. As the lower limit of the δ tracers is 0, they can be transported with the positive transport schemes, and subsistence PON and POP are implicitly advected and mixed in proportion to phytoplankton POC and added back onto the δ tracers where required: where C p , N p , and P p are POC, PON, and POP, respectively, of phytoplankton group p (phytoplankton or diazotrophs). The local rates of change of the phytoplankton tracers are then defined by sources-minus-sinks terms (S): where µ p is net relative (C-specific) growth rate (C fixation minus the sum of respiration and release of dissolved organic carbon by phytoplankton, immediately respired to DIC here), λ p leakage, M p mortality, G n p grazing by zooplankton, V N p and V P p DIN and DIP uptake, and Q N p and Q P p biomassnormalised N and P cell quotas (N : C and P : C ratios). The last term in Eq. (3) accounts for the subsistence amounts of N and P implicitly contained in C p and subtracted from δn p via Eq. (1). Leakage is the fast-recycling term parameterising the microbial loop (Keller et al., 2012). Definitions for all terms in Eqs.
We set up configurations with two representations of temperature dependence for diazotrophs: (1) configuration OPEM with the same temperature dependence as in the original UVic, and (2) configuration OPEM-H with the Eppley (1972) temperature dependence applied to both phytoplankton (subscript phy) and diazotroph (subscript dia) growth and temperature-dependent potential rates of C, N, P acquisition Zooplankton and detritus prey-specific rate of C, N, P ingestion g max , g zoo d −1 reference, actual relative rate of total ingestion M zoo m 3 (mol C) −1 d −1 zooplankton mortality µ zoo d −1 net relative growth rate ν det d −1 detritus reference decay rate C , N , P mol m −3 effective prey C, N, P concentration φ p m 3 (mol C) −1 prey-capture coefficients, p ∈ {phy, dia, det, zoo} Q N zoo , Q P zoo mol (mol C) −1 zooplankton N : C, P : C ratio R C zoo , R N zoo , R P zoo mol m −3 d −1 respiration, dissolved N, P loss r Q -stoichiometric reduction factor S g -degree of ingestion saturation X C zoo , X N zoo , X P zoo mol m −3 d −1 particulate C, N, P loss (egestion) Figure 1. Optimality-based plankton-ecosystem model (OPEM, panel a). Ordinary phytoplankton, diazotrophs, and zooplankton are represented by optimality-based physiological regulatory formulations. Ordinary phytoplankton and diazotrophs are driven by optimal allocation of cellular resources (b), balancing the benefits of nutrient assimilation and light harvesting against allocation and energetic costs (respiration, R) of these processes. The optimal allocation trades off, e.g. cellular N as defined by Q N , between the requirements for photosynthesis (green) and nutrient acquisition (blue), with an additional compartment for N 2 fixation in diazotrophs (not shown). The phosphorus quota (Q P ) controls N assimilation (see Appendix C1.2), but only Q N affects the growth rate directly (see Appendix C1.1). Zooplankton foraging (c) is optimised by balancing costs and benefits of allocating total activity (A t ) between foraging activity (A f ) and assimilation activity (A t − A f ). Both foraging and assimilation incur energy costs (c f and c a , respectively) fuelled by respiration (R). Increasing ingestion (g) reduces assimilation efficiency (E ≤ E max ), causing more particulate egestion (X). nutrient uptake, and the temperature function from Houlton et al. (2008) for N 2 fixation ( Fig. 2; see Appendix C1.3). The maximum temperature-dependent rates for diazotrophs are multiplied by 0.4 in the original UVic but not in OPEM, so that they remain below those of ordinary phytoplankton for the whole temperature range in Fig. 2. All other temperature dependencies are unchanged from the original UVic; i.e. they follow the Eppley (1972) curve (dashed red line in Fig. 2).

Zooplankton
Zooplankton foraging is described by the model of optimal current feeding (OCF; Pahlow and Prowe, 2010). The OCF is based on the idea that the animal has a certain inherent maximum total activity (A t ), which can be allocated between foraging activity (A f ) and activity for the assimilation of food (A t − A f ), so that the net relative growth rate is maximised, considering the costs of foraging and assimilation (represented by the coefficients c f and c a , respectively). While A t is a rather abstract quantity, it can be expressed as a function of the maximal ingestion rate, which is routinely determined in feeding experiments, and temperature (see Eq. C19 in Appendix C2). The OCF can represent different foraging strategies via its prey-capture coefficient (φ) and c f . Very low φ and c f ≈ 0 represent ambush feeding, whereas c f ≈ c a is representative of current feeding for intermediate φ and cruise feeding for high φ. The parameter values in OPEM and OPEM-H (Table 2) are between values determined for cruise and current feeders by Pahlow and Prowe (2010). The OCF has two more parameters than the original UVic, but since two of them can be considered constant (β = 0.2 and E max = 1; Pahlow and Prowe, 2010), the number of parameters which have to be calibrated is the same as in the original UVic. Besides its mechanistic foundation, the main advantages over the Holling type II formulation in the original UVic model are the predicted feeding threshold and variable assimilation efficiency. Assimilation efficiency is constant and a feeding threshold does not exist in the original UVic model. Temperature dependence is accounted for by multiplying the maximum ingestion rate and maintenance respiration with the temperature function as described in Keller et al. (2012) but here without the cap at 20 • C. The cap on the increase of maximum ingestion rate with grazing in the original version was deemed necessary in order to avoid inordinately high grazing in the tropics (Keller et al., 2012). It is noteworthy that this does not appear to be a problem in OPEM even though maximum ingestion rates g max are about 4-fold higher than in the original UVic version (Table 2). We attribute this to the feeding threshold in the OCF, which reduces grazing in oligotrophic regions. Since zooplankton stoichiometry is fixed (constant Q N zoo and Q P zoo ) but that of the food is variable, any excess C, N, or P must be released, assumed here in mostly dissolved form (as inorganic nutrients). For example, all the excess ingested C is respired (see Eq. C16 in Appendix C2), as also suggested by Talmy et al. (2016). To this end, we define a stoichiometric reduction factor r Q that reduces net uptake and growth of zooplankton to the uptake of the most limiting nutrient of the ingested food: where n is the effective prey concentration for nutrient element n and φ p is the prey-specific capture coefficient. The relations among the φ p effectively determine the (relative) food preferences. The sources-minus-sinks term for zooplankton biomass S(N zoo ) is expressed here in terms of nitrogen, which can easily be converted to P and C via the zooplankton's fixed stoichiometry. S(N zoo ) is the difference between net growth (µ zoo ), which is corrected for r Q (Appendix C2), and losses due to intra-guild predation (G N zoo ) and background mortality (M zoo ): Equations for µ zoo and G N zoo are given in Appendix C2. The background mortality is a quadratic closure term intended to represent losses due to viruses, predation by higher trophic levels, etc.

Detritus and dissolved pools
Mortality terms and egestion of faecal particles by zooplankton produce detritus, which is itself subject to grazing and temperature-dependent remineralisation. We consider separate C, N, and P tracers for detritus: where ν det is the detritus remineralisation rate at 0 • C. Hence, the export and remineralisation fluxes are also traced individually for C, N, and P. This applies also to alkalinity, where we assume a sulfur-to-carbon ratio of 0.023 mol S mol C −1 for organic C (Matrai and Keller, 1994). For O 2 consumption during remineralisation, we consider contributions from C and N separately. We assume −O 2 : N = 2 (the N contribution to O 2 consumption) during nitrification and calculate the respiratory quotient for C based on an O 2 : C ratio of 170 : 117 = 1.45 mol O 2 mol C −1 (Anderson and Sarmiento, 1994), corrected for the contribution of nitrification, and an average C : N = 6.625 mol C mol N −1 . Thus, we obtain the respiratory quotient for C (the C contribution) as the difference between the average O 2 : C ratio and the N contribution to O 2 consumption, i.e. 1.45 mol O 2 mol C −1 − 2 mol O 2 mol N −1 /6.625 mol C mol N −1 = 1.15 mol O 2 mol C −1 . Equation (6) does not include gains and losses from sinking detritus particles. Detritus sinking speed v sink increases with depth, reflecting the remineralisation of more slowly sinking smaller particles, leading to a dominance of fast-sinking (typically larger) particles at greater depths: where v 0 = 6 m d −1 is the sinking velocity at the surface, z is depth and a v = 0.06 d −1 the rate of increase in v sink with depth (Kriest, 2017). Dissolved inorganic C and nutrients are utilised by phytoplankton and released by phytoplankton leakage, zooplankton respiration and excretion, and detritus remineralisation, as well as via rejection of surplus elements via grazing of organic matter with elemental stoichiometries differing from that of zooplankton.

Model reference simulations
We first did a preliminary sensitivity analysis to identify sensitive model parameters. Then we set up an ensemble of 400 parameter sets, using a Latin hypercube method, and ran both of our model configurations into steady state for all parameter sets. We select two reference simulations (trade-off solutions in Part 2; Chien et al., 2020), one each from the OPEM Table 2. Parameter settings for the original and our reference OPEM and OPEM-H configurations. Parameters in bold vary within the ensembles of simulations (Chien et al., 2020). Symbol descriptions are given in Table 1.
and OPEM-H ensembles, according to two objectives: (1) we minimise a cost function under the condition that (2) we obtain realistic levels of global water-column denitrification, i.e. at least 60 Tg N yr −1 (DeVries et al., 2012). Thus, no weighting had to be applied to our objectives. The cost function quantifies the model-data misfit by a measure of the discrepancies between observed and simulated O 2 , NO − 3 , PO 3− 4 , and Chl, considering also correlations and covariances (see Part 2; Chien et al., 2020).
In the following, we describe and discuss the behaviour of the two reference simulations, which turned out to have same parameter set (Table 2). While this may be a coincidence, it has the advantage that all differences between OPEM and OPEM-H can be ascribed unequivocally to the difference in the temperature dependence of the diazotrophs. We specifically consider the models' ability to reproduce features not included in the cost function, namely the surplus nitrate with respect to the Redfield N equivalent of phosphate, termed N * = NO − 3 − 16 · PO 3− 4 + 2.9 mmol m −3 (Gruber and Sarmiento, 1997), where the constant factor 0.87 was dropped as recommended by Mills et al. (2015), and global N 2 -fixation rates and distributions within current observa-tional ranges. All our UVic model results are shown as annual averages at the end of the spin-up (i.e. after at least 10 000 years), when a seasonally cycling steady state has been reached.
We compare the predictions of our reference simulations with data from these sources: NO − 3 , PO 3− 4 , and O 2 data are from the World Ocean Atlas 2013 annual objectively analysed mean fields (WOA 2013; Garcia et al., 2013a, b). Dissolved inorganic carbon (DIC) data are from GLobal Ocean Data Analysis Project version 2 (GLODAPv2) (Key et al., 2015;Lauvset et al., 2016). Estimates of Chl (MODIS Aqua, level 3; https://oceancolor.gsfc.nasa.gov/l3/, last access: 11 November 2019; Hu et al., 2012), particulate organic carbon, and net primary and community production (POC, NPP, and NCP; Westberry et al., 2008;Li and Cassar, 2016) are based on satellite data. In situ N 2 -fixation data are from MARine Ecosystem biomass DATa (MAREDAT) (Luo et al., 2012).    (Fig. 3). The original UVic better reproduces the NO − 3 profile above 1000 m than OPEM and OPEM-H but overestimates NO − 3 below 2000 m. The DIC and PO 3− 4 pro-files from our reference simulations are very similar to those of the original UVic model (Fig. 3). Surface nitrate concentrations are generally slightly higher and more evenly distributed in OPEM and OPEM-H than in the original UVic model (Fig. 4). For most of the Atlantic, OPEM and OPEM-H are closer to the WOA 2013 data. Surface NO − 3 in the Indian Ocean is underestimated by the original UVic and overestimated by OPEM and OPEM-H. Surface patterns of N* are much closer to observations in both OPEM and OPEM-H than in the original UVic configuration (Fig. 5). However, while N* in the northern North Pacific and Arctic oceans is lower in OPEM and OPEM-H than in the original UVic, all UVic configurations still fail to reproduce the very low N* in large parts of the North Pacific and Arctic oceans (Fig. 5). While N 2 fixation is not limited to temperatures higher than 15 • C in OPEM-H, only very little N 2 fixation occurs in the high northern and southern latitudes and thus cannot explain the higher surface N* values in OPEM-H there (see Sect. 3.3 below). In our model simulations, low N* in the eastern tropical Pacific and South Atlantic result from denitrification in underlying oxygen-minimum zones (OMZs) (Landolfi et al., 2013). The original UVic configuration also displays very low N* in the Andaman Sea, whereas results of OPEM and OPEM-H are somewhat closer to the WOA 2013 data in the northern Indian Ocean (Fig. 5).
Interestingly, these differences cannot be seen in the O 2 distribution at 300 m, the depth of the OMZs, which is very similar in the Indian Ocean and eastern tropical Pacific among all our UVic simulations (Fig. 6), indicating that the carbon export and subsequent remineralisation are very similar as well. The main differences in O 2 distribution are that O 2 is slightly higher in the Arctic Ocean and slightly lower in the equatorial Pacific and northern North Pacific in both OPEM and OPEM-H compared to the original UVic (Fig. 6).
The OPEM simulations allow for a variable C : N ratio in detritus leaving the surface layers and reveal C : N ratios higher than the canonical value of 6.625 mol C (mol N) −1 , which is also the stoichiometry of zooplankton, almost everywhere between 40 • S and 40 • N in OPEM and OPEM-H (Fig. 7). Even though detritus C : N is lower in the Bay of Bengal than in the remainder of the Indian Ocean in both OPEM simulations, this feature cannot explain the lower denitrification compared to the original UVic in this area, since the C : N ratio, which determines the O 2 demand for the remineralisation of sinking detritus, remains above the original UVic value of 6.625 mol C (mol N) −1 . Rather, the lack of denitrification in the Indian Ocean in OPEM and OPEM-H ( Fig. 13 right) appears to result simply from the reduced C export in this area compared to the original UVic (Fig. 12).
Another interesting feature of the OPEM and OPEM-H simulations is their ability to reproduce, at least qualitatively, the global-scale gradient of DIN : DIP ratios in the deep ocean (Fig. 8). The WOA 2013 data indicate relatively high DIN : DIP in the deep North Atlantic, decreasing towards the Southern, Indian, and Pacific oceans. This gradient is very weak (and reversed) in the original UVic model (Fig. 8). Also, not all simulations in our OPEM and OPEM-H ensembles can reproduce this gradient, whereas other models without variable stoichiometry can (e.g. Kriest and Oschlies, 2015). Thus, reproducing the deep DIN : DIP distribution appears to mostly require a suitable model calibration. Note that deep-water N : P ratios are systematically higher in OPEM-H compared to OPEM, because of the elevated N* values in OPEM-H in high-latitude surface waters that feed the deepocean interior (Fig. 5). We interpret the surface N* distribution outside the deep-water formation regions as a consequence rather than a cause of the deep-ocean nutrient distributions, however.

Chlorophyll, primary production, and autotrophic biomass
Chlorophyll concentrations are generally more evenly distributed in OPEM and OPEM-H, which agrees better with the MODIS Aqua (level 3) satellite estimates (Hu et al., 2012) than the original UVic model, which also overestimates chlorophyll in the tropics and the Indian Ocean more pronouncedly. Only the OPEM-H simulation predicts reasonably high chlorophyll in the Arctic Ocean compared to the satellite estimates ( Fig. 9). OPEM and OPEM-H apparently overestimate surface Chl in the oligotrophic subtropical gyres compared to the satellite estimate, which may be partly explained by both the inability of the satellite to detect deep chlorophyll maxima (DCM) and the coarse vertical resolution of the UVic grid. Unlike the original UVic, OPEM and OPEM-H have variable Chl : C ratios leading to pronounced DCM in the second layer (not shown). The surface layer in the UVic grid is 50 m thick, i.e. much thicker than the surface mixed layer in the typically strongly stratified oligotrophic subtropical gyres. Thus, the model underestimates light and overestimates nutrient supply to the surface in these regions, both of which tend to raise the Chl : C ratio (Pahlow et al., 2013), so that some of the high predicted surface Chl concentrations can be understood as the manifestation of an unresolved DCM within UVic's surface layer. As discussed below, however, part of the high Chl prediction also reflects an overestimation of autotrophic biomass (POC). NPP and NCP are defined here as where µ is the net relative growth rate, λ the leakage rate representing fast remineralisation in UVic, f det (T ) · ν det · C det detritus remineralisation, and R C zoo zooplankton respiration, defined in Eq. (C16) in Appendix C2. NCP represents the net production of organic carbon after accounting for the metabolic needs of the autotrophic and heterotrophic components of the ecosystem (Ducklow and Doney, 2013). NPP in OPEM is the same as in OPEM-H (88.0 Pg C yr −1 ) and is    much higher than the estimate from Westberry et al. (2008) of 52 Pg C yr −1 , which in turn exceeds that in the original UVic model (44.3 Pg C yr −1 ). The NPP for the original UVic is lower than previously published (55 Pg C yr −1 ; Nickelsen et al., 2015) because we include λ phy in Eq. (8). The global averages predicted by the OPEM and OPEM-H simulations are slightly higher than the range of predictions from ocean colour-and model-based estimates reported by Carr et al. (2006). NPP is much more evenly distributed in OPEM and OPEM-H than in the original UVic model, but the carbon-based productivity model (CbPM) (Westberry et al., 2008) predicts an even more uniform distribution (Fig. 10). The original configuration clearly underestimates NPP in the oligotrophic gyres, whereas OPEM and OPEM-H overestimate NPP in the tropical ocean. The high predicted NPP in OPEM and OPEM-H (Fig. 10) is apparently linked to an overestimate of (1) autotrophic biomass throughout most of the World Ocean (Fig. 11) and (2) (Fig. 4). In addition, the 50 m thick surface layer in the UVic grid implies that the integrated biomass may be overestimated even more strongly than the surface POC concentration, particularly under stratified conditions. A possible explanation for the discrepancy between the OPEM and CbPM predictions may be that we do not include light affinity (α) among the list of parameters to be calibrated, because this parameter showed relatively little effect during our preliminary sensitivity analysis used to select sensitive model parameters. However, Arteaga et al. (2016) found that simple adaptive equations for α and A 0 , meant to represent adaptation to nutrient or light limitation, greatly improved predicted Chl : C compared to constant α and A 0 as applied in the present study. The use of constant parameters means that the OPEM and OPEM-H represent physiological flexibility as observed within species but do not consider variations in plankton community composition.
Comparing the patterns in NPP and surface autotrophic POC (Figs. 10 and 11) suggests a spatial correlation between deviations in these two quantities (Fig. 11, lower left). Thus, some of the NPP overestimate could result from an overestimate in POC: the predicted NPP in both OPEM and OPEM-H is 1.7 times the CbPM estimate in Fig. 10 and the average surface autotrophic POC in OPEM and OPEM-H is 1.4 and 1.7 times that of the satellite-based CbPM estimate in Fig. 11. We interpret this as indicating that the growth rates of the primary producers may be relatively well represented by their optimality-based formulation, but the model behaviour might benefit from improvements in the representation of top-down control. While the growth of the primary producers is defined by the optimality-based formulation of phytoplankton introduced here, mortality is only partly determined by the optimal current-feeding model employed to describe zooplankton behaviour. A large part of phytoplankton mortality is still due to the mortality terms of the original UVic. The importance of top-down control becomes apparent from the result that autotrophic POC is much greater than the zooplankton feeding threshold throughout most of the World ocean in OPEM and OPEM-H (contours in the right panels of Fig. 11). Thus, the feeding threshold itself appears to be reasonable compared to the satellite-derived autotrophic POC, but our zooplankton somehow fails to exert sufficient top-down control when food availability is high.
Net community production (NCP) is spatially more evenly distributed in OPEM and OPEM-H than in the original UVic model. Both the more evenly distribution and the subsequently higher global total NCP are much closer to the satellite-based estimate of Li and Cassar (2016) than the original UVic model, except in the Indian Ocean (Fig. 12). The relatively low NPP in the original UVic model appears to be connected to a correspondingly low NCP (9.3 Pg C yr −1 ), which is close to previous model predictions (clustering around 10 Pg C yr −1 ; Laws et al., 2000;Dunne et al., 2005;DeVries and Weber, 2017). The high (overestimated) NPP in OPEM and OPEM-H is associated with much higher NCP predictions (12.9 and 13.0 Pg C yr −1 , respectively), which are much closer to the satellite-based estimate of 13.5 Pg C yr −1 (Fig. 12) based on Li and Cassar (2016).

N 2 fixation and diazotrophs
N 2 -fixation rates are shown in Fig. 13. Unfortunately, our model simulations differ most strongly in the Indian Ocean, for which no data exist in the MAREDAT database of Luo et al. (2012). One of the problems we face regarding N 2 fixation is that our UVic simulations do not include benthic denitrification and hence miss the dominant oceanic fixed-N loss term (e.g. Gruber, 2004;Wang et al., 2019). Since we have run the models into steady state, N 2 fixation must balance denitrification, which in our case occurs only in the water column. Thus, our UVic simulations cannot be expected to generate realistic global rates of N 2 fixation unless watercolumn denitrification is strongly overestimated. Accordingly, our predicted N 2 -fixation rates (53.9 Tg N yr −1 in the original UVic, 71.4 Tg N yr −1 in OPEM, and 69.5 Tg N yr −1 in OPEM-H; Fig. 13) are much closer to current estimates of water-column denitrification than total N 2 fixation (≈ 70 vs. ≈ 160 Tg N yr −1 ; Wang et al., 2019). Another major difference is the much larger relative contribution of Northern Hemisphere N 2 fixation in OPEM and OPEM-H compared to the original UVic (Fig. 13, top right). The North Atlantic contributes only 4 % in the original UVic, but the 23 % and 24 % contributions in OPEM and OPEM-H, respectively, are closer to the observation-based estimate of 23 % reported by Landolfi et al. (2018), for the data from Luo et al. (2012), than any other model mentioned there.
Both OPEM and OPEM-H predict less N 2 fixation than the original UVic model in the Indian Ocean, which explains (at least partly) the differences in N* there (Fig. 5). OPEM and OPEM-H have no N 2 fixation in the northern Indian Ocean, which is an area of intense diazotrophy in the original UVic, owing the presence of diazotrophs in the original UVic and their absence in OPEM and OPEM-H in this region (Fig. 15). Other models, for example, the one of Monteiro et al. (2011), also produce high rates of N 2 fixation in the northern Indian Ocean, similar to the distribution simulated by the original UVic. In contrast, Löscher et al. (2020) recently found no evidence for significant N 2 fixation in the Bay of Bengal. Whether the qualitative change towards very little N 2 fixation also in other parts of the Indian Ocean, as simulated by both OPEM and OPEM-H, is a qualitative improvement in the representation of N 2 fixation by biogeochemical ocean models, remains to be seen. OPEM-H predicts a wider geographical range for N 2 fixation than the other UVic configurations, due to Houlton et al.'s 2008 temperature function for diazotrophy, now occurring in a few spots north of 40 • N (Fig. 13) (Westberry et al., 2008) with missing data treated as 0. Figure 11. Annually averaged distribution of surface autotrophic POC estimated from satellite data via the CbPM and predicted from the OPEM and OPEM-H UVic simulations. The contours in the right panels indicate multiples of the zooplankton feeding threshold ( th , Eq. C18); i.e. a value of 1 means that effective autotrophic POC (defined as φ phy C phy + φ dia C dia ) is equal to th . The lower left panel illustrates the relation between relative errors in vertically integrated NPP and surface autotrophic POC (δNPP and δPOC, respectively) with respect to the CbPM data. The relative errors δx are defined as δx = x model /x CbPM − 1. The solid lines show the regressions forced through the origin. The slopes of these lines are 1.064 ± 0.059 (R 2 = 0.05, OPEM) and 1.028 ± 0.024 (R 2 = 0.25, OPEM-H). The satellite-based CbPM estimate is the average for 1998-2007 (Westberry et al., 2008) with missing data treated as 0.
compared to the UVic temperature function for diazotrophs at high temperatures appears to be rather small, mostly restricted to the tropics (Fig. 13, top right) but may be the main reason for the slightly lower global N 2 fixation in OPEM-H compared to OPEM. Thus, widening the temperature range of N 2 fixation as in OPEM-H could well be a prerequisite for a more realistic representation of diazotrophy.
Comparing the distributions of simulated N* and N 2 fixation reveals a positive relation with N 2 fixation, which occurs mostly in regions with N* > 0 (Fig. 13). This pattern is very different from that in the analysis of Deutsch et al. (2007), who assumed a high PO 3− 4 demand of diazotrophs, whereas our model does not make this assumption and actually predicts that N 2 fixation can greatly increase the competitive ability of diazotrophs at low PO 3− 4 concentrations (Pahlow et al., 2013). Thus, in our models, the rise in N* due to N 2 fixation does not destroy the niche of the diazotrophs but rather creates an environment in which their ability to utilise very low PO 3− 4 concentrations allows them to persist. This ability derives from the absence of N limitation in the original UVic and from the additional N allocation towards P uptake in OPEM and OPEM-H. Pahlow et al. (2013) suggested that the coexistence of ordinary and diazotrophic phytoplankton should result in a roughly inverse relation between NO − 3 /PO 3− 4 and [PO 3− 4 ] due to the high competitive ability of diazotrophs under low NO − 3 and in particular PO 3− 4 concentrations. This inverse relation implies that N 2 fixation can occur under high NO − 3 /PO 3− 4 ratios only when [PO 3− 4 ] is low and is indeed observed in data from WOCE section A05 in the subtropical North Atlantic (Millero et al., 2000) and predicted by OPEM and OPEM-H, but not by the original UVic, for the same region (Fig. 14a). The patterns for the global surface ocean reveal a similar inverse relation for the original UVic, albeit much less constrained than for OPEM (Fig. 14b, c). In all cases, the patterns for locations with N 2 fixation are very different from those for all regions (green and blue dots in Fig. 14b, c). Whereas the pattern for the original UVic appears more similar to the pattern in the data from Luo et al. (2012) corresponding to total N 2 fixation, except where both NO − 3 and PO 3− 4 are very low (Fig. 14b), the pattern in OPEM is closer to that where N 2 fixation by Trichodesmium occurs (Fig. 14c). Thus, the representation of diazotrophy still appears to warrant further investigation. While none of our UVic configurations can explain N 2 fixation occurring at very low NO − 3 and PO 3− 4 concentrations (Fig. 14b), the physiology of N 2 fixation clearly has a strong influence on NO − 3 /PO 3− 4 and hence N* patterns, as revealed, in particu-  (Fig. 14c).
Contrary to the original UVic model, we do not apply any explicit growth-rate reduction to the diazotrophs in our OPEM simulations, but we assign a lower nutrient affinity and a higher Fe half-saturation concentration to diazotrophs (k Fe,dia > k Fe,phy , whereas k Fe,dia < k Fe,phy in the original UVic), and the model calibration yielded a higher value of the prey-capture coefficients for diazotrophs (Table 2; see also Part 2; Chien et al., 2020). Both OPEM and OPEM-H have a similar phytoplankton biomass and distribution (Fig. 15). Phytoplankton biomass (not Chl; see Fig. 9) is much more evenly distributed, and the integrated biomass is about 2.3 times as large as in the original UVic model.
Diazotrophs are implemented as facultative and their biomass is distributed very differently in all three UVic simulations (Fig. 15). In the original UVic and OPEM, the diazotroph distribution roughly matches that of N 2 fixation, whereas prominent diazotroph biomass appears at high latitudes, even in the Arctic and Antarctic oceans, in OPEM-H, mostly unassociated with N 2 fixation (see also Fig. 13). In fact, non-N 2 -fixing diazotrophs are responsible for the improved representation of Chl, NPP, and NCP in the Arctic when compared to satellite-based estimates (Figs. 9-12) in OPEM-H but also for the somewhat higher N* values at high latitudes compared to OPEM (Fig. 5).
The main reason why the facultative diazotrophs can populate the high latitudes in OPEM-H is their higher light affinity (α = 0.5 compared to 0.4 m 2 mol C W −1 (g Chl) −1 d −1 for ordinary phytoplankton), which can overwhelm the effect of the much higher food preference for diazotrophs (compare φ dia and φ phy , Table 2) under light-limited conditions. A high α for diazotrophs was also obtained by Pahlow et al. (2013). In these areas, characterised by low light and high inorganic nutrient availability, the advantage of a higher α more than compensates for the lower nutrient affinity (A 0 ) and higher N demand (Q N 0 ) of the diazotrophs. Our interpretation of this behaviour is that the diazotroph compartment in OPEM-H actually represents two functional groups: one occurring at low latitudes, representing what we usually associate with facultative diazotrophs, and one occurring at high latitudes, representing non-N 2 -fixing species adapted to low light and long periods of darkness. While the diazotrophs in OPEM-H are able to fix N 2 in high latitudes, they do not because it would reduce their net growth rate compared to utilising nitrate. The (facultative) diazotrophs occur mostly where their realised net relative growth rate exceeds that of ordinary phytoplankton ( µ > 0, µ = µ dia − µ phy ) for OPEM and OPEM-H but not for the original UVic (Fig. 13). The main reason for this discrepancy in the original UVic is the much lower food preference for diazotrophs (0.1) compared to ordinary phytoplankton (0.3) in this configuration, which partly decouples the competitive balance between the two autotrophic groups from µ.
While the occurrence of diazotrophs in the Arctic (mostly without N 2 fixation) appears helpful in view of high-latitude NPP, it is also responsible for the overestimation of N* there (Fig. 5), due to their high N : P ratios. The C : N : P of or-dinary phytoplankton in the Arctic (not shown) is close to Redfield proportions in OPEM, but this simulation fails to generate any appreciable NPP there. Although it might also be possible to explain the low N* in the Arctic with a high N : P ratio in Arctic zooplankton, we are not aware of any indication of this. Hence, phytoplankton in the Arctic appears    Table 3. Log-averaged C : N and C : P ratios for the depth ranges of the upper two layers in the UVic model. to have a low N : P ratio and cannot be represented by our facultative diazotrophs. Low phytoplankton N : P utilisation ratios in the Arctic have been reported by, e.g. Mills et al. (2015), who also inferred high rates of benthic denitrification there. Since we have no benthic denitrification and almost no N 2 fixation in our UVic simulations, it is clear that the stoichiometric imbalance between phytoplankton and zooplankton strongly affect surface N* in the Arctic. Thus, the most likely explanation of the low Arctic N* may be the combination of benthic denitrification and phytoplankton communities dominated by species with high light affinity and a low N subsistence quota.

C : N : P ratios
Simulated log-averaged particulate (i.e. the sum of phytoplankton, diazotrophs, zooplankton, and detritus) C : N and C : P ratios of both OPEM and OPEM-H are well above the canonical Redfield ratios (C : N = 6.625 mol mol −1 and C : P = 106 mol mol −1 ; Table 3) in the topmost two layers. Both simulations tend to overestimate C : N ratios in the surface layer and underestimate C : P compared to observations compiled by Martiny et al. (2014), though not as much as the uniform Redfield C : P ratio employed in the original UVic model. While the data indicate increasing C : P with depth, it is lower in the second compared to the first layer in OPEM and OPEM-H (Table 3). The increasing C : P in the data may be indicative of preferential remineralisation of P relative to C and N (e.g. Letscher and Moore, 2015), which is absent in the current UVic configurations. The decline of C : N and C : P with depth in UVic is the result of primary production with lower light and greater nutrient availability in the second layer. This effect may well be too strong in UVic, due to its coarse vertical resolution, enforcing a homogeneous vertical distribution of all biological tracers within the upper 50 m. The latitudinal patterns of the particulate C : N and C : P ratios are shown in Fig. 16. Interestingly, the simulated C : N ratios are closer to the observations in the Southern Hemisphere, while the simulated C : P ratios match better in the Northern Hemisphere. C : N ratios in the surface layer appear too high throughout, whereas those in the second layer are a lot closer to the observations, whereas C : P ratios seem to match similarly in both layers (Table 3 and Fig. 16).
Patterns of C : N ratios mirror the relation between light and nutrient limitation in our OPEM simulations, with high C : N ratios indicating strong nutrient limitation, which is also generally observed in phytoplankton culture experiments (Pahlow et al., 2013). Thus, one possible explanation for the too-high particulate C : N ratios in the surface layer could be that too little nutrients reach the surface ocean at subtropical northern latitudes. This is consistent with toolow rates of NPP being predicted around 20 • N (Fig. 10), where the overestimation in surface C : N ratios is strongest (Fig. 16). The lower C : N ratios at high latitudes (60 • S and 60 • N) in OPEM-H reflect the dominance of (non-N 2 -fixing) diazotrophs there in this simulation.
The relatively high C : N ratios throughout most of the surface layer also largely explain the lower export efficiency, as indicated by the much higher NPP estimate (Fig. 10) relative to NCP (Fig. 12) in OPEM and OPEM-H compared to the original UVic. Since the average particulate C : N and C : P ratios are much greater in OPEM and OPEM-H than the (Redfield) C : N and C : P ratios of the zooplankton, the excess C is released in dissolved form (as CO 2 ) by the zooplankton according to Eq. (C16). Thus, consumption of particles with elevated C : N and/or C : P relative to the zooplankton lowers the export efficiency. While particulate C : P agrees much better with the observations than C : N, it is still on average well above the (Redfield) C : P ratio of the zooplankton, which implies that a better match of surface particulate C : N alone might not reconcile the relative magnitudes of NPP and NCP in OPEM and OPEM-H with the satellitederived estimates. Both the high surface C : N and low P : C in midlatitude regions might result from the underestimation of N 2 fixation, due to the lack of benthic denitrification. Enhanced N 2 fixation would add fixed N to the surface ocean, partly releasing phytoplankton from N limitation and intensifying P limitation, and could thus bring C : N and C : P ratios closer to the observations. Further promising approaches in this respect may be the consideration of preferential remineralisation, which could allow enhanced N assimilation due to additional P availability, or allowing for variable stoichiometry in zooplankton (e.g. Talmy et al., 2014).
The C : N and C : P ratios of sinking particles (detritus) in OPEM and OPEM-H are greater than those of total particulate matter (Fig. 7), because the C : N : P ratio of zooplankton is 106 : 16 : 1 but that of its food is larger. Zooplankton respire the excess C in the food, thereby reducing the average particulate C : N : P, whereas the detritus pool is fed not only by zooplankton egestion but also by the phytoplankton and diazotroph mortality terms with relatively high C : N : P ratios. The magnitude of this effect is modulated by the zooplankton assimilation efficiency (E zoo ) as this determines the fraction of particulate egestion. In regions with high E zoo ≈ 1 (Fig. 17), almost no particles are egested, whereas for E zoo ≈ 0.5 about half of the ingested food is lost to detritus. The relatively low assimilation efficiencies in the Arctic between 90 • E and 120 • W in OPEM-H compared to OPEM in Fig. 17 result from the availability of food, as OPEM-H is the only simulation with any appreciable Figure 16. Zonally averaged particulate C : N and C : P ratios for the depth ranges of the two topmost layers of UVic for 5 • latitude bands. Lines are predictions from the OPEM and OPEM-H simulations, and circles represent data from Martiny et al. (2014). POC < 0.01 mmol m −3 , PON < 1 µmol m −3 , and POP < 0.1 µmol m −3 were removed from the observations prior to calculating the ratios. Observed ratios were mapped onto the UVic grid by taking the median of all available data for each grid cell and then zonal log averages were calculated. NPP (Fig. 10) and hence biomass in this region (Fig. 15), and E zoo is inversely related to ingestion in OPEM and OPEM-H. Food availability exceeds the zooplankton feeding threshold in this region only for OPEM-H (contours in Fig. 17).

Conclusions
The above description of the model behaviour highlights some of the improvements of our optimality-based (OPEM, OPEM-H) compared to the original biogeochemistry in the UVic model. Some of these may also be possible with the original UVic with improved parameters, e.g. the deep-ocean N : P distribution (Fig. 8) or a better global NCP (Fig. 12), as these vary strongly among our different parameter sets tested during the calibration process of OPEM and OPEM-H (Chien et al., 2020). Others are simply impossible to reach with a fixed-stoichiometry model, e.g. the distribution of C : N and C : P ratios in particulate matter (Fig. 16). Apparently, our optimality-based biology has a certain internal rigidity (Krishna et al., 2019), preventing us from tuning the OPEM simulations so that, e.g. global NPP, NCP, and N 2fixation distributions can simultaneously be reproduced very well with the same parameter settings. We thus try to use the resulting, and often systematic, model-data discrepancies in the behaviour of OPEM and OPEM-H as a magnifying glass on model deficiencies to identify avenues for future biogeochemical model development.
A similar difference in low-latitude NPP pattern as between the CbPM and OPEM predictions can be seen on the Ocean Productivity website (O'Malley, 2017) as resulting from the use of a polynomial (Behrenfeld and Falkowski, 1997) vs. an exponential (Eppley, 1972) temperature func-tion, as also applied in the UVic model. The CbPM does not have a direct temperature dependence, and Taucher and Oschlies (2011) found that omission of direct temperature effects on biotic processes did not reduce the ability of the UVic model to reproduce observed tracer distributions. Mechanistically, temperature effects might well be subdued under light-limiting conditions, since photochemical reactions are less temperature sensitive than most other biochemical processes. The wider temperature range for diazotrophy in OPEM-H allows for N 2 fixation north of 40 • N, which has been observed recently in the western North Atlantic (Mulholland et al., 2019). Therefore, investigating temperature effects could be a promising approach towards more realistic NPP and N 2 -fixation rates.
Environmental constraints on diazotrophy in our UVic simulations suffer from the absence of benthic denitrification, as mentioned above. In addition, preferential P remineralisation could be important for a better representation of N 2 fixation (Monteiro and Follows, 2012). For example, Fernández-Castro et al. (2016) found that preferential P remineralisation is essential for reproducing observed N 2fixation rates at Bermuda Atlantic Time-series Site (BATS), particularly when atmospheric deposition of fixed N is also considered. Thus, preferential P remineralisation may not only be important for improving the vertical distribution of particulate C : P (Fig. 16) but also for the simulation of diazotrophy. According to Fernández-Castro et al. (2016), this phenomenon could also be a prerequisite for realistically accounting for the effects of atmospheric deposition of nutrients into the surface ocean.
The similarity in the spatial patterns of NPP and surface autotrophic POC, also as they compare to satellite-derived estimates, suggests that the growth of primary producers might be relatively well described but further developments in the representation of top-down control by zooplankton, but also by higher trophic levels or viruses, may be another promising route towards a better resolution of plankton biogeochemical processes.
Besides temperature and top-down effects, the distributions of NPP and particulate C : N ratios are also strongly affected by light and nutrient affinity (model parameters α and A 0 ). The use of fixed settings in these parameters may be responsible for both overestimating NPP at low latitudes ( Fig. 10) and preventing ordinary phytoplankton from growing in the Arctic Ocean (Fig. 15), as indicated by the growth of facultative (but mostly non-N 2 -fixing) diazotrophs there in the OPEM-H simulation. The biotic compartments of the OPEM configurations have been shown to match the observed behaviour of at least some phytoplankton and zooplankton species (Pahlow and Prowe, 2010;Pahlow et al., 2013). Thus, the failure to obtain a better fit to the observed NPP distribution may reflect a certain rigidity, brought about by attempting to represent plankton communities by a globally uniform parameter set, i.e. one and the same combination of one phytoplankton, one diazotroph, and one zooplankton species. As mentioned above, Arteaga et al. (2016) achieved a strong improvement in model behaviour by replacing α and A 0 with a trade-off represented by opposite linear functions of light and nutrient limitation. Since our cost function does not appear to be very sensitive to α, we interpret these findings as an indication that the regional variability of α may be more important for the model behaviour than its global average. Similar formulations could be introduced, e.g. to represent species sorting (Norberg, 2004;Smith et al., 2016), possibly responsible for regional and local variations in α and A 0 . Whether variations in these two parameters suffice, e.g. to explain the low N* in the Arctic, remains to be seen. The approach might have to be extended to further parameters for a more realistic representation of different phytoplankton and zooplankton communities Su et al., 2018). Nevertheless, it is clear from Fig. 5 that N* in the surface ocean is very sensitive to plankton physiology (subsistence quotas), which could greatly complicate inferring regional balances of N 2 fixation and denitrification from N* or similar quantities (e.g. Mills et al., 2015).

Appendix A: Bug fixes applied to all configurations
UVic has already contained code intended to reduce the occurrence of negative concentrations by setting all sink terms to 0 once a concentration drops below a certain threshold. This mechanism was made partly ineffective, however, by passing positive values to the biogeochemical subroutine (npzd_src), even when the actual tracer concentration was negative, so that the negative concentration was not detected, or too late, and sink terms could still apply. This was corrected by passing the actual tracer values to the npzd_src subroutine.
The dynamic Fe model (Nickelsen et al., 2015) injects atmospheric Fe deposition directly into the surface layer, which we consider as a bug, as this bypasses the surface-flux mechanism built into UVic. Correcting this bug also reduces the occurrence of negative Fe concentrations.

Appendix B: Preventing negative concentrations in OPEM
One of the main problems for implementing our variablestoichiometry formulation in UVic's finite-difference code is the occurrence of negative concentrations in UVic. Negative concentrations occur predominantly as a result of the semiimplicit vertical mixing scheme when applied to steep vertical gradients (with smaller contributions arising from advection, the explicit isopycnal mixing scheme, and high-latitude filtering), as revealed by detailed inspection of the model's behaviour. Since the vertical gradients related to the biotic tracers in OPEM are generally much steeper, at least in the upper three layers of the ocean grid, negative concentrations can become much larger and more widespread in OPEM than in the original UVic. Inside its biogeochemical module, UVic deals with negative concentrations by preventing, at every time step and in every grid box, any fluxes out of negative tracer compartments, as mentioned above. UVic also applies a flux-corrected central-differencing scheme for tracer advection (FCT, applied here also in the vertical) in order to prevent generation of negative concentrations. Negative concentrations are also generated in the main biogeochemical module of UVic (subroutine npzd_src), due to the long time steps (we use 0.5 times the physical time step of 30 h and, if this would generate negative tracer concentrations, subcycle with 0.25 times the physical time step) and the Euler scheme used for calculating the sources-minus-sinks terms.
For many cases (parameter settings), phytoplankton and/or diazotrophs can end up negative everywhere in OPEM, compromising our calibration procedure, which depends on the reliability of simultaneous evaluation of simulation ensembles (see Sect. 2.4 and Part 2; Chien et al., 2020). We have addressed the problem in OPEM by limiting the biological tracer fluxes of the subcycled biological time step at every grid box, so that not more than 90 % of any tracer is removed within any grid box during one time step. In order to counter the generation of negative concentrations by advection and vertical mixing, we also modify the physical transport of all particulate tracers and dissolved iron as follows: the sourcesminus-sinks terms of the biogeochemical module are applied before calculating advective and diffusive fluxes, so that diffusion is the only remaining source of negative concentrations. In all cases where the sum of all diffusive fluxes (D) would remove more of a tracer than is present in a grid cell after applying advective fluxes (T ), we calculate a correction factor, f D = −T /(D × t), where t is the time step, which is then multiplied with all outward diffusive fluxes to ensure a non-negative tracer concentration. This flux limitation does not affect tracer conservation. Since limiting the flux out of one grid cell reduces the flux into the neighbouring cell, this procedure is applied recursively until nonnegative concentrations are guaranteed everywhere. Whenever high-latitude filtering (Kvale et al., 2017) results in negative concentrations, we multiply positive changes T + by a factor f filt = T filt <0.1T (0.1T − T filt )/ T + and hence allow filtering-induced reductions by at most 90 %, where T filt is the (possibly negative) result of the high-latitude filter.
Appendix C: Optimality-based process descriptions

C1 Phytoplankton and diazotrophs
Please note that we omit the subscripts phy and dia in this subsection.

C1.1 Optimal-growth regulation
Our optimality-based formulations use allocation factors to allocate energy and other resources between light harvesting and nutrient acquisition at each grid point and time step, such that net growth of phytoplankton is maximised. The rates of net relative growth (µ), nutrient uptake (V N and V P ), and N 2 fixation (F N ) in the OGM (optimal-growth model) are given by the optimality-based chain model of Pahlow et al. (2013), modified here to allow for temperature dependence and Fe limitation and to avoid outgrowing the P subsistence quota during transition towards P limitation. Net relative growth rate is the difference between C fixation (V C ) and the sum of respiration (R) and extra DIC release (r DIC ; see below) to prevent outgrowing the P subsistence quota. The chain model idea is based on the roles of N and P in a phytoplankton cell, where P is mainly needed for N assimilation and N drives all other biochemical rates (Ågren, 2004), including growth. Thus, the optimal regulation can be described in terms of two conceptual levels, with the lower level consisting of the nutrient-uptake apparatus and the chloroplast, and the upper level being the whole cell. Within the nutrient-uptake apparatus, cellular N is allocated between N and P uptake so as to maximise N assimilation (see Sect. C1.2 below). Since the role of P is restricted to the nutrient-uptake apparatus in this model, we can ignore P in the formulation of the optimal allocation scheme at the whole-cell level: We collect all N-independent gain and loss terms in µ * : where the allocation factors f C and f V ensure optimal allocation of cellular N between C fixation and nutrient uptake, respectively (see Pahlow et al., 2013, for derivation), f (T ) is temperature dependence, L day is day length, V C 0 the temperature-and Fe-dependent maximum potential rate for C processing, α the light-absorption coefficient (light affinity),θ the Chl : C ratio of the chloroplast, I irradiance, ζ Chl and ζ N the costs of Chl synthesis and N assimilation, R Chl the cost of Chl synthesis and maintenance, R Chl M the cost of Chl maintenance, and S I the depth-and time-averaged light saturation of the photosynthetic apparatus.θ is obtained as its balanced-growth optimum according to Pahlow et al. (2013), modified for numerical stability: is the minimum light intensity for photosynthesis (see Pahlow et al., 2013), and W 0 is an approximation of the 0 branch of Lambert's W function (W 0 ) for very large arguments x, applied here only for x > e 700 (whence the relative error < 10 −7 ) to prevent numeric overflows in Fortran's exp function. S I is calculated assuming a triangular light cycle and uniform light attenuation within a grid cell: where I 0 and I ( z) are the mean daytime light intensities at the top and bottom of the current grid cell of height z, is the light-attenuation coefficient, Ei is the exponentialintegral function, and the factor of 2 converts the mean to the maximum irradiance in the triangular light cycle. As in the original UVic code, we assume that ∝ N phy + N dia + absorption by seawater, since chlorophyll is not a tracer. Equations (C6) and (C7) apply only for I > I min . Thus, for I 0 > I min > I ( z), Eq. (C7) is applied to the part of the grid cell where I > I min and then multiplied with z * / z, where I ( z * ) = I min . In effect, this means that S I > 0 occurs only in the upper 240 m (the top three layers) of the UVic grid.

C1.2 Optimal uptake kinetics
DIN and DIP uptake and N 2 fixation are defined as products of allocation factors, setting the size of the respective cellular compartment and the rate of uptake normalised to the size of that compartment ( V ). V is defined in Eq. (C9) via optimal uptake kinetics (Pahlow, 2005;Smith et al., 2009). The nutrient-uptake compartment, responsible for DIN and DIP uptake and N 2 fixation, contains fraction f V of the cellular N resources, of which fraction f N is available for DIN uptake, leaving f V (1 − f N ) for DIP uptake: where A 0 is nutrient affinity and f F the allocation for N 2 fixation within the nutrient-uptake compartment. The allocation factor f F is implemented as a switch, so that the facultative diazotrophs either fix N 2 or utilise DIN (see Pahlow et al., 2013, for derivation). The dependence of V max and F N on Q P introduces a chain of limitations, where the P quota limits N uptake and N limits all other processes. Extra DIC release (r DIC ) during transition towards severe P limitation prevents outgrowing of the P subsistence quota (Q P 0 ): where the first term limits r DIC to conditions of declining Q P and the second term states that r DIC > 0 occurs only for Q P < 2Q P 0 . Equation (C11) is an admittedly rather arbitrary measure to stabilise the OGM, but it did result in reasonable rates of DOC production in a previous study (Fernández-Castro et al., 2016).

C1.3 Temperature and Fe limitation
Temperature and Fe limitation are implemented by V C 0 (T ) = V N 0 (T ) = f p (T ) · S Fe · V 0 , V P 0 (T ) = f p (T ) · V 0 , F N 0 (T ) = f nfix (T ) · S Fe · F 0 p ∈ {phy, dia} (C12) λ phy = λ 0,phy · f phy (T ) M dia = M 0,dia · f dia (T ), where V 0 is the potential-rate parameter, F 0 the potential rate of N 2 fixation, f p (T ) the group-specific temperature dependence of nutrient uptake and photosynthesis, f dia (T ) the temperature dependence of N 2 fixation, and S Fe the Fe limitation term.

C2 Zooplankton
Net growth (µ zoo ) is described in terms of total (A t ; see Eq. C19 below) and foraging activity (A f ), and corrected for r Q : µ zoo = (E zoo · g zoo − R * zoo ) · r Q , g zoo = A f · S g , where C zoo = 6.625·N zoo and N zoo are zooplankton POC and PON, µ zoo net relative growth rate, G N zoo predation on zooplankton, M zoo (quadratic) mortality, Q N zoo N : C ratio, g zoo relative ingestion rate, E zoo and E max actual and maximal assimilation efficiency, X C zoo egestion, R * zoo and R C zoo minimal (uncorrected for r Q ) and actual respiration, R n zoo metabolic N and P losses, β digestion coefficient, c a and c f cost of assimilation and foraging coefficients, and R M zoo maintenance respiration. The same relation between dissolved and particulate losses applies for N and P as for C in Eq. (C17). Equations (C14)-(C16) define the benefits (g zoo ) and costs (E zoo and R * zoo ) of foraging, whence the optimal foraging activity is obtained as where W −1 is the −1 branch of Lambert's W function and th is the feeding threshold. A t is a function of the maximal ingestion rate (g max ) and temperature: The predation rates for individual prey types are G C p = φ p C p C · g zoo · C zoo , C zoo = N zoo Q N zoo , G N p = G C p · Q N p , G P p = G C p · Q P p , p ∈{phy, dia, det, zoo}. (C20) Equations (4) and (C14)-(C17) stipulate that most of the excess C, N, or P rejected to maintain homeostasis is released in dissolved inorganic form (see Eqs. C14 and C16). This is because the actual growth rate µ zoo is obtained as the product of r Q and the potential growth rate, i.e. that obtained for food with the same stoichiometry as the zooplankton in Eq. (C14), and respiration R C zoo is then derived from µ zoo in Eq. (C16), whereas egestion X C zoo is not affected by r Q in Eq. (C14). Since the relation of dissolved and particulate N and P losses follows that for C (X n zoo in Eq. C14), a stoichiometric imbalance between zooplankton and its food increases dissolved losses for N and P as well.
Code availability. The University of Victoria Earth System Climate Model version 2.9 is available at http://www.climate.uvic. ca/model/ (last access: 11 March 2020) ("updated model" version). The code for the original model and OPEM is available at https://doi.org/10.3289/SW_1_2020 (Pahlow, 2020). The instructions needed to reproduce the model results described in this article are in the Supplement.
Supplement. The supplement related to this article is available online at: https://doi.org/10.5194/gmd-13-4663-2020-supplement. DIN : DIP gradients in their model simulations. The manuscript benefitted strongly from the reviews by David Talmy, Emily Zakem, an anonymous referee, and Andrew Yool.
Financial support. Markus Pahlow was supported by Deutsche Forschungsgemeinschaft (DFG) by SFB754 (Sonderforschungsbereich 754 "Climate-Biogeochemistry Interactions in the Tropical Ocean"; http://www.sfb754.de, last access: 8 September 2020) and as part of Priority Programme 1704 (DynaTrait). Markus Pahlow and Chia-Te Chien were supported by the BMBF-funded Deutsches Zentrum für Luft-und Raumfahrt (PalMod). Lionel A. Arteaga was partially supported by the 2014 miniproposal award from the Integrated School of Ocean Sciences (ISOS) Kiel and by NASA (award no. NNX17AI73G).
The article processing charges for this open-access publication were covered by a Research Centre of the Helmholtz Association.
Review statement. This paper was edited by Andrew Yool and reviewed by David Talmy, Emily Zakem, and one anonymous referee.