Significance

The future of the terrestrial carbon (C) sink has tremendous consequences for society and the rate of climate change, but is highly uncertain. The sensitivity of interannual variability in the C sink to climate drivers can help elucidate the mechanisms driving the C sink. Here, we test the statistical strength of major climate drivers of the C sink and find that nighttime tropical temperatures are most strongly associated with the global C sink from 1959–2010, likely acting through their effect on respiration. The temperature-mediated sensitivity of the global C sink to respiration highlights that tropical C stores may be vulnerable to robust projected increases in tropical nighttime temperatures.

Abstract

The terrestrial biosphere is currently a strong carbon (C) sink but may switch to a source in the 21st century as climate-driven losses exceed CO2-driven C gains, thereby accelerating global warming. Although it has long been recognized that tropical climate plays a critical role in regulating interannual climate variability, the causal link between changes in temperature and precipitation and terrestrial processes remains uncertain. Here, we combine atmospheric mass balance, remote sensing-modeled datasets of vegetation C uptake, and climate datasets to characterize the temporal variability of the terrestrial C sink and determine the dominant climate drivers of this variability. We show that the interannual variability of global land C sink has grown by 50–100% over the past 50 y. We further find that interannual land C sink variability is most strongly linked to tropical nighttime warming, likely through respiration. This apparent sensitivity of respiration to nighttime temperatures, which are projected to increase faster than global average temperatures, suggests that C stored in tropical forests may be vulnerable to future warming.

Terrestrial ecosystems have been a substantial net sink of anthropogenic carbon (C) emissions since the 1960s (1⇓⇓–4), but the terrestrial C sink could switch to a C source in the 21st century, resulting in a positive C cycle-climate feedback that would accelerate global surface warming with potentially major consequences for the biosphere (5⇓–7). The interannual variability of the terrestrial C sink can help constrain our understanding of C/climate feedbacks and identify regions and mechanisms of the terrestrial C cycle that are most sensitive to climate parameters, shedding light on the future of the sink and its possible transition to a source (8). Currently, several major drivers have been shown to be correlated with the interannual variability of the terrestrial C sink, including (i) tropical temperature, which is tightly coupled to interannual variability in the atmospheric growth rate (AGR) of CO2 (8, 9); (ii) tropical drought stress, including major droughts in the Amazon (10⇓–12), which has been suggested to underlie increasing sensitivity of the AGR to tropical temperature over the period from 1959–2010 (13); (iii) temperature and precipitation variability in semiarid regions (14, 15); and (iv) average minimum daily (hereafter “nighttime”) temperatures, which studies of several local field sites in the tropics have found play a major role in interannual productivity (16⇓–18).

Determining the mechanism underlying the interannual variability of the terrestrial C sink, including the relative roles of precipitation vs. temperature stress and their effects on gross primary productivity (GPP) vs. total respiration (both autotrophic and heterotrophic; R), is critical to predict the sink’s future and to improve Earth system models. Here, we quantify changes in the interannual variability of the terrestrial C sink over the past half-century and then statistically evaluate four hypotheses that the variability of the terrestrial sink is most strongly influenced by (i) tropical mean temperature, (ii) tropical precipitation, (iii) precipitation and temperature in semiarid regions, and (iv) nighttime tropical temperatures. We combine multiple simulations from an atmospheric mass balance of the land C sink [net ecosystem exchange (NEE)] from 1959 to 2010, remote sensing-modeled datasets of vegetation greenness and GPP from 1982 to 2010, and global gridded climate datasets to constrain globally the fundamental equation NEE = GPP − R and the relative sensitivities of each component to temperature and precipitation. We draw on a combination of model selection and partial correlation analysis to provide relative likelihood estimates of each driver and to account for covariation between predictor variables (e.g., tropical mean temperature vs. nighttime temperature).

Results and Discussion

We calculated a global mass balance (Fig. 1 A–E), such that changes in the atmosphere (AGR) must be the result of emissions from fossil fuels (EF) and land use (EL) and uptake by land (NEE) and ocean (SO) from 1959 to 2010. We first estimated time-varying uncertainties for EF and EL, ocean uptake (SO), and the atmospheric C increase. Each of these flux estimates was derived from multiple datasets. The errors and estimates were then used in a Monte Carlo simulation to generate 4,500 alternative trajectories of NEE. Then, using a moving window of 20 y and multiple alternate methods (SI Appendix, Fig. S1), we found that the variance in the AGR has increased over the period from 1959–2010 (99% of linear slopes positive; F test between 1959–1984 and 1985–2010; F = 0.5, P = 0.001) (Fig. 1 F and K) and that neither the emissions fluxes (Fig. 1 L and M) nor the SO flux (Fig. 1N) can account for this increase in variance. We caution, however, that the true variability of fossil fuel, land use, and ocean fluxes is likely underestimated due to data and model limitations.

Increasing variance in atmospheric growth rate of CO2 is driven by increasing variance of the terrestrial C sink from 1959 to 2010. Annual components of the global C budget (PgC⋅y−1; shading is ±1 SD): AGR of CO2 (A), anthropogenic EF (B), anthropogenic EL (C), SO (D), and terrestrial C uptake (E). (F–J) Twenty-year moving windows of variance of each component. (K–O) Density of linear slopes of the variance for each flux across all Monte Carlo simulations, with the percentage of slopes that were positive (as shown by shading).

Thus, we can attribute this increasing variability with high confidence to C uptake by terrestrial ecosystems (Fig. 1O; 96% of slopes positive; F test between 1959–1984 and 1985–2010; F = 0.4, P = 0.009), whose interannual variability is widely considered to be driven by sensitivities to climate, such as El Niño-Southern Oscillation (ENSO)–mediated shifts in precipitation and temperature patterns (9, 19, 20). Although there is general consensus that interannual variability in the global C cycle is dominated by terrestrial processes (20, 21), we demonstrate here that this variability has been increasing over the past five decades. Recently reported increasing sensitivity of NEE to tropical temperatures (13) may partially explain the increased variability, although we explore this hypothesis in more detail below. Follow-up analyses revealed that this increasing variability was not due to increasing mean sink (SI Appendix, Fig. S2), increasing coverage of the CO2 observing network (SI Appendix, Fig. S3), or volcanic eruptions in the 1990s (SI Appendix, Fig. S4). However, variance in NEE peaked during the 1990s due, in part, to the greater frequency and intensity of El Niño events in the 1980s and 1990s, which have subsequently waned in the 2000s (Fig. 1J and SI Appendix, Fig. S4). Nevertheless, an increase in NEE variance is still observed when ENSO is removed, albeit much smaller (SI Appendix, Figs. S1 and S4).

Examining the climate drivers of interannual variability using both model selection and partial correlation analysis, we found that the interannual variability of terrestrial NEE was most sensitive (SI Appendix, Figs. S5 and S6) to nighttime temperature anomalies averaged over the terrestrial tropics (Tmin-Tropics) (R2 = 0.31, P < 0.0001; Fig. 2). Comparing models’ relative likelihood with Akaike’s information criterion (AIC), tropical nighttime temperatures were fourfold to fivefold more likely to explain interannual NEE variability than were tropical mean (0.25 relative likelihood) or tropical maximum (0.22 relative likelihood) temperature, and more than 100-fold more likely than the other climate drivers tested (SI Appendix, Fig. S6). When using both forward and backward stepwise model selection to allow for models with up to 12 independent variables, the model with only Tmin-Tropics and lagged semiarid precipitation was selected as the most parsimonious (R2 = 0.47, P < 0.0001), with coefficients indicating a much stronger role of Tmin-Tropics than of semiarid precipitation (SI Appendix, Table S1). Partial correlation analyses produced similar results to model selection (SI Appendix, Fig. S5), and all results were robust regardless of the SO estimates used (SI Appendix, Table S2). Moreover, the sensitivity of the terrestrial sink to tropical nighttime temperature has increased over time (SI Appendix, Fig. S5), indicating that nighttime temperatures are most likely driving the previously reported correlations between the AGR and tropical temperature (13). Spatially averaged tropical precipitation was not significantly correlated with NEE using either Pearson correlation coefficients (r = 0.19, P = 0.18) or partial correlation coefficients accounting for other correlated variables (r = 0.11, P = 0.76). This lack of correlation indicates that tropical moisture availability alone does not explain interannual variation in NEE, although temperature and precipitation interactions may still be important.

Respiration appears to mediate the impact of tropical nighttime temperature on the interannual variability of NEE. (A) Partial regression of NEE vs. tropical maximum (Tmax) temperatures, where Tmin temperatures have been removed. (B) Partial regression of NEE vs. Tmin temperatures, where Tmax temperatures have been removed. (C) Anomaly of detrended global GPP compared with the anomaly of detrended Tmin temperatures from 1982 to 2010. (D) Anomaly of detrended global respiration (R) compared with the anomaly of detrended Tmin temperatures from 1982 to 2010. The dashed line indicates the best fit of ordinary least squares regression.

Consistent with recent findings, our results show that the variance in GPP has been increasing in arid ecosystems (15) (SI Appendix, Fig. S7). However, we find decreasing variance in GPP across all other climate zones, resulting in decreasing variance in GPP at the global scale (SI Appendix, Fig. S7). Decreasing variance in GPP provides ancillary evidence that R is driving increasing variance in NEE. These findings are in contrast to recent work that suggested increasing variance of GPP, particularly in semiarid ecosystems, is the main driver of observed variance in NEE (15). Differences in the magnitude and direction of variance in satellite vegetation indices may explain this discrepancy and should be the focus of future research (SI Appendix). The GPP remote-sensing algorithm has known limitations (22). Although it has reasonable skill in capturing interannual variability in semiarid and evergreen broadleaf tropical forest biomes (23) central to our analyses here, it has difficulties in capturing interannual variability in many biomes and the required climatic data inputs could drive observed trends. Thus, we tested the robustness of our results in two ways. First, we verified that all trends in GPP were apparent in the raw satellite data [fraction of photosynthetically active radiation (FPAR); SI Appendix, Fig. S7]. Second, we determined that the trends were robust to the full uncertainty in climate inputs to GPP (SI Appendix, Fig. S7). Model selection using the AIC revealed that tropical nighttime temperatures were the only significant variable in explaining interannual patterns in R (SI Appendix, Fig. S8 and Table S3). Thus, the correlation between respiration and nighttime tropical temperatures was found to be stronger than mean (relative likelihood of 0.45 compared with the nighttime temperature model) or maximum (relative likelihood of 0.39) tropical temperatures, although these differences are less robust than the NEE model differentiation (black lines in SI Appendix, Fig. S6).

Although high tropical daytime temperatures may sometimes reduce photosynthesis rates, the greater temperature sensitivity of the global C sink is expected to act through respiration (24). Furthermore, nighttime warming patterns will differentially affect respiration more than photosynthesis, thereby reducing C uptake (25, 26). Nighttime temperatures have risen faster than daytime temperatures globally and in the tropics (SI Appendix, Fig. S9). Long-term measurements from old-growth tropical forests in Costa Rica show that nighttime temperatures are a major factor determining interannual differences in aboveground productivity, explaining 73% of the interannual variance in tree growth when combined with dry season vapor pressure deficit (16, 17). Previous analyses of nighttime warming at a global scale have only considered the sensitivity of C uptake (i.e., GPP) (26), whereas we consider both C uptake and loss here. Changing disturbance regimes, particularly fire, may also play a role. We found that fire emissions are only weakly associated with interannual NEE variability (SI Appendix, Figs. S10 and S11), but limitations in long-term disturbance and fire data prevent a robust analysis of how they may be affecting variability in the terrestrial C cycle.

In conclusion, we have shown that the variance in the AGR has increased from 1959 to 2010 due to increased interannual fluctuations in net terrestrial C uptake. Tropical nighttime temperatures are the strongest statistical driver of this observed increase in interannual variability of net terrestrial C uptake and total respiration. The central role of elevated nighttime temperatures confirms the global importance of a mechanism previously observed only in a local field study (16, 17) and may reduce the terrestrial C sink in light of nighttime warming projections (25, 27) because C gains in tropical ecosystems may be offset by greater respiratory losses as a result of nighttime warming. Respiration-driven losses in forest C constitute one major scenario through which the terrestrial C sink could switch to a source. Because tropical forests account for 33% of the annual primary productivity of the terrestrial biosphere (28) and are currently a significant C sink (excluding ELs) (4), the C balance of tropical forests is of heightened concern. Furthermore, because C residence times in tropical forests are relatively low (13), these forests could also respond relatively rapidly to directional changes, including warming nighttime temperatures.

Although considerable research has focused on the temperature sensitivity of the high-latitude C balance to warming (29), our results highlight that tropical C dynamics are also likely sensitive to warming. If respiratory losses of C eclipse current C gains due to CO2 fertilization and secondary regrowth from land use (30), forest loss could undermine efforts to offset emissions through reducing emissions from deforestation and degradation policies and lead to a strong C/climate feedback. Our results suggest that the previously identified emergent relationship between temperature and the AGR is most likely due to respiration losses; thus, respiratory processes should be prioritized in the improvement of Earth system models and their future climate predictions. In addition, establishment of long time series CO2 measurements at sufficient spatial density in tropical regions (4) and additional constraints on GPP measurements through solar-induced chlorophyll fluorescence (31) would greatly reduce uncertainties in interannual variability of regional and global C budgets. Increased understanding of the drivers of mean and interannual variability in terrestrial C uptake can ultimately help identify the causes and magnitudes of C cycle feedbacks to anthropogenic climate change in the 21st century.

Methods

C Cycle Inversion.

We constructed an annual global C budget via standard methods (1, 32) of mass balance where anthropogenic EF and EL emission fluxes to the atmosphere are in balance with the atmospheric growth rate AGR and uptake of sinks in the ocean (SO) and land (NEE):EF+EL=AGR+SO+NEE.[1]

As the least well-constrained component of the global C budget, the land C sink is calculated as a residual of the other terms. The interannual variability in this inversion-calculated residual lines up relatively well with the interannual variability from process-based dynamic global vegetation models (33).

The AGR was determined using measurements of atmospheric CO2 concentrations from the Global Greenhouse Gas Network (data source: www.esrl.noaa.gov/gmd/ccgg/). Following standard C cycle inversion methods, we calculated the annual AGR as the average of a year’s December and January concentrations minus the average December and January concentrations from the previous year, converting concentrations to mass [petagrams of C (PgC)] using the molar mass conversion factor 2.124 PgC⋅ppm−1. This standard method assumes that CO2 in the atmosphere is relatively well mixed on time scales longer than 1 y, although we considered a temporal autocorrelation of error in our propagation of uncertainty in calculating the AGR (discussed below). December/January concentrations are commonly used because they accurately capture annual atmospheric growth rate (r = 0.93, P < 0.0001; SI Appendix, Fig. S12), whereas other months can be strongly influenced by widely documented changes in the seasonality of atmospheric CO2 concentrations (2), which are likely a signal of agricultural productivity and not necessarily the terrestrial C sink (34, 35). Because measurement accuracy and precision are quite high for CO2 concentrations (36), the majority of the uncertainty in calculating the annual AGR derives from spatial heterogeneities in the uneven sampling network (37).

To assess uncertainty in the AGR, we performed 100 bootstrap simulations of the AGR from the global sampling network. We resampled with replacement of 40 sites from the marine boundary layer to create each bootstrap simulation, with the only constraint being that at least one site must come from each of the North Atlantic, North Pacific, Arctic, Antarctic, and tropics regions. To account for the documented negative interannual autocorrelation in errors in the AGR (37), we applied autocorrelated noise equivalent to the autocorrelated uncertainty (i.e., same SD, same autocorrelation coefficients) observed across the 100 marine boundary layer sampling sites from 1980 to 2010 to the bootstrapped simulations from 1959 to 1979. This method yielded 100 bootstrapped simulations that should reasonably account for spatial and temporal uncertainty in the calculation of the AGR from the global observing network over time. To ensure that the increasing coverage of the observing network was not driving the rising interannual variability, we verified that the trends in variability were also evident in the AGR of the two longest measurement sites: Mauna Loa (SI Appendix, Fig. S3) and the South Pole (not shown).

EF were calculated from three independent global inventory datasets: British Petroleum, the Carbon Dioxide Information and Analysis Center, and the Emission Database for Global Atmospheric Research [data available at the Global Carbon Project (GCP): www.globalcarbonproject.org/carbonbudget/]. These inventories use energy consumption statistics to calculate emissions to the atmosphere from fossil fuel combustion and cement production. Although errors in these inventories are generally thought to be relatively small, EF constitute the largest flux in the C budget; thus, accounting for potential uncertainty and errors is critical to evaluating and constraining the global C cycle (38). The datasets differ in their treatment of international transport, gas flaring, cement production, and other areas; accounting practices of different institutions and countries provide an additional component of uncertainty (39).

To account for this uncertainty, we performed 100 bootstrap simulations per inventory. Although the true value of spatial errors across countries is largely unknown, there is general consensus that errors are likely smaller in Organization for Economic Cooperation and Development (OECD) nations than in non-OECD nations (39). We account for spatial uncertainty and errors in reporting practices by assigning a 5% potential error in emissions from OECD nations and a 10% potential error in emissions from non-OECD nations. Similar to the AGR, errors in EF are expected to exhibit positive temporal autocorrelation (i.e., erroneously high reported emissions will remain erroneously high until that error is spotted and then retroactively corrected for all previous years) (39). Thus, we added temporally autocorrelated random noise with an AR1 coefficient of 0.95 to our bootstrap simulations. This approach gives a persistence of autocorrelated errors of around 20 y in simulations.

We used three estimates of EL: a bookkeeping method based on deforestation patterns and biomass data (40) and two vegetation model-based estimates [Land surface Processes and eXchanges model from the University of Bern (LPX-Bern) (41) and Integrated Science Assessment Model (ISAM) (42)] (data available at www.globalcarbonproject.org/carbonbudget/). ELs contain perhaps the largest uncertainty of any estimate in our analysis because they include a large number of poorly constrained processes, different combinations of which are included in each estimate. The standard bookkeeping method is considered a benchmark product and includes the flux to the atmosphere of clearing of primary forest and conversion to agriculture based on deforestation statistics, satellite estimates of deforestation, and biomass maps of deforested areas (40). The LPX-Bern model-based estimate is based on historical land-use maps and also includes the abandonment of land and regrowth of secondary forest (41). The ISAM model-based estimate includes historical land use and environmental responses such as nitrogen cycling, which influences secondary growth and uptake (42). Although other model-based datasets of EL exist, they rely on similar boundary conditions of historical land-use change maps; thus, we believe these three datasets likely encompass the entire range and uncertainty of CO2 emissions estimates for scenarios of land-use and land-cover change.

To account for uncertainty in EL, we performed 100 bootstrap simulations for each of the EL datasets, using uncertainty provided by the GCP. Similar to EF, positive autocorrelation of error is expected in EL because the input data used to derive some components of EL (e.g., deforestation statistics from the United Nations Food and Agriculture Organization) are released every 5 y. Thus, we used the same autocorrelated noise approach as in EF above to account for this known autocorrelation in the bootstrap simulations.

Ocean C uptake (SO) was determined using two independent methods with differing degrees of observational constraints. First, we used a data assimilation method where annual SO in an oceanic pulse response model was constrained with the oceanic inventory of anthropogenic CO2 from the 1990s (43), historical fluxes from 1980 to 2009 estimated from surface pCO2 measurements, and the historical atmospheric CO2 concentration (44). The method relies on Markov chain Monte Carlo sampling of the model parameter space, and we took 1,000 samples from the Markov chain as samples of the SO under uncertainty (45). The second method used the multimodel average of six process-based atmosphere-ocean general circulation models (AOGCMs) contained in the GCP data archive (data available at www.globalcarbonproject.org/carbonbudget/), which have been constrained to match the 1990–1999 observation-driven uptake estimates of 2.2 PgC⋅y−1 (33). The AOGCMs we used were Nucleus for European Modelling of the Ocean (NEMO-PlankTOMS) (46), Laboratoire des Sciences du Climat et de l'Environnement (LSCE) (47), Community Climate System Model - Biogeochemical Elemental Cycling (CCSM-BEC) (48), Model for Interdisciplinary Research On Climate - HAMburg Ocean Carbon Cycle (MICOM-HAMOCC) (49), Max Plank Institute Ocean Model - HAMburg Ocean Carbon Cycle (MPIOM-HAMOCC) (50), and Princeton (51).

To account for uncertainty in SO in the GCP models, we used the published 1σ uncertainty (33) of 0.5 PgC⋅y−1 to construct 1,000 bootstrap simulations around the model mean. As with the emissions fluxes, autocorrelation of uncertainty is expected in ocean fluxes; we thus followed the same method as above, with an autocorrelation coefficient of 0.9. The data assimilation technique used for constraining the pulse response model generates estimates of the autocorrelation coefficient and total error from the residuals with the ocean flux history. We constructed bootstrapped noise for the data-constrained SO with an estimated autocorrelation coefficient of 0.29 and total error of 0.51 PgC⋅y−1, following the assumption that the residuals from the ocean flux time series represented short-term variability not resolved by the pulse response model.

To calculate the residual land C sink with a full accounting for uncertainty through the global C budget, we used a Monte Carlo approach to construct a set of emissions simulations based on a factorial combination of the three EF and three EL datasets (3 × 3 = 9 dataset combinations × 500 runs each = 4,500 simulations). We then subtracted the 100 simulations of the AGR to generate the simulations of the net C sink (ocean + land). Finally, to constrain the residual land sink, we constructed two datasets of 4,500 simulations each by subtracting the estimated ocean C uptake (random subsampling with replacement of the 1,000 runs) from each SO dataset. For all subsequent analyses, we analyzed these two land C sink datasets separately, although all analyses were robust to dataset choice (discussed below).

Climate, Volcanic, and Solar Radiation Data.

We downloaded global average and gridded climate datasets to assess the sensitivity of the land C sink and its interannual variability to climate variables between 1959 and 2010. We acquired gridded mean annual temperature data from the Hadley Center’s Climatic Research Unit (HadCRUT4) (52), the National Aeronautics and Space Administration’s (NASA) Goddard Institute for Space Studies (53), and the National Oceanic and Atmospheric Administration’s (NOAA) National Climatic Data Center (54). We acquired gridded annual precipitation data from the NOAA’s Global Precipitation Climatology Center (55). We acquired gridded monthly mean maximum and monthly mean minimum (nighttime) temperatures from the Hadley Center’s Climatic Research Unit Time Series (CRU TS 3.21) (56). From each dataset, we extracted and calculated the annual global average and zonal averages for the tropics (25°S to 25°N), Northern extratropics (25.5°N to 90°N), and boxes covering the semiarid regions presented by Poulter et al. (14) (three boxes: South America: 50°S 80°W to 20°S 40°W, Southern Africa: 35°S 10°E to 10°S 40°E, and Australia: 40°S 115°E to 10°S 155°E). We calculated the 1-y lags from the semiarid regions as well, because lags may be important in these areas (14). We drew volcanic aerosol emissions from NASA’s Goddard Institute for Space Studies dataset of optical aerosol thickness at 550 nm (57). Global average solar radiation data were 10.7 cm of solar flux data, provided by the National Research Council of Canada. Details about GPP and fire emissions analyses are provided in SI Appendix.

Statistical Analyses.

For all analyses except calculation of the coefficient of variation of the land C sink, all variables were first detrended with a linear trend. To analyze the change in variability of each component of the C cycle, we used a 20-y moving window across all Monte Carlo simulations of each flux (discussed above) and calculated the variance of each simulation in each window. We plotted the trajectory of each simulation’s variance trend and calculated the percentage of these linear slopes that were positive (as assessed with ordinary least squares regression) and the statistical significance of the mean simulation with ordinary least squares regression. For the land C sink (NEE), we performed this analysis on both datasets (one estimated using the GCP biogeochemistry ocean models and the other estimated using the data assimilation-constrained SO). The increase in NEE variance was robust across both datasets (SI Appendix, Fig. S13). We repeated this analysis with differing window lengths ranging from 15 to 25 y and found the increase in NEE variability to be robust to window length (95–97% of slopes positive in all window lengths). We further validated the robustness of the increase in variance by using a Fisher test that compares the variance of two groups. We performed this test on the first and last 20 y and 25 y of the record and found that the variances were significant (20-y windows: F = 0.4, P = 0.03; 25-y windows: F = 0.4, P = 0.02). We performed the same moving window analyses on GPP and FPAR estimated from Moderate Resolution Imaging Spectroradiometer (MODIS), using 5- and 10-y moving windows (SI Appendix, Fig. S7). All NEE, GPP, and FPAR variables met the assumptions of normality.

To test the robustness of this increase in variance in NEE to method and volcanic signals, we used wavelet transformations on the mean detrended NEE across the Monte Carlo simulations and the residuals of mean detrended NEE after regressing against volcanic aerosol in the current year and a 1-y lag. Similar to a fast Fourier transform, a wavelet transform decomposes a time series into intensity at each frequency, centered around each year (58). We used a Morlet wavelet and a Mexican hat wavelet, which provide better resolution in the frequency and time domains, respectively (58). We found that the variability of NEE increased in all cases, shown as hotter colors (higher intensity) over time in the wavelet power spectra (SI Appendix, Fig. S1).

To test whether increasing variability was a product of increasing mean, we repeated the moving window analyses detailed above but, instead, did not detrend the time series of NEE and calculated the coefficient of variation (defined as SD/mean) in each window. We found that the coefficient of variation also increased over time (SI Appendix, Fig. S2), indicating that the SD (and thus variance) increased faster than the mean. We further found via partial regression that the variance of NEE increased when removing the volcanic signal and ENSO signals (SI Appendix, Fig. S4), indicating that the Pinatubo eruption and strong 1997/1998 ENSO event were not the sole drivers of the increased variance.

To assess the sensitivity of interannual variation in NEE to climate, we used two independent techniques to determine the most important variables among a suite of highly correlated explanatory climate variables. We first used model selection based on the AIC, considered a standard technique for determining the most parsimonious model for a statistical relationship (59⇓–61). All variables were first detrended and converted to Z-scores, with the volcanic signal removed from NEE via partial regression. These climate variables met the assumptions of normality. For each of the 4,500 simulations of NEE, we calculated the AIC of single linear regressions of NEE vs. each of 12 climate variables. Climate variables included in stepwise model selection were as follows: Tmin-Tropics (nighttime); average tropical temperatures (Tave-Tropics); maximum tropical temperatures (Tmax-Tropics); global average temperature (Tave-Globe); average temperature of latitudes >30°N (Tave-North); maximum and minimum temperatures at those latitudes (Tmax-North and Tmin-North, respectively); precipitation for the globe (Precip-Globe), tropical regions (Precip-Tropics), and northern extratropics (Precip-North); semiarid regions’ current year’s temperature (SaT) and 1-y lag (SaT1); and semiarid regions’ current year’s precipitation (SaP) and 1-y lag (SaP1). Per standard model selection methods, we then calculated the relative likelihood of each model compared with the best model (lowest AIC) as Relative LikelihoodModelI = exp[(AICMin − AICModelI)/2].

To determine the most parsimonious model, we then performed stepwise model selection both “forward” (starting with the null model and adding explanatory variables) and “backward” (starting with the complete model and removing explanatory variables). To remove covarying variables, we used a standard model selection technique that involved creating a correlation matrix between explanatory variables, identifying variables correlated at r > 0.5, and comparing each independently with NEE. The variable that explained more of the variance in NEE was retained. All stepwise model selection routines culminated with mean monthly minimum temperatures in the tropics (Tmin-Tropics) and 1-y lagged semiarid precipitation as the most parsimonious model. We performed model selection with mean NEE calculated from both SO datasets (GCP models and data assimilation uptake) and found that the final model was the same in both cases. We tested all models to ensure that explanatory variables and residuals were normally distributed and met the assumptions of linear regression.

We further tested the robustness of the stepwise model selection using partial correlation analysis. Partial correlation analysis is used to determine the correlations between two given variables (in our case, terrestrial NEE and a given climate variable) while accounting for other correlated variables. It has been used widely to examine the sensitivity of C cycle fluxes to climate parameters (9, 26, 28). All variables were first detrended and converted to Z-scores. We performed partial regressions between each pair of variables while accounting for all other variables. We performed partial correlations with NEE calculated from both SO datasets (i.e., GCP models, data assimilation method) and found that the partial correlations exhibited the same qualitative pattern; thus, we only display the GCP model uptake-based land NEE in SI Appendix, Fig. S5. We further examined the strength of the partial correlation of NEE and tropical nighttime temperatures over time, accounting for mean monthly maximum tropical temperature, solar radiation, and tropical precipitation in one model and global mean annual temperature, solar radiation, and tropical precipitation in a second model, in moving 20-y windows. We performed this analysis on the 4,500-member ensemble of NEE estimated from the GCP models of SO.

All statistics were performed in the R computing environment (58). Wavelet analysis was performed using the “biwavelet” package (62). Climate variables in NetCDF were imported and calculated using the “RnetCDF” package (63). Stepwise model selection was performed using the “MASS” package (64). Partial correlations were performed using the “ggm” package (65).

Acknowledgments

We thank Pekka Kauppi and the Finnish Society of Sciences and Letters, and the Carbon Mitigation Initiative at Princeton University for sponsoring the workshop that led to this paper. W.R.L.A. was supported, in part, by US National Science Foundation (NSF) MacroSystems Biology Grant Division of Environmental Biology EF-1340270, NSF RAPID Grant DEB-1249256, NSF Early Concept Grants for Exploratory Research (EAGER) Grant 1550932, and a National Oceanic and Atmospheric Administration (NOAA) Climate and Global Change postdoctoral fellowship, administered by the University Corporation of Atmospheric Research.

Footnotes

↵1To whom correspondence may be addressed. Email: anderegg{at}princeton.edu or pacala{at}princeton.edu.

↵2Present address: Institute on the Environment, University of Minnesota, St. Paul, MN 55108.

(1994) Evidence for interannual variability of the carbon cycle from the National Oceanic and Atmospheric Administration/Climate Monitoring and Diagnostics Laboratory global air sampling network. J Geophys Res Atmos99(D11):22831–22855

Blood-sucking sand flies from disparate global regions have a predilection for feeding on the marijuana plant (Cannabis sativa), and the findings hint at a potential avenue for controlling sand flies, which can transmit leishmaniasis.