The dominant role of semi-arid ecosystems in the trend and variability of the land CO2 sink

See allHide authors and affiliations

Science  22 May 2015:
Vol. 348, Issue 6237, pp. 895-899
DOI: 10.1126/science.aaa1668

The difference is found at the margins

The terrestrial biosphere absorbs about a quarter of all anthropogenic carbon dioxide emissions, but the amount that they take up varies from year to year. Why? Combining models and observations, Ahlström et al. found that marginal ecosystems—semiarid savannas and low-latitude shrublands—are responsible for most of the variability. Biological productivity in these semiarid regions is water-limited and strongly associated with variations in precipitation, unlike wetter tropical areas. Understanding carbon uptake by these marginal lands may help to improve predictions of variations in the global carbon cycle.

Science, this issue p. 895


The growth rate of atmospheric carbon dioxide (CO2) concentrations since industrialization is characterized by large interannual variability, mostly resulting from variability in CO2 uptake by terrestrial ecosystems (typically termed carbon sink). However, the contributions of regional ecosystems to that variability are not well known. Using an ensemble of ecosystem and land-surface models and an empirical observation-based product of global gross primary production, we show that the mean sink, trend, and interannual variability in CO2 uptake by terrestrial ecosystems are dominated by distinct biogeographic regions. Whereas the mean sink is dominated by highly productive lands (mainly tropical forests), the trend and interannual variability of the sink are dominated by semi-arid ecosystems whose carbon balance is strongly associated with circulation-driven variations in both precipitation and temperature.

Since the 1960s, terrestrial ecosystems have acted as a substantial sink for atmospheric CO2, sequestering about one-quarter of anthropogenic emissions in an average year (1). This ecosystem service, which helps mitigate climate change by reducing the rate of increase of atmospheric greenhouse gases, is due to an imbalance between the uptake of CO2 through gross primary production (GPP, the aggregate photosynthesis of plants) and the release of carbon to the atmosphere by ecosystem respiration (Reco) and other losses, including wildfires (Cfire). The net carbon flux (net biome production, NBP = GPP – Reco – Cfire) results from the small imbalance between the much larger uptake and release fluxes. Consequently, small fractional variations in either of these fluxes can cause substantial absolute variations in net carbon exchange with the atmosphere. These variations account almost entirely for year-to-year variations around the overall trend in atmospheric concentrations of CO2 (2, 3).

Modeling studies suggest a large uncertainty of the future magnitude and sign of the carbon sink provided by terrestrial ecosystems (48). Robust projections are crucial to assessments of future atmospheric CO2 burdens and associated climate change, and are therefore central to the effectiveness of future mitigation policies. Reducing the uncertainty of these projections requires better knowledge of the regions and processes governing the present sink and its variations. Inventories suggest that since the beginning of industrialization, the majority of carbon sequestered by the terrestrial biosphere has accumulated in forest ecosystems of the tropics and temperate zones (9). However, the relative contributions of ecosystems of different, climatically distinct, regions to variations in the land sink on interannual to multidecadal time scales are not well characterized. Here, we investigated relative regional contributions to the mean sink, to its trend over recent decades, and to the interannual variability (IAV) around the trend.

We used LPJ-GUESS (1012), a biogeochemical dynamic global vegetation model, to simulate the geographic pattern and time course of NBP. LPJ-GUESS explicitly accounts for the dependency of plant production and downstream ecosystem processes on the demography (size structure) and composition of simulated vegetation. We forced the model with historical climate (13) and CO2 concentrations, accounting for emissions from land use change and carbon uptake due to regrowth after agricultural abandonment (14). We compared the results to an ensemble of nine ecosystem and land surface model simulations from the TRENDY model intercomparison project (12, 15) (hereinafter TRENDY models; table S1). The TRENDY ensemble is similarly based on historical climate and CO2 but uses a static 1860 land use mask.

Global NBP, as simulated by LPJ-GUESS, shows strong agreement (r2 = 0.62) with the Global Carbon Project (GCP) estimate of the net land CO2 flux—an independent, bookkeeping-based estimate derived as the residual of emissions, atmospheric growth, and ocean uptake of CO2 (1) (Fig. 1A). TRENDY models do not account for land use change. Relative to the GCP land flux estimate, they consequently predict a higher average NBP but similar interannual variation. Moreover, the offset between the TRENDY model ensemble mean and the GCP land flux estimate is comparable to the GCP estimate of mean land use change emission flux for the period 1982–2011 (fLUC).

Fig. 1 Global and regional NBP mean, trend, and variations (1982–2011).

(A) Global NBP from LPJ-GUESS (red line) and GCP land flux time series (black line) with ±0.8 Pg C uncertainty range (shaded gray area). TRENDY models mean (blue line) and first and third quartiles (shaded blue area) are plotted on a separate axis with a time-invariant offset corresponding to the time period average GCP fLUC estimate (1.2 P Pg C year-1). (B) Tropical forest NBP. LPJ-GUESS (red line) includes emissions from land use change. TRENDY models average (blue line) and first and third quartiles of the ensemble (shaded blue area) do not include emissions from land use change. (C) NBP of semi-arid ecosystems from LPJ-GUESS (including land use change emissions) and TRENDY models (excluding land use change emissions); colors and shading as in (B). (D) Contribution of land cover classes to global mean NBP (1982–2011) (mean NBP of land cover class as a proportion of mean global NBP). Horizontal lines in box plots show, from top to bottom, 95th, 75th, 50th, 25th, and 5th percentiles. (E) Contribution of land cover classes to global NBP trend (land cover class NBP trend as a proportion of global NBP trend). (F) Contribution of land cover classes to global NBP IAV (Eq. 1).

We divided the global land area into six land cover classes, following the MODIS MCD12C1 land cover classification (12, 16): tropical forests (Fig. 1B), extratropical forest, grasslands and croplands (here combined), semi-arid ecosystems (Fig. 1C), tundra and arctic shrub lands, and sparsely vegetated lands (areas classified as barren) (figs. S1 and S2).

When the global terrestrial CO2 sink (average NBP) and its trend (1982–2011) are partitioned among land cover classes, we find that tropical forests account for the largest fraction (26%, 0.33 Pg C year−1) of the average sink over this period (1.23 Pg C year−1) (Fig. 1D). In contrast, we find that semi-arid ecosystems dominate the positive global CO2 sink trend (57%, 0.04 Pg C year−2; global, 0.07 Pg C year−2) (Fig. 1E). The TRENDY model ensemble shows a consistent pattern, with tropical forests dominating the mean sink (median 24%) and semi-arid ecosystems dominating the trend (median 51%). The predominance of semi-arid ecosystems in explaining the global land sink trend is consistent with widespread observations of woody encroachment over semi-arid areas (17) and increased vegetation greenness inferred from satellite remote sensing over recent decades (1719). Likewise, a recent study attributes the majority of the record land sink anomaly of 2011 to the response of semi-arid ecosystems in the Southern Hemisphere, Australia in particular, to an anomalous wet period; the study further postulates a recent increase in the sensitivity of carbon uptake to precipitation for this region, which is attributed to vegetation expansion (20).

We further partitioned IAV in global NBP among land cover classes according to the contribution of individual regions (grid cells or land cover classes) to global NBP IAV (12). To this end, we adopted an index that scores individual geographic locations according to the consistency, over time, with which the local NBP flux resembles the sign and magnitude of global NBP (fig. S4): Embedded Image (1)where xjt is the flux anomaly (departure from a long-term trend) for region j at time t (in years), and Xt is the global flux anomaly, so that Xt = ∑ jxjt. By this definition fj is the average relative anomaly xjt/Xt for region j, weighted with the absolute global anomaly |Xt|. Regions receiving higher and positive average scores are inferred to have a larger contribution in governing global NBP IAV, as opposed to regions characterized by smaller or negative (counteracting) scores (fig. S3). The index we adopt does not characterize the variability of ecosystems of different regions, as, for example, the standard deviation would do (fig. S5); rather, it enables a comparison of their relative importance (contribution) in governing global IAV.

Semi-arid ecosystems were found to account for the largest fraction, 39%, of global NBP IAV, exceeding tropical forest (19%), extratropical forest (11%; all forest, 30%), and grasslands and croplands (27%) (Fig. 1F). The TRENDY model ensemble shows a similar partitioning, with semi-arid ecosystems accounting for 47% (median; tropical forests, 28%; extratropical forest, 6%; all forest, 35%). The overall contributions per land cover class are the sum of both positive and negative contributions that result from differences in phase between IAV of individual grid cells compared with global IAV (fig. S4). The extent to which negative contributions reduce the overall land cover class contributions is minor for all regions except grasslands and croplands (fig. S6) (LPJ-GUESS, –13%; TRENDY median, –13%) because the latter are distributed widely across climate zones, and because both climate variations and the sensitivity of NBP to climate variations differ among regions.

To partition the global NBP IAV among component fluxes (GPP, Reco, Cfire) and among land cover classes, we applied Eq. 1. We found that global NBP IAV is most strongly associated with variation in GPP; interannual GPP anomalies contribute 56% of the global NBP IAV in LPJ-GUESS and a median of 90% in the TRENDY model ensemble. Comparing different land cover classes, the GPP anomalies of semi-arid ecosystems alone contribute 39% in LPJ-GUESS and a median of 65% in the TRENDY model ensemble to global NBP IAV (fig. S7). Semi-arid vegetation productivity thus emerges clearly as the single most important factor governing global NBP IAV.

We used two complementary methods to attribute the variability in GPP—as the inferred primary driver of global NBP IAV—to its environmental drivers. First, we analyzed simulation results from LPJ-GUESS, linking output GPP anomalies to variability in the climatic input data. Second, we used a time-resolved gridded global GPP product derived from upscaled flux tower measurements (12, 21) (hereafter, empirical GPP product). This product uses an empirical upscaling of flux measurements and is thus entirely independent of the modeled GPP in our study.

The three main climatic drivers—temperature (T), precipitation (P), and shortwave radiation (S)—are interdependent and correlated. To account for the combined effects of these drivers, we adopted an analysis of GPP variations from an “impact perspective” (2224): We first identified GPP anomalies and then extracted their climatic covariates. The primary challenge of such an analysis on an annual scale is to target climate indices that adequately characterize the “period of climatic influence” (e.g., growing season average, annual averages, minima or maxima of a given climatic forcing). To overcome this challenge, we used semiannual time series of climate drivers constructed via an optimization procedure that weights monthly anomalies of a given climate variable (T, P, or S), accounting for time lags of up to 24 months while making no additional prior assumptions as to the period of influence (12). For each GPP event, we extracted climatic covariates as z scores of the semiannual climatic drivers.

We evaluated the climatic covariates of GPP anomalies for semi-arid ecosystems from the empirical GPP product and modeled by LPJ-GUESS, focusing on T and P, and found similar responses of GPP to climate with both approaches across all latitude bands (Fig. 2, A and B). Negative GPP anomalies in semi-arid ecosystems are mainly driven by warm and dry (low rainfall) climatic events in most latitudes, suggestive of drought. By contrast, positive GPP anomalies are dominated by cool and wet conditions. Averaging the distributions over latitudes (Fig. 2, A and B) and extracting the climatic covariates per percentile of the GPP distributions shows that GPP varies with climatic conditions on a straight line in T-P space (fig. S8), with a stronger covariation with P than with T. This implies that the full GPP distributions are driven by similar climatic patterns—that is, anomalies that differ in size and sign covary with corresponding differences in size and sign in the drivers. GPP extremes (the tails of the distribution of GPP among years) covary with El Niño–Southern Oscillation (ENSO) across all latitudes (Fig. 2, C and D). In both the model and the empirical GPP product, GPP anomalies are more strongly associated with the positive phase of ENSO (El Niño) than with the negative phase (La Niña); the sign of the relationship varies with latitude. Positive ENSO tends to coincide with negative GPP anomalies in the tropics (30°S to 20°N) and with positive GPP anomalies north of 20°N.

Fig. 2 Climatic covariates of semi-arid ecosystem GPP variations.

(A) Distribution by latitude of the empirical GPP product anomalies normalized by average standard deviation of GPP in semi-arid lands. The distribution is colored according to average local climatic covariates per latitude zone and distribution bin. (B) LPJ-GUESS GPP distribution calculated and colored as in (A). (C) Covariation of the multivariate ENSO index [MEI (31, 32)] anomalies with the empirical GPP product. (D) Covariation of MEI and modeled GPP anomalies per latitudinal zone. Note that the figure shows the covariates of latitudinal average local GPP anomalies, and not the average covariates based on GPP IAV contribution to NBP IAV.

The agreement between climatic covariates of the data-based empirical GPP product and modeled GPP alongside the comparatively robust pattern of the covariation with climate suggests that GPP IAV for semi-arid ecosystems is mediated by climate. Because ENSO covaries with a considerable portion of the GPP distribution, we infer that ENSO is the dominating mode of global circulation variations driving GPP IAV over semi-arid ecosystems. Recent modeling studies have found that extreme El Niño events could become more common under climate change (25), which, together with an increased atmospheric demand for water associated with global warming, might exacerbate the impact of El Niño events over semi-arid ecosystems and further increase the role of semi-arid regions in driving global NBP IAV (2628).

We repeated the calculation of climatic covariates to simulated NBP for LPJ-GUESS and each of the TRENDY models. The resulting maps of covariates in T-P space are shown as average covariates of negative NBP extremes (Fig. 3, A and B) and positive NBP extremes (Fig. 3, C and D). In general, semi-arid ecosystems stand out as regions in which strong CO2 uptake events are consistently associated with cool and moist conditions, and strong CO2 release events with warm and dry conditions. In tropical forests, NBP covaries with both T and P as in semi-arid regions, but also with T alone. In high latitudes, wet or warm and wet conditions lead to negative NBP extremes, whereas dry or warm and dry conditions tend to lead to positive extremes, although the spatial heterogeneity of the covariates is large in this region (Fig. 3).

Fig. 3 Climatic covariates of NBP extremes.

(A) Climatic covariates of LPJ-GUESS negative NBP extremes (1st to 10th percentiles). (B) Mean climatic covariates of TRENDY models’ negative NBP extremes (1st to 10th percentiles). (C) Covariates of LPJ-GUESS positive NBP extremes (90th to 99th percentiles). (D) Mean climatic covariates of TRENDY models’ positive NBP extremes (90th to 99th percentiles).

Our approach offers detailed spatial and temporal disaggregation of drivers and responses, which is important when analyzing drivers or covariates of global NBP IAV because of the high temporal and spatial variability in P (figs. S9 to S11). Using four upscaling levels with increasing spatial and temporal disaggregation [ranging from land surface mean P and T to semiannual P and T, averaged according to the spatial origin of each year’s global NBP anomaly (eqs. S5 and S6)], we found that P and NBP IAV become more correlated at higher levels of disaggregation. At the highest disaggregation level, P is almost as strongly correlated with NBP IAV as T, suggesting a strong influence of soil moisture variations on global NBP IAV (28). This strong increase in P correlations with disaggregation resolves an apparent conflict between our findings and those of studies using regionally averaged drivers that emphasize the role of T in governing IAV in atmospheric CO2 (2830). For semi-arid ecosystems, T correlations with NBP IAV are slightly stronger than P correlations with NBP IAV (Fig. 4B), partly because of an asymmetric distribution of P and/or an asymmetric response of NBP to P IAV (fig. S12). The correlation of tropical forest P with NBP IAV increases when we use the semiannual drivers, which suggests the importance of accounting for time lags and the “period of climatic influence” of P variations (12), but P correlations with NBP IAV are still weaker than T correlations with NBP IAV (Fig. 4C).

Fig. 4 Correlations between climatic drivers IAV (P and T) and global NBP IAV (mean of all 10 models).

(A) Global P and T correlations to global NBP IAV. From black to white and left to right, bars represent annual P and T IAV correlations to global NBP IAV with increasing spatial and temporal disaggregation of P and T while averaging to global time series. Black bars represent averaged global land surface P and T weighted by grid cell area; dark gray bars represent P and T weighted by 30-year average contribution to global NBP IAV (Eq. 1 and fig. S4); light gray bars represent averaged P and T weighted by each year’s contributions, thus accounting for the difference in the spatial distribution of contributions between years (eqs. S5 and S6); white bars represent semiannual climate drivers averaged to global time series using the annual spatial contributions (as for light gray bars), thereby accounting for the “period of climatic influence” and time lags of up to 24 months. (B) Correlations between P and T IAV and global NBP IAV for semi-arid ecosystems. Weights, where applicable, are based on contributions to global NBP IAV as in (A) but with P and T averaged over semi-arid ecosystems only. (C) Correlations between P and T IAV and global NBP IAV for tropical forest. Weights, where applicable, are based on contributions to global NBP IAV as in (A) but with P and T averaged over tropical forest only.

Our analysis provides evidence that semi-arid ecosystems, largely occupying low latitudes, have dominated the IAV and trend of the global land carbon sink over recent decades. Semi-arid regions have been the subject of relatively few targeted studies that place their importance in a global context. Our findings indicate that semi-arid regions and their ecosystems merit increased attention as a key to understanding and predicting interannual to decadal variations in the global carbon cycle.

Supplementary Materials

Materials and Methods

Figs. S1 to S12

References (3356)

References and Notes

  1. See supplementary materials on Science Online.
  2. Acknowledgments: This paper is dedicated to the memory of Michael Robin Raupach (1950–2015), whose scientific integrity and novel contributions leave a long-lasting legacy in the field of carbon cycle sciences. The MODIS MOD12C1 land cover product was obtained through the online Data Pool at the NASA Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota ( Supported by the Royal Physiographic Society in Lund (Birgit and Hellmuth Hertz Foundation), Swedish Research Council grant 637-2014-6895, and the Mistra-SWECIA program (A. Ahlström); EC FP7 grant LUC4C (603542) (A. Arneth); OCE Distinguished Visiting Scientist to the CSIRO Ocean and Atmosphere Flagship, Canberra (B.S.); EC FP7 grant EMBRACE (282672) (A. Arneth, M.R., and B.D.S.); the Australian Climate Change Science Program (J.G.C.); NSF grant AGS 12-43071, U.S. Department of Energy grant DE-SC0006706, and NASA LCLUC program grant NNX14AD94G (A.K.J.); the Environmental Research and Technology Development Fund (S-10) of the Ministry of Environment of Japan (E.K.); CSIRO strategic research funds (Y.P.W.); and NOAA grants NA10OAR4310248 and NA09NES4400006 and NSF grant AGS-1129088 (N.Z.). This study is a contribution to the Lund Centre for Studies of Carbon Cycle and Climate Interactions (LUCCI) and the strategic research areas MERGE and BECC.
View Abstract

Stay Connected to Science

Navigate This Article