Interannual global carbon cycle variations linked to atmospheric circulation variability

One of the least understood temporal scales of global carbon cycle (C-cycle) dynamics is its interannual variability (IAV). This variability is mainly driven by variations in the local climatic drivers of terrestrial ecosystem activity, which in turn are controlled by large-scale modes of atmospheric variability. Here, we quantify the fraction of global C-cycle IAV that is explained by large-scale atmospheric circulation variability, which is quantified by spatiotemporal sea level pressure (SLP) fields. C-cycle variability is diagnosed from the global detrended atmospheric CO₂ growth rate and the land CO₂ sink from 16 dynamic global vegetation models and two atmospheric inversions in the Global Carbon Budget 2018. We use a regularized linear regression model, which represents a statistical learning technique apt to deal with the large number of atmospheric circulation predictors (p≥800, each representing one pixel-based time series of SLP anomalies) in a relatively short observed record (n<60 years). We show that boreal winter and spring SLP anomalies allow predicting IAV in the atmospheric CO₂ growth rate and the global land sink, with Pearson correlations between reference and predicted values between 0.70 and 0.84 for boreal winter SLP anomalies. This is comparable to or higher than that of a similar model using 15 traditional teleconnection indices as predictors. The spatial patterns of regression coefficients of the model based on SLP fields show a predominant role of the tropical Pacific and over Southeast Asia extending to Australia, corresponding to the regions associated with the El Niño–Southern Oscillation variability. We also identify another important region in the western Pacific, roughly corresponding to the West Pacific pattern.

We further evaluate the influence of the time series length on the predictability of IAV and find that reliable estimates of global C-cycle IAV can be obtained from records of 30–54 years. For shorter time series (n<30 years), however, our results show that conclusions about CO₂ IAV patterns and drivers need to be evaluated with caution. Overall, our study illustrates a new data-driven and flexible approach to model the relationship between large-scale atmospheric circulation variations and C-cycle variability at global and regional scales, complementing the traditional use of teleconnection indices.

1 Introduction

The global carbon cycle (C-cycle) varies at multiple timescales ranging from minutes to a millennium (Ciais et al., 2013). Quantifying and understanding the patterns of variability in the C-cycle and their drivers are crucial to better understanding C-cycle dynamics and better constraining future climate projections (Cox et al., 2013; Friedlingstein et al., 2014). Primarily driven by the land sink dynamics (Piao et al., 2020), interannual variability (IAV) is one of the most uncertain and poorly understood terms of the global C-cycle in the observational period (Friedlingstein et al., 2019).

A fundamental challenge is that variability in land–atmosphere carbon exchange is affected in complex ways by large-scale atmosphere circulation modes but also land use change, forced climate changes, and direct physiological CO₂ effects on ecosystems, among others (IPCC, 2013). Separating these effects is difficult because of the large uncertainties associated with some processes, for example land use change (Friedlingstein et al., 2022) or processes not directly observable at global scale (e.g., photosynthesis or heterotrophic respiration) (Schimel et al., 2015; Basile et al., 2020). The second challenge is that the land sink, as a balance of carbon uptake and release, responds differently to variations in the climatic forcing (Jung et al., 2017; Piao et al., 2020). This makes it hard to attribute induced land sink IAV to specific drivers, which is crucial for process understanding (Jung et al., 2017; Humphrey et al., 2018, 2021; Wang et al., 2022). Last, the limited length of observational records may hamper robust statistical analysis (IPCC, 2013): the longest continuous observations of atmospheric CO₂ at the South Pole and Mauna Loa observatory only exist from 1958 onwards (Dlugokencky and Tans, 2019).

Land biospheric CO₂ uptake results from the net balance of carbon uptake from gross primary productivity, release from multiple respiration terms, and disturbance-induced fluxes such as fires, amongst other smaller terms (IPCC, 2013). Most of these processes are primarily driven by temperature, water, and radiation availability (Jung et al., 2017). These meteorological drivers are, in turn, modulated by large-scale modes of atmospheric circulation on multiple timescales, such as the El Niño–Southern Oscillation (ENSO; Gu and Adler, 2011) and the Pacific Decadal Oscillation (PDO; Newman et al., 2016). These climate variability modes are generated within coupled atmosphere–ocean systems (Ghil, 2002) and considered irreducible noise in climate projections (Madden, 1976; Schneider and Kinter, 1994; Deser et al., 2012). Because these modes typically interact and affect weather dynamics in regions beyond those where they emerge, such modes are collectively referred to as teleconnections (IPCC, 2013). Bacastow (1976) showed that ENSO is highly correlated with annual variations in observed atmospheric CO₂ measured at the South Pole and Mauna Loa, Hawaii. Keeling et al. (1995) attributed these correlations to the ENSO impact on the biospheric sink. In addition to ENSO, Zhu et al. (2017) showed that the PDO and the Atlantic Multidecadal Oscillation (AMO; Enfield et al., 2001; Rayner et al., 2003) may also influence global terrestrial ecosystem carbon fluxes and that other modes of variability in the Northern Hemisphere also have local impacts on carbon cycling (Zhu et al., 2017).

A common approach to diagnose the impacts of natural climate variability is to use ensembles of Earth system model simulations with perturbed initial conditions to quantify the impacts of natural climate variability at decadal to millennium scales (Frölicher et al., 2013). However, the inherently chaotic atmosphere, in combination with model structural uncertainty, implies large uncertainties for future projections (Deser et al., 2020). In addition, Earth system model projections can be compromised by limited representation of the full complexity of physical processes involved, lack of observational constraints, and high computational demands when aiming to resolve high resolutions (Randall et al., 2007; IPCC, 2013).

Statistical approaches are a simplified but effective way to reveal physical processes in observations (von Storch, 1995). A traditional approach consists of evaluating relationships between the variables of interest (e.g., CO₂ time series) and teleconnection indices (Bacastow, 1976; Bastos et al., 2013; Zhu et al., 2017). As a simple representation of the large-scale atmospheric circulation modes, teleconnection indices are extracted mainly from sea surface temperature or atmospheric anomalies (Kumar and Hoerling, 1997; IPCC, 2013). Such indices are an effective way to reduce the complexity of the spatiotemporal variability in multiple variables (Stenseth et al., 2003; Wills et al., 2017), but they may not be able to capture spatial variations in the large-scale atmospheric circulation modes themselves.

Recently, Sippel et al. (2019) applied ridge regression, a regularized linear regression method (Hastie et al., 2009; Friedman et al., 2010), to quantify the component of precipitation and temperature variability driven by atmospheric variations based on sea level pressure (SLP) fields rather than teleconnection indices. Their approach allowed them to robustly infer the main spatiotemporal patterns of atmospheric variability influencing these two climate variables. On the one hand, including a field of circulation-based predictors, avoids considering predefined assumptions about their spatial configurations as they are common to teleconnection indices, while compensating for relatively short historical records. The regularization approach, on the other hand, allows overcoming overfitting and multicollinearity issues due to short time series and a very large number of spatial predictors.

In this study, we adopt the ridge regression approach in Sippel et al. (2019), aiming to quantify the fraction of global C-cycle IAV influenced by large-scale atmospheric circulation variability. We use observation-based time series of global atmospheric CO₂ growth rate (AGR) and land CO₂ surface fluxes from atmospheric inversions and dynamic global vegetation models (DGVMs), as well as the land sink estimated as a residual of other terms in the Global Carbon Budget 2018 (Le Quéré et al., 2018). We first evaluate and compare the predictive skill of the ridge regression model when using SLP as predictors versus commonly used traditional teleconnection indices (Sect. 3.1). Next, we analyze and discuss the global C-cycle sensitivity to atmospheric circulation variability from various latitudinal domains of SLP anomaly fields (Sect. 3.2). Finally, we evaluate the sensitivity of the results to the length of the time series (Sect. 3.3) by comparing the fraction of C-cycle IAV that can be explained by large-scale atmospheric circulation variability based on these datasets with that of a very long time series (4000 years) of land CO₂ fluxes simulated by the Community Earth System Model (CESM).

2 Data and methods

2.1 CO₂ datasets for the recent past

We select the CO₂ time series datasets from the Global Carbon Budget 2018 version 1.0 (Le Quéré et al., 2018): the atmospheric CO₂ growth rate (AGR), the land sink from models (SL_DGVMs), the residual land sink (SL_Resid), and the land sink from two atmospheric inversions.

In the Global Carbon Budget 2018 (Le Quéré et al., 2018), the global CO₂ balance is calculated based on the carbon emissions from fossil fuel (FF) (Boden et al., 2017; UNFCCC, 2018; Peters et al., 2011) and land use change (FLUC) (Houghton and Nassikas, 2017; Hansis et al., 2015), the AGR (Dlugokencky and Tans, 2018), and the carbon uptake by the ocean sink (SO) and the land sink (SL) (references for individual models of SO and SL can be found in Table 4 of Le Quéré et al., 2018).

The difference of annual atmospheric CO₂ in a given year and the previous year (Ballantyne et al., 2012; Dlugokencky and Tans, 2018; Le Quéré et al., 2018) corresponds to the AGR, which is based on direct observations. The AGR is based on the average of well-mixed CO₂ measurements at multiple global stations from the US National Oceanic and Atmospheric Administration Earth System Research Laboratory (NOAA ESRL) (Dlugokencky and Tans, 2018).

FF emissions are based on inventories, while FLUC, SL, and SO are estimated by models (the latter two are denoted by SL_DGVMs and SO in Eq. 1), all of which contain uncertainties (Le Quéré et al., 2018). The total emissions from FF and FLUC minus AGR should equal the SO and SL_DGVMs (Eq. 1). Due to uncertainties in modeled land and/or ocean sinks or in land use estimations (Bastos et al., 2020; Hauck et al., 2020), the budget cannot be balanced, and thus an imbalance term (IMB) is introduced to the budget.

$$\begin{array}{}\text{(1)}& \mathrm{FF}+\mathrm{FLUC}-\mathrm{AGR}-\mathrm{SO}={\mathrm{SL}}_{\mathrm{DGVMs}}+\mathrm{IMB}={\mathrm{SL}}_{\mathrm{Resid}}\end{array}$$

The annual land sink of CO₂ (SL_DGVMs) is the average net biome production (NBP) simulated by 16 dynamic global vegetation models (DGVMs) forced with historical CO₂ and changing climate (Le Quéré et al., 2018). The residual land CO₂ (SL_Resid) is calculated from emissions, AGR, and ocean sinks, as described in Eq. (1). SL_Resid corresponds to the balance of the fossil fuel and land use change emissions as well as the sinks in the atmosphere and ocean, and it provides an alternative estimate of the global land sink.

The time series of AGR_R, SL_DGVMs, and SL_Resid in the Global Carbon Budget 2018 (Le Quéré et al., 2018) are provided at annual time steps over the period 1959–2017. In the following analysis, we invert the AGR time series (AGR_R for reversed AGR, i.e., $-\mathrm{1}\times \mathrm{AGR}$) for sign consistency with the land sink datasets used (defined as a positive flux from the atmosphere to the land).

Additionally, we use the globally aggregated net atmosphere to land CO₂ flux (positive sign as a sink in the biosphere) estimated from two atmospheric CO₂ inversions in the Global Carbon Budget 2018 (Le Quéré et al., 2018): the Jena CarboScope SL_CarboScope (Rödenbeck, 2005; Rödenbeck et al., 2018) and the Copernicus Atmosphere Monitoring Service inversion SL_CAMS (Chevallier et al., 2005), which cover the periods 1976–2017 and 1979–2017, respectively. Here we use the global annual CO₂ fluxes of these two inversions adjusted for fossil fuel emissions and lateral fluxes from Bastos et al. (2020). The period common to the CO₂ time series (1980–2017) is selected.

2.1.1 Sea level pressure

We use global monthly mean SLP fields from ERA5 reanalysis with the spatial resolution of $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ (Bell et al., 2020) at monthly time steps covering the period 1950–1978 (Bell et al., 2020) and 1979–present (Hersbach et al., 2019). The period common to other datasets of 1958–2017 is selected here.

2.1.2 Teleconnection indices

In addition to SLP fields, we select 15 teleconnection indices from the atmosphere–ocean variability, Northern Hemisphere, and Southern Hemisphere.

Three important atmosphere–ocean coupled variability modes influence global climate and the C-cycle: the El Niño-Southern–Oscillation (SOI), the Pacific Decadal Oscillation (PDO), and the Atlantic Multidecadal Oscillation (AMO) (Zhu et al., 2017).

In the Northern Hemisphere, the most relevant indices are the Arctic Oscillation (AO), the North Atlantic Oscillation (NAO), the Pacific North American pattern (PNA), the East Atlantic (EA), East Atlantic/West Russia (EAWR), the Scandinavian pattern (SCAND), polar–Eurasia (polarEA), and the West Pacific (WP). These indices are calculated and provided by the Climate Prediction Centre (CPC) of the National Oceanic and Atmospheric Administration (NOAA) (CPC, 2008). Detailed information on calculation procedures is described in NOAA CPC (2008) and Barnston and Livezey (1987).

In the Southern Hemisphere, important indices are the Antarctic Oscillation (AAO), the tropical Atlantic dipole (TAD), the Dipole Mode Index (DMI) of the Indian Ocean Dipole, and the Trans Polar Index (TPI).

The teleconnection indices used here have been summarized in Table 1. All the indices are provided as monthly means and selected for the period of 1958–2017 except the AAO, which is available for 1979–2017 only.

(McBride and Nicholls, 1983; Ropelewski and Jones, 1987)(Mantua et al., 1997; Mantua and Hare, 2002)(Mantua et al., 1997)(Rayner et al., 2003; Enfield et al., 2001)(Enfield et al., 2001)(Higgins et al., 2000, 2002)(Barnston and Livezey, 1987)Barnston and Livezey, 1987Chen and Van den Dool, 2003(Barnston and Livezey, 1987)(Barnston and Livezey, 1987)(Barnston and Livezey, 1987)(Barnston and Livezey, 1987)(Barnston and Livezey, 1987)(Mo, 2000)(CPC, 2008)Enfield et al., 1999Reboita et al., 2021(Saji and Yamagata, 2003; Reboita et al., 2021)(Pittock, 1980, 1984; Jones et al., 1999)Henley et al. (2015)

Table 1 Teleconnection indices.

2.1.3 Long-term pre-industrial control simulations for statistical benchmarking

Here we select the SLP fields and global net biome production fields (NBP) from simulations by the Community Earth System Model (CESM) version 1.2.2 (in the B1850C5CN configuration), which has been used by Stolpe et al. (2019). This experiment corresponds to a 4000-year control run. The simulation was run at an atmospheric resolution of $\mathrm{1.9}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$ using the Community Atmosphere Model version 5 (CAM5.3; Neale et al., 2012) with 30 vertical levels. The model consists of fully coupled atmosphere, ocean, sea ice, and land surface components (Hurrell et al., 2013; Meehl et al., 2013b) and did not include dynamic vegetation. This simulation includes no external forcing, so it is ideal to analyze patterns driven by internal variability.

2.2 Data pre-treatment

For all historical datasets (CO₂ time series, SLP fields, and teleconnection indices), we first remove years corresponding to volcanic eruptions (1963, 1982, 1983, 1991, 1992). We then pre-treat the datasets as follows.

2.2.1 Trend removal

The long-term trend of CO₂ time series is removed by locally weighted scatterplot smoothing (LOWESS) (Cleveland et al., 1991) of the annual time series with a fixed window size that is 25 % of the longer interval (1959–2017) and 45 % for the shorter period (1980–2017). For monthly teleconnection indices, we first calculate the seasonal mean values of DJF, MAM, JJA, and SON; we then remove the seasonal long-term trends by applying the LOWESS as for the CO₂ time series and further include DJF and MAM combined (DJF+MAM) as treated in SLP (as described below).

2.2.2 Spatial and temporal aggregation

The monthly mean SLP fields are area-weighted and aggregated to ${\mathrm{2}}^{{}^{\circ }}\times \mathrm{2}{}^{\circ }$, $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$, and $\mathrm{9}{}^{\circ }\times \mathrm{9}{}^{\circ }$ spatial resolution, and the seasonal cycle is removed by subtracting the monthly mean values for each pixel. We then aggregate SLP values in seasonal means for December of the previous year to February (DJF), March–May (MAM), June–August (JJA), and September–November (SON) of each given year and further consider DJF and MAM combined (DJF+MAM), so the number of pixel-based time series (predictors) in DJF+MAM is double DJF. Note that a large fraction of the pixel-based time series of seasonal SLP anomalies show no long-term trend, and the predicted differences between LOWESS detrended and non-detrended SLP are small. Here we keep the analysis of SLP anomalies with no LOWESS detrending. We refer to DJF and MAM as boreal winter and boreal spring.

For the CESM simulations, the SLP fields are originally provided at $\mathrm{1.9}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$ spatial resolution at monthly mean time steps, which we then resample to $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ spatial resolution. Annual mean NBP is calculated from the monthly fields. NBP and SLP fields are selected for the simulation period of 1000–5000 years.

2.3 Statistical analysis

The overall goal is to characterize annual variations in the global C-cycle that can be explained by large-scale atmospheric circulation variability. Here, the pixel-based time series of SLP anomalies are used as predictors (p≥800) of CO₂ time series (n≤ 54 years) in a linear regression model. However, the small sample size relative to the large number of predictors (n<p) can cause severe overfitting problems and result in unstable predictions (Hastie et al., 2009). Moreover, the existing spatial correlations among the neighboring pixels of SLP anomalies might cause multicollinearity among the predictors (von Storch and Zwiers, 1999). The potential multicollinearity problem results in unstable ridge regression coefficients in least square estimation, making it difficult to diagnose the most sensitive spatial patterns of predictors (von Storch and Zwiers, 1999).

Sippel et al. (2019) applied ridge regression to avoid these overfitting and multicollinearity problems. Ridge regression is a regularized linear regression, the fundamental principle of which is to introduce a constraint (hyper-parameter λ) to regularize the varying regression coefficients in least squares estimation (Hastie et al., 2009; Friedman et al., 2010). The regularized variance comes with a compromise of biased predictions and is addressed as the bias–variance trade-off (Hastie et al., 2009). When selecting the best hyper-parameter λ, this trade-off is considered to achieve stable (low variance) while slightly biased predictions (Hastie et al., 2009).

Model performance is evaluated by the R², Pearson's correlation R, and mean squared error of the original CO₂ time series against predicted values. Pearson's correlation R is selected as the main measure of predictability, and the significance P<0.05 is selected. Given the relatively short period (n<60), here we use leave-one-out cross-validation to achieve optimal model training and testing. For each training and testing group splitting, we select the training group as all years excluding 3 consecutive years and the middle year of those 3 is then selected as the test sample. We exclude the preceding and following years to reduce the potential influence of temporal autocorrelation.

Figure 1 Schematic representation of the statistical approach and model design, with an example of the selected time interval of 1959–2017. (a) Fundamental principle of ridge regression. (b) Model training and validation under ridge regression leave-one-out cross-validation. (c) Training and testing grouping through leave-one-out cross-validation.

A schematic description of the workflow from model training, validation, and selection with the selected interval 1959–2017 is shown in Fig. 1. In Fig. 1a, the global maps represent an example of the spatial distribution of SLP with a resolution of $m{}^{\circ }\times m{}^{\circ }$, where m varies in [2, 5, 9]. Each pixel corresponds to a time series of SLP, so that p predictors (x_ik for i in 1∼n years and k in 1∼p, where p depends on the spatial resolution selected) with n=54 time steps are defined. Each predictor is assigned a coefficient ω_k to collectively predict CO₂ time series, with length n=54. The cost function (Harrington, 2012) is the sum of all the squared errors of original y_i minus estimated $y_i^{\prime }$ ($x_i^T\mathit{\omega }+b$, ω represents the vector of all ω_k). At the same time, the constraint function (Harrington, 2012) suppresses the coefficient variations under a regularized range defined by the hyper-parameter λ. Figure 1b shows the model training and validation processes. In this study, SLP is aggregated to $\mathrm{9}{}^{\circ }\times \mathrm{9}{}^{\circ }$, so p=800. Note that to reduce the heavy computational load, when conducting the spatial and temporal sensitivity study as described in Sect. 2.4, the range of λ is lower and with smaller steps than the range and step shown in Fig. 1b. For example, when only selecting the tropical domain of SLP, the number of predictors is much less than 800. So λ is selected in the range [10, 1000] with a step of 50. When using teleconnection indices instead of SLP anomalies, the predictors are equal to or less than 15, and the range of λ is selected from [1, 200] with a step of 2.

The first step is to divide datasets into training and testing groups, as shown in Fig. 1c. The grouped training datasets are then used for model training to tune the best λ through 5-fold cross-validation. The best λ that achieves the optimal prediction (highest R²) is then selected by the model to predict test datasets. The model then starts another iteration with training and testing grouping. Figure 1c describes the leave-one-out training and testing grouping. The years before and after the selected test year are removed for each grouping to reduce the impact of temporal autocorrelations in CO₂ time series.

Ridge regression leave-one-out cross-validation is performed using the Python package Scikit-learn “Ridge”, and λ is tuned by Scikit-learn “RidgeCV” (Pedregosa et al., 2011). The global maps (Figs. 3, A4, A11, and A12) are plotted by Cartopy (Met Office, 2010–2015).

2.4 Experimental design

1. Preliminary dependency tests. To evaluate the robustness of the results for different characteristics of the datasets and methodological choices, we perform several preliminary tests. (1) A resolution dependency test is performed to evaluate the sensitivity of results to the SLP spatial resolution under $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$, $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$, and $\mathrm{9}{}^{\circ }\times \mathrm{9}{}^{\circ }$. (2) A seasonality dependency test is performed to evaluate the dependence of results on the definition of particular seasons, with each season being the combination of 3 consecutive months (from November last year to July in the given year). (3) Temporal autocorrelation of the CO₂ time series is performed to ensure no significant trend remains in the detrended CO₂ time series.

   Here we directly use the results following the preliminary dependency test (Appendix A). The spatial resolution does not influence the results considerably (see Fig. A1), and therefore we select $\mathrm{9}{}^{\circ }\times \mathrm{9}{}^{\circ }$ SLP spatial resolution given its smaller number of grid points. The seasonal dependency test shows that DJF and MAM are seasonal combinations more representative of boreal winter and spring (see Fig. A2). JJA and SON are found to have lower or no predictability to CO₂ time series; therefore, we limit our results to DJF and MAM. As shown in Fig. A3, the temporal autocorrelation of all CO₂ time series is mostly less than 0.4 with lag ranging from 1 to 35 years. With a lag of 1 year, absolute values of autocorrelation are below 0.2 so that we can exclude strong temporal autocorrelation effects.

2. Model training and evaluation. We evaluate the predictability of annual CO₂ time series using SLP anomalies, teleconnection indices, and SOI independently in the DJF, MAM, and DJF+MAM seasons with the approach described above. We compare Pearson's correlation of observations and predicted values for different periods (1959–2017 and 1980–2017, Fig. 2) and show the corresponding ridge regression coefficient distribution maps (Fig. 3a).

3. Spatial sensitivity study. We evaluate the predictability of historical annual CO₂ time series using DJF SLP anomalies under different spatial domains in the periods of 1959–2017 and 1980–2017. Then we use a 30-year sliding window with annual AGR_R to depict how the predictability under various SLP domains evolves in the period 1959–2017.

4. Temporal sensitivity study. We evaluate the predictability of annual CO₂ time series AGR_R, SL_DGVMs, and CESM using DJF and MAM SLP anomalies under different time intervals. Sliding windows are employed at time intervals of 15, 20, 30, and 40 years for historical datasets and CESM, as well as 100, 500, and 2000 years for CESM only. For the interval of 100 and 500 years, we use the sliding window of a 50-year step and a 500-year step for the 2000-year interval. The intervals shorter than 100 years are all in the 1-year step. We also evaluate the error rate of the model in each sliding window of 15-, 20-, 30-, and 40-year lengths. The error rate is calculated by the number of invalid predictions that have significance P>0.05 divided by the number of total predictions within a given window.

5. Comparison with empirical orthogonal function analysis. We compare the predictability of AGR_R by applying ridge regression and a generalized linear model that uses the principal components of SLP fields estimated by empirical orthogonal function (see, e.g., von Storch and Zwiers, 1999) decomposition as predictors. We first select three different DJF SLP spatial domains (global, 18^∘ N–18^∘ S, and 18^∘ N–72^∘ S). For each spatial domain, we select the first 10 components of SLP anomalies. When selecting the global SLP field, the first 10 components explain up to 75 % of the SLP variance (see Fig. A14). We then reconstruct SLP fields based on an increasing number of components from 1 (based on component 1) to 10 components. These components are used as predictors of AGR_R in a simple linear regression (also with leave-one-out cross-validation). The periods 1959–2017 and 1980–2017 are included (see Figs. A13 and A14).

Figure 2 (a) Standardized annual observed and modeled CO₂ time series over the period 1959–2017 (AGR_R in black, SL_Resid in red, and SL_DGVMs in yellow) and in the period 1980–2017 (SL_CAMS in light blue and SL_CarboScope in dark blue). The CO₂ time series have all been detrended as described in Sect. 2. Note that the AGR_R, SL_Resid, and SL_DGVMs in the period 1980–2017 are detrended data based on their relevant period; compared with detrended data based on 1959–2017, the difference is negligible. (b) Pearson correlation of predicted vs. observed and modeled CO₂ time series based on the ridge regression with SLP fields (ρ_SLP) or teleconnection indices (ρ_Tele) as predictors. Additionally, the Pearson correlation of predicted vs. observed and modeled CO₂ time series by linear regression is based on the single predictor of the SOI index (ρ_SOI). SLP fields, teleconnection indices, and SOI are aggregated for different seasons: DJF, MAM, and DJF + MAM. Panel (b) shows results for 1959–2017 and panel (c) for 1980–2017. Note that in panel (c), the ρ_SLP of SL_CAMS using MAM SLP as a predictor has significance P=0.09, and all others have significance P<0.05.

Figure 3 (a) Distribution of ridge regression coefficients of SLP with the time series of AGR_R (top row), SL_Resid (center row), and SL_DGVMs (bottom row) in DJF (left column) and MAM (right column) based on SLP fields in the period 1959–2017. (b) Distribution of ridge regression coefficients of teleconnection indices with AGR_R, SL_Resid, and SL_DGVMs. Both ω_SLP and ω_Tele are the mean of the n=54 run ridge regression coefficients.

3 Results and discussion

3.1 Global IAV patterns

In this section, we test the predictability of the global C-cycle IAV by using different predictors: global SLP fields by ridge regression, teleconnection indices by ridge regression, and the single SOI index by simple linear regression. In this context, each detrended annual CO₂ time series (including AGR_R, SL_Resid, SL_DGVMs, SL_CAMS, and SL_CarboScope) is predicted by the above predictors over DJF, MAM, and DJF + MAM separately (Fig. 2). The predictability is evaluated by the Pearson correlation (ρ) between the original and predicted detrended annual CO₂ time series. ρ_SLP, ρ_Tele, and ρ_SOI represent the predictability by using different predictors of global SLP fields, teleconnection indices, and SOI, respectively. Accordingly, the relevant ridge regression coefficients with the above different predictors are represented as ω_SLP, ω_Tele, and ω_SOI.

First, we find that the detrended annual CO₂ time series are generally consistent with each other except that the SL_DGVMs shows slight deviation (Fig. 2a). We find two anomalous years (1987 and 1998), which show deviations larger than 2 standard deviations in most CO₂ time series, both signifying apparent AGR increases and subsequent lower land sink (Fig. 2a). These two years correspond to strong El Niño events, which are usually associated with below-average land CO₂ uptake (Keeling et al., 1995; van der Werf et al., 2004; Bonan, 2016).

Figure 4 Heat map of predictability with CO₂ time series over various SLP latitude domains in DJF. Each heat map contains 5×5 squares, and each square represents one domain of SLP. For example, the square 36^∘ N–72^∘ S is the domain of SLP extending from 36^∘ N to 72^∘ S. All latitudinal domains include the tropical area (18^∘ N–18^∘ S). The top right square thus represents global-scale SLP. ρ_SLP of AGR_R, SL_Resid, SL_DGVMs, SL_CAMS, and SL_CarboScope in 1980–2017 is shown here.

Global SLP and teleconnection indices show comparable predictive skill of global C-cycle IAV in winter, while teleconnection indices have higher predictive skill in spring (Fig. 2b). In both periods, the value of ρ_SLP (except SL_DGVMs) is higher in DJF (0.51–0.70) than in MAM (0.29–0.60). On the other hand, the values of ρ_Tele are higher in MAM (0.57–0.79) than in DJF (0.53–0.70). The relatively low predictive skill of global SLP anomalies compared to teleconnection indices might result from (1) limited sample size (less than 60 years) and a large number of predictors (p=800) for ridge regression training with global SLP anomalies. But for teleconnection indices and SOI, the predictive skills are much less influenced by the limited sample size due to their limited predictors (p≤15 for teleconnection indices and p=1 for SOI). As we increase the sample size to over 100 years, the predictive skill of SLP anomalies increases considerably, as is shown in temporal sensitivity study (Fig. 6), and (2) the predictive skill of SLP anomalies in explaining global C-cycle IAV can be reduced in domains with large local rather than global impacts of atmospheric variations to land carbon sinks (Jung et al., 2017). In such domains, the SLP anomalies might show a strong relationship to local C-cycle variations but a weaker link to global C-cycle variations. Selecting the SLP domains with a higher contribution to the global C-cycle variability could improve the predictability, as is shown by the analyzes of the sensitivity of the results to the SLP spatial domains (Fig. 4).

We find that the ridge regression using the full SLP fields as predictors yields generally higher predictability than that of linear regression based on principal components of SLP (see Figs. A13 and A14). In winter and in the SLP domain of 18^∘ N to 72^∘ S, the linear regression based on the leading principal components of SLP shows lower predictability that the ridge regression for different numbers of components retained (see Fig. A13, ρ<0.7 in 1959–2017 and ρ<0.75 in 1980–2017). The lower predictive skill of the linear regression based on the leading components of SLP might be due to (1) principal component analysis capturing the main variances of the SLP field but not necessarily that of the main patterns that influence CO₂ IAV and the fact that, (2) mathematically, ridge regression and principal component analysis are deeply connected. Principal component analysis cuts off all components with small variance beyond a certain threshold, while ridge regression shrinks them, which allows for information in low-variance components to be used for prediction (van Wieringen, 2021). Thus, ridge regression might reveal hidden components of SLP variability that are nevertheless important for CO₂ IAV.

Compared to the predictive skill of teleconnection indices, which includes a set of 14 teleconnection indices for the period 1959–2017 and 15 for the period 1980–2017 as predictors, the predictive skill of SOI is slightly lower or similar in both seasons, with 0.53–0.67 in DJF and 0.60–0.74 in MAM (Fig. 2). This is consistent with the dominant role of ENSO in driving global C-cycle IAV, with other modes showing fewer contributions. Such interpretation requires caution as the indices cannot fully represent the complex atmospheric dynamics.

The predictive skill of the combined winter and spring global SLP anomalies reveals the different seasonal responses of global C-cycle IAV to large-scale atmospheric circulation variability (Fig. 2). The predictive skill of SLP and teleconnection indices in DJF+MAM is within the values for DJF and MAM for most datasets and slightly higher than the best-performing season for the predictive skill of SOI.

The predictive skill of SLP to SL_Resid is similar to SL_DGVMs in MAM and slightly higher than SL_DGVMs in DJF. The difference in the predictive skill of SLP to SL_Resid and SL_DGVMs in DJF may due to the fact that (1) compared to SL_Resid, land sink IAV simulated by DGVMs is less sensitive to DJF climate forcing (Bastos et al., 2018), and (2) SL_Resid implicitly includes the variability from land use change as well as ocean sink variations (Dufour et al., 2013; DeVries et al., 2017; Friedlingstein et al., 2019).

Figure 5 Predictability of AGR_R with DJF SLP over various latitude domains. A 30-year sliding window in the period 1959–2017 with a 1-year step is created. The starting and end year of each interval is labeled on the top of each heat map. Here we only show the results of every second starting year; the full results are in Fig. A9.

We next compare the spatial patterns of the ridge regression coefficients of SLP and teleconnection indices in the period of 1959–2017; the results of the period 1980–2017 can be found in Figs. A4 and A5. The spatial patterns of the ω_SLP are similar for the three CO₂ time series: positive coefficients over the eastern tropical Pacific Ocean, negative coefficients from Southeast Asia extending to Australia that together are roughly consistent with ENSO, and negative coefficients from the western Pacific (Fig. 3a). In winter the positive coefficients over the eastern tropical Pacific are higher than in other regions, which are influenced by El Niño and La Niña, respectively (Monahan, 2001; Hsieh, 2004; Rodgers et al., 2004; Schopf and Burgman, 2006; Sun and Yu, 2009; Yu and Kim, 2011): El Niño induces negative SLP anomalies over the eastern Pacific and positive SLP anomalies over the western Pacific (see King et al., 2020, Fig. 5). The results are consistent with the land sink being negatively driven by ENSO in winter: strong El Niño, decreased land sink; strong La Niña, increased land sink. In spring, the area over the central and western tropical Pacific also shows stronger coefficients and likely corresponds to a mix of different modes, such as the ENSO, the western Pacific teleconnection, and the Interdecadal Pacific Oscillation, all showing strong coefficients in Fig. 3b (SOI, WP, and TPI indices). In Fig. A11 we show the anomalies in temperature and precipitation associated with these patterns, as well as those in NBP from the two atmospheric inversions (see Fig. A12). Generally, the temperature anomalies over the tropics show negative correlations with annual land sink (SLP-driven AGR_R) in both winter (as high as −0.85) and spring (as high as −0.73), while weaker but positive correlations are found in Eurasia. Tropical precipitation anomalies show roughly positive correlations in winter (as high as 0.73) and in spring (as high as 0.67). This pattern indicates that AGR_R is generally higher for cooler and wetter conditions over the tropics and Southern Hemisphere semi-arid regions in both seasons, which result in increased NBP (see Fig. A12), as well as cooler but predominantly drier conditions over Eurasia, which result in a complex pattern of NBP anomalies (see Fig. A12). These results are consistent with the strong ENSO fingerprint on the IAV of global CO₂ atmospheric growth rate and global land sink, e.g., as pointed out by Piao et al. (2020) and with the importance of southern semi-arid ecosystems (Ahlström et al., 2015), for IAV in the global land sink.

The higher ridge regression coefficients of teleconnection indices are consistent with the high-sensitivity domains of the ridge regression coefficients of SLP corresponding to the patterns of ENSO and WP (Fig. 3b). Our results show that C-cycle IAV reveals a high positive sensitivity to SOI (0.18 to 0.21) and negative sensitivity to TPI (−0.24 to −0.28) and WP (−0.09 to −0.15) in winter. High sensitivities are also found for DMI in winter (negative) and AMO in spring (negative) (Fig. 3b). We find that the global C-cycle IAV is very sensitive to TPI as well as SOI, which is not so obvious for the spatial patterns of SLP ridge regression coefficients. TPI is a hemispheric-scale index defined as the pressure anomaly differences between the locations Hobart (43^∘ S, 147^∘ E) and Stanley (52^∘ S, 58^∘ W) (Pittock, 1980, 1984). We find the TPI to be strongly anticorrelated with SOI in winter and spring (−0.89 and −0.85, respectively). This might indicate an amplification of ENSO impacts on C-cycle IAV due to large-scale atmospheric circulation variability in the Southern Hemisphere.

However, the observed patterns of the ridge regression coefficients in teleconnection indices and SLP need to be interpreted with caution, since these patterns are not necessarily independent from each other. For example, the area from Southeast Asia extending to Australia corresponds to a region influenced by several large-scale atmospheric circulation modes: ENSO, the Indian Ocean Dipole (IOD), and the Southern Annular Mode (SAM) (Cleverly et al., 2016). Interactions between these modes have been shown to modulate the occurrence of drought and extreme precipitation in semi-arid areas of Australia and thus induce large interannual variability in gross primary productivity in the region (Cleverly et al., 2016).

Compared to other CO₂ datasets, the predictability of SL_DGVMs is higher when using SLP fields in MAM as predictors rather than DJF (Fig. 2b). Moreover, ω_SLP of SL_DGVMs exhibits distinct spatial patterns, especially in winter; ω_SLP for SL_DGVMs shows higher values in the southern Pacific than over the tropical Pacific region (Fig. 3a). Compared to historical results, predictability of SL_DGVMs is lower by using winter SLP as predictors in ridge regression rather than by using spring SLP, and the spatial patterns of the ridge regression coefficients for SL_DGVMs are slightly different. These differences might be an indication of shortcomings of DGVMs in simulating the sensitivity of land sink to climatic drivers.

The general match of spatial patterns of the ridge regression coefficients using SLP and the teleconnection indices as predictors of global C-cycle IAV indicates that SLP can capture the spatial distribution of the atmospheric patterns that influence IAV, with the advantage of being more flexible than teleconnection indices, since it does not require predefined definitions. However, the short sample size and the large number of predictors for ridge regression training hinder the performance of SLP anomalies, especially the lower predictability when using SLP anomalies rather than teleconnection indices in spring. Reducing the number of predictors (smaller domains of SLP anomalies) or increasing the sample size (longer time interval) for ridge regression training could improve the predictive skill of SLP anomalies. Therefore, in the next subsections, we conduct a spatial and temporal sensitivity study of the global C-cycle to SLP anomalies.

3.2 Sensitivity to the SLP domains

Here, we test the sensitivity of the global C-cycle IAV predictability by using different spatial domains of SLP as predictors in ridge regression. The SLP domains are selected over different latitudinal bands in DJF and MAM separately (Fig. 4). We find improved predictability in both seasons when selecting smaller spatial domains (particularly including the tropics to high latitudes of the Southern Hemisphere) rather than global SLP anomalies. MAM fields show lower predictability in general and are less sensitive to the spatial domain considered. In the following, we show the results for DJF; the results for MAM can be found in Fig. A7. Here we only show the results of the period 1980–2017. The results of the period 1959–2017 show a similar trend (see Fig. A6).

Consistent with previous studies (Zeng et al., 2005; Piao et al., 2020), the tropical domain corresponds to higher predictability for all datasets, but stronger predictability is found for regions extending from the tropics to the Southern Hemisphere (Fig. 4). Including the Northern Hemisphere regions results in lower predictability. The domain 18^∘ N–72^∘ S shows the highest predictability, with ρ_SLP of 0.81 for AGR_R and 0.85 for SL_Resid in 1980–2017.

The results for net atmosphere–land fluxes estimated by atmospheric inversions are consistent with those of AGR_R, with ρ_SLP of 0.70 for SL_CAMS and 0.76 for SL_CarboScope in the same domain of 18^∘ N–72^∘ S. The values of ρ_SLP of SL_DGVMs are systematically lower than the other datasets, independently of the domain.

The weaker values of predictability when extending SLP domains from the tropics to the Northern Hemisphere (Fig. 4) might be due to the local rather than global impacts of large-scale atmospheric circulation variability in the Northern Hemisphere to land sink IAV. Additional explanations include the fact that carbon fluxes are weaker in the winter Northern Hemisphere so that large-scale atmospheric circulation variability exerts a weaker influence in the global land sink and that there are strong compensatory effects of gross primary productivity versus terrestrial ecosystem respiration in the Northern Hemisphere in response to water and temperature variations (Jung et al., 2017; Wang et al., 2022).

The increasing predictability of global C-cycle IAV when using SLP fields extending from the tropics to the Southern Hemisphere (Fig. 4) is likely due to the strong contribution of semi-arid regions in the Southern Hemisphere extratropics to the global sink through their drought–wet anomalies (Poulter et al., 2014; Ahlström et al., 2015). The drought–wet anomalies in these regions are controlled by large-scale atmospheric circulation variability in the Southern Hemisphere (ENSO and ENSO related modes) and to the interactions between ENSO and other large-scale atmospheric circulation modes in the Southern Hemisphere, such as the synergistic effects from ENSO, IOD, and the SAM on Australia C-cycle variability (Cleverly et al., 2016).

3.3 Sensitivity to the temporal domains

Because of multidecadal variability in the climate system, it is possible that the relationships found for short intervals are not stable. In order to investigate whether these results depend on the temporal domain considered, we additionally analyze the influence of different temporal domains (30-year interval sliding window) on the predictability of global C-cycle IAV by selecting different temporal domains and performing the ridge regression separately for the global and the tropical domains (Fig. 5).

Results show stronger predictability of AGR_R confined to the tropics of SLP domains in earlier periods and an intensification of this predictability for the SLP domains extending to the Southern Hemisphere over the study period. In some periods, the tropics and southern extratropics show the highest values of ρ_SLP, for example in 1978–2011 and 1984–2015 (Fig. 5). There is, however, high temporal variability in the predictability and the most relevant spatial domain, with other periods showing higher global coherence (e.g., 1968–2001). It is unclear whether these temporal variations occur randomly due to internal variability in the climate system or are influenced by external forcing. Potential explanations for this pattern include trends found in SLP variability over the Pacific and southern Atlantic (Schneider et al., 2012; IPCC, 2013; Roxy et al., 2019) or enhanced sensitivity of C-cycle variability to climatic drivers, particularly in semi-arid areas, under progressive climate change (Wang et al., 2014; Poulter et al., 2014).

Understanding and attributing these changes to given processes are beyond the scope of this study, but these results highlight the importance of the temporal domain when analyzing IAV in the global C-cycle. Since the observed CO₂ time series are short and cover only limited temporal domains, results are likely to be affected by multidecadal internal climate variability, in addition to external forcing. Moreover, the data-driven ridge regression method to quantify circulation-induced global C-cycle variability uses a large number of predictors, while only relatively short time series are available for training, which may negatively affect the model's performance. Therefore, we further test the sensitivity of the predictability to the length of the time series (Fig. 6). We test the predictability of global C-cycle IAV for different lengths of the temporal domain: 15, 20, 30, and 40 years for the datasets in the Global Carbon Budget 2018 as well as CESM simulations and 100, 500, and 2000 years for CESM only.

Figure 6 Predictability of AGR_R, SL_DGVMs, and CESM NBP under various time intervals. The predictability of AGR_R and SL_DGVMs is within the period 1959–2017, with a 1-year step sliding window of 15, 20, 30, and 40 years. The predictability of CESM is in the period 1000–5000 and covers extra intervals of 100, 500, and 2000 years. The distribution of the predictability under each sliding window with SLP in DJF global, MAM global, DJF tropical, and MAM tropical is shown. Tropical domains are 18^∘ N–18^∘ S for SLP as predictors in predicting AGR_R and SL_DGVMs, as well as 20^∘ N–20^∘ S for SLP as predictors in predicting CESM. Note that the mean values are shown as black dots.

The box plots in Fig. 6 show the distribution of predictability calculated for multiple time intervals, with each time interval using a sliding window over the whole period of the respective time series of AGR_R, SL_DGVMs, and CESM. The spread of predictability provides an indication of internal variability in the predictability of global C-cycle IAV due to the choice of temporal domain and the uncertainty in the ridge regression fit for a large number of predictors and comparatively small number of training samples.

We find that the longer the time interval the higher the mean predictability and the smaller the variation, i.e., the less dependent the results are on the temporal domain considered (Fig. 6). However, the mean value tends to be lower than the median for intervals shorter than 30 years and similar to the median for longer intervals. The lower mean is influenced by some domains with very low or even negative predictability from invalid predictions in shorter time intervals.

The mean predictability of AGR_R in winter for the global domain increases from 0.45 to 0.61 from 15 to 40 years, respectively, while the spread (maximum ρ_SLP − minimum ρ_SLP) decreases from 1.04 to 0.18. Predictability for SL_DGVMs is consistent with those of AGR_R, with systematically lower mean predictability for winter and higher for spring, but similar spread in both. The mean value of predictability for CESM in SLP DJF over the global domain increases from 0.42 (15 years) to 0.57 (40 years) and to 0.84 (2000 years), and the spread decreases from 1.72 to 0.72 and to 0.008, respectively.

At global scale, the predictive skill of SLP anomalies with AGR_R and with models from SL_DGVMs and CESM is different in winter and spring (Fig. 6). Predictability of AGR_R is higher with winter SLP (ρ_SLP is 0.61 in 40-year DJF), but predictability of SL_DGVMs and CESM is higher with spring SLP (ρ_SLP in 40-year MAM is 0.70 and 0.65, respectively).

When limiting the SLP domain to the tropics, results follow the same patterns as those at the global scale, but with better predictive skill: AGR_R shows the highest mean predictability of 0.74 for intervals of 40 years in winter (Fig. 6). SL_DGVMs shows the highest mean predictability, with ρ_SLP of 0.73 for 40 years in spring (Fig. 6), a result that is very similar to those of CESM for the same temporal length (ρ_SLP=0.70). We find that with different time intervals, tropical SLP in winter leads to higher predictability of AGR_R than global SLP fields. The spring tropical SLP only shows slightly higher predictive skill than global SLP for SL_DGVMs and CESM. AGR_R is highly influenced by winter tropical SLP, while SL_DGVMs and CESM are more sensitive to spring tropical SLP (Fig. 6). This is consistent with the results of Fig. 4 even when different time intervals are considered. This might be due to the fact that Earth system models, as CESM used here, have been found not to reliably simulate the seasonal timing of ENSO occurrence (Sheffield et al., 2013). The predictive skill of SLP in seasonal and domain differences between observation-based data and models (SL_DGVMs and CESM) could be used as an indicator to reveal their different driving mechanisms.

We evaluate the error rate of valid predictions within each time interval sliding window in Fig. 6 (i.e., the fraction of predictions in one time interval sliding window with ρ_SLP of significance P>0.05) for AGR_R, SL_DGVMs, and CESM (see Fig. A10). We find that at least 30-year-long intervals are needed for robust prediction of global C-cycle IAV. For periods shorter than 30 years, the rate of invalid predictions (in the sense given above) can be higher than 40 % for most datasets and SLP domains and seasons. It is worth noting that even when predicting AGR_R by using winter SLP in the tropical domain, the error rate can still be as low as 13 % in a 15-year interval. From the 15- to 30-year interval, the error rate of AGR_R is reduced from 0.4 to 0 in winter global and 0.13 to 0 in winter tropical. The error rate decreases to less than 0.16 in a 30-year interval except that AGR_R decreases to 0.24 in spring global, which also matches the relative low predictability of spring SLP anomalies in the period 1980–2017 with different spatial domains (see Fig. A7b). All error rates are reduced to almost zero in a 40-year interval.

4 Conclusions

The major objective of this study is to explore the relationship between SLP anomalies (as a proxy for large-scale atmospheric circulation variability) and the global C-cycle IAV. Specifically, our goals are (1) to investigate the skill of SLP to predict global C-cycle IAV using ridge regression as well as to compare with traditional teleconnection indices and (2) to establish statistical links at different spatiotemporal scales between large-scale atmospheric circulation variability and global C-cycle IAV.

First, we find that boreal winter and spring SLP anomalies allow predicting IAV of atmospheric CO₂ growth rate and IAV of the global land sink, with correlations between predicted and reference values between 0.70 and 0.84 when using winter SLP fields. This is comparable to or higher than the predictive skill of a similar model using 15 teleconnection indices as predictors. The spatial patterns of the ridge regression coefficients reveal a strong influence of large-scale atmospheric circulation variability on global C-cycle IAV, particularly in the El Niño–Southern Oscillation and western Pacific domains. Second, the comprehensive spatiotemporal sensitivity analysis indicates an increasing sensitivity of global C-cycle IAV to large-scale atmospheric circulation variability during boreal winter in the Southern Hemisphere extratropics in the recent decades. This increased sensitivity may be influenced by internal climate variability or by enhanced sensitivity of global C-cycle variability to externally forced changes but requires further research. Finally, we find that time series of at least 30 years are needed for robust predictability of global C-cycle IAV. For shorter time series, predictability is highly dependent on the particular period considered and thus largely due to statistical artifacts of random internal climate variability in the fitting process.

Overall, ridge regression using seasonal SLP fields as predictors of global C-cycle variability provides a novel and efficient data-driven approach for detecting the relationship of large-scale atmospheric circulation variability to global C-cycle variability. Compared to teleconnection indices, this approach requires no predefined spatial configurations, is more flexible to the particular domain considered, and shows equal or higher predictability of global C-cycle IAV. This method allows quantifying the contribution of atmospheric dynamical processes in driving variability in the C-cycle at global and regional scales and may further be useful for attributing observed changes to internal climate variability versus anthropogenic climate change.

Appendix A

Figure A1 ρ_SLP of AGR_R, SL_Resid, and SL_DGVMs under different DJF SLP resolution ($\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$, $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$, $\mathrm{9}{}^{\circ }\times \mathrm{9}{}^{\circ }$) by ridge regression leave-one-out cross-validation in the period of 1959–2017.

Figure A2 ρ_SLP of AGR_R, SL_Resid, and SL_DGVMs under different seasonal SLP (with different month combination) by ridge regression leave-one-out cross-validation in the period of 1959–2017. Each combination represents NDJ (November, December, and January), DJF (December, January, and February), JFM (January, February, and March), MAM (March, April, and May), and AMJ (April, May, and June).

Figure A3 Time series autocorrelations of pre-treated CO₂ time series in the period (a) 1959–2017 for AGR_R, SL_Resid, and SL_DGVMs. (b) 1980–2017 included two more inversions: SL_CAMS and SL_CarboScope. The shaded areas are the 95 % confidence interval of the calculated autocorrelation under different lags.

Figure A4 Distribution of ω_SLP with the time series of AGR_R (top row), SL_Resid (second row), SL_DGVMs (third row), SL_CAMS (fourth row), and SL_CarboScope (last row) in DJF (left column) and MAM (right column) based on SLP fields in the period 1980–2017. ω_SLP represents the mean of the n=34 run ridge regression coefficients. Note that ρ_SLP of SL_CAMS MAM has P>0.5.

Figure A5 Distribution of ω_Tele with the time series of AGR_R, SL_Resid, SL_DGVMs, SL_CAMS, and SL_CarboScope based on teleconnection indices in the period of 1980–2017.

Figure A6 Heat map of ρ_SLP with CO₂ time series over various SLP latitude domains in DJF. Each heat map contains 5×5 squares, and each square represents one domain of SLP. For example, the square 36^∘ N–72^∘ S is the domain of SLP extending from 36^∘ N to 72^∘ S. All latitudinal domains include the tropical area (18^∘ N–18^∘ S). The top right square thus represents global-scale SLP. ρ_SLP of AGR_R, SL_Resid, and SL_DGVMs in 1959–2017.

Figure A7 Heat map of ρ_SLP with CO₂ time series over various SLP latitude domains in MAM. Each heat map contains 5×5 squares, and each square represents one specific domain of SLP. For example, the square 36^∘ N–72^∘ S is the domain of SLP extending from 36^∘ N to 72^∘ S. All latitudinal domains include the tropical area (18^∘ N–18^∘ S). The top right square thus represents global-scale SLP. The time series is from (a) 1959–2017 for AGR_R, SL_Resid, and SL_DGVMs. (b) 1980–2017 included two more inversions: SL_CAMS and SL_CarboScope.

Figure A8 Heat map of ρ_SLP with CO₂ time series over various SLP latitude domains in DJF+MAM. Each heat map contains 5×5 squares, and each square represents one specific domain of SLP. For example, the square 36^∘ N–72^∘ S is the domain of SLP extending from 36^∘ N to 72^∘ S. All latitudinal domains include the tropical area (18^∘ N–18^∘ S). The top right square thus represents global-scale SLP. The time series is from (a) 1959–2017 for AGR_R, SL_Resid, and SL_DGVMs. (b) 1980–2017 included two more inversions: SL_CAMS and SL_CarboScope.

Figure A9 Heat map of ρ_SLP of AGR_R with DJF SLP over various latitude domains. A 30-year sliding window in the period of 1959–2017 with a 1-year step is created. The starting and end year of each interval is labeled on the top of each heat map.

Figure A10 Error rate of the sliding window in 15-, 20-, 30-, and 40-year intervals. For each sliding window, the error rate is calculated by the number of invalid predictions (with significance P>0.05 in ρ_SLP of CO₂ time series) divided by the number of total predictions. With SLP in DJF global, MAM global, DJF tropical, and MAM tropical as predictors, the error rates of ρ_SLP of (a) AGR_R, (b) SL_DGVMs, and (c) CESM are plotted.

Figure A11 The spatial distribution of Pearson correlations between global SLP-predicted AGR_R (one time series) and global pixel-based land temperature and precipitation anomalies (both from the CRU TS4.05 Harris et al., 2020, monthly dataset, aggregated to annual mean temperature and annual sum precipitation, detrended by LOWESS) in the period 1980–2017. The top panel shows the spatial distribution of Pearson correlations between pixel-based land temperature anomalies and global SLP-predicted AGR_R for DJF (left) and MAM (right), and the bottom panel shows Pearson correlations of land precipitation anomalies.

Figure A12 The spatial distribution of Pearson correlations between global DJF–MAM SLP-predicted AGR_R (one time series) with pixel-based annual sum NBP variation (LOWESS detrended) from the atmospheric inversion CarboScope s76 (Bastos et al., 2020; Chevallier et al., 2005; Rödenbeck et al., 2003) (upper panel) and CAMS (Bastos et al., 2020; Chevallier et al., 2005; Rödenbeck et al., 2003) (lower panel) in the period 1980–2017.

Figure A13 Comparison of linear regression based on SLP principal components to ridge regression based on SLP fields in the period of (a) 1959–2017 and (b) 1980–2017. The label of the y axis shows the Pearson correlation between the predicted and original CO₂ time series. The solid lines represent the results by using the linear regression based on a different number of SLP principal components as predictors. The dashed lines represent the results by using ridge regression. Differently colored lines represent different SLP spatial domains selected.

Figure A14 The explained variance of SLP principal components extracted from global DJF SLP fields in the period of 1959–2017 (a) and 1980–2017 (b).