Seasonal to decadal spatiotemporal variations of the global ocean carbon sink.

The global ocean has absorbed approximately 30% of anthropogenic CO2  since the beginning of the industrial revolution. However, the spatiotemporal evolution of this important global carbon sink varies substantially on all timescales and has not yet been well evaluated. Here, based on a reconstructed observation-based product of surface ocean pCO2 and air-sea CO2  flux (the MPI-SOMFFN method), we investigated seasonal to decadal spatiotemporal variations of the ocean CO2  sink during the past three decades using an adaptive data analysis method. Two predominant variations are modulated annual cycles and decadal fluctuations, which account for approximately 46% and 25% of all extracted components, respectively. Although the whole summer to non-summer seasonal difference pattern is determined by the Southern Ocean, the non-summer CO2  sink at mid-latitudes in both hemispheres shows an increasing trend (a total increase of approximately 1.0 PgC during the period 1982-2019), while it is relatively stable in summer. On decadal timescales for the global ocean carbon sink, unlike the weakening decade (1990-1999) and the reinvigoration decade (2000-2009) in which the Southern Ocean plays the dominant role, the reinforcement decade (2010-2019) is mainly the result from the weakening source effect in the equatorial Pacific Ocean. Our results suggest that except for the Southern Ocean's role in the global ocean carbon sink, the strengthening non-summer's sink at mid-latitudes in both hemispheres and the decadal or longer timescales of equatorial Pacific Ocean dynamics should be fully considered in understanding the oceanic carbon cycle on a global scale.

INTRODUCTION

The global ocean has absorbed approximately one‐third of the carbon emitted by fossil and land‐use change since the beginning of the industrial revolution, significantly regulating the growth patterns of the atmospheric CO2 and the associated climate change (Friedlingstein et al., 2020; Landschützer et al., 2014; Sabine et al., 2004). The uptake of CO2 by the ocean and the outgassing of CO2 from the ocean are spatiotemporally uneven. The strength of the global ocean carbon sink in each ocean region is determined by chemical, biological, and physical processes. There are considerable large discrepancies in CO2 uptake among oceans regarding their seasonal to decadal variations, which are linked to the internal variability of climate systems (Ilyina, 2016; McKinley et al., 2020).

The Southern Ocean, south of 35°S, accounts for approximately 40% of global oceanic uptake of anthropogenic CO2 (Landschützer et al., 2015). The carbon sink in the Southern Ocean is quite sensitive to physical climate variability, varying substantially on all timescales but showing more distinct seasonal and decadal variations. Seasonal fluctuations are induced mainly by the strong compensating effects of ocean biology, mixing, and warming/cooling (Gruber et al., 2019). Decadal variabilities show a period of weakening throughout much of the 1990s and then a period of rebound after 2000 (Landschützer et al., 2015). Although mechanisms are not fully captured, the weakening period during the 1990s should be caused by a southward shift of the westerlies that enhanced the upwelling and outgassing of CO2. The rebound decade in the 2000s was probably driven by cooling in the Pacific sector, enhanced stratification in the Atlantic and Indian Ocean sectors, and a weakening of meridional overturning circulation (Gruber et al., 2019; Landschützer et al., 2016; Pérez et al., 2013).

The equatorial ocean plays an important role in the global carbon cycle, as it is the ocean's largest natural source of CO2 to the atmosphere, with estimations of its annual contribution ranging between 0.6 and 1.5 Pg (Feely et al., 2006). Equatorial waters are comparatively active in exchanging CO2 with the atmosphere in the context of long‐term changes and exhibit a large amount of spatial and temporal variability due to interannual and decadal variabilities, especially in the central and eastern Pacific Ocean (Feely et al., 2006; Sutton et al., 2014). The large interannual variability in air–sea CO2 flux within the tropical Pacific is thought to play a dominant role in the interannual variability in the global oceanic CO2 uptake (Takahashi et al., 2009; Wanninkhof et al., 2013). This variability has been proven to be associated with El Niño‐Southern Oscillation (ENSO) (Yasunaka et al., 2019). During El Niño events, equatorial CO2 outgassing in the tropics weakens due to the weakening upwelling and retreat of the cold tongue (Ishii et al., 2014; Strutton et al., 2008). On decadal timescales, decadal changes in the fugacity of CO2 have been proven to be correlated with the sea surface temperature, wind speed, and regime shifts of the Pacific Decadal Oscillation (PDO) (Feely et al., 2006).

Even though these main mechanisms are known, there are considerable uncertainties regarding seasonal to decadal variations of the global ocean carbon sink (Ilyina, 2016) due to limited observational data. In addition, ocean models fail to simulate low‐frequency variabilities, such as decadal variations (Ishii et al., 2014; Lenton et al., 2013; Wanninkhof et al., 2013). Although atmospheric CO2 inversions show lagerer variability than ocean inversions, the general patterns of the timing and magnitude are rather incongruent with those of ocean models (Peylin et al., 2013).

To quantify long‐term variations of the global ocean carbon sink, the surface ocean CO2 measurement community created the Surface Ocean CO2 Atlas (SOCAT), with the latest version of SOCAT (version 2021) including 30.6 million observations of global oceans and coastal seas ranging from 1957 to 2020 (Bakker et al., 2016). Although data coverage has increased substantially in recent years, there are only a few regions (about 1.4%) in the global ocean for which sea surface pCO2 has been monitored for decades at the same location (Bakker et al., 2016; Bates et al., 2015; Landschützer et al., 2016). Thus, a series of data interpolation and extrapolation‐based methods, together with advances in remote sensing and machine learning techniques (Land et al., 2019), have been proposed to fill data gaps in the majority of the global ocean (Jones et al.,2015; Landschützer et al., 2014; Rödenbeck et al., 2014; Sasse et al., 2013; Zeng et al., 2014). To date, the most established approaches are the SOM‐FFN method of Landschützer et al. (2013), the CSIR‐ML6 method of Gregor et al. (2019), the CMEMS‐FFNNv2 method of Denvil‐Sommer et al. (2019), the Jena‐MLS method of Rödenbeck et al. (2014), and the NIES‐FNN method of Zeng et al. (2014). Although each of these methods has its own strength and weakness, the majority of methods fall into an appropriate root mean squared error of approximately 10 μatm when compared with SOCAT (Gregor et al., 2019). A series of previous works using the so‐called SOM‐FFN method have pointed out seasonal, interannual, and decadal variations in the global ocean carbon sink; however, substantial uncertainties remain (Landschützer et al., 2016). Particularly, one needs attention is that the level of variabilities at different timescales in the global ocean carbon sink, and mechanisms that drive these variations in the global ocean carbon sink.

Here, we further investigate the spatiotemporal evolutions from seasonal to decadal timescales using an adaptive data analysis method based on the reconstructed product of the surface ocean pCO2 via the SOM‐FFN method (Landschützer et al., 2013). We extracted the seasonal to decadal variability from these air–sea flux maps and then examined the contributions from high‐ to low‐frequency components. Emphasis is placed on the spatiotemporal variability of seasonal and decadal fluctuations, as well as the main drivers of these fluctuations.

DATA AND METHODOLOGY

The data used in this study are an updated observation‐based mapped estimate of the air–sea CO2 flux taken from SOCAT (Bakker et al., 2016; Landschützer et al., 2014). Landschützer et al. (2013, 2014) developed a two‐step neural network mapping approach to reconstruct basin‐wide monthly maps of the sea surface partial pressure of CO2 (pCO2) flux at a resolution of 1° × 1°. The SOM‐FFN approach successfully overcame the low spatiotemporal density of surface carbon measurements. Unlike other data interpolation algorithms (e.g. Nakaoka et al., 2013; Rödenbeck et al., 2014; Sasse et al., 2013; Takahashi et al., 2009), the SOM‐FFN method can produce nearly bias‐free estimates on a global scale while retaining a reasonable amount of fine‐scale structure in the observations (Landschützer et al., 2014). The version used here is the latest version (2020v) ranging from 1982 to 2019, which is publicly available at https://www.ncei.noaa.gov/data/oceans/ncei/ocads/data/0160558/MPI_SOM‐FFN_v2020/. For a more detailed description of this reconstruction product, please refer to Landschützer et al. (2013).

The multi‐dimensional ensemble mode decomposition (MEEMD) method is an adaptive analysis method that was developed to extract spatiotemporal signals from high‐ to low‐frequency timescales (Wu et al., 2009). MEEMD was developed based on empirical mode decomposition (EMD) (Huang et al., 1998) and ensemble empirical mode decomposition (EEMD) (Wu & Huang, 2009). Unlike almost all previous signal analysis methods, EMD/EEMD/MEEMD methods are adaptive and appropriate for both linear and nonlinear time series analysis.

The spatiotemporal data x(s, t) are sliced into temporal data series x(t) at any spatial grid point s by EEMD. All sliced data in our study can be interpreted as a series of amplitude‐frequency modulated oscillatory components (intrinsic mode function, IMF) taken from high‐frequency to low‐frequency timescales and having a residue trend R as follows: where SemiA(t) is the high‐frequency component with a time scale of several months to half a year, MAC (t) is the modulated annual cycle with a quasi‐annual period, in which amplitude and frequency change with time, InterA (t) is the interannual timescale component with a period of several years, Dec (t) is the decadal component with a period of approximately 10 years, and R(t) is the residual trend. Unlike a linear fitted trend, the time‐varying trend obtained by EEMD is an intrinsically fitted monotonic function or a function in which there can be at most one extremum within the whole data span (Wu et al., 2007, 2011). For MEEMD, we then piece together identical timescale components of data series from all grids to form a temporal evolution of the spatially coherent structure of that timescale. It should be noted that the decomposed IMFs and R(t) by EEMD and MEEMD require no a priori function form (basis). All the time‐varying components and residue trends have low sensitivity to the addition of new data and would be included in the cycles with the extension of the data length (Wu et al., 2009, 2011).

RESULTS AND DISCUSSION

Global air–sea CO2 flux in different timescales

From the integrated global air–sea CO2 flux time series (Figure 1), we can see substantial seasonal and decadal variabilities (Figure 1a). The global ocean carbon sink remained stable before 1990 and weakened through much of the 1990s; then, the ocean sink started to strengthen substantially during the 2000s, followed by a slightly increased level of carbon uptake during the 2010s. These decadal fluctuations can be well captured by the time‐varying trend with decadal fluctuations (Figure 1a, red). Here, we further decomposed the original global integrated ocean carbon sink into several IMFs ranging from high‐frequency to low‐frequency: semi‐annual (SemiA), modulated annual cycle (MAC), interannual (interA), and decadal (dec) functions, all of which are followed by a residual trend. These components account for approximately 22% (SemiA), 46% (MAC), 7% (InterA) and 25% (Dec) of the total extracted components, respectively. As expected, the MAC and decadal fluctuations are the two main variabilities for the whole variations, which is also consistent with the previous studies (Landschützer et al., 2014, 2016).

FIGURE 1

Time series of the integrated global air–sea CO2 flux and the components extracted by EEMD. (a) The lines show the integrated global air–sea CO2 flux derived from the SOM‐FFN estimate (black, updated from Landschützer et al., 2016), the 12‐month running mean (blue), and the extracted time‐varying trend with decadal fluctuations (red); (b–f) the extracted IMFs from high‐frequency to low‐frequency: semi‐annual (b), MAC (c), interannual (d), decadal (e), and residual trend (f). EEMD, ensemble empirical mode decomposition; IMF, intrinsic mode function; MAC, modulated annual cycle

As the global patterns of ocean carbon sinks are spatiotemporally different, we further calculated the proportion of all the extracted components to the total variance for each grid during the time period 1982 to 2019 (Figure 2). Compared with other components, semiannual variability is more homogeneously in terms of spatial distribution, with the proportion of most regions ranging between 20% and 25% (Figure 2a). Large proportions of semiannual values are sporadically distributed among the West Indian Ocean, part of the North Atlantic Ocean, and along the Antarctic Front. MAC exhibits the largest variance, varying from 20% to 80%, and the regions having the largest proportions are mainly found at temperate latitudes (10°–40° for both hemispheres) (Figure 2b). The spatial distribution of the proportion of interannual changes is relatively simple and is mainly concentrated in the equatorial Pacific Ocean (Figure 2c), which is believed to be driven by ENSO (Landschützer et al., 2014). The decadal variability of the air–sea flux is not driven by a single region alone (e.g., the Southern Ocean); instead, this variability is largely distributed in the equatorial Pacific Ocean and high‐latitude regions (Figure 2d).

FIGURE 2

Proportions of semiannual (a), MAC (b), interannual (c), and decadal (d) components to total variance (Note that the range of color bars is different). MAC, modulated annual cycle

EMD/EEMD/MEEMD provides a novel time‐frequency‐energy representation of any given time series. Contrary to the previous traditional methods, they can adaptively extract components from high‐ to low‐frequency at different timescales without any “shape function” assumptions. The MEEMD can well preserve both the temporal and spatial localities of the data and thus have advantages when examining both the static climate mode and the changing spatial structure of climate variability (Ji et al., 2014). Figure 2 clearly shows that the information extracted from the data does reflect the physical process during the given time period from 1982 to 2019. Due to the combination of the spatial and temporal localities of MEEMD, the subsequent evolution of the system cannot alter the reality that has already happened. Therefore, MEEMD is an effective method for facilitating comparisons of the global ocean carbon sink at different timescales and providing a visual impression of the percentages of changes.

Seasonality

Previous studies have pointed that the winter‐to‐summer difference has increased substantially in recent decades with the rate increasing at an average of 2.2 ± 0.4 μatm per decade from 1982 to 2015 poleward of 10° latitude, inducing stress in ocean ecosystems and fisheries (Doney et al., 2009; Hauck & Völker, 2015; Landschützer et al., 2018; McNeil & Sasse, 2016). However, from Figure 1c, in contrast with the mean seasonal cycle, we can see that the amplitude and frequency of the modulated annual cycle change with time. In particular, the amplitude of the MAC shows large seasonal differences during the 1990s, a weakening trend during the 2000s, and then a substantial increasing trend during the 2010s. It means that the seasonal difference in the air–sea CO2 flux does not increase at all times or everywhere.

To clarify the spatiotemporal distribution of seasonal difference, we first checked the whole pattern of the seasonal difference from the global integrated air–sea CO2 flux time series (Figure 3). The air–sea fluxes in July, August, and September are generally similar, which represent the boreal‐summer variations. Flux changes in other months are relatively consistent, which can be represented as the boreal non‐summer seasonal variations. The variations in these two periods show decadal fluctuations, with the difference between summer and non‐summer being large during the 1990s and decreasing and stabilizing during the 2000s. After 2010, unlike changes during the 1990s and 2000s, the changes between summer and non‐summer are inconsistent; the sink in boreal‐summer remained unchanged, whereas the sink in boreal non‐summer strengthened over time. It should be noted that the monthly changes corresponding to boreal summer and non‐summer are also relatively concentrated and consistent.

FIGURE 3

Time series of the global integrated air–sea CO2 flux for each month. The black line represents the average flux for July, August, and September. The red line is the average flux for the rest of the months

We further examine the spatial distribution of the boreal summer minus boreal non‐summer (SMNS for short) of air–sea CO2 flux. It shows a zonally distributed structure (Figure 4e). At the moderate and low latitudes (Figure 4a,c) because the ocean sink of CO2 is a negative value, the SMNS at northern mid‐latitudes is positive and the southern hemisphere is negative. The seasonal flux has a negative maximum in non‐summer seasons, which is due to an increase in the dissolution rate of CO2 as a result of a decrease in sea surface temperature (Keppler et al., 2020). It is worth noting that the ocean sink remains stable in summer and that in non‐summer seasons is continuously strengthening at mid‐latitudes in both hemispheres. However, the strengthening non‐summer trend and seasonal difference at mid‐latitudes in both hemispheres also show distinct patterns. The non‐summer sink at latitudes 10°–45°N strengthens more linearly, and the corresponding months are concentrated. At latitudes 10°–45°S, the strengthening mainly occurs after 2000 for boreal‐summer seasons, with the contributing months being scattered. This strengthening of the seasonal SMNS CO2 flux is consistent with previous studies (Landschützer et al., 2018), which is probably due to increasing atmospheric CO2. One spontaneous question is that why only the sink in non‐summer seasons for both hemispheres strengthened, which needs further investigation.

FIGURE 4

Time series of air–sea CO2 flux for each month (unit: PgC/year), (a) 10°–45°N, (b) 45°S−65°N, (c) 10°–40°S, (d) 40°–65°S. The black and red lines in each panel represent the averaged boreal summer (July, August, and September) and non‐summer seasons (the rest of the months). (e) The mean boreal SMNS of air–sea CO2 flux for each 1° × 1° pixel. Positive SMNS are marked in red, while negative SMNS are marked in blue. SMNS, summer minus boreal non‐summer

In contrast, the air–sea CO2 flux at high latitudes (Figure 4b,d) shows a negative maximum in boreal summer, leading to a negative SMNS in north high latitudes and a positive value in the Southern Ocean. Gregor et al. (2018) argued that biological activity drives Southern Ocean carbon sink variability in austral summer, and wind stress in austral winter. In austral winter, intensified westerly winds enhanced outgassing in upwelling regions, causing the carbon sink to weaken during this season (Keppler, 2020). It is intriguing that the seasonal difference in the Southern Ocean is the main characteristic of that in global integrated air–sea CO2 flux. The seasonal difference in the Southern Ocean is not always strengthened with increasing atmospheric CO2, but shows decadal fluctuations.

With the help of the EEMD method, we can extract the modulated annual cycle with the amplitude and frequency varying with time. We compare the MACs in these three decades: 1990s, 2000s, and 2010s (Figure 5). There is no obvious shift in the peak value, whereas the magnitude of the maximum in the 2010s increased by 0.15 PgC/year compared with the 2000s and 1990s, and the duration of trough time increased during the 2010s. It indicates that the processes controlling the air–sea CO2 flux during the 2010s are probably different from those during the other two decades. A detailed discussion of decadal variations can be found in Section 3.4.

FIGURE 5

The average MACs for three decades, the 1990s (blue), 2000s (red), and 2010s (yellow). MAC, modulated annual cycle

Interannual variability

The equatorial ocean, as the ocean's largest natural source of CO2 to the atmosphere, plays an important role in the global carbon cycle (Feely et al., 2006). The annual contribution of CO2 to the atmosphere from the oceanic equatorial belt is estimated to be between 0.6 and 1.5 PgC (Takahashi et al., 2002). In particular, the central and eastern Pacific is a major source of CO2 to the atmosphere during non‐El Niño periods.

Our extracted interannual variability is mainly located in the equatorial Pacific Ocean (Figure 2c), with a proportion ranging from 20% to 30%. The interannual air–sea CO2 flux in the equatorial Pacific Ocean (especially in the central and eastern Pacific Ocean) is negatively modulated by ENSO events (Figure 6) (R = −.79), which is consistent with previous studies (Feely et al., 2004; Takahashi et al., 2002). During El Niño periods, the weakened upwelling results in less CO2‐rich water near the surface and leads to a near cessation of the outgassing of CO2. On the contrary, during La Niña events, the unusually strong upwelling of cold, CO2‐rich water leads to anomalously strong outgassing (Feely et al., 2006). From Figure 6 (black dotted line), we can see that the interannual variability in air–sea CO2 flux in the equatorial Pacific Ocean mainly occurs in the central and eastern Pacific Ocean, such that a one‐degree change in Nino3.4 SST anomaly corresponds to 0.08 PgC/year interannual variation.

FIGURE 6

The interannual variability of air–sea CO2 flux in the equatorial Pacific Ocean (120°E–90°W, 10°S–10°N) (blue line) and the Niño 3.4 SST anomaly (orange). The black dotted line and the grey dashed line represents the contributions of the central/eastern equatorial Pacific Ocean and the western equatorial Pacific Ocean, respectively

Decadal fluctuations

Previous studies have noticed that the global ocean carbon sink weakened during the 1990s with a minimum uptake of only approximately −0.8 PgC/year in 2000 and thereafter strengthened considerably to a level of more than −2.0 PgC/year in 2010 (Landschützer et al., 2016). After 2010, changes in the ocean sink remained relatively stable, increasing slightly to approximately −2.3 PgC/year in 2019 (Figure 1a). CO2 uptake by the Southern Ocean (<35°S) was believed to contribute to decadal global ocean carbon sink variations for the first two decades (Gruber et al., 2019). The question is what the role of the Southern Ocean has been during the last decade.

As the time‐varying trend with decadal fluctuations extracted from EEMD could perfectly capture low‐frequency variations throughout the entire time series (Figure 1a), we further examine the decadal spatial distribution by extracting the time‐varying trend with decadal fluctuations at each grid point (Figure 7). Since the integrated global air–sea CO2 flux shows obvious decadal changes from the weakening decade during the 1990s, the reinvigoration period during the 2000s and the reinforce period during the 2010s, and since the changes in each decade as a whole are relatively consistent, we use the difference between the end point and the starting point of each decade to calculate the magnitude change for the corresponding decade. Figure 7 shows the magnitude change of the air–sea CO2 flux over three decades: the 1990s, 2000s, and 2010s.

FIGURE 7

Decadal magnitude change of air–sea CO2 flux (a–c), westerly wind anomaly (d–f), and sea surface temperature (g–i) over three decades: the 1990s (left), 2000s (middle), and 2010s (right), respectively

Air–sea CO2 flux varied substantially over these three decades. During the weakening decade in the 1990s, the weakening ocean sink resulted primarily from the outgassing of respired CO2 in the high latitude regions, especially in the Southern Ocean (south of 40°S) (Figure 7a). In the reinvigoration period during the 2000s (Landschützer et al., 2015), even though the CO2 source region in the equatorial Pacific Ocean increased, the Southern Ocean, as well as the high latitude and sub‐polar regions in the Northern Hemisphere, contributed to this reinvigoration of the carbon sink (Figure 7b). During the recent decade of the 2010s, the whole global carbon sink remained stable and showed a slight increase (Figure 7c). This slight increase in the carbon sink was due to the weakening of the source region in the equatorial Pacific Ocean, as the Southern Ocean (especially south of 60°S) showed a slight weakening trend.

Although the potential drivers of carbon uptake variability in the Southern Ocean are still debated (Bronselaer et al., 2020; DeVries et al., 2017; Landschützer et al., 2015; Le Quéré et al., 2007), the interannual and decadal variability of the Southern Ocean carbon sink was believed to be closely correlated with changes in surface winds (Kepper & Landschützer, 2019). We further compare the decadal ocean sink fluctuations with the westerly wind anomaly and SST variations during these three decades (Figure 7d–i). The changes in westerly winds are consistent with the decadal fluctuations of the carbon sink in the Southern Ocean. Strengthened westerly winds lead to enhanced outgassing in the Southern Ocean (especially south of 60°S) which is probably attributed to stagnation in the 1990s (Figure 7d). In addition, weakened westerly winds during the 2000s have also been linked to enhanced CO2 uptake during that period (Figure 7b,e). Unlike the 1990s and 2000s, strengthened westerly wind anomalies were mainly concentrated in the western equatorial Pacific Ocean, leading to the warming in the central and eastern Pacific Ocean (Figure 7f,i). Therefore, this process weakened CO2 emissions in the equatorial Pacific Ocean (Figure 7c).

Based on the above analysis, decadal fluctuations are unevenly spatially distributed, with the following four regions predominating: the northern Pacific Ocean (NP, 40°–70°N, 120°W–120°E), the northern Atlantic Ocean (NA, 40°–80°N, 70°W–30°E), the equatorial Pacific Ocean (EP, 10°S–10°N, 180°W–80°W), and the Southern Ocean (SO, 35°–90°S, 180°W–180°E). The sum of these four regions could capture the decadal fluctuations of the entire dataset well (Figure 8), although there remains a discrepancy of 0.6–0.8 PgC/year. Analyzing the relative contributions of these four regions to global air–sea CO2 flux can provide initial insights into the processes driving these decadal variations (Figure 9). During the 1990s, the magnitude of change in decadal global air–sea CO2 flux reached approximately 0.5 PgC/decade with the Southern Ocean accounting for approximately 70%. During the first decade of the 2000s, the magnitude of change reversed to less than −1.0 PgC/decade with the Southern Ocean accounting for approximately 45%. A striking difference is the relative importance of the Pacific Ocean to total flux during the 2010s. The equatorial Pacific Ocean has replaced the Southern Ocean as the main contributor, accounting for approximately 35% of total decadal change (−0.35 PgC/decade).

FIGURE 8

Decadal fluctuations of the regional mean air–sea CO2 flux in the northern Pacific Ocean (NP), the northern Atlantic Ocean (NA), the equatorial Pacific Ocean (EP), the Southern Ocean (SO), and global oceans

FIGURE 9

Decadal differences in air–sea CO2 flux for different ocean regions over three decades. The value is the difference between the end and the beginning of a decadal period of the time‐varying trend with the decadal fluctuations extracted by EEMD. Global, EP, NP, NA, and SP are the global ocean, equatorial Pacific Ocean, the northern Pacific Ocean, the northern Atlantic Ocean, and the Southern Ocean, respectively. EEMD, ensemble empirical mode decomposition

For the decades 1990s and the 2000s, with the Southern Ocean dominating the decadal changes, DeVries et al. (2017) suggested that the ocean carbon cycle in the Southern Ocean was primarily driven by meridional overturning circulation. During the 2010s, accompanied by the regime shift in PDO during the early 2010s, the central and eastern Pacific Ocean showed warm signals during the decade 2010–2019 (Figure 7g,i); thus, carbon emissions in the central Pacific Ocean decreased. The tropics can also influence climate variations at high latitudes through atmospheric teleconnections and thus change regional air–sea CO2 flux (Gruber et al., 2019; Karoly, 1989). The interannual to decadal timescale signals present in the tropics can rectify the high‐frequency variations in the high latitudes to generate low‐frequency (decadal) variations. However, such teleconnections between the tropics and high latitudes through atmospheric bridges are not well established and need further investigation.

CONCLUSION AND DISCUSSION

In this study, we revisit the seasonal to decadal variations of the global ocean carbon sink based on a reconstructed observation‐based product of the surface ocean pCO2 and air–sea CO2 flux (the MPI‐SOMFFN method) using an adaptive data analysis method. The global ocean carbon sink for atmospheric CO2 varies on all spatial and temporal scales. Two predominant variations are modulated annual cycles and time‐varying decadal variations, which account for approximately 46% and 25% of all extracted components, respectively. The dominant mode of variability is the MAC at mid‐latitudes and decadal in tropical and high‐latitude regions.

The global average SMNS can reach 0.5–1.0 PgC/year, with the whole SMNS pattern being dominated by the Southern Ocean. Even though air–sea CO2 flux at mid‐latitudes remains relatively stable in boreal summer, it shows an approximately linear increasing trend (a total increase of 1.0 PgC during the period 1982–2019) for the boreal non‐summer CO2 sink. Compared with mid‐latitudes in the Southern Hemisphere, the corresponding periods of the non‐summer sink at northern mid‐latitudes are more concentrated. This strengthening SMNS of CO2 flux is believed to be correlated with increasing atmospheric CO2, with non‐summer seasons for both hemispheres dominating.

On decadal timescales for the global ocean carbon sink, unlike the weakening decade (1990–1999) and the reinvigoration decade (2000–2009), during which the Southern Ocean played the dominant role, the reinforcement decade (2010–2019) mainly results from the weakening source effect in the equatorial Pacific Ocean. The equatorial Pacific Ocean, which replaced the Southern Ocean as the major contributor, coincides with the phase transition of the PDO. However, the mechanisms driving different periods need further investigation. It remains a fascinating and largely unresolved puzzle how the interaction of all processes can generate such dominant decadal variabilities in the global ocean carbon sink. Our results suggest that except for the Southern Ocean's role in the global ocean carbon sink, the strengthening non‐summer's sink at mid‐latitudes and the decadal or longer timescales of the equatorial Pacific Ocean dynamics should also be fully considered in understanding the oceanic carbon cycle on a global scale.

It should be noted that our results are only based on one dataset using the SOM‐FFN method (Landschützer et al., 2019). Different methods have their own strengths and weaknesses (Denvil‐Sommer et al., 2019; Gregor & Gruber, 2021; Gregor et al., 2019; Rödenbeck et al., 2014; Takahashi et al., 2014 etc). Surface Ocean CO2 Mapping (SOCOM) collated 14 methods in an intercomparison of “gap‐filling” methods (Rödenbeck et al., 2015). They found that the MPI‐SOMFFN (Landschützer et al., 2016) and the Jena‐MLS methods (Rödenbeck et al., 2014) present better interannual variability compared to other methods in the comparison (Gregor et al., 2019). In data‐poor regions (e.g., the Southern Ocean) and the mesoscale intensive ocean regions, they found little agreement between these methods on seasonal and decadal timescales (Ritter et al., 2017), with the amplitude of decadal variability being overestimated by 31% in the Southern Ocean (Gloege et al., 2021). Gregor et al. (2019) developed an ensemble average of six machine‐learning models (CSIR‐ML6) to analyze this problem. The intercomparison among these products is needed for further investigation. Additionally, to improve observational reconstructions of the global ocean carbon sink in the future, additional observations will be critical.

CONFLICT OF INTEREST

We declare that we have no conflicts of interest.