Accelerating rates of Arctic carbon cycling revealed by long-term atmospheric CO2 measurements.

The contemporary Arctic carbon balance is uncertain, and the potential for a permafrost carbon feedback of anywhere from 50 to 200 petagrams of carbon (Schuur et al., 2015) compromises accurate 21st-century global climate system projections. The 42-year record of atmospheric CO2 measurements at Barrow, Alaska (71.29 N, 156.79 W), reveals significant trends in regional land-surface CO2 anomalies (ΔCO2), indicating long-term changes in seasonal carbon uptake and respiration. Using a carbon balance model constrained by ΔCO2, we find a 13.4% decrease in mean carbon residence time (50% confidence range = 9.2 to 17.6%) in North Slope tundra ecosystems during the past four decades, suggesting a transition toward a boreal carbon cycling regime. Temperature dependencies of respiration and carbon uptake suggest that increases in cold season Arctic labile carbon release will likely continue to exceed increases in net growing season carbon uptake under continued warming trends.

INTRODUCTION

The Arctic has experienced unprecedented changes in recent decades (), including rapid and diverse changes in Arctic ecosystems such as shrub cover expansion, enhanced vegetation productivity (Arctic greening), longer growing season, and permafrost thaw (–). These ecosystem changes influence rates of carbon cycling and the net Arctic carbon balance (). For instance, increasing carbon release due to permafrost thaw or soil warming could significantly increase atmospheric CO2 concentrations, whereas increasing carbon uptake due to enhanced Arctic vegetation productivity could significantly reduce atmospheric CO2. In either case, Arctic carbon-climate feedbacks will drive global atmospheric CO2 concentrations and alter Earth’s climate trajectory ().

Long-term atmospheric CO2 surface measurements reveal an increase in the CO2 seasonal amplitude in high latitudes (), although it is unclear whether this increase is due to enhanced carbon uptake linked to increased productivity () or to accelerated decomposition of soil organic matter (). With vast amounts of carbon stored in both the active layer and the underlying permafrost, an increase in soil carbon turnover could irreversibly transform the Arctic from a long-term sink to a long-term source (–). We hypothesize that climate change has increased both productivity and respiration and decreased the mean carbon residence time within Arctic ecosystems. Placing observational constraints on Arctic carbon residence time is therefore key to understanding the evolution of Arctic carbon balance and disentangling changes in productivity and respiration.

Direct measurements of pan-Arctic or regional northern high-latitude terrestrial carbon exchange are extremely challenging (). Here, we place observational constraints on the tundra ecosystem carbon balance using the long-term atmospheric CO2 record from Barrow, Alaska (, ). We constrain the northern Alaska region (or “North Slope”) land-atmosphere CO2 exchange [that is, net ecosystem exchange (NEE)] for the period 1974–2015 by deconvolving background-influenced CO2 variability from the Barrow CO2 time series (see Materials and Methods). We treat monthly CO2 anomalies (ΔCO2) as a proxy for North Slope regional-scale NEE. Following Commane et al. (), we assume that background variations are predominantly influenced by large-scale ecosystems above 60°N, including tundra, boreal forest, and other Arctic ecosystems. We evaluate the long-term changes in intra- and interannual variations of ΔCO2 to characterize the long-term carbon dynamics over the Arctic regions. Here, negative (positive) ΔCO2 values indicate regional carbon uptake (release) relative to background CO2.

RESULTS AND DISCUSSION

We find a 0.04 ± 0.022 parts per million (ppm)/year (P < 0.05) increase in the seasonal amplitude of ΔCO2 (calculated as the annual maximum-to-minimum ΔCO2 difference) for the period 1974–2015 (Fig. 1, A and B). The increase in seasonal amplitude emerges mainly from an increasing early cold season (September and October) ΔCO2 (+0.031 ± 0.02 ppm/year, P < 0.05), combined partially with a decreasing early summer (June and July) ΔCO2 (−0.01 ± 0.01 ppm/year, P < 0.5) throughout the 42-year record. The observed amplification of seasonal ΔCO2 variability is consistent with observed changes in CO2 seasonal amplitude in northern mid-latitudes (). Previous studies attribute ΔCO2 amplitude changes to an increase in summertime carbon uptake by enhanced vegetation productivity () due to tundra greening (), shrubification (), a longer growing season (), and northward migration of the boreal forest (). However, these studies do not explicitly account for background CO2 variability, which is dominated by CO2 transport from low latitudes and imparts a phase delay in the seasonal CO2 signal (, ). Our results indicate that local changes in seasonal CO2 amplitude are mainly due to early cold season respiration and that nongrowing season CO2 fluxes are an increasingly important component of the annual carbon balance (, ).

Fig. 1

Changes in seasonal variations of ΔCO2 (CO2local − CO2background) for the last four decades (1974–2015).

(A) Monthly mean ΔCO2 for the period 1974–1983 and that for the period 2006–2015. (B) Changes in annual amplitude of ΔCO2 (maximum ΔCO2 − minimum ΔCO2). (C) Changes in annual ΔCO2 (sum of monthly ΔCO2). Gray dots in (A) show daily values of ΔCO2. The asterisk indicates statistical significance at the 95% confidence level.

We find that variations in early cold season ΔCO2 are strongly correlated to the soil temperature, especially in upper soil layers (0 to 40 cm; P < 0.01) (fig. S4). The positive relationship between near-surface soil temperature and ΔCO2 is consistent with anticipated soil carbon losses due to warming (). In addition to early cold season warming, positive correlation between summertime (June, July, and August) averaged normalized difference vegetation index and early cold season (September to October) ΔCO2 shows a direct contribution of Arctic greening and increased productivity to early cold season ΔCO2 (fig. S5). This is consistent with a previous field study over the Arctic region (), supporting the hypothesis that the increase in respiration of recent labile soil organic matter is proportional to the increase in carbon input by enhanced vegetation productivity (, ).

In accordance with trends in seasonal ΔCO2, annual ΔCO2 exhibits a positive trend over the past 40 years (+0.132 ± 0.064 ppm/year, P < 0.01; Fig. 1C). Assuming that annual ΔCO2 is related to annual NEE, this result suggests that, despite the compensation between the increasing trends in carbon uptake (for example, negative trends in ΔCO2 during summer) and loss (for example, positive trend in ΔCO2 during early cold season), the mean carbon balance has significantly changed over time, given a stasis of seasonal changes in the wind speeds and boundary layer heights. Therefore, on the basis of 1974–2015 annual ΔCO2 trends, our results indicate that carbon release linked to increased soil thaw is increasing. However, given changes in the seasonal boundary layer height, it is difficult to assess whether soil carbon losses are outpacing carbon uptake due to enhanced vegetation productivity.

The carbon residence (or turnover) time is an emergent ecosystem property that diagnoses the interplay of climate, carbon fluxes, soil, and vegetation (, ). Shorter residence times reflect ecosystems that can respond more rapidly to changing climate (, , ). Residence time depends on the intrinsic stability of the carbon stock and the environmentally limited (that is, soil temperature and moisture) rates of plant and microbial respiration (, ). We evaluate integrated dynamics of northern Alaska carbon cycling by combining observed ΔCO2 and an ecosystem carbon balance model to retrieve the decadal trends in mean carbon residence time over Alaska. We use a Bayesian approach with observed ΔCO2 anomalies and an ecosystem carbon balance model (see Materials and Methods) and find a 99% probability that carbon residence time decreased throughout the study period. Mean carbon residence time for 2004–2013 is 13.4% lower than that for 1979–1988 (50% confidence interval, 9.2 to 17.6%, Fig. 2). In comparison, the Multi-scale Synthesis and Terrestrial Model Intercomparison Project (MsTMIP) ensemble of terrestrial biosphere models indicates a mean decrease of 8.8% (50% confidence interval, 5.0 to 12.8%; Fig. 2) in residence time for the same period (). Although process-based models exhibit a weaker response, they do indicate a consistent decrease in carbon residence time across Arctic regions. In addition, our estimate for northern Alaska (50% confidence interval, 23 to 782 years; median, 139 years) is larger but broadly consistent with gridded residence time estimates () (range, 29 to 180 years; median, 55 years) and roughly an order of magnitude greater than tropical ecosystem mean residence times (, , ).

Fig. 2

Retrieved changes in carbon residence time based on the difference between 1979–1988 and 2004–2013 10-year periods.

The vertical red line indicates the average of retrieved changes in carbon residence time. The blue line indicates mean (solid line) and range (shading) in the equivalent residence time change estimates from the MsTMIP model ensemble.

Our finding of a 13.4% decrease in the carbon residence time suggests that Alaska’s North Slope tundra is becoming more boreal. A decrease in high-latitude residence time indicates an increased sensitivity of Arctic carbon cycling to warming and, therefore, an enhanced role of Arctic ecosystems in the variability of atmospheric CO2. Although several processes can contribute to residence time changes, we hypothesize that residence time changes are occurring as a result of changes in dead organic carbon turnover rates and/or trends in the live biomass/dead organic carbon ratio. Carbon held in permafrost is effectively permanently stored, with a minuscule decomposition rate; therefore, permafrost thaw transfers carbon to a far more active phase and effectively reduces the average carbon residence time. Shrubification and forest migration would likely lead to higher carbon allocation to long-lived woody carbon pools (), therefore leading to a likely increase or no change in the residence time of live biomass. However, biomass accumulation due to shrubification would lead to a higher live biomass/dead organic carbon ratio, leading to a decrease in mean ecosystem carbon residence time, since live biomass residence times are roughly an order of magnitude lower than those of dead organic carbon matter ().

Our analysis also indicates a 54% probability that soil (heterotrophic) respiration is more temperature-sensitive than net carbon uptake (fig. S3), resolving a key uncertainty in the response of Arctic ecosystems to climate change. The enhanced response of soil decomposition and respiration to warming and deepening of the soil active layer implies net ecosystem carbon loss under elevated temperatures in the 21st century. To gain confidence in the sensitivity of carbon respiration and uptake to climatic variability, the effects of permafrost mobilization and priming (), biomass accumulation, and subsequent shifts in live biomass inputs into dead organic carbon pools need to be explored further; however, these processes will likely contribute to a sustained or increasing future net carbon loss from Arctic ecosystems. Soil moisture and precipitation trends will likely influence the rate of permafrost degradation, decomposition rates, and the CO2:CH4 of heterotrophically respired carbon (); these carbon-water relationships remain uncertain and highlight the need for an improved understanding of the integrated long-term carbon dynamics and carbon balance sensitivity to hydrological variability in the Arctic system ().

Observationally constrained carbon residence time estimates can help narrow uncertainties of carbon cycle predictions () and improve understanding of Arctic carbon-climate feedbacks. Absolute CO2 flux constraints based on past and current CO2 observation systems will ultimately help resolve the processes controlling long-term variations in terrestrial carbon exchanges. We anticipate that regional top-down constraints on ecosystem carbon cycling—including atmospheric CO2 flux constraints (), such as satellite, airborne, and tower solar-induced fluorescence (, ), and land-surface data constraints ()—will together provide an enhanced process understanding of Arctic terrestrial carbon cycling sensitivity and vulnerability to long-term climate trends.

MATERIALS AND METHODS

Regional land attributions of observed CO2 anomalies (ΔCO2)

We used half-hourly CO2 concentration measurements at Barrow, Alaska, from the National Oceanic and Atmospheric Administration (NOAA) Earth System Research Laboratory archive for the period 1974–2015 (). We discriminated between local CO2 (CO2local) and background CO2 concentrations (CO2background) based on wind direction, wind speed, and time of day (, ). Sweeney et al. () applied the method to separate land sources of CH4 from the clean air sector. Commane et al. () verified this method to estimate land sources of observed CO2 at Barrow by using the Hybrid Single-Particle Lagrangian Integrated Trajectory Model driven by meteorological data from NOAA. Locally influenced CO2 concentrations (CO2local) were derived on the basis of 12:00 p.m. to 4:00 p.m. half-hourly CO2 measurements where wind speed exceeded 3 m/s and wind direction was between 150° and 210°. Background CO2 concentrations (CO2background) were derived on the basis of 12:00 p.m. to 4:00 p.m. half-hourly CO2 measurements where wind speed exceeded 3 m/s and wind direction was between 0° and 90°. Finally, we defined local attributions of CO2 changes as ΔCO2 = CO2local − CO2background. More details in this method are described in Commane et al. () and Sweeney et al. ().

We selected data from 12:00 p.m. to 4:00 p.m. during the day because boundary layer height is the highest during this time. A higher boundary layer expands the surface area of influence significantly, resulting in a larger volume of air that is mixed from the surface to the top of the boundary layer. Thus, air that may have been at the surface a day ago but has escaped the nighttime boundary layer can be reincorporated during this time of the day. As the boundary layer stratifies during other times of the day, the area of influence is significantly reduced, making the change from background more local. It should also be noted that, because boundary layer is higher between 12:00 p.m. and 4:00 p.m. than any other time of the day, it will also be more diluted by the overlying free tropospheric air that sits over the boundary layer.

Daily measurements show a mix of negative and positive anomalies (ΔCO2) over the first and last decades of our study (1974–1983 and 2006–2015, respectively; Fig. 1). However, on average, monthly ΔCO2 are predominantly positive over the course of the year. This pattern is mainly due to the limited duration of the drawdown period and variability of that period from year to year, making the monthly average positive despite many negative daily observations. It also should be noted that the delayed growing season drawdown in tundra (summer) relative to southern boreal and temperate forests (spring) will result in pan-Arctic “background” CO2 significantly depleted relative to that observed at Barrow (, ).

Model-data fusion

We retrieved the change in ecosystem carbon residence time by assimilating June and September 1979–2015 ΔCO2 into a single-pool ecosystem carbon balance model. We used a Bayesian model-data fusion approach () to optimize five model parameters relating to carbon uptake, initial ecosystem carbon stock, turnover rate, and the temperature sensitivities of carbon uptake and respiration. We used a Metropolis-Hastings Markov Chain Monte Carlo algorithm to sample 5000 model parameter sets. We evaluated the model results against a process-based terrestrial biosphere model ensemble (), flux tower–based estimates of carbon uptake (), and residence time retrievals (, ).

In the following sections, we describe the optimization of a first-order ecosystem carbon balance model scheme for the period 1979–2015. The temporal change in mean ecosystem C residence time—and its associated uncertainty—is derived on the basis of the optimized model output.

Ecosystem carbon balance model

We expressed the total active terrestrial ecosystem carbon C within the northern Alaska region (longitude, 170°W to 150°W; latitude, 65°N to 70°N) at time t + 1 aswhere NPP is the net primary productivity, R is the heterotrophic respiration, and Δt is the time step. We represent the first-order sensitivity of NPP to temperature (T), leaf area index (LAI), and global radiation (S) as followswhere p1 is a baseline productivity parameter and p2 is an exponential temperature dependence parameter. S and T (2-m air temperature) data were derived from ERA-interim 1979–2015 1° × 1° reanalysis fields. LAI was obtained from Zhu et al. (). R was derived as a function of temperature, wherewhere p3 is a baseline turnover rate and p4 is an exponential temperature dependence parameter. At each time step, the NEE isOur model is a first-order representation of the integrated ecosystem carbon dynamics and their sensitivity to leaf area, solar radiation, and temperature. We noted that our NPP derivation assumes (i) no hydrological limitations on carbon uptake in Arctic ecosystems, (ii) a linear relationship with respect to LAI and S, and (iii) respiration dependence on the 2-m air temperature (instead of surface or subsurface temperature). We determined the viability of our approach by comparing northern Alaska–optimized monthly NPP against a range of process-based terrestrial biosphere models (, ). The optimization of model parameters and the comparison of our model with process-based models are described in the following section. Our model assumes that carbon losses due to disturbance are a minimal component of the Alaska North Slope carbon balance: Land management and land-use change are minimal components of the Arctic ecosystem carbon cycling, while mean C fire emissions (based on the global fire emission database, version 4) () account for <0.19 grams of carbon (gC) m−2 month−1 for the period 2001–2013, roughly two orders of magnitude smaller than mean uptake and respiration fluxes (see fig. S2).

Model optimization

We optimized model parameters based on 1979–2013 monthly ΔCO2 observations; we henceforth referred to model parameters—p1, p2, p3, p4, and C0—as x and to monthly ΔCO2 observations as O. We used Bayesian inference to derive the probability distribution of x relative to observational constraints O, p(x|O), wherep(x) is the prior probability of x, and p(O|x) is the likelihood of x given O. We used a Markov Chain Monte Carlo approach () to extract 5000 samples of x. In addition to broad uniform prior parameter ranges (table S1), we (i) prescribed a 10 gC m−2 month−1 constraint on mean annual NPP with a factor of 2 uncertainty, which is broadly consistent with the MsTMIP [BG1 simulations, comparison in fig. S2; see Materials and Methods (Comparison with process-based terrestrial biosphere model ensemble) for MsTMIP experiment details], and (ii) ensured that |NEE| is statistically within 1 SD of mean NPP. Hence, for each Monte Carlo iteration, the prior probability of parameter sample x, p(x), is evaluated aswhere prange(x) is 1 if x is within the prior parameter ranges (or zero otherwise), and are the mean model NPP and NEE for the period 1979–2013, respectively (derived as shown in Eqs. 2 and 4), and NPP0 = 10 gC m−2 month−1. We defined the likelihood of x given observational constraints ΔCO2 aswhere nee and a are the vectors of monthly NEE and ΔCO2 standardized anomalies for month m, respectively. For each m, standardized anomalies nee (a) were derived by removing mean 1979–2013 NEE (ΔCO2). We hence assumed a linear relationship between ΔCO2 and NEE for each month. By standardizing monthly anomalies independently, we did not assume that monthly ΔCO2 anomalies were comparable. For each month, the model-data residual (nee − a) was weighed by uncertainty σ. We weighed the model-data residuals (nee – a) by σ = 10 SD(a), where SD(a) represents the SD of the bracketed vector. We noted that the uncertainty and error covariance structure of a and the representation error of nee are largely unknown: To avoid overconstraining the model, we chose to prescribe σ = 10 SD(a), since we found that prescribing σ < 10 SD(a) did increase model confidence but without a substantial reduction in model-data residuals.

Standardized observations and optimized standardized model NEE are shown in fig. S1. The details of the Metropolis-Hastings Markov Chain Monte Carlo algorithm are described in Bloom and Williams () and references therein. Both model NPP and heterotrophic respiration peak in July, which is consistent with half (four of eight) of the MsTMIP terrestrial biosphere models (fig. S2).

Retrieved residence time change and temperature sensitivities

We derived the percentage of residence time change (Δτ) as followswhere τstart and τend correspond to the mean ecosystem C residence times based on the 1979–1988 and 2004–2013 10-year periods. The residence time is derived in accordance with Bloom et al. (), where each τ is derived aswhere and are the 10-year mean C stocks and heterotrophic respirations, respectively, for the corresponding τ time spans. The range of Δτ values is derived by repeating Eq. 8 for each model parameter vector sample x. On the basis of Δτ outcomes for all accepted parameter vectors, we found that the probability of Δτ < 0 is 99%. The distribution of the optimized temperature sensitivities of NPP and R (parameters p2 and p4; table S1) is shown in fig. S3. On the basis of all accepted parameter vector samples x, we found that the probability of p4 > p2 is 54%.

Comparison with process-based terrestrial biosphere model ensemble

We compared the ecosystem carbon balance model NPP, heterotrophic respiration, and retrieved mean residence times with the MsTMIP terrestrial biosphere model ensemble outputs [V1.0 MsTMIP outputs ()]. To compare model C residence times (see section S4), we limited our comparison to BG1 simulation models in MsTMIP, which provided both heterotrophic respiration and total soil carbon outputs (table S2): We noted that since τ is dependent on total modeled C stocks and fluxes during a certain time period, it provides a complexity-independent first-order metric of the rate at which C is cycled through the terrestrial biosphere. The BG1 simulations include time-varying climate, nitrogen deposition, atmospheric CO2, and land-use history: For the sake of brevity, we refer the reader to Huntzinger et al. () for individual model details.