MESMO 3 : Flexible phytoplankton stoichiometry and refractory DOM 1 2

Abstract. 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, a 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. Also, we describe the semi-labile and refractory dissolved organic pools of C, N, P, and Fe that can be enabled in MEMSO 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 30
Here we document the development of the third version of the Minnesota Earth System 31 Model for Ocean biogeochemistry (MESMO 3). As described for the first two versions 32 (Matsumoto et al., 2008, MESMO is based on the non-modular version of the Grid 33 ENabled Integrated Earth (GENIE) system model (Lenton et al., 2006;Ridgwell et al., 2007). 34 The computationally efficient ocean-climate model of Edwards  where I is the seasonally variable solar short-wave irradiance in W m -2 . Light is attenuated 140 exponentially from the ocean surface with a 20 m depth scale. 141 142 Nutrient uptake in Equation 1 has a dependence on zml, which is diagnosed using the st 143 density gradient criterion (Levitus, 1982). Following the Sverdrup (1953)  Irradiance (I): 160 Equations 5 and 6 are the power-law equations that calculate the change in P:C and N:C for 163 fractional changes in environmental drivers relative to the reference P:C and N:C, 164 respectively (Matsumoto et al., 2020;Tanioka and Matsumoto, 2017

Production of POM and DOM 227
In the top 100 m of the model domain, where phytoplankton P uptake occurs (i.e., Γ ! > 0, 228 see section 2.1), NPP is immediately routed to POM and DOM pools (Figure 1). The 229 production fluxes of POM, DOMsl, and DOMr from NPP are given as Jprod. Here we write the 230 equations in terms of P, which is the master nutrient variable: 231 232 The term fDOM denotes the fraction of NPP that is routed to DOM as opposed to POM. 234 Likewise, fDOMr is the fraction of DOM that is routed to DOMr as opposed to DOMsl. The 235 value of fDOMr is not well known but estimated to be ~1% , which we 236 tentatively adopt in MESMO 3. If DOMr is not selected in the model run, fDOMr = 0. In 237 previous versions of MESMO, fDOM was assigned a constant value of 0.67. In reality, a large 238 variability is observed locally for this ratio, ranging from 0.01-0.2 in temperate waters to In MEMSO 3, a new DOM production pathway below the production layer is available as an 250 option. In previous MESMO versions, sinking POM was respired in the water column with 251 the loss of O2 directly to the dissolved inorganic forms (i.e., POC-->DIC, POP-->PO4, and 252 POP-->NO3). In the new "deep POC split" pathway, sinking POM is simply broken down into 253 DOM without the loss of O2 as in the production layer ( Figure 1). If DOMr is selected in the 254 model, the broken-down POM is further routed to both DOMsl and DOMr according to

Production of CaCO3 and opal by eukaryotes 262
In MESMO 2, opal production was associated with the large PFT and CaCO3 production was 263 associated with the "small" PFT. We recognize that coccolithophorids and diatoms, which 264 are the producers of these biogenic tests, are both eukaryotes. Therefore, in MEMSO 3, we 265 associate both CaCO3 and opal production with the POM production by the same eukaryote allowing competition between diatoms and non-siliceous phytoplankton within the same 271 If RNPG is more positive, the index indicates that nitrate-dependent growth exceeds silica-279 dependent growth. Thus, non-Si phytoplankton are more competitive, and this leads to 280 higher CaCO3 production. On the other hand, a more negative RNPG implies that silica 281 limitation for diatoms is relieved, leading to enhanced diatom growth and reduced CaCO3 26 303 V @LV is the base remineralization rate, parameter k W expresses the temperature sensitivity 304 of remineralization, and K O 2 is half-saturation constant for oxygen-dependent 305 remineralization. When the sediment model is not coupled, any POM that reaches the 306 seafloor dissolves completely to its inorganic form and is returned to the overlying water. All forms of DOMr also remineralize at the same rate in MESMO 3. In total, there are three 314 optional, additive sinks of DOMr in the model: slow background decay, photodegradation, 315 and degradation via hydrothermal vents ( Figure 1). Observations clearly indicate that the 316 14 C age of deep ocean DOCr is 10 3 years (e.g., Druffel et al., 1992), much older than DI 14 C. 317 Also, the deep ocean DOCr concentration decreases modestly along the path of the deep 318 water from the deep North Atlantic to the deep North Pacific (Hansell and Carlson, 1998). 319 Thus, it is understood that there is a slow DOMr background decay in the deep ocean. We 320 represent this ubiquitous process with tbg, which is the inverse of the background decay 321 time scale, estimated to be ~16,000 years . 322

323
Observations to date indicate that photodegradation is a major sink of DOMr (e.g., Mopper interior into the euphotic zone into more labile forms of DOM. We represent 326 photodegradation with tphoto, the inverse of the decay time scale, estimated to be ~70 years 327 (Yamanaka and Tajika, 1997). This occurs only in the surface. The three DOMr sinks are not mutually exclusive. They can thus be combined to yield the 346 total DOMr remineralization rate: 347 where SWflux_local is the mass of seawater that circulates through the vents in each grid box 350 in Dt, and SWgrid is the total mass of seawater in the same grid box. 351 The amount of O2 respired as a result of these POM and DOM remineralization processes is 353 related to the organic carbon pools by the respiratory quotients of POC and DOC, r 2L ( :@L? 354 and r 2L ( :`L? , respectively. These are molar ratios of O2 consumed per unit organic carbon 355 respired. They are variable and calculated from the ambient POM and DOM concentration 356

Remineralization of CaCO3 and opal 360
Remineralization of CaCO3 and opal particles occurs as they sink through the water column 361 and remains the same as in MESMO 2. Key parameter values are given in Table 2d. 362 Remineralization of CaCO3 is a function of temperature similar to that of particulate organic 363 matter remineralization but without oxygen dependency. The temperature dependence 364 term kR modifies the base remineralization rate V ?Q?L & : 365 Opal remineralization in the water column follows Ridgwell et al. (2002). The rate of opal 368 remineralization R abQc is given by the product of normalized dissolution rate (r abQc ), base 369 opal dissolution rate (k abQc ), and opal concentration [opal]: 370

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

Conservation of inorganic nutrients 396
The time rate of change of the inorganic nutrients have organic carbon production as sink 397 terms and remineralization as source terms. The production terms (Jprod) are zero below 398 the upper ocean production layer. Nutrients have a unit of mol element kg -1 in the model. 399 where FixN is the nitrogen fixation rate and I DL & is the nitrate dependency term in 416 quadratic Michaelis-Menten kinetics form with the half-saturation constant K D ( . See Table  417 2e for the values related to the N cycle. 418 419 Water-column denitrification is formulated in an approach similar to that of the original 420 GENIE model (Ridgwell et al., 2007), in which 2 moles of NO3 are converted to 1 mole of N2 421 and liberating 2.5 moles of O2 as a byproduct: 422 Denitrification takes place in grid boxes, in which O2 concentration is below a threshold 425 concentration (O 9,gFh ) and is stimulated if the total global inventory of NO3 relative to PO4 is 426 high. In other words, denitrification can effectively act as negative feedback to nitrogen 427 fixation. The threshold O2 concentration (O 9,gFh ) takes the minimum of the hard-bound O2 428 threshold concentration (O 9,ij!k ) and the NO3 to PO4 ratio, scaled by a parameter k`. The 429 parameters O 9,ij!k and k` are calibrated to give the global denitrification rate of roughly 430   Table 3). It is a challenge for MESMO and other coarse resolution models to 495 simulate narrow dynamical features such as equatorial upwelling and reproduce 496 biogeochemical features with sharp gradients. The spatial pattern of POC export that drives 497 this surface nutrient pattern is similar in the two models ( Figure S2). In the 1D global 498 profile, there is a marked improvement in the subsurface distribution of O2 in MESMO 3 499 over MESMO 2. Whereas the depth of the oxygen minimum was ~300 m in MEMOS 2, it is 500 ~1000 m in both MESMO 3 and the World Ocean Atlas (Figure S3). At 1000 m, the O2 501 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 503 1000 m is too well-oxygenated in MESMO 2. We believe that the improved match in the O2 504 minimum depth would help simulate denitrification in the correct depth range, and there is 505 a modest improvement in the data-model O2 mismatch in terms of RMSE in MESMO 3 over 506 MESMO 2 (Table 3). The deepening of the O2 minimum was achieved largely by increasing 507 the particle sinking speed, which tends to strengthen the biological pump and deplete the 508 surface nutrients. This also helps MESMO 3 preserves MESMO 2's surface Si(OH)4 depletion 509 in much of the world ocean except in the North Pacific and Southern Ocean ( Figure S4). DONsl, and DOFesl. We note only that the surface DOCsl concentration of 58 µmol kg -1 and 517 DOC export production of 1.4 Pg C yr -1 in MESMO 3 are higher than in MESMO 2 (24 µmol 518 kg -1 and 0.4 Pg C yr -1 , respectively). The higher surface concentration is due to the longer tsl 519 in MESMO 3 (Table 2d). The global average of the temperature dependent fDOM in MESMO 520 3 is 0.69, which is slightly higher than the spatially uniform value of 0.67 in MESMO 2. 521 522

Novel features of MESMO 3 523
An important new feature of MESMO 3 is the representation of the primary producers by 524 three PFTs (Figure 2). The eukaryotes are characterized by the highest maximum growth 525 rate and high half-saturation constants. Thus, the eukaryotes are more dominant than the 526 other PFTs in the more eutrophic waters of the equatorial and polar regions (Figure 2a). 527 The cyanobacteria have smaller half-saturation constants and thus are more dominant in 528 the oligotrophic subtropical gyres (Figure 2c). The diazotrophs do not have NO3 limitation 529 but have the lowest maximum growth rate. Thus it is much lower in abundance than the Outside the Southern Ocean, the eukaryotes are primarily limited by Si(OH)4 ( Figure 2b). 535 As far as organic carbon is concerned, we consider the eukaryotes to basically represent 536 diatoms, which are arguably the most important agent of organic carbon export. In this 537 context, the widespread silica limitation for eukaryotes would be consistent with the 538 notions that silica uptake by diatoms should be limited in ~60% of the world surface ocean 539 The global pattern of the mean C:P uptake ratio in the production layer is shown in Figure  558 4. Consistent with observations , the simulated C:P ratio of the 559 phytoplankton community is elevated in the oligotrophic subtropical gyres and low in the 560 eutrophic polar waters (Figure 4a). The community C:P ratio exceeds 200 in the gyres and 561 reaches as low as 40-50 in the Southern Ocean. The community C:P has contributions from taxonomic effects (i.e., the shift in the community composition changes the weighting of 564 each PFT's C:P ratio). Figure 4b shows that the community C:P is high in oligotrophic gyres 565 because cyanobacteria and to a lesser extent diazotrophs dominate the community, and 566 their C:P ratio is high. Conversely, the community C:P is low in the polar waters because the 567 eukaryotes dominate and their C:P ratio is low. For both eukaryotes and cyanobacteria, 568 their C:P is high in oligotrophic subtropical gyres because PO4 is low (Figure 4c and d). This 569 physiological effect is larger in eukaryotes than cyanobacteria because the former has 570 greater sensitivity (i.e., larger sensitivity factor +, ! +:4 , Equation 5 , Table 2b). However, the 571 cyanobacteria PFT's C:P ratio has an additional sensitivity to temperature (i.e., % +:4 ≠ 0) 572 that elevates their C:P in the lower latitudes. We do not show the C:P ratio for diazotrophs 573 because it is very similar to that of cyanobacteria (Figure 4b, d). 574

575
In order to gain more insights into the spatial patterns of the C:P ratio (Figure 4) Table 2b). The C:N-PO4 correlation exists, simply because the nutrients are 585 well correlated. Similarly, because temperature and photosynthetically active radiation 586 (PAR) tend to be correlated via latitude, the stoichiometry has a similar correlation to the 587 two drivers. For example, cyanobacteria C:P has a strong correlation with both 588 temperature and PAR (Figure 5j, 4n), but only the temperature is a real driver. Figure 5  589 indicates which are the dominant drivers of the C:N:P ratio in MESMO 3. For the eukaryote 590 Figure 6 shows the community C:P and C:N ratios plotted against the four environmental 595 drivers. Unlike Figure 5, which reflected the individual PFT's physiological response, Figure  596 6 includes the effect of taxonomy as well. Still, the effects of PO4 and temperature are 597 clearly visible on the community C:P ratio. Both low [PO4] and warmer waters are found in 598 the lower latitudes, so the P frugality and temperature effects are additive. The effect of 599 NO3 on the community C:N ratio is also very clear, but the effect of PAR is not as clear. Thus for the key DOM remineralization model parameters (Table 2d)  the surface DOCt is too high because DOCr is too high, while DOCsl is not unreasonable 672 ( Figure 7). DOCr is too high because there is too much DOCr production (e.g., fDOMr=1% is 673 too large), there is too little DOCr degradation (e.g., one of the DOM decay mechanisms is 674 too slow; Equation 28 and Table 2d), or some combination of both. For example, fDOMr is a 675 key parameter that is not well constrained by observations. Had we used 0.2% instead of 676 1% for fDOMr, the global mean surface DOCt drops to 76 µmol kg -1 (red line, Figure 7), 677 consistent with observations. For achieving a better surface DOCt pattern, we may need a 678 different formulation of fDOM that is, for example, negatively related to nutrient 679 concentrations so that fDOM increases in the oligotrophic subtropical gyres (Roshan and Another lesson from the DOM modeling exercise is that it is important to simulate DOPr 683 reasonably well in order to preserve the favorable results we achieved in MESMO 3 with 684 respect to biological production and the phytoplankton C:N:P ratio. We find that in the 685 experiment LV, the global mean DOPr concentration becomes steady at 0.45 µmol-P kg -1 . In 686 observations, the mean DOCr is about 40 µmol-C kg -1 , and the DOCr:DOPr ratio is estimated 687 to be ~1370:1 (Letscher and Moore, 2015), so DOPr concentration should only be roughly 688 0.03 µmol-P kg -1 . Thus, the simulated DOPr=0.45 µmol-P kg -1 is an order of magnitude too 689 high. Because there is more P in the form of DOPr in LV, the oceanic inventory of PO4 690 declines, causing a nearly 10% drop in export production compared to the standard 691 The modeled global depth-integrated N2 fixation is 101 Tg N yr -1 (Table 3) where nitrate limits eukaryotes and cyanobacteria (Figure 2b, d), can be explained by the 714 healthy growth of diazotrophs, which is not limited by N. In the subtropical and tropical 715 Atlantic and the Indian Ocean, high N2 fixation is driven by elevated C:P and N:P ratio 716 (Figure 4), exemplified by low phosphate availability and warm surface temperature. This 717 spatial pattern agrees with a recent inverse model study (Wang et al., 2019), which showed 718 an elevated N2 fixation rate in subtropical gyres. 719 720 Global water-column denitrification is 101 Tg N yr -1 (Table 3)