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 Biogeochemical 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 component of the biological carbon pump, to an existing ecosystem model. This new model combines previously published parameterizations of a diatom functional type, opal production and export with a novel, temperature-dependent dissolution scheme. Modelled steady-state biogeochemical rates, carbon and nutrient distributions are similar to those found in previous model versions. The new model performs well against independent ocean biogeochemical indicators and captures the large-scale features of the marine silica cycle to a degree comparable to similar Earth system models. Furthermore, it is computationally efficient, allowing both fully coupled, long-timescale transient simulations and “offline” transport matrix spinups. We assess the fully coupled model against modern ocean observations, the historical record starting from 1960 and a business-as-usual atmospheric CO₂ forcing to the year 2300. The model simulates a global decline in net primary production (NPP) of 1.4 % having occurred since the 1960s, with the strongest declines in the tropics, northern midlatitudes and Southern Ocean. The simulated global decline in NPP reverses after the year 2100 (forced by the extended RCP8.5 CO₂ concentration scenario), and NPP returns to 98 % of the pre-industrial rate by 2300. This recovery 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 calcite export in 2300 that varies only slightly from pre-industrial rates. Diatom distribution moves southward in our simulations, following the receding Antarctic ice front, but diatoms are outcompeted by calcifiers across most of their pre-industrial Southern Ocean habitat. Global opal export production thus drops to 75 % of its pre-industrial value by 2300. Model nutrients such as phosphate, silicate and nitrate build up along the Southern Ocean particle export pathway, but dissolved iron (for which ocean sources are held constant) increases in the upper ocean. This different behaviour of iron is attributed to a reduction of low-latitude NPP (and consequently, a reduction in both uptake and export and particle, including calcite scavenging), an increase in seawater temperatures (raising the solubility of particulate iron) and stratification that “traps” the iron near the surface. These results are meant to serve as a baseline for sensitivity assessments to be undertaken with this model in the future.

1 Introduction

It has become apparent in recent decades that the representation of elemental cycles of nutritive elements (N, P, Si, Fe) is important in order to simulate critical biogeochemical feedbacks (Heinze et al., 2019). Here, we describe an updated representation of ocean biogeochemistry (the Kiel Marine Biogeochemical Model, version 3; KMBM3) embedded in an Earth system model, the University of Victoria Earth System Climate Model version 2.9, second-level updates (UVic ESCM; Weaver et al., 2001; Eby et al., 2009), to which we have added a silicate cycle.

Silicate (considered here as any anion that combines silicon, Si, and between two and four oxygen molecules) is important for simulating ocean biogeochemistry for several reasons. It regulates the growth of diatoms, which are phytoplankton responsible for 40 % of the particulate carbon export in the modern ocean (Jin et al., 2006). Fluctuations in the relative abundance of diatoms might therefore affect air–sea CO₂ exchange, with a globally significant influence on atmospheric CO₂ concentrations over millennial and longer timescales (e.g. Matsumoto et al., 2002; Renaudie, 2016). High-latitude diatoms with low N : P ratios have been demonstrated to exert global control on the nitrogen budget (Weber and Deutsch, 2012) and are speculated to exert seasonal control on carbonate chemistry through competition with calcifiers (Merico et al., 2006). Silicic acid (Si(OH)₄ and SiO(OH)${}_{\mathrm{3}}^{-}$) distributions are tightly controlled by diatom primary production (biogenic silica or opal; SiO₂) and subsequent dissolution (e.g. only about 34 % global average ocean silicic acid is preformed, as opposed to biologically regenerated; Holzer et al., 2014). First attempts to constrain the marine silicate cycle demonstrate that representations of the silica cycle are well suited to parameter optimization (Holzer et al., 2014; Pasquier and Holzer, 2017). Due to diatoms' important role in biological carbon export, evaluating simulated silicate and diatoms presents a unique opportunity to constrain ocean carbon cycling.

In the following sections, we describe our new model of explicit diatoms and silicate cycling in the ocean. As far as we are aware, this model is unique among ecosystem models embedded within Earth system models in that it combines established parameterizations of a diatom functional type (i.e. fast growth rates and inefficient nutrient utilization, Tréguer et al., 2018) with a novel temperature-dependent opal dissolution parameterization and bottom-up and top-down competition for resources with explicit, fully independent calcifying phytoplankton and nitrogen fixing diazotrophs. These assumptions bias KMBM3 diatoms towards large species in the well-mixed high latitudes and early seasonal succession. Modelled silicate boundary conditions include hydrothermal silicate inputs and a benthic transfer function to approximate sediment sequestration. A ballast model applied to the calcifiers (Kvale et al., 2015b) further differentiates KMBM3, as does the ability to operate in an “offline” parameter optimization framework (Kvale et al., 2017; Yao et al., 2019).

2 Model description

2.1 Physics

UVic ESCM version 2.9 is a coarse-resolution (1.8^∘ × 3.6^∘ × 19 ocean depth layers) ocean–atmosphere–biosphere–cryosphere–geosphere model. Model structure and physics are described in Weaver et al. (2001), Meissner et al. (2003) and Eby et al. (2009). We use the second-level updates of version 2.9. Additionally, an ideal age tracer (Koeve et al., 2015) is merged into this version of the model to allow for explicit tracing of water mass age.

2.2 Biogeochemistry

Development of the Kiel Marine Biogeochemistry Model has progressed in recent years. Beginning with a base composition of nitrate, phosphate, dissolved inorganic carbon (DIC) and alkalinity as well as general phytoplankton, diazotrophs and one zooplankton functional type (Schmittner et al., 2005, 2008), Keller et al. (2012) updated zooplankton grazing and included a seasonally cycling iron mask (KMBM1). To this version, prognostic iron tracers were implemented by Nickelsen et al. (2015) (KMBM2), and prognostic CaCO₃ and calcifying phytoplankton were added by Kvale et al. (2015b). Somes et al. (2013) included benthic denitrification and Muglia et al. (2017) included hydrothermal iron in MOBI, a separate branch of the UVic ESCM biogeochemistry. Kvale and Meissner (2017) examined the sensitivity of primary production in the model to changes in the light attenuation parameter following a correction in the calculation of light availability (Partanen et al., 2016). To the Keller et al. (2012) version, Kvale et al. (2017) added offline transport matrix method (TMM) capability, allowing for rapid matrix-based computation of model mean ocean tracer states and coupling to alternative physical models (Khatiwala et al., 2005; Khatiwala, 2007). The current model description unites divergent code, assimilating the Nickelsen et al. (2015) and Kvale et al. (2015b, 2017) models as well as the benthic denitrification code of Somes et al. (2013) and hydrothermal iron of Muglia et al. (2017) and expands the code to include explicit diatoms and a silicic acid tracer.

With up to 17 model tracers, KMBM3 biogeochemistry is now fairly complex. An abbreviated description is provided in the following sections, with only the relevant model developments described here. All simulations use the fully coupled UVic ESCM framework.

We introduce a new phytoplankton functional type to which have been assigned key physiological attributes of diatoms. This new model version therefore contains phytoplankton functional types “slow-growing low-latitude phytoplankton (LP)”, “slow-growing phytoplankton which can fix nitrogen when necessary, including diazotrophs (DZ)”, “moderate-growth phytoplankton, including calcifiers (CP)” and “fast-growing diatoms (DT)”. Zooplankton contribute to CaCO₃ production. Prognostic tracers include nitrate, phosphate, oxygen, DIC, alkalinity, iron, CaCO₃, silicic acid and particulate forms of iron and organic detritus. The new model schematic is shown in Fig. 1, and state variables are given in Table 1. Diatoms compete with the other functional types for light and nutrients. The diatoms' growth is also limited by dissolved silicic acid availability; they implicitly produce opal that instantly remineralizes back into dissolved silicic acid throughout the water column. To ensure absolute model conservation of silica, any opal that is lost due to implicit burial at the seafloor is replaced by external sources (prescribed atmospheric dust deposition, sediment release via a benthic transfer function, prescribed hydrothermal vent release and river fluxes that can be set to compensate for any remainder; see below for details). The phytoplankton and detritus production and remineralization are linked to nutrients through fixed Redfield stoichiometry using a base unit of mmol nitrogen m⁻³.

Figure 1 The latest biogeochemical model structure for KMBM3. Previously unpublished features are shown in orange.

Table 1 KMBM3 state variables.

In the following model description, notation will generally follow the symbols used in Kvale et al. (2015b), with the abbreviations LP, CP, DT and DZ representing the phytoplankton functional types, and “Z” representing zooplankton when a distinction is necessary. The most important biogeochemical model parameters are listed in Tables 2 to 5. The model description here covers only the most relevant equations and equations that have changed in this newest version; please see Kvale et al. (2015b), Nickelsen et al. (2015), Keller et al. (2012) and Schmittner et al. (2005, 2008) for original references and a complete description of the other equations.

Table 2 Miscellaneous KMBM3 parameters.

2.2.1 Sources and sinks

Tracer concentrations (C) vary according to

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \left[C\right]}{\partial t}}=T+S+B,\end{array}$$

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 evaporation–precipitation correction) fluxes and hydrothermal inputs, where relevant.

2.2.2 Phytoplankton

Phytoplankton (X representing all types except diazotrophs) biomass (in units of mmol nitrogen m⁻³) source and sink terms are

$$\begin{array}{}\text{(2)}& S\left(X\right)=J_{X}X-G_X-{\mathit{\mu }}_X^{*}X-m_{X}X,\end{array}$$

where growth (J), mortality (m) and fast recycling (μ^*) rates 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.

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) in degrees Celsius, an e-folding temperature parameter T_b and iron limitation (u_Fe) modifying a growth parameter (a). This parameterization assumes sufficient iron is required for the utilization of other nutrients (Galbraith et al., 2010; Keller et al., 2012; Nickelsen et al., 2015).

$$\begin{array}{}\text{(3)}& J_{\max }=a{\exp }^{\frac{T}{T_b}}{u}_{\mathrm{Fe}}\end{array}$$

As in earlier versions of the model, diazotroph growth rate is calculated following ordinary phytoplankton (now LP) and then multiplied with an additional handicap.

Nickelsen et al. (2015) assigned a constant iron half saturation (k_Fe) to diazotrophs, but in ordinary phytoplankton this parameter varied as a function of biomass (X) to implicitly represent different cell sizes in the model. Because KMBM3 contains multiple phytoplankton functional types, we revert to a prescribed but unique k_Fe for all functional types and calculate iron limitation as

$$\begin{array}{}\text{(4)}& u_{\mathrm{Fe}}={\displaystyle \frac{\left[\mathrm{Fe}\right]}{k_{\mathrm{Fe}}+\left[\mathrm{Fe}\right]}}.\end{array}$$

The maximum potential growth rate is then multiplied by the minimum nutrient limitation term (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 functional 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 following Liebig's law of the minimum:

$$\begin{array}{}\text{(5)}& {\displaystyle }{u}_{{\mathrm{NO}}_{\mathrm{3}}^{-}}={\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]}{k_{\mathrm{N}}+\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]}}\text{(6)}& {\displaystyle }{u}_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}}={\displaystyle \frac{\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]}{k_{\mathrm{P}}+\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]}}\text{(7)}& {\displaystyle }{u}_{\mathrm{Si}}={\displaystyle \frac{\left[\mathrm{Si}\right]}{k_{\mathrm{Si}}+\left[\mathrm{Si}\right]}}.\end{array}$$

Silica limitation uses the empirical fit to experimental data scaling of k_Si in mmol Si m⁻³ from Aumont et al. (2003):

$$\begin{array}{}\text{(8)}& k_{\mathrm{Si}}=\mathrm{0.8}\mathrm{mmol}\mathrm{Si}{\mathrm{m}}^{-\mathrm{3}}+\mathrm{7.2}\mathrm{mmol}\mathrm{Si}{\mathrm{m}}^{-\mathrm{3}}{\displaystyle \frac{\left[\mathrm{Si}\right]}{k_{\mathrm{Si}}^{*}+\left[\mathrm{Si}\right]}},\end{array}$$

with a $k_{\mathrm{Si}}^{*}$ value adopted from Aumont et al. (2003) of 30 mmol Si m⁻³.

The potential growth rate under limited light availability (J_I) but replete nutrients is calculated as

$$\begin{array}{}\text{(9)}& J_I={\displaystyle \frac{J_{\max }{\mathit{\alpha }}^{\mathrm{chl}}\mathit{\theta }I}{(J_{\max }^{\mathrm{2}}+({\mathit{\alpha }}^{\mathrm{chl}}\mathit{\theta }I{)}^{\mathrm{2}}{)}^{\frac{\mathrm{1}}{\mathrm{2}}}}}.\end{array}$$

Nickelsen et al. (2015) introduced this parameterization to the model and discussed it at length; interested readers are recommended to read their Sect. 2.3.2. In the above equation, α^chl is the initial slope of the photosynthesis versus irradiance (I) curve in chlorophyll units:

$$\begin{array}{}\text{(10)}& {\mathit{\alpha }}^{\mathrm{chl}}={\mathit{\alpha }}_{\min }^{\mathrm{chl}}+({\mathit{\alpha }}_{\max }^{\mathrm{chl}}-{\mathit{\alpha }}_{\min }^{\mathrm{chl}})u_{\mathrm{Fe}},\end{array}$$

and θ is a Chl : C ratio (noting, chlorophyll is not a prognostic variable, Nickelsen et al., 2015):

$$\begin{array}{}\text{(11)}& \mathit{\theta }={\mathit{\theta }}_{\min }+({\mathit{\theta }}_{\max }-{\mathit{\theta }}_{\min })u_{\mathrm{Fe}}.\end{array}$$

The approximation of θ is simplistic and neglects other factors, e.g. irradiance, which can affect the ratio. For simplicity, the same maximum and minimum values of α^chl and θ are used for all functional types.

Light attenuation by algae biomass and by coccoliths is included in the calculation of available irradiance at each depth level:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}I& =I_{z=\mathrm{0}}\mathrm{PAR}{e}^{-k_{\mathrm{w}}\stackrel{\mathrm{̃}}{z}-k_{\mathrm{c}}{\int }_{\mathrm{0}}^{\stackrel{\mathrm{̃}}{z}}(\mathrm{LP}+\mathrm{CP}+\mathrm{DZ}+\mathrm{DT})\mathrm{d}z-k_{{\mathrm{CaCO}}_{\mathrm{3}}}{\int }_{\mathrm{0}}^{\stackrel{\mathrm{̃}}{z}}\left[{\mathrm{CaCO}}_{\mathrm{3}}\right]\mathrm{d}z}\\ & \cdot (\mathrm{1}+a_{\mathrm{i}}(e^{-k_{\mathrm{i}}(h_{\mathrm{i}}+h_{\mathrm{s}})}-\mathrm{1}\left)\right),\end{array}\end{array}$$

where PAR stands for the photosynthetically available radiation, k_w, k_c, $k_{{\mathrm{CaCO}}_{\mathrm{3}}}$ and k_i are the light attenuation coefficients for water, all phytoplankton, CaCO₃ and ice, $\stackrel{\mathrm{̃}}{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.

The actual growth rate (J_X) of phytoplankton is taken to be the minimum of the growth functions described above:

$$\begin{array}{}\text{(13)}& {\displaystyle }{J}_{\mathrm{LP}}=\min (J_I,J_{\max }{u}_{{\mathrm{NO}}_{\mathrm{3}}^{-}},J_{\max }{u}_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}}),\text{(14)}& {\displaystyle }{J}_{\mathrm{CP}}=\min (J_I,J_{\max }{u}_{{\mathrm{NO}}_{\mathrm{3}}^{-}},J_{\max }{u}_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}}),\text{(15)}& {\displaystyle }{J}_{\mathrm{DT}}=\min (J_I,J_{\max }{u}_{{\mathrm{NO}}_{\mathrm{3}}^{-}},J_{\max }{u}_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}},J_{\max }{u}_{\mathrm{Si}}),\text{(16)}& {\displaystyle }{J}_{\mathrm{DZ}}=\min (J_I,J_{\max }{u}_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}}).\end{array}$$

In Eq. (2), two loss terms other than predation are considered. Non-grazing mortality is parameterized using a linear mortality rate (m). Temperature-dependent fast remineralization is a loss term used to implicitly account for the microbial loop and dissolved organic matter cycling, and is parameterized using a temperature dependency multiplied by a constant (${\mathit{\mu }}_{\mathrm{0}}^{*}$) (Schmittner et al., 2008):

$$\begin{array}{}\text{(17)}& {\mathit{\mu }}^{*}={\mathit{\mu }}_{\mathrm{0}}^{*}{\exp }^{\frac{T}{T_b}}.\end{array}$$

With this formulation, increasing seawater temperature increases the return of nutrients to the upper ocean.

2.2.3 Zooplankton

Changes in zooplankton population (Z) are calculated as the total ingested food (phytoplankton, zooplankton and organic detritus) scaled with a growth efficiency coefficient (ϖ) minus mortality. In addition to mortality from higher trophic level predation calculated with a quadratic mortality function (m_ZZ²), intra-guild predation is represented with a separate term (G_Z) (Keller et al., 2012; Kvale et al., 2015b; Nickelsen et al., 2015).

$$\begin{array}{}\text{(18)}& \begin{array}{rl}S\left(Z\right)& =\mathit{\varpi }(G_{\mathrm{LP}}+G_{\mathrm{CP}}+G_{\mathrm{DZ}}+G_{\mathrm{DT}}+G_{{\mathrm{Detr}}_{\mathrm{tot}}}+G_{\mathrm{Z}})\\ & -m_{\mathrm{Z}}{Z}^{\mathrm{2}}-G_{\mathrm{Z}}.\end{array}\end{array}$$

Zooplankton grazing (G) follows Keller et al. (2012). Relevant parameters are listed in Table 4. Grazing of each food source is calculated using a Holling II function, where a calculated maximum zooplankton grazing rate (${\mathit{\mu }}_Z^{\max }$) is reduced by a scaling that is weighted by a relative food preference (ψ_X, where “X” stands for any of the food sources and the sum of all preferences must equal 1), the total prey population and a half-saturation constant for zooplankton ingestion (k_z):

$$\begin{array}{}\text{(19)}& G_X={\mathit{\mu }}_Z^{\max }Z{\displaystyle \frac{{\mathit{\psi }}_{X}X}{\begin{array}{c}{\mathit{\psi }}_{\mathrm{LP}}\mathrm{LP}+{\mathit{\psi }}_{\mathrm{CP}}\mathrm{CP}+{\mathit{\psi }}_{\mathrm{DZ}}\mathrm{DZ}+{\mathit{\psi }}_{\mathrm{DT}}\mathrm{DT}\\ +{\mathit{\psi }}_{\mathrm{Det}}{\mathrm{Detr}}_{\mathrm{tot}}+{\mathit{\psi }}_{\mathrm{Z}}Z+k_z\end{array}}}.\end{array}$$

The calculated maximum potential grazing rate is a function of a maximum potential grazing rate at 0 ^∘C (${\mathit{\mu }}_Z^{\mathit{\theta }}$), temperature and oxygen, where grazing activity is capped when temperatures exceed 20 ^∘C (Keller et al., 2012):

$$\begin{array}{}\text{(20)}& {\mathit{\mu }}_Z^{\max }={\mathit{\mu }}_Z^{\mathit{\theta }}\max (\mathrm{0},r_{\mathrm{sox}}^{{\mathrm{O}}_{\mathrm{2}}}{b}^{c\cdot min(\mathrm{20}{}^{\circ }\mathrm{C},T)}).\end{array}$$

Grazing is also reduced under suboxic conditions ($r_{\mathrm{sox}}^{{\mathrm{O}}_{\mathrm{2}}}$):

$$\begin{array}{}\text{(21)}& r_{\mathrm{sox}}^{{\mathrm{O}}_{\mathrm{2}}}=\mathrm{0.5}\left(\tanh \right(\left[{\mathrm{O}}_{\mathrm{2}}\right]-\mathrm{8}\mathrm{mmol}{\mathrm{m}}^{-\mathrm{3}})+\mathrm{1}),\end{array}$$

where O₂ is dissolved oxygen in mmol m⁻³.

Table 3 KMBM3 phytoplankton production and mortality parameters.

Table 4 KMBM3 zooplankton parameters. Temperature-dependent parameter values are given for 0 ^∘C.

2.2.4 Organic detritus

As was introduced in Kvale et al. (2015b), organic carbon detritus sources and sinks are split into “free” (Detr_free) and “ballast” (Detr_bal) pools using a fixed ratio (R_{bal : tot}). Ballast detritus is formed from the CaCO₃-protected portion 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₃ dissolution (${\mathit{\lambda }}_{{\mathrm{CaCO}}_{\mathrm{3}}}$):

$$\begin{array}{}\text{(22)}& {\displaystyle }{\mathrm{Detr}}_{\mathrm{tot}}={\mathrm{Detr}}_{\mathrm{bal}}+{\mathrm{Detr}}_{\mathrm{free}}\text{(23)}& {\displaystyle }\begin{array}{rl}S\left({\mathrm{Detr}}_{\mathrm{bal}}\right)& =(\mathrm{1}-\mathit{\gamma })(G_{{\mathrm{Detr}}_{\mathrm{bal}}}+(G_{\mathrm{Z}}+G_{\mathrm{CP}}\left)R_{\mathrm{bal}:\mathrm{tot}}\right)\\ & +(m_{\mathrm{Z}}Z+m_{\mathrm{CP}}\mathrm{CP})R_{\mathrm{bal}:\mathrm{tot}}-G_{{\mathrm{Detr}}_{\mathrm{bal}}}\\ & -{\displaystyle \frac{R_{\mathrm{bal}:\mathrm{tot}}{\mathit{\lambda }}_{{\mathrm{CaCO}}_{\mathrm{3}}}\left[{\mathrm{CaCO}}_{\mathrm{3}}\right]}{R_{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{POC}}{R}_{\mathrm{C}:\mathrm{N}}}}-w_{\mathrm{C}}{\displaystyle \frac{\partial {\mathrm{Detr}}_{\mathrm{bal}}}{\partial z}},\end{array}\end{array}$$

where γ is the food assimilation efficiency, $R_{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{POC}}$ is a fixed production ratio of CaCO₃ to organic detritus, R_C : N is a Redfield molar ratio, and w_C is the sinking speed of calcite. Free detritus is described as

$$\begin{array}{}\text{(24)}& \begin{array}{rl}S\left({\mathrm{Detr}}_{\mathrm{free}}\right)& =(\mathrm{1}-\mathit{\gamma })(G_{\mathrm{Z}}+G_{\mathrm{DZ}}+G_{\mathrm{DT}}+G_{{\mathrm{Detr}}_{\mathrm{free}}}\\ & +G_{\mathrm{Z}}(\mathrm{1}-R_{\mathrm{bal}:\mathrm{tot}})+G_{\mathrm{CP}}(\mathrm{1}-R_{\mathrm{bal}:\mathrm{tot}}))\\ & +m_{\mathrm{Z}}{Z}^{\mathrm{2}}(\mathrm{1}-R_{\mathrm{bal}:\mathrm{tot}})+m_{\mathrm{DZ}}\mathrm{DZ}+m_{\mathrm{DT}}\mathrm{DT}\\ & +m_{\mathrm{LP}}\mathrm{LP}+m_{\mathrm{CP}}\mathrm{CP}(\mathrm{1}-R_{\mathrm{bal}:\mathrm{tot}})-{\mathit{\mu }}_{\mathrm{D}}{\mathrm{Detr}}_{\mathrm{free}}\\ & -G_{{\mathrm{Detr}}_{\mathrm{free}}}+{\displaystyle \frac{R_{\mathrm{bal}:\mathrm{tot}}{\mathit{\lambda }}_{{\mathrm{CaCO}}_{\mathrm{3}}}\left[{\mathrm{CaCO}}_{\mathrm{3}}\right]}{R_{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{POC}}{R}_{\mathrm{C}:\mathrm{N}}}}\\ & -w_{\mathrm{D}}{\displaystyle \frac{\partial {\mathrm{Detr}}_{\mathrm{free}}}{\partial z}},\end{array}\end{array}$$

where μ_D is the detrital remineralization rate and w_D is the sinking speed of organic detritus.

2.2.5 Detrital iron

Iron in detritus follows Nickelsen et al. (2015) but with the additional contribution of calcifiers and diatoms:

$$\begin{array}{}\text{(25)}& \begin{array}{rl}S\left({\mathrm{Detr}}_{\mathrm{Fe}}\right)& =R_{\mathrm{Fe}:\mathrm{N}}\left(\right(\mathrm{1}-\mathit{\gamma }\left)\right(G_{\mathrm{LP}}+G_{\mathrm{DZ}}+G_{\mathrm{DT}}\\ & +G_{{\mathrm{Detr}}_{\mathrm{tot}}}+G_{\mathrm{Z}}+G_{\mathrm{CP}})+m_{\mathrm{LP}}\mathrm{LP}+m_{\mathrm{DZ}}\mathrm{DZ}\\ & +m_{\mathrm{DT}}\mathrm{DT}+m_{\mathrm{Z}}{Z}^{\mathrm{2}}+m_{\mathrm{CP}}\mathrm{CP}-G_{{\mathrm{Detr}}_{\mathrm{tot}}})\\ & +[\mathrm{Fe}{]}_{\mathrm{orgads}}+[\mathrm{Fe}{]}_{{\mathrm{orgads}}_{\mathrm{ca}}}+[\mathrm{Fe}{]}_{\mathrm{col}}\\ & -{\mathit{\mu }}_{\mathrm{D}}{\mathrm{Detr}}_{\mathrm{Fe}}-w_{\mathrm{D}}{\displaystyle \frac{\partial {\mathrm{Detr}}_{\mathrm{Fe}}}{\partial z}}.\end{array}\end{array}$$

R_Fe : N is a fixed iron to nitrogen ratio that converts total organic nitrogen detritus sources into iron units. The remineralization and sinking of detrital iron occurs at the same rates as free organic detritus. Two scavenging processes are considered: adsorption of dissolved iron onto organic detritus ([Fe]_orgads):

$$\begin{array}{}\text{(26)}& [\mathrm{Fe}{]}_{\mathrm{orgads}}={\mathrm{kFe}}_{\mathrm{org}}[{\mathrm{Fe}}^{\prime }\left]\right({\mathrm{Detr}}_{\mathrm{free}}{R}_{\mathrm{C}:\mathrm{N}},M_{\mathrm{C}}{)}^{\mathrm{0.58}}\end{array}$$

and particulate CaCO₃ ($[\mathrm{Fe}{]}_{{\mathrm{orgads}}_{\mathrm{ca}}}$):

$$\begin{array}{}\text{(27)}& [\mathrm{Fe}{]}_{{\mathrm{orgads}}_{\mathrm{ca}}}=k{\mathrm{Fe}}_{\mathrm{ca}}[{\mathrm{Fe}}^{\prime }\left]\right(\left[{\mathrm{CaCO}}_{\mathrm{3}}\right]M_{{\mathrm{CaCO}}_{\mathrm{3}}}{)}^{\mathrm{0.58}},\end{array}$$

calculated as a function of free organic detritus and carbonate particles converted to carbon units, respectively, a scavenging rate constant (kFe_org and kFe_ca), iron available for scavenging ([Fe^′]) 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. (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.

2.2.6 Calcite

As in Kvale et al. (2015b), the source and sink terms for CaCO₃ include both phytoplankton calcifier and zooplankton sources from grazing and mortality, and losses from dissolution and sinking:

$$\begin{array}{}\text{(28)}& \begin{array}{rl}S\left(\right[{\mathrm{CaCO}}_{\mathrm{3}}\left]\right)& =\left(\right(\mathrm{1}-\mathit{\gamma }\left)\right(G_{\mathrm{CP}}+G_{\mathrm{Z}})+m_{\mathrm{CP}}\mathrm{CP}+m_{\mathrm{Z}}Z)\\ & \cdot R_{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{POC}}{R}_{\mathrm{C}:\mathrm{N}}-{\mathit{\lambda }}_{{\mathrm{CaCO}}_{\mathrm{3}}}\left[{\mathrm{CaCO}}_{\mathrm{3}}\right]\\ & -w_{\mathrm{C}}{\displaystyle \frac{\partial \left[{\mathrm{CaCO}}_{\mathrm{3}}\right]}{\partial z}},\end{array}\end{array}$$

where w_C is the sinking speed of calcite, and the dissolution-specific rate (${\mathit{\lambda }}_{{\mathrm{CaCO}}_{\mathrm{3}}}$) is calculated using a parameterization from Aumont et al. (2003) and Kvale et al. (2015b).

Calcite associated with living plankton is calculated separately (but not traced) as

$$\begin{array}{}\text{(29)}& S[{\mathrm{CaCO}}_{\mathrm{3}}{]}_{\mathrm{liv}}=(S\left(\mathrm{CP}\right)+S\left(Z\right))R_{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{POC}}{R}_{\mathrm{C}:\mathrm{N}}.\end{array}$$

2.2.7 Opal

The production of biogenic opal is calculated as a function of the diatom grazing and linear mortality loss terms:

$$\begin{array}{}\text{(30)}& \mathrm{Pr}\left(\mathrm{Opal}\right)=\left(\right(\mathrm{1}-\mathit{\gamma })G_{\mathrm{DT}}+m_{\mathrm{DT}}\mathrm{DT})R_{\mathrm{Opal}:\mathrm{POC}}{R}_{\mathrm{C}:\mathrm{N}},\end{array}$$

where R_{Opal : POC} is a production ratio that varies as

$$\begin{array}{}\text{(31)}& \begin{array}{rl}& R_{\mathrm{Opal}:\mathrm{POC}}=\\ & R_{\mathrm{Opal}:\mathrm{POC},\mathrm{0}}\left(min\left({\displaystyle \frac{\left[\mathrm{Si}\right]}{k_{\mathrm{Si}}}},\mathrm{1}\right)\left(\mathrm{4}\text{–}\mathrm{3}min\left({\displaystyle \frac{\left[\mathrm{Fe}\right]}{k_{\mathrm{Fe}DT}}},\mathrm{1}\right)\right)\right).\end{array}\end{array}$$

This parameterization was introduced by Aumont et al. (2003) and yields an average surface opal : free detritus export value of around 1 Si : C across the Southern Ocean, using a fixed average ratio (R${}_{\mathrm{Opal}:\mathrm{POC},\mathrm{0}}$) 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. In linking opal production to diatom mortality, our model essentially focuses on the importance of its export pathway in establishing vertical gradients. Opal that is not exported is effectively assumed to be recycled to dissolved silicate within the surface layer.

Dissolution of opal in the water column is calculated by assuming instantaneous sinking of the vertically integrated production, where the local flux of opal (F_Opal in mol Si m⁻²) is distributed down the water column using the e-folding temperature parameterization (unitless), scaled by a dissolution rate constant (λ_Opal, in day⁻¹), which is divided by a sinking rate (w_Opal, in metres day⁻¹):

$$\begin{array}{}\text{(32)}& \mathrm{Di}\left(\mathrm{Opal}\right)=F_{\mathrm{Opal}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}=F_{\mathrm{Opal}}\left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{Opal}}}{w_{\mathrm{Opal}}}}{\exp }^{\frac{T}{T_b}}\right).\end{array}$$

Dissolution is presented in units of mol Si m⁻³ d⁻¹. This parameterization results in greater dissolution at warm temperatures and is similar to the instant-sinking-and-dissolution function applied to model calcite (Schmittner et al., 2008). Opal dissolution in the water column is primarily controlled by temperature and silicate saturation, with higher temperatures and lower silicate saturations leading to faster dissolution (Sarmiento and Gruber, 2006; Ridgwell et al., 2002). Secondary drivers include the aluminium content and surface area of diatom shells (Van Cappellen et al., 2002), and the presence or absence of organic coatings, where the loss of the coating increases opal dissolution rates (Bidle and Azam, 1999). We do not consider silicate or aluminium concentration dependencies explicitly in this model because we do not explicitly calculate thermodynamic dissolution of opal (nor do we simulate aluminium), but the temperature dependency exerts an influence on dissolution of the same sign (faster dissolution in warmer, i.e. shallower and undersaturated, water). We find this parameterization offers improved model fit to World Ocean Atlas silica distributions relative to other parameterizations that we tested, e.g. the temperature-dependent parameterization of Gnanadesikan (1999) or the temperature- and oxygen-dependent parameterization of Enright et al. (2014). The Gnanadesikan (1999) parameterization yields lower dissolution rates at low temperatures than the Enright et al. (2014) parameterization, which is similarly formulated but includes additional oxygen scaling. The Enright et al. (2014) oxygen scaling is not justified in their model description, but it has the effect of increasing Si dissolution rates in the deep ocean (exacerbating the over-estimation 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 area. Our temperature scaling has the effect of raising dissolution rates at the surface. Greater dissolution rates at the surface may be necessary to compensate for the low vertical resolution of the model. Global average implicit opal export, dissolution and seafloor loss rate profiles are plotted in Fig. 2 for model year 2014. Consistent with the estimates of Tréguer et al. (2021), opal cycling in the open ocean largely occurs in the shallower depths, with opal export following a power law decrease with increasing depth. The global ocean average opal dissolution rate at 1000 m depth is only half that at the ocean surface, and the average flux of opal at the seafloor at 500 m depth is only half that at the surface. Shallow shelf areas are hence an important and relatively fast region of turnover for our simulated silicate cycle.

Figure 2 Implicit opal export (a), water column dissolution (b) and seafloor burial (c) at model year 2014. Rates are globally integrated in the latitudinal and longitudinal planes but depth information is retained. Thus, opal export, dissolution and seafloor flux are reported in Tmol Si flux per metre of depth per year.

2.2.8 Particle sinking

Detritus (Schmittner et al., 2005), calcite (Kvale et al., 2015b) and iron (Nickelsen et al., 2015) particles are exported from the surface with a sinking speed (w) that increases linearly (wdc, wdd) with depth (z; Berelson, 2001) for calcite and ballasted organic detritus:

$$\begin{array}{}\text{(33)}& w_{\mathrm{C}}=w_{C,\mathrm{0}}+wdc\times z,\end{array}$$

and free detritus and associated iron:

$$\begin{array}{}\text{(34)}& w_{\mathrm{D}}=w_{D,\mathrm{0}}+wdd\times z.\end{array}$$

Alternative parameterizations exist and their effects on fluxes and model performance make for interesting comparisons (e.g. Kriest and Oschlies, 2011; Cael and Bisson, 2018), but we do not explore them here. The initial surface sinking speeds of particles are assigned different values to represent the denser structure of CaCO₃ relative to that of organic detritus. Ballasted detritus sinks at the CaCO₃ speed, but once it enters the free pool it uses the organic detrital sinking speed and remineralization rate. Any organic detritus reaching the sediments is dissolved back in to the water column to ensure conservation of carbon and phosphate in the model domain. 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 instantly at the seafloor). Iron detritus reaching the sediments is lost from the ocean, unless bottom water oxygen falls below 5 mmol m⁻³, whereupon detrital iron reaching the bottom is 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):

$$\begin{array}{}\text{(35)}& [\mathrm{Fe}{]}_{\mathrm{sed}}=R_{\mathrm{Fe}:{\mathrm{P}}_{\mathrm{sed}}}{F}_{\mathrm{POP}}{\exp }^{\frac{T}{T_b}},\end{array}$$

where $R_{\mathrm{Fe}:{\mathrm{P}}_{\mathrm{sed}}}$ 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 ([Si]_sed) according to criteria laid out in Sarmiento and Gruber (2006):

$$\begin{array}{}\text{(36)}& \begin{array}{rl}& \mathrm{For}F\left(\mathrm{Opal}\right)>\mathrm{2}\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}[\mathrm{Si}{]}_{\mathrm{sed}}=\mathrm{0.3}F(\mathrm{Opal})\\ & \mathrm{For}F\left(\mathrm{Opal}\right)\le \mathrm{2}\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}[\mathrm{Si}{]}_{\mathrm{sed}}=\mathrm{0.05}F(\mathrm{Opal}).\end{array}\end{array}$$

Opal flux, F(Opal), is calculated as depth-integrated opal production (Pr(Opal)) less the depth-integrated 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 subgrid 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 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 paper.

2.2.9 Dissolved inorganic tracers

Ocean nutrient sources and sinks (in concentration/time units) follow

$$\begin{array}{}\text{(37)}& {\displaystyle }\begin{array}{rl}S\left(\right[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\left]\right)& =({\mathit{\mu }}_{\mathrm{D}}{\mathrm{Detr}}_{\mathrm{free}}+{\mathit{\mu }}_{\mathrm{LP}}^{*}\mathrm{LP}+{\mathit{\mu }}_{\mathrm{CP}}^{*}\mathrm{CP}+{\mathit{\mu }}_{\mathrm{DT}}^{*}\mathrm{DT}\\ & +(\mathit{\gamma }-\mathit{\varpi })(G_{\mathrm{LP}}+G_{\mathrm{CP}}+G_{\mathrm{DT}}+G_{\mathrm{DZ}}\\ & +G_{{\mathrm{Detr}}_{\mathrm{free}}}+G_{\mathrm{Z}})-J_{\mathrm{LP}}\mathrm{LP}-J_{\mathrm{CP}}\mathrm{CP}\\ & -J_{\mathrm{DT}}\mathrm{DT}-J_{\mathrm{DZ}}\mathrm{DZ})R_{\mathrm{P}:\mathrm{N}}\end{array}\text{(38)}& {\displaystyle }\begin{array}{rl}S\left(\right[{\mathrm{NO}}_{\mathrm{3}}^{-}\left]\right)& =({\mathit{\mu }}_{\mathrm{D}}{\mathrm{Detr}}_{\mathrm{free}}+{\mathit{\mu }}_{\mathrm{LP}}^{*}\mathrm{LP}+{\mathit{\mu }}_{\mathrm{CP}}^{*}\mathrm{CP}+{\mathit{\mu }}_{\mathrm{DT}}^{*}\mathrm{DT}\\ & +(\mathit{\gamma }-\mathit{\varpi })(G_{\mathrm{LP}}+G_{\mathrm{CP}}+G_{\mathrm{DT}}+G_{\mathrm{DZ}}\\ & +G_{{\mathrm{Detr}}_{\mathrm{free}}}+G_{\mathrm{Z}})-J_{\mathrm{LP}}\mathrm{LP}-J_{\mathrm{CP}}\mathrm{CP}\\ & -J_{\mathrm{DT}}\mathrm{DT}-u_{{\mathrm{NO}}_{\mathrm{3}}}{J}_{\mathrm{DZ}}\mathrm{DZ}\left)\right(\mathrm{1}-\mathrm{0.8}{R}_{\mathrm{O}:\mathrm{N}}{r}_{\mathrm{sox}}^{{\mathrm{NO}}_{\mathrm{3}}}),\end{array}\end{array}$$

where R_P : N and R_O : N are Redfield molar ratios and $u_{{\mathrm{NO}}_{\mathrm{3}}}$ is nitrate availability (diazotrophs use nitrate when available). In suboxic water (less than 5 mmol m⁻³), oxygen consumption is replaced by denitrification ($r_{\mathrm{sox}}^{{\mathrm{NO}}_{\mathrm{3}}^{-}}$):

$$\begin{array}{}\text{(39)}& r_{\mathrm{sox}}^{{\mathrm{NO}}_{\mathrm{3}}^{-}}=\max (\mathrm{0},\mathrm{0.5}(\mathrm{1}-\tanh \left(\right[{\mathrm{O}}_{\mathrm{2}}]-\mathrm{5}{\mathrm{mmol}\mathrm{m}}^{-\mathrm{3}})\left)\right).\end{array}$$

There are no additions of phosphate, nitrate or oxygen along the boundary (with the exception of air/sea gas exchange in the case of oxygen; Keller et al., 2012).

Dissolved iron includes sources and sinks of particulate iron mentioned above, as well as prescribed dust deposition ([Fe]_dust) as in Nickelsen et al. (2015) and hydrothermal iron ([Fe]_hydr) (Muglia et al., 2017) boundary terms (also in concentration/time units):

$$\begin{array}{}\text{(40)}& {\displaystyle }\begin{array}{rl}S\left(\right[\mathrm{Fe}\left]\right)& =R_{\mathrm{Fe}:\mathrm{N}}({\mathit{\mu }}_{\mathrm{D}}{\mathrm{Detr}}_{\mathrm{tot}}+{\mathit{\mu }}_{\mathrm{LP}}^{*}\mathrm{LP}+{\mathit{\mu }}_{\mathrm{CP}}^{*}\mathrm{CP}\\ & +{\mathit{\mu }}_{\mathrm{DT}}^{*}\mathrm{DT}+(\mathit{\gamma }-\mathit{\varpi })(G_{\mathrm{LP}}+G_{\mathrm{CP}}+G_{\mathrm{DT}}\\ & +G_{\mathrm{DZ}}+G_{{\mathrm{Detr}}_{\mathrm{tot}}}+G_{\mathrm{Z}})-J_{\mathrm{LP}}\mathrm{LP}-J_{\mathrm{CP}}\mathrm{CP}\\ & -J_{\mathrm{DT}}\mathrm{DT}-J_{\mathrm{DZ}}\mathrm{DZ})-[\mathrm{Fe}{]}_{\mathrm{orgads}}\\ & -[\mathrm{Fe}{]}_{{\mathrm{orgads}}_{\mathrm{ca}}}-[\mathrm{Fe}{]}_{\mathrm{col}}+{\mathit{\mu }}_{\mathrm{D}}{\mathrm{Detr}}_{\mathrm{Fe}}\end{array}\text{(41)}& {\displaystyle }B\left(\right[\mathrm{Fe}\left]\right)=[\mathrm{Fe}{]}_{\mathrm{sed}}+[\mathrm{Fe}{]}_{\mathrm{dust}}+[\mathrm{Fe}{]}_{\mathrm{hydr}}.\end{array}$$

The Nickelsen et al. (2015) model applied the iron dust flux to the surface ocean after the biological routine and before the mixing routine. This resulted in low model sensitivity to iron dust inputs, because the dust was mixed away faster than the biological processes 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 time step 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₃ (Kvale et al., 2015b).

$$\begin{array}{}\text{(42)}& {\displaystyle }\begin{array}{rl}S\left(\right[\mathrm{DIC}\left]\right)& =S\left(\right[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\left]\right)R_{\mathrm{C}:\mathrm{P}}+{\mathit{\lambda }}_{{\mathrm{CaCO}}_{\mathrm{3}}}\left[{\mathrm{CaCO}}_{\mathrm{3}}\right]\\ & -S[{\mathrm{CaCO}}_{\mathrm{3}}{]}_{\mathrm{liv}}-((\mathrm{1}-\mathit{\gamma })(G_{\mathrm{Z}}+G_{\mathrm{CP}})\\ & +m_{\mathrm{Z}}Z+m_{\mathrm{CP}}\mathrm{CP})R_{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{POC}}{R}_{\mathrm{C}:\mathrm{N}}\end{array}\text{(43)}& {\displaystyle }\begin{array}{rl}S\left(\right[\mathrm{Alk}\left]\right)& =-S\left(\right[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\left]\right)R_{\mathrm{N}:\mathrm{P}}+\mathrm{2}\left({\mathit{\lambda }}_{{\mathrm{CaCO}}_{\mathrm{3}}}\right[{\mathrm{CaCO}}_{\mathrm{3}}]\\ & -S\left[{\mathrm{CaCO}}_{\mathrm{3}}{]}_{\mathrm{liv}}\right)-\mathrm{2}\left(\right(\mathrm{1}-\mathit{\gamma }\left)\right(G_{\mathrm{Z}}+G_{\mathrm{CP}})\\ & +m_{\mathrm{Z}}Z+m_{\mathrm{CP}}\mathrm{CP})R_{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{POC}}{R}_{\mathrm{C}:\mathrm{N}}\end{array}\end{array}$$

With the exception of air–sea gas exchange of CO₂, there are no boundary additions of DIC or alkalinity when the sediment model is deactivated (as is presented here).

Dissolved silicic acid tracer sources, sinks and boundary terms (all in concentration/time units) follow

$$\begin{array}{}\text{(44)}& {\displaystyle }S\left(\right[\mathrm{Si}\left]\right)=\mathrm{Di}\left(\mathrm{Opal}\right)-\mathrm{Pr}\left(\mathrm{Opal}\right)\text{(45)}& {\displaystyle }B\left(\right[\mathrm{Si}\left]\right)=[\mathrm{Si}{]}_{\mathrm{riv}}+[\mathrm{Si}{]}_{\mathrm{dust}}-[\mathrm{Si}{]}_{\mathrm{sed}}+[\mathrm{Si}{]}_{\mathrm{hydr}},\end{array}$$

where discharge from rivers ([Si]_riv) is used as a budget balancing term to compensate for any remainder in the other external sources and sinks: wind-borne dust ([Si]_dust), hydrothermal silicate ([Si]_hydr) and loss to the sediments ([Si]_sed). River inputs of silicate, alkalinity and DIC are scaled against seasonally variable river flow using the standard UVic ESCM version 2.9 O_rivflux.F forcing file. 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 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) (0.6 Tmol Si yr⁻¹). This estimate has recently been revised upward to 1.7 ± 0.8 Tmol Si yr⁻¹ (Tréguer et al., 2021); thus, our hydrothermal contribution is under-estimated. Model sensitivity to both sources will be explored in future tests.

3 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₂ concentrations, agricultural land cover, volcanic radiative forcing, sulfate aerosol and chlorofluorocarbon (CFC) concentrations, changes in land ice and solar forcing (Machida et al., 1995; Battle et al., 1996; Etheridge et al., 1996, 1998; Flückiger et al., 1999, 2004; Ferretti et al., 2005; Meure et al., 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 6 lists key biogeochemical properties diagnosed by the model for a modern climate, as well as corresponding properties diagnosed from KMBM1 (Keller et al., 2012). We compare model output at the year 2004, although data sources reflect a range of collection and publication dates. Net primary production (NPP) is similar to previous model versions (53.60 Pg C yr⁻¹, e.g. compared to 54.33 Pg C yr⁻¹ in Keller et al., 2012), and still within the literature range of 44–78 Pg C yr⁻¹ (e.g. Carr et al., 2006; Jin et al., 2006). Calcite production is similar to and nitrogen fixation rates are somewhat improved with respect to earlier model versions, though they remain too low, as do deep particle fluxes of particulate organic carbon (POC) and shallow and deep calcite (also referred to as PIC, for particulate inorganic carbon). Diagnosed surface opal production is also too low and below the range of a recent estimate (Tréguer et al., 2021). Accordingly, deep ocean (2 km depth) opal fluxes are also too low (only about 50 % of the observationally based deep ocean estimates of Tréguer and De La Rocha, 2013). However, upper ocean (130 m) opal flux is about 9 % larger than the estimate provided by Tréguer et al., 2021. This may be due to the model under-estimating upper-ocean dissolution (10.3 Tmol Si yr⁻¹ in our model compared to 124 Tmol Si yr⁻¹ estimated in the surface mixed layer by Tréguer et al., 2021). There is too little dissolution of opal throughout the water column (74 % of the observational estimate in Tréguer et al., 2021, but the ratio of total water column dissolution to biogenic production is too high, with a ratio of 0.95 compared to 0.67 calculated by Tréguer et al., 2021). The calculated river flux is 2.38 Tmol Si yr⁻¹; this is lower than the Tréguer et al. (2021) estimate of 8.1 Tmol Si yr⁻¹. Overall, an increase in production and subsequent processes is required to improve the model agreement with independently estimated fluxes. Low bias in fluxes is also found in calcite production and nitrogen fixation, as well as in deep particle fluxes, which is due at least in part to the model bias of a too-sluggish deep circulation (Kriest et al., 2020). The calculated flux of silicate through the seafloor is only 70 % of that estimated by Tréguer et al. (2021). An under-estimation of seafloor flux is potentially corrected with an increase in pelagic opal production rates or an decrease in benthic return. Benthic fluxes are potentially improved with newer models of benthic transfer (e.g. Dale et al., 2021) and will be explored in the future. A recent review suggests the global silicate cycle sources and sinks are roughly balanced, but uncertainties are still substantial (Tréguer et al., 2021). It should be noted that the global silica cycle may be unbalanced between sources and sinks at present and observations of global input and output rates are both affected by anthropogenic perturbation and have large uncertainties (Tréguer and De La Rocha, 2013). Atmospheric and hydrothermal silica inputs are fixed in these simulations, but each source and sink of silicate produce a unique spatial distribution in the water column. These terms are slated for automated calibration in the future (e.g. Yao et al., 2019; Kriest, 2017; Kriest et al., 2017); thus, our hand-tuned simulations are meant to serve as a baseline for future improvements to the model. Also note that despite the deficiencies described above, KMBM3 demonstrates reasonable performance against a suite of state-of-the-art Earth system models (described below).

Table 5 KMBM3 particle export–production parameters.

Table 6 Globally integrated diagnosed biogeochemical properties at the year 2004. Corresponding values are also given using KMBM1 (Keller et al., 2012).

^a Carr et al. (2006), Jin et al. (2006). ^b Smith and Gattuso (2011), Smith and Mackenzie (2016). ^c Tréguer et al. (2021). ^d Luo et al. (2014). ^e All particle fluxes from Honjo et al. (2008) unless noted. ^f Luo et al. (2014). ^g All biomass estimates from Buitenhuis et al. (2013). ^h Picophytoplankton. ⁱ Pteropods.

Total global phytoplankton biomass is on the low end, but within the range, of previous estimates (0.55 Pg C, virtually unchanged from Keller et al., 2012, compared to 0.5–2.4 Pg C from Buitenhuis et al., 2013; Table 6). KMBM3 explicitly represents only a fraction of the ecological complexity found in the real ocean, which causes our biomass estimates for the phytoplankton functional types to look dissimilar to the Buitenhuis et al. (2013) biomass estimates (e.g. with low-latitude mixed phytoplankton (a model-specific category, LP) having a biomass only 11 % of the lowest picophytoplankton estimate and both diazotrophs (DZ) and calcifiers (CP) having more biomass than the upper observational estimate). These functional types proved particularly difficult to tune by hand, with small variations in parameter values causing extinction. An over-estimate of calcifiers is compensated by the low PIC : POC production ratio (0.07), which is meant to represent 7 % of the phytoplankton class having a PIC : POC production ratio of 1 (compare to a ratio for Emiliania huxleyi of 0.51–2.30 from Paasche, 2001). Over-estimated diazotroph biomass results in an increase (0.04 Pg N yr⁻¹) in nitrogen fixation compared to earlier model versions, which is improved with respect to, but still lower than, the observational independent estimate. In KMBM3, diazotrophs use preformed nitrate when available. Thus in our modelling context, this functional type can be considered “slow-growing phytoplankton capable of fixing nitrogen when necessary”. Constraints on this functional type will be explored in the future. Diatom biomass estimates are too low and outside the Buitenhuis et al. (2013) range, although the functional type is represented in the same relative proportion to low-latitude phytoplankton (LP) as is reported for picophytoplankton to diatoms by Buitenhuis et al. (2013).

Looking next at spatial distributions of biological rates, Fig. 3 compares KMBM3 NPP at year 2014 to the Westberry et al. (2008) model applied to MODIS (NASA, 2018) climatology from 2012–2018. As is repeatedly found in KMBM (see plot of the Keller et al., 2012 model NPP), open-ocean primary production rates are generally too high, particularly in the eastern equatorial Pacific and northern Indian Ocean. These very-high-production zones compensate in the globally integrated rate estimate for the diffuse and more widespread production calculated from satellites. The spatial pattern in NPP has changed somewhat from earlier model versions, with a continued under-estimate of NPP occurring within gyres and a new over-estimate of NPP in the 40–60^∘ N and S ranges.

Figure 3 Comparison of model NPP output at 2014 (a Keller et al., 2012; b this model) and satellite-derived NPP climatology (NASA, 2018), 2012–2018 (c) in gC m⁻² d⁻¹.

Spatial biases in NPP are also held in surface calcite concentrations. Figure 4 compares KMBM3 CaCO₃ averaged between 2004–2014 with the 2002–2018 climatology data from MODIS (NASA, 2018). Model calcite is too high in the North Pacific and Southern Ocean between 40 and 60^∘ S, and too low south of 60^∘ S. Upwelling zones and the Indian Ocean have an over-estimated role in calcite production, and coastal calcite production is either not resolved or under-estimated.

Figure 4 Comparison of model average 2004–2014 (a) and satellite CaCO₃ climatology 2002–2018 (b) in mmol CaCO₃ m⁻³. Data product is scaled by the model grid in the z direction, and in both plots only the upper 20 m are represented, where uniform coccolith concentration is assumed (Balch and Utgoff, 2009).

Phytoplankton biomass is compared to the MARine Ecosystem DATa (MAREDAT) datasets (Leblanc et al., 2012; Luo et al., 2012; O'Brien, 2012) in Fig. 5. Globally integrated diatom biomass is too low compared to observational estimates but the diatom geographical distributions agree roughly with the very sparse Leblanc et al. (2012) dataset (highest concentrations in the Southern Ocean, North Pacific and North Atlantic). Like the globally integrated biomass estimate, calcifiers (CP) are universally over-estimated compared to O'Brien (2012). As explained above, this functional 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). Mixed low-latitude phytoplankton (LP) is a model-specific category with no clear analogue in the MAREDAT dataset.

Figure 5 Annually averaged and depth-integrated model biomass at the year 2014 (top row). Zonally averaged and depth-integrated model biomass plotted with the equivalent MAREDAT compilations (bottom row; Leblanc et al., 2012; O'Brien, 2012; Luo et al., 2012; no data compilation of “mixed” phytoplankton, a model-specific category). All units are mmol C m⁻². Black lines represent MAREDAT; red lines are the model output.

We compare CMIP6 model output relative to KMBM3 model output due to the very low normalized correlation (0.04–0.15) and normalized standard deviation (0.10–0.22) of all models against the sparse Leblanc et al. (2012) diatom biomass dataset. At the time of writing, available CMIP6-simulated annual average diatom biomass shows diverse quantities (maximum concentrations from 0.0035 to 0.03 mol C m⁻³) and spatial distributions ranging from global maximum concentrations at the Equator (CanESM5-CanOE, Swart et al., 2019), to shallow seas and coastlines (IPSL-CM6A-LR, Boucher et al., 2018), to the high latitudes (CMCC-ESM2, EC-Earth3-CC, GFDL-ESM4, CESM2; Lovato et al., 2021; EC-Earth Consortium, 2021; Krasting et al., 2018; Danabasoglu, 2019, respectively). Figure 6 contextualizes these diatom biomass estimates relative to the KMBM3 diatom biomass, simulated at year 2014. KMBM3 is most closely correlated with GFDL-ESM4. KMBM3 has a near-perfect standard deviation and second-closest correlation with CMCC-ESM2, a higher-resolution fully coupled Earth system model with full representation of the global carbon cycle and ocean biogeochemistry. The CMCC-ESM4, KMBM3, CESM2 and CanESM5-CanOE models all fall outside the cluster in correlation and standard deviation of EC-Earth3-CC, IPSL-CM6A-LR and GFDL-ESM4. EC-Earth3-CC and IPSL-CM6A-LR both utilize the PISCES biogeochemical model (and are therefore expected to show similar correlation and standard deviation relative to a reference model), but GFDL-ESM4 utilizes COBALT (an independent biogeochemical model lineage). As stated above, phytoplankton biomass is difficult to tune for (particularly when multiple functional types are represented) due to the under-constrained parameter space, the theoretical nature of phytoplankton functional categories and the sparsity of gridded, annually averaged biomass datasets. Therefore, a wide range in model biomass estimates is expected across the CMIP6 ensemble.

Figure 6 Taylor diagram (Taylor, 2001) of CMIP6 annual average diatom biomass concentrations (in mol C m⁻³) at the year 2014 referenced to KMBM3. The distance to the origin represents the normalized standard deviation. Normalized correlation with KMBM3 is read from the azimuthal position. Perfect agreement with KMBM3 is a normalized standard deviation of 1 and a normalized correlation of 1. Normalization against the model, rather than observations, is shown due to the models having greater similarity with each other than with the very sparse observations of diatom biomass currently available. Models are KMBM3 (red circle), GFDL-ESM4 (green star; Krasting et al., 2018), CanESM5-CanOE (yellow star; Swart et al., 2019), CESM2 (light blue square; Danabasoglu, 2019), CMCC-ESM2 (orange triangle; Lovato et al., 2021), EC-Earth3-CC (blue diamond; EC-Earth Consortium, 2021) and IPSL-CM6A-LR (purple hexagon; Boucher et al., 2018). With the exception of KMBM3, the data for all models were obtained by mining the CMIP6 database (https://esgf-node.llnl.gov/search/cmip6/, last access: 5 February 2021) using the following search terms: CMIP/phydiat/historical/annual output. Data were accessed between 1–5 February 2021.

KMBM3 interior ocean particle flux performance is encouraging. Figure 7 compares model POC, PIC and opal fluxes at 2 km depth to the Honjo et al. (2008) data compilation. The model shows a general under-estimate with respect to the observations for all deep ocean particle fluxes, particularly for intermediate rates. Root mean square error is improved beyond Kvale et al. (2015a) for PIC (101.4 compared to 147.1 mmol C m⁻² yr⁻¹) and POC (96.7 versus 98.0 mmol C m⁻² yr⁻¹ in Kvale et al., 2015a). Root mean square error for opal flux at 2 km depth is 223.0 mmol Si m⁻² yr⁻¹. However, the model captures large-scale features of high flux rates in the north west Pacific and Southern Ocean, as well as moderate flux rates in the Indian and North Atlantic basins. KMBM3 demonstrates moderate skill at simulating annual average deep ocean opal fluxes relative to CMIP6 models (Fig. 8), with two of six models having a lower correlation with the Honjo et al. (2008) data compared to KMBM3 and four of six having a lower, or equal, normalized standard deviation. Two of six have an over-estimated bias in opal flux, relative to KMBM3. The GFDL-ESM4 performs best against the Honjo et al. (2008) dataset of the model simulations examined here.

Figure 7 Comparison of model POC (a, d), CaCO₃ (b, e) and opal (c, f) fluxes at the year 2004 and 2 km depth with Honjo et al. (2008) data. Red points are this model (KMBM3); blue points are the Keller et al. (2012) model (KMBM1). Statistics are given for KMBM3 only.

Figure 8 Taylor diagram (Taylor, 2001) of CMIP6 annual average opal flux (in mol Si m⁻² yr⁻²) at 2 km depth, for the year 2014, normalized to the Honjo et al. (2008) dataset. The distance to the origin represents the normalized standard deviation. Normalized correlation with observations is read from the azimuthal position. Perfect agreement with observations is a normalized standard deviation of 1 and a normalized correlation of 1. Models are KMBM3 (red circle), GFDL-ESM4 (green star; Krasting et al., 2018), CESM2 (light blue square; Danabasoglu, 2019), MPI-ESM1.2-LR (pink pentagon; Wieners et al., 2019), CMCC-ESM2 (orange triangle; Lovato et al., 2021), EC-Earth3-CC (blue diamond; EC-Earth Consortium, 2021) and IPSL-CM6A-LR (purple hexagon; Boucher et al., 2018). With the exception of KMBM3, the data for all models were obtained by mining the CMIP6 database (https://esgf-node.llnl.gov/search/cmip6/) using the following search terms: CMIP/expsi/historical/annual output. Data were accessed between 1–5 February 2021.

Particle flux rates, alongside the model's ocean circulation, impact ocean tracer distributions. Figure 9 shows carbon and nutrient profiles, globally averaged and for each basin, compared to Global Ocean Data Analysis Project (GLODAP) (Key et al., 2015; Lauvset et al., 2016) and World Ocean Atlas (WOA) (Garcia et al., 2014a, b) data. Nutrients are generally too high at depth, and oxygen is too low, partly as a result of over-estimated deep ocean POC flux rates in the Southern Ocean. Global root mean squared error is listed in the figure for each model dissolved inorganic tracer. All tracers except DIC, phosphate and oxygen show improvement with respect to the Keller et al. (2012) model version. Alkalinity is primarily affected by PIC flux and attenuation rates and shows relatively good agreement (global RMSE value of 0.34 µmol kg⁻¹) with observations. Silicic acid has a global root mean squared error of 0.388 mmol Si m⁻³. Concentrations are too high in the Atlantic by as much as 50 mmol Si m⁻³ and too low in the Indian Ocean by as much as 50 mmol Si m⁻³.

Figure 9 Model global and basin average nutrient and carbon depth profiles at the year 2014 (KMBM3; red lines), compared to model output from the Keller et al. (2012) version (KMBM1; blue lines) and observational data (black lines; Garcia et al., 2014a, b; Key et al., 2015; Lauvset et al., 2016). Global root mean square error is given below each column.

Greater basin and surface detail in carbon and nutrient concentrations is displayed in Figs. 10–15. Interior ocean alkalinity is too high by as much as 50 µmol kg⁻¹ and surface alkalinity too low by as much as 50 µmol kg⁻¹, in the Atlantic, as it was with previous model versions (Eby et al., 2009; Keller et al., 2012; Kvale et al., 2015a). Alkalinity in the interior and surface Pacific is too high by as much as 20 µmol kg⁻¹ along the Equator, whereas previously it was over-estimated at depth and under-estimated at the surface (Kvale et al., 2015a). The physical circulation has not changed since previously published model versions. Differences might be partly due to shifts in phytoplankton biogeography and calcification rates. Note, however, that surface alkalinity is also subject to evaporation–precipitation effects that might have also changed since previous model versions. The Indian Ocean interior alkalinity is under-estimated by as much as 60 µmol kg⁻¹. Deep ocean DIC is not improved with respect to past model versions (Fig. 11), though (at least compared to Kvale et al., 2015a) this discrepancy might be partly attributed to the use of a transient year 2014 model output that includes an anthropogenic signal, instead of the pre-industrial spinup. The deep North Pacific shows a low bias of up to 40 µmol kg⁻¹, similar to previous versions. Surface DIC anomalies are also similar to previous model versions, with DIC being too low in the western Pacific (a consequence of too-high export production in the model and physical biases) and too high in the surface Southern Ocean. As with earlier model versions, deep ocean phosphate concentrations (Fig. 12) are also generally too high, particularly in the Southern Ocean-sourced deep and intermediate water masses in the Indian and Pacific sectors. Nitrate, however, is too low in these water masses (Fig. 13), although in basin average the depth profiles of KMBM1 and KMBM3 are nearly identical for the Southern Ocean, Pacific and Indian basins below 3000 m. Oxygen anomalies mirror nutrient biases (Fig. 14), with oxygen being up to 50 mmol m⁻³ too low along the Southern Ocean particle export pathway and up to 50 mmol m⁻³ too high in the subsurface tropical ocean. The oxygen bias in the deep water masses might be addressed with better tuning of the biological production, export and flux parameters, but the tropical ocean deficiencies will require improvements to both the biogeochemistry as well as the physics (Oschlies et al., 2017).

Figure 10 Model (left column) year 2014 alkalinity (µmol kg⁻¹) averaged by basin compared to GLODAP (right column; Key et al., 2015; Lauvset et al., 2016). Regions are as follows: Atlantic (top row), Pacific (second row), Indian (third row) and global surface (bottom row).

Figure 11 Model (left column) year 2014 DIC (µmol kg⁻¹) averaged by basin compared to GLODAP (right column; Key et al., 2015; Lauvset et al., 2016). Regions are as follows: Atlantic (top row), Pacific (second row), Indian (third row) and global surface (bottom row).

Figure 12 Model (left column) year 2014 PO${}_{\mathrm{4}}^{\mathrm{3}-}$ (mmol m⁻³) averaged by basin compared to WOA (right column; Garcia et al., 2014b). Regions are as follows: Atlantic (top row), Pacific (second row), Indian (third row) and global surface (bottom row).

Figure 13 Model (left column) year 2014 NO${}_{\mathrm{3}}^{-}$ (mmol m⁻³) averaged by basin compared to WOA (right column; Garcia et al., 2014b). Regions are as follows: Atlantic (top row), Pacific (second row), Indian (third row) and global surface (bottom row).

Figure 14 Model (left column) year 2014 O₂ (mmol m⁻³) averaged by basin compared to WOA (right column; Garcia et al., 2014a). Regions are as follows: Atlantic (top row), Pacific (second row), Indian (third row) and global surface (bottom row).

Silicic acid distributions are reasonably well captured by the model (Figs. 15, 16). KMBM3 demonstrates the second-highest correlation with gridded observations (Garcia et al., 2014b) of dissolved silicic acid at year 2014 compared to CMIP6 model output (Fig. 16). Two of seven models, MPI-ESM1.2-LR and CMCC-ESM2, have a greater difference in normalized standard deviation than KMBM3. Both correlation and standard deviation are generally better and show less inter-model spread for dissolved silicic acid than for other metrics. KMBM3 simulates the deep Southern Ocean maximum to within 20 mmol Si m⁻³, with a slight high bias in the deep Atlantic, Pacific and especially Indian sectors (Fig. 15). This spatial weighting to the Indian sector might represent the high bias in diatom biomass in this basin (also clearly seen in the opal flux plots). North Atlantic silicic acid gradients are well represented, while a low bias is apparent in the deep North Pacific of around 20 mmol Si m⁻³. A low bias is also simulated in the surface North Pacific, which possibly suggests deficiencies in the circulation within and between regional marginal seas (Nishioka et al., 2020). The North Pacific appears to perform well with respect to deep opal flux (Fig. 7), but small biases in fluxes can compound with time in dissolved nutrient fields. It might also be that sedimentary processes that influence deep water silicic acid concentrations would be important to resolve in this region. Our low-bias result in the deep North Pacific is interesting because the origins of the silicic acid-rich deep ocean plume in this region are still under debate (Tréguer and De La Rocha, 2013).

Figure 15 Model (left column) year 2014 Si (mmol m⁻³) averaged by basin compared to WOA (right column; Garcia et al., 2014b). Regions are as follows: Atlantic (top row), Pacific (second row), Indian (third row) and global surface (bottom row).

Figure 16 Taylor diagram (Taylor, 2001) of CMIP6 annual average dissolved silicic acid distribution (in mol Si m⁻³) for year 2014, normalized to the Garcia et al. (2014b) dataset. The distance to the origin represents the normalized standard deviation. Normalized correlation with observations is read from the azimuthal position. Perfect agreement with observations is a normalized standard deviation of 1 and a normalized correlation of 1. Models are KMBM3 (red circle), GFDL-ESM4 (green star; Krasting et al., 2018), CESM2 (light blue square; Danabasoglu, 2019), MPI-ESM1.2-LR (pink pentagon; Wieners et al., 2019), NorESM2-LM (brown plus; Seland et al., 2019), CMCC-ESM2 (orange triangle; Lovato et al., 2021), EC-Earth3-CC (blue diamond; EC-Earth Consortium, 2021) and IPSL-CM6A-LR (purple hexagon; Boucher et al., 2018). With the exception of KMBM3, the data for all models were obtained by mining the CMIP6 database (https://esgf-node.llnl.gov/search/cmip6/) using the following search terms: CMIP/si/historical/annual output. Data were accessed between 1–5 February 2021.

Seasonal succession in functional types follows the general progression of zonal maxima in diatoms preceding calcifiers by a few weeks (albeit, in separate zonal ranges, Fig. 17). 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 Ocean. In the lower 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.

Figure 17 Model year 2015 zonally averaged surface NPP and phytoplankton biomass.

4 Model assessment under climate change

In addition to the historical model forcings described earlier, from the year 2005 to 2300 the simulations are forced using increasing CO₂ and non-CO₂ greenhouse gas concentrations, projected changes to the fraction of the land surface devoted to agricultural uses (calculated to the year 2100 by Hurtt et al., 2011, and then held constant after) and changes in the direct effect of sulfate aerosols following “business-as-usual” RCP8.5 scenario (RCP8.5, Riahi et al., 2007; Meinshausen et al., 2011). The wind fields continue to be geostrophically adjusted to air temperature anomalies.

4.1 Historical changes (1964–2014)

We next explore model trends from the 1960s to the 2010s and compare to available data from this period. Significant changes are simulated to have already occurred in NPP and phytoplankton biomass, with most of the change over this period occurring since the 1980s (Fig. 18). Global NPP is simulated to have declined 1.4 % between 1964 and 2014, with the strongest declines in the tropics, northern midlatitudes and Southern Ocean due to increasing thermal stratification. In KMBM3, that decline is dominated by the simulated loss of diatom biomass (a 24.4 % decline) due to their high nutrient requirements, largely in the Southern Ocean (Fig. 19). Diazotrophs are simulated to have experienced a 2.0 % loss in biomass globally over this period. Calcifying phytoplankton and low-latitude non-calcifying phytoplankton are simulated to have experienced a 0.4 % and 1.5 % loss, respectively, in net biomass. Increases in calcifying phytoplankton are simulated to have occurred in the Arctic and in the Southern Ocean and decreases are simulated to have occurred in the middle latitudes, where LP have increased their biomass. Warming and a lengthening growing season (but increased stratification; see Arctic temperature trend in Fig. 20) in the Arctic benefit calcifiers. Expansion of coccolithophores into the Arctic has been observed over recent decades (Neukermans et al., 2018). In KMBM3, the Southern Ocean is also simulated to have experienced increasingly favourable conditions for calcifiers, with diatom biomass declining as calcifiers increase. Note, however, that this model does not resolve ocean acidification effects on calcification.

Figure 18 Modelled historical changes in zonally averaged and depth-integrated NPP, and phytoplankton biomass, from 1964 to 2014.

Figure 19 Modelled changes in phytoplankton biomass (mmol C m⁻³) from 1964 to 2014 (top row) and from 2014 to 2294 (bottom row).

KMBM3 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. KMBM3 also does not simulate a large decline prior to the 1950s (Fig. 21). Wernand et al. (2013) found no historical trend in chlorophyll in their Forel–Ule ocean colour scale proxy from 1899 to 2000. However, both chlorophyll reconstructions (and KMBM3) suggest strong regional variation in the historical trends. KMBM3 also simulates a smaller decline in global NPP than an earlier modelling effort, which obtained a 6.5 % decline between 1960–2006 (Laufkötter et al., 2013). Their simulation similarly resulted in large declines in production in the low latitudes, which they attributed to warming and stratification over the historical period. KMBM3 simulates strong warming in the low latitudes (up to 0.5^∘ warming in the zonal mean between 1964 and 2014; Fig. 20). Zonal mean trends in idealized natural radiocarbon (simulated without bomb or Suess effects, also Fig. 20) are also positive in the upper ocean, suggesting enhanced stratification and reduced vertical mixing, particularly in the tropical and subtropical Pacific and Indian Ocean basins. Laufkötter et al. (2013) also found strong changes in phytoplankton biogeography in the Southern Ocean and North Atlantic, which they attributed to increasing zonal wind stress and vertical mixing, and surface freshening inducing stratification, respectively. However, our models do not agree in the response of diatoms and high-latitude calcifiers, for which their model simulates high-latitude increases 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 functional types. However, the Southern Ocean in their model showed no clear trend in phytoplankton community composition between 1998–2012 (Rousseaux and Gregg, 2015). The model we present here is fully competition driven. Variable particle sinking and remineralization rates, and explicit CaCO₃ ballasting, further differentiate KMBM3. All of these factors may contribute to the different behaviour reported here.

Figure 20 Modelled change in temperature, radiocarbon and ideal age profiles by major ocean basins from 1964 to 2014.

Figure 21 Major biogeochemical fluxes in the model, 1776 to 2300. Panel (a) is total net primary production, (b) is total particle fluxes at 130 m depth, (c) is total particle fluxes at 2 km depth. Blue lines in the middle and right panels are PIC, green lines are POC, and red lines are opal. Solid lines are from this model (KMBM3); dashed lines are from the Keller et al. (2012) version (KMBM1).

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 increased in the Indian and Pacific sectors of the Southern Ocean and in the Arctic. Trends in deep ocean carbon particle export in KMBM3 generally follow trends in calcifier (CP) biomass (Fig. 19), as CaCO₃ (PIC) ballasting contributes significantly to deep export (Kvale et al., 2015a, b, 2019). Declining POC export in the North Atlantic is also simulated to have occurred, though in this region the trend is driven by an increase in diatoms, which are less efficient exporters of organic carbon. Deep ocean opal export is simulated to have declined almost globally, with the largest loss in the Southern Ocean. Opal dissolution 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.8 % (POC), 2.2 % (PIC) and 6.2 % (opal).

Unfortunately, there is no comparable historical reconstruction of deep particle fluxes over this time period. Just as with NPP, Laufkötter et al. (2013) simulated a larger decline (8 %) in export production (POC) between 1960 and 2006, which they calculated at 100 m depth. They also simulated the largest declines in the Indian Ocean and west-central Pacific basins, driven by declines in NPP, which in turn were driven by declines in nutrient availability owing to increasing stratification. The models disagree with respect to trends in the North Atlantic, with KMBM3 producing a decline in POC and PIC deep export (due to a decrease in calcifiers), and only a small increase in opal (due to an increase in diatoms). Both Laufkötter et al. (2013) and KMBM3 simulate a historical increase in diatoms in this region; the difference in export trends can be explained by which functional type (calcifiers or diatoms) is more efficient at POC export, with calcifiers being the more efficient carbon exporter in our formulation. Likewise, KMBM3 simulates different POC export trends in the Southern Ocean as Laufkötter et al. (2013), due to differences in model structure. Increasing calcifiers in the Indian sector in KMBM3 increase POC and PIC export there, while in Laufkötter et al. (2013) diatoms regionally increase.

Changes in carbon and nutrient profiles between 1964 and 2014 (Fig. 22) reflect a combination of physical–chemical uptake of anthropogenic CO₂ and the changes in ocean circulation on tracer accumulation (Figs. 20 and 24). Dissolved inorganic carbon (DIC) is simulated to have increased by more than 10 mmol C m⁻³ in zonal mean over this time period in both the Southern Ocean and North Atlantic. Both regions, but particularly the Southern Ocean, are thought to be the primary regions for anthropogenic carbon uptake into the ocean interior (e.g. Khatiwala et al., 2009). The low latitudes are primary regions of carbon storage (Frölicher et al., 2015), and strong (greater than 25 mmol C m⁻³) increases in DIC are also seen here, in all ocean basins, in the zonal mean profiles.

Figure 22 Modelled change in carbon and nutrient profiles by major ocean basins from 1964 to 2014.

Declines in phosphate and nitrate are simulated in the upper ocean in all basins, with the largest zonal mean declines of up to 0.1 mmol P m⁻³ and 0.5 mmol N m⁻³). Silicic acid increases by up to 20 mmol Si m⁻³ in the central North Pacific surface (Fig. 23) due to a regional shift from DT to CP dominance. A recent reconstruction suggests a positive silicate concentration trend in the upper 300 m in this region over the 1970–2017 time period (Stramma et al., 2020), although the estimated magnitude is lower (0.013±0.013 mmol Si m⁻³ yr⁻¹). Dissolved iron shows an increasing trend in the tropical surface despite using constant atmospheric sources. Surface and subsurface (300 m depth) concentration trends (Fig. 23) reveal spatial heterogeneity, with declining surface and increasing subsurface phosphate occurring in the subarctic North Pacific (a result of increasing particle export) and increasing concentrations along the Humboldt Current (driven by a regional reduction in nitrate; not shown). These results compare favourably with surface declines in phosphate recorded in the North Pacific between 1961 and 2012 (Yasunaka et al., 2016). They also compare favourably with the recent data compilation of Stramma et al. (2020), which shows a decline in surface nutrients, and an increase in subsurface nutrients, in the subarctic Pacific. In KMBM3, this pattern is produced by the replacement of diatoms with calcifiers, who are more efficient exporters of nutrients. The UVic ESCM lacks a fully dynamic atmosphere model and therefore does not simulate multi-decadal oscillations in climate, which have been implicated in recent Pacific interior nutrient trends (Stramma et al., 2020). This may be why KMBM3 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 KMBM3 does not simulate. This may be because KMBM3 does not include anthropogenic sources of nitrate (summarized by Stramma et al., 2020).

Figure 23 Modelled change in nutrients at the surface (a–c) and 300 m depth (d–e) from 1964 to 2014.

Significant increases in North Atlantic deep ocean concentrations (below 3000 m depth) occur in phosphate, nitrate and silicic acid (Fig. 22). A decline in maximum meridional overturning (MOC) of about 1 Sv is apparent over this time, and the water has warmed more than 0.3^∘ (Fig. 24). Warming of the Gulf Stream increases particle remineralization rates, thereby raising nutrient concentrations along the North Atlantic Deep Water pathway.

Figure 24 Major physical changes in the model from 1776 to 2300.

Southern Ocean-sourced intermediate and deep water is simulated to have increased nutrient concentrations from 1964 to 2014, with a positive trend in phosphate and nitrate outcropping along the Antarctic ice margin, in qualitative agreement with the limited observations of an increasing trend in the Indian sector from 1965–2008 (Iida et al., 2013). The Southern Ocean is simulated to have experienced subsurface 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 declines in NPP and export production upstream of the Southern Ocean 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; Fig. 20 shows a decline in ideal age (greater than 5 years in the zonal mean) in Pacific Intermediate water and the Indian Ocean, which leads to less particle remineralization (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 remineralization, resulting in higher nutrient concentrations.

Oxygen is simulated to have declined in all ocean basins, with the exception of the subsurface low latitudes. It has been previously estimated that the global ocean lost more than 2 % of its oxygen since the 1960s (Schmidtko et al., 2017). KMBM3 simulates a 0.7 % decline in total oxygen content from 1964–2014, which is an under-estimate 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 remineralization) contributes to the simulated decline in oxygen. In the upper ocean, warming has reduced oxygen solubility, lowering near-surface concentrations.

4.2 Long-term future changes (2014–2294)

More significant changes in ocean biogeochemistry are still to come, if applied boundary forcing assumptions hold over the next centuries. Figures 25 to 28 extend the previous analysis to the year 2294, with respect to 2014 biogeochemistry. Critically, spatial patterns in NPP trends reverse, with strong increases, in places exceeding 80 gC m⁻² yr⁻¹ in the zonal mean, in the Southern Ocean (Fig. 25). This trend is dominated by the increase in calcifiers (CP) and, to a lesser geographical extent, diatoms (DT) (Fig. 19). Calcifiers also continue their historical expansion into the Arctic (again, noting adverse physiological effects of ocean acidification are not simulated in this model), while diazotrophs are simulated to significantly increase in the middle latitudes after year 2100. Low-latitude phytoplankton (LP) also broadly increase in biomass between 20^∘ S and 60^∘ N, though the net NPP trend in the low latitudes remains negative.

Figure 25 Modelled changes in zonally averaged and depth-integrated NPP, and phytoplankton biomass, from 2014 to 2294.

Figure 26 Modelled change in temperature, radiocarbon and ideal age profiles by major ocean basins from 2014 to 2296.

This new ecosystem model responds differently to forcing than previous versions. Kvale et al. (2015a) compared the responses of the KMBM1 biogeochemistry (no calcifiers) to two versions of calcifier model, also integrated to 2300 using RCP8.5 forcing, and found the application of calcifiers eliminated the global reduction of NPP found in the Keller et al. (2012) version until around the year 2100. However, the introduction of diatoms and iron, and reorganization of the phytoplankton community structure, produces an even larger decline in global NPP to 2100 (close to 5 Pg C yr⁻¹) in KMBM3 than found in KMBM1 (less than 1 Pg C yr⁻¹). The difference appears to be in the low-latitude response, where calcifiers with low nutrient requirements in the previous version maintained NPP despite stratification, partly by supporting diazotroph nitrogen fixation through efficient ballast removal of surface nitrate (Kvale et al., 2019). In the current model, calcifiers have a lower biomass in the low latitudes and do not establish this symbiotic relationship with diazotrophs to the same extent. As a consequence, diazotrophs decline more strongly than previously. After 2100, NPP increases abruptly (Fig. 21). Rising NPP over the long-term is a long-standing feature of the KMBM biogeochemical model formulation and occurs because of the acceleration of nutrient recycling by the temperature-sensitive microbial loop in the low and middle latitudes (Schmittner et al., 2008). In this latest model version, this increase in NPP is much smaller (less than 4 Pg C yr⁻¹ by 2300, compared to 11–13 Pg C yr⁻¹ 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 attenuation, is difficult to assess, but it confirms the finding of Fu et al. (2016) of a strong dependence on phytoplankton community structure in model response to climate change.

Trends in phytoplankton biomass and productivity can be explained by the physical changes in the model over this time (Figs. 24 and 26). Maximum MOC declines from 18 to 11 Sv over the period 1800 to 2200 before starting to increase after 2200. Slowing overturning helps to accelerate surface warming, and zonal mean temperatures in the upper North Atlantic rise more than 5^∘ by 2294. Globally increased radiocarbon in the upper ocean suggests widespread increased stratification, as well as a more complete separation between upper and lower water masses globally (as previously reported by Kvale et al., 2018, with only slightly different forcing conditions). Ideal age trends similarly show the lengthening Southern Ocean-sourced deep water pathway that extends to ventilate all ocean basins from the south, replacing the shoaled northern ventilation pathways (also described in detail in Kvale et al., 2018). Ideal ages increase more than 300 years in the deep southern Pacific and Atlantic basins, and in the deep North Atlantic, where North Atlantic Deep 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 functional types (DT, CP) thrive, while in the lower latitudes it is the more efficient nutrient consumers (DZ and LP) which benefit.

Long-term particle export trends generally follow the historical trend but with increasing magnitude (Figs. 21 and 27). 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 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 KMBM3 is larger than in previous versions (about 2.0 Pg C yr⁻¹ by 2100, rather than 1.5 Pg C yr⁻¹).

Figure 27 Modelled change in particle fluxes at 2 km depth, from 2014 to 2294.

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 (Fig. 27). Strong decreases in calcifiers, and associated deep carbon export, occur in the Indian Ocean and North Atlantic. Opal export declines by more than 100 mmol Si m⁻² yr⁻¹ 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 (Fig. 28), particularly iron, and a short growing season favours diatoms over calcifiers.

Figure 28 Modelled change in carbon and nutrient profiles by major ocean basins from 2014 to 2294.

The historical trend in carbon and nutrients is similarly extended, with continuing increases in DIC in the upper ocean (as atmospheric CO₂ 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 (Fig. 28). Oxygen declines along the Southern Ocean–abyssal global ocean pathway due to both warming and increasing particle remineralization, which is also responsible for the increasing nutrient concentrations in the deep ocean. Decreasing phosphate and nitrate concentrations in the subsurface tropical ocean basins are a product of declining particle remineralization there, brought about by both warming, which shoals remineralization 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⁻³ in the zonal mean) concentrations in the upper ocean was previously described by Nickelsen et al. (2015). They attributed the increase to stratification “trapping” aerosol iron near the surface. However, the regions showing the greatest increases in 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 in Southern Ocean primary production that couples with a reorganization of deep ocean circulation to produce a long term “nutrient trapping” effect in Southern Ocean-sourced interior water masses (e.g. Moore et al., 2018; Kvale et al., 2019). Near-surface increases in iron, and decreases in nitrate, phosphate and silicic acid, have also been observed to 2100 in a comparison of nine 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 (unlike in KMBM3). Phytoplankton community composition and export formulation were 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.

5 Conclusions

Our paper describes a new model of the marine silicate cycle (KMBM3), evaluates its performance against previous KMBM versions and other Earth system models, 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 with respect to previous versions of KMBM and is mechanistically more realistic, with the added complexity of iron, calcite and silicate merged into a single model code. Furthermore, its representation of the silicate cycle has similar performance with other state-of-the-art biogeochemical models.

Simulations using our new model suggest diatoms have been, and will continue to be, the losers as the Earth system warms. 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 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.

Several aspects of KMBM3, 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 optimization 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 look forward to further refinements and the many applications of this model to come.