MESMO 3: Flexible phytoplankton stoichiometry and refractory dissolved organic matter

We describe the third version of Minnesota Earth System Model for Ocean biogeochemistry (MESMO 3), an Earth system model of intermediate complexity, with a dynamical ocean, dynamic–thermodynamic sea ice, and an energy–moisture-balanced atmosphere. A major feature of version 3 is the flexible C : N : P ratio for the three phytoplankton functional types represented in the model. The flexible stoichiometry is based on the power law formulation with environmental dependence on phosphate, nitrate, temperature, and light. Other new features include nitrogen fixation, water column denitrification, oxygen and temperaturedependent organic matter remineralization, and CaCO3 production based on the concept of the residual nitrate potential growth. In addition, we describe the semi-labile and refractory dissolved organic pools of C, N, P, and Fe that can be enabled in MESMO 3 as an optional feature. The refractory dissolved organic matter can be degraded by photodegradation at the surface and hydrothermal vent degradation at the bottom. These improvements provide a basis for using MESMO 3 in further investigations of the global marine carbon cycle to changes in the environmental conditions of the past, present, and future.


Introduction
Here we document the development of the third version of the Minnesota Earth System Model for Ocean biogeochemistry (MESMO 3). As described for the first two versions (Matsumoto et al., 2008, MESMO is based on the non-modular version of the Grid ENabled Integrated Earth (GENIE) system model (Lenton et al., 2006;Ridgwell et al., 2007). The computationally efficient ocean-climate model of Edwards and Marsh (Edwards and Marsh, 2005) forms the core of GENIE's physical model. MESMO is an Earth system model of intermediate complexity (EMIC), which occupies a midpoint in the continuum of climate models that span high-resolution, comprehensive coupled models on one end and box models on the other (Claussen et al., 2002). MESMO has a 3D dynamical ocean model on a 36 × 36 equal-area horizontal grid with 10 • increments in longitude and uniform in the sine of latitude. There are 16 vertical levels. It is coupled to a 2D energy moisture-balanced model of the atmosphere and a 2D dynamic and thermodynamic model of sea ice. Thus, MESMO retains important dynamics that allow for simulations of transient climate change, while still being computationally efficient.
Since the first version, MESMO has continued to be developed chiefly for investigations of ocean biogeochemistry (Table 1). Briefly, in MESMO 1 the main improvements over the predecessor GENIE focused on the biological production and remineralization, as well as on the uptake of natural radiocarbon ( 14 C) and anthropogenic transient tracers (Matsumoto et al., 2008). The net primary production (NPP) in MESMO 1 occurred in the top two vertical levels, representing the surface 100 m, and depended on temperature, nutrients, light, and mixed-layer depth (MLD). The nutrient dependence was based on the Michaelis-Menten uptake kinetics of phosphate (PO 4 ), nitrate (NO 3 ), and aqueous CO 2 . The limiting nutrient was determined by Liebig's rule of the minimum relative to the fixed uptake stoichiometry of C : N : P = 117 : 16 : 1. A single generic phytoplankton functional type (PFT) carried out NPP, which was split between particulate organic matter (POM) and dissolved organic matter (DOM) in a globally constant ratio of 1 : 2. The semi-labile form of the dissolved organic carbon (DOC) was the only form of DOM simulated in MESMO 1. The POM flux across the 100 m level defined the export production. The vertical flux of POM was driven by a fixed rate of sinking and a temperature-dependent, variable remineralization rate.
The main aim of MESMO 2 was a credible representation of the marine silica cycle . To this end, the set of limiting nutrients (P, N, and C) in MESMO 1 was augmented to include iron (Fe) and silicic acid (Si(OH) 4 ) in MESMO 2 ( Table 1). The stable isotope of Si ( 30 Si) was also added as a state variable. The Fe cycle included an aeolian flux of Fe, complexation with organic ligand, and particle scavenging of free Fe. The scavenged Fe that reached the seafloor was removed from the model domain. This burial flux of Fe balanced the aeolian flux at steady state. In addition, a new PFT was added in MESMO 2 chiefly to represent diatoms. This new "large" PFT was limited by Si and characterized by a high maximum growth rate and large halfsaturation constants for the nutrient uptake kinetics. It represented fast and opportunistic phytoplankton that do well under nutrient replete conditions. In comparison, the "small" PFT was characterized by a lower maximum growth rate and smaller half-saturation constants and outperformed the large PFT in oligotrophic subtropical gyres. CaCO 3 production was associated with the "small" PFT in MESMO 2. The addition of Fe, Si, and the large PFT in MESMO 2 allowed it to have a Fe-dependent, variable Si : N uptake ratio (Hutchins and Bruland, 1998;Takeda, 1998), which is critical for simulating important features of the global ocean Si distribution.
MESMO 1 and 2 were assessed and calibrated by multiobjective tuning and extensive model-data comparisons of transient tracers (anthropogenic carbon, CFCs), deep ocean 14 C, and nutrients (Matsumoto et al., 2008. These versions have been employed successfully in a number of studies of global distributions of carbon and carbon isotopes under various conditions of the past, present, and future (Cheng et al., 2018;Lee et al., 2011;Matsumoto et al., 2010Matsumoto et al., , 2020Matsumoto and McNeil, 2012;Matsumoto and Yokoyama, 2013;Sun and Matsumoto, 2010;Tanioka and Matsumoto, 2017;Ushie and Matsumoto, 2012). In addition, MESMO 1 and 2 have participated in model intercomparison projects Cao et al., 2009;Eby et al., 2013;Weaver et al., 2012;Zickfeld et al., 2013).
In this contribution, we describe the third and latest version of MESMO with a number of substantial biogeochemical model modifications and new features that bring MESMO up to date with the evolving and accumulating knowledge of the ocean biogeochemical cycle (Table 1). There is no change in the physical model between MESMO 3 and MESMO 2. The most significant new feature of MESMO 3 over the previous versions is the power law formulation of flexible phytoplankton C : N : P ratio. Other new features include additional PFT diazotrophs that carry out N fixation, water col-umn denitrification, the dependence of organic matter remineralization on the dissolved oxygen (O 2 ) and temperature, and CaCO 3 production based on the concept of the residual nitrate potential growth. In addition, we describe the semilabile DOM for P, N, and Fe (DOP sl , DON sl , and DOFe sl ) and the refractory DOM for C, P, and N (DOC r , DOP r , and DON r ), which can be activated as an optional feature in MESMO 3. Some of these features have been described separately in different publications (Matsumoto et al., 2020;Matsumoto and Tanioka, 2020;Matsumoto, 2017, 2020a). This work consolidates the descriptions of all these features in a single publication.

Model description
Here we present the full set of biogeochemical equations of MESMO 3 and the key model parameters (Table 2). We describe only the biogeochemical source and sink terms and omit the physical (advective and diffusive) transport terms that are calculated by the ocean circulation model. We discuss the production terms first, followed by remineralization terms, and finally the conservation equations that incorporate both terms.

Phytoplankton nutrient uptake
NPP occurs in the top two vertical levels of the ocean domain above the fixed compensation depth (z c ) of 100 m. Key parameter values are given in Table 2a. Nutrient uptake by phytoplankton type i ( i ) depends on the optimal nutrient uptake timescale (τ i ), nutrients, temperature (T ), irradiance (I ), and mixed-layer depth (z ml ): Subscript i refers to PFT (i = 1: eukaryotes; i = 2: cyanobacteria; i = 3: diazotrophs). The nutrient dependence F N,i is given by Liebig's law of minimum combined with Michaelis-Menten uptake kinetics of limiting nutrients: PO 4 , NO 3 , CO 2 (aq), total dissolved iron (sum of free iron and ligand-bound iron: FeT = Fe +FeL), and Si(OH) 4 : where K X,i is the half-saturation concentration of nutrient X for PFT i. Only eukaryotes (i = 1) are limited by Si(OH) 4 . LG stands for large (diatoms) and SM stands for small. MESMO 3 PFTs are as follows: Eu stands for eukaryotes, Cy stands for cyanobacteria, and Dz stands for diazotrophs. OM stands for organic matter. RNPG stands for residual nitrate potential growth. T stands for temperature. PAR stands for photosynthetically available radiation. fDOM stands for the fraction of NPP routed to dissolved organic matter (DOM). The two types of DOM are semi-labile (DOC, DOP, DON, and DOFe) and refractory (DOCr, DOPr, and DONr). Carbon isotopes ( 12 C, 13 C, and 14 C) are calculated separately for DOC and DOCr. The run ID is 210310m for the MESMO 3 experiment LVR and 210310o for the experiment LVR with fDOM r = 0.2 %. Diazotrophs (i = 3) are not limited by NO 3 . Nutrient uptake is based on the master nutrient variable P , and all other nutrient uptake is related to by the uptake stoichiometry Q X,i , where X is N, Fe, Si, or C. For example, Q C,i = 1 [P:C] i for PFT i. Thus, Q C,i is numerically equivalent to C : P for PFT i, but we write the equations in terms of P : C for numerical stability and convenience. The Q X,i ratios represent the flexible phytoplankton uptake stoichiometry and are described more fully in Sect. 2.2.
The temperature dependence F T of Eq. (1) is given by which is analogous to the commonly used Q 10 = 2 relationship. Light limitation F I of Eq. (1) is described by a hyperbolic function: where I is the seasonally variable solar short-wave irradiance in W m −2 . Light is attenuated exponentially from the ocean surface with a 20 m depth scale. Nutrient uptake in Eq.
(1) has a dependence on z ml , which is diagnosed using the σ t density gradient criterion (Levitus, 1982). Following the Sverdrup (1953) model of the spring bloom, Eq. (1) allows for the shoaling of z ml relative to z c to enhance nutrient uptake.

Phytoplankton uptake stoichiometry
As noted above, all nutrients and O 2 are related to the main model currency P by Q X,i . We describe three different, mutually exclusive formulations in this section. The standard formulation is the power law model (Matsumoto et al., 2020;Tanioka and Matsumoto, 2017). The other two (linear model and optimality-based model of stoichiometry) are alternative formulations that have been coded, and the user can activate them (one at a time) in place of the power law formulation. However, the alternative formulations are not calibrated. Key parameter values are given in Table 2b for the power law formulation.

Power law model of stoichiometry
The uptake P : C and N : C ratios are calculated using the power law formulation as a function of ambient concentrations of phosphate [PO 4 ], nitrate [NO 3 ], temperature (T ),   = 0; thus, the environmental driver PO 4 does not drive the N : C ratio). The reference ratios are in ‰ so that [P : C] 0 = 11.6 ‰ (i.e., C : P = 86.2) for eukaryotes, and [P : C] 0 = 6.3 ‰ (i.e., C : P = 158.7) for cyanobacteria and diazotrophs. The reference ratio [N : C] 0 = 151.0 ‰ for all PFTs (i.e., C : N : = 106 : 16) is the Redfield ratio. and irradiance (I ). [ Equations (5) and (6) are the power law equations that calculate the change in P : C and N : C for fractional changes in environmental drivers relative to the reference P : C and N : C, respectively (Matsumoto et al., 2020;Tanioka and Matsumoto, 2017). The exponents are the sensitivity factors determined by a meta-analysis (Tanioka and Matsumoto, 2020a). Subscript "0" indicates the reference values (Table 2b). We have hard bounds for the calculated P : C and N : C ratios to be within 26.6 < C : P < 546.7 and 2 < C : N < 30 as observed (Martiny et al., 2013).
The P : C and N : C ratios from Eqs. (5) and (6) can then be converted to Q N,i and Q C,i for use in Eq. (2).

Linear model of stoichiometry by Galbraith and Martiny
A much simpler, alternative formulation for P : C and N : C is the model of Galbraith and Martiny (2015), where P : C is a linear function of [PO 4 ] (in µM) and N : C is a Holling type 2 functional form with a frugality behavior only at very low [NO 3 ] (in µM). The same P : C and N : C values are applied to all three PFTs.

Optimality-based model of stoichiometry
The optimality-based model of phytoplankton growth is based on the chain model, which connects the cellular P, N, and C acquisition via a chain of limitations, where the P quota limits N assimilation and the N quota drives carbon fixation Oschlies, 2009, 2013). Resource allocations of cellular P, N, and C among different cellular compartments are derived from balancing energy gain from gross carbon fixation and energy loss due to nutrient acquisition and light harvesting. The optimality-based model by  computes C : N and C : P as a function of nutrient availability (PO 4 and NO 3 ), irradiance, and day length. Temperature dependence was added by Arteaga et al. (2014) following the simple logarithmic temperature dependence on maximum nutrient uptake rate of Eppley (1972). Different versions of this optimality-based model have previously been successfully implemented in global ocean biogeochemical models, such as the Pelagic Interactions Scheme for Carbon and Ecosystem Studies (PISCES) (Kwiatkowski et al., 2018(Kwiatkowski et al., , 2019 and the University of Victoria Earth System Model (UVic) Pahlow et al., 2020). However, as we are not describing any results in this paper, we will only mention here that there is an option to calculate C : N : P using this stoichiometry model in MESMO 3. The full description of the optimality-based stoichiometry model and its parameter calibration is presented specifically for the UVic model elsewhere Pahlow et al., 2020).

Stoichiometry of iron and silica
Iron uptake stoichiometry Q Fe,i is calculated as a function of FeT following the power law formulation of Ridgwell (2001). Key parameter values are given in Table 2c.
[Fe : For all PFTs, the power law exponent s Fe:C in Eq. (12) is −0.65. The allowable Fe : C ratio is bounded at the low end by the hard-bound minimum Fe : C of 1 : 220 000. The scaling constant or [C : Fe] ref,i is set differently for PFTs, with eukaryotes having a higher base [C : Fe] ref,i than cyanobacteria and diazotrophs (115 623 : 1 and 31 805 : 1, respectively). The high end of the allowable Fe : C ratio is bounded by [C : Fe] min,i (i.e., maximum Fe : C) of 15 000 : 1 for eukaryotes and 20 000 : 1 for cyanobacteria or diazotrophs. These parameters directly follow Ridgwell (2001), who fitted power law functions to the experimental data (Sunda and Huntsman, 1995).
Silica uptake stoichiometry by eukaryotes Q Si is a power law of FeT concentration and increases with a decrease in [FeT] (Brzezinski, 2002). The power law exponent s Si:N is set to 0.7. The Si : N ratio is limited to a maximum of 18 and a minimum of 1.
O 2 liberated by phytoplankton during photosynthesis per PO 4 consumed (Q −O 2, i ) is calculated from the uptake C : P and N : P ratios (Tanioka and Matsumoto, 2020b):

Production of POM and DOM
In the top 100 m of the model domain, where phytoplankton P uptake occurs (i.e., i > 0, see Sect. 2.1), NPP is immediately routed to POM and DOM pools (Fig. 1). The production fluxes of POM, DOM sl , and DOM r from NPP are given as Jprod. Here we write the equations in terms of the master nutrient variable P: The term fDOM denotes the fraction of NPP that is routed to DOM as opposed to POM. Likewise, fDOM r is the fraction of DOM that is routed to DOM r as opposed to DOM sl . The value of fDOM r is not well known but estimated to be ∼ 1 % (Hansell, 2013), which we tentatively adopt in MESMO 3. If DOM r is not selected in the model run, fDOM r = 0. In previous versions of MESMO, fDOM was assigned a constant value of 0.67. In reality, a large variability is observed locally for this ratio, ranging from 0.01-0.2 in temperate waters to 0.1-0.7 in the Southern Ocean (Dunne et al., 2005;Henson et al., 2011;Laws et al., 2000). In MESMO 3, fDOM is calculated as a function of the ambient temperature following Laws et al. (2000): This formulation gives low export efficiency (i.e., high fDOM) in the warmer regions compared to the colder highlatitude regions. Locally, we impose fixed fDOM upper and lower bounds of 0.96 and 0.28, respectively, as estimated from a previous study (Dunne et al., 2005).
In MESMO 3, a new DOM production pathway below the production layer is available as an option. In previous MESMO versions, sinking POM was respired in the water column with the loss of O 2 directly to the dissolved inorganic forms (i.e., POC → DIC, POP → PO 4 , and PON → NO 3 ). In the new "deep POC split" pathway, sinking POM is simply broken down into DOM without the loss of O 2 as in the production layer (Fig. 1). If DOM r is selected in the model, the broken-down POM is further routed to both DOM sl and DOM r according to fDOM r . If not, all of the broken down POM is converted to DOM sl . Thus, when the deep POC split is activated, the presence of DOM in the deep ocean can be accounted for by in situ production of DOM and DOM r in addition to DOM transport from the surface. Thus, the deep POC split pathway offers an alternative means to control deep ocean DOM distribution.

Production of CaCO 3 and opal by eukaryotes
In MESMO 2, opal production was associated with the "large" PFT, and CaCO 3 production was associated with the "small" PFT. We recognize that coccolithophorids and diatoms, which are the producers of these biogenic tests, are both eukaryotes. Therefore, in MESMO 3, we associate both CaCO 3 and opal production with the POP production by the same eukaryote PFT (Jprod POP1 ): The concept of the residual nitrate potential growth (RNPG) (Balch et al., 2016) is useful in allowing competition between diatoms and non-siliceous phytoplankton within the same PFT (Matsumoto et al., 2020). Typically, in the real ocean, non-Si phytoplankton are able to grow faster and dominate the community if Si concentration is low and diatom growth is Si limited. Otherwise, diatoms are more competitive, as they have higher intrinsic growth rates. The RNPG index recasts the ambient concentrations of NO 3 and Si(OH) 4 into potential algal growth rates: If RNPG is more positive, the index indicates that nitratedependent growth exceeds silica-dependent growth. Thus, non-Si phytoplankton are more competitive, and this leads to higher CaCO 3 production. On the other hand, a more negative RNPG implies that silica limitation for diatoms is relieved, leading to enhanced diatom growth and reduced CaCO 3 production. The RNPG index is incorporated in the calculation of the rain ratio r CaCO 3 :POC presented in Eq. (20) as follows: Equation (23) indicates the base rain ratio r CaCO 3 :POC 0 (set to 0.30) is also modified by the carbonate ion saturation state by η (set to 1.28) and temperature (see Ridgwell et al., 2007, and references therein): K sp is the solubility product of CaCO 3 . The temperature dependency of CaCO 3 formation (k T ,CaCO 3 ) is similar to that of Moore et al. (2004) where warmer temperatures favor the growth of carbonate-bearing phytoplankton. In the new model, DOM r can be activated. DOM r is produced from POM breakdown, which can occur in the production layer or throughout the water column in the "deep POC split". Possible DOM r remineralization mechanisms are the slow background degradation that occurs everywhere, thermal degradation in hydrothermal vents, and photodegradation at the surface. See the text for details.

Remineralization of POM and DOM
Once produced, both POM and DOM undergo remineralization throughout the water column. Key remineralization parameter values are given in Table 2d. Previously, POM remineralization had a temperature dependence and decayed exponentially with depth (Yamanaka et al., 2004). In MESMO 3, we incorporate an additional dependency on dissolved oxygen following Laufkötter et al. (2017): V POM is the base remineralization rate, k R expresses the temperature sensitivity of remineralization, and K O 2 is halfsaturation constant for oxygen-dependent remineralization. When the sediment model is not coupled, any POM that reaches the seafloor dissolves completely to its inorganic form and is returned to the overlying water.
In MESMO 3, all forms of semi-labile DOM remineralize at the same rate. It is represented by τ sl , the inverse of the timescale of DOM sl decay, which has been estimated previously to be ∼ 1.5 years (Hansell, 2013): All forms of DOM r also remineralize at the same rate in MESMO 3. In total, there are three optional, additive sinks of DOM r in the model: slow background decay, photodegradation, and degradation via hydrothermal vents (Fig. 1). Observations clearly indicate that the 14 C age of deep-ocean DOC r is 10 3 years (e.g., Druffel et al., 1992), much older than DI 14 C. In addition, the deep ocean DOC r concentration decreases modestly along the path of the deep water from the deep North Atlantic to the deep North Pacific (Hansell and Carlson, 1998). Thus, it is understood that there is a slow DOM r background decay in the deep ocean. We represent this ubiquitous process with τ bg , which is the inverse of the background decay timescale, estimated to be ∼ 16 000 years (Hansell, 2013).
Observations to date indicate that photodegradation is a major sink of DOM r (e.g., Mopper et al., 1991). This process is believed to convert DOM r that is upwelled from the ocean interior into the euphotic zone into more labile forms of DOM. We represent photodegradation with τ photo , the inverse of the decay timescale, estimated to be ∼ 70 years (Yamanaka and Tajika, 1997). This occurs only in the surface.
Finally, observations of DOM emanating from different types of hydrothermal vents indicate that they have variable impacts on the deep-sea DOM r (Lang et al., 2006). However, the off-axis vents circulate far more seawater through the fractured oceanic crust than the high-temperature and diffuse vents and are thus believed to determine the overall impact of the vents on the deep-sea DOM r as a net sink (Lang et al., 2006). Here we assume simply that seawater that circulates through the vents loses all DOM r (i.e., 1/τ vent < t, where t is the biogeochemical model time step of 0.05 year). This means that the more seawater circulates through the vents, the more DOM r is removed: the total removal rate depends on the vent flux of seawater H flux . We implement the vent degradation of DOM r in MESMO 3 by first identifying the wet grid boxes located immediately above known midocean ridges. We then distribute the annual global H flux of 4.8 × 10 16 kg yr −1 (Lang et al., 2006) equally among those ridge-associated grid boxes. The grid cells contain a mass of seawater much greater than the mass that circulates through vents in t (10 21 kg vs. 10 13 kg). Therefore, the seawater mass in the vent grid cells that does not circulate through the vents in t is subject only to background degradation in MESMO 3.
The three DOM r sinks are not mutually exclusive. They can thus be combined to yield the total DOM r remineralization rate: where SW flux_local is the mass of seawater that circulates through the vents in each grid box in t, and SW grid is the total mass of seawater in the same grid box. The amount of O 2 respired as a result of these POM and DOM remineralization processes is related to the organic carbon pools by the respiratory quotients of POC and DOC, r −O 2 :POC and r −O 2 :DOC , respectively. These are molar ratios of O 2 consumed per unit organic carbon respired. They are variable and calculated from the ambient POM and DOM concentration (Tanioka and Matsumoto, 2020b):

Remineralization of CaCO 3 and opal
Remineralization of CaCO 3 and opal particles occurs as they sink through the water column and remains the same as in MESMO 2. Key parameter values are given in Table 2d. Remineralization of CaCO 3 is a function of temperature similar to that of particulate organic matter remineralization but without oxygen dependency. The temperature-dependence term k R modifies the base remineralization rate V CaCO 3 : Opal remineralization in the water column follows Ridgwell et al. (2002). The rate of opal remineralization R opal is given by the product of normalized dissolution rate (r opal ), base opal dissolution rate (k opal ), and opal concentration [opal].
Si(OH) 4 eq (34) r opal is a function of temperature (T ) and the degree of undersaturation (u opal ), which in turn is calculated from the ambient Si(OH) 4 and Si(OH) 4 at equilibrium. The equilibrium concentration is a function of ambient temperature: log 10 Si(OH) 4 eq = 6.44 − 968 T (K) .
Without the sediment module of MESMO activated, both CaCO 3 and opal particles that reach the seafloor are completely dissolved back to inorganic forms.

Conservation of organic matter and biogenic tests
The time rate of change of the biogenic organic matter and tests are given by the sum of the production terms (i.e., sources) and the remineralization terms (i.e., sinks). The circulation-related transport terms are omitted as noted above, but the vertical transport due to particle sinking is included here. The sinking speed w is the same for all particles. The sum of POM i of all the PFTs give the total POM concentrations.
The time rate of change of CaCO 3 and opal is expressed in much the same way as POM.
The DOM pools have the production and remineralization terms without the particle sinking term.

Conservation of inorganic nutrients
The time rate of change of the inorganic nutrients have organic carbon production as sink terms and remineralization as source terms. The production terms (J prod ) are zero below the upper-ocean production layer. Nutrients have a unit In Eq. (51), Fix N is the N fixation carried out by diazotrophs, and Den N is the water column denitrification. There is an air-sea gas exchange term F gas in Eqs. (52) and (56) for gaseous CO 2 and O 2 , respectively. In Eq. (53), alkalinity increases with decreasing nitrate concentrations and increasing CaCO 3 dissolution. Equation (54) contains R POMFe , which is an iron source that represents remineralization of the Fe scavenged by sinking particles. These terms are explained in the following sections.

Prognostic nitrogen cycle
Biological production by diazotrophs is stimulated when the ambient NO 3 is low. Nitrogen fixed by diazotrophs during their growth is added to the marine NO 3 pool. The prognostic nitrogen fixation model employed here is similar to that used in the HAMOCC biogeochemical module (Paulsen et al., 2017): where Fix N is the nitrogen fixation rate and I NO 3 is the nitrate dependency term in quadratic Michaelis-Menten kinetics form with the half-saturation constant K NO 3 _Nfix . See Table 2e for the values related to the N cycle. Water column denitrification is formulated in an approach similar to that of the original GENIE model (Ridgwell et al., 2007), in which 2 mol of NO 3 are converted to 1 mol of N 2 and liberating 2.5 mol of O 2 as a byproduct: Denitrification takes place in grid boxes, in which O 2 concentration is below a threshold concentration (O 2,def ) and is stimulated if the total global inventory of NO 3 relative to PO 4 is high. In other words, denitrification can effectively act as negative feedback to nitrogen fixation. The threshold O 2 concentration (O 2,def ) takes the minimum of the hard-bound O 2 threshold concentration (O 2,crit ) and the NO 3 to PO 4 ratio, scaled by a parameter k D . The parameters O 2,crit and k D are calibrated to give the global denitrification rate of roughly 100 Tg N yr −1 , which balances the total nitrogen fixation rate in the model.

Prognostic iron cycle
The iron cycle in MESMO 3 remains the same as in MESMO 2. Key parameter values are given in Table 2e. The two species of dissolved iron (Fe and FeL) are partitioned according to the following equilibrium relationship: where [L] is the ligand concentration and K ligand is the conditional stability constant. The sum of ligand and FeL is set at a constant value of 1 nmol kg −1 everywhere. Iron is introduced into the model domain by a constant fraction (3.5 wt %) of aeolian dust deposition at the surface (F in ) following the prescribed modern flux pattern (Mahowald et al., 2006) with constant solubility (β): Particle-scavenged iron POM Fe (note the difference from POFe) is produced below the productive layer when sinking POM scavenges Fe to sinking POM: where τ sc and K o are empirical parameters that determine the strength of scavenging. Remineralization of Fe scavenged to POM (POM Fe ) is identical in form to that of POM remineralization: The conservation equation of the particle scavenged iron is thus expressed as follows: Any scavenged iron that escapes remineralization in the water column reaching the seafloor is removed from the model domain in order to keep the total Fe inventory at a steady state.

Air-sea gas exchange
The air-sea gas exchange formulation remains the same as in MESMO 2 and follows Ridgwell et al. (2007). It is the function of gas transfer velocity, the ambient dissolved gas concentration, and saturation gas concentration. The flux of CO 2 and O 2 gases across the air-sea interface is given by where k is the gas transfer velocity, [CO 2 ] sat and [O 2 ] sat are saturation concentrations, and A is the fractional ice-covered area that is calculated by the physical model. Gas transfer velocity k is a function of wind speed (u) following Wanninkhof (1992) where Sc is the Schmidt number for a specific gas:

Results and discussion
All new results from MESMO 3 presented here are from the steady state. On a single computer core at the Minnesota Supercomputing Institute, it takes approximately 1 h to complete 1000 years of MESMO 3 simulation. The "standard" MESMO 3 has the power law model of flexible stoichiometry but no DOM r . The results from the standard model (hereafter just MESMO 3) are presented in Sect. 3.1, and the results from the DOM r -enabled model are presented in Sect. 3.2. In Table 3, we summarize and compare key biogeochemical diagnostics from MESMO 3 against those from MESMO 2 and available observational constraints. The global NPP, as well as global export production of POC and opal, are comparable or somewhat lower in MESMO 3 than MESMO 2. We relied on experience to calibrate MESMO 3 with the primary goal of reasonably simulating the phytoplankton community composition and C : N : P ratio (e.g., abundant cyanobacteria and high ratio for all PFTs in oligotrophic gyres). We tried to improve or at least preserve the gains that we achieved in earlier versions of MESMO in terms of the global distributions of PO 4 , NO 3 , O 2 , Si(OH) 4 , and FeT (Supplement Figs. S1, S2, S3, S4, and S5). Overall, there is a stronger nutrient depletion in the new model. For example, the surface PO 4 and NO 3 in MESMO 3 are now sufficiently depleted in the subtropical gyres but are too low in the eastern equatorial Pacific when compared to the World Ocean Atlas (Fig. S1; see RMSE in Table 3). It is a challenge for MESMO and other coarse-resolution models to simulate narrow dynamical features such as equatorial upwelling and reproduce biogeochemical features with sharp gradients. The spatial pattern of POC export that drives this surface nutrient pattern is similar in the two models (Fig. S2). In the 1D global profile, there is a marked improvement in the subsurface distribution of O 2 in MESMO 3 over MESMO 2. Whereas the depth of the oxygen minimum was ∼ 300 m in MESMO 2, it is ∼ 1000 m in both MESMO 3 and the World Ocean Atlas (Fig. S3). At 1000 m, the O 2 minimum is located in the far North Pacific in MESMO 3, whereas in the World Ocean Atlas it occurs in both the Northeast Pacific and the Arabian Sea. In contrast, the world ocean at 1000 m is too well oxygenated in MESMO 2. We believe that the improved match in the O 2 minimum depth would help simulate denitrification in the correct depth range, and there is a modest improvement in the data-model O 2 mismatch in terms of RMSE in MESMO 3 over MESMO 2 ( Table 3). The deepening of the O 2 minimum was achieved largely by increasing the particle sinking speed, which tends to strengthen the biological pump and deplete the surface nutrients. This also helps MESMO 3 preserve MESMO 2's surface Si(OH) 4 depletion in much of the world ocean except in the North Pacific and Southern Ocean (Fig. S4). This is a feature captured by Si * < 0 (Si * =[Si(OH) 4 ]-[NO 3 ]) in observations (Sarmiento et al., 2004) and simulated previously by MESMO 2 and now MESMO 3. Finally, surface FeT is also depleted more strongly in MESMO 3 over MESMO 2, except the North Atlantic, where aeolian deposition of dust from the Sahara maintains a steady Fe supply (Fig. S5).
In MESMO 3, we made no effort to calibrate all the semi-labile pools of DOM: DOC sl , DOP sl , DON sl , and DOFe sl . We note only that the surface DOC sl concentration of 58 µmol kg −1 and DOC export production of 1.4 Pg C yr −1 in MESMO 3 are higher than in MESMO 2 (24 µmol kg −1 and 0.4 Pg C yr −1 , respectively). The higher surface concentration is due to the longer τ sl in MESMO 3 (Table 2d). The global average of the temperature-dependent fDOM in MESMO 3 is 0.69, which is slightly higher than the spatially uniform value of 0.67 in MESMO 2.

Novel features of MESMO 3
An important new feature of MESMO 3 is the representation of the primary producers by three PFTs (Fig. 2). The eukaryotes are characterized by the highest maximum growth rate and high half-saturation constants. Thus, the eukaryotes are more dominant than the other PFTs in the more eutrophic waters of the equatorial and polar regions (Fig. 2a). The cyanobacteria have smaller half-saturation constants and a NPP for MESMO 2 was unavailable as a model output and is therefore estimated from POC and fDOM = 0.67. b NPP (in terms of C) is needed in the calculation of the PFT abundance. The root-mean-square error (RMSE) of the simulated P, N, Si, and O 2 distributions from MESMO 2 and 3 was calculated relative to the World Ocean Atlas 2018 (WOA18) gridded data . The model-data comparison is made in the top 100 m for nutrients and below 100 m for O 2 . WOA18 was regridded to the MESMO 3 grid to calculate the RMSE. References for independent constraints are as follows: (1) global NPP (Carr et al., 2006), (2) global POC export (DeVries and Weber, 2017), (3) global DOC export assumed to be 20 % of total carbon export (Hansell et al., 2009;Roshan and DeVries, 2017), (4) global opal , (5) global CaCO 3 export (Berelson et al., 2007), (6) global N fixation and denitrification rates (Landolfi et al., 2018), (7) uptake C : N : P ratio is based on POM measurements (Martiny et al., 2013), (8) export C : N : P ratio is assumed to equal the subsurface remineralization ratio (Anderson and Sarmiento, 1994), and (9) deep O 2 from WOA18 below 100 m .
thus are more dominant in the oligotrophic subtropical gyres (Fig. 2c). The diazotrophs do not have NO 3 limitation but have the lowest maximum growth rate. Thus it is much lower in abundance than the other two PFTs generally, and outcompeted in transient blooms and thus excluded in higher latitudes (Fig. 2e).
Figure 2 also indicates that all three PFTs show Fe limitation in the Southern Ocean. Outside the Southern Ocean, the eukaryotes are primarily limited by Si(OH) 4 (Fig. 2b). As far as organic carbon is concerned, we consider the eukaryotes to basically represent diatoms, which are arguably the most important agent of organic carbon export. In this context, the widespread silica limitation for eukaryotes would be consis- tent with the notions that silica uptake by diatoms should be limited in ∼ 60 % of the world surface ocean (Ragueneau et al., 2000) and that much the world ocean thermocline is filled with silica-depleted water (Si * < 0 as noted above). On the other hand, the cyanobacteria are largely limited by NO 3 outside the Southern Ocean (Fig. 2d). The diazotrophs are limited by iron in much of the world ocean except in the Atlantic basin (Fig. 2f), where surface PO 4 is strongly depleted in both observations (Mather et al., 2008) and in our model (Fig. S1). Figure 3 illustrates the influence of the RNPG index, which was implemented in MESMO 3 to allow for the effect of competition between diatoms and coccolithophores within the same PFT (Eqs. 22 and 23). The eukaryote NPP (Fig. 3a) is effectively split into two parts: one is associated with diatoms and opal production (Fig. 3b), and the other is associated with coccolithophores and CaCO 3 production (Fig. 3c). According to the RNPG index, opal production is simulated more in the higher latitudes of the Southern Ocean and the North Pacific, where surface [Si(OH) 4 ] is abundant. Elsewhere, CaCO 3 production is relatively large. The decoupling is prominent in the North Indian Ocean. Note that the spatial pattern of CaCO 3 production is quite different in MESMO 3 (Fig. 3c) and MESMO 2 (Fig. 3d) because CaCO 3 production was associated in MESMO 2 with the "small" PFT, which corresponds to the cyanobacteria PFT in MESMO 3.
The global pattern of the mean C : P uptake ratio in the production layer is shown in Fig. 4. Consistent with ob- Figure 3. Eukaryote production in MESMO 3 and CaCO 3 export in MESMO 2. In MESMO 3, eukaryote NPP (a) is linked to both opal export (b) and CaCO 3 export (c) but the two export productions are differentiated by the residual nitrate potential growth (RNPG). Compare CaCO 3 export in MESMO 3 (c) to MESMO 2 (d) (unit: mol m −2 yr −1 ).
servations (Martiny et al., 2013), the simulated C : P ratio of the phytoplankton community is elevated in the oligotrophic subtropical gyres and low in the eutrophic polar waters (Fig. 4a). The community C : P ratio exceeds 200 in the gyres and reaches as low as 40-50 in the Southern Ocean. The community C : P has contributions from both physiological effects (i.e., environmental drivers acting on each PFT's C : P ratio) and taxonomic effects (i.e., the shift in the community composition changes the weighting of each PFT's C : P ratio). Figure 4b shows that the community C : P is high in oligotrophic gyres because cyanobacteria and to a lesser extent diazotrophs dominate the community and their C : P ratio is high. Conversely, the community C : P is low in the polar waters because the eukaryotes dominate and their C : P ratio is low. For both eukaryotes and cyanobacteria, their C : P is high in oligotrophic subtropical gyres because PO 4 is low (Fig. 4c and d). This physiological effect is larger in eukaryotes than cyanobacteria because the former has greater sensitivity (i.e., larger sensitivity factor s P:C PO 4 ; see Eq. 5 and Table 2b). However, the cyanobacteria PFT's C : P ratio has an additional sensitivity to temperature (i.e., s P:C T = 0) that elevates their C : P in the lower latitudes. We do not show the C : P ratio for diazotrophs because it is very similar to that of cyanobacteria (Fig. 4b, d).
In order to gain more insights into the spatial patterns of the C : P ratio (Fig. 4), we examined the relationships between the C : P and C : N ratios and the four possible environmental drivers for eukaryotes and cyanobacteria ( Fig. 5; again, diazotrophs are not shown). The red plots show that there is a causal relationship between the ratios and the drivers as formulated by the power law model (Eqs. 5 and 6). The black plots show the absence of a causal relationship. For example, the C : P ratio of both eukaryotes and cyanobacteria is strongly correlated with PO 4 because there is a causal relationship (Fig. 5a, b, shown in red). Similarly, the C : N ratio of the same two PFTs have a strong correlation with PO 4 (Fig. 5c, d in black), but there is actually not a causal relationship (i.e., s N:C PO 4 = 0, Table 2b). The C : N-PO 4 correlation exists simply because the nutrients are well correlated. Similarly, because temperature and photosynthetically active radiation (PAR) tend to be correlated via latitude, the stoichiometry has a similar correlation to the two drivers. For example, cyanobacteria C : P has a strong correlation with both temperature and PAR (Figs. 5j, 4n), but only the temperature is a real driver. Figure 5 indicates which are the dominant drivers of the C : N : P ratio in MESMO 3. For the eukaryote C : P ratio, it is PO 4 . For the cyanobacteria C : P ratio, the important drivers are temperature and PO 4 . For the C : N ratio for , and diazotroph C : P (blue). In addition, (b) shows the mean range of observed C : P ratio binned by latitude (Martiny et al., 2013). both eukaryotes and cyanobacteria, NO 3 is more important than PAR. Figure 5 also serves to remind us that correlation does not indicate causation. Figure 6 shows the community C : P and C : N ratios plotted against the four environmental drivers. Unlike Fig. 5, which reflected the individual PFT's physiological response, Fig. 6 includes the effect of taxonomy as well. Still, the effects of PO 4 and temperature are clearly visible on the community C : P ratio. Both low [PO 4 ] and warmer waters are found in the lower latitudes, so the P frugality and temperature effects are additive. The effect of NO 3 on the community C : N ratio is also very clear, but the effect of PAR is not as clear. Thus, the overall physiological effects seen in the PFTspecific C : N : P are obvious in the community C : N : P ratio.

DOM r -enabled MESMO 3
In MESMO 2, DOC sl was a standard state variable. In MESMO 3, other forms of DOM are available as options. They are the semi-labile forms of DOM, i.e., DOP sl , DON sl , and DOFe sl , and the refractory forms of DOM, i.e., DOC r , DOP r , and DON r . MESMO 3 is not yet calibrated with re-spect to all the DOM variables, but here we demonstrate their potential use in future biogeochemical investigations by presenting steady state DOM results from the model experiment LV (experiment ID: 210310m). In this run, all three sinks of DOM r are activated: slow background decay, photodegradation, and degradation in hydrothermal vents.
The experiment name LV stands for "literature values". In LV, we use the literature values for the key DOM remineralization model parameters (Table 2d) and fDOM r = 0.01 (Hansell, 2013). All other model parameter values in the LV run are identical to the standard MESMO 3 model ( Table 2). The black lines in Fig. 7 show the global mean vertical profiles of the total DOC (DOC t = DOC sl + DOC r ) with solid lines and DOC r with a dashed line. Qualitatively, the simulated profiles are consistent with the observations, showing a near-uniform DOC r concentration and a DOC sl profile that rapidly attenuates with depth in the top few hundred meters (Hansell, 2013). However, the simulated values reach 130 µmol kg −1 at the surface, which is approximately twice the observations. More typically, the observed DOC r is 30-40 µmol kg −1 , and the observed DOC sl attenuates with depth from 30-40 µmol kg −1 near the surface. So their sum, which is represented by DOC t , is approximately 60-80 µmol kg −1 at the surface in observations. Figure 8 adds a lateral perspective to Fig. 7. The rapid DOC t attenuation in the vertical is strong in the lower latitudes where stratification is generally stronger. The transport of DOC sl from the surface to deeper waters is evident in the high latitudes of the North Atlantic and the Southern Ocean. The DOC t change in the deep ocean is limited.
Observations of deep-ocean DOC t indicate a reduction by 29 % or 14 µmol kg −1 from the deep North Atlantic to the deep North Pacific (Hansell and Carlson, 1998). Figure 8 shows that the deep-ocean DOC t gradient in LV is approximately 10 µmol kg −1 from 70-75 µmol kg −1 in the North Atlantic to < 65 µmol kg −1 in the North Pacific.
The horizontal DOC t distributions from the LV run can also be compared to a global extrapolation based on an artificial neural network (ANN) of the available DOC t data Figure 6. Scatterplots of surface ocean community C : P and C : N vs. environmental drivers in MESMO 3. (Roshan and DeVries, 2017). At the surface, the extrapolation indicates higher DOC t concentrations in the subtropical gyres (Fig. 9a), while our simulation does not clearly delineate the gyres (Fig. 9c). In our model, fDOM is temperaturedependent and strongly controls the production of DOM. The surface DOC t is thus more elevated in the lower latitudes. Interestingly, the ANN study diagnosed higher rates of DOM production in the subtropical gyres. Since the oligotrophic subtropical gyres have low NPP, the diagnosis would thus suggest that somehow fDOM is higher in the gyres. At depth, both the extrapolated and simulated DOC t show a gradual decline in concentrations from the North Atlantic to the North Pacific (Fig. 9b, d). The highest deep DOC t in the LV run is seen just south of Greenland, where convection occurs in the model.
Finally, we show that the deep-ocean radiocarbon aging is larger in DIC than in DOC t in the model (Fig. 10). The North Pacific-North Atlantic 14 C gradient is roughly −100 ‰ for DIC and −70 ‰ for DOC t . The oldest DOC t 14 C is approximately −430 ‰ in the North Pacific. If 14 C decay were the only mechanism of change along the path of the deepwater circulation, the 14 C gradient should be quite similar between DIC and DOC t , which are both dissolved phases and transported passively by the same circulation. The one potentially important difference is that the addition of the relatively young DI 14 C and DO 14 C to the deep ocean by the "deep POC split" (see Sect. 2.3). The addition impacts DOC t 14 C more than DIC 14 C because DOC t is 2 orders of magnitude lower in concentration than DIC.
In observations, the aging of DIC and DOC t is reportedly similar in the Antarctic Bottom Water (below 4000 m) of the deep Pacific (Druffel et al., 2019). This may be explained by the fact that there would not be much deep POC split occurring so deep in the ocean. The North Pacific-North Atlantic 14 C gradient, accounting for thermonuclear bomb 14 C, may be as large as −100 ‰ for DOC t (about −550 ‰ in the deep Pacific and −456 ‰ in the deep Atlantic) (Druffel et al., 2019). This gradient is not rigorously determined because there is not enough data to do an objective analysis. Therefore, the equivalent 14 C gradient for DIC cannot Figure 7. Global mean vertical profiles of DOC from the DOM Renabled MESMO 3. DOC t (DOC sl +DOC r , black line) and DOC r (dashed black line) from the LV run. The red line is DOC t after reducing fDOM r from 1 % in LV to 0.2 % (Experiment 210310o) (unit: µmol kg −1 ). be determined. However, the DIC 14 C endmember values by inspection (about −250 ‰ in the deep Pacific and −70 ‰ in the deep Atlantic) (Matsumoto and Key, 2004) indicate a clearly larger 14 C gradient for DIC than DOC t as simulated by the experiment LV.
One lesson from the data-LV run mismatch in the overall DOC t concentration (Fig. 7) and surface DOC t pattern ( Fig. 9) is that the parameter values from the literature do not fully capture the DOC cycle and/or that MESMO 3 is still lacking some important DOC process. In LV, the surface DOC t is too high because DOC r is too high, while DOC sl is not unreasonable (Fig. 7). DOC r is too high because there is too much DOC r production (e.g., fDOM r = 1 % is too large), there is too little DOC r degradation (e.g., one of the DOM decay mechanisms is too slow; Eq. 28 and Table 2d), or some combination of both. For example, fDOM r is a key parameter that is not well constrained by observations. Had we used 0.2 % instead of 1 % for fDOM r , the global mean surface DOC t drops to 76 µmol kg −1 (red line, Fig. 7), consistent with observations. For achieving a better surface DOC t pattern, we may need a different formulation of fDOM that is, for example, negatively related to nutrient concentrations so that fDOM increases in the oligotrophic subtropical gyres (Roshan and DeVries, 2017).
Another lesson from the DOM modeling exercise is that it is important to simulate DOP r reasonably well in order to preserve the favorable results we achieved in MESMO 3 with respect to biological production and the phytoplankton C : N : P ratio. We find that in the experiment LV, the global mean DOP r concentration becomes steady at 0.45 µmol-P kg −1 . In observations, the mean DOC r is about 40 µmol-C kg −1 and the DOC r : DOP r ratio is estimated to be ∼ 1370 : 1 (Letscher and Moore, 2015), so DOP r concentra-tion should only be roughly 0.03 µmol-P kg −1 . Thus, the simulated DOP r = 0.45 µmol-P kg −1 is an order of magnitude too high. Because there is more P in the form of DOP r in LV, the oceanic inventory of PO 4 declines, causing a nearly 10 % drop in export production compared to the standard MESMO 3. In LV, the decline in the surface ocean PO 4 that accompanies the change in the PO 4 inventory acts on the phytoplankton physiology (i.e., P effect on C : P in Eq. 5), which leads to a large rise in the global mean phytoplankton community C : P export ratio from 113 : 1 to 127 : 1. The implementation of preferential remineralization of DOP (and DON) over DOC (Letscher and Moore, 2015) is one way to deal with the problem of too high DOP r concentrations.

Large-scale patterns of N 2 fixation and denitrification
The modeled habitat of diazotrophs is concentrated in tropical and subtropical waters between 40 • S and 40 • N and limited by iron (Fig. 1e, f). Most noticeably in the North Pacific subtropical gyre, diazotrophs constitute ∼ 40 % of total NPP. The latitudinal extent of diazotrophs is mainly determined by surface nitrate availability and physical factors such as surface temperature and irradiance. Low nitrate availability in subtropical gyres gives diazotrophs a competitive advantage over small cyanobacteria. Warm temperature and high irradiance are also critical physical factors that drive the growth of diazotrophs in the model. The modeled global depth-integrated N 2 fixation is 101 Tg N yr −1 (Table 3), and this value falls well within the range of observational and geochemical constraints of 80-200 Tg N yr −1 (Landolfi et al., 2018). In MESMO 3, N 2 fixation occurs in the North Pacific and mid-to-low latitudes of the Atlantic basin (Supplement Fig. S6), where diazotrophs are generally more abundant (Fig. 2e). The elevated N 2 fixation rate in the North Pacific, where nitrate limits eukaryotes and cyanobacteria (Fig. 2b, d), can be explained by the healthy growth of diazotrophs, which is not limited by N. In the subtropical and tropical Atlantic and the Indian Ocean, high N 2 fixation is driven by a elevated C : P and N : P ratio (Fig. 4), exemplified by low phosphate availability and warm surface temperature. This spatial pattern agrees with a recent inverse model study , which showed an elevated N 2 fixation rate in subtropical gyres.
Global water column denitrification is 101 Tg N yr −1 (Table 3) and is equal to the global N 2 fixation because the model has reached steady state. Denitrification is restricted to the subpolar North Pacific (Fig. S6), where sub-surface oxygen concentration is significantly depleted (Fig. S3d). Enhanced denitrification in this region is in qualitative agreement with a previous modeling study (Bianchi et al., 2018), which showed the anaerobic niche due to particle microenvironments can significantly expand the hypoxic expanses in the North Pacific. However, the extent of denitrification in our model does not include the eastern equatorial Pacific   and northern Indian oceans, which are important hotspots for denitrification (Codispoti, 2007;Deutsch et al., 2007). This issue is typical of coarse-resolution global ocean biogeochemistry models that lack spatial resolution in reproducing intense upwelling (Marchal et al., 1998;Najjar et al., 1992;Yamanaka and Tajika, 1997).
Finally, the ratio of the global inventories of NO 3 and PO 4 in MESMO 3 is just about 16 at steady state, consistent with observations (Gruber and Sarmiento, 1997). One key model parameter in this regard is the nitrate uptake half-saturation constant of diazotrophs, K NO 3 _Nfix in Eq. (58). A large value of K NO 3 _Nfix will make it hard for diazotrophs to obtain fixed N from NO 3 , which would facilitate N 2 fixation and pushes up the global N/P ratio. With a smaller value of K NO 3 _Nfix , diazotrophs will more easily uptake NO 3 , thus depressing N 2 fixation and lowering the global N/P ratio.

Conclusions
MESMO 3, the third and latest version of MESMO, is comprehensively described here. With a fully flexible C : N : P ratio in three PFTs, a prognostic N cycle, and more mechanistic schemes of organic matter production and remineralization, MESMO 3 reflects the evolving and accumulating knowledge of the ocean biogeochemistry. The model thus remains an effective tool for investigations of the global biogeochemical cycles, especially over long timescales, given the model's computational efficiency. In particular, MESMO 3 holds promise for studying the marine DOM cycle. The optional features of MESMO 3 include the semi-labile and refractory pools of C, P, N, and Fe. The fact that the literature values regarding the present marine DOM cycle are unable to simulate key observations indicates an opportunity for MESMO 3 to contribute to an improved understanding of the marine DOM cycle.
Code and data availability. Model results presented in this study are archived and available with the code. The complete code of MESMO version 3.0 and the results presented here are available at GitHub https://github.com/gaia3intc/mesmo.git (last access: 26 April 2021) and have the following DOI: https://doi.org/10.5281/zenodo.4403605 (Matsumoto, 2020).
Author contributions. KM, TT, and JZ developed the model code. KM performed the simulations, carried out the analyses, and archived the model code and results. KM and TT wrote the paper.