Sensitivity of atmospheric CO₂ to regional variability in particulate organic matter remineralization depths

The concentration of CO₂ in the atmosphere is sensitive to changes in the depth at which sinking particulate organic matter is remineralized: often described as a change in the exponent “b” of the Martin curve. Sediment trap observations from deep and intermediate depths suggest there is a spatially heterogeneous pattern of b, particularly varying with latitude, but disagree over the exact spatial patterns. Here we use a biogeochemical model of the phosphorus cycle coupled with a steady-state representation of ocean circulation to explore the sensitivity of preformed phosphate and atmospheric CO₂ to spatial variability in remineralization depths. A Latin hypercube sampling method is used to simultaneously vary the Martin curve independently within 15 different regions, as a basis for a regression-based analysis used to derive a quantitative measure of sensitivity. Approximately 30 % of the sensitivity of atmospheric CO₂ to changes in remineralization depths is driven by changes in the subantarctic region (36 to 60^∘ S) similar in magnitude to the Pacific basin despite the much smaller area and lower export production. Overall, the absolute magnitude of sensitivity is controlled by export production, but the relative spatial patterns in sensitivity are predominantly constrained by ocean circulation pathways. The high sensitivity in the subantarctic regions is driven by a combination of high export production and the high connectivity of these regions to regions important for the export of preformed nutrients such as the Southern Ocean and North Atlantic. Overall, regionally varying remineralization depths contribute to variability in CO₂ of between around 5 and 15 ppm, relative to a global mean change in remineralization depth. Future changes in the environmental and ecological drivers of remineralization, such as temperature and ocean acidification, are expected to be most significant in the high latitudes where CO₂ sensitivity to remineralization is also highest. The importance of ocean circulation pathways to the high sensitivity in subantarctic regions also has significance for past climates given the importance of circulation changes in the Southern Ocean.

1 Introduction

Sinking particles of organic matter transfer 5–10 Pg C per year from the upper ocean to the ocean interior (Henson et al., 2011) as part of a process known as the biological pump. As these particles sink, they are remineralized through bacterial- and zooplankton-related activity, releasing the carbon and nutrients back into solution at depth. Vertical fluxes of particulate organic carbon (POC) in the water column have historically been described by the Martin curve, a power-law function that describes the rapid decrease in flux (F_z) from a maximum value at depth z₀, nominally the base of the mixed layer, to a small asymptotic value in deep waters (Fig. 1; Eq. 1; Martin et al., 1987):

$$\begin{array}{}\text{(1)}& F_z=F_{z_{\mathrm{0}}}{\left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)}^{-b}.\end{array}$$

The dimensionless exponent in the power law (b) describes whether organic matter is remineralized predominantly at shallower depths (larger values of b; e.g. b=1.6) or deeper in the water column (smaller values of b; e.g. b=0.4) (Fig. 1). The exponent itself parameterizes the rate at which POC sinks through the water column (units of m d⁻¹) and the rate at which it is remineralized (units of d⁻¹; Kriest and Oschlies, 2008; Lam et al., 2011). In this paper we use the term “remineralization depth”, defined as a depth at which a defined % of POC has been remineralized. Previously, this has been defined as an e folding depth: the depth at which ∼63 % of POC has been remineralized (Kwon et al., 2009; although note that the Martin curve is not exponential).

Figure 1 The normalized water column distribution of particulate fluxes defined using the Martin curve. As a comparison, Martin curves are shown with the exponent found by Martin et al. (1987) (b=0.858) and minimum and maximum exponent values used in this study based on sediment trap data compilations (b=0.4, b=1.6; Henson et al., 2012; Marsay et al., 2015; Guidi et al., 2015). All curves have export depth (z₀, Eq. 1) of 120 m.

Ocean biogeochemical models predict that the concentration of CO₂ in the atmosphere is sensitive to changes in a globally uniform remineralization depth. Kwon et al. (2009) showed that a deepening of the remineralization depth globally of 24 m (from b=1.0 to 0.9) redistributed dissolved inorganic carbon (DIC) from the intermediate waters to the deep ocean, leading to a reduction in atmospheric CO₂ of between 10 and 27 ppm. The drawdown was also associated with a decrease in the global mean concentration of preformed nutrients in the ocean interior (nutrients that are not utilised by biology in the surface ocean and enter the ocean interior via circulation; Ito and Follows, 2005). Kwon et al. (2009) found that an increase in respired carbon in the deep ocean was balanced by a reduction in preformed nutrients exported in the North Atlantic. Deepening of the POC remineralization depth could also drive dissolution of calcium carbonate (CaCO₃) in ocean sediments, ultimately drawing down more CO₂ over millennial timescales (Roth et al., 2014). The potential impact of remineralization depth changes on atmospheric CO₂ is therefore a highly relevant component of the marine carbon cycle for both past and current changes in climate (Riebesell et al., 2009; Hülse et al., 2017; Meyer et al., 2016).

Analyses of global sediment trap observations suggest there is a spatially heterogeneous pattern of remineralization depths in the modern ocean that varies particularly with latitude. A synthesis of observations from deep sediment traps (>1500–2000 m; Henson et al., 2012) suggests that POC fluxes in high latitudes attenuate faster with depth (shallower remineralization depth: b=1.6) than in low latitudes, where a greater proportion of POC is transported to depth (deeper remineralization depth: b=0.4); see Fig. 1). However, POC fluxes measured using neutrally buoyant sediment traps at shallower depths (<1000 m) suggest the inverse of this latitudinal pattern (Marsay et al., 2015); see also Weber et al. (2016). A recent compilation of sediment trap data and profiles of particle size distributions observed in the water column highlight additional intra-basin variability in b (e.g. shallower remineralization in the eastern equatorial Pacific than in the western equatorial Pacific) and inter-basin variability (e.g. deeper remineralization in the Atlantic and Indian basins compared to the Pacific) (Guidi et al., 2015). The uncertainty in the spatial variability of remineralization depths presents a challenge for determining which mechanisms may be responsible for changes in remineralization depths and how these might drive future or past changes in remineralization (e.g. Boyd, 2015). Additionally, this also presents a challenge for biogeochemical models, which are beginning to resolve the mechanisms that are potentially responsible for these spatial patterns, such as particle-size-dependent sinking rates (DeVries et al., 2014), temperature-dependent remineralization (John et al., 2014), and oxygen dependence (Laufkötter et al., 2017).

A key question in light of the observed spatial variability in remineralization depths and the associated uncertainty in spatial patterns is what the sensitivity of atmospheric CO₂ concentrations to spatial variability in remineralization depths is. Kwon et al. (2009) further quantified the sensitivity of atmospheric CO₂ to basin-scale changes in remineralization depths by perturbing them in each basin individually, finding that the Pacific, Southern Ocean (defined as >40^∘ S), Atlantic, and Indian oceans contributed 38 %, 22 %, 21 %, and 19 % of the total CO₂ drawdown, respectively (Kwon et al., 2009). The variability in CO₂ sensitivity between basins matched the variability in the magnitude of export production integrated over the basins and basin area, suggesting that no one region was more significant when varying the globally uniform remineralization depth (Kwon et al., 2009). However, this basin-scale analysis does not resolve the sensitivity of atmospheric CO₂ occurring at the resolution suggested by observations, i.e. a latitudinal and within-basin scale, or at the resolution of ecological and biogeochemical variability (Longhurst, 1998; Fay and McKinley, 2014). Additionally, the analysis does not allow for the identification of potential interactions and feedbacks between regions when remineralization depths are changing simultaneously.

Here we aim to address these issues by performing a global sensitivity analysis of regionally varying remineralization depths. To this end, we use a transport matrix (a steady-state, computationally efficient representation of ocean transport) derived from the MIT global circulation model (MITgcm) with a model of phosphorus and carbon cycling where the ocean is divided into 15 regions in which remineralization depths can change independently. Remineralization depths are perturbed simultaneously using Latin hypercube sampling and sensitivity quantified using regression analysis,.

2 Methods

2.1 Model description

We provide a brief description of the model here and a more detailed description in Appendix A. The approach to quantifying sensitivity used here relies on the ability to run an ensemble of model experiments. To make this approach feasible we use the “transport matrix method” (Khatiwala et al., 2005; Khatiwala, 2007), a steady-state, computationally efficient representation of ocean transport and climate derived from a global circulation model. We use monthly mean transport matrices derived from the 2.8^∘ global configuration of the MIT global circulation model (MITgcm) with 15 vertical levels (Khatiwala et al., 2005; Khatiwala, 2007). These specific matrices have been previously applied to model biogeochemistry (Kriest et al., 2012; Kriest and Oschlies, 2015).

The biogeochemical model used here is a model of the marine phosphorus and carbon cycle that resolves phosphate (PO₄), dissolved organic phosphorus (DOP), dissolved inorganic carbon (DIC), total alkalinity, and atmospheric CO₂. Following Kwon et al. (2009), we calculate the production of organic matter using either a nutrient-restoring scheme, where [PO₄] is restored to monthly observations of [PO₄] (Garcia et al., 2014) using a timescale of 30 d (Eq. A3) or using constant export production, where export production is fixed to that of the control run, unless local nutrients fall below zero. These two schemes represent two end-member scenarios, strictly within the context of this model, where organic matter production either depends entirely on macronutrient concentrations and can increase with higher nutrient fluxes (restoring) or is limited by other factors such as light or micronutrients (constant export). The remineralization of particulate organic matter is parameterized using the Martin curve (Eq. 1). To further facilitate a large number of experiments for the sensitivity analysis, we use the model to define a statistical relationship between preformed PO₄ (P_pre) and atmospheric CO₂ and run the model with a phosphorus cycle only (see Sect. A2.3).

2.2 Experiment design

Defining regions

We define a set of oceanic regions to approximately encapsulate the large-scale variability in biogeochemistry and patterns of remineralization depths observed in sediment trap studies. We define regions by lines of latitude and basins, similar to the approach used by air–sea flux inversion studies, e.g. Gloor et al. (2001) and Mikaloff Fletcher et al. (2006). All 15 regions are defined based on a partitioning by Gloor et al. (2001), with some minor changes (Fig. 2a). The assigned regions broadly correspond with major features in observed surface [PO₄], such as higher concentrations in upwelling regions and lower concentrations in the nutrient-depleted gyres (Fig. 2b). This suggests the regions should be a reasonable analogue for an alternative approach that captures key spatial variability in ecology and biogeochemistry by defining regions using vertical mixing, mixed layer depths, sea ice, and sea surface temperature (Longhurst, 1998; Sarmiento et al., 2004; Henson et al., 2010; Fay and McKinley, 2014). The regions are also comparable to the ocean biomes defined in previous biological pump studies (e.g. Weber et al., 2016; Pasquier and Holzer, 2016).

Figure 2 (a) Location and names of the 15 regions defined on the model grid based on Gloor et al. (2001). Boundaries are at 58^∘ S, 36^∘ S, 13^∘ S, 13^∘ N, and 36^∘ N. The equatorial Pacific is split at 98.75^∘ E, following Mikaloff Fletcher et al. (2006). Each region can be assigned a value of b that is independent of other regions. (b) Location of regions superimposed on the annual mean surface [PO₄] from World Ocean Atlas 2013 (Garcia et al., 2014) regridded to the model grid.

2.3 Experiments

We perform a set of experiments to explore the sensitivity of atmospheric CO₂ to regional variability in b with the aim of quantitatively ranking the sensitivity of atmospheric CO₂ to remineralization depth changes in each region (e.g. Pianosi et al., 2016):

1. Control run. A pre-industrial control run is set up with export production diagnosed by restoring it to observed surface [PO₄] values, a globally uniform Martin exponent of 1.0 and initialized with globally uniform tracer concentrations. Atmospheric CO₂ is restored to 278 ppm. A globally uniform Martin exponent gives the lowest volume-weighted root-mean-square misfit compared to annual mean World Ocean Atlas 2013 [PO₄] observations (Garcia et al., 2014), as found in other studies using the same MITgcm transport matrices (Kriest et al., 2012). The control run is spun-up from uniform initial conditions for 5000 years.

2. Global sensitivity. The Martin curve is varied globally, i.e. all regions are assigned the same exponent, between 0.4 and 1.6 in 0.1 increments, based on the range of spatial variability observed in the modern ocean (Henson et al., 2012; Marsay et al., 2015; Guidi et al., 2015). Each experiment is run for 3000 years continuing from the control run using the nutrient-restoring scheme to predict export production (Eq. A3) and freely evolving atmospheric CO₂.

3. Regional sensitivity. Latin hypercube sampling, a stratified random procedure that provides an efficient way of sampling high dimensional parameter space (McKay et al., 1979), is used to vary the Martin curves in every region simultaneously. Values of b are sampled from a uniform distribution ranging from 0.4 to 1.6 using the “lhsdesign” function in MATLAB with “maximin” sampling (an additional constraint that helps reduce clustering of samples, by maximizing the minimum distance between points, in order to give a well-spread distribution of points across the parameter space). The range of b used centres around b=1.0, as also used for the control run. We generate a Latin hypercube ensemble with 200 experiments, balancing the need for higher sampling resolution of the parameter space and total computational time. Each experiment is run for 3000 years continuing on from the control run with a phosphorus cycle only. Changes in atmospheric CO₂ are inferred from changes in preformed PO₄.

Both the global and regional sensitivity experiments are repeated with constant-export production that is diagnosed from the control run. All output is diagnosed from the last full simulation year.

2.4 Sensitivity analysis

We use multiple linear regression analysis to derive the sensitivity of atmospheric CO₂ to changes in b in each region (k) where the fitted coefficients (β_k) give a quantitative measure of the sensitivity (e.g. Pianosi et al., 2016):

$$\begin{array}{}\text{(2)}& {\mathrm{CO}}_{\mathrm{2}}={\mathit{\beta }}_{\mathrm{0}}+\sum _k{\mathit{\beta }}_k{b}_k.\end{array}$$

3 Results

3.1 Sensitivity of CO₂ to regional variability in remineralization depths

To quantify the sensitivity of CO₂ to regional changes in b, we fit linear regression models (Eq. 2) to the results of the Latin hypercube ensembles. The resulting regression models explain a large proportion of the variability between CO₂ and b (R²=0.88 and 0.90 for the constant-export and nutrient-restoring ensembles, respectively). Residuals of the regression models showed no significant bias versus the regression output (Fig. S5 in Supplement), suggesting that a linear model was appropriate. Although the relationship between CO₂ and a globally uniform remineralization depth is non-linear (e.g. Fig. 6), the relationship is near linear around the observed global mean in the centre of the range of b tested (see also, Kwon et al., 2009). Overall, the absence of a strongly non-linear relationship suggests the use of a linear regression model is appropriate (Pianosi et al., 2016).

When b is varied as a globally uniform parameter from 0.4 to 1.6, atmospheric CO₂ varies from 197 to 347 ppm (range of 150 ppm) for the constant-export scheme and from 257 to 288 ppm (range of 31 ppm) for the nutrient-restoring scheme, consistent with previous model experiments (Kwon et al., 2009). Figure 3 shows how the sensitivity of CO₂ to changes in b varies as a function of region. CO₂ is most sensitive to changes in b occurring in the subantarctic regions, with CO₂ being most sensitive to changes in the Indian Ocean sector of the subantarctic (Fig. 3, Table 1). The Southern Ocean and subtropical gyres, with the exception of the gyre in the North Pacific, are consistently the regions where b has the smallest impact on CO₂. Other regions, including the equatorial Indian Ocean, equatorial Pacific, and North Pacific have an intermediate sensitivity. As a region, the subantarctic is responsible for ∼30 % of the CO₂ sensitivity, comparable to the Pacific basin-scale sensitivity (Table 1).

Table 1 Key metrics and sensitivity estimates for each region and basin. Representative metrics including area (percentage of global area), region-integrated mean annual POP export (percentage of global POP export), and mean preformed [PO₄] (percentage of global mean) are taken from the control run. CO₂ sensitivity estimates for both the nutrient-restoring and constant-export ensembles are given as ΔCO₂∕Δb (β_k from Eq. 2) with 95 % confidence intervals and as a percentage (${\mathit{\beta }}_k/{\sum }_k{\mathit{\beta }}_k$).

^a For comparison with the Southern Ocean estimate, defined as >40^∘ S in Kwon et al. (2009). ^b SubPol-PAC, SubPol-ATL and SubPol-IND.

Figure 3 Regional sensitivity (ΔCO₂∕Δb) of atmospheric CO₂ (ppm) to changes in Martin curves for the constant-export scheme (a, b) and nutrient-restoring scheme (c, d). The sensitivity value reflects the increase in CO₂ (preformed PO₄) for an increase in b (shallower remineralization). Error bars in panels (a) and (c) are 95 % confidence intervals for the linear regression coefficients. Atmospheric CO₂ is inferred from modelled preformed PO₄ using empirical relationships in Fig. A1 in Appendix A.

As with the globally uniform changes in b, the magnitudes of regional sensitivities are smaller when run using nutrient-restoring uptake as opposed to a constant-export scheme because export production is able to convert any increase in surface nutrient and carbon fluxes back into organic matter, limiting any change in CO₂ fluxes. However, the relative sensitivity ranked across regions remains similar, as shown by expressing β_k as a percentage of ∑_kβ_k (Table 1). Therefore, the regional patterns in Fig. 3 are not sensitive to assumptions about the response of nutrient uptake to the redistribution of nutrients. This suggests that whilst the absolute magnitude of CO₂ sensitivity to changes in b is related to global export production, it is not driven by local changes in export production specific to any region(s).

Kwon et al. (2009) demonstrated that the sensitivity of CO₂ to basin-scale changes in b correlated with the magnitude of export production in each basin. Similarly, we find a general positive correlation between sensitivity and regional export production (r=0.79, p<0.01 for constant export; r=0.47, p=0.07 for restoring uptake), as measured by the mean annual average export production across the 200 ensemble runs (Fig. 4). The correlation is much weaker with nutrient-restoring uptake compared to the constant-export production. Intuitively, regions with lower export production, i.e. that contribute less to the inventory of regenerated PO₄, have a smaller impact on the balance between preformed and regenerated nutrients and therefore on atmospheric CO₂. Whilst CO₂ is generally more sensitive to remineralization depths in regions with higher export production, sensitivity varies across regions with similar export production. For example, the sensitivity for the temperate North Pacific (NTemp-PAC; Fig. 4a) ($\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}/\mathrm{\Delta }b=\mathrm{13.72}$, export production = 1.4±0.11 Tmol P yr⁻¹) is approximately double that of the subpolar region of the South Pacific (STemp-PAC; Fig. 4a) ($\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}/\mathrm{\Delta }b=\mathrm{5.72}$, export production = 1.4±0.17 Tmol P yr⁻¹). There are no apparent relationships between the variability of export production across the ensemble in each region, as shown by the horizontal error bars, and sensitivity (Fig. 4). This further supports the finding that the response of export production to changes in nutrient distributions are not an important factor in the sensitivity of CO₂ to regional changes in b.

Figure 4 Relationship between regional CO₂ sensitivity and annual export of PO₄ in each region using (a) the constant-export scheme and (b) the nutrient-restoring scheme. Annual export is shown as the mean of the 200 ensemble experiments with ±1 standard deviation error bars.

The variability in sensitivity not explained by the magnitude of POC export is likely a function of how changing remineralization depths interact with ocean circulation. To quantify this effect we calculate the mean preformed PO₄ in the ocean interior derived from each region individually ($\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}$; see Appendix A2.3) and repeat the sensitivity regression analysis. In order to compare regression coefficients from different regions with significantly different magnitudes of $\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}$ we first normalize the concentrations ($\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}$) to the range of variability across the 200 experiments:

$$\begin{array}{}\text{(3)}& \widehat{\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}-min\left(\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}\right)}{max\left(\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}\right)-min\left(\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}\right)}}.\end{array}$$

The new regression analysis (Eq. 4) now predicts the contribution of changing b in all regions to the change in $\widehat{\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}}$ derived from a single region, rather than globally (Fig. 5). The regression analysis is repeated for each region.

$$\begin{array}{}\text{(4)}& \widehat{\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}}={\mathit{\beta }}_{\mathrm{0}}^r+\sum _k{\mathit{\beta }}_k^{\mathrm{region}}{b}_k\end{array}$$

Figure 5 shows the relative sensitivity of $\widehat{\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}}$ to changes in b. R² ranges from 0.82 to 0.97, suggesting that, overall, the linear regression models are appropriate. The regression coefficients specific to a single region are collected from across the 15 regression results in each column of Fig. 5. By definition, each row of Fig. 5 shows the effect of changing b in the corresponding region on $\widehat{\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}}$ from all other regions.

Figure 5 Sensitivity of steady-state normalized mean preformed [PO₄] exported from each region ($\widehat{\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}}$). $\widehat{\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}}$ from each region is expressed as a function of b using linear regression. $\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}$ is normalized to the range of values within each region in the ensemble to account for large differences in preformed $\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}$ between regions. The regression coefficients are arranged such that each row shows the impact of changing b in that region on $\widehat{\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}}$ across other regions. Panels (a) and (b) show the results for the constant-export and nutrient-restoring schemes, respectively.

The sensitivity analysis for each region on an individual basis shows that changes in b in the subantarctic regions have large impacts on $\widehat{\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}}$ across regions globally (Fig. 5). In particular, these regions have a particular effect on the $\widehat{\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}}$ export in the Southern Ocean and in the Atlantic as a basin, comparable in magnitude to the local changes in $\widehat{\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}}$. Changes in b in the equatorial upwelling regions of the Pacific and Indian Oceans also have a large global effect but with a more pronounced local effect. These features are more pronounced with nutrient-restoring uptake (Fig. 5b). The Southern Ocean and North Atlantic regions are those with the highest $\widehat{\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}}$ export across the ensemble (Table 1), consistent with previous findings about the global importance of these regions for preformed PO₄ (DeVries et al., 2012; Pasquier and Holzer, 2016). This suggests the larger sensitivity of CO₂ to changes in b in the subantarctic regions is due to the way in which the ocean circulation connects these regions to the Southern Ocean and North Atlantic. In contrast, changing b in the Southern Ocean and North Atlantic has a relatively minimal effect on $\widehat{\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}^{\mathrm{region}}\right]}}$ (and by inference CO₂) both locally and globally.

3.2 Regional versus global sensitivity

Lastly, we explore whether the spatial patterns in sensitivity (Fig. 3) are significant on a global scale. Global average values of b are calculated for each of the 200 Latin hypercube samples using an area-weighted mean and compared against experiments where b is perturbed uniformly across regions (Fig. 6). The amount of organic matter reaching the deep ocean is a non-linear function of b, following from the fact that the Martin curve represents the scenario of a fixed remineralization rate and an increasing sinking rate (Kriest and Oschlies, 2008; Cael and Bisson, 2018; see Supplement). Therefore, larger values of b, i.e. shallower remineralization, have disproportionally more weight when calculating a global arithmetic mean of spatially variable b values (Fig. S6). To account for this, we find the equivalent e folding depths for each Latin hypercube sample, which form a skewed distribution due to higher occurrence of shallower remineralization, calculate the area-weighted geometric mean e folding depth for each sample, and re-arrange again for b (see Supplement for details).

The relationship between CO₂ and the globally averaged b values from the sensitivity experiments closely matches the relationship between CO₂ and globally uniform b for the both constant-export and nutrient-restoring schemes (Fig. 6). Note that b in each region is varied within the full parameter range but because Latin hypercube sampling varies all parameters across their parameter range simultaneously the global mean does not reach the highest and lowest global b values. The average regionally varying b values vary within $\sim \pm \mathrm{15}$ ppm of the globally uniform experiments with constant-export experiments and $\sim \pm \mathrm{5}$ ppm for the nutrient-restoring experiments, comparable to the change in CO₂ for a globally uniform change in b of ∼0.2.

Figure 6 Comparison of CO₂ sensitivity when b is varied as a globally uniform parameter (solid lines) and when b is varied regionally in the Latin hypercube samples and calculated as an area-weighted geometric mean of e folding depths converted back to b to correct for non-linearities in the Martin curve. Runs using the constant-export scheme are shown in black and nutrient-restoring runs are shown in grey.

4 Discussion

Sediment trap observations reveal significant spatial variability in remineralization depths. Here we have quantified the sensitivity of atmospheric CO₂ to regional changes in remineralization depths and show that CO₂ is most sensitive to changes in the subantarctic regions. Much of the observed spatial variability varies across latitudes (Henson et al., 2012; Guidi et al., 2015; Marsay et al., 2015; Weber et al., 2016). Additionally, the mechanisms potentially driving these patterns are also likely to vary on a latitudinal basis, with changes in related environmental properties in response to anthropogenic CO₂ emissions affecting the high latitudes in particular: temperature changes (Kirtman et al., 2013) affecting temperature-dependent remineralization rates, a reduction in carbonate saturation state with ocean acidification (Orr et al., 2005) affecting ballasting and changes in plankton community composition, and cell size (Lefort et al., 2015) affecting aggregation dynamics and particle sinking velocities. Additionally, this is a consideration for changes in remineralization depths occurring in past climates (Meyer et al., 2016). This suggests that the spatial patterns in CO₂ sensitivity could be significant when considering the impact of remineralization depth changes.

Changes in the air–sea balance of carbon are commonly related to changes in preformed nutrients. Because of the inefficient utilization of upwelled nutrients in the Southern Ocean, this region has been identified as key to setting the efficiency of the biological pump (Ito and Follows, 2005; DeVries et al., 2012). Our results show that this is also key for the sensitivity of CO₂ to regional variability in remineralization depths because of upwelling in subantarctic regions (Fig. 5). This relationship has implications when invoking changes in the efficiency of the biological pump in past climates such as the Last Glacial Maximum (LGM). Processes that increase the utilization of nutrients in the Southern Ocean, such as iron fertilization, and processes that reduce the delivery of nutrients to the Southern Ocean, such as increased stratification, have been implicated in the drawdown of atmospheric CO₂ during the LGM (Sigman et al., 2010). Any changes in stratification will also impact the sensitivity of CO₂ to any additional changes in remineralization depths, such as from changes in ballasting minerals and/or temperature-dependent remineralization (Chikamoto et al., 2012). In comparison, processes such as iron fertilization will not impact on this sensitivity. Because the spatial patterns of CO₂ sensitivity to regional changes in remineralization are predominantly constrained by ocean circulation pathways, this also suggests that the sensitivity may change with a reorganization of ocean circulation as suggested for the LGM (Sigman et al., 2010).

The Martin curve is a commonly used parameterization of the remineralization of particulate organic matter with depth in marine biogeochemical models and is commonly applied with a globally uniform exponent (b) (Hülse et al., 2017). However, the Martin curve used in this way has potential limitations: it is an empirical and static parameterization that does not represent the mechanisms affecting remineralization and sinking rates and it does not capture spatial variability in remineralization observed in sediment trap data (Henson et al., 2012; Marsay et al., 2015; Guidi et al., 2015). In our sensitivity analysis, we have shown that CO₂ has a similar sensitivity to the global mean change in b, as compared to a globally uniform change in b with an uncertainty of ±5–15 ppm, equivalent to a change in b of ∼0.2 (Fig. 6). Kwon et al. (2009) suggest a decrease of 0.3 from the modern remineralization depth is sufficient to explain the increase in deep ocean nutrient concentrations during the Last Glacial Maximum. For the 21st century, Laufkötter et al. (2017) predict a decrease in POC export at 500 m by 2100 under RCP8.5 in response to temperature- and oxygen-dependent remineralization, equivalent to a decrease in b of ∼0.25. As such, the global mean change in potential future and past changes in remineralization depth may be larger than the uncertainty associated with spatial variability. This has potentially useful implications for modelling the remineralization of particulate organic matter fluxes. Models resolving the various processes that affect remineralization rates and sinking velocities have recently been developed (Jokulsdottir and Archer, 2016; Cram et al., 2018); however, the requirements to model processes such as particle aggregation can be computationally expensive, limiting their application to 1-D models (Jokulsdottir and Archer, 2016; Cram et al., 2018) or to offline models (DeVries et al., 2014). A globally uniform change in b informed by these models could then be used to calculate the impact on atmospheric CO₂ if the change in b is greater than 0.2. However, we note that the modern global mean b is subject to uncertainty associated with under-sampled spatial variability.

Our results are dependent on the use of transport matrices derived from one global circulation model. Whilst this model has been widely applied to study biogeochemistry previously, it is subject to a number of caveats. The ocean model predicts significantly larger outcrops of dense water in the Southern Ocean compared to observations (see Fig. S3 in Supplement, also Duteil et al., 2013), leading to deep-water formation occurring at latitudes around 50^∘ S (Fig. S4 in Supplement). The volumetric fraction of water in the ocean interior derived from the subantarctic is also higher (26 %) compared to data-constrained estimates (18 %, Table S1: Khatiwala et al., 2012). As such, the sensitivity estimates for the subantarctic may be overestimated. This is also consistent with the higher sensitivity compared to the basin-scale analysis of Kwon et al. (2009), who found that the Southern Ocean (>40^∘ S) contributed 22 % of the global CO₂ sensitivity, compared with 36 % in this study (>38^∘ S, see Table 1). However, our results have key similarities, including absolute and relative magnitudes of regional preformed PO₄ export, to other studies using alternative steady-state circulation states (DeVries et al., 2012; Pasquier and Holzer, 2016). As such, our results should be broadly reproducible with other models. A disadvantage to using a steady-state circulation is that we cannot quantify impact of the CO₂–climate feedback on ocean circulation and atmospheric CO₂. Studies exploring the simultaneous effects of warming temperatures on circulation and biology in response to anthropogenic CO₂ emissions show that changes in circulation could be as important as biological changes (Cao and Zhang, 2017; Yamamoto et al., 2018). Quantifying the regional sensitivity with a dynamic ocean is therefore an important focus for future work. Ratios of CaCO₃ to POC vary latitudinally, and could also therefore modify our sensitivity results. Segschneider and Bendtsen (2013) found important feedbacks involving interactions between calcifiers and silicifiers in an marine ecosystem model when exploring temperature-dependent remineralization rates in the 21st century. Future model experiments including a representation of plankton ecosystems would therefore help explore the impact of CaCO₃ export on regional sensitivity patterns.

5 Conclusions

We have presented a sensitivity analysis that quantifies the sensitivity of atmospheric CO₂ to regional variability in particulate organic carbon remineralization depths. CO₂ is most sensitive to changes in remineralization depths occurring in the subantarctic regions, particularly in the Indian Ocean sector. As a whole, the subantarctic regions have a sensitivity similar to that of the Pacific basin despite the smaller area and levels of production. Sensitivity patterns are in part a function of the magnitude of export production in each region and the physical circulation pathways specific to each region. Whilst the overall magnitude of CO₂ sensitivity to regional changes is dependent on the magnitude and response of export production to changes in nutrients, the relative spatial patterns in sensitivity are predominantly constrained by ocean circulation pathways. We also find that the regional variability adds ±5–15 ppm uncertainty to global mean changes in remineralization depths. The regional patterns in sensitivity could be significant if a number of processes that potentially drive changes in remineralization depths, including temperature-dependent remineralization rates and plankton community structure, vary predominantly in the high latitudes. However, this uncertainty is similar to the change in CO₂ for a globally uniform change in b of ∼0.2, meaning that larger changes in b could be reliably approximated by a globally uniform b as commonly used in biogeochemical models.

Appendix A: Model description

The Latin hypercube sampling approach used relies on the ability to run an ensemble of model experiments. To make this approach feasible we use the transport matrix method (Khatiwala et al., 2005; Khatiwala, 2007).

A1 Steady-state ocean circulation model

The matrix used here is the 2.8^∘ global configuration of the MIT model with 15 vertical levels driven by seasonally cycling fluxes of momentum, heat, and freshwater (publicly available from http://kelvin.earth.ox.ac.uk/spk/Research/TMM/TransportMatrixConfigs; last access: 29 July 2019). Seasonally varying ocean circulation is calculated at each time step by linearly interpolating between monthly mean matrices. An advantage of using transport matrices is that the time step can be made longer to reduce computational expense (Khatiwala, 2007). Here we extend the circulation time step to 3.8 d.

A2 Biogeochemical model

The biogeochemical model represents the cycle of phosphorus and carbon in the ocean with four dissolved tracers, PO₄, dissolved organic phosphorus (DOP), dissolved inorganic carbon (DIC), and total alkalinity. The biogeochemical model has the same time step as the ocean circulation model (3.8 d).

A2.1 Phosphorus cycle

PO₄ and DOP are governed by the following equations:

$$\begin{array}{}\text{(A1)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{\mathrm{PO}}_{\mathrm{4}}}{\mathrm{d}t}}=\mathbf{A}{\mathrm{PO}}_{\mathrm{4}}-J_{\mathrm{up}}+J_{\mathrm{POP}}+J_{\mathrm{DOP}},\text{(A2)}& {\displaystyle }& {\displaystyle \frac{\mathrm{dDOP}}{\mathrm{d}t}}=\mathbf{A}\mathrm{DOP}+v\cdot J_{\mathrm{up}}-J_{\mathrm{DOP}},\end{array}$$

where A denotes the transport matrix calculation of ocean transport and J denotes biogeochemical source and sink terms.

The uptake of PO₄ during production of organic matter occurs in the euphotic zone, here defined as the base of the upper two grid boxes (120 m). Following Kwon et al. (2009) we calculate the production of organic matter using either a nutrient-restoring scheme or a constant-export scheme. The nutrient-restoring scheme restores surface concentrations of PO₄ to observed [PO₄] with a restoring timescale (τ=30 d Najjar et al., 2007) and is scaled by the fraction of sea ice present (F_seaice, as monthly average fields from the original global circulation model):

$$\begin{array}{}\text{(A3)}& J_{\mathrm{up}}={\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}max\left(\right({\mathrm{PO}}_{\mathrm{4}}-{\mathrm{PO}}_{\mathrm{4},\mathrm{obs}}),\mathrm{0})(\mathrm{1}-F_{\mathrm{seaice}}).\end{array}$$

Organic matter production in the constant-export scheme is fixed to that of the experiment defined as the control run, unless surface [PO₄] is depleted below zero, in which case J_up is set to zero at that time step. The control run is defined as having the run with the lowest root-mean-square misfit compared to annual mean World Ocean Atlas [PO₄] observations.

A fixed fraction (v=0.66) of the organic matter production integrated across the upper two grid boxes is routed directly to dissolved organic phosphorus (DOP) and remineralized back to PO₄ in a first-order reaction with a decay rate of κ throughout the water column:

$$\begin{array}{}\text{(A4)}& J_{\mathrm{DOP}}=\mathit{\kappa }\mathrm{DOP}.\end{array}$$

The remaining fraction of organic matter production ($\mathrm{1}-v=\mathrm{0.34}$) is integrated across the upper two grid boxes and exported as particulate organic phosphorus (POP) at the base of the of the second grid box in the vertical (120 m). The remineralization of POP is parameterized with the Martin curve (Eq. 1). POP that has reached the sediment is remineralized fully in the lowermost grid box of the water column, maintaining a closed system with respect to [PO₄]. As such, there is no sediment component in this model.

Figure A1 Relationship between preformed PO₄ and atmospheric CO₂ in the biogeochemical model using nutrient-restoring and constant-export schemes when the Martin curve is varied globally between 0.4 and 1.6. A quadratic function ($\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}={\mathit{\beta }}_{\mathrm{1}}\mathrm{\Delta }{{\mathrm{PO}}_{\mathrm{4},\mathrm{pre}}}^{\mathrm{2}}+{\mathit{\beta }}_{\mathrm{2}}\mathrm{\Delta }{\mathrm{PO}}_{\mathrm{4},\mathrm{pre}}+{\mathit{\beta }}_{\mathrm{3}}$) is fitted to the combined data with non-linear least squares is shown with 95 % confidence intervals. The coefficients for the two fits are β₁=66.59, ${\mathit{\beta }}_{\mathrm{2}}=-\mathrm{22.48}$ and β₃=187.09 for constant-export schemes and β₁=42.80, ${\mathit{\beta }}_{\mathrm{2}}=-\mathrm{6.69}$ and β₃=207.95 for nutrient-restoring schemes.

A2.2 Carbon cycle

The uptake of nutrients and remineralization of particulate and dissolved organic phosphorus are related to dissolved inorganic carbon and alkalinity via Redfield ratios of $\mathrm{1}:\mathrm{16}:\mathrm{116}$ ($\mathrm{P}:\mathrm{N}:\mathrm{C}$). Carbonate chemistry parameters are computed from dissolved inorganic carbon and alkalinity using the method described by Follows et al. (2006). The method provides a simplified but accurate solution with computational efficiency. The air–sea gas exchange of CO₂ is calculated as per Orr et al. (2017).

A2.3 Preformed PO₄ and atmospheric CO₂

Changes in atmospheric CO₂ due to changes in the biological pump can be directly related to the inventory or average concentration of preformed PO₄ (P_pre) if total nutrient concentrations are conserved (Ito and Follows, 2005; Marinov et al., 2008). This provides a way of relating changes in our model of the phosphorous cycle to changes in atmospheric CO₂ without simulating a relatively computationally expensive carbon cycle. The distribution of annual mean [P_pre] for each run is calculated by splitting the transport matrices into “interior” matrices (A^I) and “exterior” matrices (B) for both the explicit and implicit matrices (subscripts e and i, respectively; see Khatiwala, 2007). The annual mean surface [PO₄] from the end of a simulation is set as a boundary condition and solved for the interior distribution of P_pre:

$$\begin{array}{}\text{(A5)}& P_{\mathrm{pre}}=({\mathbf{A}}_{\mathrm{i}}^I{\mathbf{A}}_{\mathrm{e}}^I-I{)}^{-\mathrm{1}}(({\mathbf{A}}_{\mathrm{i}}^I{\mathbf{B}}_{\mathrm{e}}+{\mathbf{B}}_{\mathrm{i}}){\mathrm{PO}}_{\mathrm{4}}).\end{array}$$

The global mean concentration of P_pre ($\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}\right]}$) is related to CO₂ using a empirical quadratic function (Eq. A6) fitted to a series of experiments where the Martin curve is varied globally (see below). The function is derived from a non-linear least-squares regression (Fig. A1). The resulting regression fit (details in the caption of Fig. A1) is used to estimate changes in CO₂:

$$\begin{array}{}\text{(A6)}& {\mathrm{CO}}_{\mathrm{2}}=({\mathit{\beta }}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}\right]}}^{\mathrm{2}}+{\mathit{\beta }}_{\mathrm{2}}\stackrel{\mathrm{‾}}{\left[P_{\mathrm{pre}}\right]}+{\mathit{\beta }}_{\mathrm{3}}).\end{array}$$

The approach to predict CO₂ from regional changes in b is evaluated in the Supplement and shows the approach is robust to factors such as the slow timescale for the air–sea gas exchange of CO₂ (e.g. Eggleston and Galbraith, 2018).