Articles | Volume 14, issue 4
Geosci. Model Dev., 14, 2265–2288, 2021
Geosci. Model Dev., 14, 2265–2288, 2021

Model description paper 30 Apr 2021

Model description paper | 30 Apr 2021

MESMO 3: Flexible phytoplankton stoichiometry and refractory dissolved organic matter

MESMO 3: Flexible phytoplankton stoichiometry and refractory dissolved organic matter
Katsumi Matsumoto1, Tatsuro Tanioka1,a, and Jacob Zahn1 Katsumi Matsumoto et al.
  • 1Department of Earth & Environmental Sciences, University of Minnesota, Minneapolis, Minnesota, USA
  • acurrent affiliation: Department of Earth System Science, University of California Irvine, Irvine, California, USA

Correspondence: Katsumi Matsumoto (


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 temperature-dependent 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.

1 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, 2013), 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 (14C) 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 (PO4), nitrate (NO3), and aqueous CO2. 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.

Table 1Summary of MESMO development.

PFT stands for phytoplankton functional types. MESMO 2 PFTs are as follows: 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 (12C, 13C, and 14C) 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 fDOMr=0.2 %.

Download Print Version | Download XLSX

The main aim of MESMO 2 was a credible representation of the marine silica cycle (Matsumoto et al., 2013). 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 (30Si) 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 half-saturation 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. CaCO3 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 multi-objective tuning and extensive model–data comparisons of transient tracers (anthropogenic carbon, CFCs), deep ocean Δ14C, and nutrients (Matsumoto et al., 2008, 2013). 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., 2010, 2020; Matsumoto 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 (Archer et al., 2009; Cao et al., 2009; Eby et al., 2013; Joos 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 column denitrification, the dependence of organic matter remineralization on the dissolved oxygen (O2) and temperature, and CaCO3 production based on the concept of the residual nitrate potential growth. In addition, we describe the semi-labile DOM for P, N, and Fe (DOPsl, DONsl, and DOFesl) and the refractory DOM for C, P, and N (DOCr, DOPr, and DONr), 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; Tanioka and Matsumoto, 2017, 2020a). This work consolidates the descriptions of all these features in a single publication.

2 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.

Table 2MESMO 3 biogeochemical model parameters values: (a) phytoplankton nutrient uptake; (b) power law model of flexible C:N:P stoichiometry; (c) iron uptake stoichiometry; (d) parameters related to POM, DOM, CaCO3, and opal; (e) nitrogen and iron cycles.

* Sensitivity factors not listed in Table 2b have a value of zero (e.g., sPO4N:C=0; thus, the environmental driver PO4 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.

Download XLSX

2.1 Phytoplankton nutrient uptake

NPP occurs in the top two vertical levels of the ocean domain above the fixed compensation depth (zc) 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 (zml):

(1) Γ i = 1 τ i F N , i F T F I max 1 , z c z ml .

Subscript i refers to PFT (i=1: eukaryotes; i=2: cyanobacteria; i=3: diazotrophs). The nutrient dependence FN,i is given by Liebig's law of minimum combined with Michaelis–Menten uptake kinetics of limiting nutrients: PO4, NO3, CO2 (aq), total dissolved iron (sum of free iron and ligand-bound iron: FeT= Fe+FeL), and Si(OH)4:

(2) F N , i = min PO 4 PO 4 + K PO 4 , i PO 4 , NO 3 NO 3 + K NO 3 , i NO 3 Q N , i - 1 , CO 2 aq CO 2 aq + K CO 2 , i CO 2 aq Q C , i - 1 , FeT FeT + K FeT , i FeT Q Fe , i - 1 , Si OH 4 Si OH 4 + K Si OH 4 Si OH 4 Q Si - 1 ,

where KX,i is the half-saturation concentration of nutrient X for PFT i. Only eukaryotes (i=1) are limited by Si(OH)4. Diazotrophs (i=3) are not limited by NO3. Nutrient uptake Γ is based on the master nutrient variable P, and all other nutrient uptake is related to Γ by the uptake stoichiometry QX,i, where X is N, Fe, Si, or C. For example, QC,i=1P:Ci for PFT i. Thus, QC,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 QX,i ratios represent the flexible phytoplankton uptake stoichiometry and are described more fully in Sect. 2.2.

The temperature dependence FT of Eq. (1) is given by

(3) F T = T C + 2 T C + 10 ,

which is analogous to the commonly used Q10=2 relationship. Light limitation FI of Eq. (1) is described by a hyperbolic function:

(4) F T = I I + 20 ,

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 zml, 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 zml relative to zc to enhance nutrient uptake.

2.2 Phytoplankton uptake stoichiometry

As noted above, all nutrients and O2 are related to the main model currency P by QX,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.

2.2.1 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 [PO4], nitrate [NO3], temperature (T), 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 QN,i and QC,i for use in Eq. (2).


2.2.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 [PO4] (in µM) and N:C is a Holling type 2 functional form with a frugality behavior only at very low [NO3] (in µM). The same P:C and N:C values are applied to all three PFTs.


2.2.3 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 (Pahlow et al., 2013; Pahlow and 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 Pahlow et al. (2013) computes C:N and C:P as a function of nutrient availability (PO4 and NO3), 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, 2019) and the University of Victoria Earth System Model (UVic) (Chien et al., 2020; 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 (Chien et al., 2020; Pahlow et al., 2020).

2.2.4 Stoichiometry of iron and silica

Iron uptake stoichiometry QFe,i is calculated as a function of FeT following the power law formulation of Ridgwell (2001). Key parameter values are given in Table 2c.


For all PFTs, the power law exponent sFe: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 QSi is a power law of FeT concentration and increases with a decrease in [FeT] (Brzezinski, 2002). The power law exponent sSi:N is set to 0.7. The Si:N ratio is limited to a maximum of 18 and a minimum of 1.


O2 liberated by phytoplankton during photosynthesis per PO4 consumed (Q-O2,i) is calculated from the uptake C:P and N:P ratios (Tanioka and Matsumoto, 2020b):

(15) Q - O 2 , i = 1.1 Q C , i + 2 Q N , i .

2.3 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, DOMsl, and DOMr from NPP are given as Jprod. Here we write the equations in terms of the master nutrient variable P:


Figure 1Schematic diagram of DOM cycling in MESMO 2 vs. MESMO 3. In the new model, DOMr can be activated. DOMr is produced from POM breakdown, which can occur in the production layer or throughout the water column in the “deep POC split”. Possible DOMr 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.


The term fDOM denotes the fraction of NPP that is routed to DOM as opposed to POM. Likewise, fDOMr is the fraction of DOM that is routed to DOMr as opposed to DOMsl. The value of fDOMr is not well known but estimated to be  1 % (Hansell, 2013), which we tentatively adopt in MESMO 3. If DOMr is not selected in the model run, fDOMr=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):

(19) fDOM = 1.0 - min 0.72 , max 0.04 , 0.62 - 0.02 T C .

This formulation gives low export efficiency (i.e., high fDOM) in the warmer regions compared to the colder high-latitude 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 O2 directly to the dissolved inorganic forms (i.e., POC DIC, POP PO4, and PON NO3). In the new “deep POC split” pathway, sinking POM is simply broken down into DOM without the loss of O2 as in the production layer (Fig. 1). If DOMr is selected in the model, the broken-down POM is further routed to both DOMsl and DOMr according to fDOMr. If not, all of the broken down POM is converted to DOMsl. 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 DOMr in addition to DOM transport from the surface. Thus, the deep POC split pathway offers an alternative means to control deep ocean DOM distribution.

2.4 Production of CaCO3 and opal by eukaryotes

In MESMO 2, opal production was associated with the “large” PFT, and CaCO3 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 CaCO3 and opal production with the POP production by the same eukaryote PFT (JprodPOP1):


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 NO3 and Si(OH)4 into potential algal growth rates:

(22) RNPG = NO 3 NO 3 + K NO 3 , 1 - Si OH 4 Si OH 4 + K Si OH 4 .

If RNPG is more positive, the index indicates that nitrate-dependent growth exceeds silica-dependent growth. Thus, non-Si phytoplankton are more competitive, and this leads to higher CaCO3 production. On the other hand, a more negative RNPG implies that silica limitation for diatoms is relieved, leading to enhanced diatom growth and reduced CaCO3 production. The RNPG index is incorporated in the calculation of the rain ratio rCaCO3:POC presented in Eq. (20) as follows:

(23) r CaCO 3 : POC = r 0 CaCO 3 : POC Ω - 1 η min 1 , max 0.1 , RNPG k T , CaCO 3 .

Equation (23) indicates the base rain ratio r0CaCO3:POC (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):


Ksp is the solubility product of CaCO3. The temperature dependency of CaCO3 formation (kT,CaCO3) is similar to that of Moore et al. (2004) where warmer temperatures favor the growth of carbonate-bearing phytoplankton.

2.5 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):

(26) R POM i = V POM e k R T O 2 [ O 2 ] + K O 2 POM i .

VPOM is the base remineralization rate, kR expresses the temperature sensitivity of remineralization, and KO2 is half-saturation 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 DOMsl decay, which has been estimated previously to be  1.5 years (Hansell, 2013):

(27) R DOM sl = τ sl DOM sl .

All forms of DOMr also remineralize at the same rate in MESMO 3. In total, there are three optional, additive sinks of DOMr in the model: slow background decay, photodegradation, and degradation via hydrothermal vents (Fig. 1). Observations clearly indicate that the 14C age of deep-ocean DOCr is 103 years (e.g., Druffel et al., 1992), much older than DI14C. In addition, the deep ocean DOCr 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 DOMr 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 DOMr (e.g., Mopper et al., 1991). This process is believed to convert DOMr 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 DOMr (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 DOMr as a net sink (Lang et al., 2006). Here we assume simply that seawater that circulates through the vents loses all DOMr (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 DOMr is removed: the total removal rate depends on the vent flux of seawater Hflux. We implement the vent degradation of DOMr in MESMO 3 by first identifying the wet grid boxes located immediately above known mid-ocean ridges. We then distribute the annual global Hflux of 4.8×1016 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 (1021 kg vs. 1013 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 DOMr sinks are not mutually exclusive. They can thus be combined to yield the total DOMr remineralization rate:

(28) R DOM r = τ bg + τ photo + τ vent SW flux_local SW grid DOM r ,

where SWflux_local is the mass of seawater that circulates through the vents in each grid box in Δt, and SWgrid is the total mass of seawater in the same grid box.

The amount of O2 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-O2:POC and r-O2:DOC, respectively. These are molar ratios of O2 consumed per unit organic carbon respired. They are variable and calculated from the ambient POM and DOM concentration (Tanioka and Matsumoto, 2020b):


2.6 Remineralization of CaCO3 and opal

Remineralization of CaCO3 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 CaCO3 is a function of temperature similar to that of particulate organic matter remineralization but without oxygen dependency. The temperature-dependence term kR modifies the base remineralization rate VCaCO3:

(31) R CaCO 3 = V CaCO 3 e k R T CaCO 3 .

Opal remineralization in the water column follows Ridgwell et al. (2002). The rate of opal remineralization Ropal is given by the product of normalized dissolution rate (ropal), base opal dissolution rate (kopal), and opal concentration [opal].


ropal is a function of temperature (T) and the degree of under-saturation (uopal), 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:

(35) log 10 Si OH 4 eq = 6.44 - 968 T K .

Without the sediment module of MESMO activated, both CaCO3 and opal particles that reach the seafloor are completely dissolved back to inorganic forms.

2.7 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 POMi of all the PFTs give the total POM concentrations.


The time rate of change of CaCO3 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.


2.8 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 (Jprod) are zero below the upper-ocean production layer. Nutrients have a unit of mol element kg−1 in the model.


In Eq. (51), FixN is the N fixation carried out by diazotrophs, and DenN is the water column denitrification. There is an air–sea gas exchange term Fgas in Eqs. (52) and (56) for gaseous CO2 and O2, respectively. In Eq. (53), alkalinity increases with decreasing nitrate concentrations and increasing CaCO3 dissolution. Equation (54) contains RPOMFe, which is an iron source that represents remineralization of the Fe scavenged by sinking particles. These terms are explained in the following sections.

2.9 Prognostic nitrogen cycle

Biological production by diazotrophs is stimulated when the ambient NO3 is low. Nitrogen fixed by diazotrophs during their growth is added to the marine NO3 pool. The prognostic nitrogen fixation model employed here is similar to that used in the HAMOCC biogeochemical module (Paulsen et al., 2017):


where FixN is the nitrogen fixation rate and INO3 is the nitrate dependency term in quadratic Michaelis–Menten kinetics form with the half-saturation constant KNO3_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 NO3 are converted to 1 mol of N2 and liberating 2.5 mol of O2 as a byproduct:

(59) 2 NO 3 - + 2 H + 2.5 O 2 + N 2 + H 2 O .

Denitrification takes place in grid boxes, in which O2 concentration is below a threshold concentration (O2,def) and is stimulated if the total global inventory of NO3 relative to PO4 is high. In other words, denitrification can effectively act as negative feedback to nitrogen fixation. The threshold O2 concentration (O2,def) takes the minimum of the hard-bound O2 threshold concentration (O2,crit) and the NO3 to PO4 ratio, scaled by a parameter kD. The parameters O2,crit and kD 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.


2.10 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:

(62) K ligand = FeL Fe L ,

where [L] is the ligand concentration and Kligand 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 (Fin) following the prescribed modern flux pattern (Mahowald et al., 2006) with constant solubility (β):

(63) S Fe = β F in .

Particle-scavenged iron POMFe (note the difference from POFe) is produced below the productive layer when sinking POM scavenges Fe to sinking POM:

(64) J Fe = - τ sc K o POC 0.58 Fe ,

where τsc and Ko are empirical parameters that determine the strength of scavenging. Remineralization of Fe scavenged to POM (POMFe) is identical in form to that of POM remineralization:

(65) R POM Fe = V POM e k R T O 2 [ O 2 ] + K O 2 POM Fe .

The conservation equation of the particle scavenged iron is thus expressed as follows:

(66) POM Fe t = J Fe - z w POM Fe - R POM Fe .

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.

2.11 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 CO2 and O2 gases across the air–sea interface is given by


where k is the gas transfer velocity, CO2sat and O2sat 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:

(69) k = 0.31 u 2 S c 660 - 0.5 .
3 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 DOMr. The results from the standard model (hereafter just MESMO 3) are presented in Sect. 3.1, and the results from the DOMr-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.

Table 3Key Biogeochemical diagnostics

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 O2 distributions from MESMO 2 and 3 was calculated relative to the World Ocean Atlas 2018 (WOA18) gridded data (Garcia et al., 2018, 2019). The model–data comparison is made in the top 100 m for nutrients and below 100 m for O2. 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 (Dunne et al., 2007), (5) global CaCO3 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 O2 from WOA18 below 100 m (Garcia et al., 2019).

Download Print Version | Download XLSX

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 PO4, NO3, O2, 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 PO4 and NO3 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 O2 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 O2 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 O2 minimum depth would help simulate denitrification in the correct depth range, and there is a modest improvement in the data–model O2 mismatch in terms of RMSE in MESMO 3 over MESMO 2 (Table 3). The deepening of the O2 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]-[NO3]) 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: DOCsl, DOPsl, DONsl, and DOFesl. We note only that the surface DOCsl 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.

3.1 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 thus are more dominant in the oligotrophic subtropical gyres (Fig. 2c). The diazotrophs do not have NO3 limitation but have the lowest maximum growth rate. Thus it is much lower in abundance than the other two PFTs generally, and out-competed in transient blooms and thus excluded in higher latitudes (Fig. 2e).

Figure 2NPP-based surface phytoplankton functional type (PFT) abundance and nutrient limitation in MESMO 3. Fractional abundance and nutrient limitation for eukaryotes (a, b), cyanobacteria (c, d), and diazotrophs (e, f).

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 consistent 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 NO3 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 PO4 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 CaCO3 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, CaCO3 production is relatively large. The decoupling is prominent in the North Indian Ocean. Note that the spatial pattern of CaCO3 production is quite different in MESMO 3 (Fig. 3c) and MESMO 2 (Fig. 3d) because CaCO3 production was associated in MESMO 2 with the “small” PFT, which corresponds to the cyanobacteria PFT in MESMO 3.

Figure 3Eukaryote production in MESMO 3 and CaCO3 export in MESMO 2. In MESMO 3, eukaryote NPP (a) is linked to both opal export (b) and CaCO3 export (c) but the two export productions are differentiated by the residual nitrate potential growth (RNPG). Compare CaCO3 export in MESMO 3 (c) to MESMO 2 (d) (unit: mol m−2 yr−1).

The global pattern of the mean C:P uptake ratio in the production layer is shown in Fig. 4. Consistent with observations (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 PO4 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 sPO4P:C; see Eq. 5 and Table 2b). However, the cyanobacteria PFT's C:P ratio has an additional sensitivity to temperature (i.e., sTP:C0) 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).

Figure 4Uptake C:P ratio in the top 100 m in MESMO 3: (a) phytoplankton community C:P, (b) zonal mean C:P of all three PFTs and phytoplankton community, (c) eukaryote C:P, and (d) cyanobacteria C:P. The colors in (b) indicate community C:P (black), eukaryote C:P (red), cyanobacteria C:P (green), and diazotroph C:P (blue). In addition, (b) shows the mean range of observed C:P ratio binned by latitude (Martiny et al., 2013).

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 PO4 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 PO4 (Fig. 5c, d in black), but there is actually not a causal relationship (i.e., sPO4N:C=0, Table 2b). The C:N-PO4 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 PO4. For the cyanobacteria C:P ratio, the important drivers are temperature and PO4. For the C:N ratio for both eukaryotes and cyanobacteria, NO3 is more important than PAR. Figure 5 also serves to remind us that correlation does not indicate causation.

Figure 5Scatterplots of surface ocean eukaryote and cyanobacteria C:P and C:N vs. environmental drivers in MESMO 3. Columns show the following data, from left to right, eukaryote C:P, cyanobacteria C:P, eukaryote C:N, and cyanobacteria C:N. Rows show the following data, from top to bottom, PO4, NO3, temperature, and PAR. Red indicates the causal relationship according to the power law formulation of flexible C:N:P ratio. PAR stands for photosynthetically active radiation in W m−2.


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 PO4 and temperature are clearly visible on the community C:P ratio. Both low [PO4] and warmer waters are found in the lower latitudes, so the P frugality and temperature effects are additive. The effect of NO3 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 PFT-specific C:N:P are obvious in the community C:N:P ratio.

Figure 6Scatterplots of surface ocean community C:P and C:N vs. environmental drivers in MESMO 3.


3.2 DOMr-enabled MESMO 3

In MESMO 2, DOCsl 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., DOPsl, DONsl, and DOFesl, and the refractory forms of DOM, i.e., DOCr, DOPr, and DONr. MESMO 3 is not yet calibrated with respect 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 DOMr 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 fDOMr=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 (DOCt = DOCsl + DOCr) with solid lines and DOCr with a dashed line. Qualitatively, the simulated profiles are consistent with the observations, showing a near-uniform DOCr concentration and a DOCsl 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 DOCr is 30–40 µmol kg−1, and the observed DOCsl attenuates with depth from 30–40 µmol kg−1 near the surface. So their sum, which is represented by DOCt, is approximately 60–80 µmol kg−1 at the surface in observations.

Figure 7Global mean vertical profiles of DOC from the DOMR-enabled MESMO 3. DOCt (DOCsl+DOCr, black line) and DOCr (dashed black line) from the LV run. The red line is DOCt after reducing fDOMr from 1 % in LV to 0.2 % (Experiment 210310o) (unit: µmol kg−1).


Figure 8 adds a lateral perspective to Fig. 7. The rapid DOCt attenuation in the vertical is strong in the lower latitudes where stratification is generally stronger. The transport of DOCsl from the surface to deeper waters is evident in the high latitudes of the North Atlantic and the Southern Ocean. The DOCt change in the deep ocean is limited.

Figure 8Global depth–latitude transect of DOCt from the DOMR-enabled MESMO 3 LV run. Transects are north–south along 25 W in the Atlantic, east–west along 60 S in the Southern Ocean, and north–south along 165 E in the Pacific (unit: µmol kg−1).


Observations of deep-ocean DOCt 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 DOCt 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 DOCt distributions from the LV run can also be compared to a global extrapolation based on an artificial neural network (ANN) of the available DOCt data (Roshan and DeVries, 2017). At the surface, the extrapolation indicates higher DOCt concentrations in the subtropical gyres (Fig. 9a), while our simulation does not clearly delineate the gyres (Fig. 9c). In our model, fDOM is temperature-dependent and strongly controls the production of DOM. The surface DOCt 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 DOCt show a gradual decline in concentrations from the North Atlantic to the North Pacific (Fig. 9b, d). The highest deep DOCt in the LV run is seen just south of Greenland, where convection occurs in the model.

Figure 9Assessment of surface and deep-ocean DOCt from the DOMR-enabled MESMO 3 LV run. Data-derived DOCt distributions in the top 100 m (a) and 2500–4000 m (b). Model-simulated DOCt distributions in the top 100 m (c) and 2500–4000 m (d). Data-derived DOCt are from Roshan and DeVries (Roshan and DeVries, 2017) (unit: µmol kg−1).

Finally, we show that the deep-ocean radiocarbon aging is larger in DIC than in DOCt in the model (Fig. 10). The North Pacific–North Atlantic Δ14C gradient is roughly 100 ‰ for DIC and 70 ‰ for DOCt. The oldest DOCtΔ14C is approximately 430 ‰ in the North Pacific. If 14C decay were the only mechanism of change along the path of the deepwater circulation, the Δ14C gradient should be quite similar between DIC and DOCt, 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 DI14C and DO14C to the deep ocean by the “deep POC split” (see Sect. 2.3). The addition impacts DOCtΔ14C more than DIC Δ14C because DOCt is 2 orders of magnitude lower in concentration than DIC.

Figure 10Δ14C of deep-ocean DIC (a) and DOCt (b) from the DOMR-enabled MESMO 3 LV run. Vertical average over 2500–4000 m water depth (unit: ‰).

In observations, the aging of DIC and DOCt 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 Δ14C gradient, accounting for thermonuclear bomb 14C, may be as large as 100 ‰ for DOCt (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 Δ14C gradient for DIC cannot be determined. However, the DIC Δ14C endmember values by inspection (about 250 ‰ in the deep Pacific and 70 ‰ in the deep Atlantic) (Matsumoto and Key, 2004) indicate a clearly larger Δ14C gradient for DIC than DOCt as simulated by the experiment LV.

One lesson from the data–LV run mismatch in the overall DOCt concentration (Fig. 7) and surface DOCt 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 DOCt is too high because DOCr is too high, while DOCsl is not unreasonable (Fig. 7). DOCr is too high because there is too much DOCr production (e.g., fDOMr=1 % is too large), there is too little DOCr degradation (e.g., one of the DOM decay mechanisms is too slow; Eq. 28 and Table 2d), or some combination of both. For example, fDOMr is a key parameter that is not well constrained by observations. Had we used 0.2 % instead of 1 % for fDOMr, the global mean surface DOCt drops to 76 µmol kg−1 (red line, Fig. 7), consistent with observations. For achieving a better surface DOCt 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 DOPr 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 DOPr concentration becomes steady at 0.45 µmol-P kg−1. In observations, the mean DOCr is about 40 µmol-C kg−1 and the DOCr:DOPr ratio is estimated to be 1370:1 (Letscher and Moore, 2015), so DOPr concentration should only be roughly 0.03 µmol-P kg−1. Thus, the simulated DOPr=0.45 µmol-P kg−1 is an order of magnitude too high. Because there is more P in the form of DOPr in LV, the oceanic inventory of PO4 declines, causing a nearly 10 % drop in export production compared to the standard MESMO 3. In LV, the decline in the surface ocean PO4 that accompanies the change in the PO4 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 DOPr concentrations.

3.3 Large-scale patterns of N2 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 N2 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, N2 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 N2 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 N2 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 (Wang et al., 2019), which showed an elevated N2 fixation rate in subtropical gyres.

Global water column denitrification is 101 Tg N yr−1 (Table 3) and is equal to the global N2 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 NO3 and PO4 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, KNO3_Nfix in Eq. (58). A large value of KNO3_Nfix will make it hard for diazotrophs to obtain fixed N from NO3, which would facilitate N2 fixation and pushes up the global N/P ratio. With a smaller value of KNO3_Nfix, diazotrophs will more easily uptake NO3, thus depressing N2 fixation and lowering the global N/P ratio.

4 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 (last access: 26 April 2021) and have the following DOI: (Matsumoto, 2020).


The supplement related to this article is available online at:

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.


Numerical modeling and analysis were carried out using resources at the University of Minnesota Supercomputing Institute.

Financial support

This research has been supported by the National Science Foundation (grant no. OCE-1827948).

Review statement

This paper was edited by Andrew Yool and reviewed by two anonymous referees.


Anderson, L. A. and Sarmiento, J. L.: Redfield ratios of remineralization determined by nutrient data analysis, Global Biogeochem. Cy., 8, 65–80,, 1994. 

Archer, D., Eby, M., Brovkin, V., Ridgwell, A., Cao, L., Mikolajewicz, U., Caldeira, K., Matsumoto, K., Munhoven, G., Montenegro, A., and Tokos, K.: Atmospheric Lifetime of Fossil Fuel Carbon Dioxide, Annu. Rev. Earth Pl. Sc., 37, 117–134,, 2009. 

Arteaga, L., Pahlow, M., and Oschlies, A.: Global patterns of phytoplankton nutrient and light colimitation inferred from an optimality-based model, Global Biogeochem. Cy., 28, 648–661,, 2014. 

Balch, W. M., Bates, N. R., Lam, P. J., Twining, B. S., Rosengard, S. Z., Bowler, B. C., Drapeau, D. T., Garley, R., Lubelczyk, L. C., Mitchell, C., and Rauschenberg, S.: Factors regulating the Great Calcite Belt in the Southern Ocean and its biogeochemical significance, Global Biogeochem. Cy., 30, 1124–1144,, 2016. 

Berelson, W. M., Balch, W. M., Najjar, R., Feely, R. A., Sabine, C., and Lee, K.: Relating estimates of CaCO3 production, export, and dissolution in the water column to measurements of CaCO3 rain into sediment traps and dissolution on the sea floor: A revised global carbonate budget, Global Biogeochem. Cy., 201, 1–15,, 2007. 

Bianchi, D., Weber, T. S., Kiko, R., and Deutsch, C.: Global niche of marine anaerobic metabolisms expanded by particle microenvironments, Nat. Geosci., 11, 263–268,, 2018. 

Brzezinski, M. A.: A switch from Si(OH)4 to NO3- depletion in the glacial Southern Ocean, Geophys. Res. Lett., 29, 1564,, 2002. 

Cao, L., Eby, M., Ridgwell, A., Caldeira, K., Archer, D., Ishida, A., Joos, F., Matsumoto, K., Mikolajewicz, U., Mouchet, A., Orr, J. C., Plattner, G.-K., Schlitzer, R., Tokos, K., Totterdell, I., Tschumi, T., Yamanaka, Y., and Yool, A.: The role of ocean transport in the uptake of anthropogenic CO2, Biogeosciences, 6, 375–390,, 2009. 

Carr, M. E., Friedrichs, M. A. M., Schmeltz, M., Noguchi Aita, M., Antoine, D., Arrigo, K. R., Asanuma, I., Aumont, O., Barber, R., Behrenfeld, M., Bidigare, R., Buitenhuis, E. T., Campbell, J., Ciotti, A., Dierssen, H., Dowell, M., Dunne, J., Esaias, W., Gentili, B., Gregg, W., Groom, S., Hoepffner, N., Ishizaka, J., Kameda, T., Le Quéré, C., Lohrenz, S., Marra, J., Mélin, F., Moore, K., Morel, A., Reddy, T. E., Ryan, J., Scardi, M., Smyth, T., Turpie, K., Tilstone, G., Waters, K., and Yamanaka, Y.: A comparison of global estimates of marine primary production from ocean color, Deep-Sea Res. Pt. II, 53, 741–770,, 2006. 

Cheng, H., Edwards, R. L., Southon, J., Matsumoto, K., Feinberg, J. M., Sinha, A., Zhou, W., Li, H., Li, X., Xu, Y., Chen, S., Tan, M., Wang, Q., Wang, Y., Ning, Y., Lawrence Edwards, R., Southon, J., Matsumoto, K., Feinberg, J. M., Sinha, A., Zhou, W., Li, H., Li, X., Xu, Y., Chen, S., Tan, M., Wang, Q., Wang, Y., and Ning, Y.: Atmospheric 14C/12C changes during the last glacial period from hulu cave, Science, 362, 1293–1297,, 2018. 

Chien, C.-T., Pahlow, M., Schartau, M., and Oschlies, A.: Optimality-based non-Redfield plankton–ecosystem model (OPEM v1.1) in UVic-ESCM 2.9 – Part 2: Sensitivity analysis and model calibration, Geosci. Model Dev., 13, 4691–4712,, 2020. 

Claussen, M., Mysak, L., Weaver, A., Crucifix, M., Fichefet, T., Loutre, M. F., Weber, S., Alcamo, J., Alexeev, V., Berger, A., Calov, R., Ganopolski, A., Goosse, H., Lohmann, G., Lunkeit, F., Mokhov, I., Petoukhov, V., Stone, P., and Wang, Z.: Earth system models of intermediate complexity: Closing the gap in the spectrum of climate system models, Clim. Dynam., 18, 579–586,, 2002. 

Codispoti, L. A.: An oceanic fixed nitrogen sink exceeding 400 Tg N a−1 vs the concept of homeostasis in the fixed-nitrogen inventory, Biogeosciences, 4, 233–253,, 2007. 

Deutsch, C., Sarmiento, J. L., Sigman, D. M., Gruber, N., and Dunne, J. P.: Spatial coupling of nitrogen inputs and losses in the ocean, Nature, 445, 163–167,, 2007. 

DeVries, T. and Weber, T.: The export and fate of organic matter in the ocean: New constraints from combining satellite and oceanographic tracer observations, Global Biogeochem. Cy., 31, 535–555,, 2017. 

Druffel, E. R. M., Williams, P. M., Bauer, J. E., and Ertel, J. R.: Cycling of dissolved and particulate organic matter in the open ocean, J. Geophys. Res., 97, 15639,, 1992. 

Druffel, E. R. M., Griffin, S., Wang, N., Garcia, N. G., McNichol, A. P., Key, R. M., and Walker, B. D.: Dissolved Organic Radiocarbon in the Central Pacific Ocean, Geophys. Res. Lett., 46, 5396–5403,, 2019. 

Dunne, J. P., Armstrong, R. A., Gnanadesikan, A., and Sarmiento, J. L.: Empirical and mechanistic models for the particle export ratio, Global Biogeochem. Cy., 19, GB4026,, 2005. 

Dunne, J. P., Sarmiento, J. L. ,and Gnanadesikan, A.: A synthesis of global particle export from the surface ocean and cycling through the ocean interior and on the seafloor, Global Biogeochem. Cy., 21, 1–16,, 2007. 

Eby, M., Weaver, A. J., Alexander, K., Zickfeld, K., Abe-Ouchi, A., Cimatoribus, A. A., Crespin, E., Drijfhout, S. S., Edwards, N. R., Eliseev, A. V., Feulner, G., Fichefet, T., Forest, C. E., Goosse, H., Holden, P. B., Joos, F., Kawamiya, M., Kicklighter, D., Kienert, H., Matsumoto, K., Mokhov, I. I., Monier, E., Olsen, S. M., Pedersen, J. O. P., Perrette, M., Philippon-Berthier, G., Ridgwell, A., Schlosser, A., Schneider von Deimling, T., Shaffer, G., Smith, R. S., Spahni, R., Sokolov, A. P., Steinacher, M., Tachiiri, K., Tokos, K., Yoshimori, M., Zeng, N., and Zhao, F.: Historical and idealized climate model experiments: an intercomparison of Earth system models of intermediate complexity, Clim. Past, 9, 1111–1140,, 2013. 

Edwards, N. R. and Marsh, R.: Uncertainties due to transport-parameter sensitivity in an efficient 3-D ocean-climate model, Clim. Dynam., 24, 415–433,, 2005. 

Eppley, R. W.: Temperature and phytoplankton growth in the sea, Fish. Bull., 70, 1063–1085, 1972. 

Galbraith, E. D. and Martiny, A. C.: A simple nutrient-dependence mechanism for predicting the stoichiometry of marine ecosystems, P. Natl. Acad. Sci. USA, 112, 8199–8204,, 2015. 

Garcia, H., Weathers, K. W., Paver, C. R., Smolyar, I., Boyer, T. P., Locarnini, R. A., Zweng, M. M., Mishonov, A. V., Baranova, O. K., Seidov, D., and Reagan, J. R.: World Ocean Atlas 2018. Volume 4: Dissolved Inorganic Nutrients (phosphate, nitrate and nitrate+nitrite, silicate), NOAA Atlas NESDIS, 84, 35 pp., 2018. 

Garcia, H. E., Weathers, K., Paver, C. R., Smolyar, I., Boyer, T. P., Locarnini, R. A., Zweng, M. M., Mishonov, A. V., Baranova, O. K., Seidov, D., and Reagan, J. R.: World Ocean Atlas 2018, Volume 3: Dissolved Oxygen, Apparent Oxygen Utilization, and Oxygen Saturation, NOAA Atlas NESDIS, 3, 38 pp., 2019. 

Gruber, N. and Sarmiento, J. L.: Global patterns of marine nitrogen fixation and denitrification, Global Biogeochem. Cy., 11, 235–266,, 1997. 

Hansell, D. A.: Recalcitrant Dissolved Organic Carbon Fractions, Ann. Rev. Mar. Sci., 5, 421–445,, 2013. 

Hansell, D. A. and Carlson, C. A.: Deep-ocean gradients in the concentration of dissolved organic carbon, Nature, 395, 263–266,, 1998. 

Hansell, D. A., Carlson, C. A., Repeta, D. J., and Schlitzer, R.: Dissolved Organic Matter in the Ocean: A Controversy stimulates new insights, Oceanography, 22, 202–211,, 2009. 

Henson, S. A., Sanders, R., Madsen, E., Morris, P. J., Le Moigne, F., and Quartly, G. D.: A reduced estimate of the strength of the ocean's biological carbon pump, Geophys. Res. Lett., 38, L04606,, 2011. 

Hutchins, D. A. and Bruland, K. W.: Iron-limited growth and Si:N ratios in a costal upwelling regime, Nature, 393, 561–564, 1998. 

Joos, F., Roth, R., Fuglestvedt, J. S., Peters, G. P., Enting, I. G., von Bloh, W., Brovkin, V., Burke, E. J., Eby, M., Edwards, N. R., Friedrich, T., Frölicher, T. L., Halloran, P. R., Holden, P. B., Jones, C., Kleinen, T., Mackenzie, F. T., Matsumoto, K., Meinshausen, M., Plattner, G.-K., Reisinger, A., Segschneider, J., Shaffer, G., Steinacher, M., Strassmann, K., Tanaka, K., Timmermann, A., and Weaver, A. J.: Carbon dioxide and climate impulse response functions for the computation of greenhouse gas metrics: a multi-model analysis, Atmos. Chem. Phys., 13, 2793–2825,, 2013. 

Kwiatkowski, L., Aumont, O., Bopp, L., and Ciais, P.: The Impact of Variable Phytoplankton Stoichiometry on Projections of Primary Production, Food Quality, and Carbon Uptake in the Global Ocean, Global Biogeochem. Cy., 516–528,, 2018. 

Kwiatkowski, L., Aumont, O., and Bopp, L.: Consistent trophic amplification of marine biomass declines under climate change, Glob. Chang. Biol., 25, 218–229,, 2019. 

Landolfi, A., Kähler, P., Koeve, W., and Oschlies, A.: Global marine N2 fixation estimates: From observations to models, Front. Microbiol., 9, 1–8,, 2018. 

Lang, S. Q., Butterfield, D. A., Lilley, M. D., Paul Johnson, H., and Hedges, J. I.: Dissolved organic carbon in ridge-axis and ridge-flank hydrothermal systems, Geochim. Cosmochim. Ac., 70, 3830–3842,, 2006. 

Laufkötter, C., John, J. G., Stock, C. A., and Dunne, J. P.: Temperature and oxygen dependence of the remineralization of organic matter, Global Biogeochem. Cy., 31, 1038–1050,, 2017. 

Laws, E. A., Falkowski, P. G., Smith, W. O., Ducklow, H. W., and McCarthy, J. J.: Temperature effects on export production in the open ocean, Global Biogeochem. Cy., 14, 1231–1246,, 2000. 

Lee, S. Y., Chiang, J. C. H., Matsumoto, K. and Tokos, K. S.: Southern Ocean wind response to North Atlantic cooling and the rise in atmospheric CO2: Modeling perspective and paleoceanographic implications, Paleoceanography, 26, 1–16,, 2011. 

Lenton, T. M., Williamson, M. S., Edwards, N. R., Marsh, R., Price, A. R., Ridgwell, A. J., Shepherd, J. G., and Cox, S. J.: Millennial timescale carbon cycle and climate change in an efficient Earth system model, Clim. Dynam., 26, 687–711,, 2006. 

Letscher, R. T. and Moore, J. K.: Preferential remineralization of dissolved organic phosphorus and non-Redfield DOM dynamics in the global ocean: Impacts on marine productivity, nitrogen fixation, and carbon export, Global Biogeochem. Cy., 29, 325–340,, 2015. 

Levitus, S.: Climatological Atlas of the World Ocean, NOAA Professional Paper 13, 173 pp., 1982. 

Mahowald, N. M., Muhs, D. R., Levis, S., Rasch, P. J., Yoshioka, M., Zender, C. S., and Luo, C.: Change in atmospheric mineral aerosols in response to climate: Last glacial period, preindustrial, modern, and doubled carbon dioxide climates, J. Geophys. Res.-Atmos., 111,, 2006. 

Marchal, O., Stocker, T. F. and Joos, F.: A latitude-depth, circulation-biogeochemical ocean model for palaeoclimate studies. Development and sensitivities, Tellus B, 50B, 290–316,, 1998. 

Martiny, A. C., Pham, C. T. A. A., Primeau, F. W., Vrugt, J. A., Moore, J. K., Levin, S. A., and Lomas, M. W.: Strong latitudinal patterns in the elemental ratios of marine plankton and organic matter, Nat. Geosci., 6, 279–283,, 2013. 

Mather, R. L., Reynolds, S. E., Wolff, G. A., Williams, R. G., Torres-Valdes, S., Woodward, E. M. S., Landolfi, A., Pan, X., Sanders, R., and Achterberg, E. P.: Phosphorus cycling in the North and South Atlantic Ocean subtropical gyres, Nat. Geosci., 1, 439–443,, 2008. 

Matsumoto, K.: gaia3intc/mesmo: Release of MESMO v3.0 (Version v3.0), Zenodo,, 2020. 

Matsumoto, K. and Key, R. M.: Natural radiocarbon distribution in the deep ocean, in Global Environmental Change in the Ocean and on Land, edited by: Shiyomi, M., Kawahata, H., Koizumi, H., Tsuda, A., and Awaya, Y., 45–58, Terrapub, Tokyo, available at: (last access: 23 April 2021), 2004. 

Matsumoto, K. and McNeil, B.: Decoupled response of ocean acidification to variations in climate sensitivity, J. Climate, 26, 1764–1771,, 2012. 

Matsumoto, K. and Tanioka, T.: Shifts in regional production as a driver of future global ocean production stoichiometry, Environ. Res. Lett., 15, 124027,, 2020. 

Matsumoto, K. and Yokoyama, Y.: Atmospheric Δ 14 C reduction in simulations of Atlantic overturning circulation shutdown, Global Biogeochem. Cy., 27, 296–304,, 2013. 

Matsumoto, K., Tokos, K. S., Price, A. R., and Cox, S. J.: First description of the Minnesota Earth System Model for Ocean biogeochemistry (MESMO 1.0), Geosci. Model Dev., 1, 1–15,, 2008. 

Matsumoto, K., Tokos, K., Chikamoto, M., and Ridgwell, A.: Characterizing post-industrial changes in the ocean carbon cycle in an Earth system model, Tellus B, 62, 296–313,, 2010. 

Matsumoto, K., Tokos, K., Huston, A., and Joy-Warren, H.: MESMO 2: a mechanistic marine silica cycle and coupling to a simple terrestrial scheme, Geosci. Model Dev., 6, 477–494,, 2013. 

Matsumoto, K., Rickaby, R. and Tanioka, T.: Carbon Export Buffering and CO2 Drawdown by Flexible Phytoplankton C:N:P Under Glacial Conditions, Paleoceanogr. Paleocl., 35, 1–22,, 2020. 

Moore, J. K., Doney, S. C., and Lindsay, K.: Upper ocean ecosystem dynamics and iron cycling in a global three-dimensional model, Global Biogeochem. Cy., 18, 1–21,, 2004. 

Mopper, K., Zhou, X., Kieber, R. J., Kieber, D. J., Sikorski, R. J., and Jones, R. D.: Photochemical degradation of dissolved organic carbon and its impact on the oceanic carbon cycle, Nature, 353, 60–62,, 1991. 

Najjar, R. G., Sarmiento, J. L., and Toggweiler, J. R.: Downward transport and fate of organic matter in the ocean: Simulations with a general circulation model, Global Biogeochem. Cy., 6, 45–76,, 1992. 

Pahlow, M. and Oschlies, A.: Chain model of phytoplankton P, N and light colimitation, Mar. Ecol. Prog. Ser., 376, 69–83,, 2009. 

Pahlow, M. and Oschlies, A.: Optimal allocation backs droop's cell-quota model, Mar. Ecol. Prog. Ser., 473, 1–5,, 2013. 

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

Pahlow, M., Chien, C.-T., Arteaga, L. A., and Oschlies, A.: Optimality-based non-Redfield plankton–ecosystem model (OPEM v1.1) in UVic-ESCM 2.9 – Part 1: Implementation and model behaviour, Geosci. Model Dev., 13, 4663–4690,, 2020. 

Paulsen, H., Ilyina, T., Six, K. D., and Stemmler, I.: Incorporating a prognostic representation of marine nitrogen fixers into the global ocean biogeochemical model HAMOCC, J. Adv. Model. Earth Sy., 9, 438–464,, 2017. 

Ragueneau, O., Tréguer, P., Leynaert, A., Anderson, R., Brzezinski, M., DeMaster, D., Dugdale, R., Dymond, J., Fischer, G., François, R., Heinze, C., Maier-Reimer, E., Martin-Jézéquel, V., Nelson, D., and Quéguiner, B.: A review of the Si cycle in the modern ocean: recent progress and missing gaps in the application of biogenic opal as a paleoproductivity proxy, Global Planet. Change, 26, 317–365,, 2000. 

Ridgwell, A.: Glacial-interglacial perturbations in the global carbon cycle, PhD thesis, Univ. East Anglia, Norwich, U. K. Ridgwell, A. J., U. Edwards, 134 pp., 2001. 

Ridgwell, A., Hargreaves, J. C., Edwards, N. R., Annan, J. D., Lenton, T. M., Marsh, R., Yool, A., and Watson, A.: Marine geochemical data assimilation in an efficient Earth System Model of global biogeochemical cycling, Biogeosciences, 4, 87–104,, 2007. 

Ridgwell, A. J., Watson, A. J., and Archer, D. E.: Modeling the response of the oceanic Si inventory to perturbation, and consequences for atmospheric CO 2, Global Biogeochem. Cy., 16, 19-1–19-25,, 2002. 

Roshan, S. and DeVries, T.: Efficient dissolved organic carbon production and export in the oligotrophic ocean, Nat. Commun., 8, 2036,, 2017. 

Sarmiento, J. L., Gruber, N., Brzezinski, M. A., and Dunne, J. P.: High-latitude controls of thermocline nutrients and low latitude biological productivity, Nature, 427, 56–60,, 2004. 

Sun, X. and Matsumoto, K.: Effects of sea ice on atmospheric pCO2: A revised view and implications for glacial and future climates, J. Geophys. Res., 115, G02015,, 2010. 

Sunda, W. G. and Huntsman, S. A.: Iron uptake and growth limitation in oceanic and coastal phytoplankton, Mar. Chem., 50, 189–206,, 1995. 

Sverdrup, H. U.: On the conditions for the vernal blooming of phytoplankton, J. Cons. Perm. Int. Pour l'Exploration La Mer, 18, 287–195, 1953. 

Takeda, S.: Influence of iron availability on nutrient consumption ratio, Nature, 393, 774–777, 1998. 

Tanioka, T. and Matsumoto, K.: Buffering of Ocean Export Production by Flexible Elemental Stoichiometry of Particulate Organic Matter, Global Biogeochem. Cy., 31, 1528–1542,, 2017. 

Tanioka, T. and Matsumoto, K.: A meta-analysis on environmental drivers of marine phytoplankton C : N : P, Biogeosciences, 17, 2939–2954,, 2020a. 

Tanioka, T. and Matsumoto, K.: Stability of Marine Organic Matter Respiration Stoichiometry, Geophys. Res. Lett., 47, 1–10,, 2020b. 

Ushie, H. and Matsumoto, K.: The role of shelf nutrients on glacial-interglacial CO2: A negative feedback, Global Biogeochem. Cy., 26, 1–10,, 2012. 

Wang, W.-L., Moore, J. K., Martiny, A. C., and Primeau, F. W.: Convergent estimates of marine nitrogen fixation, Nature, 566, 205–211,, 2019. 

Wanninkhof, R.: Relationship between wind speed and gas exchange over the ocean, J. Geophys. Res., 97, 7373–7382,, 1992. 

Weaver, A. J., Sedlá, J., Eby, M., Alexander, K., Crespin, E., Fichefet, T., Philippon-berthier, G., Joos, F., Kawamiya, M., Matsumoto, K., Steinacher, M., Tachiiri, K., Tokos, K., Yoshimori, M., and Zickfeld, K.: Stability of the Atlantic meridional overturning circulation: A model intercomparison, Geophys. Res. Lett., 39, 1–7,, 2012. 

Yamanaka, Y. and Tajika, E.: Role of dissolved organic matter in the marine biogeochemical cycle: Studies using an ocean biogeochemical general circulation model, Global Biogeochem. Cy., 11, 599–612,, 1997. 

Yamanaka, Y., Yoshie, N., Fujii, M., Aita, M. N., and Kishi, M. J.: An Ecosystem Model Coupled with Nitrogen-Silicon-Carbon Cycles Applied to Station A7 in the Northwestern Pacific, J. Oceanogr., 60, 227–241,, 2004. 

Zickfeld, K., Eby, M., Weaver, A. J., Alexander, K., Crespin, E., Edwards, N. R., Eliseev, A. V., Feulner, G., Fichefet, T., Forest, C. E., Friedlingstein, P., Goosse, H., Holden, P. B., Joos, F., Kawamiya, M., Kicklighter, D., Kienert, H., Matsumoto, K., Mokhov, I. I., Monier, E., Olsen, S. M., Pedersen, J. O. P., Perrette, M., Philippon-Berthier, G., Ridgwell, A., Schlosser, A., Von Deimling, T. S., Shaffer, G., Sokolov, A., Spahni, R., Steinacher, M., Tachiiri, K., Tokos, K. S., Yoshimori, M., Zeng, N., and Zhao, F.: Long-Term climate change commitment and reversibility: An EMIC intercomparison, J. Climate, 26, 5782–5809,, 2013. 

Short summary
MESMO is a mathematical model that represents the essential components of the Earth, such as the global ocean, atmosphere, and sea ice. It is used to study the global climate and the global carbon cycle. We describe the third version of MESMO. A novel feature of the new version is its mathematical representations of the chemical composition of marine phytoplankton and the marine dissolved organic matter, which are both recognized as important for the global ocean carbon cycle.