Biological carbon pump estimate based on multidecadal hydrographic data

Data

Observational concentrations of DIP, DIC, ALK and O₂ were downloaded from the Global Ocean Data Analysis Project website, the second version (GLODAPv2 (ref. ¹⁵); Extended Data Fig. 1). DOC observations were from ref. ¹⁶. The data were then binned into the ocean circulation inverse model (OCIM) model grid that has a horizontal resolution of 2° × 2° and vertically 24 layers, with finer resolution in the upper ocean and coarser resolution in the deep ocean. The updated DOC compilation contains 25,869 valid data points (Extended Data Fig. 2) after binning to our model grid and has better coverage than previously widely used ones (ref. ⁵¹) that only had 14,034 valid data points in the model grid. The DOC dataset has a slight seasonal bias, with more samples collected in the summer season. However, we think that the influence is minor because: (1) a substantial proportion of the total DOC is composed of refractory DOC; unlike labile and semi-labile DOC, refractory DOC does not exhibit strong seasonality owing to its long residence time in the ocean; (2) we used tracer data from the full-water depth to constrain our model parameters. The deep ocean experiences lesser seasonal variability compared with the surface ocean. Therefore, using full-water-depth data helps to anchor the stability of the inversion. Two NPP products, carbon-based NPP from Sea-viewing Wide Field-of-view Sensor (SeaWiFS CbPM)⁵² and CAFE, were downloaded from https://doi.org/10.6084/m9.figshare.19074521. The NPP products were interpolated and averaged by Nowicki et al.²⁴ to the same model grid as used in this study. The climatological ocean temperature and silicate are from World Ocean Atlas 2018 (refs. ^53,54). The projected temperature at 2099 was obtained from a CESM-BGC model prediction under the RCP8.5 scenario⁵⁵. The historical atmospheric pCO₂ data were obtained from ref. ⁵⁶ for the period from 1850 to 2015 and were downloaded from https://scrippsco2.ucsd.edu/data/atmospheric_co2/primary_mlo_co2_record.html (ref. ⁵⁷) for the period from 2016 to 2020.

Biogeochemical inverse model

A schematic of the structure of the biogeochemical model is shown in Fig. 1. The model couples the cycling of phosphorus (P), carbon (C) and oxygen (O). The phosphorus model is the base model that provides a biological uptake rate (γ(r), in which r is a position coordinate) in P units (G ≡ (γ[DIP])), which is then converted to a DIC uptake rate in the carbon model by incorporating a C:P ratio (r_C:P). In the P-cycle model, the DIP assimilation rate is modelled using a spatial pattern obtained from satellite-derived NPP (mg C m⁻² day⁻¹) and a gridded surface DIP climatology as follows

in which NPP₀ and [DIP]₀ are 1 mmol C m⁻² day⁻¹ and 1 μM that are functioned to remove dimensions of NPP and DIP; α and β are adjustable parameters that are constrained in the inversion; r_C:P is the C:P ratio that is used to convert NPP from C unit to P unit and modelled according to Galbraith and Martiny⁵⁸ (r_C:P = (0.006 + 0.0069[DIP]_obs)⁻¹). z and z_C are water depth and the euphotic zone depth, respectively. Photosynthesis is assumed to occur only in the euphotic zone and to be zero below. The euphotic zone is defined as the top two model layers (73.4 m).

Phosphorus model

The phosphorus model considers four explicit tracers: dissolved inorganic phosphorus (DIP), dissolved semi-labile organic phosphorus (DOP), dissolved labile organic phosphorus (DOP_l) and particulate organic phosphorus (POP). We assign an e-folding remineralization time (1/κ_l) of 12 h for DOP_l so that it quickly cycles in the upper ocean, with little being transported below the euphotic zone. We use a parameter δ to allocate production to labile pools. The remaining production (total production less production to DOP_l) is allocated to DOP and POP. The factions σ_P and (1 − σ_P − δ) of the production allocated, respectively, to DOP and POP are determined by estimating the parameter σ_P through our Bayesian inversion procedure. The advective-diffusive transport of dissolved tracers (DIP, DOP and DOP_l in the P model; DIC, semi-labile dissolved organic carbon (DOC), labile dissolved organic carbon (DOC_l), refractory dissolved organic carbon (DOC_r) and ALK in the C model; and O₂ in the O model) is computed using the OCIM tracer transport matrix, \({\bf{T}}\left[{\rm{C}}\right]\equiv \nabla \cdot \left(\vec{U}\left[{\rm{C}}\right]-{\rm{K}}\nabla \left[{\rm{C}}\right]\right)\), in which \(\vec{U}\) is the velocity vector and K is the diffusive term. Τ represents the climatological mean circulation of the ocean. The OCIM tracer transport matrix is constrained using salinity, temperature, sea-surface height, CFC-11, CFC-12, ¹⁴C, ³He etc. (see DeVries and Holzer⁵⁹ for details). We neglect the advective-diffusive transport of particulate tracers (POP in the P model and PIC and POC in the C model) so that particulate tracers are transported only vertically. The vertical transport of POP is modelled using a sinking flux divergence operator (\({{\bf{F}}}_{{\rm{POP}}}\equiv \nabla \cdot (\vec{w}[{\rm{POP}}])\)), in which \(\vec{w}\) is the sinking speed of POP and is directed downward. We choose a sinking speed that increases linearly with depth and a constant dissolution rate, κ_P = (1/30) days⁻¹, so that the attenuation of the vertical flux of POP follows a power-law function, F(z) = F(z₀)(z/z₀)^−b, in which F(z) and F(z₀) are fluxes at a depth of z and z₀, respectively⁶⁰. A sensitivity test with κ_P = (1/60) days⁻¹ suggests that the choice of κ_P does not markedly influence our results. The exponent b for the P model (C model in the following section) is defined in the following way (ref. ⁴²), b(P) = b_PθT + b_P, in which b_Pθ and b_P are two adjustable parameters and T is the average temperature of the model euphotic zone. The initial guess of b_Pθ is set to zero, thereby avoiding any intentional imposition of temperature dependence. The optimization process determines both the sign and magnitude of b_Pθ. The governing equations for the phosphorus cycle are as follows:

in which κ_dP is the DOP remineralization rate constant that is a function of temperature defined using a Q₁₀ function (κ_dP = κ_PθQ₁₀^(T−30)/10), in which T is water temperature from World Ocean Atlas 2018 (ref. ⁵⁴). κ_Pθ and Q₁₀ are optimized in the inversion. κ_l is the e-folding remineralization time of DOP_l, which is fixed at κ_l = (1/12) h⁻¹. We tested the sensitivity to a smaller κ_l = (1/24) h⁻¹ and found that the choice of κ_l did not substantially change the fittings to the tracers but could alter the export flux of labile organic matter. We therefore include different κ_l values in the uncertainty analysis (see the ‘Uncertainty analysis’ section). κ_g is prescribed to (1/10⁶) years⁻¹ and is used to set the global mean DIP concentration to the observed global mean concentration ([DIP]_obs). κ_P is a prescribed POP remineralization rate constant (κ_P = (1/30) days⁻¹). A sensitivity test shows that increases or decreases in the fraction of DOP_l production (δ in equations (2) and (3)) does not alter the fit to the observational data nor does it change the inferred export fluxes of POC and semi-labile DOC. We therefore set δ to be zero in the first-round optimization.

Carbon model

The carbon model explicitly simulates seven tracers: DIC, DOC_l, DOC, DOC_r, POC, PIC and ALK (Fig. 1b). The DIP assimilation rate G is converted to the DIC assimilation rate by incorporating a C:P ratio (r_C:P) that is allowed to vary spatially according to the modelled DIP concentration, r_C:P = (cc[DIP] + dd)⁻¹, in which cc and dd are estimated as part of the inversion. As in the P model, we set the allocation to DOC_l to be zero in the first round of optimization. Subsequently, we prescribe the difference between satellite NPP and model organic carbon production as the production for labile DOC_l, so that our model production matches satellite NPP exactly. The fraction, σ_C, of the organic carbon production allocated to POC and DOC pools is estimated as part of the inversion and does not need to be the same as the fraction σ_P allocated to the POP and DOP pools. A further adjustable parameter, η, is used to control the fraction of DOC that is transferred to the refractory pool by bacterial reworking. The remaining DOC fraction (1 − η) is remineralized back to DIC. The e-folding decay times of DOC_r (κ_ur and κ_dr for the upper and deeper ocean, respectively) are estimated as part of the inversion. POC sinks and is gradually remineralized to DIC in the water column. The downward transport of POC is modelled using a flux divergence operator (F_POC), which is formulated in the same way as the POP sinking flux-divergence operator F_POP with independent adjustable parameters b_Cθ and b_C that are determined as part of the inversion (Extended Data Table 1). Unlike DIP, DIC experiences sea-to-air gas exchange at the surface. This gas exchange is modelled according to the method used for phase 2 of the Ocean Carbon-Cycle Model Intercomparison Project (OCMIP-2)⁶¹ using a recalibrated piston velocity (see the next section). Also, freshwater precipitation and evaporation can greatly affect surface ocean DIC and ALK concentrations. Precipitation will dilute, whereas evaporation will concentrate their concentrations. A virtual flux according to OCMIP-2 (ref. ⁶¹) is applied to model for the effects of precipitation and evaporation on DIC and ALK (F_vDIC[DIC]_s and F_vALK[ALK]_s, in which [DIC]_s and [ALK]_s are the mean surface-ocean concentrations of DIC and ALK, respectively).

Production of PIC is modelled to be proportional to the production of POC using two adjustable parameters, r_Si and r_RR, that are estimated in the inversion. The parameter r_Si adjusts PIC production according to silicate concentration in the surface ocean in linear form (R_RR = r_Si[SiO₄⁴⁻] + r_RR). The downward transport of PIC is modelled using a flux divergence operator (F_PIC), which generates a PIC flux profile that follows an exponential function F_PIC(z) = F₀exp((z − z₀)/d), in which d is the PIC dissolution length scale, whose value is estimated as part of the inversion (Extended Data Table 1). Compared with a power-law function, an exponential function with a length scale on the order of several thousand metres leads to a much smaller CaCO₃ dissolution rate in the shallow water in which CaCO₃ is supersaturated⁶². Every mole of PIC production consumes two moles of ALK. By contrast, the dissolution of one mole of PIC releases two moles of ALK (equation (3)). From the perspective of carbon, photosynthesis and remineralization of organic matter do not change alkalinity. However, in the processes of photosynthesis and remineralization, chemical forms of nitrogen change, which influences alkalinity so that a mole of organic carbon production increases alkalinity by r_N:C moles, whereas a mole of organic carbon remineralization decreases alkalinity by r_N:C moles. The governing equations for carbon cycling are as follows:

Anthropogenic DIC

To use DIC observations to constrain our inverse model, we have to take into account the changing DIC concentration owing to the invasion of anthropogenic CO₂ into the ocean. To obtain a self-consistent estimate of the anthropogenic carbon signal, we performed a time-dependent simulation using equation (3). Starting from an assumed steady state, we time-stepped our carbon-cycle model forward in time from 1850 to 2020, using an implicit trapezoid-rule time-integration scheme for all terms except for the gas exchange, for which we used an explicit Euler forward scheme. In this calculation, we prescribed the surface SST according to a time-dependent reanalysis product (ref. ⁶³). The transient integration was carried out with a time step size of Δt = 2 months. The atmospheric pCO₂ was prescribed according to ref. ⁵⁶ from 1850 to 2015 and according to ref. ⁵⁷ from 2016 to 2020. We also simulated δ¹⁴C to better calibrate the air–sea gas-exchange velocity as described below. The atmospheric δ¹⁴C was prescribed according to ref. ⁶⁴ for the period from 1850 to 2015 and according to ref. ⁶⁵ from 2016 to 2020. To produce the initial conditions, we assumed that the system was in steady state in 1850 and used Newton’s method to find the steady state.

To calibrate the air–sea gas exchange parameterization, we re-optimized the scaling factor in the OCIM2 gas-exchange scheme by minimizing the misfit between our modelled δ¹⁴C and the GLODAPv2 δ¹⁴C data. See Extended Data Fig. 5a,b for the number of observations as a function of time. To compute the misfit, we sampled our model at the location and times of the bottle measurements in the GLODAPv2 database. Our calibration method followed an iterative two-step process in which we first optimized the air–sea gas exchange through a series of transient carbon-cycle simulations. After obtaining the optimal air–sea gas exchange, we subtracted the excess anthropogenic DIC from the GLODAPv2 measurements to produce an estimate of the natural background DIC for the year 1850. The resulting DIC data and optimal gas-exchange velocity were then used for the optimization of the biogeochemistry model (see the ‘Parameter estimation’ section). The optimized biogeochemical model was then used to produce an updated initial condition for the transient carbon-cycle simulation and a re-optimization of the air–sea gas-exchange velocity. We repeated this two-step process until we obtained self-consistent estimates of: (1) the optimal biogeochemical parameter values (Extended Data Table 1); (2) the biogeochemical state; (3) the scaling factor for the air–sea gas–transfer velocity, a = 0.234 cm h⁻¹ (m s⁻¹)⁻²; (4) transient DIC; and (5) the transient δ¹⁴C signal including the combined effects of radioactive decay, the Suess effect and the bomb radiocarbon signal. Extended Data Figure 5c shows a time series of the excess anthropogenic DIC concentration averaged over the top 100 m of the water column and for the water column below 100 m. By 2020, the vertical DIC gradient is reduced by 20%.

Oxygen model

Oxygen production is modelled by applying a ratio of oxygen to carbon (r_O:C) to the DIC assimilation rate (Gr_C:P). The ratio r_O:C is optimized in the process of inversion. We convert the DOC and POC remineralization rates (ηκ_dC[DOC] + κ_r[DOC_r] + κ_l[DOC_l] + κ_p[POC]) to an oxygen consumption rate using the same r_O:C ratio and gradually shut down oxygen consumption as the oxygen concentration falls below the critical value (O_crit = 5 mmol l⁻¹) using a hyperbolic equation (R([O₂]) = 0.5 + 0.5tanh[([O₂] − O_crit)/[O₂]₀]), in which [O₂]₀ (1 mmol l⁻¹) is used to remove the O₂ dimension. Sea-to-air O₂ flux (F_O2) is modelled according to OCMIP-2 (ref. ⁶¹):

in which the matrix R is a diagonal matrix, whose elements are given by R([O₂]).

Parameter estimation

The 21 adjustable parameters of the model (Extended Data Table 1) were estimated using a Bayesian inversion method. In this approach, the solutions to our model equations define the tracer fields as implicit functions of the adjustable parameters, which we then compare with the observations to construct a likelihood function. We obtain the P, C and O fields by finding the steady-state solutions of the governing equations for the P, C and O models (equations (2)–(4)). Because the governing equations for the P model are linear, their steady-state solution can be obtained efficiently by direct matrix inversion after setting the time derivatives in equation (2) to zero. We fix the atmospheric CO₂ concentration at the preindustrial level (278 ppm) to compute the preindustrial sea-to-air CO₂ flux (F_CO2). The steady-state solution for the C model is solved using Newton’s method because of nonlinearity in F_CO2. The governing equation for O is also nonlinear because of the hyperbolic function (R) that turns off oxygen consumption when oxygen concentration is critically low. We solve the oxygen equations using Newton’s method.

To find the most probable parameter values, we minimize the negative logarithm of the posterior probability function, which is equivalent to minimizing the negative log-likelihood because we log-transformed our parameters (except of the slopes of exponent b) so that they have flat priors:

in which the e_X are column vectors whose elements are given by the difference between modelled and observed concentrations, e_X = H_X[X_mod] − [X_obs], in which the X label denotes the specific tracer and H_x is a rectangular matrix that picks out the model grid boxes that have observations of tracer X. Because there are no measurements that precisely separate DOC into different pools according to their lability, we sum all three pools in the model (DOC, DOC_l and DOC_r) and compare the sum to observations. For DIC, we subtracted our estimated anthropogenic DIC from the bottle measurements in the GLODAPv2 database according to the location and time of measurement (see the ‘Anthropogenic DIC’ section). W_X is a precision matrix for tracer X and is defined in the following way:

in which V_X is a diagonal matrix with the fractional volumes of the model grid boxes (V = diag(ΔV_i/Σ_iΔV_i), in which the subscript i is the index of the grid boxes that have at least one observation) and \({\sigma }_{{\rm{x}}}^{2}\) is the spatial variance of the observations, that is,

with the spatial mean given by

in which V_X is a diagonal matrix with the grid-box volumes and the subscript X represents the grid boxes that have observations of tracer X. The bold 1 represents a column vector. The transpose turns it into a row vector. Thus, the numerator yields the volume integral of X_obs and the denominator yields the total volume.

The optimization is conducted using MATLAB’s fminunc function, which is computationally efficient because we can supply hand-coded first and second derivatives of the objective function with respect to the adjustable parameters. The optimization generally takes fewer than 100 iterations. The most probable model parameter values are presented in Extended Data Table 1. Parameter error bars that correspond to ±1 standard deviations are calculated using Laplace’s approximation as described in ref. ⁶⁶.

Calculation of carbon flux

The two-dimensional vertical-flux field (f_POC) is calculated by vertically integrating POC remineralization below the euphotic zone (\({f}_{{\rm{P}}{\rm{O}}{\rm{C}}}={\sum }_{i=1}^{z}{\kappa }_{{\rm{p}}}{\rm{P}}{\rm{O}}{{\rm{C}}}_{i}{\Delta V}_{i}{M}_{i}\), in which i represents the index for deep grid, ΔV_i is the volume of the ith grid box and M_i is a mask that is set to 1 below the euphotic zone and 0 elsewhere). In our model, POC is transported vertically in the water column and is not advected to neighbouring grids. This approximation is appropriate for the coarse horizontal resolution of our model. The non-advective-diffusive flux below the first two layers is calculated on the basis of the following power-law function (also known as the Martin curve function), f_POC(z) = f_POC(z₀)(z/z₀)^−b, in which z₀ is the euphotic zone depth and z is the depth at which non-advective-diffusive flux is calculated. The exponent, b, depends linearly on the surface water temperature⁴². The estimated slope and intercept (b_C and b_Cθ) for this linear relationship are presented in Extended Data Table 1. The advective-diffusive fluxes of labile and semi-labile organic carbon are calculated using an adjoint method, which tracks the export and subsequent remineralization of DOC, as described by ref. ³¹. Only DOC respired below the depth of the euphotic zone is counted as export. The flux of TOC is the sum of the non-advective-diffusive flux and fluxes from labile DOC and semi-labile DOC. We ignore the export of refractory DOC because of its negligible contribution.

To compare the non-advective-diffusive flux to CMIP6 models at their consensus reference depth (about 100 m), we scale our estimated non-advective-diffusive flux using a power-law function with our optimized temperature-dependent b exponents.

To compare our export estimates to the geochemical ANCP estimates that are calculated at the base of spatially varying mMLDs obtained from a CESM simulation³⁵, we estimated the export fluxes of POC and semi-labile DOC to the depth z = mMLD. For cases in which mMLD is deeper than our euphotic zone depth, we scaled the fluxes using the power-law function with our optimized b exponents. For cases in which mMLD is above the base of the euphotic zone, we did not apply the power-law scaling because it tends to amplify errors. In those cases, we used the export flux at the base of the model’s euphotic zone for the comparison. Note that the contribution of labile DOC is ignored when scaling flux down to mMLD owing to the short e-folding decay time (12 h or 24 h).

Uncertainty analysis

The uncertainty analysis is conducted in two ways. First, we use a Monte Carlo method whereby an ensemble of parameter values is drawn from a multivariate normal distribution whose mean is given by our estimated most probable parameter values and whose covariance matrix is given by the inverse of the matrix of second partial derivatives of the negative logarithm of the posterior probability distribution, that is, by the Hessian matrix. For each ensemble member, we solve the steady-state model equations and calculate the organic carbon fluxes. However, the parameters are so well constrained that their uncertainties are small, and the flux uncertainties calculated this way are small. Second, because the DIP uptake model is constructed with two different satellite NPP products (SeaWiFS CbPM and CAFE), and the 21 adjustable parameters are optimized for each NPP field (Extended Data Table 1), the influence of NPP fields on export fluxes are much larger than that of parameter uncertainties. Also, the e-folding remineralization time of labile DOC is prescribed at 12 h and 24 h. We, therefore, report flux uncertainties estimated from the results based on different initial NPP fields and on different labile DOC e-folding decay timescales. The distributions of the standard deviation of key outputs are illustrated in Extended Data Fig. 9.

DOC sequestration time

To calculate the DOC sequestration time, we injected unit DOC pulses in the model euphotic zone and tracked this DOC as it was transported by the circulation, respired into DIC (according to the timescale given in Extended Data Table 1) and then transported back to the surface, where it was rapidly removed with a loss frequency of (1 day)⁻¹. We then spatially integrated the removal rate for each DOC pulse to obtain residence-time distributions for the DOC exported from the surface of each water column.

Sequestration-time-partitioned distribution functions

To compute the sequestration-time-partitioned distribution functions, we use the three-dimensional organic carbon respiration rate to construct a Dirac δ-function pulse of labelled regenerated inorganic carbon. The resulting tracer field is then transported using the circulation model until it is removed in the 36.1-m-thick surface layer of the model using a loss frequency of (1/500) year⁻¹. We integrate the system forward in time for 10,000 years, by which time all of the regenerated-carbon pulse has left the system. We use a second-order-accurate trapezoidal integration rule starting with a time-step size of less than 10⁻⁴ years and gradually increase it to 10 years by the end of the simulation. A sequestration-time density distribution function is obtained by globally integrating the loss rate and the cumulative distribution function is then obtained by integrating the density function for progressively longer times. To obtain the cumulative sequestration-time distribution for the stock of regenerated DIC, we first integrate the tracer field over the whole volume of the ocean and then integrate the resulting stock for progressively longer sequestration times. By year 10,000, the resulting integral is equal to the global inventory of regenerated DIC.