Ocean biogeochemistry in the Canadian Earth System Model version 5.0.3: CanESM5 and CanESM5-CanOE

The ocean biogeochemistry components of two new versions of the Canadian Earth System Model (CanESM) are presented and compared to observations and other models. CanESM5 employs the same ocean biology model as CanESM2, whereas CanESM5-CanOE (Canadian Ocean Ecosystem model) is a new, more complex model developed for CMIP6, with multiple food chains, flexible phytoplankton elemental ratios, and a prognostic iron cycle. This new model is described in detail and the outputs (distributions of major tracers such as oxygen, dissolved inorganic carbon, and alkalinity, the iron and nitrogen cycles, plankton biomass, and historical trends in CO₂ uptake and export production) compared to CanESM5 and CanESM2, as well as to observations and other CMIP6 models. Both CanESM5 models show gains in skill relative to CanESM2, which are attributed primarily to improvements in ocean circulation. CanESM5-CanOE shows improved skill relative to CanESM5 for most major tracers at most depths. CanESM5-CanOE includes a prognostic iron cycle, and maintains high-nutrient/low-chlorophyll conditions in the expected regions (in CanESM2 and CanESM5, iron limitation is specified as a temporally static “mask”). Surface nitrate concentrations are biased low in the subarctic Pacific and equatorial Pacific, and high in the Southern Ocean, in both CanESM5 and CanESM5-CanOE. Export production in CanESM5-CanOE is among the lowest for CMIP6 models; in CanESM5, it is among the highest, but shows the most rapid decline after about 1980. CanESM5-CanOE shows some ability to simulate aspects of plankton community structure that a single-species model can not (e.g. seasonal dominance of large cells) but is biased towards low concentrations of zooplankton and detritus relative to phytoplankton. Cumulative ocean uptake of anthropogenic carbon dioxide through 2014 is lower in both CanESM5-CanOE (122 PgC) and CanESM5 (132 PgC) than in observation-based estimates (145 PgC) or the model ensemble mean (144 PgC).

1 Introduction

The Canadian Centre for Climate Modelling and Analysis has been developing coupled models with an interactive carbon cycle for more than a decade (Christian et al., 2010; Arora et al., 2011). The Canadian Earth System Model version 5 (CanESM5, Swart et al., 2019a) is an updated version of CanESM2 (Arora et al., 2011), with a new ocean model based on the Nucleus for European Modelling of the Ocean (NEMO) system version 3.4. The ocean biogeochemistry modules were developed in house. CanESM5 uses the same ocean biology model as CanESM1 (Christian et al., 2010) and CanESM2 (Arora et al., 2011), the Canadian Model of Ocean Carbon (CMOC; Zahariev et al., 2008). An additional model was developed for CMIP6, called the Canadian Ocean Ecosystem model (CanOE). The biological components of CanOE are of substantially greater complexity than CMOC, including multiple food chains, flexible phytoplankton elemental ratios, and a prognostic iron (Fe) cycle. The two coupled models are known as CanESM5 and CanESM5-CanOE, respectively.

The reasons for developing both models are, firstly, to evaluate the effect of changes in ocean circulation between CanESM2 and CanESM5 on ocean biogeochemistry by running the new climate model with the same ocean biogeochemistry, and secondly because CanOE is substantially more expensive computationally (19 tracers vs. 7, so the total computation time to integrate the ocean model with biogeochemistry is approximately double). Most CMIP6 experiments were run with CanESM5 only, as ocean biogeochemistry is not central to their purpose. Additional tracers requested by the Ocean Model Intercomparison Project – Biogeochemistry (OMIP-BGC) including abiotic and natural dissolved inorganic carbon (DIC), DI¹⁴C, CFCs, and SF₆ (see Orr et al., 2017) were run only in CanESM5. The CMIP6 experiments published for CanESM5-CanOE are listed in Table S1 in the Supplement.

CMOC is a nutrient–phytoplankton–zooplankton–detritus (NPZD) model with highly parameterized representations of phytoplankton Fe limitation, dinitrogen (N₂) fixation and denitrification, and calcification and calcite dissolution (Zahariev et al., 2008; Fig. S1 in the Supplement). CanESM1 and CanESM2 did not include oxygen; CanESM5 includes oxygen as a purely “downstream” tracer that does not affect other biogeochemical processes. In CanESM5-CanOE, denitrification is prognostic and dependent on the concentration of oxygen. Among the less satisfactory aspects of CMOC biogeochemistry are, firstly, that Fe limitation is specified as a static “mask” that does not change with climate (it is calculated from the present-day climatological distribution of nitrate, based on the assumption that regions without iron limitation will have complete drawdown of surface nitrate at some point in the year), and secondly, that denitrification is parameterized so that nitrogen (N) is conserved within each vertical column, i.e. collocated with N₂ fixation in tropical and subtropical open-ocean regions (Zahariev et al., 2008; Riche and Christian, 2018). This latter simplification produced excessive accumulations of nitrate in Eastern Boundary Current (EBC) regions where most denitrification occurs. CMOC also has a tendency to produce rather stark extremes of high and low primary and export production (Zahariev et al., 2008), a well-known problem of NPZD models (Armstrong, 1994; Friedrichs et al., 2007). Our intent in developing CanOE was to alleviate, or at least reduce, these biases, by including multiple food chains, a prognostic Fe cycle, and prognostic denitrification. Dinitrogen fixation is still parameterized, but the CanOE parameterization includes Fe (but not P) limitation, whereas in CMOC N₂ fixation tends to grow without bound in a warming ocean as CMOC does not include P or Fe limitation (Riche and Christian, 2018).

In this paper, we present a detailed model description for CanOE and an evaluation of both CanESM5 and CanESM5-CanOE relative to observational data products and other available models. CMOC has been well described previously (Zahariev et al., 2008) and the details are not reiterated here. In some cases, CanESM2 results are also shown to illustrate which differences in the model solutions arise largely from the evolution of the physical climate model, and which are specifically associated with different representations of biogeochemistry. An overall evaluation of the CanESM5 climate including the physical ocean is given in Swart et al. (2019a). Here, we focus on biogeochemical variables, and we have evaluated model performance in three main areas: (1) the distribution of major tracers like oxygen, DIC, and alkalinity, and the resulting saturation state for CaCO₃ minerals, (2) the iron cycle and its interaction with the nitrogen cycle, and (3) plankton community structure and the concentration and export of particulates. We first address the major chemical species that are common to both models (and almost all other Earth system models) to determine whether a more complex biology model measurably improves skill, and whether the updated circulation model improves skill relative to CanESM2. Then we examine the areas where our two models differ: the presence of a prognostic iron cycle and multiple food chains in CanOE. More specifically, does CanESM5-CanOE reproduce the geographic distribution of high-nutrient/low-chlorophyll (HNLC) regions? Does the large phytoplankton/large zooplankton food chain become dominant under nutrient-rich conditions, and how does having multiple detrital size classes affect particle flux and remineralization length scale? Following this model evaluation, we present historical trends in ocean anthropogenic CO₂ uptake, export production, and total volume of low-oxygen waters over the historical (1850–2014) experiment. Possible future changes under Shared Socioeconomic Pathway experiments will be addressed in subsequent publications.

2 Model description

CanESM5 (Swart et al., 2019a) is an updated version of CanESM2 (Arora et al., 2011), with an entirely new ocean. The atmosphere model has the same T63 horizontal resolution, and contains some important improvements in atmospheric physics (Swart et al., 2019a). The land surface (Canadian Land Surface Scheme) and terrestrial carbon cycle (Canadian Terrestrial Ecosystem Model) models are substantially the same as in CanESM2 with minor modifications as described by Arora et al. (2020). The CanESM5 ocean is based on the NEMO modelling system version 3.4, with a horizontal resolution of 1^∘, telescoping to 1/3^∘ in the tropics, and 45 vertical levels ranging in thickness from ∼ 6 m near the surface to ∼ 250 m in the deep ocean (Swart et al., 2019a). All physical climate model components are the same in CanESM5 and CanESM5-CanOE. There are no feedbacks between biology and the physical ocean model, so the physical climate of CanESM5 and CanESM5-CanOE is identical in experiments with prescribed atmospheric CO₂ concentration.

The NEMO system is a publicly available archive of codes based on the OPA (Océan PArallelisé) ocean model (Madec and Imbard, 1996; Guilyardi and Madec, 1997) and the Tracers in Ocean Paradigm (TOP) module for tracer advection and mixing. Our ocean biogeochemistry modules are built within TOP, using NEMO v3.4.1, but have also been implemented in NEMO 3.6 for regional downscaling applications (Holdsworth et al., 2021).

Carbon chemistry is based on the Best Practices Guide (Dickson et al., 2007) and the OMIP-BGC data request (Orr et al., 2017) and is identical in CanESM5 and CanESM5-CanOE. All calculations are done on the total scale and the recommended formulae for the equilibrium constants are employed. The carbon chemistry solver was run for a fixed number of iterations (10 in the surface layer and 5 in the subsurface layers in CanESM5-CanOE). CanESM5 does not solve the carbon chemistry equations in the subsurface layers. OMIP-BGC formulations for CO₂ and O₂ solubility and gas exchange are employed. It is important to note here that the carbon chemistry and gas exchange formulations used in CanESM2 (and other CMIP5 models) are slightly different than those used in CMIP6. However, this difference is of little functional significance; i.e. it will have a negligible impact on the distribution of [CO${}_{\mathrm{3}}^{--}$] compared to the differences in DIC and alkalinity distribution. The initialization fields for nitrate, DIC and alkalinity were also different in CanESM2. This will affect the total ocean inventory of DIC but not the spatial distribution if the model is well equilibrated.

The CanOE biology model is based on the cellular regulation model of Geider et al. (1998). There are two phytoplankton size classes, and each group has four state variables: C, N, Fe, and chlorophyll. Photosynthesis is decoupled from cell production and photosynthetic rate is a function of the cell's internal N and Fe quotas. Each functional group has a specified minimum and maximum N quota and Fe quota, and nutrient uptake ceases when the maximal cell quota is reached. Chlorophyll synthesis is a function of N uptake and increases at low irradiance. There are also two size classes each of zooplankton and detritus. Small zooplankton graze on small phytoplankton, while large zooplankton graze on both large phytoplankton and small zooplankton. Small detritus sinks at 2 m d⁻¹ and large detritus at 30 m d⁻¹ (in CanESM5, there is a single detrital pool with a sinking rate of 8 m d⁻¹). Model parameters and their values are listed in Table 1. A schematic of the model is shown in Fig. 1.

Figure 1 Schematic of the CanOE biology model. Model currencies including chlorophyll (Chl) are indicated by coloured boxes except oxygen (O₂) and carbonate (CaCO₃). Arrows indicate flows of carbon (C), nitrogen (N), and iron (Fe) between compartments containing small (S) and large (L) phytoplankton (P), zooplankton (Z), and detritus (D) components; counterflows of oxygen are not shown.

Table 1 Ecosystem model parameters.

2.1 Photosynthesis and phytoplankton growth

For simplicity and clarity, the equations are shown here for a single phytoplankton species and do not differ structurally for small and large phytoplankton. Some parameter values differ for the two phytoplankton groups; all parameter values are listed in Table 1.

Temperature dependence of photosynthetic activity is expressed by the Arrhenius equation:

$$\begin{array}{}\text{(1)}& T_{\mathrm{f}}=\exp \left(-{\displaystyle \frac{E_{\mathrm{ap}}}{R}}\left({\displaystyle \frac{\mathrm{1}}{T}}-{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{ref}}}}\right)\right),\end{array}$$

where E_ap is an enzyme activation energy that corresponds approximately to that of RuBisCo (see Raven and Geider, 1988), R is the gas constant (8.314 J mol⁻¹ K⁻¹), and temperature T and reference temperature T_ref are in Kelvin. Maximal rates of nutrient (either N or Fe but generically referred to here with the superscript X) uptake are given by

$$\begin{array}{}\text{(2)}& V_{\max }^X=V_{\mathrm{ref}}^X{T}_{\mathrm{f}}{\left({\displaystyle \frac{Q_{\max }^X-Q^X}{Q_{\max }^X-Q_{\min }^X}}\right)}^{\mathrm{0.05}},\end{array}$$

where $V_{\max }^X$ is the maximal uptake rate in mg of nutrient X per mg of cell C, X can represent N or Fe, Q is the nutrient cell quota and Q_min and Q_max its minimum and maximum values, and $V_{\mathrm{ref}}^X$ is a (specified) basal rate at T=T_ref and Q=Q_min. These maximum rates are then reduced according to the ambient nutrient concentration, i.e.

$$\begin{array}{}\text{(3a)}& V^{\mathrm{N}}=V_{\max }^{\mathrm{N}}\left(L_{{\mathrm{NH}}_{\mathrm{4}}}+\left(\mathrm{1}-L_{{\mathrm{NH}}_{\mathrm{4}}}\right)L_{{\mathrm{NO}}_{\mathrm{3}}}\right),\end{array}$$

where $L_{\mathrm{NH}\mathrm{4}}=\frac{N_{\mathrm{a}}}{K_{\mathrm{Na}X}+N_{\mathrm{a}}}$ and $L_{\mathrm{NO}\mathrm{3}}=\frac{N_{\mathrm{i}}}{K_{\mathrm{Ni}X}+{\mathrm{N}}_{\mathrm{i}}}$, with N_i and N_a indicating nitrate and ammonium respectively, and

$$\begin{array}{}\text{(3b)}& V^{\mathrm{Fe}}=V_{\max }^{\mathrm{Fe}}\left({\displaystyle \frac{\mathrm{Fe}}{K_{\mathrm{Fe}X}+\mathrm{Fe}}},\right)\end{array}$$

where X indicates large or small phytoplankton (Table 1). The maximal carbon-based growth rate is given by

$$\begin{array}{}\text{(4)}& {\mathrm{P}}_{\max }^C={\mathrm{P}}_{\mathrm{ref}}^C{T}_{\mathrm{f}}\min \left\{{\displaystyle \frac{Q^{\mathrm{N}}-Q_{\min }^{\mathrm{N}}}{Q_{\max }^{\mathrm{N}}-Q_{\min }^{\mathrm{N}}}},{\displaystyle \frac{Q^{\mathrm{Fe}}-Q_{\min }^{\mathrm{Fe}}}{Q_{\max }^{\mathrm{Fe}}-Q_{\min }^{\mathrm{Fe}}}}\right\},\end{array}$$

where $P_{\mathrm{ref}}^{\mathrm{C}}$ is the rate at the reference temperature T_ref under nutrient-replete conditions (Q=Q_max). The light-limited growth rate is then given by

$$\begin{array}{}\text{(5)}& {\mathrm{P}}_{\mathrm{phot}}^C={\mathrm{P}}_{\max }^C\left(\mathrm{1}-e^{-{\mathit{\alpha }}_{\mathrm{chl}}E{\mathit{\theta }}_{\mathrm{C}}/{\mathrm{P}}_{\max }^C}\right),\end{array}$$

where E is irradiance and θ_C is the chlorophyll-to-carbon ratio. The rate of chlorophyll synthesis is

$$\begin{array}{}\text{(6)}& {\mathit{\rho }}_{\mathrm{chl}}={\mathit{\theta }}_{\max }^{\mathrm{N}}{\displaystyle \frac{{\mathrm{P}}_{\mathrm{phot}}^C}{E{\mathit{\alpha }}_{\mathrm{chl}}\mathit{\theta }}}.\end{array}$$

These rates are then used to define a set of state equations for phytoplankton carbon (C_p), nitrogen (N_p), iron (Fe_p), and chlorophyll (M).

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\displaystyle \frac{{\mathrm{dC}}_{\mathrm{p}}}{\mathrm{d}t}}& =({\mathrm{P}}_{\mathrm{phot}}^C-\mathit{\zeta }{V}_{\mathrm{N}}){\mathrm{C}}_{\mathrm{p}}-(G+{\mathrm{C}}_{\mathrm{XS}})-m_{\mathrm{1}}{\mathrm{C}}_{\mathrm{p}}\\ & -m_{\mathrm{2}}{\mathrm{C}}_{\mathrm{p}}^{\mathrm{2}}-k_{\mathrm{XU}}{\mathrm{C}}_{\mathrm{INTR}},\end{array}\end{array}$$

where ζ is the respiratory cost of biosynthesis, G is the grazing rate (Eq. 12), C_XS is the excess (above the ratio in grazer biomass) carbon in grazing losses (see Eq. 16a below), m₁ and m₂ are coefficients for linear and quadratic nonspecific mortality terms, C_INTR is the concentration of intracellular carbohydrate carbon in excess of biosynthetic requirements, and k_XU is a rate coefficient for its exudation to the environment. The nonspecific mortality terms are set to 0 below 0.01 mmol C m⁻³ to prevent biomass from being driven to excessively low levels in the high latitudes in winter; linear mortality terms can result in biomass declining to levels from which recovery would take much longer than the brief Arctic summer (Hayashida, 2018). The full equations for phytoplankton N, Fe, and chlorophyll are

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{{\mathrm{dN}}_{\mathrm{p}}}{\mathrm{d}t}}={\displaystyle \frac{V^{\mathrm{N}}}{Q_{\mathrm{N}}}}-(G+m_{\mathrm{1}}{\mathrm{C}}_{\mathrm{p}}+m_{\mathrm{2}}{\mathrm{C}}_{\mathrm{p}}^{\mathrm{2}})R_{\mathrm{NC}}-{\mathrm{N}}_{\mathrm{XS}}\text{(9)}& {\displaystyle \frac{{\mathrm{dFe}}_{\mathrm{p}}}{\mathrm{d}t}}={\displaystyle \frac{V^{\mathrm{Fe}}}{Q_{\mathrm{Fe}}}}-(G+m_{\mathrm{1}}{\mathrm{C}}_{\mathrm{P}}+m_{\mathrm{2}}{\mathrm{C}}_{\mathrm{p}}^{\mathrm{2}})R_{\mathrm{FeC}}-{\mathrm{Fe}}_{\mathrm{XS}}\text{(10)}& {\displaystyle \frac{\mathrm{dM}}{\mathrm{d}t}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{chl}}{V}^{\mathrm{N}}}{{\mathit{\theta }}_{\mathrm{C}}}}M-(G+m_{\mathrm{1}}{\mathrm{C}}_{\mathrm{p}}+m_{\mathrm{2}}{\mathrm{C}}_{\mathrm{p}}^{\mathrm{2}}){\mathit{\theta }}_{\mathrm{C}}-k_{\mathrm{dgr}}M,\end{array}$$

where k_dgr is a rate coefficient for nonspecific losses of chlorophyll e.g. by photooxidation, in addition to losses to grazing and other processes that also affect C_p, N_p, and Fe_p. N_XS and Fe_XS are remineralization of “excess” (relative to grazer or detritus ratios) N or Fe and are defined below (Eq. 16).

2.2 Grazing and food web interactions

Grazing rate depends on the phytoplankton carbon concentration, which most closely represents the food concentration available to the grazer (Elser and Urabe 1999; Loladze et al. 2000). Zooplankton biomass is also in carbon units. State equations for small and large zooplankton are

$$\begin{array}{}\text{(11a)}& {\displaystyle \frac{\mathrm{d}{Z}_s}{\mathrm{d}t}}=\mathit{\lambda }{G}_s-(R+G_{\mathrm{Z}}+m_{\mathrm{1}s}{Z}_s+m_{\mathrm{2}s}{Z}_s^{\mathrm{2}})\end{array}$$

$$\begin{array}{}\text{(11b)}& {\displaystyle \frac{\mathrm{d}{Z}_{\mathrm{L}}}{\mathrm{d}t}}=\mathit{\lambda }{G}_{\mathrm{L}}-(R+m_{\mathrm{1}\mathrm{L}}{Z}_{\mathrm{L}}+m_{\mathrm{2}\mathrm{L}}{Z}_{\mathrm{L}}^{\mathrm{2}}),\end{array}$$

where

$$\begin{array}{}\text{(12a)}& G_s=G_{\mathrm{so}}(\mathrm{1}-e^{-a_s{\mathrm{C}}_{\mathrm{ps}}})Z_s\end{array}$$

$$\begin{array}{}\text{(12b)}& G_{\mathrm{L}}=G_{\mathrm{L}\mathrm{0}}(\mathrm{1}-e^{-a_l({\mathrm{C}}_{\mathrm{pl}}+Z_s)})Z_{\mathrm{L}},\end{array}$$

for small and large zooplankton, respectively, G_Z is grazing of small zooplankton by large zooplankton, R is respiration, and m₁ and m₂ are non-grazing mortality rates. Large zooplankton grazing is divided into grazing on large phytoplankton and small zooplankton in proportion to the relative abundance of each:

$$\begin{array}{}\text{(13a)}& G_{\mathrm{P}}=G_{\mathrm{L}}{\displaystyle \frac{{\mathrm{P}}_l}{{\mathrm{P}}_l+Z_s}}\end{array}$$

$$\begin{array}{}\text{(13b)}& G_{\mathrm{Z}}=G_{\mathrm{L}}{\displaystyle \frac{Z_s}{{\mathrm{P}}_l+Z_s}}.\end{array}$$

Zooplankton biomass loss to respiration is given by

$$\begin{array}{}\text{(14)}& R=\max \left\{r_z{T}_{\mathrm{f}}Z-{\mathrm{C}}_{\mathrm{XS}},\mathrm{0}\right\},\end{array}$$

and uses the same activation energy as photosynthesis. Respiration (R) is assumed to consume only carbon and not result in catabolism of existing biomass when “excess” carbon is available in the prey. In addition, conservation of mass must be maintained by recycling to the dissolved pool grazer consumption of elements in excess of biosynthetic requirements when grazer and prey elemental ratios differ. In the case where the nutrient quota (relative to carbon) exceeds the grazer fixed ratio, the excess nutrient is remineralized to the dissolved inorganic pool. In the case where the nutrient quota is less than the grazer ratio, the grazer intake is reduced to what can be supported by the least abundant nutrient (relative to the grazer biomass ratio) and excess carbon is remineralized. For the case of two nutrients (in this case N and Fe), it is necessary to define

$$\begin{array}{}\text{(15)}& G^{\prime }=G\min \left\{{\displaystyle \frac{{\mathrm{N}}_{\mathrm{P}}}{{\mathrm{C}}_{\mathrm{P}}}}{R}_{\mathrm{CN}},{\displaystyle \frac{{\mathrm{Fe}}_{\mathrm{P}}}{{\mathrm{C}}_{\mathrm{P}}}}{R}_{\mathrm{CFe}},\mathrm{1}\right\},\end{array}$$

where G is equal to G_S (Eq. 12a) for small zooplankton and G_P (Eq. 13a) for large zooplankton, and R_XY indicates the fixed ratio of element X to element Y in grazer biomass. The “excess” carbon available for respiration is

$$\begin{array}{}\text{(16a)}& {\mathrm{C}}_{\mathrm{XS}}=G{}^{\prime }max\left\{{\displaystyle \frac{{\mathrm{C}}_{\mathrm{P}}}{{\mathrm{N}}_{\mathrm{P}}}}{R}_{\mathrm{NC}}-\mathrm{1},{\displaystyle \frac{{\mathrm{C}}_{\mathrm{P}}}{{\mathrm{Fe}}_{\mathrm{P}}}}{R}_{\mathrm{FeC}}-\mathrm{1},\mathrm{0}\right\},\end{array}$$

and the excess nutrients remineralized to their inorganic pools are

$$\begin{array}{}\text{(16b)}& \begin{array}{rl}{\mathrm{N}}_{\mathrm{XS}}& =G{}^{\prime }\max \left\{{\displaystyle \frac{{\mathrm{N}}_{\mathrm{P}}}{{\mathrm{C}}_{\mathrm{P}}}}-R_{\mathrm{NC}},\mathrm{0}\right\}\mathit{\epsilon }\\ & +G{}^{\prime }\max \left\{R_{\mathrm{NC}}\left({\displaystyle \frac{{\mathrm{N}}_{\mathrm{P}}}{{\mathrm{Fe}}_{\mathrm{P}}}}{R}_{\mathrm{FeN}}-\mathrm{1}\right),\mathrm{0}\right\}\left(\mathrm{1}-\mathit{\epsilon }\right)\end{array}\end{array}$$

$$\begin{array}{}\text{(16c)}& \begin{array}{rl}{\mathrm{Fe}}_{\mathrm{XS}}& =G{}^{\prime }\max \left\{{\displaystyle \frac{{\mathrm{Fe}}_{\mathrm{P}}}{{\mathrm{C}}_{\mathrm{P}}}}-R_{\mathrm{FeC}},\mathrm{0}\right\}\mathit{\epsilon }\\ & +G{}^{\prime }\max \left\{R_{\mathrm{FeC}}({\displaystyle \frac{{\mathrm{Fe}}_{\mathrm{P}}}{{\mathrm{N}}_{\mathrm{P}}}}{R}_{\mathrm{NFe}}-\mathrm{1}),\mathrm{0}\right\}(\mathrm{1}-\mathit{\epsilon }),\end{array}\end{array}$$

where $\mathit{\epsilon }=\frac{\max \left\{{\mathrm{C}}_{\mathrm{xs}},\mathrm{0}\right\}}{{\mathrm{C}}_{\mathrm{xs}}+\mathrm{\Delta }}$ is a switch to prevent double counting in cases where one of the terms is redundant (the excess relative to the least abundant element is included in the other term) but would otherwise be nonzero (Δ is a constant equal to 10⁻¹⁵, to prevent division by zero). For three elements, there are 3! = 6 possible cases: for N greater or less than C_PR_NC, Fe may be either in excess relative to both C and N, deficient relative to both, or in excess relative to one but not the other (Table 2).

Table 2 Cases where the “excess” terms are nonzero. These terms are always greater than or equal to zero, and always zero when the phytoplankton elemental ratio is equal to the grazer biomass ratio. A plus (+) sign indicates that a specific term is positive. N₁ and N₂, Fe₁ and Fe₂ indicate the first and second terms in Eq. (16b) and (16c). R_NC is the grazer $\mathrm{N}/\mathrm{C}$ (Redfield) ratio.

2.3 Organic and inorganic pools

There are two pools of detritus with different sinking rates but the same fixed elemental ratios. Detrital C/N/Fe ratios are the same as zooplankton, so zooplankton mortality or grazing of small zooplankton by large zooplankton produce no “excess”. Phytoplankton mortality, and defecation by zooplankton grazing on phytoplankton, produces excess nutrient or excess C that needs to be recycled into the inorganic pool in a similar fashion as outlined above for the assimilated fraction of grazing on phytoplankton.

The conservation equations for detrital C are

$$\begin{array}{}\text{(17a)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{D}_s}{\mathrm{d}t}}& =m_{\mathrm{1}}({\mathrm{C}}_{\mathrm{ps}}+Z_s)+m_{\mathrm{2}}({\mathrm{C}}_{\mathrm{ps}}^{\mathrm{2}}+Z_{\mathrm{S}}^{\mathrm{2}})\\ & -r_{\mathrm{1}}{D}_s{T}_g-w_s{\displaystyle \frac{\mathrm{d}{D}_s}{\mathrm{d}z}}\end{array}\end{array}$$

$$\begin{array}{}\text{(17b)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{D}_l}{\mathrm{d}t}}& =m_{\mathrm{1}}({\mathrm{C}}_{\mathrm{pl}}+Z_{\mathrm{L}})+m_{\mathrm{2}}({\mathrm{C}}_{\mathrm{pl}}^{\mathrm{2}}+Z_{\mathrm{L}}^{\mathrm{2}})\\ & -r_{\mathrm{2}}{D}_l{T}_g-w_l{\displaystyle \frac{\mathrm{d}{D}_l}{\mathrm{d}z}},\end{array}\end{array}$$

where T_g is an Arrhenius function for temperature dependence of remineralization and w is the sinking speed. The conservation equations for inorganic C, N, and Fe are

$$\begin{array}{}\text{(18a)}& {\displaystyle \frac{{\mathrm{dC}}_{\mathrm{i}}}{\mathrm{d}t}}=(\mathit{\zeta }{V}^{\mathrm{N}}-{\mathrm{P}}_{\mathrm{phot}}^{\mathrm{C}}){\mathrm{C}}_{\mathrm{p}}+R+{\mathrm{C}}_{\mathrm{XS}}+(r_{\mathrm{1}}{D}_s+r_{\mathrm{2}}{D}_l)T_g\end{array}$$

$$\begin{array}{}\text{(18b)}& \begin{array}{rl}{\displaystyle \frac{{\mathrm{dN}}_{\mathrm{i}}}{\mathrm{d}t}}& =-{\displaystyle \frac{V^{\mathrm{N}}}{Q^{\mathrm{N}}}}{N}_{\mathrm{p}}\left({\displaystyle \frac{L_{\mathrm{NO}\mathrm{3}}}{L_{\mathrm{NO}\mathrm{3}}+L_{\mathrm{NH}\mathrm{4}}}}\right)\\ & +{\mathrm{N}}_{\mathrm{ox}}-{\mathrm{N}}_{\mathrm{dentr}}(\mathrm{1}-A_{\mathrm{f}})\end{array}\end{array}$$

$$\begin{array}{}\text{(18c)}& \begin{array}{rl}{\displaystyle \frac{{\mathrm{dN}}_{\mathrm{a}}}{\mathrm{d}t}}& =-{\displaystyle \frac{V^{\mathrm{N}}}{Q^{\mathrm{N}}}}{N}_{\mathrm{p}}\left({\displaystyle \frac{L_{\mathrm{NH}\mathrm{4}}}{L_{\mathrm{NO}\mathrm{3}}+L_{\mathrm{NH}\mathrm{4}}}}\right)+{\displaystyle \frac{R}{R_{\mathrm{CN}}}}+{\mathrm{N}}_{\mathrm{XS}}\\ & +(r_{\mathrm{1}}{D}_s+r_{\mathrm{2}}{D}_l)R_{\mathrm{NC}}{T}_g-{\mathrm{N}}_{\mathrm{ox}}+{\mathrm{N}}_{\mathrm{dnf}}\\ & -{\mathrm{N}}_{\mathrm{dentr}}{A}_{\mathrm{f}}\end{array}\end{array}$$

$$\begin{array}{}\text{(18d)}& {\displaystyle \frac{\mathrm{dFe}}{\mathrm{d}t}}={\displaystyle \frac{V^{\mathrm{Fe}}}{Q^{\mathrm{Fe}}}}{\mathrm{Fe}}_{\mathrm{p}}+{\displaystyle \frac{R}{R_{\mathrm{CFe}}}}+{\mathrm{Fe}}_{\mathrm{XS}}+(r_{\mathrm{1}}{D}_s+r_{\mathrm{2}}{D}_l)R_{\mathrm{FeC}}{T}_g,\end{array}$$

where N_ox is microbial oxidation of ammonium to nitrate (nitrification), N_dnf and N_dentr are sources and sinks associated with dinitrogen fixation and denitrification, and A_f is the ammonium fraction of denitrification losses, associated with anaerobic ammonium oxidation (“anammox”). The oxygen equation is essentially the inverse of Eq. (18a), with additional terms for oxidation and reduction of N, i.e.

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{{\mathrm{dO}}_{\mathrm{2}}}{\mathrm{d}t}}=-{\displaystyle \frac{{\mathrm{dC}}_{\mathrm{i}}}{\mathrm{d}t}}+\mathrm{2}{\displaystyle \frac{V^{\mathrm{N}}}{Q^{\mathrm{N}}}}{\mathrm{N}}_{\mathrm{p}}\left({\displaystyle \frac{L_{\mathrm{NO}\mathrm{3}}}{L_{\mathrm{NO}\mathrm{3}}+L_{\mathrm{NH}\mathrm{4}}}}\right)-\mathrm{2}{\mathrm{N}}_{\mathrm{ox}}.\end{array}$$

Nitrification is given by

$$\begin{array}{}\text{(20)}& {\mathrm{N}}_{\mathrm{ox}}=k_{\mathrm{NH}\mathrm{4}\mathrm{ox}}{\mathrm{N}}_{\mathrm{a}}{\displaystyle \frac{K_E}{K_E+E\left(z\right)}},\end{array}$$

where E(z) is the layer mean irradiance at depth z. Dinitrogen fixation is parameterized as an external input of ammonium dependent on light, temperature and Fe availability, and inhibited by high ambient concentrations of inorganic N:

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{\mathrm{N}}_{\mathrm{dnf}}=k_{\mathrm{dnf}}{T}_{\mathrm{dnf}}\left(\mathrm{1}-e^{-\mathrm{aE}}\right)& \left({\displaystyle \frac{\mathrm{Fe}}{K_{\mathrm{Fe}}+\mathrm{Fe}}}\right)\\ & \left({\displaystyle \frac{K_{\mathrm{NO}\mathrm{3}}}{K_{\mathrm{NO}\mathrm{3}}+{\mathrm{N}}_{\mathrm{i}}+{\mathrm{N}}_{\mathrm{a}}}}\right),\end{array}\end{array}$$

where T_dnf=max(0, 1.962(T_f−0.773)), i.e. a linear multiple of Eq. (1) that is 0 at T<20 ^∘C and unity at T=30 ^∘C. The temperature, iron, and light limitation terms are based on Pelagic Interactions Scheme for Carbon and Ecosystem Studies (PISCES) (Aumont et al,, 2015); the N-inhibition term is from CMOC (Zahariev et al., 2008) (CMOC implicitly combines nitrate and ammonium into a single inorganic N pool).

Denitrification is parameterized as a fraction of total remineralization that increases as a linear function of oxygen concentration for concentrations less than a threshold concentration O_mxd:

$$\begin{array}{}\text{(22)}& {\mathrm{N}}_{\mathrm{frxn}}=\mathrm{1}-{\displaystyle \frac{\min ({\mathrm{O}}_{\mathrm{2}},{\mathrm{O}}_{\mathrm{mxd}})}{{\mathrm{O}}_{\mathrm{mxd}}}}.\end{array}$$

Remineralization is then divided among oxygen (1−N_frxn), nitrate (0.875N_frxn), and ammonium (0.125N_frxn) assuming an average anammox contribution of 25 % (Babbin et al., 2014). We use this average ratio of anammox to classical denitrification to partition fixed N losses between NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$; the DIC sink and organic matter source associated with anammox are small and are neglected here.

2.4 Calcification, calcite dissolution, and alkalinity

In CanOE, calcification is represented by a prognostic detrital calcite pool with its own sinking rate (distinct from that of organic detritus), and calcite burial or dissolution in the sediments depends on the saturation state (100 % burial when Ω_C≥1, 100 % dissolution when Ω_C<1). Calcification is represented by a detrital calcium carbonate (CaCO₃) state variable, but no explicit calcifier groups. Detrital CaCO₃ sinks in the same fashion as detrital particulate organic carbon (POC), with a sinking rate independent of those for large and small organic detritus. Calcite production is represented as a fixed fraction of detritus production from small phytoplankton and small zooplankton mortality:

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{dCa}}{\mathrm{d}t}}& =m_{\mathrm{1}}\left({\mathrm{C}}_{\mathrm{ps}}+Z_s\right){\mathrm{P}}_{\mathrm{Ca}}+m_{\mathrm{2}}\left({\mathrm{C}}_{\mathrm{ps}}^{\mathrm{2}}+Z_{\mathrm{S}}^{\mathrm{2}}\right){\mathrm{P}}_{\mathrm{Ca}}\\ & -k_{\mathrm{Ca}}\mathrm{Ca}-w_{\mathrm{Ca}}{\displaystyle \frac{\mathrm{dCa}}{\mathrm{d}z}}.\end{array}\end{array}$$

Calcite dissolution occurs throughout the water column as a first-order process (i.e. no dependence on temperature or saturation state). Approximately 80 % of calcite produced is exported from the euphotic zone. In CanESM5-CanOE, burial in the sediments is represented as a simple “on/off” switch dependent on the calcite saturation state (zero when Ω_C<1 and 1 when Ω_C≥1). In CanESM5, calcification is parameterized by a temperature-dependent “rain ratio” (Zahariev et al., 2008), and 100 % burial of calcite that reaches the seafloor is assumed. Calcite burial in both models is balanced by an equivalent source of DIC and alkalinity at the ocean surface (in the same vertical column) as a crude parameterization of fluvial sources.

For each mole of calcite production, two moles of alkalinity equivalent are lost from the dissolved phase; the reverse occurs during calcite dissolution. There are additional sources and sinks for alkalinity associated with phytoplankton nutrient (NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$) uptake, organic matter remineralization, nitrification, denitrification, and dinitrogen fixation (Wolf-Gladrow et al., 2007; see Table S2 in the Supplement). The anammox reaction does not in itself contribute to alkalinity (Jetten at al., 2001), but there is a sink associated with ammonium oxidation to nitrite (the model does not distinguish between nitrite and nitrate). Autotrophic production of organic matter by anammox bacteria is a net source of alkalinity (Strous et al., 1998), but this source is extremely small (∼ 0.03 mol molN⁻¹) and is neglected here. Globally, the sources and sinks of alkalinity from the N cycle offset each other such that there is no net gain or loss as long as the global fixed N pool is conserved (see Sect. 2.5 below). If dinitrogen fixation and denitrification are allowed to vary freely, there will generally be a net gain or loss of fixed N and therefore of alkalinity.

2.5 External nutrient sources and sinks

External sources and sinks consist of river inputs, aeolian deposition, biological N₂ fixation, denitrification, mobilization of Fe from reducing sediments, loss of Fe to scavenging, and burial of calcium carbonate in the sediments. There is no burial of organic matter; organic matter reaching the seafloor is instantaneously remineralized. Aeolian deposition of Fe is calculated from a climatology of mineral dust deposition generated from offline (atmosphere-only) simulations with CanAM4 (von Salzen et al., 2013), with an Fe mass fraction of 5 % and a fractional solubility of 1.4 % in the surface layer. Subsurface dissolution is parameterized based on PISCESv2 (Aumont et al., 2015); the total dissolution is 6.35 %, with 22 % of soluble Fe input into the first vertical layer (see the Supplement). Iron from reducing sediments is also based on PISCES, with a constant areal flux of 1000 nmol m⁻² d⁻¹ in the first model level, declining exponentially with increasing seafloor depth (i.e. assuming that shelf sediments are the strongest source and the sediments become progressively more oxygenated with increasing seafloor depth) with an e-folding length scale of about 600 m. Scavenging of dissolved iron is first order with a high rate (2.5 d⁻¹) for concentrations in excess of 0.6 nM (Johnson et al., 1997). For concentrations below this threshold, the rate is much lower (0.001 d⁻¹) and is weighted by the concentration of organic detritus (Christian et al., 2002b), i.e.

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\mathrm{dFe}}{\mathrm{d}t}}=-{\mathrm{FeS}}_{\mathrm{Fe}\mathrm{1}}\min \left\{\left(D_{\mathrm{S}}+D_{\mathrm{L}}\right){\mathrm{P}}_{\mathrm{Fe}},\mathrm{1}\right\},\end{array}$$

where Fe is the dissolved iron concentration, D_S and D_L are the small and large detritus concentrations, S_Fe1 is the first-order scavenging rate in surface waters with abundant particulates, and P_Fe is an empirical parameter to determine the dependence on particle concentration (Table 1). The basis for this parameterization is that the rate of scavenging must depend not only on the concentration of iron but on the concentration of particles available for it to precipitate onto and assumes that detrital POC is strongly positively correlated with total particulate matter. Scavenging is treated as irreversible; i.e. scavenged Fe is not tracked and does not re-enter the dissolved phase.

N₂ fixation and denitrification vary independently in CanOE, so the global total N pool can change. Conservation is imposed by adjusting the global total N pool according to the difference between the gain from N₂ fixation and the loss to denitrification. A slight adjustment is applied to the nitrate concentration at every grid point, while preserving the overall spatial structure of the nitrate field. Adjustments are multiplicative rather than additive to avoid producing negative concentrations. This adjustment does not maintain (to machine precision) a constant global N inventory but is intended to minimize long-term drift, keeping it much smaller than the free surface error (see below). This adjustment is applied every 10 d and has a magnitude of approximately $\mathrm{7}\times {\mathrm{10}}^{-\mathrm{8}}$ of the total N.

When the total fixed N adjustment is applied, one mole of alkalinity is added (removed) per mole of N removed (added), to account for the alkalinity sources associated with N₂ fixation (creation of new NH${}_{\mathrm{4}}^{+}$) and denitrification (removal of NO${}_{\mathrm{3}}^{-}$) (Wolf-Gladrow et al., 2007; see Table S2 in the Supplement). As there is a 2 mol molN⁻¹ sink associated with nitrification, this formulation is globally conservative. As noted above, in CanOE CaCO₃ can dissolve or be buried in the sediments depending on the calcite saturation state. DIC and alkalinity lost to burial are reintroduced at the ocean surface, at the same grid point as burial occurs, providing a crude parameterization of river inputs so that global conservation is maintained (fresh water runoff contains no DIC or alkalinity). However, the OPA free surface formulation is inherently imperfect with regard to tracer conservation. Drift in total ocean alkalinity and nitrogen over time is on the order of 0.01 % and 0.03 % per thousand years, respectively.

2.6 Ancillary data

For first-order model validation, we have relied largely on global gridded data products rather than individual profile data. Global gridded data from World Ocean Atlas 2018 (WOA2018) (Locarnini et al., 2018; Zweng et al., 2018; Garcia et al., 2018a, 2018b) were used for temperature, salinity, and oxygen and nitrate concentration. DIC and alkalinity were taken from the GLODAPv2.2016b gridded data product (Key et al., 2015; Lauvset et al., 2016). Offline carbon chemistry calculations were done following the Best Practices Guide (Dickson et al., 2007) and the OMIP-BGC protocols (Orr et al., 2017), and they are identical to those used in the models except that constant reference concentrations were used for phosphate (1 µM) and silicate (10 µM).

There is no global gridded data product for Fe, but we have made use of the GEOTRACES Intermediate Data Product 2017 (Schlitzer et al., 2018), and the data compilations from Monterey Bay Aquarium Research Institute (MBARI) (Johnson et al., 1997; 2003) and North Pacific Marine Science Organization (PICES) Working Group 22 (Takeda et al., 2013). The latter two are concentrated in the Pacific, while GEOTRACES is more global. The combined data sets provide more than 10000 bottle samples from more than 1000 different locations (Fig. S10a in the Supplement) (excluding some surface transect data that involve frequent sampling of closely spaced locations along the ship track). More details about model comparison to these data compilations and the list of original references are given in the Supplement.

Satellite ocean colour estimates of surface chlorophyll were taken from the combined SeaWiFS/MODIS climatology described by Tesdal et al. (2016). Climatological satellite POC was downloaded from the NASA Ocean Color website and is based on the algorithm of Stramski et al. (2008) using MODIS Aqua data. This climatology differs slightly from the chlorophyll one in terms of years included and sensors utilized, but as only climatological concentrations are considered and each climatology covers ∼ 15 years, these differences will have negligible effect on the results presented. Satellite chlorophyll concentrations greater than 1 mg m⁻³ were excluded as these are mostly associated with coastal regions not resolved by coarse-resolution global ocean models.

Figure 2 Global distribution of oxygen (O₂) concentration in mmol m⁻³ at 400, 900, and 1400 m (rows) for CanESM5-CanOE, CanESM5, the mean for other (non-CanESM) CMIP6 models, and World Ocean Atlas 2018 (WOA2018) observations (columns). Numbers on the lower left are the mean model bias. Differences from the observation-based fields are shown in Fig. S2 in the Supplement.

Figure 3 Latitude–depth distribution (surface to 1750 m) of zonal mean oxygen concentration (O₂), oxygen concentration at saturation (O₂(sat)), and apparent oxygen utilization (AOU) in mmol m⁻³ for CanESM5-CanOE, CanESM5, the mean for other CMIP6 models, and observations (WOA2018). Note different colour scales for different rows. Numbers on the lower left are the mean model bias. Differences from the observation-based fields are shown in Fig. S2 in the Supplement.

CMIP6 model data were regridded by distance-weighted averaging using the Climate Data Operators (https://code.mpimet.mpg.de/projects/cdo/, last access: 19 May 2022) to a common grid (2×2^∘, 33 levels) to facilitate ensemble averaging. The vertical levels used are those used in GLODAP and in earlier (through 2009) versions of the World Ocean Atlas (e.g. Locarnini et al., 2010). For large-scale tracer distributions, using a 1 or 2^∘ grid makes little difference (for example, the spatial pattern correlation between CanESM5 and observed oxygen concentration at specific depths on a 1 or 2^∘ grid differs by an average of 0.0011). The years 1986–2005 of the historical experiment were averaged into climatologies or annual means, for meaningful comparison with observation-based data products. The CMIP6 historical experiment runs from 1850–2014 with atmospheric CO₂ concentration (and other atmospheric forcings) based on historical observed values. A single realization was used in each case (see Table S3); 20-year averages are used to minimize the effect of internal variability (e.g. Arguez and Vose, 2011; see Table S4). Where time series are shown, 5-year means are used.

Sampling among CMIP6 models was somewhat opportunistic, and the exact suite of models varies among the analyses presented. When we conducted a search for a particular data field, we included in the search parameters all models that published that field and repeated the search at least once for models that were unavailable the first time the search was executed. In some cases, model ensemble means excluded all but one model from a particular “family” (e.g. there are three different MPI-ESM models for which ocean biogeochemistry fields were published), as the solutions were found to be similar and would bias the ensemble mean towards their particular climate. The models used are ACCESS-ESM1-5, CESM2, CESM2-WACCM, CNRM-ESM2-1, GFDL-CM4, GFDL-ESM4, IPSL-CM6A-LR, MIROC-ES2L, MPI-ESM-1-2-HAM, MPI-ESM1-2-LR, MPI-ESM1-2-HR, MRI-ESM2-0, NorESM2-LM, NorESM2-MM, and UKESM1-0-LL. Details of which variables and realizations are used for which models are given in Table S3 in the Supplement.

3 Results

We first describe the large-scale distribution of oxygen, DIC, alkalinity, and the saturation state with respect to CaCO₃ that derives from these large-scale tracer distributions. Tracer distributions result partly from ocean circulation and partly from biogeochemical processes. An overall evaluation of the ocean circulation model is given in Swart et al. (2019a). Analyzing CanESM5 and CanESM5-CanOE (with identical circulation) as well as CanESM2 where possible (same biogeochemistry as CanESM5 but different circulation) allows us to separate the effects of physical circulation and biogeochemistry on evolving model skill with respect to large-scale tracer distributions. In subsequent sections, we address the main areas where CanESM5 and CanESM5-CanOE differ, such as the interaction of the iron and nitrogen cycles and plankton community structure. Finally, we present some temporal trends over the course of the historical experiment (1850–2014).

Figure 4 Taylor diagrams (Taylor, 2001) comparing modelled and observed distributions of oxygen at specific depths from 100 to 3500 m. The angle from the vertical indicates spatial pattern correlation. Distance from the origin indicates ratio of standard deviation in modelled vs. observed (WOA2018) fields. Red dots represent CanESM5-CanOE, blue dots CanESM5, small grey dots other CMIP6 models, and large grey dots the model ensemble mean for all CMIP6 models except CanESM5 and CanESM5-CanOE.

3.1 Distribution of oxygen

The spatial distribution of oxygen concentration ([O₂]) at selected intermediate depths (400, 900, and 1400 m) is shown in Fig. 2 for gridded data from WOA2018 and differences from that observational data product for CanESM5, CanESM5-CanOE, a model ensemble mean (MEM) of CMIP6 models (excluding CanESM5 and CanESM5-CanOE). The depths were chosen to span the depth range where low oxygen concentrations exist; these low-oxygen environments are of substantial scientific and societal interest and are sensitive to model formulation. The major features are consistent across the models. Both CanESM models as well as the MEM show elevated oxygen concentrations relative to observations, particularly in the North Pacific, the North Atlantic, and the Southern Ocean. In the Indian Ocean, both CanESM models show high oxygen concentrations in the Arabian Sea and deeper layers of the Bay of Bengal relative to observations and the MEM; these biases are somewhat smaller in CanESM5-CanOE than in CanESM5 (Fig. 2).

The ocean's oxygen minimum zones (OMZs) are mostly located in the eastern Pacific Ocean, the northern North Pacific, and the northern Indian Ocean; the spatial pattern changes with increasing depth (Fig. 2), but the OMZs are mostly located between 200 and 2000 m depth. Biases in the EBC regions are depth and model specific. CanESM5 shows particularly strong oxygen depletion at 1400 m in the eastern tropical Pacific. In the southeastern Atlantic, models tend to be biased low at the shallower depths and show somewhat more variation at greater depths (Fig. 2). Overall, [O₂] biases tend to be positive over large areas of ocean with the exception of some EBC regions, implying that models exaggerate the extent to which remineralization is concentrated in these regions. An alternate version of Fig. 2 that shows the modelled concentrations is given in Fig. S2 in the Supplement.

Figure 5 Total volume of ocean with oxygen (O₂) concentration less than (a) 6 mmol m⁻³ (mean for last 30 years of the historical experiment) and (b) 60 mmol m⁻³. Observation are from WOA2018.

Figure 6 Global distribution of aragonite saturation (Ω_A) at 400, 900, and 3500 m for CanESM5-CanOE, CanESM5, the mean for other CMIP6 models, and observations (GLODAPv2 and WOA2018). Note different colour scales for different depths. Numbers on the lower left are the mean model bias. Differences from the observation-based fields are shown in Fig. S2 in the Supplement.

The zonal mean oxygen concentration, saturation concentration, and apparent oxygen utilization (AOU) are shown in Fig. 3 for the same four cases. Again, the models generally show a positive bias in [O₂], particularly in high-latitude deep waters. The major ocean circulation features are reproduced fairly well in all cases (e.g. weaker ventilation of low-latitude subsurface waters, greater vertical extent of well-ventilated surface waters in the subtropics). The saturation concentration (a function of temperature and salinity) generally shows relatively little bias, implying that the bias in [O₂] arises mainly from remineralization and/or ventilation. AOU is lower than observed over much of the subsurface ocean. CanESM5 and CanESM5-CanOE show a high bias over much of the Northern Hemisphere that reflects the high concentrations in the North Pacific and North Atlantic (Fig. 2). The overall trend of bias with latitude in CanESM5 and CanESM5-CanOE is generally similar to the MEM, but the biases are larger. The bias in CanESM5 is generally slightly larger than in CanESM5-CanOE, except in the Arctic Ocean. Again, Fig. S2 in the Supplement includes a version of this plot that shows the modelled concentration fields.

The skill of each model with respect to the distribution of O₂ at different depths is represented by Taylor diagrams (Taylor, 2001) in Fig. 4. These diagrams allow us to assess how well the model reproduces the spatial distribution at a range of depths, because different physical and biogeochemical processes determine the distribution in different depth ranges. All of the CMIP6 models that were shown as an ensemble mean in Figs. 2 and 3 are shown individually. The large blue dots represent CanESM5, red CanESM5-CanOE, and grey the MEM; the smaller grey dots represent the individual models. CanESM5-CanOE shows slightly higher pattern correlation than CanESM5 at all depths. Both models compare favourably with the full suite of CMIP6 models, with r>0.85 for CanESM5 and r>0.9 for CanESM5-CanOE at all depths examined, and a normalized standard deviation within ±25 % of unity.

The total volume of ocean with [O₂] less than 6 mmol m⁻³ (the threshold for denitrification; Devol, 2008) and 60 mmol m⁻³ (a commonly used index of hypoxia) is shown in Fig. 5. The total volume is highly variable among models (note, however, that there are several clusters of related models with quite similar totals). CanESM5 and CanESM5-CanOE have among the lowest total volumes (i.e. the interior ocean is relatively well ventilated) and are among the nearest to the observed total. For [O₂]<60 mmol m⁻³ the bias is, nonetheless, quite large (i.e. the observed volume is underestimated by almost 50 % in both models). The volume of water with [O₂] below the denitrification threshold is overestimated in both CanESM5 and CanESM5-CanOE; CanESM5-CanOE has a much smaller total that is closer to the observed value. The bias in the spatial pattern of hypoxia (not shown) is generally similar to the bias in dissolved oxygen distribution (Fig. 2). The low-oxygen regions are generally more concentrated in the eastern tropical Pacific in the models than in observations, and the low-oxygen region in the northwest Pacific is not well reproduced in CanESM models.

3.2 Distribution of DIC, alkalinity, and CaCO₃ saturation

The spatial distribution of aragonite saturation state (Ω_A) at selected depths is shown in Fig. 6. The first two depths are the same as in Fig. 2, but a much greater depth is also included, as the length scale for CaCO₃ dissolution is greater than for organic matter remineralization. In this case, the observations are a combination of GLODAPv2 (Key et al., 2015; Lauvset et al., 2016) for DIC and alkalinity, and WOA2018 for temperature and salinity. CanESM5 and CanESM5-CanOE show an overall high saturation bias at the shallower depths, particularly in the North Atlantic, with a low bias found mainly in the eastern Pacific. The low saturation bias in the eastern tropical Pacific is substantially reduced in CanESM5-CanOE compared to CanESM5. On the other hand, CanESM5 generally does better than CanESM5-CanOE, or the MEM, at reproducing the low saturation states in the northwestern Pacific and the Bering Sea. Both CanESM models show a high saturation state bias in the North Atlantic and the well-ventilated regions of the north Pacific subtropical gyre; these biases are slightly smaller in CanESM5-CanOE. Maps of the calcite and aragonite saturation horizon (Ω=1) depth are shown in Fig. S3 in the Supplement; these generally confirm the same biases noted in Fig. 6.

Zonal mean distributions of aragonite saturation state (Ω_A), calcite saturation state (Ω_C), and carbonate ion concentration ([CO${}_{\mathrm{3}}^{--}$]) and the differences of the models from the observations are shown in Fig. 7 (Fig. S2 in the Supplement includes versions of Figs. 6 and 7 that show the modelled fields). The models generally compare well with the observations in the representation of the latitude–depth distribution of high- and low-saturation waters. CanESM5 has a high saturation bias in low-latitude surface waters that is somewhat reduced in CanESM5-CanOE. Both CanESM5 models show a high saturation bias in Northern Hemisphere intermediate (e.g. 200–1000 m) depth waters that is larger than in the MEM. This is primarily a result of low Ω in the North Atlantic Ocean (Fig. 6).

Figure 7 Latitude–depth distribution of zonal mean (surface to 1150 m) aragonite saturation state (Ω_A), calcite saturation state (Ω_C), and carbonate ion concentration ([CO${}_{\mathrm{3}}^{--}$]) in mmol m⁻³ for CanESM5-CanOE, CanESM5, the mean for other CMIP6 models, and observations (GLODAPv2 and WOA2018). Numbers on the lower left are the mean model bias. Differences from the observation-based fields are shown in Fig. S2 in the Supplement.

Taylor diagrams for a range of depths are shown for DIC in Fig. 8 and for Ω_A in Fig. 9 (for alkalinity; see Fig. S4 in the Supplement). As expected, the MEM generally compares favourably with the individual models (e.g. Lambert and Boer, 2001). CanESM5 and CanESM5-CanOE compare favourably with the full suite of CMIP6 models. CanESM5-CanOE shows a gain in skill relative to CanESM5, and both show improvement relative to CanESM2. At 400 m, CanESM2 stands out as having extremely high variance, which is mostly due to extremely high DIC concentrations occurring over a limited area in the eastern equatorial Pacific (not shown). This bias is present in CanESM5 and in CMIP6 models generally (Fig. 6) but involves much lower concentrations spread over a larger area.

Figure 8 Taylor diagrams comparing modelled and observed distributions of DIC at specific depths from 100 to 3500 m. Observations are from GLODAPv2 (Lauvset et al., 2016). Red dots represent CanESM5-CanOE, blue dots CanESM5, magenta dots CanESM2, small grey dots other CMIP6 models, and large grey dots the model ensemble mean for all CMIP6 models except CanESM5 and CanESM5-CanOE.

Figure 9 Taylor diagrams comparing modelled and observed (GLODAPv2 and WOA2018) distributions of Ω_A at specific depths from 100 to 3500 m. Symbol colours are as in Fig. 8.

3.3 N and Fe cycles

An important difference between CanESM5 and CanESM5-CanOE is the inclusion of a prognostic Fe cycle. The CMOC iron mask (Zahariev et al., 2008) was a pragmatic solution in the face of resource limitations but is inherently compromised as it can not evolve with a changing climate. The first-order test of a model with prognostic, interacting Fe and N cycles is whether it can reproduce the distribution of HNLC regions and the approximate surface macronutrient concentrations within these. CanESM5-CanOE succeeded by this standard, although the surface nitrate concentrations are biased low in the subarctic Pacific and equatorial Pacific and high in the Southern Ocean and in the global mean (Fig. 10).

Figure 10 Climatological seasonal cycle of surface nitrate concentration averaged for selected ocean regions. The thick red line represents CanESM5-CanOE, thick blue line CanESM5, thick black line observations (WOA2018), thin grey lines individual CMIP6 models, and thick grey line the model ensemble mean (excluding CanESM5 and CanESM5-CanOE). Regional boundaries are given in Table S5 and Fig. S5.

The seasonal cycle of the zonal mean surface nitrate concentration for a selection of CMIP6 models is shown in Fig. 11. CanESM5, CanESM5-CanOE, and CNRM-ESM2-1 reproduce the equatorial enrichment and the low concentrations in the tropical–subtropical latitudes fairly well. Some models either have very weak equatorial enrichment (MPI-ESM1-2-LR) or too high a concentration in the off-equatorial regions (UKESM1-0-LL, NorESM2-LM). UKESM1-0-LL has very high concentrations throughout the low-latitude Pacific, which biases the ensemble mean (Fig. 11). Figure S6 in the Supplement shows the same data as Fig. 11 but for a more limited latitude range to better illustrate model behaviour in the tropics. CanESM5, CanESM5-CanOE, and CNRM-ESM2-1 reproduce the seasonal cycle of tropical upwelling (e.g. Philander and Chao, 1991), with highest concentrations in summer.

Figure 11 Climatological seasonal cycle of zonal mean surface nitrate concentration for a selection of CMIP6 models, a model ensemble mean (MEM) excluding CanESM5 and CanESM5-CanOE, and an observation-based data product (WOA2018). An alternate version showing only latitudes <20^∘ is given in Fig. S6.

The surface distribution of dissolved iron (dFe) in various CMIP6 models is shown in Fig. 12. For Fe there is no observation-based global climatology with which to compare the model solutions (some comparisons to available profile data are shown in Fig. S10b–h in the Supplement). CanESM5-CanOE shows a similar overall spatial pattern to other models, and generally falls in the middle of the spread, particularly regarding concentrations in the Southern Ocean. Several models show extremely high concentrations in the tropical–subtropical North Atlantic (Sahara outflow region). CanESM5-CanOE, along with CNRM-ESM2-1 and CESM2, has much less elevated concentrations in this region, due to lower deposition or greater scavenging or both. CanESM5-CanOE has its lowest concentration in the eastern subtropical South Pacific, which is common to many models (Fig. 12). The area of strong surface depletion is generally more spatially restricted in CanESM5-CanOE than in other models, and surface dFe concentrations are greater over large areas of the Pacific. Both the north–south and east–west asymmetry of distribution in the Pacific is greater in CanESM5-CanOE than in most other models, some of which show the South Pacific minimum extending westward across the entire basin, and others into the Northern Hemisphere. Only in CESM2 is this minimum similarly limited to the southeast Pacific.

Figure 12 Global distribution of dissolved iron (dFe) concentration (nmol m⁻³) at the ocean surface for CanESM5-CanOE and other CMIP6 models that published this field. Concentrations exceeding 1000 nmol m⁻³ are masked in white. CanESM5 is not included because it does not have prognostic iron.

The mean depth profiles of dFe are shown in Fig. 13. Some models show more of a “nutrient-type” (increasing with depth due to strong near-surface biological uptake and subsequent remineralization) profile, some a more “scavenged-type” (maximal at the surface, declining with depth) profile (see Li, 1991; Nozaki, 2001), and others a hybrid profile (increasing downward but with a surface enrichment). CanESM5-CanOE is at the “nutrient-type” end of spectrum with a generally monotonic increase with depth to a near-constant deep-water concentration of 0.6 nM and a very slight near-surface enrichment (see also Fig. S10b, c in the Supplement).

Figure 13 Global mean depth profiles of dissolved iron concentration for CanESM5-CanOE and other CMIP6 models that published this field. GFDL-CM4 is excluded because it has very high concentrations (>2000 nmol m⁻³) near the surface. The thick red line represents CanESM5-CanOE, the thin grey lines individual CMIP6 models, and the thick grey line the model ensemble mean (excluding CanESM5-CanOE and GFDL-CM4).

Mean surface nitrate and dFe concentrations for selected ocean regions are shown in Fig. 14. CanESM5-CanOE shows concentrations that are within the range of CMIP6 models, although in some cases at the higher or lower end. Surface nitrate concentrations generally compare favourably with the observation-based climatology, but are biased low in HNLC regions other than the Southern Ocean. These biases are not necessarily a consequence of having too much or too little iron. For example, in the Southern Ocean CanESM5-CanOE has among the highest surface nitrate concentrations, but it also has some of the highest dFe concentrations, and the high nitrate bias is present in CanESM5 as well. Comparisons with the limited GEOTRACES data available suggest that near-surface dFe concentrations in the Southern Ocean are biased high rather than low in CanESM5-CanOE (not shown). One region where there does seem to be a strong correlation between surface nitrate and dFe concentrations is the western subarctic Pacific. All but two models (CNRM-ESM2-1, NorESM2-LM) fall along a spectrum from high Fe/low nitrate to low Fe/high nitrate. CanESM5-CanOE falls near the high Fe/low nitrate end of the range.

Figure 14 Mean surface nitrate (NO₃) vs. dissolved iron (dFe) concentrations in different oceans, including the major high-nutrient/low-chlorophyll (HNLC) regions. CanESM5-CanOE is shown as a red dot and other CMIP6 models as grey dots (CanESM5 is not included because it does not have iron). Observed NO₃ is shown as a vertical black line as there are no observational estimates of dFe concentration. For GFDL-CM4, nitrate is estimated as phosphate ×16. Region definitions are given in Table S5 and Fig. S5.

Surface nitrate concentrations along the Pacific Equator during the upwelling season (June–October) for CanESM5 and CanESM5-CanOE are shown in Fig. 15. The range of other CMIP6 models is not shown here because it is large and therefore adds little information (see Figs. 11 and S6). CanESM5-CanOE better represents the east-west gradient, while CanESM5 has slightly higher concentrations in the core upwelling region. Both models underestimate the highest concentrations around 100^∘ W. Although some localized maxima in this data product are due to undersampling, equatorial upwelling is strong at this location (e.g. Lukas, 2001) and the spatial coherence of the data strongly suggests that this maximum accurately reflects reality. It should be noted that CanESM5 iron limitation is calculated from a version of the same data product; however, the Fe mask is based on the minimum nitrate concentration over the annual cycle, whereas the data shown here are for the upwelling season.

Figure 15 Surface nitrate (NO₃) concentrations along the Pacific Equator (mean from 2^∘ S–2^∘ N) during the upwelling season (June–October) for CanESM5-CanOE (red), CanESM5 (blue), and WOA2018 observations (black).

3.4 Plankton biomass, detritus, and particle flux

The relative abundance of the four plankton groups is shown in Fig. 16 for a range of ocean regions. Both CanESM models mostly compare favourably with observation-based estimates of phytoplankton biomass, except in the tropics where CanESM5-CanOE has very high biomass. Both CanESM models have low phytoplankton biomass in the North Atlantic. In the North Pacific and the Southern Ocean, CanESM5-CanOE reproduces the observation-based estimates well and CanESM5 slightly less well. The general pattern is that large and small phytoplankton have similar abundance and are substantially more abundant than zooplankton.

Figure 16 Annual mean surface ocean concentration of large and small phytoplankton and zooplankton in CanESM5-CanOE (red) and of phytoplankton and zooplankton in CanESM5 (blue) for the representative ocean regions shown in Fig. 14. Observational estimates (black) are for phytoplankton biomass calculated from satellite ocean colour estimates of surface chlorophyll (SeaWiFS/MODIS; Tesdal et al. 2016), assuming a carbon-to-chlorophyll ratio of 50 g g⁻¹. Region definitions are given in Table S5 and Fig. S5.

Figure 17 Mean annual cycle of surface chlorophyll for the representative ocean regions shown in Figs. 14 and 16. CanESM5-CanOE large and small phytoplankton concentrations are shown separately and combined (red) along with CanESM5 (blue) and observational estimates (black). Region definitions are shown in Table S5 and Fig. S5.

Figure 18 Climatological surface particulate organic carbon (POC) vs. chlorophyll for CanESM5-CanOE (red) and observations (black). Data are for all ocean grid points (2×2^∘ uniform global grid) for all months of the year where observational data are available. Model POC is offset 17 mg m⁻³ for illustrative purposes. Chlorophyll concentrations >1 mg m⁻³ are excluded as they largely represent coastal areas poorly resolved by coarse-resolution global ocean models.

Part of the rationale for multiple food chains is that they better represent the way that actual plankton communities adapt to different physical ocean regimes and therefore are better able to simulate distinct ocean regions with a single parameter set (e.g. Chisholm, 1992; Armstrong, 1994; Landry et al., 1997; Friedrichs et al., 2007). The expectation is that small phytoplankton will be more temporally stable and large phytoplankton will fluctuate more strongly between high and low abundances. The mean annual cycles of surface chlorophyll largely conform to this pattern; e.g. in the North Atlantic and the western subarctic Pacific, large phytoplankton are dominant in summer and much more variable over the seasons (Fig. 17). Compared to observations, CanESM5 models underestimate the amplitude of the seasonal cycle in the North Atlantic and overestimate it in the North Pacific. CanESM5 shows a stronger and earlier North Atlantic spring bloom compared to CanESM5-CanOE; the observations are in between the two in terms of timing, and both models underestimate the amplitude (Fig. 17). In the tropics, the seasonal cycle is weak. CanESM5-CanOE in the tropical Atlantic shows the expected seasonal cycle but not the expected dominance of large phytoplankton in summer. CanESM5-CanOE generally overestimates the total near-surface chlorophyll in both the tropical Pacific and the tropical Atlantic.

Zooplankton biomass (especially microzooplankton) is also somewhat difficult to test against observations, but our model concentrations appear to be biased low. Stock et al. (2014) estimated depth-integrated biomass of phytoplankton, mesozooplankton, and microzooplankton for a range of oceanic locations in which intensive field campaigns have occurred (estimates of microzooplankton biomass are relatively sparse). They found that in most locations phytoplankton and (combined) zooplankton biomass are of comparable magnitude, whereas in CanESM5-CanOE zooplankton biomass is consistently lower (Fig. 16). The global integral biomass of mesozooplankton is about an order of magnitude less than the 0.19 PgC estimated by Moriarty and O'Brien (2013). The CanESM5 total of 0.14 Pg is relatively close to the Moriarty estimate but implicitly includes microzooplankton.

Surface chlorophyll and POC for CanESM5-CanOE and for ocean colour observational data are shown in Fig. 18 (POC in the model is the sum of phytoplankton, microzooplankton, and detrital carbon). The observations have a lower limit for POC that is not present in the model (∼ 17 mgC m⁻³), which is unsurprising given the processes neglected in the model; i.e. in regions of very low chlorophyll, there is still substantial dissolved organic carbon, bacteria that consume it, and microzooplankton that consume the bacteria and produce particulate detritus. The observational data show a fairly linear relationship at low concentrations, but with a curvature that implies a greater phytoplankton fraction in more eutrophic environments (see Chisholm, 1992). The model, by contrast, shows a fairly linear relationship over the whole range of concentrations. In other words, the phytoplankton share of POC is higher and more constant in the model than in the observations. The living biomass (phytoplankton and microzooplankton) fraction of total POC in CanOE is generally in excess of 50 % (not shown), which is implausible for a real-world oceanic microbial community (e.g. Christian and Karl, 1994) but consistent with the relatively low rates of export from the euphotic zone.

Export production for a range of CMIP6 models is shown in Fig. 19a. CanESM5-CanOE is at the low end of the range. Observations are not shown because the range of observational estimates covers the entire range of model estimates (e.g. Siegel et al., 2016). Note also that CanESM5 export is quite a bit lower than in CanESM2, which is relatively high for CMIP5 models (not shown). The difference between CanESM2 and CanESM5 is attributable primarily to different circulation, although the different initialization fields for nitrate might also play a small role. The lower rate in CanESM5-CanOE is consistent with the above results regarding plankton community structure (e.g. the concentration of detritus is generally low compared to living biomass), as well as the lower sinking rate for small detritus. The latitudinal distribution of export is shown in Fig. 19b. CanESM5 shows very high export in the midlatitudes of the Southern Ocean, similar to CanESM2 (not shown). Both CanESM5 and CanESM5-CanOE show latitudinal patterns consistent with the range of other CMIP6 models. CanESM5 has slightly greater export in the equatorial zone; in both CanESM5 and CanESM5-CanOE, the equatorial enrichment attenuates very rapidly with latitude and the rates are low in the subtropics.

Figure 19 (a) Global total export production (epc100) in PgC yr⁻¹ (b) and zonal mean export production in molC m⁻² yr⁻¹​​​​​​​ according to selected CMIP6 models (mean for 1985–2014 of historical experiment). The thick red line represents CanESM5-CanOE, thick blue line CanESM5, thin grey lines individual CMIP6 models, and thick grey line the model ensemble mean (excluding CanESM5 and CanESM5-CanOE).

3.5 Historical trends

Cumulative ocean uptake of CO₂ is shown in Fig. 20 for the historical experiment (1850–2014). CanESM models are biased low relative to observation based estimates (∼ 145 PgC; see Friedlingstein et al., 2020) and the MEM (144 PgC, Fig. 20) but fall well within the spread of CMIP6 models. Some of the difference may be attributable to differences in the way cumulative uptake is calculated in models vs. observations (Bronselaer et al., 2017), although this should apply to other CMIP6 models as well. CanESM5-CanOE has lower cumulative uptake than CanESM5 by ∼ 10 PgC. As the models were not fully equilibrated when the historical run was launched, this difference does not necessarily arise from the biogeochemical model structure; part of the difference can be attributed to differences in the spinup (see Séférian et al., 2016). The drift in the piControl experiment over the 165 years from the branching off of the historical experiment is −10.0 PgC in CanESM5-CanOE and −5.1 PgC in CanESM5 (see Table S6 in the Supplement), so drift accounts for about half (48 %) of the difference in net ocean CO₂ uptake. The spatial distribution of anthropogenic DIC is very similar between CanESM5 and CanESM5-CanOE (Fig. S7 in the Supplement). CanESM5 and CanESM5-CanOE show a high bias in near-surface DIC relative to alkalinity (a measure of the ocean's capacity to absorb CO₂) in the midlatitudes of both hemispheres (Fig. S8 in the Supplement), which may in part explain the weak uptake of CO₂.

Figure 20 Cumulative ocean uptake of carbon dioxide (CO₂) as anthropogenic dissolved inorganic carbon (AnthDIC) in PgC over the course of the historical experiment (1850–2014). Data are shown as successive 5-year means. CMIP6 mean (thick grey line) indicates ensemble mean for CMIP6 models (thin grey lines) excluding CanESM5 (blue) and CanESM5-CanOE (red). An observation-based estimate of 145±20 PgC (Friedlingstein et al., 2020) is shown for a nominal year (2014) (black).

The long-term trend in global total export production is shown in Fig. 21. The model values must be normalized in order to compare trends, since the differences among means are large compared to the changes over the historical period (Fig. 19). Such trends are difficult or impossible to meaningfully constrain with observations, but the general expectation has been that export will decline somewhat due to increasing stratification (e.g. Steinacher et al., 2010). CanESM5 shows a greater decline than most other CMIP6 models, while CanESM5-CanOE is more similar to non-CanESM models. The change in CanESM5 is geographically widespread and not concentrated in a specific region or regions: export is maximal in the tropics and the northern and southern midlatitudes (Fig. 19b) and declines over the historical period in all of these regions (Fig. S9 in the Supplement). In CanESM5-CanOE, export declines in the same regions, but the magnitude of the change is smaller, and in the Southern Ocean increases and decreases in different latitude bands largely offset each other.

Figure 21 Change in export production (epc100) over the course of the historical experiment (1850–2014), normalized to the 1850–1900 mean. Data are shown as successive 5-year means. The thick red line represents CanESM5-CanOE, thick blue line CanESM5, thin grey lines other CMIP6 models, and thick grey line the ensemble mean of non-CanESM models.

The trend in the volume of ocean water with O₂ concentration less than 6 or 60 mmol m⁻³ is shown in Fig. 22. Again, the totals are normalized to a value close to the preindustrial, as the differences among models are large (Fig. 5). For the volume with <60 mmol m⁻³, CanESM models show relatively little change; in CanESM5, the volume actually declines slightly, while in CanESM5-CanOE it increases, but the total change is <1 % in each case. As with the baseline volumes, the range among models is large, with one model showing an increase approaching 10 % of the total volume estimated for WOA2018 (Figs. 5b and 22b). For the volume with <6 mmol m⁻³ (Fig. 22a), CanESM models are among the most stable over time. In CanESM5, the volume again declines, although this is within the range of internal variability. Again, some models show fairly large excursions, but in this case none show a strong secular trend over the last half century.

Figure 22 (a) Change in total ocean volume with oxygen (O₂) concentration less than (a) 6 mmol m⁻³ and (b) 60 mmol m⁻³ over the course of the historical experiment (1850–2014), normalized to the 1850–1870 mean. Data are shown as successive 5-year means. The thick red line represents CanESM5-CanOE, thick blue line CanESM5, and thin grey lines other CMIP6 models.

4 Discussion

CanESM5 and CanESM5-CanOE are new coupled ocean–atmosphere climate models with prognostic ocean biogeochemistry. The two have the same physical climate (in experiments with specified atmospheric CO₂) and differ only in their ocean biogeochemistry components. CanESM5-CanOE has a much more complex biogeochemistry model including a prognostic iron cycle. We have presented results that assess how these two models simulate the overall distribution of major tracers like DIC, alkalinity, nitrate, and oxygen, as well as analyses of the interaction of the iron and nitrogen cycles, plankton community structure, export of organic matter from the euphotic zone, and historical trends over 1850–2014.

The overall distribution of major tracers indicates that both models do a reasonable job of simulating both biogeochemical (e.g. export and remineralization of organic matter) and physical (e.g. deep and intermediate ocean ventilation) processes. The volume of ocean with oxygen concentration below 6 or 60 µM compares favourably with other CMIP6 models (Fig. 5) and is among the most stable over historical time (Fig. 22). CanESM5-CanOE has a substantially lower volume of water with [O₂]<6 µM than CanESM5 and much closer to observation-based estimates (Fig. 5). Both models are biased slightly low in terms of historical uptake of anthropogenic CO₂, which may indicate weak Southern Ocean upwelling or too-shallow remineralization of DIC or both (Fig. 20). The spatial distribution of anthropogenic DIC is very similar between the two models (Fig. S7 in the Supplement), which is expected as it is mainly a function of the physical ocean model circulation. However, CanESM5 has higher concentrations in the main areas of accumulation, particularly the North Atlantic and the Southern Ocean. This probably indicates more efficient removal and export of “natural” DIC by the plankton, particularly in the Southern Ocean upwelling zone (Fig. 19), and deeper average remineralization, with the caveat that the preindustrial control simulations had different degrees of equilibration when the historical experiment was launched (see Séférian et al., 2016, Table S6 in the Supplement).

Analyses of phytoplankton and zooplankton biomass concentrations show that CanESM5 and CanESM5-CanOE compare somewhat favourably with available observational data but do have distinct biases. In particular, both zooplankton biomass and detrital organic matter concentration tend to be very low in CanESM5-CanOE; the total biomass of the plankton community and the standing crop of particulate organic matter are dominated by phytoplankton (e.g. Fig. 17). Regional biases differ between the two models, with CanESM5-CanOE showing excessively large phytoplankton biomass in the tropics. We note, however, that the seasonal cycle of equatorial upwelling and the formation of the equatorial Pacific HNLC are reproduced rather well by our models (e.g. Figs. 11, 15, and S6), and that CanESM5-CanOE is the first CanESM model to have genuinely simulated this as an emergent property (see Sect. 3.3). In CanESM5-CanOE, decoupling of large and small phytoplankton populations associated with seasonal upwelling or convection (see below) is observed in some regions but not others.

Global export production is biased low, particularly in CanESM5-CanOE. This is due in part to the biogeochemical model and in part to ocean circulation. CanESM5 has the same ocean biology as CanESM2 but a different physical ocean model, and global ocean export production is substantially lower in CanESM5. It is lower still in CanESM5-CanOE (Fig. 19). We note that CanESM5 performs better than CanESM2 on most metrics of physical ocean model evaluation (Swart et al., 2019a) and shows a more realistic distribution of major tracers like DIC (Fig. 8). While the range of observation-based estimates of global ocean export production is large and encompasses the full range of CMIP5 and CMIP6 models, the change between CanESM2 and CanESM5 is large. Changes in the physical ocean are not entirely independent of the biogeochemistry model even when the latter is ostensibly identical. In CanESM2 and CanESM5, iron limitation is specified as a spatially static “mask” based on the observed distribution of surface nitrate, and it is possible that in these two models ocean upwelling occurs in different places relative to the specified boundary of the region of Southern Ocean iron imitation (Fig. 3 of Zahariev et al., 2008). It is also possible that the lower export production in CanESM5-CanOE is due to low iron supply to the surface waters of the Southern Ocean, but comparisons with available observations do not suggest that this is the case. Several biases are common to CanESM5 and CanESM5-CanOE that relate to Southern Ocean upwelling (high Southern Ocean surface nitrate concentration, low export production, weak anthropogenic CO₂ uptake) and so are probably more attributable to the physical ocean model than to the Fe submodel. The difference between CanESM2 and CanESM5 bears this out.

The development of CanOE was undertaken in response to some of the most severe limitations of CanESM2. Many of the additional features that CanOE introduces were already in the models published by other centres even in CMIP5. In addition to CMOC (Zahariev et al., 2008), previous models developed by members of our group include Denman and Peña (1999; 2002), Christian et al. (2002a, b), Christian (2005), and Denman et al. (2006). Christian et al. (2002a) had a prognostic Fe cycle and multiple phytoplankton and zooplankton species but had fixed elemental ratios. Christian (2005) incorporated a cellular-regulation model but only for a single species and without Fe limitation. Christian (2005) had prognostic chlorophyll, whereas Denman and Peña (1999; 2002) and Christian et al. (2002a) used an irradiance-dependent diagnostic formulation. Christian et al. (2002a) used multiplicative (Franks et al., 1986) grazing, which creates stability in predator–prey interactions but severely limits phytoplankton biomass accumulation under nutrient-replete conditions.

One of the most important lessons from Christian et al. (2002a, b) was that when a fixed $\mathrm{Fe}/\mathrm{N}$ ratio is employed, sensitivity to this parameter is extreme. Because Fe cell quotas are far more variable than N, P, or Si quotas, treating this parameter as constant results in the specified value influencing the overall solution far more than any other parameter. CanESM5-CanOE largely succeeded in creating a prognostic Fe–N limitation model that produces HNLC conditions in the expected regions (Figs. 10, 11, 14, 15, S6), although surface nitrate concentration is low relative to observation-based estimates in some cases. External Fe sources and scavenging parameterizations will be revisited and refined in future versions. In CanESM5-CanOE, the scavenging model is very simple, with distinct regimes for concentrations greater or less than 0.6 nM; scavenging rates are very high above this threshold which causes deep-water concentrations to converge on this value. The generally nutrient-like profile suggest that in CanOE the scavenging rate is quite low for concentrations below 0.6 nM (Fig. 13; see also Fig. S10h in the Supplement). We note that the aeolian mineral dust deposition field employed here is derived from the CanESM atmosphere model; these processes are not presently interactive but could be made so in the future.

A particular issue with CanESM2 was that extremely high concentrations of nitrate occurred under the EBC upwelling regions. This error resulted from spreading denitrification out over the ocean basin so that introduction of new fixed N from N₂ fixation would balance denitrification losses within each vertical column, whereas in the real world denitrification is highly localized in the low oxygen environments under the EBCs. CanESM2 did not include oxygen, but CanESM5 incorporates oxygen as a “downstream” tracer that does not feed back on other biogeochemical processes. The incorporation of a more process-based denitrification parameterization in CanESM5-CanOE is independent of the many other processes that are present in CanESM5-CanOE but not in CanESM5: a CMOC-like model with prognostic denitrification is clearly an option. We chose not to include explicit, oxygen-dependent denitrification in CanESM5 because we wanted to maintain a CMOC-based model as close to the CanESM2 version as possible, and because oxygen would not then be a downstream tracer that does not affect other processes.

Plankton community structure in CanESM5-CanOE is somewhat biased toward high concentrations of phytoplankton, low concentrations of zooplankton and detritus, and low export (Sect. 3.4). In the development phase, a fair number of experiments were conducted with various values of the grazing rates and detritus sinking speeds. A wide range of values of these parameters was tested, with no resulting improvement in the overall results. Possibly the detrital remineralization rates are too high, although primary production is also on the low end of the CMIP6 range (not shown) and would probably decline further if these rates were decreased. The model was designed around the Armstrong (1994) hypothesis of “supplementation” vs. “replacement”; i.e. small phytoplankton and their grazers do not become much more abundant in more nutrient-rich environments but rather stay at about the same level and are joined by larger species that are absent in more oligotrophic conditions (see also Chisholm, 1992; Landry et al., 1997; Friedrichs et al., 2007). The results presented here suggest that this was partially achieved, but further improvement is possible (Fig. 17).

As to whether the gains in skill with CanESM5-CanOE justify the extra computational cost, Taylor diagrams (Figs. 4, 8, 9, and Fig. S4 in the Supplement) show a modest but consistent gain in skill at simulating the major biogeochemical species (O₂, DIC, alkalinity) across variables and depths, especially for alkalinity at middle depths (Fig. S4 in the Supplement), for which CanESM5 displays the least skill relative to other fields or depths. Other processes that are highly parameterized in CanESM5, such as calcification and CaCO₃ dissolution, were not addressed in detail in this paper but are an important factor in determining the subsurface distribution of alkalinity. Again, we emphasize that we are simulating as an emergent property of a process-based model something that is parameterized in CanESM5 (as previously noted for surface nitrate concentration in HNLC regions), and doing at least as well in terms of model skill. As a general rule, the potential for improving skill and achieving better results in novel environments (e.g. topographically complex regional domains like the Arctic Ocean and the boreal marginal seas) is expected to be greater in less parameterized, more mechanistic models (e.g. Friedrichs et al., 2007; Tesdal et al., 2016). Inclusion of a prognostic iron cycle and C/N/Fe stoichiometry also open up additional applications and scientific investigations that are not possible with CMOC.

An updated version of CanESM5 with prognostic denitrification is clearly possible. However, for the reasons discussed above, a prognostic Fe cycle with a fixed phytoplankton $\mathrm{Fe}/\mathrm{N}$ remains problematic, and the model would still have a single detritus sinking speed and remineralization length scale. We are also developing CanOE for regional downscaling applications (Hayashida, 2018; Holdsworth et al., 2021). The regional domains have complex topography and prominent continental shelf and slope, and the single remineralization length scale in CMOC may not be well suited to such an environment. The number of tracers in CanOE is not particularly large compared with other CMIP6 models. We expect to further refine CanOE and its parameterizations, evaluate it against new and emerging ocean data sets (e.g. GEOTRACES, biogeochemical Argo), and incrementally improve CMOC (which we will maintain for a wide suite of physical climate experiments for which ocean biogeochemistry is not central to the purpose). For CMIP6, we chose to keep CMOC as close to the CanESM2 version as possible. This strategy allows us to quantify how much of the improvement in model skill is due to the physical circulation, as is illustrated by greater skill with respect to DIC (Fig. 8) and alkalinity (Fig. S4 in the Supplement), particularly at intermediate depths (400–900 m). The CanESM terrestrial carbon model is also undergoing important new developments (e.g. Asaadi and Arora, 2021) and we expect CanESM to continue to offer a credible contribution to global carbon cycle studies, as well as advance regional downscaling and impact science.