External Forcing Explains Recent Decadal Variability of the Ocean Carbon Sink

The ocean has absorbed the equivalent of 39% of industrial‐age fossil carbon emissions, significantly modulating the growth rate of atmospheric CO2 and its associated impacts on climate. Despite the importance of the ocean carbon sink to climate, our understanding of the causes of its interannual‐to‐decadal variability remains limited. This hinders our ability to attribute its past behavior and project its future. A key period of interest is the 1990s, when the ocean carbon sink did not grow as expected. Previous explanations of this behavior have focused on variability internal to the ocean or associated with coupled atmosphere/ocean modes. Here, we use an idealized upper ocean box model to illustrate that two external forcings are sufficient to explain the pattern and magnitude of sink variability since the mid‐1980s. First, the global‐scale reduction in the decadal‐average ocean carbon sink in the 1990s is attributable to the slowed growth rate of atmospheric pCO2. The acceleration of atmospheric pCO2 growth after 2001 drove recovery of the sink. Second, the global sea surface temperature response to the 1991 eruption of Mt Pinatubo explains the timing of the global sink within the 1990s. These results are consistent with previous experiments using ocean hindcast models with variable atmospheric pCO2 and with and without climate variability. The fact that variability in the growth rate of atmospheric pCO2 directly imprints on the ocean sink implies that there will be an immediate reduction in ocean carbon uptake as atmospheric pCO2 responds to cuts in anthropogenic emissions.


Introduction
The ocean has absorbed the equivalent of 39% of fossil carbon emissions since 1750, significantly modulating the growth of atmospheric CO 2 and the associated climate change (Ciais et al., 2013;Friedlingstein et al., 2019;Le Quéré et al., 2018a, 2018bMcKinley et al., 2017). If emissions continue to accelerate, this sink is expected to grow and substantially mitigate atmospheric carbon accumulation for the next several centuries (Randerson et al., 2015). Yet, we lack a detailed understanding of spatial and temporal variability in the sink and its underlying mechanisms. Incomplete understanding of ocean flux variability contributes to significant uncertainty in the global anthropogenic carbon budget for recent decades; this uncertainty is equivalent to~10% of the atmospheric pCO 2 growth rate (Friedlingstein et al., 2019; Le Quéré • The reduced ocean carbon sink in the decade of the 1990s can be explained by the slowed growth rate of atmospheric CO 2 • The global sea surface temperature response to Mt Pinatubo in 1991 explains the intra-decadal pattern of the ocean carbon sink in the 1990s • There will be an immediate reduction in ocean carbon uptake as atmospheric pCO 2 responds to cuts in anthropogenic CO 2 emissions  Eddebbar, Y. A., Gloege, L., & Lovenduski, N. S. (2020). External forcing explains recent decadal variability of the ocean carbon sink. AGU Advances, 1, e2019AV000149. https://doi.org/ 10. 1029/2019AV000149 et al., 2018a, 2018b. This imbalance, in turn, limits the scientific community's ability to inform international efforts to curb fossil fuel emissions (Peters et al., 2017).
Recent studies have concluded that ocean hindcast models, which have long been used to assess the ocean carbon sink, may significantly underestimate its variability (Gruber et al., 2019a;Le Quéré et al., 2018b). These conclusions are based on comparisons to new observation-based gap-filled products that suggest substantially larger interannual to decadal variability (Landschützer et al., 2015(Landschützer et al., , 2016Rödenbeck et al., 2013Rödenbeck et al., , 2015. It is possible, however, that these gap-filled products do not accurately represent the surface ocean carbon cycle given that their raw input data cover only 1.5% of the global ocean in the last 3 decades, and only 3.5% in the most recent years . Variability may be amplified by the significant extrapolation that occurs when global full-coverage maps are produced from these very sparse data (Fay et al., 2018;Rödenbeck et al., 2015).
The ocean carbon sink of the 1990s is of particular interest. During this time, the growth of the ocean sink stalled from its expected growth (DeVries et al., 2019;Fay & McKinley, 2013;Landschützer et al., 2015;Le Quéré et al., 2007;Lovenduski et al., 2008). Using a data-constrained ocean circulation model, it has been suggested that this slow-down was caused by excess outgassing of natural carbon due to an anomalously strong overturning of the ocean's upper 1000 m (DeVries et al., 2017). Changing patterns of wind-driven circulation in the Southern Ocean have also been proposed (Gruber et al., 2019a;Landschützer et al., 2015).
Consensus has yet to be achieved on the mechanisms driving changes in the 1990s sink. It is critical that we accurately quantify and diagnose this variability so that we can better project the future ocean carbon sink and, thus, the degree to which the ocean carbon sink will continue to mitigate global climate change.

Methods
In this study, we compare the ocean carbon sink since 1980 as estimated from six ocean hindcast models, four observationally-based products, and a theoretical upper ocean box model. We supplement our analysis with results from nine ocean hindcast models run in constant climate model (DeVries et al., 2019). Supporting information provides additional methodological details.

Models and Products
Six (6) ocean hindcast models for 1980-2017 are the primary basis for this analysis (Table S1). Ocean hindcast models are gridded, three dimensional representations of the evolution of the ocean state for recent decades. These models have been forced at the surface by winds, heat and buoyancy fluxes from reanalyses of the past atmospheric state. The models we analyze are those used in the recent versions of the Global Carbon Budget (Le Quéré et al., 2018aQuéré et al., , 2018b. Four (4) observationally-based products are also utilized, chosen because they all cover 1985-2016 (Table S2). While models have full coverage in space and time, in situ observations only cover a small fraction of the global surface ocean. These observation-based products utilize gap-filling techniques to estimate values in all periods and areas not directly observed. Interpolation methods fill spatial and temporal gaps by assuming statistical relationships with neighboring or similar areas with available observations.
Ensemble means of the 4 observationally-based products and the 6 hindcast models are calculated.

Flux Analysis
Model fluxes span years 1980-2017, while observationally-based products span the 1985-2016 period (Table S1, S2). An area-weighted annual mean timeseries is calculated for each model and product for regions of interest: i) global, ii) east equatorial Pacific biome (Fay & McKinley, 2014), and iii) global, excluding the east equatorial Pacific biome. We note that individual models and products utilize different methods for the flux calculation, including wind speed products and parameterizations. The standard approach of our field is to use the mean of these estimates as the current best-estimate of the air-sea flux (DeVries et al., 2019; Le Quéré et al., 2018aQuéré et al., , 2018b.

Accounting for Open Ocean Areas Without Observationally-Based Estimates
There are large differences in the spatial extent of the observationally-based products compared to the hindcast models. This difference, when left unaccounted for, results in significant discrepancies between models and products for both global mean flux and global mean pCO 2 . We correct for spatial coverage differences by calculating the mean flux or mean pCO 2 in the area of each model where there is no coverage for a specific data product. By calculating this offset from 6 models for each individual data-product's spatial coverage, we generate a product-specific offset for the annual time series that corrects for the difference in spatial coverage. As expected, observationally-based products with more complete spatial coverage have smaller offsets (Table S2).
The only exception to this data-product correction process is the JENA product (Rödenbeck et al., 2013) because it is produced with full spatial coverage and is primarily used in this analysis on its coarser native grid.

Upper Ocean Box Model
The box model ( Figure S1) solves for the time change of Dissolved Inorganic Carbon (DIC) in single upper ocean box (Equation 1).
The first term on the right of Equation 1 is the overturning circulation (ν) acting on the surface to depth gradient of DIC. V is the volume of the global ocean box, V = A*dz. Our value for the global area removes ice-covered regions. The second term on the right of Equation 1 is the air-sea exchange of CO 2 . The rate of flux is set by a piston velocity (k w ), solubility (S o ) and density (ρ) over the depth of the box (dz = 200 m), multiplied by the ocean to atmosphere difference of pCO 2 (Wanninkhof, 2014). S o and pCO 2 ocean are calculated using full carbon chemistry (Dickson & Millero, 1987;Mehrbach et al., 1973) given inputs of temperature, salinity, alkalinity and DIC. Parameters choices are globally representative values outside the tropics (Table S3). Consistent with current understanding of the drivers of ocean uptake of anthropogenic carbon (Gruber et al., 2019b), the biological pump is assumed constant over time. This leads to our abiotic formulation. Thus, we remove from the DIC deep concentration the amount of carbon that would be vertically supplied, and then removed biologically. We take a global mean DIC deep concentration of 2320 mmol/m 3 and a biological pump component of this of 265 mmol/m 3 (Sarmiento & Gruber, 2006,  NOAA ESRL surface marine boundary layer annual mean observed xCO 2 atmosphere is used to force the model. This is the same xCO 2 atmosphere data used to force the ocean hindcast models and observationally-based products, and in conversion to pCO 2 atmosphere , the water vapor correction is applied.   (Table S4). In A, the mean flux of the observationally-based products is increased by 0.45 PgC/yr (Jacobson et al., 2007) to account for the outgassing of natural carbon supplied by rivers to the ocean. extends to several hundred meters depth in CESM LENS. The box model is time stepped with monthly resolution for 1959-2018. In all cases, the box model is spun up from 1959-1979 using observed pCO 2 atmosphere and the values presented in Table S3.
The mean uptake flux in the box model is most sensitive to the depth of the box and the global overturning rate ( Figure S3, Supporting Information). We use dz = 200 m and set other parameters to result in a mean flux and ΔpCO 2 consistent with the ocean models and observationally-based products.

Results/Discussion
Global air-sea CO 2 flux variability estimated by the ensemble means of ocean hindcast models and of observationally-based products (Gruber et al., 2019a;Landschützer et al., 2016) are highly correlated (r = 0.95) ( Figure 1a, Table S4). The decadal variability of the air-sea CO 2 flux is not driven by a single region, but instead is largely globally coherent  Sarmiento, 1993). This change was due in part to a pause of growth in fossil fuel emissions from 1989 to 1994 when fossil fuel emissions were approximately constant at 6.1 PgC/yr (Friedlingstein et al., 2019;Sarmiento et al., 2010). Increased uptake of carbon by the terrestrial biosphere also contributed significantly to this slow down (Angert et al., 2004;Brovkin et al., 2010;Sarmiento et al., 2010). Though some studies have attributed the increased land carbon sink to temperature and radiation changes caused by Mt. Pinatubo's 1991 eruption, there is not a consensus with respect to the mechanisms on land (Angert et al., 2004;Brovkin et al., 2010;Frölicher et al., 2011;Lucht et al., 2002;Sarmiento et al., 2010).
Deviations in the evolution of pCO 2 ocean from the evolution of pCO 2 atmosphere drive ΔpCO 2 changes in the observationally-based products and hindcast models (Figure 3c). Since global mean ΔpCO 2 is only~5 μatm (Figure 3c), anomalies of a few μatm in either pCO 2 atmosphere or pCO 2 ocean can significantly impact the air-sea CO 2 flux. In 1992, growth of pCO 2 ocean abruptly slowed in both the hindcast models and the observationally-based products ( Figure 3b) and ΔpCO 2 becomes more negative (Figure 3c). From 1992 through 2001, pCO 2 ocean increases more rapidly than pCO 2 atmosphere , driving a positive change in ΔpCO 2 over this period (Figure 3c). For 2002-2011 and beyond, pCO 2 ocean grows more slowly than the strongly accelerating pCO 2 atmosphere , leading to increasingly negative ΔpCO 2 (Figure 3c).
Given that global mean changes in the ocean sink ( Figure 1) are found to occur at most latitudes (Figure 2), we hypothesize an important role for external forcing. To explore these drivers, we apply the upper ocean box model ( Figure S1, Equation 1). The key processes captured by this model are (1) ventilation to the surface of deep waters with low anthropogenic carbon content, and (2) air-sea gas exchange. First, we ask: Are  Table S4). This correspondence emphasizes the critical role of variability in the growth of pCO 2 atmosphere to variability in the ocean sink. It also serves as evidence that the box model is realistically estimating the timing and magnitude of this response.
Though it is the slowed growth of pCO 2 atmosphere in the early 1990s that causes the mean 1990s sink to be only 0.1 PgC/yr larger than the sink of the 1980s (Table S5), the pattern of the sink change in the 1990s is clearly inconsistent with ocean hindcast models with variable climate (compare bold green to red dash in Figure 1b), and thus is also inconsistent with the observationally-based products (Figure 1a). An additional mechanism is required.
Volcanic eruptions of El Chichon in 1982 and Mt. Pinatubo in 1991 injected large quantities of sulfate aerosols into the stratosphere and dramatically altered global air and sea surface temperatures (Church et al., 2005). The forced response to these eruptions was a substantial oceanic uptake of carbon and oxygen for the following 2-3 years (Eddebbar et al., 2019). Modern earth system models indicate that a significant negative anomaly in global sea surface temperatures (SST) was driven by the eruptions (Eddebbar et al., 2019). For the diagnostic box model, we apply the same magnitude of forced global SST cooling estimated by these models, 0.1°C in 1982 and 0.2°C in 1991 ( Figure S2) to evaluate the impact on ocean carbon sink variability.
With this volcano-forced SST variability applied to the box model, strong coolings in 1982 and in 1991 drive a rapid drop of pCO 2 ocean (Figure 3b, solid red) and a strong uptake anomaly (Figure 1, red bold). The reduced flux that would have occurred in the early 1990s due to the slowing of pCO 2 atmosphere (Figure 1b, red dash) was overwhelmed by the rapid cooling due to Mt. Pinatubo and thus, a strong uptake pulse occurred (Figure 1b, red bold). The re-warming and excess DIC in the surface ocean that follow Mt. Pinatubo elevates pCO 2 ocean relative to pCO 2 atmosphere through 2001, leading to ΔpCO 2 becoming less negative over this period Observationally-based products mean (blue), hindcast model mean (green) and upper ocean diagnostic box model (red). The box model is forced with only pCO 2 atmosphere (dashed) and with both pCO 2 atmosphere and volcano-associated SST change (solid). Hindcast models that did not have the water vapor correction applied to their atmospheric pCO 2 when the simulation was performed are adjusted to account for that difference (see Supporting Information). Figure S6 expands on these results by including additional box model forcing scenarios.
( Figure 3c). Thus, the sink stagnates from the early to late 1990s ( Figure 1). The net effect of both forcings is that the reduced sink of the early 1990s caused by the slowed pCO 2 atmosphere growth rate is shifted to the late 1990s by the rapid cooling and slow re-warming caused by Mt. Pinatubo. In summary, the climate variability mechanism that led to a neutral 1990s intra-decadal trend of the ocean carbon sink (DeVries et al., 2019) contains a major contribution from the ocean's response and recovery from Pinatubo cooling, i.e. a response to external forcing from volcanos.
CO 2 fluxes from the diagnostic box model forced with both observed pCO 2 atmosphere and the volcanoes' impacts on sea surface temperature are highly correlated to the ensemble mean of the observationally-based products for their overlapping periods (r = 0.89, 1985-2016) and hindcast models (r = 0.92, 1980-2017) (Figure 1a, Table S4). The simplicity of these global mechanisms and the strong correspondence of the resulting decadal variability to the products and the models supports the conclusion that global air-sea CO 2 flux variability since 1980 has been significantly driven by external forcing from (1) the changing pCO 2 atmosphere growth rate and (2) in the 1990s, the surface ocean temperature effects of Mt.
Since ocean carbon uptake is enhanced with the eruption of large volcanos, the effect on pCO 2 atmosphere would ideally be modeled interactively. Unfortunately, land carbon sink uncertainties preclude this. The sea surface temperature effects of Pinatubo led to an increased ocean sink of approximately 0.5 PgC/yr (Figure 1b), but estimates of the land sink anomaly at this time are much larger and uncertain, ranging 1-2 PgC/yr (Angert et al., 2004;Sarmiento, 1993;Sarmiento et al., 2010). In fact, the post-Pinatubo period is one of maximum uncertainty in the post-1960 Global Carbon Budget (Friedlingstein et al., 2019;Peters et al., 2017). Soon after the eruption of Mt. Pinatubo, Sarmiento (1993) noted the coincident slowdown in growth of pCO 2 atmosphere and reported that 13 C records at that time suggested a terrestrial driver, while O 2 /N 2 records suggested an oceanic driver. A modern reconsideration of 13 C and O 2 /N 2 records may lead to better understanding of this partitioning.
This analysis illustrates that externally forced variability played an important role in recent decadal variability of the ocean carbon sink. However, the total climate variability in any variable is the sum of forced variability caused by drivers external to the system, and internal variability due to system dynamics (Deser et al., 2012). We have recently reviewed the many previous studies on mechanisms of ocean carbon sink variability (McKinley et al., 2017). Processes discussed have included the variable upper ocean circulation, wind and circulation patterns in the Southern Ocean, and modes of coupled atmosphere/ocean variability in both hemispheres (DeVries et al., 2017(DeVries et al., , 2019Gruber et al., 2019a;Landschützer et al., 2015Landschützer et al., , 2019Lovenduski et al., 2007). These analyses have focused on variability internal to the ocean or associated with coupled atmosphere/ocean modes (McKinley et al., 2017), but have not been able to comprehensively explain the global-scale decadal variability. Here, we illustrate that the observed changes can, to first order, be attributed to two forcings external to the ocean.
Previous studies have also typically focused on a single model or a single observationally-based product. However, for the best estimate of the real ocean's flux variability, it is common practice to use the ensemble mean of ocean models and/or of observationally-based products (DeVries et al., 2019;Friedlingstein  Le Quéré et al., 2018a, 2018b, which is our approach. Only the internal variability that is represented in most ensemble members will be preserved in the ensemble average. Averaging damps internal variability of the individual members and thus amplifies the common forced component (Deser et al., 2012;McKinley et al., 2016McKinley et al., , 2017. Our results illustrate that the current best estimate of the real ocean's flux variability, based on this ensemble average, can be explained largely with forced mechanisms. However, because we do not have enough information to determine which member of the ensemble best approximates the ocean's true internal variability, this current best estimate potentially underestimates the full impact of internal variability in the carbon sink of the real ocean. What is the range of magnitude of internal variability that may be occurring in addition to the forced variability that we identify? Individual observationally-based products have a range of detrended flux variability from 0.14-0.30 PgC/yr, while the ensemble mean variability is 0.19 PgC/yr (1σ for 1985PgC/yr (1σ for -2016. For the hindcast models, the range is 0.10-0.20 PgC/yr and the ensemble mean variability is 0.11 PgC/yr for the same years. Our estimate of the amplitude of externally-forced variability from the box model is 0.14 PgC/yr ( Figure 1). By this measure, externally forced variability as estimated by the box model is approximately equal to the total amplitude common to the hindcast models, but is only about 70% of the variability common to the products. On top of this, the individual models and products suggest a wide range of additional internal variability. In future studies, separation of the forced component of ocean carbon sink variability driven by changing pCO 2 atmosphere and volcanos from the total variability in individual models and products should help to clarify the patterns, magnitudes, and physical and biogeochemical mechanisms of internal variability in the real ocean. For diagnostic (Friedlingstein et al., 2019;Le Quéré et al., 2018a, 2018b and predictive purposes (Randerson et al., 2015) it is critical to also determine which model and observationally-based estimates best represent both the internal and forced variability of the real ocean.
Though our box model is sufficient to represent the global-mean behavior of the externally-forced ocean carbon sink in recent decades, other mechanisms may increase in importance in the future. As climate changes have increased impact on ocean physics and biogeochemistry, feedbacks on the carbon sink of increasing magnitude can be expected. Future reduction in the overturning circulation, or increased re-emergence of waters already carrying a high anthropogenic carbon load would reduce the sink. A weaker biological pump would also damp net ocean carbon uptake (Kwon et al., 2009). The reduced buffer capacity of the surface ocean should grow in importance over time, particularly under high emission scenarios (Fassbender et al., 2017). As mitigation of CO 2 emissions occurs, the growth rate of pCO 2 atmosphere will slow. With this reduced external forcing, the imprint of internal variability on the sink should become more evident. Improved understanding of both internal and external mechanisms is essential to continue to accurately diagnose the evolving ocean carbon sink, and to improve model-based predictions.

Conclusions
We have shown that externally forced variability is sufficient to explain a significant portion of current model and observationally-based best-estimates of the recent decadal variability of the global ocean carbon sink (Figure 1a). The reduced ocean carbon sink in the decade of the 1990s was driven by a slowed growth rate of pCO 2 atmosphere . The intra-decadal timing of the slowed growth rate in the 1990s was due to the surface ocean temperature response to the Mt. Pinatubo eruption in 1991. Volcano-driven cooling first led to an anomalously large sink, and then as surface ocean temperatures recovered, pCO 2 ocean was elevated causing the sink to slow. In the box model, only this SST response is needed to replicate the behavior of the observationally-based products and the ocean hindcast models (Figure 1), but it would be of great value to perform a deeper analysis of the upper ocean response to Mt. Pinatubo with future studies. From 2001 on, the recovery of the global ocean carbon sink is attributable to the enhanced growth rate of pCO 2 atmosphere (Figure 4).
Implications for the future ocean carbon sink are several. First, we note the relative importance of external vs internal drivers of ocean sink change can be expected to change, and thus both must be understood. Regarding external forcing, future large volcanic eruptions cannot be predicted, and it is difficult to predict the detailed future of pCO 2 atmosphere . Thus, these are now identified as additional sources of uncertainty in decadal predictions and long-term projections ( timescales on which additional human interventions in the climate system, such as solar radiation management or nuclear conflict, would mimic these externally forced changes and modify the ocean carbon sink should be considered (Lauvset et al., 2017;Lovenduski et al., 2020). Finally, since the changing growth rate of pCO 2 atmosphere is the primary driver of recent variability in the ocean carbon sink, the ocean sink should be expected to slow as reductions in the pCO 2 atmosphere growth rate occur in response to climate change mitigation efforts (Peters et al., 2017). It is important that this critical feedback on the atmospheric CO 2 content be accurately estimated and accounted for in policy making.

Data Availability Statement
Ocean hindcast models with real climate are available from https://www. . Funding from many countries and agencies has supported the collection of surface ocean pCO 2 data, for development of ocean hindcast models and observationally-based products, and for international coordination. G.A. M., A.R.F and L.G. were supported by NASA NNX17AK19G and by Columbia University. G.A.M., A.R.F. and N.S.L. were supported by National Science Foundation OCE-1948664 and OCE-1558225. N.S.L. was also supported by NSF OCE-1752724, NSF PLR-13009540, and the Open Philanthropy Project. This work would not be possible without the efforts of many scientists who have collected surface ocean pCO 2 data and contributed it to the SOCAT database, and to the developers of the observationally-based products based on these data. We thank also the scientists who have contributed their ocean hindcast model results to the Global Carbon Project. Leadership from the International Ocean Carbon Coordinating Project (IOCCP) and the Global Carbon Project (GCP) has been essential to the success of these efforts. This is Lamont Doherty Earth Observatory contribution #8400.