Explicit silicate cycling in the Kiel Marine Biogeochemistry Model, version 3 (KMBM3) embedded in the UVic ESCM version 2.9

. We describe and test a new model of biological marine silicate cycling, implemented in the Kiel Marine Biogeo-chemical Model version 3 (KMBM3), embedded in the University of Victoria Earth System Climate Model (UVic ESCM) version 2.9. This new model adds diatoms, which are a key aspect of the biological carbon pump, to an existing ecosystem model. The new model performs well against important ocean biogeochemical indicators and captures the large-scale features of the marine silica cycle. Furthermore it is computationally efﬁcient, allowing both fully-coupled, long-timescale transient 5 simulations, as well as “ofﬂine” transport matrix spinups. We assess the fully-coupled model against modern ocean observations, the historical record since 1960, and a business-as-usual atmospheric CO 2 forcing to the year 2300. The model simulates a global decline in net primary production (NPP) of 1.8% having occurred since the 1960s, with the strongest declines in the tropics, northern mid-latitudes, and Southern Ocean. The simulated global decline in NPP reverses after the year 2100 (forced by the extended RCP 8.5 CO 2 concentration scenario), and NPP returns to 96% of the pre-industrial rate by 2300. This recovery 10 is dominated by increasing primary production in the Southern Ocean, mostly by calcifying phytoplankton. Large increases in calcifying phytoplankton in the Southern Ocean offset a decline in the low latitudes, producing a global net

, and Schmittner et al. (2008) for original references and a complete description of the other equations.

90
Tracer concentrations (C) vary according to: with T including all transport terms (advection, diffusion, and convection), S representing all source and sink terms, and B representing air-sea interface boundary (including virtual) fluxes.

95
Phytoplankton (X representing all types except diazotrophs) biomass source and sink terms are: where growth rate (J), mortality (m), and fast recycling (µ * ) terms are described below, and losses to zooplankton grazing (G) are described in the zooplankton equations. The diazotroph equation is similar except that it does not include a linear mortality loss, only a loss to fast recycling. 100 As in Keller et al. (2012), the maximum possible growth rate of phytoplankton (J max ) is a modified Eppley curve (Eppley, 1972), and is a function of seawater temperature (T ), an e-folding temperature parameter T b , and iron availability (u Fe ).
As in earlier versions of the model, diazotroph growth rate is calculated following general phytoplankton (now LP) and then assigned an additional handicap.
105 Nickelsen et al. (2015) assigned a constant iron half saturation (k Fe ) to diazotrophs but in general phytoplankton this parameter varied as a function of biomass (X) to implictly represent different cell sizes in the model. Because our model contains multiple phytoplankton types we revert to a prescribed but unique k Fe for all phytoplankton types and calculate iron availability as: 110 The maximum potential growth rate is then multiplied by a nutrient availability (u) for nitrate, phosphate, and dissolved silicic acid (the latter for diatoms only) to calculate potential growth under nutrient limitation but replete light, where k N , and k P are fixed half saturation constants unique to each phytoplankton type (and k P is calculated from 16 P to 1 N Redfield stoichiometry). These equations are applied to obtain maximum possible growth rates as a function of temperature and nutrients: u Si = Si 120 Silica uptake uses the (Aumont et al., 2003) scaling of k Si in mol Si m −3 : The potential growth rate under limited light availability (J I ) but replete nutrients is calculated as: where α chl is the initial slope of the photosynthesis versus irradiance (I) curve in chlorophyll units: and θ is a Chl:C ratio (Nickelsen et al., 2015): For simplicity, the same maximum and minimum values of α chl and θ are used for all phytoplankton types.
Light attenuation by coccoliths is included in the calculation of available irradiance at each depth level: where PAR stands for the photosynthetically available radiation, k w , k c , k CaCO3 , and k i are the light attenuation coefficients for water, all phytoplankton, CaCO 3 , and ice,z is the effective vertical coordinate, a i is the fractional sea ice cover, and h i and h s are calculated sea ice and snow cover thickness. Opal generated by diatoms is not explicitly traced and is therefore not included in the underwater light field.
In Equation 2, two loss terms other than predation are considered. Mortality from old age or disease is parameterised using a linear mortality rate (m). Temperature-dependent fast remineralisation is a loss term used to implicitly account for the microbial 145 loop and dissolved organic matter cycling, and is parameterised using a temperature dependency multiplied by a constant (µ * 0 ):

Zooplankton
Changes in zooplankton population (Z) are calculated as the total available food (phytoplankton, zooplankton, and organic 150 detritus) scaled with a growth efficiency coefficient ( ) minus mortality. In addition to old age and disease calculated with a quadratic mortality function (m Z Z 2 ), zooplankton mortality also encompasses losses from higher trophic level predation (G Z ).
Zooplankton grazing (G) follows Keller et al. (2012). Relevant parameters are listed in Table 3. Grazing of each food source is 155 calculated using a Holling II function, where a calculated maximum zooplankton grazing rate (µ max Z ) is reduced by a scaling that is weighted by a food preference (ψ X , where "X" stands for any of the food sources), the total prey population and a half saturation constant for zooplankton ingestion (k z ): The calculated maximum potential grazing rate is a function of a maximum potential grazing rate at 0 • C (µ θ Z ), temperature, 160 and oxygen, where grazing activity is capped when temperatures exceed 20 • C: Grazing is also reduced under hypoxic conditions (r O2 sox ): where O 2 is dissolved oxygen in mmol m −3 . 165 6 https://doi.org/10.5194/gmd-2020-235 Preprint. Discussion started: 9 September 2020 c Author(s) 2020. CC BY 4.0 License.

Organic Detritus
As was introduced in Kvale et al. (2015b), organic carbon detritus sources and sinks are split into "free" (Detr f ree ) and "ballast" (Detr bal ) pools using a fixed ratio (Rbal : tot). Ballast detritus is formed of the CaCO 3 -protected portions of calcifying phytoplankton (CP), and zooplankton. This protected portion does not interact with nutrient pools directly, and instead transfers from the "ballast" to the "free" detrital pool at the rate of CaCO 3 dissolution (λ CaCO3 ): 170 where γ is the food assimilation efficiency, R CaCO3:POC is a fixed production ratio of CaCO 3 to organic detritus, R C:N is a

175
Redfield molar ratio, and w C is the sinking speed of calcite. Free detritus is described as: where µ D is the detrital remineralisation rate and w D is the sinking speed of organic detritus.

Detrital Iron
Iron in detritus follows Nickelsen et al. (2015), but with the additional contribution of calcifiers and diatoms: R Fe:N is a fixed iron to nitrogen ratio that converts total organic carbon detritus sources into iron units. The remineralisation and sinking of detrital iron occur at the same rates as free organic detritus. Two scavenging processes are considered: adsorption of dissolved iron onto organic detritus (Fe orgads ): 190 and particulate CaCO 3 (Fe orgads ca ): calculated as a function of free organic detritus and carbonate particles, respectively, a scavenging rate constant (kFe org and kFe ca ), iron available for scavenging (Fe prime ), and inorganic precipitation of dissolved iron into colloidal material (Fe col ), which is calculated as a linear function independent of organic particle concentration and unchanged from Nickelsen et al. 195 (2015). While calcite is known to be a powerful scavenger of iron and other trace metals (Olsson et al., 2014), the strength of plankton-derived calcite in scavenging dissolved iron is unquantified. A test of this sensitivity is planned at a later stage, so for now the calcite scavenging parameter value is set equal to that of organic detritus.

Calcite
As in Kvale et al. (2015b), the source and sink terms for CaCO 3 include both phytoplankton calcifier and zooplankton sources 200 from grazing and mortality, and losses from dissolution and sinking: where w C is the sinking speed of calcite, and the dissolution specific rate (λ CaCO3 ) is calculated using a parameterisation from Calcite associated with living plankton is calculated separately as:

Opal
The production of biogenic opal is calculated as a function of the diatom grazing and linear mortality loss terms: where R Opal:POC is a production ratio that varies as: This parameterisation was introduced by Aumont et al. (2003) and yields an average surface opal:free detritus export value of around 1 across the Southern Ocean, using a fixed average ratio (R Opal:POC0 ) of 0.5. Production of lithogenic opal occurs mostly on land (Tréguer and De La Rocha, 2013), so its contribution to marine silicate cycling is included simplistically via the dissolved silicate river flux calculation. 215 We present a novel opal dissolution parameterisation based on evidence that opal dissolution rates increase when organic coatings are stripped away by strong bacterial activity (Sarmiento and Gruber, 2006). We approximate an exponential flux function and apply our e-folding temperature parameterisation to represent microbially enhanced dissolution: 8 https://doi.org/10.5194/gmd-2020-235 Preprint. Discussion started: 9 September 2020 c Author(s) 2020. CC BY 4.0 License.
We find this parameterisation offers simplicity and consistency with other temperature-dependent rates in our model, as well as adding mechanistic realism, and improving model fit to World Ocean Atlas silica distributions (relative to other parameterisations that we tested, e.g. the temperature-dependent parameterisation of Gnanadesikan (1999) or the temperature and oxygen-dependent parameterisation of Enright et al. (2014)). The Gnanadesikan (1999) parameteristion yields lower dissolution rates at low temperatures than the Enright et al. (2014) parameterisation, which is similarly formulated but which includes an additional oxygen scaling. The Enright et al. (2014) oxygen scaling is not justified in their model description, but it has the 225 effect of increasing Si dissolution rates in the deep ocean (exacerbating the overestimation of Si dissolution in this region by the Gnanadesikan (1999) scaling described in Ridgwell et al. (2002)) and decreasing Si dissolution rates (to a lesser extent) in the near-surface. Our temperature scaling has the effect of raising dissolution rates at the surface.

Particle Sinking
Detritus, calcite, and iron particles are exported from the surface with a sinking speed (w) that increases linearly (wdd, wdc) 230 with depth (z) for calcite: and free detritus and associated iron: The initial surface sinking speeds of particles are assigned different values to represent the denser structure of CaCO 3 relative 235 to that of organic detritus. Ballasted detritus sinks at the CaCO 3 speed, but once it enters the free pool it uses the organic detrital sinking speed and remineralisation rate. Any organic detritus reaching the sediments is dissolved back in to the water column to ensure conservation of carbon and phosphate. Calcite particles that reach the seafloor enter the sediment model (Kvale et al., 2015b) if it is active, though we do not use it here (in this case the particles dissolve). Iron detritus reaching the sediments is lost from the ocean, unless bottom water oxygen falls below 5 mmol m −3 , whereupon detrital iron reaching the bottom is 240 dissolved back into the water column (Nickelsen et al., 2015). By definition, iron is not mass-conserving in the model, but is formulated with open boundary conditions (atmosphere and sediments) -hence it is not just the oxygen threshold that can cause a loss or gain of marine iron. In the model the oxygen threshold is only met along coastal boundaries in the north Pacific under modern conditions. A sedimentary release of dissolved iron (Fe sed ) is prescribed as in Nickelsen et al. (2015):

245
where R Fe:Psed is a ratio of iron released from sediment in proportion to the flux of phosphorus in the organic detritus (F POP , which includes free and ballast detritus) reaching the bottom.
Opal reaching the seafloor is either returned to the silicic acid tracer or considered lost to the sediments (bsi) according to criteria laid out in Sarmiento and Gruber (2006):

250
Opal flux, Exp(Opal), is calculated as opal production (Pr(Opal)) less the instantaneous dissolution term (Di(Opal)). The model conserves silica when river fluxes are allowed to compensate for the net change due to the ocean sedimentary sink, hydrothermal input, and atmospheric deposition.
The sub-grid bathymetric scaling that was introduced for the iron model in Nickelsen et al. (2015) is extended here to apply to all particle fluxes reaching the bottom ocean grid. This scaling feature calculates a sedimentary exchange factor based on the 255 proportion of each grid cell that falls outside of the model grid's depth according to a high-resolution bathymetric dataset. This scaling is used to account for high-relief bathymetric features, such as ridges and troughs, in the sediment transfer functions.
We do not use the sediment model in this manuscript.

Dissolved Inorganic Tracers
Ocean nutrient sources and sinks follow: where R P:N and R O:N are Redfield molar ratios and u N is a Michaelis-Menten uptake rate. In suboxic water (less than 5 mmol m −3 ), oxygen consumption is replaced by denitrification (r Dissolved iron includes sources and sinks of particulate iron mentioned above, as well as prescribed dust deposition as in 270 Nickelsen et al. (2015) and hydrothermal iron (Fe hydr ) (Muglia et al., 2017) boundary terms: The Nickelsen et al. (2015) model applied the iron dust flux to the surface ocean after the biological routine and before the 275 mixing routine. This resulted in low model sensitivity to iron dust inputs, because the dust was mixed away before the biology had a chance to access it. In this version, the iron dust flux is added prior to the biological routine (which operates on a shorter timestep than ocean mixing), and results in a greater biological sensitivity to iron dust flux.
DIC and alkalinity tracer sources and sinks are a function of sources and sinks of prognostic CaCO 3 (Kvale et al., 2015b).
Dissolved silicic acid tracer sources, sinks, and boundary terms follow: where deposition from rivers (Si riv ) is used as a budget balancing term to compensate for any remainder in the other external sources and sinks: windborne dust (Si dust ), hydrothermal silicate (Si hydr ), and loss to the sediments (Si sed ). Dust deposition from the atmosphere is prescribed using an interpolated monthly pre-industrial dust flux derived from the NCAR's Community Climate System Model (Mahowald et al., 2006). The silica content of the dust is derived from maps (Zhang et al., 2015), with a 290 global average dust solubility of 3% assumed to produce annual bioavailable fluxes in agreement with estimates (Sarmiento and Gruber, 2006;Tréguer and De La Rocha, 2013). Observations demonstrate wide spatial variability in the solubility parameter (Tréguer and De La Rocha, 2013) so our single value represents a simplification. Silicate from hydrothermal sources are prescribed using a static mask scaled from the hydrothermal iron mask (Muglia et al., 2017) using a Fe:Si ratio to obtain the estimated annual total contribution from hydrothermal sources in Sarmiento and Gruber (2006).

Model Assessment in a Modern Climate
The model is spun up to equilibrium (greater than 15,000 years) at year 1765 boundary conditions, then forced with historical data for comparison to observational datasets. Forcing includes historical atmospheric CO 2 concentrations, agricultural land cover, volcanic radiative forcing, sulphate aerosol and CFC concentrations, changes in land ice and solar forcing (Machida et al., 1995;Battle et al., 1996;Etheridge et al., 1996Etheridge et al., , 1998Flückiger et al., 1999Flückiger et al., , 2004Ferretti et al., 2005;Meure et al., 300 2006). Solar insolation at the top of the atmosphere, wind stress, and wind fields vary seasonally (Kalnay et al., 1996), and wind fields are geostrophically adjusted to air temperature anomalies (Weaver et al., 2001). Table 5 lists key biogeochemical properties diagnosed by the model for a modern climate, as well as corresponding properties Phytoplankton biomass is compared to the MAREDAT datasets (Leblanc et al., 2012;Luo et al., 2012;O'Brien, 2012) in North Atlantic), though the model appears to over-estimate the North Atlantic biomass. Like the globally integrated biomass estimate, calcifiers (CP) are universally over-estimated compared to O'Brien (2012). As explained above, this phytoplankton type is meant to represent a variety of moderate growth, moderate nutrient affinity phytoplankton types, including those that calcify. Therefore it is not unexpected to have an over-estimate of the biomass. Diazotroph biomass is primarily concentrated in the tropics, which agrees spatially but not in magnitude with the limited data of Luo et al. (2012).

355
Interior ocean particle flux performance is encouraging. Figure  Seasonal succession in phytoplankton types follows the general progression of zonal maxima in diatoms preceeding calcifiers by a few weeks ( Figure 13). This succession is due to the higher nutrient requirements, and faster growth rates, of the diatoms, which are able to take advantage of winter mixing early in the growing season. Once surface nutrient concentrations start to decline, the calcifying phytoplankton become relatively more successful. This pattern is most pronounced in the Southern

405
Ocean. In the low latitudes, diazotroph biomass peaks earlier in the growing season than the low-latitude, non-calcifying phytoplankton (LP). Diazotrophs have a growth handicap with respect to LP, but their ability to fix nitrate gives them an advantage in more stratified summer conditions. This fixed nitrate is then used by the LP in the winter months.
In addition to the historical model forcings described earlier, from year 2005 to 2300 the simulations are forced using increasing 410 CO 2 and non-CO 2 greenhouse gas concentrations, projected changes to the fraction of the land surface devoted to agricultural uses (calculated to year 2100 by Hurtt et al. (2011), and then held constant after), and changes in the direct effect of sulphate aerosols following "business-as-usual" RCP scenario 8.5 (RCP8.5, Riahi et al., 2007;Meinshausen et al., 2011). The wind fields continue to be geostrophically adjusted to air temperature anomalies. where non-calcifying phytoplankton have increased their biomass. Warming, and a lengthening growing season (but increased stratification; see Arctic temperature trend in Figure 16) in the Arctic benefits calcifiers. Expansion of coccolithophores into 425 the Arctic has been observed over recent decades (Neukermans et al., 2018). The Southern Ocean is also simulated to have experienced increasingly favourable conditions for calcifiers, with diatom biomass declining as calcifiers increase.
Our model simulates a slower global decline in productivity and biomass than the scarce, and controversial, satellite and in-situ chlorophyll record reconstructed by Boyce et al. (2010), who calculated a 1% per year decline in chlorophyll. Our model also does not simulate a large decline prior to the 1950s ( Figure 17). Wernand et al. (2013) found no historical trend in diatoms. Our results are more similar to those of Rousseaux and Gregg (2015), who combined a model with satellite data to reconstruct ocean surface changes from 1998-2012. They simulated a global decline in diatoms over this period, which they attributed to an increase in nutrient limitation and photosynthetically available radiation, which favours other phytoplankton types. However, the Southern Ocean in their model showed no clear trend in phytoplankton community composition between 445 1998-2012 (Rousseaux and Gregg, 2015). The model we present here is fully competition-driven. Variable particle sinking and remineralisation rates, and explicit CaCO 3 ballasting, further differentiate our model. All of these factors may contribute to the different behaviour reported here.
Trends in biomass and NPP affect deep particle fluxes. Low latitude declines in primary production result in less deep ocean particle export in the western Pacific and Indian Ocean basins, while deep POC and PIC export is simulated to have slightly is temperature-dependent, therefore regional warming is almost certainly contributing to the reduction of deep ocean opal flux, but the loss of diatom biomass is the major driver of this trend in the Southern Ocean. Global losses of particle export across 2 km depth between 1964 and 2014 are calculated at 2.6% (POC), 2.0% (PIC), and 6.1% (opal).
Unfortunately there is no comparable historical reconstruction which have been implicated in recent Pacific interior nutrient trends (Stramma et al., 2020). This may be why our model does a poor job reproducing observations of nutrient trends observed in the central Pacific (Stramma et al., 2020). Also, increases in nitrate are found in the observed record, which our model does not simulate. This may be because our model does not include anthropogenic sources of nitrate (summarized by Stramma et al., 2020).
Surface nutrient concentrations are simulated to have increased slightly or show no trend in the North Atlantic (Figure 19),  (Iida et al., 2013). The Southern Ocean is simulated to have experienced sub-surface increases in phosphate in the regions also experiencing increases in silicic acid and iron; this is due to the less efficient export of particles by diatoms (locally increasing over this time), relative to calcifiers.
Taken together, these results suggest low latitude declines in NPP and export production have produced a decline in down-500 stream nutrient accumulation in northern hemisphere intermediate water masses that are particularly apparent in the North Pacific and north Atlantic regions. In the Southern Ocean, declines in NPP and export production upstream have introduced excess nutrients to the basin, raising the nutrient concentrations in Antarctic Bottom and Intermediate water masses despite declining Southern Ocean NPP. These production and export effects have been exacerbated by physical changes in the circulation; Figure 16 shows a decline in ideal age (greater than 5 years in the zonal mean) in Pacific Intermediate water and the 505 Indian Ocean, which leads to less particle remineralisation (and hence, lower nutrient concentrations) there. Increasing water mass age in Southern Ocean-sourced intermediate and deep water masses (upwards of 30 years in the deep Pacific), likewise has the effect of producing more complete particle remineralisation, resulting in higher nutrient concentrations.
Oxygen is simulated to have declined in all ocean basins, with the exception of the sub-surface low latitudes. It has been previously estimated that the global ocean lost more than 2% of its oxygen since the 1960s ). Our 510 model simulates a 0.6% decline in total oxygen content from 1964-2014, which is an underestimate resulting from physical biases in our model (i.e., deficiencies in simulating low latitude ventilation), though biogeochemical deficiencies might also be contributing (Oschlies et al., 2017). In the deep ocean, the ageing of water masses (and associated more complete particle remineralisation) contributes to the simulated decline in oxygen. In the upper ocean, warming has reduced oxygen solubility, lowering near-surface concentrations. ). In this latest model version, this increase in NPP is much smaller (less than 6 Pg C y −1 by 2300, compared to 11-13 Pg C y −1 in previous model versions). Again, this reduced sensitivity in the low latitudes is due to the reduction of low-latitude ballast-forming calcifiers in the new version. The response of calcifiers outside of the low latitudes to ocean changes is also changed from previous versions, in that our model now simulates increases in biomass in the Southern Ocean as well as the Arctic. Whether these differences are competition effects, or due to other model changes, such as the correction of light atten- Water formation has declined. The net effect of these physical changes is an overall decline in low-latitude productivity (driven by increased nutrient limitation) but a strong increase in Southern Ocean productivity, with a faster biogeochemical connection between the surface ocean south of the Polar Front, and the abyssal basins (see the improved ventilation in the radiocarbon plots). At the poles, fast-growing, nutrient-demanding phytoplankton types (DT, CP) thrive, while in the lower latitudes it is 555 the more efficient nutrient consumers (DZ and LP) who benefit.
Long-term particle export trends generally follow the historical trend, but with increasing magnitude (Figures 17 and 23).
Globally integrated particle fluxes decline, and remain suppressed with respect to pre-industrial rates, for POC and opal. PIC surface export rates change very little and deep export rates increase with climate forcing as a response to increasing surface calcifier POC export fluxes (e.g., Kvale et al., 2015a). The production of both opal and PIC are scaled against their respective 560 plankton types' POC production, so it is expected that PIC and opal fluxes follow the POC export production trend. Just as with NPP, the POC export production decline in our new model is larger than in previous versions (about 1.0 Pg C y −1 by 2100, rather than 0.8 PgC y −1 ).
POC and PIC fluxes increase where calcifying phytoplankton biomass also increases; south of 40 S, in the eastern equatorial Pacific upwelling zone, and along the Kuroshio Current into the North Pacific ( Figure 23). Strong decreases in calcifiers, and 565 associated deep carbon export, occur in the Indian Ocean and North Atlantic. Opal export declines by more than 100 mmol Si m −2 y −1 both south of 40 S and north of 40 N (with only small changes in the Arctic). However, opal export (and diatoms) increase south of 60 S, where increased nutrients (Figure 24), particularly iron, and a short growing season favours diatoms over calcifiers.
The historical trend in carbon and nutrients is similarly extended, with continuing increases in DIC in the upper ocean (as 570 atmospheric CO 2 continues to enter the ocean), declines in low latitude upper ocean nutrients phosphate, nitrate, and silicic acid (due to decreasing resupply from the deep ocean), increases in the deep ocean in the same nutrients, and widespread declines in oxygen ( Figure 24). Oxygen declines along the Southern Ocean-abyssal global ocean pathway due to both warming and increasing particle remineralisation, which is also responsible for the increasing nutrient concentrations in the deep ocean.
Decreasing phosphate and nitrate concentrations in the sub-surface tropical ocean basins are a product of declining particle 575 remineralisation there, brought about by both warming, which shoals remineralisation and increases respiration rates, and a shift to less efficiently exporting phytoplankton (LP).
The striking trend in dissolved iron that emerges in these future projections of strongly increasing (more than 80 nmol m −3 in the zonal mean) concentrations in the upper ocean was previously described by Nickelsen et al. (2015). They attributed dissolved iron are also the regions experiencing both strong declines in NPP (and hence, lower iron uptake) and strong declines in particle export (and hence, less particulate iron scavenging and removal). The loss of calcifiers in the Indian Ocean and central Pacific particularly increases iron concentrations there, because of the dual effect of reduced POC and PIC scavenging of iron.
Our future simulation results broadly agree with other long-term simulations in the sustained, and significant, increase 585 in Southern Ocean primary production that couples with a reorganisation of deep ocean circulation to produce a long term "nutrient trapping" effect in Southern Ocean-sourced interior water masses (Moore et al., 2018;Kvale et al., 2019). Nearsurface increases in iron, and decreases in nitrate, phosphate, and silicic acid, have also been observed to 2100 in a comparison of 9 other earth system models by Fu et al. (2016). These same models also simulate weak to strong increases in diatoms in the Southern Ocean to 2100, though in most, if not all, of them, diatoms are the most efficient exporters of carbon and nutrients

590
(unlike in our model). Phytoplankton community composition and export formulation was discussed by Fu et al. (2016) to be of critical importance in determining trends in NPP, nutrients, and particle export over the coming century, thus a diversity of model formulations benefits our understanding of how the global ocean ecosystem might change in the future.

Conclusions
Our manuscript describes a new model of the marine silicate cycle, evaluates its performance against previous KMBM versions 595 as well as key biogeochemical data derived from observations of the ocean, and compares long-term ecosystem projections to similar models available in the literature. We find our new model shows general improvement in the representation of nutrients and particle fluxes and is mechanistically more realistic, with the added complexity of iron, calcite, and silicate merged into a single model code.
Simulations using our new model suggest diatoms have been, and will continue to be, the losers as the earth system warms.

600
Their high nutrient requirements prove a disadvantage as the upper ocean stratifies, and small gains in productivity provided by sea ice retreat cannot compensate for the fact that their southern bound is ultimately limited geographically. Calcifying phytoplankton with more moderate nutrient requirements are the big winners across the high latitudes, while in the tropics slow-growing, less nutrient-hungry phytoplankton are projected to thrive. From a deep ocean carbon sequestration perspective, the loss of diatom export production is of transient importance, as the calcifying phytoplankton increase their role in carbon 605 export, efficiently sinking organic carbon as well as carbonate.
Our simulations also reveal the past may not accurately portray future trends, as evidenced by simulated historical declines in NPP in the Southern Ocean that reverse as conditions become more favourable for calcifiers. Significant and rapid increases in dissolved iron in the low latitude tropical ocean is another potential biogeochemical "surprise", still to come, if anthropogenic emissions of carbon follow the present trajectory.

610
Several novel aspects of our model, including iron scavenging by calcite, silicate source and sink strengths, and different zooplankton grazing preferences are slated for further study. The impact of variable stoichiometry is another important potential aspect of biogeochemical modelling that is not explored here. More complete parameter assessment is planned in the context of offline parameter optimisation and model calibration experiments (e.g., Kriest, 2017;Kriest et al., 2017;Yao et al., 2019) in the future, as is merging this new biogeochemical model into the latest UVic ESCM version 2.10 (Mengis et al., 2020). We 615 look forward to further refinements, and the many applications of this model to come.     29 https://doi.org/10.5194/gmd-2020-235 Preprint. Discussion started: 9 September 2020 c Author(s) 2020. CC BY 4.0 License.