Optimal parameters for the ocean's nutrient, carbon, and oxygen cycles compensate for circulation biases but replumb the biological pump

Accurate predictive modeling of the ocean's global carbon and oxygen cycles is challenging because of uncertainties in both biogeochemistry and ocean circulation. Advances over the last decade have made parameter optimization feasible, allowing models to better match observed biogeochemical fields. However, does fitting a biogeochemical model to observed tracers using a circulation with known biases robustly capture the inner workings of the biological pump? Here we embed a mechanistic model of the ocean's coupled nutrient, carbon, and oxygen cycles into two circulations for the current climate. To assess the effects of biases, one circulation (ACCESS-M) is derived from a climate model and the other from data assimilation of observations (OCIM2). We find that parameter optimization compensates for circulation biases at the expense of altering how the biological pump operates. Tracer observations constrain pump strength and regenerated inventories for both circulations, but ACCESS-M export production optimizes to twice that of OCIM2 to compensate for ACCESS-M having lower sequestration efficiencies driven by less efficient particle transfer and shorter residence times. Idealized simulations forcing complete Southern Ocean nutrient utilization show that the response of the optimized system is sensitive to the embedding circulation. In ACCESS-M, Southern Ocean nutrient and dissolved inorganic carbon (DIC) trapping is partially short circuited by unrealistically deep mixed layers. For both circulations, intense Southern Ocean production deoxygenates Southern-Ocean-sourced deep waters, muting the imprint of circulation biases on oxygen. Our findings highlight that the biological pump's plumbing needs careful assessment to predict the biogeochemical response to ecological changes, even when optimally matching observations.

1 Introduction

The ocean's nutrient, carbon, and oxygen cycles are of central importance for the climate and the fertility of the ocean. The cycling rates and patterns are shaped by the subtle interplay between the ocean circulation and the generation, transport, and respiration of organic matter (the biological pump) as well as by air–sea gas exchange. Building robust predictive models of the ocean's biological pump poses a formidable challenge because of the myriad of biogeochemical processes that must be parameterized. Current prognostic earth-system models are computationally expensive, which prohibits systematic parameter space exploration. Relatedly, the long timescales of the global ocean circulation make it expensive to spin up earth-system models by brute-force time stepping, and differences among models have been shown to reflect differences in spin-up strategy (Séférian et al., 2016), compounding the parametric uncertainties. As a result of these computational challenges, simulations with earth-system models have made widely varying predictions of future ocean biogeochemistry (Bopp et al., 2013; Cocco et al., 2013; Henson et al., 2022).

To reduce the parametric uncertainties, biogeochemical parameters can be objectively determined by minimizing the quadratic mismatch between model-predicted and observed tracer distributions. Systematic parameter optimization is made possible by embedding the biogeochemical model in climatological steady flow and efficiently solving the generally nonlinear equations of the system directly for steady state using Newton-type implicit solvers (e.g., Kwon and Primeau, 2006). This approach exploits the matrix representation of the discretized advective–diffusive flux-divergence operator (the “transport matrix”; e.g., Khatiwala et al., 2005; Primeau, 2005; Chamberlain et al., 2019) and has been applied to a number of biogeochemical cycles embedded in data-assimilated ocean circulation models (e.g., Primeau et al., 2013; DeVries, 2014; Teng et al., 2014; Holzer et al., 2014; Pasquier and Holzer, 2017; DeVries and Weber, 2017; Wang et al., 2019, to cite a few).

When the circulation is data-assimilated to provide a realistic representation of the ocean's advective–eddy-diffusive transport, optimizing biogeochemical parameters is a natural strategy for obtaining robust representations of the ocean's biogeochemical cycles. However, Kriest and Oschlies (2015) demonstrated that optimal parameters for the Model of Oceanic Pelagic Stoichiometry (MOPS; Kriest and Oschlies, 2015) differ depending on the circulation model that is used. This raises the following questions: if the circulation is taken from a climate-model simulation with known biases in the ocean's physical state, do optimized biogeochemical parameters still provide a reliable estimate of the ocean's biogeochemistry, and are the simulated responses of the system to either biogeochemical or physical perturbations robust?

To answer these questions, we develop a relatively simple model (dubbed PCO2 here) of the coupled phosphorous, carbon, and oxygen cycles and contrast the properties of its biological pump and its response to perturbations depending on whether the model is optimized for a data-assimilated circulation or for a climate-model-derived circulation. The PCO2 model was constructed with a particular focus on capturing the coupling between oxygen and organic-particle respiration. We use a mechanistic formulation of nutrient uptake rather than an observation-based parameterization so that biological production can fully respond to the embedding circulation (and also to make PCO2 suitable for exploring climate-change scenarios in future studies). PCO2 improves on the well-established approach of the Ocean Carbon-Cycle Model Intercomparison Project (OCMIP; Najjar et al., 1992) by not constraining the functional form of the particle-flux-divergence profile to be either a power law or an exponential. Instead, we explicitly model sinking biogenic particles and let them react with the ambient oxygen in a simple temperature-dependent model of microbial respiration (DeVries and Weber, 2017; Laufkötter et al., 2017) to determine the flux divergence and respiration rates mechanistically.

We focus on a decadal-mean circulation derived from the ACCESS1.3 climate model (Chamberlain et al., 2019; Holzer et al., 2020) for the 1990s. To assess PCO2's biogeochemistry and biological pump when embedded and optimized in this climate-model circulation, we compare with PCO2 embedded and optimized in a data-assimilated ocean circulation (OCIM2; DeVries and Holzer, 2019). OCIM2 has been optimized so that its transport produces tracer fields that are as close as possible to observations. Thus, OCIM2 provides a realistic reference circulation with minimal biases. To explore the effect of optimized biogeochemical parameters on the response to biogeochemical perturbations, we consider idealized perturbations in which biological production in the Southern Ocean is intensified to cause a nearly complete nutrient drawdown. Similar perturbations have been used previously to quantify the Southern Ocean's key role in supplying the rest of the ocean with preformed nutrients (e.g., Sarmiento et al., 2004; Marinov et al., 2006; Holzer and Primeau, 2013) and to illustrate Southern Ocean nutrient trapping (Primeau et al., 2013). Here, we establish how the response to Southern Ocean nutrient drawdown differs depending on whether biogeochemistry is optimized for the OCIM2 or ACCESS circulations.

We find that the biogeochemical model can be optimized to fit the observed phosphate, dissolved inorganic carbon, oxygen, and total alkalinity distributions with reasonable fidelity for both data-assimilated and climate-model-derived circulations. However, the biological pump operates very differently for the OCIM2 and ACCESS-M circulations, largely because of differences in the sequestration time of regenerated organic matter. The differences in the biological pump in turn produce significant differences in the response of the system to imposed Southern Ocean nutrient drawdown.

2 Methods

2.1 Biogeochemical model

We model the ocean's coupled phosphorus, carbon, and oxygen cycles using mechanistic representations of nutrient uptake, particle transport, and respiration. Depending on whether oxygen is prescribed by observations or explicitly modeled, we refer to the model as the PC or PCO2 model, respectively.

Biological production is approximated as requiring only phosphate as a nutrient, and the production of organic carbon is keyed to phosphate uptake. For simplicity, dissolved organic phosphorus (DOP) is deemed to be not bioavailable, although it has been shown to be utilized in oligotrophic regions (Letscher et al., 2016). Nitrate, silicic acid, and iron limitations are not explicitly modeled either, but we do parameterize the effect of denitrification on respiration as described below. We justify our approximations a posteriori by the fidelity of the modeled tracers to observations.

We model four distinct phosphorus pools: dissolved inorganic phosphorus (DIP, which is phosphate, PO₄), semi-labile DOP, and fast and slowly sinking particulate organic phosphorus (POP_f and POP_s). The steady-state model equations for these tracers are

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& \mathcal{T}\left[\mathrm{DIP}\right]=-U_{\mathrm{P}}+R_{\mathrm{DOP}}+R_{\mathrm{POP}}+J_{\mathrm{DIP}}^{\mathrm{geo}}\\ & \mathcal{T}\left[\mathrm{DOP}\right]=\mathit{\sigma }{U}_{\mathrm{P}}-R_{\mathrm{DOP}}+D_{\mathrm{POP}}\\ & {\mathcal{S}}_{\mathrm{f}}\left[{\mathrm{POP}}_{\mathrm{f}}\right]={\mathit{\sigma }}_{\mathrm{f}}(\mathrm{1}-\mathit{\sigma })U_{\mathrm{P}}-R_{{\mathrm{POP}}_{\mathrm{f}}}-D_{{\mathrm{POP}}_{\mathrm{f}}}\\ & {\mathcal{S}}_{\mathrm{s}}\left[{\mathrm{POP}}_{\mathrm{s}}\right]=\left(\mathrm{1}-{\mathit{\sigma }}_{\mathrm{f}}\right)(\mathrm{1}-\mathit{\sigma })U_{\mathrm{P}}-R_{{\mathrm{POP}}_{\mathrm{s}}}-D_{{\mathrm{POP}}_{\mathrm{s}}},\end{array}\end{array}$$

with transport terms (circulation or gravitational settling) on the left and sources and sinks on the right. Specifically, $\mathcal{T}\left[X\right]=\mathrm{\nabla }\cdot \left(\mathit{u}\right[X\left]\right)-\mathrm{\nabla }\cdot \left(\mathbf{K}\mathrm{\nabla }\right[X\left]\right)$ is the flux divergence of dissolved tracer X due to advection (velocity u) and eddy diffusion (diffusivity tensor K). Similarly, ${\mathcal{S}}_k\left[{\mathrm{POP}}_k\right]={\partial }_z\left(w_k\right[{\mathrm{POP}}_k\left]\right)$ is the flux divergence of POP_k sinking with speed w_k (where k= f or k= s). Note that here particles are only subjected to gravitational sinking; including their advective–diffusive transport does not significantly change the solutions but greatly increases computational cost. A fraction σ of the phosphorus uptake U_P is allocated to DOP, a fraction σ_f (1−σ) is allocated to POP_f, and the remainder is allocated to POP_s, all of which are remineralized back into the DIP pool (through R_DOP and $R_{\mathrm{POP}}=R_{{\mathrm{POP}}_{\mathrm{f}}}+R_{{\mathrm{POP}}_{\mathrm{s}}}$). The global phosphate inventory is constant in our model and prescribed by weakly restoring DIP to the observed global mean $\stackrel{\mathrm{‾}}{\left[\mathrm{DIP}\right]}=\mathrm{2.17}\mathrm{\mu }\mathrm{M}$ via $J_{\mathrm{DIP}}^{\mathrm{geo}}=(\stackrel{\mathrm{‾}}{\left[\mathrm{DIP}\right]}-[\mathrm{DIP}\left]\right)/{\mathit{\tau }}_{\mathrm{geo}}$ with “geological” timescale τ_geo=1 Myr. Without prescribing the total amount of phosphate in this way there would be no unique solution to steady-state Eq. (1). (This contrasts with time-stepped models where, in the absence of external sources and sinks, the total amount of phosphate is set by the initial conditions.) Remineralization of organic phosphorus and respiration of organic carbon are modeled as having the same specific (i.e., per molecule) rates so that remineralization preserves the C : P ratio of organic-matter production (discussed below). We now briefly describe how the remineralization and respiration rates are modeled; details of the phosphorus uptake rate per unit volume U_P are provided in Appendix A1.

The remineralization of particulate organic matter (POM; either POP or particulate organic carbon, POC) is known to have relatively simple dependencies on oxygen and temperature that are parameterized following previous work (e.g., Laufkötter et al., 2017; DeVries and Weber, 2017; Dinauer et al., 2022):

$$\begin{array}{}\text{(2)}& R_{{\mathrm{POM}}_k}={\mathit{\gamma }}_k{q}_{\mathrm{10}}^{\frac{T-T_{\mathrm{ref}}}{\mathrm{10}\mathrm{K}}}{\displaystyle \frac{max\left(\left[{\mathrm{O}}_{\mathrm{2}}\right],\left[{\mathrm{O}}_{\mathrm{2}}^{\lim }\right]\right)}{max\left(\left[{\mathrm{O}}_{\mathrm{2}}\right],\left[{\mathrm{O}}_{\mathrm{2}}^{\lim }\right]\right)+K_{{\mathrm{O}}_{\mathrm{2}}}}}\left[{\mathrm{POM}}_k\right],\end{array}$$

where $T_{\mathrm{ref}}=\mathrm{20}{}^{\circ }\mathrm{C}$ and $\left[{\mathrm{O}}_{\mathrm{2}}^{\lim }\right]=\mathrm{5}\mathrm{\mu }\mathrm{M}$ defines the oxygen limit below which water is deemed anoxic. Note that Eq. (2) differs from previous parameterizations in that we include the effect of microbes switching to nitrate for organic-matter oxidization (denitrification) by explicitly disallowing respiration rates to decline in anoxic waters through the max function in Eq. (2). This means that respiration continues when [O₂] falls below $\left[{\mathrm{O}}_{\mathrm{2}}^{\lim }\right]$ at the same per molecule rate as would occur if [O₂] was equal to $\left[{\mathrm{O}}_{\mathrm{2}}^{\lim }\right]$ but without utilizing oxygen (note that [O₂] can fall below $\left[{\mathrm{O}}_{\mathrm{2}}^{\lim }\right]$ in spite of this being explicitly disallowed in the oxygen tracer Eq. (5) discussed below because we smooth step functions for differentiability as described in Appendix B3). To limit unrealistic POM accumulation in the bottom grid boxes under anoxic conditions, a small fraction of POM is dissolved into dissolved organic matter (DOM) at rate $D_{\mathrm{POM}}=\left[\mathrm{POM}\right]/{\mathit{\tau }}_{\mathrm{POM}}$ with τ_POM=1 year. Remineralization of dissolved organic matter (DOM; either DOP or dissolved organic carbon, DOC) is simply approximated as $R_{\mathrm{DOM}}=\left[\mathrm{DOM}\right]/{\mathit{\tau }}_{\mathrm{DOM}}$ with a globally uniform timescale τ_DOM=2 years.

The steady-state model equations for dissolved inorganic carbon (DIC), DOC, fast and slow POC (POC_f and POC_s), and particulate inorganic carbon (PIC, which is CaCO₃) are

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& \mathcal{T}\left[\mathrm{DIC}\right]=-\left(\mathrm{1}+r_{\mathrm{PIC}}(\mathrm{1}-\mathit{\sigma })\right)U_{\mathrm{C}}+R_{\mathrm{DOC}}\\ & +R_{\mathrm{POC}}+D_{\mathrm{PIC}}+J_{\mathrm{DIC}}^{\mathrm{atm}}\\ & \mathcal{T}\left[\mathrm{DOC}\right]=\mathit{\sigma }{U}_{\mathrm{C}}-R_{\mathrm{DOC}}+D_{\mathrm{POC}}\\ & {\mathcal{S}}_{\mathrm{f}}\left[{\mathrm{POC}}_{\mathrm{f}}\right]={\mathit{\sigma }}_{\mathrm{f}}(\mathrm{1}-\mathit{\sigma })U_{\mathrm{C}}-R_{\mathrm{POCf}}-D_{\mathrm{POCf}}\\ & {\mathcal{S}}_{\mathrm{s}}\left[{\mathrm{POC}}_{\mathrm{s}}\right]=\left(\mathrm{1}-{\mathit{\sigma }}_{\mathrm{f}}\right)(\mathrm{1}-\mathit{\sigma })U_{\mathrm{C}}-R_{\mathrm{POCs}}-D_{\mathrm{POCs}}\\ & {\mathcal{S}}_{\mathrm{PIC}}\left[\mathrm{PIC}\right]=r_{\mathrm{PIC}}(\mathrm{1}-\mathit{\sigma })U_{\mathrm{C}}-D_{\mathrm{PIC}},\end{array}\end{array}$$

where ${\mathcal{S}}_{\mathrm{PIC}}\left[\mathrm{PIC}\right]={\partial }_z\left(w_{\mathrm{PIC}}\right[\mathrm{PIC}\left]\right)$ is the flux divergence of PIC sinking with speed w_PIC. The uptake rate of carbon per unit volume $U_{\mathrm{C}}=r_{\mathrm{C}:\mathrm{P}}{U}_{\mathrm{P}}$ is keyed to phosphate uptake U_P using the stoichiometric C : P ratio r_C : P, parameterized here in terms of [DIP] (Galbraith and Martiny, 2015, see also Appendix A2). Uptake of DIC results in the production of DOC, POC_f, and POC_s in the same proportions as the corresponding phosphorus tracers (determined by σ and σ_f; see Eq. 1). For OCIM2, we account for the effect of precipitation and evaporation on DIC with “virtual fluxes” (Murnane et al., 1999) as described in the OCMIP protocol (Najjar and Orr, 1999). The ACCESS-M matrix captures the flux divergence due to water exchange with the atmosphere directly. The carbonate pump is keyed to the soft-tissue pump via the rain ratio $r_{\mathrm{PIC}}=\mathrm{PIC}:\mathrm{POC}$, and PIC dissolution is parameterized as $D_{\mathrm{PIC}}=\left[\mathrm{PIC}\right]/{\mathit{\tau }}_{\mathrm{PIC}}$ with τ_PIC=1 d.

In Eq. (3), $J_{\mathrm{DIC}}^{\mathrm{atm}}$ is the DIC source–sink term due to air–sea CO₂ exchange, parameterized in terms of surface winds and sea-ice fractions using the formulation of Wanninkhof (1992) with prescribed preindustrial atmospheric pCO₂=278 µatm (we optimize our model for preindustrial conditions, assuming negligible changes in circulation since preindustrial times). For OCIM2, we use the National Centers for Environmental Prediction (NCEP) reanalysis for the ice fraction and 6-hourly surface winds, while for ACCESS-M we use the corresponding quantities as simulated by the ACCESS climate model. From these fields, 6-hourly gas-exchange coefficients for CO₂ and O₂ are computed (to capture gustiness) and time averaged to form an annual mean climatology. The effective partial pressure of CO₂ in seawater needed for air–sea CO₂ exchange is calculated from the equilibrium carbonate chemistry using the MATLAB CO2SYS function (Lewis and Wallace, 1998; van Heuven et al., 2011).

The sinking speeds of the biogenic particles (w_s and w_f for POM and w_PIC for PIC) are constructed from globally uniform reference sinking speeds ($w_{\mathrm{s}}^{*}$, $w_{\mathrm{f}}^{*}$, and $w_{\mathrm{PIC}}^{*}$) that are multiplied with a dimensionless in situ viscosity factor α_μ to account for slower terminal velocities in colder (and to a lesser degree in more saline) waters. α_μ depends on seawater viscosity and on the density difference between POM and ambient seawater (Taucher et al., 2014, Appendix A3).

The concentration of modeled total alkalinity (TA) obeys (Murnane et al., 1999)

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\mathcal{T}\left[\mathrm{TA}\right]& =\mathrm{2}\left(D_{\mathrm{PIC}}-r_{\mathrm{PIC}}(\mathrm{1}-\mathit{\sigma })U_{\mathrm{C}}\right)\\ & -\mathrm{21.8}\left(R_{\mathrm{DOP}}+R_{{\mathrm{POP}}_{\mathrm{f}}}+R_{{\mathrm{POP}}_{\mathrm{s}}}-U_{\mathrm{P}}\right)+J_{\mathrm{TA}}^{\mathrm{geo}}.\end{array}\end{array}$$

The first term represents sources and sinks of TA due to the cycling of carbonate, and the second term contains the contributions from nitrate, phosphate, and sulfate (for further details see Wolf-Gladrow et al., 2007). We approximate the TA inventory as being constant and hence set the global mean [TA] to 2420 µM via the $J_{\mathrm{TA}}^{\mathrm{geo}}$ term, analogous to what we do for phosphate. For OCIM2, virtual fluxes are again used to account for concentration or dilution of TA by evaporation and precipitation.

The concentration of dissolved oxygen [O₂] is set by air–sea gas exchange, organic-matter respiration, and phytoplankton photosynthesis. In steady state, modeled [O₂] obeys

$$\begin{array}{}\text{(5)}& \begin{array}{rl}\mathcal{T}\left[{\mathrm{O}}_{\mathrm{2}}\right]=& r_{{\mathrm{O}}_{\mathrm{2}}:\mathrm{C}}[U_{\mathrm{C}}-\left(R_{\mathrm{DOC}}+R_{{\mathrm{POC}}_{\mathrm{f}}}+R_{{\mathrm{POC}}_{\mathrm{s}}}\right)\\ & \mathrm{\Theta }\left(\left[{\mathrm{O}}_{\mathrm{2}}\right]-\left[{\mathrm{O}}_{\mathrm{2}}^{\lim }\right]\right)]+J_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{atm}},\end{array}\end{array}$$

where $r_{{\mathrm{O}}_{\mathrm{2}}:\mathrm{C}}$ is the O₂ : C stoichiometric ratio of organic matter approximated as globally uniform. The Heaviside (step) function Θ switches oxygen consumption off when [O₂] falls below $\left[{\mathrm{O}}_{\mathrm{2}}^{\lim }\right]=\mathrm{5}\mathrm{\mu }\mathrm{M}$, consistent with the parameterization of anaerobic POM respiration in Eq. (2). The air–sea exchange rate $J_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{atm}}$ for oxygen is parameterized similarly to that for CO₂ (Wanninkhof, 1992) using the coefficients for oxygen solubility and Schmidt number as tabulated by Wanninkhof (2014).

2.2 Steady-state ocean circulation models

The nonlinear coupled partial differential Eqs. (1), (3), (4), and (5) are discretized on the model grid, and the three-dimensional tracer fields are organized into column vectors. Linear operators such as 𝒯 and 𝒮 then become sparse matrices, usually referred to as transport matrices, especially when referring to advection–diffusion. The discretized steady-state equations are coupled nonlinear algebraic equations that are solved using Newton's method, requiring order 10 iterations (Appendix B).

2.2.1 OCIM2

The Ocean Circulation Inverse Model version 2 (OCIM2; DeVries, 2014; DeVries and Holzer, 2019) provides a data-assimilated advection–eddy-diffusion transport matrix. The OCIM2 data assimilation uses the ventilation tracers CFC-11, CFC-12, radiocarbon, and ³He, in addition to sea-level height and air–sea heat and fresh-water fluxes. OCIM2 has a horizontal resolution of $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$ and 24 depth levels, with layer thicknesses that increase from 36 m at the surface to 634 m for the deepest layer. The OCIM2 transport operator is an N×N sparse matrix with $N\approx \mathrm{2}\times {\mathrm{10}}^{\mathrm{5}}$ and 3×10⁶ nonzero elements. OCIM2 arguably provides the most realistic estimate of the ocean's climatological steady-state transport and is thus a natural reference against which to assess biases in climate-model-derived transport for the current state of the ocean. For detailed analyses of the OCIM2 circulation, including ideal mean age (mean time since surface contact) and mean re-exposure time (mean time until next surface contact), see the work by DeVries and Holzer (2019).

2.2.2 ACCESS-M

As a climate-model-based estimate of the ocean's advection–diffusion operator, we use a slightly modified version of the “preferred” ACCESS1.3 transport matrix of Chamberlain et al. (2019). This matrix was built from the decadal-mean volume fluxes (resolved plus parameterized) for the 1990s from the ACCESS1.3 “historical” runs (Bi et al., 2013a), with nominal $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ horizontal resolution (finer in latitude near the Equator) and 50 depth levels with layer thicknesses that increase from 10 m at the surface to 335 m for the deepest layer. For detailed analyses of the circulation captured by the ACCESS1.3 transport matrix, including ideal mean age and mean re-exposure time, see the work by Chamberlain et al. (2019) and Holzer et al. (2020). This ACCESS1.3 matrix has a size of N×N with $N\approx \mathrm{2.7}\times {\mathrm{10}}^{\mathrm{6}}$ and 1.8×10⁷ nonzero entries, which is an order of magnitude larger than the OCIM2 matrix.

The tripolar grid of ACCESS1.3 results in a more complex sparsity pattern that slows the factorization of the Jacobian in Newton's method (Appendix B1). We therefore coarse-grain the ACCESS1.3 matrix by lumping together 2×2 nearest horizontal neighbors (similar to the “lump-and-spray” approach of Bardin et al. (2014)), which results in about 16× faster factorization. This coarse-graining reduces the maximum ideal mean age in the Pacific, but we compensate by reducing the interior background diffusivity from 0.3 to 0.1 cm² s⁻¹ to match OCIM2 and to retain the original ideal mean age. We refer to the resulting matrix/circulation model as ACCESS-M, which is of size N×N with $N\approx \mathrm{7}\times {\mathrm{10}}^{\mathrm{5}}$ and 5×10⁶ nonzero entries. We emphasize that high resolution is not important for transport matrices built from the output of an ocean model because the model's volume fluxes already contain the mean effects of processes resolved (and parameterized) at a higher resolution in the parent circulation model.

The most important difference between ACCESS-M and OCIM2 for simulating biogeochemistry stems from differences in how the mixed layer is modeled. Both matrix models use mean annual maximum mixed-layer depth (MLD). However, while OCIM2 specifies MLD from observational analyses (de Boyer Montégut et al., 2004), ACCESS-M uses the MLD of the parent ocean model. Overall, the ACCESS-M MLD is deeper than observed (roughly 1.5–3 times in the subtropical gyres) and has important unrealistic features. In the Weddell and Ross seas, the ACCESS-M MLD reaches all the way to the sea floor, and the deep winter mixed layers of the North Atlantic and Nordic seas occupy a much larger area than observed (Appendix C). The deep Southern Ocean mixed layers are due to unrealistic open-ocean convection (Bi et al., 2013a; Heuzé et al., 2013) and are a key ACCESS-M feature that imprints on the biological pump and affects its responses to perturbations (see below). Furthermore, ACCESS-M is simply built from time-averaged model volume fluxes, while OCIM2 has a steady transport that is optimized to yield propagated tracer concentrations that are as close as possible to observations. ACCESS-M therefore inherits documented circulation and thermodynamic biases from the parent ocean model (Marsland et al., 2013; Bi et al., 2013b, 2020).

2.3 Parameter optimization and tracer data

We optimize the PC and PCO2 model parameters by minimizing an objective (“cost”) function that measures the quadratic mismatch with observed DIP, DIC, O₂, and TA concentrations and penalizes deviations from a plausible range of values for each parameter. Details on the objective function and optimization procedure are provided in Appendix B4.

For DIP observations we use gridded annual mean phosphate from the World Ocean Atlas 2018 (Garcia et al., 2019). Gridded O₂, DIC, and TA observations are taken from the Global Data Analysis Project (Key et al., 2015; Lauvset et al., 2016, GLODAP v2;). We optimize our models for preindustrial conditions, assumed to be reasonably well represented by the OCIM2 and ACCESS-M circulations and by the observational DIP, TA, and O₂ climatologies. For DIC, we subtract an estimate of anthropogenic DIC as propagated from the reconstructed atmospheric CO₂ time history since 1720 using the data-assimilated OCIM2 as done by Holzer et al. (2021b). Observed tracer concentrations are interpolated onto the grid of each circulation model, and grid cells without observations are ignored in the objective function.

3 Results

We now focus on the optimized steady state of the PCO2 model and how it differs depending on whether PCO2 is embedded in the OCIM2 or ACCESS-M circulation. To examine the sensitivity of optimized model parameters to model complexity, we will also consider the PC model, for which O₂ concentrations are prescribed by observations.

3.1 Fidelity to observed fields

To quantify how well the optimized PCO2 model matches observations, we first examine the volume-weighted modeled–observed joint probability density functions (PDFs), which are essentially binned model-versus-observation scatter plots. These are shown in Fig. 1a–h together with globally averaged depth profiles (Fig. 1j–l) for DIP, O₂, DIC, and TA as obtained with either the OCIM2 or ACCESS-M circulations. Overall, there is good agreement (tight clustering of the PDFs around the 1:1 line) with volume-weighted root-mean-square errors (RMSEs) that are around 20 %–30 % of the observed standard deviation for OCIM2 and around 40 %–50 % for ACCESS-M.

Figure 1 (a–h) Joint probability density functions (PDFs) of the modeled and observed concentrations for DIP, O₂, DIC, and TA as optimized for the PCO2 model embedded in the OCIM2 (a–d) and ACCESS-M (e–h) circulations. The darker the colors the denser the PDF such that n % of the data lie outside of the nth percentile contour. The volume-weighted root-mean-square error (RMSE) is indicated in each panel along with its ratio (%) to the spatial mean and standard deviation (SD) of the observations. (i–l) Corresponding simulated and observed global-mean vertical profiles (because of interpolation, the observed profiles depend slightly on the grid used).

The optimized OCIM2-embedded PCO2 model compares well to other objectively optimized models of the P, C, and O₂ cycles (e.g., Primeau et al., 2013; Pasquier and Holzer, 2016; Holzer, 2022) as quantified by similar RMSEs (PDF panels of Fig. 1). However, unlike these other models, PCO2 has interactive oxygen providing mechanistic respiration and remineralization. Dissolved oxygen, with its large dynamic range from near zero to about 300 µM, is the tracer that has the largest mismatch with observations at an RMSE of 34 % of the spatial standard deviation from the global mean. The global mean vertical [O₂] profile matches the observations above 1000 m but progressively overestimates deeper concentrations, reaching a high bias of about 15 µM at 4000 m (Fig. 1j).

The optimized ACCESS-M-embedded PCO2 model performs worse for every tracer, with RMSEs that are larger by a factor of 1.4–2.4 (Fig. 1e–h). In contrast to the OCIM2 PCO2 model, the ACCESS-M PCO2 model underestimates oxygen at low concentrations. This could in part be due to ACCESS-M's finer low-latitude resolution, which may allow for more rapid nutrient supply, POM production, and in turn higher oxygen utilization rates. Despite the large local mismatches visible in the joint PDFs, the horizontally averaged ACCESS-M PCO2 tracer profiles fit the observations reasonably well (Fig. 1i–l), with the exception of [O₂] between 400 m to 1500 m depth where underestimates by up to 30 µM (the largest profile mismatch across models and tracers) indicate that ACCESS-M PCO2 has under oxygenated intermediate waters.

The basin zonal means of Figs. 2 and 3 show striking differences between the OCIM2 and ACCESS-M PCO2 oxygen and DIC fields. In the Southern Ocean, ACCESS-M PCO2 strongly overestimates [O₂] (by up to 80 µM) in the same region where ACCESS-M has unrealistic deep mixing (Fig. C1). This overestimate turns out to not be due to increased preformed oxygen (Fig. S3), which is similar for OCIM2 and ACCESS-M. Instead, the unrealistically deep vertical mixing dramatically reduces O₂ residence times for ACCESS-M such that total oxygen is closer to preformed oxygen in the Southern Ocean than is the case for OCIM2. This occurs despite the larger ACCESS-M respiration-rate coefficients (γ_s and γ_f) and lower q₁₀ (Table 1), presumably because of the relatively low organic-matter production in the polar Southern Ocean (Fig. 4 below). In the mid- and low-latitude Atlantic, the OCIM2 PCO2 generally overestimates oxygen especially in the oxygen minimum zones (OMZs), while ACCESS-M PCO2 underestimates oxygen, especially in the thermocline (by up to 80 µM). The underestimates of ACCESS-M PCO2 are consistent with its generally strengthened export production (discussed with Fig. 4 below), producing more organic matter and hence having higher oxygen demand than OCIM2 PCO2. We find that ACCESS-M mode and intermediate waters have a weaker preformed oxygen supply (by an order of 40 µM, Supplement Fig. S3), which also contributes to the large underestimates. In the Pacific, both OCIM2 and ACCESS-M generally underestimate O₂ in low latitudes and overestimate it elsewhere, but the underestimate for ACCESS-M is roughly twice that for OCIM2, again consistent with increased oxygen demand and under-ventilated mode and intermediate waters. In the Indian Ocean, the mismatches are similar in pattern but of larger amplitude for ACCESS-M PCO2.

Figure 2 (a–c) Basin zonal means of [O₂] from PCO2 embedded in OCIM2 for the Atlantic (a, left), Pacific (b, center), and Indian Ocean (c, right). (d–f) Model-observation difference. (g–l) As (a)–(f) but for ACCESS-M. Light gray indicates missing observations (see main text for details). For all zonal means shown in this work, the Atlantic basin excludes the Gulf of Mexico and the Caribbean, and the Pacific basin excludes the Sea of Japan so that the averages are more cleanly interpretable.

Figure 3 As Fig. 2 but for [DIC].

The zonal-mean [O₂] and [DIC] mismatches in Figs. 2 and 3 approximately mirror each other, with a Pearson correlation coefficient of about −0.6. To the extent that O₂ and DIC have realistic air–sea exchange, this anticorrelation is consistent with higher oxygen corresponding to reduced oxygen utilization and hence reduced DIC production. The details of the mismatch with observations are also influenced by errors in air–sea exchange, but the prominent mirroring of the O₂ and DIC mismatches suggests that errors in oxygen utilization are the dominant driver. Regionally in the North Pacific, the overall anticorrelation does not hold for ACCESS-M, suggesting that other factors play a role there.

When O₂ is prescribed from observations (PC model) rather than explicitly modeled, the mismatch improves for most tracers, despite oxygen not being self consistent. Specifically, we find that, relative to PCO2, the RMSEs of the PC model for DIP and DIC improve by 15 % and 3.1 % for OCIM2 and by 1.1 % and 5.8 % for ACCESS-M. For TA, the mismatch improves by 1.8 % for ACCESS-M but slightly degrades by 0.5 % for OCIM2.

3.2 Parameter sensitivity to circulation and model complexity

How sensitive are the values of the optimized biogeochemical parameters to whether we use the OCIM2 or ACCESS-M circulation, and how much do they depend on whether we prescribe the oxygen concentration from observations or model [O₂] self consistently? Recently, Kriest et al. (2020) showed that different circulations generally require different parameter values to best match observations. While Kriest et al. (2020) demonstrated this in the context of the Model of Oceanic Pelagic Stoichiometry (MOPS; Kriest and Oschlies, 2015) for five circulations, here we address the question for two other circulations (OCIM2 and ACCESS) using an entirely different model of biogeochemistry (PCO2). In addition, we investigate how sensitive optimized parameter values are to model complexity in the sense of whether oxygen is prescribed (PC) or modeled (PCO2). The optimized values of the PC and PCO2 parameters for both the OCIM2 and ACCESS-M cases are collected in Table 1.

Table 1 Biogeochemical model parameters. “Opt.” indicates if the parameter was optimized or not. Δ_bgc and Δ_circ are the sensitivities to model complexity and circulation (Appendix E).

The maximal phytoplankton concentration p_max and half-saturation constant K_DIP control nutrient uptake and are thus key for each modeled tracer. While p_max shows sensitivity to both circulation and complexity, it is more sensitive to circulation as quantified by the mean relative standard deviations of Δ_circ(p_max)=62 % and Δ_bgc(p_max)=35 %. (See Appendix E for definitions of Δ_bgc and Δ_circ.) By contrast, K_DIP is less sensitive to both circulation and complexity, lying in a range of 2.0–3.1 µM across models.

The carbonate pump is controlled by the rain ratio r_PIC, the fraction 1−σ_DOM of production allocated to POM, the sinking-speed parameter $w_{\mathrm{PIC}}^{*}$, and the PIC dissolution timescale τ_PIC. These parameters are sensitive to circulation with the OCIM2-embedded PCO2 exporting more PIC to greater depth: for OCIM2, $r_{\mathrm{PIC}}(\mathrm{1}-{\mathit{\sigma }}_{\mathrm{DOM}})=\mathrm{2.3}\mathit{\%}$ and $w_{\mathrm{PIC}}/{\mathit{\tau }}_{\mathrm{PIC}}=\mathrm{3500}\mathrm{m}$, while for ACCESS-M, $r_{\mathrm{PIC}}(\mathrm{1}-{\mathit{\sigma }}_{\mathrm{DOM}})=\mathrm{0.77}\mathit{\%}$ and $w_{\mathrm{PIC}}/{\mathit{\tau }}_{\mathrm{PIC}}=\mathrm{2100}\mathrm{m}$. The rain ratio is the most sensitive, with Δ_circ(r_PIC)=87 %.

Key to the strength of the biological pump are DOM and POM exports, which are controlled by a number of optimized parameters: the fraction σ_DOM of production allocated to DOM, the fraction σ_f of POP allocated to fast-sinking particles, and the POC respiration-rate amplitudes γ_f and γ_s, which are themselves dependent on temperature and oxygen via q₁₀, $K_{{\mathrm{O}}_{\mathrm{2}}}$, T_ref, and $\left[{\mathrm{O}}_{\mathrm{2}}^{\lim }\right]$. Each of these parameters is strongly dependent on circulation, with Δ_circ ranging roughly from 35 % to 110 % (notably γ_f and γ_s have Δ_circ≥80 %). These large sensitivities, despite identical observational constraints, show that the biological pump operates differently in the OCIM2 and ACCESS-M circulations.

Given that parameters can be expected to be least biased when optimized for the data-assimilated OCIM2 circulation, how well does ACCESS-M PCO2 match observations when solved with OCIM2-optimized parameters? With optimized OCIM2 parameters the fidelity of the ACCESS-M tracers to observations is strongly degraded. Oxygen is most affected (RMSE doubles to 64 µM) and particularly unrealistic in the Pacific where the tropical upper ocean and the old waters of the deep Pacific are strongly deoxygenated (Supplement Figs. S4–S6). Near the tropical surface this occurs because an increased fraction of production, which is largely unaffected, is routed to DOC (the OCIM2-optimized σ is more than twice the ACCESS-M-optimized σ). POM respiration, however, is shifted to a greater depth because the respiration amplitudes (γ_f and γ_s) are reduced relative to their ACCESS-M-optimized values, which allows POM to sink deeper (greater transfer efficiency). As a result of deeper POM respiration (relative to the optimized state), oxygen is stripped out of Antarctic bottom water (AABW), which greatly expands the volume of hypoxic waters in the Pacific. For DIP, DIC, and TA, we find RMSEs of 0.32, 54, and 36 µM, respectively, which is 30 %–60 % worse than the ACCESS-M-optimized fit. While non-optimal parameters by definition degrade the fit to observations, these large increases in mismatch with observations underline the central role of circulation biases.

3.3 Biological pump

Can the optimized PCO2 model robustly predict the patterns and strength of the ocean's biological pump regardless of whether we use the OCIM2 or ACCESS-M estimates of the current ocean state? To address this question, we consider a number of simple metrics of the biological pump and contrast the OCIM2 and ACCESS-M cases.

3.3.1 POC flux and export production

A commonly used metric of the biological carbon pump is the POC flux through a given depth horizon. Although the POC flux through the base of the euphotic zone is arguably more robust (Buesseler and Boyd, 2009), here we simply consider the POC flux through 100 m depth and then consider export production (referenced to the base of the euphotic zone), which is a more robust comprehensive metric of export.

Figure 4 shows maps and global zonal integrals of the 100 m POC flux. The geographic patterns of the OCIM2 and ACCESS-M 100 m POC fluxes are broadly similar, but the globally integrated flux of 22 PgC yr⁻¹ for ACCESS-M is 3 times larger than the 7.4 PgC yr⁻¹ flux for OCIM2. Relative to OCIM2, the ACCESS-M POC flux is too large in the subtropical gyres, indicating too much production fueled by excessive phosphate supply. The ACCESS-M POC flux is also larger in the subpolar oceans, particularly in the Pacific and Indian sectors of the Southern Ocean and in the North Atlantic along the Gulf Stream trajectory. These differences are likely due to the ACCESS-M model's deeper mixed layers (Fig. C1).

Because carbon can be exported in both particulate and dissolved form, a more comprehensive measure of export is the export production, that is, the rate of organic-matter production in a given euphotic water column that results in respired DIC anywhere in the aphotic ocean (Primeau et al., 2013, Appendix D). Figure 4 also shows maps and global zonal integrals of export production. Globally integrated export production, which includes export of DOM, is 16 PgC yr⁻¹ for OCIM2 and 36 PgC yr⁻¹ for ACCESS-M, considerably larger than the 100 m POC fluxes. The geographic pattern of export production is generally similar to that of the 100 m POC flux, but for OCIM2 subpolar export is much larger than low-latitude export in terms of export production than in terms of the 100 m POC flux.

3.3.2 Pump strength and regeneration pathways

A simple metric of the strength of the biological pump is the fraction of the global phosphate inventory that is regenerated, $E_{\mathrm{P}}\equiv \langle \left[{\mathrm{DIP}}_{\mathrm{reg}}\right]\rangle /\langle \left[\mathrm{DIP}\right]\rangle$, where the angle brackets denote a global volume-weighted integral. We define DIP_reg here as DIP that was remineralized in the aphotic ocean and has not been in contact with the euphotic zone since (Appendix D2). By contrast, preformed DIP is transported out of the euphotic zone without passing through the biological pump. E_P was introduced (as P^*) by Ito and Follows (2005) as a metric of pump efficiency, but regional variations of C : P in POM complicate this interpretation prompting alternate directly carbon-based metrics of pump efficiency (e.g., Holzer et al., 2021b). Remarkably, despite the biases of the ACCESS circulation, the OCIM2- and ACCESS-M-embedded PCO2 have almost identical pump strengths of E_P≈44 % and 43 %. These values are within the 39 %–50 % range of previous inverse models (DeVries et al., 2012; Primeau et al., 2013; Pasquier and Holzer, 2016; Holzer et al., 2021b) but above the original 36 % estimate by Ito and Follows (2005), which was based on apparent oxygen utilization (AOU).

Does the biological pump operate in the same way for OCIM2 and ACCESS-M? Phosphate can be regenerated through three mechanisms in our model: remineralization of POP_f, POP_s, or DOP. Similarly, DIC can be regenerated through the respiration of POC_f, POC_s, and DOC, and additionally through the dissolution of PIC (carbonate pump). To quantify the importance of each pathway, we partitioned regenerated DIP and DIC using a Green function approach (Appendix D2). The pie charts of Fig. 5 show that the dominant contribution comes from biogenic particles, accounting for roughly 73 %–78 % of regenerated DIP (and hence E_P) and 82 %–84 % of regenerated DIC for both circulations. For regenerated DIC, PIC dissolution makes a sizable contribution of 25 % for OCIM2 and 23 % for ACCESS-M, consistent with the very deep dissolution of PIC (exponential profile with e-folding length of 3490 m for OCIM2 and 2060 m for ACCESS-M).

To better understand how the biological pump sets the size of the regenerated DIP and DIC pools, it is useful to think about the bulk sequestration time of the regenerated pool and the corresponding export rates. We therefore write the regenerated inventory (for a given mechanism) as the product of the corresponding globally integrated export production and the corresponding bulk sequestration time (bulk sequestration/residence time is simply defined here as the ratio of inventory over rate and is thus equal to the regeneration-weighted mean water re-exposure time). The bulk sequestration times and corresponding export productions for each regeneration mechanism are plotted as boxes in Fig. 5, with the height of the box being the sequestration time, the length being the export production, and the area being proportional to the regenerated inventory. Despite similar POM-regenerated pools, the export production rates from POM are roughly 3× larger for ACCESS-M than for OCIM2, compensating for 3× shorter sequestration times (this is the case for all POM types, whether it be POP or POC, slow or fast). For ACCESS-M, strong export rates (wide boxes) are due to rapid uptake (large p_max) and deep mixed layers, while short sequestration times (short boxes) are due to rapid (large γ_s or γ_f) and thus shallow respiration. This is a striking example of how parameter optimization can change the inner workings of the biological pump to compensate for transport biases. We hasten to add, however, that we do not use POM measurements as a constraint on the model so that the relative contributions due to slow and fast POM are likely model specific.

The smaller optimized PIC : POC ratio for ACCESS-M (r_PIC=1.02 % compared to 6.28 % for OCIM2; Table 1) compensates for ACCESS-M's larger carbon export, resulting in ACCESS-M and OCIM2 having similar PIC exports (0.82 and 0.73 PgC yr⁻¹). We note that while the value of r_PIC=1.02 % is optimal for ACCESS-M, it is unrealistically small compared to other estimates that range from roughly 3 % to 12 % (Sarmiento et al., 2002; Jin et al., 2006; Kwon et al., 2022). The sequestration times of the PIC-regenerated DIC pools are also similar for the two circulations (670 years for OCIM2, 540 years for ACCESS-M) despite widely different PIC sinking speeds (and thus dissolution depths given the fixed dissolution timescale τ_PIC). This points to compensation due to subtle differences in the regeneration-weighted water re-exposure times. Overall, the carbonate pump contributes about a quarter of the global regenerated DIC inventory, regardless of circulation. The robustness of PIC export and PIC-regenerated DIC sequestration times and inventories across circulations is likely due to the alkalinity constraint, which acts to adjust the PIC pump to match TA observations.

Figure 4 Maps and global zonal integrals of the 100 m POC flux and the carbon export production out of the euphotic zone for PCO2 embedded in OCIM2 (a, b, c) and ACCESS-M (d, e, f). Global integrals are indicated on Asia.

Figure 5 (a) DIP_reg contributions from POP_f, POP_s, and DOP (red, orange, blue) for the PCO2 model embedded in OCIM2 represented as both a pie chart and a bar chart. Pump strength as a percentage is indicated above the pie chart. Each bar represents the regenerated DIP inventories in terms of the corresponding export rate (width) × bulk re-exposure time (height). (b) As (a) but for ACCESS-M. (c–d) As (a)–(b) but for DIC, including the additional contributions from PIC (gray).

Figure 5 shows that DOM remineralization makes a substantial contribution to the regenerated DIP and DIC inventories. The sequestration times of DOM-regenerated DIC and DIP are only a few decades, as expected given the short 2-year e-folding time for semilabile DOM in our model. The sequestration time of DOM-regenerated DIP or DIC for OCIM2 is ∼ 1.8× larger than for ACCESS-M, but its export rate is ∼ 1.5× smaller, giving roughly comparable DOM-regenerated DIC pools for both circulations. The larger DOC export for ACCESS-M is consistent with its larger nutrient and carbon uptake, in turn consistent with its deeper mixed layer supplying more nutrients. We emphasize that DOP and DOC are modeled very simply here with a single uniform lifetime τ_DOM and that we did not use any DOM observational constraints (which would require multiple DOM tracers with a spectrum of labilities). Thus, while our diagnostics demonstrate that DOM can be an important contributor to export production, the specific values of the DOM-driven export obtained here should not be considered to be accurate for the real ocean. With OCIM2, DOM accounts for roughly 50 % of the export production, while recent work places the DOM contribution to carbon export at around 20 % (Letscher et al., 2015).

3.3.3 POC transfer efficiency

The different ways in which OCIM2 and ACCESS-M PCO2 achieve optimum fits to the observations are also manifest in the models' particle dynamics, examined here in terms of the POC transfer efficiency. The efficiency of POC transfer from depth z₁ to a deeper depth z₂ is simply the ratio $\mathrm{\Phi }\left(z_{\mathrm{2}}\right)/\mathrm{\Phi }\left(z_{\mathrm{1}}\right)$, where Φ(z) is the POC flux at depth z. The transfer efficiency is a convenient and observable metric of POC flux attenuation with depth. High transfer efficiency corresponds to lower respiration rates and hence to particles surviving to greater depth.

Figure 6 shows maps of the POC transfer efficiency from the base of the euphotic zone (z₁=z_eu) to 500 m deeper ($z_{\mathrm{2}}=z_{\mathrm{eu}}+\mathrm{500}\mathrm{m}$) together with [O₂] averaged over the z₁ to z₂ water column. For OCIM2, the transfer efficiencies of both slow and fast particles have patterns that have a strong inverse correlation with the low oxygen concentrations of the Pacific OMZs. At a given temperature, respiration rates are modulated by the [O₂] Michaelis–Menten factor in Eq. (2), so that lower oxygen and respiration rates result in higher transfer efficiency, as expected. The fast-sinking POC_f achieves a 500 m transfer efficiency of 0.75 in the global mean, with local values over 0.85 in the Pacific OMZs. The slowly sinking POC_s has more time to respire over a given depth range and hence has a transfer efficiency of only 0.08 in the global mean, reaching around 0.25 in the Pacific OMZs. The transfer efficiencies are elevated by around 0.05 in high latitudes because of lower respiration in colder waters, as parameterized by the q₁₀ term in Eq. (2).

Figure 6 The transfer efficiency of POC_s (a, b) and POC_f (c, d) sinking from the base of the euphotic zone at depth z_eu to depth z_eu+500 m together with the oxygen concentration averaged over the transferred depth range (e, f) for OCIM2 (a, c, e) and ACCESS-M (b, d, f).

We note in passing that the reduced respiration in cold waters competes with increased seawater viscosity, which slows sinking particles down (smaller viscosity factor; see Appendix A3). The slower sinking allows for respiration to act over a longer period of time, compensating for the lower respiration rates. Depending on the value of q₁₀, this compensation could in principle erase the temperature dependence of respiration. However, for both models, parameters optimize such that the compensation is only partial, with the effect of reduced respiration dominating the effect of increased viscosity. The compensation is stronger for OCIM2 PCO2, for which the viscosity effect is empirically equivalent to dividing the temperature (in ^∘C) by roughly a factor of 2.4 in the q₁₀ term, compared to a corresponding factor of only about 1.3 for ACCESS-M PCO2.

The spatial patterns of the transfer efficiency for both fast and slow POC are markedly different for OCIM2 and ACCESS-M. The different patterns are a consequence of the different optimal respiration parameters $K_{{\mathrm{O}}_{\mathrm{2}}}$ and q₁₀. For both slow and fast POC, the highest transfer efficiencies for ACCESS-M occur in subpolar and polar waters because of the much greater sensitivity to temperature (about twice as large a value of q₁₀ compounded with weaker viscosity compensation). In terms of contributions to the regenerated DIC inventory, we note that the deeper POC respiration in the ACCESS-M Southern Ocean is compensated for in part by the shorter re-exposure times of about 200 years (Holzer et al., 2020) compared to up to 700 years for OCIM2 (DeVries and Holzer, 2019). For ACCESS-M, the temperature dependence dominates the oxygen dependence, with $K_{{\mathrm{O}}_{\mathrm{2}}}$ being 2.5× smaller than for OCIM2. Compared to the OCIM2 PCO2 model, oxygen in ACCESS-M must therefore drop to 2.5× lower concentrations for the same reduction in respiration, which is a bit more likely to occur because of ACCESS-M's lower OMZ oxygen concentrations (Figs. 1 and 2). As a result, ACCESS-M transfer efficiencies still show enhancement in the OMZs by about 0.3 for fast POC and only 0.03 for slow POC.

3.4 Response to Southern Ocean nutrient drawdown

Given optimal parameters for both embedding circulations, how robust is PCO2's response to perturbations? Motivated by previous studies that explored the importance of the iron-limited Southern Ocean as a source of preformed nutrients to the rest of the ocean (e.g., Sarmiento et al., 2004; Marinov et al., 2006; Primeau et al., 2013; Holzer and Primeau, 2013; Holzer et al., 2021b), we perturb the system by forcing nearly complete nutrient utilization south of 30^∘ S. This is accomplished by adding a DIP uptake rate of the form $\left[\mathrm{DIP}\right]/{\mathit{\tau }}^{*}$ with ${\mathit{\tau }}^{*}=$ 0.1 d. In the following, we contrast the ensuing responses of the OCIM2- and ACCESS-M-embedded nutrient, carbon, and oxygen cycles.

Our idealized perturbation increases carbon uptake south of 30^∘ S by 260 % for OCIM2 and by 360 % for ACCESS-M. This achieves nearly complete nutrient utilization south of 30^∘ S, which redistributes DIP (phosphate) globally because the total amount of phosphate is conserved in our formulation. Phosphate becomes “trapped” in the Southern Ocean (Primeau et al., 2013), reducing nutrient concentrations north of 30^∘ S, where biological production is reduced by 25 % for OCIM2 and by 30 % for ACCESS-M. The dramatic production increase in the Southern Ocean cranks up the global biological pump strength E_P to almost 90 % for both circulations, similar to findings of Primeau et al. (2013).

To visualize the global redistribution of nutrients, Fig. 7 shows the basin zonal averages of the DIP response. For OCIM2 PCO2, intense Southern Ocean nutrient trapping is evident with [DIP] increases of up to 1 µM at depth. The depletion of surface nutrients south of 30^∘ S deprives Antarctic Intermediate and Mode Waters (AAIW and AAMW) of the preformed DIP that they supply in the unperturbed state to the rest of the ocean (e.g., Sarmiento et al., 2004). The abyssal branch of the overturning circulation (Holzer et al., 2021a) extends elevated Southern Ocean DIP concentrations into the abyssal Pacific and Indian oceans. In the Atlantic, the nutrient trapping is more confined to high southern latitudes, with North Atlantic Deep Water (NADW) still supplying up to 0.5 µM of preformed DIP in the zonal mean (Supplement Figs. S1 and S2).

Figure 7 Atlantic, Pacific, and Indian Ocean zonal-mean [DIP] responses to complete nutrient drawdown south of 30^∘ S. (a–c) Perturbed [DIP] for OCIM2 PCO2. (d–f) Corresponding difference between perturbed [DIP] and base (unperturbed) [DIP]. (g–l) As (a)–(f) but for ACCESS-M PCO2.

For ACCESS-M, the Southern Ocean nutrient trapping is less intense than for OCIM2. The contrast with OCIM2 is particularly striking in the Atlantic sector, where increases in DIP barely reach a quarter of the OCIM2 DIP response. The reason for this contrast lies in ACCESS-M's unrealistically deep mixed layers in the Weddell and Ross seas (Appendix C), which are not present for OCIM2. The nutrient trapping in OCIM2 occurs by POP being regenerated deep in upwelling Circumpolar Deep Water (CDW), which returns DIP with the surface Ekman divergence to regions of high production where the POP flux to depth is maintained. In ACCESS-M, this sinking-POP–upwelling-DIP trapping loop is short circuited in the regions of unrealistic deep mixing. In these regions, DIP regenerated at depth is quickly mixed throughout the water column and utilized again instead of being slowly returned to the surface by upwelling CDW. Because DOP is not bioavailable in our model, phosphorus in the high-mixing regions is in effect siphoned out of the deep-mixing regions as DOP, with little remaining as DIP (outside the deep-mixing regions, regenerated DIP is trapped by the same mechanism as for OCIM2). The zonal-mean DIP depletion due to the short-circuit siphon is visible in Fig. 8g and k south of 60^∘ S in the Atlantic and Pacific, with weaker depletion in the Pacific where deep mixing occupies a smaller fraction of the basin. The weaker Southern Ocean nutrient trapping for ACCESS-M results in a correspondingly weaker response north of 60^∘ S: above 1000 m, ACCESS-M has smaller decreases in preformed DIP carried northwards by surface currents, AAMW, and AAIW, and at depth ACCESS-M has smaller increases in regenerated DIP extending northwards with the abyssal overturning of the Pacific and Indian oceans (Figs. S1 and S2).

The DIC response shown in the basin zonal means of Fig. 8 is the result of both Southern Ocean nutrient trapping and changes in air–sea CO₂ exchange. Similar to the DIP response, DIC trapped in the Southern Ocean propagates northwards at depth with AABW for both OCIM2 and ACCESS-M. The DIC response is generally larger than would be expected from the DIP response using the C : P stoichiometry of POM (which for zero DIP saturates at about 167 molC molP⁻¹; Galbraith and Martiny, 2015). The extra DIC is supplied by CO₂ ingassing driven by the strong surface DIC drawdown, which changes the Southern Ocean from net CO₂ outgassing to net CO₂ ingassing. The global DIC inventory increases by roughly 7 % for both circulations. For both OCIM2 and ACCESS-M, the Southern Ocean CO₂ ingassing weakens the decrease of preformed DIC (due to intensified uptake) that is propagated via AAMW and AAIW such that the total DIC response is dominated by the response of regenerated DIC. As for DIP, Southern Ocean DIC trapping is more pronounced for OCIM2 than for ACCESS-M, which is again a consequence of ACCESS-M's unrealistic deep mixing in the Weddell and Ross seas.

Figure 8 Atlantic, Pacific, and Indian Ocean zonal-mean [DIC] responses to complete nutrient drawdown south of 30^∘ S. (a–c) Perturbed [DIC] for OCIM2 PCO2. (d–f) Corresponding difference between perturbed [DIC] and base (unperturbed) [DIC]. (g–l) As (a)–(f) but for ACCESS-M PCO2.

The zonal mean [O₂] response quantified in Fig. 9 shows less sensitivity to the choice of circulation. For both OCIM2 and ACCESS-M, intensified Southern Ocean production dramatically deoxygenates the ocean, driving [O₂] in Southern-Ocean-sourced water masses (SAMW, AAIW, AABW) to near zero (the prominent oxygen plume that can be seen in the deep south Indian Ocean (Fig. 9c and i) is fed by CDW propagating eastward from the Atlantic). This deoxygenation is driven by strongly increased respiration, which balances the strongly increased Southern Ocean organic-matter production. Strongly increased dissolved organic matter propagates northward from the Southern Ocean with SAMW, AAIW, and AABW, shaping the oxygen response seen in Fig. 9. Increased photosynthetic oxygen production in the Southern Ocean increases [O₂] near the surface, most of which is quickly lost to the atmosphere and thus not able to meet the greater oxygen demand at depth. Outside of the Southern Ocean, oxygen weakly increases near the surface (by up to 50 µM) and in northern NADW consistent with decreased production north of 30^∘ S and decreased respiration in NADW, which also manifested in decreased DIC (Fig. 8d, j).

Figure 9 Atlantic, Pacific, and Indian Ocean zonal-mean [O₂] responses to complete nutrient drawdown south of 30^∘ S. (a–c) Perturbed [O₂] for OCIM2 PCO2. (d–f) Corresponding difference between perturbed [O₂] and base (unperturbed) [O₂]. (g–l) As (a)–(f) but for ACCESS-M PCO2.

Because the overall oxygen response is dominated by increased respiration driving [O₂] to near zero in much of the ocean, the difference in the trapping mechanisms between OCIM2 and ACCESS-M is not as pronounced in Fig. 9. For both circulations, most of the deep Pacific, Indian Ocean, and Southern Ocean become OMZs ([O₂]<20 µM), as the global ocean oxygen inventory is reduced by about 60 %. Nevertheless, for ACCESS-M, stronger oxygen decreases in the deep Southern Ocean, and weaker vertical gradients south of 60^∘ S still show the effect of the rapid deep mixing in ACCESS-M. The more rapid vertical exchange with the surface oxygen supply in the Ross and Weddell seas for ACCESS-M prevent the Atlantic and Pacific south of 60^∘ S from being as strongly deoxygenated as in OCIM2.

It is worth noting that the response to Southern Ocean nutrient drawdown is completely dominated by the circulation. Indeed, solving ACCESS-M PCO2 with OCIM2 optimal parameters results in responses that are nearly the same as those shown in Figs. 7–9.

4 Discussion

This study was motivated by the challenges posed in using ocean circulations from climate models to capture the workings of the biological pump and its effect on the ocean's oxygen distribution. In particular, how do biases in a circulation model for the current state of the ocean affect our ability to match observations, and if model parameters are optimized to match observations as well as possible, how do circulation biases affect the response of the biological pump to perturbations? The answers to these questions are important for assessing predictions for the future biogeochemical state of the ocean.

To address these issues, we built a model (PCO2) of the coupled nutrient, carbon, and oxygen cycles. The mechanistic nutrient and carbon uptake and the simpler treatment of DOC that this affords are the key differences between PCO2 and the SIMPLE-TRIM model of DeVries and Weber (2017), which otherwise share essentially the same formulation of POM respiration. The fully interactive oxygen of PCO2 is the key difference with OCMIP-style models (Najjar et al., 2007) for which POM flux-divergence profiles are prescribed and organic carbon passes through the semilabile DOC pool before being respired (e.g., Holzer, 2022).

We modeled a minimal set of biogeochemical tracers (PO4, POP, DOP, DIC, POC, DOP, PIC, O₂, TA), in part because of the greater computational demands of the ACCESS-M circulation even when coarse-grained to nominal $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$ horizontal resolution. In particular, we used only a single semilabile pool of DOC, as opposed to separate labile, semilabile, and recalcitrant pools (e.g., DeVries and Weber, 2017; Kwon et al., 2022). For this single DOC pool, remineralization is modeled using a simple fixed 2-year e-folding time because on the one hand we lack quantitative estimates of its biological and photochemical degradation and on the other hand the neglect of labile and refractory DOC pools justifies a simpler representation of semilabile DOC. By the same token, for simplicity neither DOC nor DOP are bioavailable in our model, although in the real ocean DOP provides phosphorus to P-limited phytoplankton in highly oligotrophic regions (e.g., Letscher and Moore, 2015).

The absence of an explicit nitrogen cycle means that we cannot make detailed statements about how denitrification might be affected by circulation changes, but the basic effect of organic-matter oxidization continuing in anoxic regions is parameterized. We do not model dissolved iron either, but PCO2 captures the large-scale patterns of production because uptake parameters are optimized against observed nutrient distributions, which are shaped by all the processes in the real ocean. We justify these approximations a posteriori by the high quality of the fit to the observations for the data-assimilated OCIM2 circulation.

For most parameters the variation of the optimal value with circulation is larger than the variation with model complexity (meaning prescribed versus modeled oxygen here, i.e., PC vs. PCO2), underlining the all-important control of transport on ocean biogeochemistry. Our findings also show that caution is necessary when comparing parameter values among models. Unless the circulation is free from biases and the formulation of a given process can be justified from fundamental biology and chemistry, parameter optimization is not the same as the estimation of fundamental parameters, the value of which could in principle be measured directly. Instead, optimized model parameters take on values that tend to compensate for shortcomings of the circulation and biogeochemical formulation. This generally changes the inner workings of the biological pump with export production and transfer efficiency adjusting to the circulation model's re-exposure times. While only two model complexities have been examined here (PC and PCO2), our main results should apply to more complex biogeochemical models, although our detailed quantitative findings are of course model specific.

Even when optimized with the data-assimilated OCIM2 circulation, significant biases in the biogeochemical tracers remain, pointing to model deficiencies. Remarkably, for the OCIM2 circulation, the remaining biases in the oxygen distribution are similar to those of a much simpler OCMIP-style model of oxygen also embedded in OCIM2 (Holzer, 2022). This points to potential issues with the OCIM2 circulation, air–sea exchange, and/or the parameterization of the oxygen dependence of microbial respiration. An important caveat that must be kept in mind is that the covariance between biological production, air–sea exchange, and seasonally varying circulation is not captured with our steady circulation models. In particular, the models specify the mixed-layer depth to be the climatological annual maximum, which could over-oxygenate high-latitude regions consistent with the OCIM2-optimized state having too much oxygen in the mid-depth subpolar North Pacific. Note that using a time-stepped forward model would have the benefit of capturing seasonal and inter-annual variability but would otherwise likely lead to qualitatively similar results at steeply increased computational cost.

Remaining model biases could potentially be reduced by using additional observational constraints. For example, POC observations (e.g., Dinauer et al., 2022) could be used to better constrain particle dynamics, although these are currently only available for a very sparse set of stations. Also, our results show that DOC transport is an important pathway for carbon export, suggesting potential value from using DOC observations as constraints. However, typically, total DOC is measured, not just the semi-labile fraction, so that modeling all DOC pools becomes necessary, which would increase computational cost considerably. One could also try to constrain nutrient uptake with satellite phytoplankton measurements, but this would entail using different functional classes (e.g., Pasquier and Holzer, 2017) and hence again lead to greater model complexity, and the larger set of parameters would make the optimization more costly.

The matrix formulation of PCO2 not only allowed for efficient steady-state solutions, and hence parameter optimization, but also enabled us to diagnose the inner workings of the biological pump. For example, partitioning regenerated DIP according to which mechanisms produced it is generally not computationally feasible for forward models (for forward models, regenerated DIP is typically either approximately inferred from apparent oxygen utilization (e.g., Ito and Follows, 2005) or computed as a residual from preformed DIP (e.g., Marinov et al., 2008), neither of which allow further partitioning according to remineralization mechanism). Here we were able to calculate this partitioning efficiently for steady state as detailed in Appendix D2. Similarly, export production is computationally prohibitive for forward models but readily available in the matrix formulation using an adjoint approach (Appendix D1).

5 Conclusions

To explore the effects of climate-model circulation biases on the global biological pump, we embedded a steady-state model of the ocean's nutrient, carbon, and oxygen cycles (PCO2) in the ACCESS-model-derived decadal-mean ocean circulation for the 1990s and contrasted the results with PCO2 embedded in the data-assimilated OCIM2 circulation. The differences between the OCIM2 and ACCESS cases in optimized biogeochemical parameters and in their responses to Southern Ocean nutrient drawdown lead us to the following main conclusions.

With optimized parameters, the PCO2 model is able to match the observed DIP, DIC, O₂, and TA fields with reasonable fidelity for both circulations, despite some strong biases in the ACCESS circulation. However, the fit for the ACCESS circulation is not as good as for OCIM2, with RMSEs that are roughly 1.4–2.4 times larger. Neither circulation captures all the features of the observed O₂ distribution. In OMZs, the oxygen concentration is overestimated for OCIM2 and underestimated for ACCESS-M, which points to biases in both models (possibly in both biogeochemistry and circulation) that are not compensated by parameter optimization. Modeling O₂, as compared to prescribing it from observations, increases the RMSEs for the other tracers regardless of the circulation.

The parameter values optimized for the realistic data-assimilated OCIM2 circulation are not optimal for the ACCESS-M-embedded biogeochemistry. Optimal parameter values vary by up to a factor of 7 between OCIM2 and ACCESS-M, and using OCIM2 parameters for ACCESS-M degrades the fit to observations by 30 %–60 %. This is in agreement with the findings of Kriest et al. (2020) that “one size does not fit all”. Circulation is a key control on biogeochemistry, with optimal parameter values varying more with circulation than with the complexity of the biogeochemistry model (Matear and Holloway, 1995).

Despite fitting observed tracers reasonably well, the optimized biological pump operates differently in the two circulations. This manifests in the ACCESS-M export production being roughly 2 times larger and its 100 m POC flux being 3 times larger than for OCMI2, which has a 16.4 PgC yr⁻¹ export production and a 7.4 PgC yr⁻¹ 100 m POC flux. Despite these large export differences, the biological pump strength (quantified by the regenerated fraction of the phosphate inventory) is robust at 43 %–44 % across embedding circulations. About 30 %–50 % of the global production is exported as DOC, which contributes less than ∼ 20 % of the regenerated DIC inventory for both OCIM2 and ACCESS-M. The remaining 50 %–70 % of the carbon export is carried by POC, with PIC contributing only a few percent.

Widely different exports with similar pump strengths are reconciled by differences in sequestration times (DeVries et al., 2012; Holzer et al., 2021b). We find that DOC- and POC-regenerated DIC takes roughly 3 times as long to return to the euphotic zone for OCIM2 than for ACCESS-M so that OCIM2 has a higher sequestration efficiency: a smaller export rate acts over a longer time resulting in similarly sized pools of respired carbon. For the carbonate pump (PIC), deep dissolution leads to a sequestration time of roughly 600 years and accounts for almost a quarter of the regenerated DIC inventory for both circulations.

Differences in particle dynamics shape differences in the biological pump. Globally, POC is respired deeper in OCIM2 compared to ACCESS-M, but regionally the largest differences in transfer efficiency occur in OMZs and at high latitudes through the oxygen and temperature dependence of respiration. For OCIM2, respiration is optimized to have a weak temperature dependence but a strong oxygen dependence, enhancing transfer efficiency primarily in OMZs. For ACCESS-M, temperature dominates variations in respiration, enhancing transfer efficiency mostly at high latitudes. In the ACCESS-M Southern Ocean, deeper remineralization is counteracted by much shorter deep re-exposure times (less than 200 years compared to up to 700 years for OCIM2), resulting in similar global pump strengths.

Despite PCO2 fitting observed tracers, the differences in the biological pump drive differences in the response to Southern Ocean nutrient drawdown. For OCIM2, strongly stimulated Southern Ocean production leads to intense nutrient trapping and increased carbon sequestration in the deep oceans. For ACCESS-M, Southern Ocean nutrient trapping is partially short circuited by rapid vertical mixing in the unrealistically deep mixed layers of the Weddell and Ross seas, where the intensified surface production siphons DIP out of the entire water column. The DIC response is broadly similar to the DIP response, but outside of the Southern Ocean DIC is not as depleted, and there is enhanced DIC leakage with mode/intermediate waters due to enhanced C : P stoichiometric ratios and CO₂ ingassing in the Southern Ocean. Southern Ocean nutrient drawdown leads to nearly complete deoxygenation of Southern-Ocean-sourced water masses: for both circulations, strongly increased POC and DOC production leads to sharply increased oxygen demand that cannot be met by increased ocean photosynthesis. Because [O₂] is almost driven to zero over much of the deep ocean, differences in the [O₂] responses between the two circulations are only pronounced south of ∼60^∘ S, where the rapid deep mixing in ACCESS-M provides better oxygenation.

Our findings show that optimizing biogeochemical parameters to match observed tracers does not guarantee a robust representation of the biological pump. Biases in the circulation influence how the biological pump operates and its response to perturbations, even when parameters are optimized to match biogeochemical tracer fields. It is thus imperative to quantify the inner workings of the biological pump in ocean models to assess the response of the carbon and oxygen cycles to climate change. We hope that our work will lead to future ocean models being assessed not just in terms of the fidelity of their physical variables but also in terms of key biogeochemical metrics.

Appendix A: Biogeochemistry model details

A1 Biological phosphate uptake

Following Pasquier and Holzer (2017), phosphate uptake U_P is parameterized as a function of temperature, light, and nutrient availability:

$$\begin{array}{}\text{(A1)}& U_{\mathrm{P}}\equiv {\displaystyle \frac{p_{\max }}{{\mathit{\tau }}_U}}{e}^{{\mathit{\kappa }}_{T}T}{\left({\displaystyle \frac{I}{I+K_I}}\right)}^{\mathrm{2}}{\left({\displaystyle \frac{\left[\mathrm{DIP}\right]}{\left[\mathrm{DIP}\right]+K_{\mathrm{DIP}}}}\right)}^{\mathrm{2}}.\end{array}$$

In Eq. (A1), T is the Celsius temperature and I is the photosynthetically available radiation (PAR). PAR is prescribed throughout the water column from the ACCESS1.3 model runs for both the ACCESS-M- and OCIM2-embedded PCO2 models (PAR is therefore not coupled to the plankton concentration, precluding any effects from self-shading). The main difference with the work of Pasquier and Holzer (2017) is that explicit iron and silicate limitation are not included for simplicity. Phosphate uptake is modeled to have exponential temperature dependence following Eppley (1972), who tuned ${\mathit{\kappa }}_T=\mathrm{0.063}{\mathrm{K}}^{-\mathrm{1}}$, which was later validated statistically (e.g., Bissinger et al., 2008). Light availability and nutrient limitation are parameterized as Monod factors, the square of which controls the strength of phosphate utilization following the logistic phytoplankton growth model used by Pasquier and Holzer (2017) and broadly inspired by Galbraith et al. (2010). The optimized parameter p_max represents the phytoplankton concentration for nutrient and light replete condition (unit Monod factors) and $T=\mathrm{0}{}^{\circ }\mathrm{C}$, while τ_U is a nominal uptake timescale set to 30 d. To avoid unrealistic nutrient trapping due to under-resolved circulation that occurs in some marginal seas for our models, we zero out production in the Sea of Japan, the Gulf of Mexico, and the Red Sea.

A2 Uptake C : P stoichiometry

The C : P stoichiometry of biological production has been shown to have strong regional deviations from the 106 : 1 Redfield value (e.g., Teng et al., 2014). Here we model the C : P of biological production to be identical to the C : P ratio of POM in the surface ocean, which is known to be strongly correlated with surface [DIP]. Galbraith and Martiny (2015) showed that the P : C of surface POM can be fit by the linear [DIP] dependence

$$\begin{array}{}\text{(A2)}& r_{\mathrm{P}:\mathrm{C}}=m\left[\mathrm{DIP}\right]+b,\end{array}$$

with slope $m=\mathrm{6.9}\mathrm{mmol}{\mathrm{mol}}^{-\mathrm{1}}\mathrm{\mu }{\mathrm{M}}^{-\mathrm{1}}$ and intercept $b=\mathrm{6.0}\mathrm{mmol}{\mathrm{mol}}^{-\mathrm{1}}$.

By constraining the parameters m and b of Eq. (A2) using an OCIM2-embedded inverse model of the carbon cycle, including semi-labile and recalcitrant DOP and DOC pools and a detailed representation of PIC and riverine carbon inputs, Kwon et al. (2022) recently provided an independent verification that a linear P : C dependence on [DIP] provides a good fit to observed tracers for values of m and b that are consistent with the fits of Galbraith and Martiny (2015). Phytoplankton frugality in very nutrient-poor regions has been hypothesized to be modeled better by a power-law dependence of P : C on [DIP] (Tanioka and Matsumoto, 2017; Matsumoto et al., 2020), but Kwon et al. (2022) show that their inverse model is able to fit observations equally well regardless of whether a linear or power-law relationship is used. We therefore use the simpler linear relationship of Eq. (A2) with $r_{\mathrm{C}:\mathrm{P}}=\mathrm{1}/r_{\mathrm{P}:\mathrm{C}}$ and the values of m and b from Galbraith and Martiny (2015).

A3 Viscosity effect on sinking speeds

Here, we follow a similar approach to that of Taucher et al. (2014) and define a viscosity factor α_μ that multiplies a constant reference sinking speed ($w_{\mathrm{f}}^{*}$, $w_{\mathrm{s}}^{*}$, or $w_{\mathrm{PIC}}^{*}$) to obtain the local sinking speed (of POM_f, POM_s, or PIC). The factor α_μ is given in terms of temperature T and salinity S by

$$\begin{array}{}\text{(A3)}& {\mathit{\alpha }}_{\mathit{\mu }}(S,T)={\displaystyle \frac{\mathit{\mu }(S,\mathrm{0}{}^{\circ }\mathrm{C})}{\mathit{\mu }(S,T)}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{sw}}(S,T)}{{\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{sw}}(S,\mathrm{0}{}^{\circ }\mathrm{C})}},\end{array}$$

where the first term represents the effect of dynamic viscosity μ, and the second term the effect of changing buoyancy (ρ_p is the particle density and ρ_sw is seawater density). For ρ_p we follow Taucher et al. (2014) and set it to a constant value of ${\mathit{\rho }}_{\mathrm{p}}=\mathrm{1060}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, representing an average across the literature (Logan and Hunt, 1987; Bach et al., 2012). For ρ_sw(S,T) we use the MATLAB Gibbs-SeaWater (GSW) Oceanographic Toolbox (IOC et al., 2010). For dynamic seawater viscosity, μ(S,T), we use the equation of Sharqawy et al. (2010).

Appendix B: Computational methods

B1 Steady-state solver

Equations (1), (3), (4), and (5) are discretized and collected into a nonlinear system of equations F(x)=0, where x is the concatenated vector of all the tracers. This nonlinear system is efficiently solved using Newton or quasi-Newton methods for root finding, which iteratively update the state vector via ${\mathit{x}}_{i+\mathrm{1}}={\mathit{x}}_i-{\mathbf{J}}_i^{-\mathrm{1}}{\mathit{F}}_i$, where F_i=F(x_i), and usually J_i=J(x_i) is the Jacobian of F at x_i. In practice, the Jacobian factors are not updated at every iteration to save computational resources (Kelley, 2003). Additionally, to reduce the memory required for factorization, we divide the system into smaller subsystems for P, C, and O₂ and then solve the subsystems iteratively until the entire system has converged. Specifically, we first solve the P subsystem (DIP, POP_f, POP_s, DOP), then the O₂ equation, then the C subsystem (DIC, POC_f, POC_s, DOC, PIC, TA), and repeat until F(x_i)≈0.

B2 Positivity

In practice, many equations of the PCO2 model are modified to ensure that some variables remain positive. This is useful, for example, in Monod factors such as that of O₂ in Eq. (2), to avoid catastrophic cancellation (e.g., if $\left[{\mathrm{O}}_{\mathrm{2}}\right]\approx -K_{{\mathrm{O}}_{\mathrm{2}}}$ numerically). Hence, where positivity of a variable X is required, we replace X with the differentiable approximation to $max(X,\mathrm{0})$ given by

$$\begin{array}{}\text{(B1)}& max(X,\mathrm{0})\approx X_{\mathrm{0}}\log (\mathrm{1}+e^{X/X_{\mathrm{0}}}),\end{array}$$

where X₀ is carefully chosen for every variable X to balance smoothness against accuracy (the larger X₀ the smoother but worse the approximation). Specific values used are X₀=0.1 µM for DIP, 10 µM for DIC and TA, and 1 µM for O₂.

B3 Smoothness

Because we are using Newton-type solvers to find the steady state of the tracer equations, we replace discontinuous or non-differentiable equations with smooth and differentiable approximations. For example in Eq. (5) where we abruptly turn off oxygen utilization for $\left[{\mathrm{O}}_{\mathrm{2}}\right]<\left[{\mathrm{O}}_{\mathrm{2}}^{\lim }\right]$, we approximate the Heaviside function by

$$\begin{array}{}\text{(B2)}& \mathrm{\Theta }\left(X\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}+\tanh (X/X_{\mathrm{0}})\right),\end{array}$$

where X₀ controls the smoothness, and the same values of X₀ are used as in Eq. (B1).

B4 Objective function

Our goal is to minimize the mismatch of the steady-state solution to Eqs. (1), (3), (4), and (5) with the corresponding observations (DIP, DIC, TA, and O₂). We measure the difference of tracer X with its observed values X^obs using the volume-weighted quadratic mismatch metric

$$\begin{array}{}\text{(B3)}& f_X={\displaystyle \frac{\int \mathrm{d}V{\left(X-X^{\mathrm{obs}}\right)}^{\mathrm{2}}}{\int \mathrm{d}V{\left(X^{\mathrm{obs}}-\stackrel{\mathrm{‾}}{X^{\mathrm{obs}}}\right)}^{\mathrm{2}}}},\end{array}$$

where the denominator is the spatial variance of X^obs, which provides a convenient scale for normalizing the misfit.

The objective function $\widehat{f}\left(\mathit{p}\right)$ to be minimized is then simply defined as the sum of the mismatch metrics for each constraining tracer field as

$$\begin{array}{}\text{(B4)}& \widehat{f}\left(\mathit{p}\right)=f_{\left[\mathrm{DIP}\right]}+f_{\left[\mathrm{DIC}\right]}+f_{\left[\mathrm{TA}\right]}+f_{\left[{\mathrm{O}}_{\mathrm{2}}\right]}+c.\end{array}$$

$\widehat{f}\left(\mathit{p}\right)$ is a function of the parameters p because [DIP], [DIC], [TA], and [O₂] are taken from the p-dependent steady-state solution. $\widehat{f}\left(\mathit{p}\right)$ also includes a small penalty c for the parameters, which prevents unrealistic values. In practice, we use MATLAB's unconstrained minimizer, fminunc, to find an optimal p. We note in passing that the minimum of the objective function $\widehat{f}$ determined in this way is not guaranteed to be the global minimum given the complex nature of $\widehat{f}$. However, during the course of this research, we optimized many versions of our biogeochemical model and found that they all converged to a similar minimum.

For the parameter penalty c, we prescribe strict bounds (a and b) on each optimizable parameter p, such that p remains in (a,b). We calculate the penalty as $c=\frac{\mathit{\omega }}{\mathrm{2}}{\sum }_p{\widehat{p}}^{\mathrm{2}}$, where the weight ω=0.001 ensures that the parameter penalty cost is smaller than the tracer cost, and $\widehat{p}=\log \left(\frac{p-a}{b-p}\right)$ transforms p from the interval (a,b) to $\widehat{p}$ on the interval$(-\mathrm{\infty },+\mathrm{\infty })$. The penalty for each parameter can be understood as a measure of its negative log-likelihood given a logit-normal distribution on (a,b) as its prior (with mean 0 and standard deviation 1 for its logit). The specified parameter ranges (a,b) are collected in Table B1.

Table B1 Permissible ranges for optimized parameters.

Appendix C: Key circulation characteristics

Figure C1 shows the mixed-layer depths (MLDs) used in OCIM2 and ACCESS-M. For both transport matrices a large vertical diffusivity of 0.1 m² s⁻¹ is used within the mixed layer. The MLD patterns are similar across models except at high latitudes. Most strikingly, near the Weddell and Ross seas the ACCESS-M MLDs reach the sea floor, while the observation-based MLDs used in OCIM2 are only a few hundred meters in these regions. The unrealistically deep mixed layer in the Weddell Sea was reported for the 500 years ACCESS1.3 benchmark run of Bi et al. (2013b), although that run did not have the deep mixed layer in the Ross Sea that was present in the ACCESS1.3 runs on which ACCESS-M is based.

Figure C1 Annual maximum mixed-layer depth as used in the OCIM2 (a) and ACCESS-M (b) transport matrices.

Appendix D: Biogeochemical diagnostic computations

D1 Export production

Export production via a given carbon species (POC_f, POC_s, DOC, or PIC), referred to here as a specific export pathway, is calculated using a Green function approach. Below, we detail the computation for POC_f as an example; the calculation for the other pathways is similar. We first replace the nonlinear processes with equivalent linear terms. For POC_f we have

$$\begin{array}{}\text{(D1)}& {\mathcal{S}}_{\mathrm{f}}\left[{\mathrm{POC}}_{\mathrm{f}}\right]={\mathit{\sigma }}_{\mathrm{f}}\left(\mathrm{1}-\mathit{\sigma }\right)U_{\mathrm{C}}-{\mathcal{R}}_{\mathrm{f}}\left[{\mathrm{POC}}_{\mathrm{f}}\right]-{\mathcal{D}}_{\mathrm{f}}\left[{\mathrm{POC}}_{\mathrm{f}}\right],\end{array}$$

where the local respiration and dissolution rates of Eq. (3) have been recast in terms of rate coefficients ℛ_f and 𝒟_f, diagnosed from the full nonlinear solution as ${\mathcal{R}}_{\mathrm{f}}=R_{{\mathrm{POC}}_{\mathrm{f}}}/\left[{\mathrm{POC}}_{\mathrm{f}}\right]$ and ${\mathcal{D}}_{\mathrm{f}}=D_{{\mathrm{POC}}_{\mathrm{f}}}/\left[{\mathrm{POC}}_{\mathrm{f}}\right]$. Note that with these coefficients Eqs. (D1) and (3) have the same solution. We may think of POC_f in Eq. (D1) as a linear labeling tracer that is attached to the actual POC_f and participates in nonlinear processes in proportion to the POC_f concentration. Linear labeling tracers have been very useful in a number of contexts (e.g., Holzer et al., 2014; Pasquier and Holzer, 2018; Holzer and DeVries, 2022).

Denoting ${\mathcal{A}}_{\mathrm{f}}={\mathcal{S}}_{\mathrm{f}}+{\mathcal{R}}_{\mathrm{f}}+{\mathcal{D}}_{\mathrm{f}}$, the Green function $G(\mathit{r},t|{\mathit{r}}^{\prime },t^{\prime })$ that is solution to

$$\begin{array}{}\text{(D2)}& {\partial }_{t}G+{\mathcal{A}}_{\mathrm{f}}G=\mathit{\delta }(t-t^{\prime })\mathit{\delta }(\mathit{r}-{\mathit{r}}^{\prime })\end{array}$$

gives us the POC_f contribution at location r from carbon uptake at r^′ through the convolution

$$\begin{array}{}\text{(D3)}& \left[{\mathrm{POC}}_{\mathrm{f}}\right]\left(\mathit{r}\right|{\mathit{r}}^{\prime })=\underset{-\mathrm{\infty }}{\overset{t}{\int }}\mathrm{d}{t}^{\prime }G(\mathit{r},t|{\mathit{r}}^{\prime },t^{\prime }){\mathit{\sigma }}_{\mathrm{f}}(\mathrm{1}-\mathit{\sigma })U_{\mathrm{C}}\left({\mathit{r}}^{\prime }\right).\end{array}$$

In steady-state discretized matrix form, the dt^′ integral of the Green function becomes the inverse matrix of the steady-state operator 𝒜_f so that

$$\begin{array}{}\text{(D4)}& \left[{\mathbf{POC}}_{\mathbf{f}}\right]={\mathbf{A}}_{\mathrm{f}}^{-\mathrm{1}}{\mathit{\sigma }}_{\mathrm{f}}(\mathrm{1}-\mathit{\sigma })\mathrm{diag}\left({\mathit{U}}_{\mathrm{C}}\right){\mathbf{V}}^{-\mathrm{1}},\end{array}$$

where diag(U_C) is a matrix with U_C along its diagonal and where we have multiplied on the right by V⁻¹ to obtain the contribution per unit r^′ volume.

The global export per unit volume due to production at r^′ is then obtained by integrating the respiration rate of [POC_f] over r in the aphotic domain Ω_a:

$$\begin{array}{}\text{(D5)}& {\mathit{\varphi }}_{{\mathrm{POC}}_{\mathrm{f}}}\left({\mathit{r}}^{\prime }\right)=\underset{{\mathrm{\Omega }}_{\mathrm{a}}}{\int }\mathrm{d}{\mathit{r}}^{\mathrm{3}}{\mathcal{R}}_{\mathrm{f}}\left(\mathit{r}\right)\left[{\mathrm{POC}}_{\mathrm{f}}\right]\left(\mathit{r}\right|{\mathit{r}}^{\prime }).\end{array}$$

In matrix form this becomes

$$\begin{array}{}\text{(D6)}& {\mathit{\varphi }}_{{\mathrm{POC}}_{\mathrm{f}}}={\mathbf{\Omega }}_{\mathrm{a}}^{\mathsf{T}}\mathbf{V}{\mathbf{R}}_{\mathrm{f}}{\mathbf{A}}_{\mathrm{f}}^{-\mathrm{1}}{\mathit{\sigma }}_{\mathrm{f}}(\mathrm{1}-\mathit{\sigma })\mathrm{diag}\left({\mathit{U}}_{\mathrm{C}}\right){\mathbf{V}}^{-\mathrm{1}}.\end{array}$$

For computational efficiency, we take the transpose and calculate

$$\begin{array}{}\text{(D7)}& {\mathit{\varphi }}_{{\mathrm{POC}}_{\mathrm{f}}}^{\mathsf{T}}={\mathit{\sigma }}_{\mathrm{f}}(\mathrm{1}-\mathit{\sigma }){\mathbf{V}}^{-\mathrm{1}}\mathrm{diag}\left({\mathit{U}}_{\mathrm{C}}\right){\mathbf{A}}_{\mathrm{f}}^{-\mathsf{T}}{\mathbf{R}}_{\mathrm{f}}\mathbf{V}{\mathbf{\Omega }}_{\mathrm{a}}.\end{array}$$

D2 Regenerated and preformed DIP, DIC, and O₂

Preformed concentrations are obtained by propagating euphotic concentrations into the aphotic interior without any interior sources or sinks. Preformed concentrations are thus conveniently calculated by solving

$$\begin{array}{}\text{(D8)}& \mathcal{T}\left[X_{\mathrm{pre}}\right]=\mathrm{\Theta }\left(z-z_{\mathrm{eup}}\right)\left(\left[X\right]-\left[X_{\mathrm{pre}}\right]\right)/{\mathit{\tau }}_{\mathrm{0}},\end{array}$$

where X can denote DIP, DIC, or O₂, X_pre its preformed part, and where the term on the right clamps the preformed concentration to the total concentration in the euphotic zone with rapid timescale τ₀=1 s (there is no sensitivity to the exact value of τ₀ as long as it is much smaller than the relevant transport timescales). In matrix form, Eq. (D8) becomes

$$\begin{array}{}\text{(D9)}& {\mathit{X}}_{\mathrm{pre}}={\left(\mathbf{T}+{\mathbf{M}}_{\mathrm{0}}\right)}^{-\mathrm{1}}{\mathbf{M}}_{\mathrm{0}}\mathit{X},\end{array}$$

where T is the transport matrix, X is the vector of simulated DIP, DIC, or O₂ concentrations, and M₀ is a matrix with the mask vector for the euphotic zone divided by τ₀ along the diagonal.

Conversely, regenerated tracers are obtained by labeling them at regeneration and removing them upon entry into the euphotic layer. They can thus be calculated by solving

$$\begin{array}{}\text{(D10)}& \mathcal{T}\left[X_{\mathrm{reg}}\right]=R-\mathrm{\Theta }\left(z-z_{\mathrm{eup}}\right)\left[X_{\mathrm{reg}}\right]/{\mathit{\tau }}_{\mathrm{0}},\end{array}$$

where the regenerated tracer X_reg is clamped to zero in the euphotic zone, and R is the corresponding regeneration rate per unit volume (e.g., for POC_f-regenerated DIC we use $R=R_{{\mathrm{POC}}_{\mathrm{f}}}$). In matrix form, Eq. (D10) becomes

$$\begin{array}{}\text{(D11)}& {\mathit{X}}_{\mathrm{reg}}={\left(\mathbf{T}+{\mathbf{M}}_{\mathrm{0}}\right)}^{-\mathrm{1}}\mathit{R},\end{array}$$

where R is the vector of the discretized R.

Appendix E: Sensitivity to BGC model and circulation

For a given metric or parameter X_b,c that depends on biogeochemical (BGC) model b and circulation c, we define its sensitivity to the choice of BGC model as

$$\begin{array}{}\text{(E1)}& \mathrm{\Delta }{X}_{\mathrm{bgc}}={\displaystyle \frac{\mathrm{1}}{n_{\mathrm{circ}}}}\sum _c{\displaystyle \frac{\mathrm{1}}{{\stackrel{\mathrm{‾}}{X}}_c}}\sqrt{{\displaystyle \frac{\mathrm{1}}{n_{\mathrm{bgc}}}}\sum _b{\left(X_{b,c}-{\stackrel{\mathrm{‾}}{X}}_c\right)}^{\mathrm{2}}},\end{array}$$

where n_circ=2 is the number of circulations used (OCIM2 and ACCESS-M), n_bgc=2 is the number of BGC models (PC and PCO2), and ${\stackrel{\mathrm{‾}}{X}}_c$ is the mean of X_b,c across BGC models at fixed circulation c. Sensitivity to circulation ΔX_circ is given by interchanging b and c in Eq. (E1).