Optimality-Based Non-Redfield Plankton-Ecosystem Model (OPEM v1.0) in the UVic-ESCM 2.9. Part I: 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.0 with variable C:N:P stoichiometry in the University of Victoria ESM (UVic, Eby et al., 2009; Weaver et al., 2001) and the behaviour 5 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 fixed-stoichiometry 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 NO3 – :PO4 – ratio, and partially explains observed patterns of particulate C:N:P 10 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 15 may be represented by locally-varying parameters based on suitable trade-offs. The similarity in the overestimation of NPP and surface autotrophic 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* (NO3 – − 16 ·PO4 – + 2.9mmol m−3), it still fails to reproduce observed N* in the Arctic, 20 possibly related to a mis-representation 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.

1 Introduction 25 Earth system models (ESMs) are routinely used for simulating both the possible future development and the past of our climate system (e.g., IPCC, 2013;Hülse et al., 2017;Keller et al., 2018;Park et al., 2019). While different ESMs agree to some extent in their predictions, they usually also encompass 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.,30 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). 35 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 40 variable C:N:P in phytoplankton in CMIP5 (Bopp et al., 2013) and CMIP6 (Arora et al., 2019) 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 half-saturation concentration of nitrate use varies systematically with nitrate concentration and suggested that optimal uptake kinetics (Pahlow, 2005) may be more 45 appropriate than the commonly-used Michaelis-Menten kinetics for simulating phytoplankton nutrient uptake. Zooplankton foraging behaviour can be characterized 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 50 (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 requirements for light harvesting and nutrient acquisition in phytoplankton and for foraging and 55 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 60 shown their ability to describe ecosystem behaviour in 0D and 1D modelling studies (e.g., Fernández-Castro et al., 2016;Su et al., 2018), and to predict patterns of phytoplankton nutrient and light colimitation 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 3D 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 65 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 in the following, Eby et al., 2009;Weaver et al., 2001). Owing to its coarse spatiotemporal resolution, UVic is a practical choice when working on long time scales (e.g., Niemeyer et al., 2017) and/or when many 70 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 (nutrients-phytoplankton-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 75 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 optimality-based 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 II (Chien et al., 2020). 80 2 Optimality-based plankton in the UVic model The 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 85 UVic, except that we use constant half-saturation 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 90 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 Foraging Assimilation 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 (panel 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 N2 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 (panel C) is optimised by balancing costs and benefits of allocating total activity (At) between foraging activity (Af) and assimilation activity (At − Af). Both foraging and assimilation incur energy costs (cf and ca, respectively) fuelled by respiration (R). Increasing ingestion (g) reduces assimilation efficiency (E ≤ Emax), causing more particulate egestion (X).

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 spatio-temporal resolution of UVic and augmented with temperature and iron effects (see equations 95 provided below). Owing to the relatively long time step, the model does not resolve the dynamics of photo-acclimation and we therefore describe the Chl:C ratio of the chloroplast by its balanced-growth optimum. 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 C, N, P (POC, PON, 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 100 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 105 OGM has only 8 parameters (maximum rate V 0 , nutrient affinity A 0 , costs of N assimilation ζ 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 6 parameters to be calibrated.    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 FCT technique (flux-corrected transport, 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 115 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: 6 where C p , N p , P p are POC, PON, 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 biomass-normalised N and P cell quotas (N:C and P:C ratios). The last term in (3) 125 accounts for the subsistence amounts of N and P implicitly contained in C p and subtracted from δn p via (1). Leakage is the fast-recycling term parametrising the microbial loop (Keller et al., 2012). Definitions for all terms in Eqs. (2) and (3)

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, 140 it can be expressed as a function of the maximal ingestion rate, which is routinely determined in feeding experiments, and temperature (see Eq. C18 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 the original UVic, but since two 145 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-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 150 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 155 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. C15 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 are the prey-specific capture coefficients. 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 170 temperature-dependent remineralisation. We consider separate C, N, and P tracers for detritus: where ν det is the detritus remineralization 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 175 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 . Eq. (6) does not include gains and losses 180 from sinking detritus particles. Detritus sinking speed v sink increases with depth, reflecting the gradual disappearance of smaller particles during sinking, according to 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).

185
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 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 200 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 4 3 -+ 2.9 mmol m −3 (Gruber and Sarmiento, 1997;Mills et al., 2015), and global N 2 -fixation rates and distributions within current observational 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.
205 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. Interestingly, these differences cannot be seen in the O 2 distribution at 300 m, the depth of the OMZs, which is very similar 230 in the Indian Ocean and eastern tropical Pacific among all our UVic simulations (Fig. 6), indicating that the carbon export and subsequent remineralization is 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 235 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 240 (Fig. 13 bottom) 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 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 245 models without variable stoichiometry can (e.g., Kriest and Oschlies, 2015). Thus, reproducing the deep DIN:DIP distribution appears to require mostly 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 deep ocean interior (Fig. 5). We interpret the surface N* distribution outside the deep-water formation regions as a consequence, rather than a cause the of deep-ocean nutrient distributions, however. 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).
Global net primary production is defined here as where µ is the net relative growth rate and λ the leakage rate representing fast remineralisation in UVic. 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)  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 Comparing the patterns in NPP and surface autotrophic POC (Figs. 10 and 11) suggests a spatial correlation between de-285 viations in these two quantities (Fig. 11, lower left). Thus, some of the NPP overestimate could result from an overestimate  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 290 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 295 reasonable compared to the satellite-derive 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 300 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 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
305 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  the UVic temperature function for diazotrophs at high temperatures appears to be rather small, 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 330 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 4 3demand 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 4 3concentrations (Pahlow et al., 2013). Thus, in our 335 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 4 3concentrations 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. 3− ] 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 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 345 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 4 3are 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 4 3concentrations (Fig. 14B), the physiology of N 2 fixation clearly has a strong influence on NO 3 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 355 for diazotrophs (Table 2, see also Part II, 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 360 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 affin-365 ity (α = 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 370 functional groups, one occurring in 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. 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 orignal UVic is the much lower food preference for diazotrophs (0.1) compared to ordinary phytoplankton (0.3) in this 375 configuration, which partly decouples the competitive balance between the two autotrophic groups from ∆µ.
While the occurrence of diazotrophs in the Arctic appears helpful in view of high-latitude NPP, they are also responsible for the overestimation of N* there (Fig. 5), owing to their high N:P ratios. The C:N:P of ordinary 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 380 of any indication of this. Hence, phytoplankton in the Arctic appears 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 385 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, 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 390 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 395 layer. This effect may well be too strong in UVic, owing to its coarse vertical resolution, enforcing a homogeneous vertical distribution of all biological tracers within the upper 50 m. Table 3. Log-averaged C:N and C:P ratios for the depth ranges of the upper two layers in the UVic model. 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, 400 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 too low rates of NPP being predicted 405 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 the surface layer also largely explains 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 410 ratios of the zooplankton, the excess C is released in dissolved form (as CO 2 ) by the zooplankton according to Eq. (C15).
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 satellite-derived estimates. Both the high surface C:N and low P:C in mid-415 latitude regions might result from the underestimation of N 2 fixation, owing 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 425 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 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). Redfield 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 calculated.

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 435 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 2 -fixation 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 440 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 function, 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 445 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 2 fixation rates at 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) 455 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 460 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.

465
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 470 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 indicating 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 475 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 485 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  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 Section 2.4 and Part II, Chien et al., 2020). We have addressed the problem in OPEM by limiting the biological tracer fluxes of the sub-cycled 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 510 also modify the physical transport of all particulate tracers and dissolved iron as follows: The sources-minus-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.

515
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 non-negative 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 = Tfilt<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.

520
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 525 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 out-growing 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 dissolved inorganic C (DIC) release (r DIC , see below) to 530 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 Section C1.2 below). Since the role of P is restricted 535 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 545 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. S I is calculated assuming a triangular light cycle and constant 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 exponential-integral function, and the factor 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. Eqs. (C5) and (C6) apply only for I > I min , where I min = ζ Chl R Chl M f (T )/(αL day ) is the minimum light intensity for photosynthesis (see Pahlow et al., 2013). Thus, for I 0 > I min > I(∆z), (C6) is applied to the part 555 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 3 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 normalized to the size of that compartment ( V ). V is defined in Eq. C8 via optimal uptake 560 kinetics (Pahlow, 2005;Smith et al., 2009). The size of 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 utilize 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 570 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 . Eq. (C10) 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 λ phy = λ 0,phy · f phy (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 580 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. C18 below) and foraging activity (A f ), and corrected for r Q : 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, 590 β 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 (C16). Eqs. (C13)-(C15) 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 Lambert's W-function and Π th is the feeding threshold. A t is a function of the maximal ingestion rate (g max ) 595 and temperature: The predation rates for individual prey types are Eqs. (4) and (C13)-(C16) stipulate that most of the excess C, N, or P rejected to maintain homeostasis is released in dissolved 600 inorganic form (see Eqs. C13 and C15). 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. (C13), and respiration R C zoo is then derived from µ zoo in Eq. (C15), whereas egestion X C zoo is not affected by r Q in Eq. (C13). Since the relation of dissolved and particulate N and P losses follows that for C (X n zoo in Eq. C13), a stoichiometric imbalance between zooplankton and its food increases dissolved losses for N and P as well.
The manuscript has benefitted strongly from the reviews by D. Talmy