Large loss of CO2 in winter observed across the northern permafrost region.

Recent warming in the Arctic, which has been amplified during the winter1-3, greatly enhances microbial decomposition of soil organic matter and subsequent release of carbon dioxide (CO2)4. However, the amount of CO2 released in winter is highly uncertain and has not been well represented by ecosystem models or by empirically-based estimates5,6. Here we synthesize regional in situ observations of CO2 flux from arctic and boreal soils to assess current and future winter carbon losses from the northern permafrost domain. We estimate a contemporary loss of 1662 Tg C yr-1 from the permafrost region during the winter season (October through April). This loss is greater than the average growing season carbon uptake for this region estimated from process models (-1032 Tg C yr-1). Extending model predictions to warmer conditions in 2100 indicates that winter CO2 emissions will increase 17% under a moderate mitigation scenario-Representative Concentration Pathway (RCP) 4.5-and 41% under business-as-usual emissions scenario-RCP 8.5. Our results provide a new baseline for winter CO2 emissions from northern terrestrial regions and indicate that enhanced soil CO2 loss due to winter warming may offset growing season carbon uptake under future climatic conditions.

Air and soil temperatures in the Arctic are increasing rapidly, with the most severe climate amplification occurring in autumn and winter[1,2]. Although warmer soils decompose more quickly, thus releasing more CO2 into the atmosphere, microbial respiration is known to occur even under extremely cold winter conditions (e.g., down to ~ −20°C) in unfrozen microsites that can persist at sub-zero soil temperatures[7]. This production and release of CO2 in winter is expected to increase substantially as soils continue to warm and thaw under a warming climate[4,8].

However, it remains highly uncertain how much CO2 is currently emitted from the permafrost region during winter[9] and how much these emissions might increase in the future[8,10]. Many ecosystem models are not well adapted to simulate respiration from high latitude soils[5] and may greatly underestimate present and future winter CO2 emissions[6]. Given the limitations in current models, lack of satellite and airborne CO2 data for the Arctic during winter[11], and gaps in spatial coverage of Arctic air monitoring networks[12], in situ CO2 flux observations provide the most direct insight into the state of winter CO2 emissions across the northern permafrost domain.

Studies of winter respiration indicate that the amount of CO2 released during cold periods depends greatly on vegetation type[13], availability of labile carbon substrates[14,15,16], non-frozen soil moisture[4,7,15,17,18], microbial community composition and function[19], and snow depth[15, 20, 21]. However, knowledge of the influence of these drivers on the rates and patterns of winter CO2 flux on a regional scale remains limited[6, 9].

Here we present a new compilation of in situ CO2 winter flux data for the northern permafrost domain (Fig. 1, Supplementary Information (SI) Table 1) to examine the drivers and magnitude of winter respiration in the Arctic. We define the winter period as October through April—months when the landscape is generally covered by snow and photosynthesis is negligible [22,23]. The dataset represents more than 100 high latitude sites and comprises more than 1,000 aggregated monthly fluxes. We examined patterns and processes driving winter CO2 emissions and scaled fluxes to the permafrost domain using a boosted regression tree (BRT) machine learning model based on hypothesized drivers of winter CO2 flux. Environmental and ecological drivers (e.g., vegetation type and productivity, soil moisture, and soil temperature) obtained from satellite remote sensing and reanalysis data were used to estimate regional winter CO2 emissions for contemporary (2003–2017) climatic conditions. We estimated winter fluxes through 2100 using meteorological and carbon cycle drivers from ensembles of Earth System Model (ESM) outputs for RCP 4.5 and RCP 8.5[24].

Fig. 1.

Distribution of in situ data included in this winter CO2 flux synthesis.

(a) Locations of in situ winter CO2 flux data (yellow circles) in this synthesis include (b) upland and wetland sites in boreal and tundra biomes located (c) within the northern permafrost region[41]. Violin plots (b,c) depict magnitude and distribution density (width; dots are monthly aggregated data) of in situ data used in our machine-learning model.

Soil temperature had the strongest influence on winter CO2 emissions, with fluxes measured at soil temperatures down to −20°C (Fig. 2a), in line with results from lab incubations (Fig. 2b), demonstrating that microbial respiration may occur in unfrozen microsites that persist at sub-zero bulk soil temperatures[18]. Diffusion of stored CO2 produced during the non-frozen season may have driven some of the emissions measured in winter, but the magnitude of this contribution is unclear. Winter CO2 emissions increased by a factor of 2.9 (95% credible interval (CI) = 2.1, 4.2) per 10°C soil temperature increase (i.e., Q10) for in situ fluxes and by a factor of 8.5 (CI= 5.0, 14.5) for CO2 release from low temperature lab incubations. Differences between in situ and lab Q10s may reflect site-level differences in environmental drivers other than temperature (in situ and lab sites were not fully overlapping), experimental design differences (e.g., less restricted diffusion in the lab), or variation in the depth of in situ CO2 production, which can occur throughout the soil profile, relative to the depth of recorded temperature, which tended to be closer to the soil surface (~ 10 cm).

Fig. 2.

Effect of soil temperature on CO2 release from soils.

(a) Relationships between in situ soil temperature (~ 10 cm average depth) and CO2 fluxes and (b) temperature and CO2 released from lab incubations. Shading represents the standard deviation of an exponential model, which, for in situ fluxes, was fit to mean CO2 flux from each sample location (symbols shown with standard error). Note that the different soil temperature scales between panels reflect data ranges.

Air and soil temperatures had the strongest influence on winter flux with a combined relative influence (RI) of 32%. Vegetation type (15% RI), leaf area index (LAI; 11%), tree cover (TC; 10%), and previous summer’s gross primary productivity (GPP; 8.5%) also influenced winter CO2 emissions (SI Fig. 1). Along with warmer air and soil temperatures in winter and corresponding increases in CO2 loss, summer GPP has also been increasing in some parts of the northern permafrost region[25]. The positive relationship between GPP and winter CO2 emissions suggests that increased CO2 uptake during the growing season may be offset, in part, by winter CO2 emissions.

Another important driver of winter respiration was unfrozen water content, which is a function of soil temperature and texture, as finer textured soils contain more unfrozen water than coarse soils for a given sub-zero temperature[26]. Indirect measurements of unfrozen water availability confirm its importance: soils with low sand and high clay content, which tend to have greater unfrozen microsites, were characterized by higher CO2 flux rates. While snow cover is a key driver of winter flux through its impact on ground temperature[27], remote sensing estimates of snow cover were not significant predictors in the model; this may be a result of high uncertainty in regional snow products or because snow depth and density, which are difficult to determine from space using currently available satellite technology[28], have a greater influence on ground temperatures than snow presence alone.

Using our model to assess winter flux for the terrestrial permafrost domain, we estimate approximately 1662 Tg C winter−1 released under current climatic conditions (2003–2017), with a corresponding uncertainty of 813 Tg C winter−1 (SI Methods). We observed no temporal trends in winter CO2 flux during this 15-year period (p > 0.1), which corresponded with the lack of a significant circumpolar trend in the reanalysis winter air or soil temperature data used as model inputs (p > 0.1). Although we did not observe region-wide trends during the past 15 years, atmospheric CO2 enhancements for Alaska[8] and site-level studies from Alaskan tundra[29,30] showed recent increases in winter emissions, which are already shifting some tundra regions from an annual carbon sink to a source.

Our flux estimates are twofold higher than a previous estimate derived from in situ measurements reported in the Regional Carbon Cycle Assessment and Processes (RECCAP) tundra and northern boreal domain[10], which was based on a much smaller dataset (< 20 site-years for winter data). The RECCAP study reported fluxes of 24 – 41 g C m−2 winter−1 from in situ data, compared to 64 g C m−2 winter−1 estimated here for the RECCAP region and 98 g C m−2 winter−1 for the full permafrost domain (SI Fig. 2). Our estimate of winter flux agrees more closely with the RECCAP atmospheric inversion estimate (27–81 g C m−2 winter−1), providing some closure between bottom-up and top-down assessments[6,12].

We then compared our permafrost region flux estimates to winter net ecosystem exchange (NEE) outputs from five process-based terrestrial models and from FluxCom, a global machine-learning NEE product[31]. Our winter CO2 flux estimate was generally higher than estimates from these models, which ranged from 377 Tg C winter−1 for FluxCom and from 503 to 1301 Tg C for the process models (mean: 1008 Tg C winter−1; SI Fig. 3). Similar variation in carbon budget estimates from terrestrial models has been reported elsewhere for high latitude regions[5], which reflects considerable differences in model parameterization of soil temperature, unfrozen water, and substrate effects on CO2 production under winter conditions. Some process-based models may underestimate winter CO2 emissions by shutting down respiration at sub-zero soil temperatures[32] or because they are unable to capture small-scale processes that influence winter flux, such as talik formation and shrub-snow interactions that are more likely to be captured by in situ measurements.

Combining growing season NEE (−687 to −1647 Tg C season−1) and winter NEE derived from the process-based terrestrial models described above results in an estimated annual NEE of −351 to 514 Tg C yr−1 (−555 for FluxCom; SI Table 2). Because our winter emissions estimate was higher than these process models, we expect that annual CO2 losses may also be higher. For example, if we account for growing season NEE using the process model estimates, this would yield an average annual CO2 emission of 646 Tg C yr−1 (range of 15 to 975) from the permafrost region, based on our estimate of winter CO2 flux.

Our assessment of future winter emissions—obtained by forcing the BRT model with environmental conditions from CMIP5 ESM outputs[2]—showed significant increases in winter CO2 emissions under both climate scenarios (p < 0.001, Fig. 3); however, emissions were substantially lower with climate mitigation in RCP 4.5 than with RCP 8.5. Compared to current winter emissions (2003–2017), there was a 17% projected increase in winter CO2 flux under RCP 4.5 by 2100 (to 1950 Tg C yr−1) and a 41% increase under RCP 8.5 by 2100 (to 2345 Tg C yr−1) (Fig. 4).

Fig. 3.

Pan-Arctic winter CO2 emissions under current and future climate scenarios.

(a) Average annual winter (October - April) CO2 emissions estimated for the permafrost region for the baseline years 2003–2017. Cumulative winter CO2 fluxes under (b) RCP 4.5 and (c) RCP 8.5 scenarios over an 80-year period (2017–2057 and 2057–2097). Fluxes are reported on an annual basis (g CO2-C m−2 yr−1).

Fig. 4.

Projected annual CO2 emissions during the winter for the northern permafrost region.

Solid lines represent BRT modeled results through 2100 under RCP 4.5 (blue solid line) and RCP 8.5 (red solid line), with bootstrapped 95% confidence intervals indicated by shading. For reference, CMIP5 ensemble respiration for RCP 4.5 and 8.5 are also shown (dashed lines).

The present-day continuous permafrost zone experienced the strongest positive trend in winter CO2 emissions under both climate scenarios (p < 0.001); however, accounting for differences in area, the largest rate of change in winter CO2 emissions occurred across the discontinuous zone (SI Table 3) where soils have warmed rapidly and permafrost has diminished in recent years[33]. The differences in projected changes in winter CO2 emissions among permafrost zones may reflect the influence of latitudinal variation in environmental and ecological variables, including tree cover, dominant vegetation, and soil organic matter content and composition[34].

Increased winter CO2 emissions from our data-driven BRT model were largely driven by changes in soil and air temperatures, which both increased by 0.04°C yr−1 under RCP 4.5, and increased by 0.08°C yr−1 for soil and 0.1°C yr−1 for air under RCP 8.5 (SI Fig. 4). Vegetation leaf area and GPP, both of which were positively related to winter CO2 flux, also significantly increased through 2100.

From 2018 to 2100, we estimated a cumulative winter flux of 150 Pg C for RCP 4.5 and 162 Pg C for RCP 8.5. This represents an additional 15 Pg C for RCP 4.5 and 27 Pg C for RCP 8.5 emitted as a result of climate change, when compared to the estimated 135 Pg of C that would be emitted through 2100 if current (2003–2017) climatic conditions remained constant. These losses are comparable to 70% of the current permafrost-region near-surface (0–30cm) soil carbon pool[35]. These projected increases are substantially lower than projections from CMIP5 ESMs, in which winter CO2 emissions from ecosystem respiration for the permafrost region (1753 ± 1066 Pg C yr−1 for 2003–2005) were projected to increase in 2100 by 37% and 86% under RCP 4.5 (2482 ± 1403 Pg C yr−1) and 8.5 (3473 ± 1731 Pg C yr−1), respectively (Fig. 4). Our data-driven BRT model may provide more conservative estimates because current in situ observations may not adequately reflect future environmental responses to substantially warmer winter conditions. However, it is also possible that the ESMs are missing stabilizing drivers and mechanisms that might provide negative feedbacks to winter CO2 emissions. Hence, we stress the importance of addressing current uncertainties in process-model estimates of both growing season and winter CO2 exchange. Given the data limitations during the winter, there is a particular need for long-term monitoring of winter CO2 exchange in permafrost regions to provide key insights into processes that may enhance or mitigate change. As most of the CMIP5 models do not currently include a permafrost component, these data are critical for improving pan-arctic carbon cycle simulations.

Some of the projected winter CO2 emissions could be offset by plant carbon uptake, which is expected to increase as plants respond favorably to warming and CO2 fertilization[36,37]. In addition, our modeled results do not explicitly account for CO2 uptake during the shoulder seasons (early and late winter period, e.g., October and April), which can occur, even under the snowpack[22,23,38] and which may increase with climate warming[22]. Our model projections also did not incorporate all changes expected under future climates, such as changes in permafrost distribution, delayed seasonal freeze-up, increased fire frequency, changes in snow cover and distribution, thermokarst frequency and extent, and landscape-level hydrologic changes (e.g., lake drainage).

The CO2 emissions reported here are only part of the winter carbon budget, which also includes significant CH4 emissions from land[17,39] and CO2 and CH4 emissions from inland waters[40]. Recent data-derived estimates of high-latitude terrestrial winter CH4 emissions range from 1.6 Tg C yr−1 (land area > 60°N)[39] to 9 Tg C yr−1 for arctic tundra[17]. Similar to winter CO2 emissions, process models significantly underestimated the fraction of annual CH4 emissions released during the winter[39].

To reduce uncertainty in estimates of current and future emissions, we recommend increased spatial and temporal coverage, and coordination and standardization of in situ winter measurements, improvements to regional snow density products, and development of remote sensing active sensors that can detect high resolution (< 20 km) changes in atmospheric CO2 concentrations during periods of low to no sunlight, which is a key constraint on efforts to monitor changes in permafrost region carbon cycling. Current rates of winter CO2 emissions may be offsetting CO2 uptake by vegetation across the permafrost region. Circumpolar winter CO2 emissions will likely increase in the near future as temperatures continue to rise; however, this positive feedback on global climate can be mitigated with a reduction of global anthropogenic greenhouse gas emissions.

Methods

Data overview

We compiled a dataset of in situ winter season (Oct-April) CO2 emissions and potential driving variables from sites within the northern permafrost zone[41]. The synthesized dataset included 66 published studies and 21 unpublished studies (SI Table 1) conducted at 104 sites (i.e., sample areas with unique geographic coordinates) and in 152 sampling locations (i.e., different locations within a site as distinguished by vegetation type, landscape position, etc.). Sites spanned boreal and tundra landcover classes (SI Fig. 5, SI Table 4) in continuous permafrost (n=69), discontinuous (n=24), and isolated/sporadic (n=11) permafrost zones (Fig. 1). Data were aggregated at the monthly level; however, the number of measurements per month varied among studies. The dataset included more than 1,000 site-month flux measurements. We also extracted CO2 data from incubations of permafrost-region soils (SI Table 5) to compare their temperature response functions (Q10) with Q10 derived from the synthesized in situ flux data. Further details of data extraction and Q10 calculations can be found in the Supplementary Methods.

Data extraction, geospatial data

We extracted data from regional gridded geospatial products including climatological data, soil temperature and moisture, snow water equivalent, soil carbon stocks and texture, permafrost status, vegetation cover, proxies of vegetation growth and productivity (e.g., enhanced vegetation index, EVI; leaf area index, LAI; gross primary productivity, GPP). See Supplementary Methods for further description and data sources. All geospatial data were re-gridded to the National Snow and Ice Data Center Equal Area Scalable Earth (EASE) 2.0 format[42] at a 25-km spatial resolution prior to the CO2 flux upscaling and simulations.

Boosted regression tree analysis

We used boosted regression tree analysis (BRT) to model drivers of winter CO2 emissions and to upscale emissions to the northern permafrost region under current and future climate scenarios. The BRT model was fit in R[43] using ‘gbm’ package version 2.1.1[44], and using code adapted from[45]. The BRT model was fitted with the following metaparameters: Gaussian error distribution, bag-fraction (i.e., proportion of data used in each iteration) of 0.5, learning rate (contribution of each tree to the final model) of 0.005, and a tree complexity (maximum level of interactions) of two. We used 10-fold cross-validation (CV) to determine the optimal number of trees to achieve minimum predictive error and to fit the final model to the data.

We used geospatial data as explanatory variables in our BRT model (See Supplementary Methods for full description of input data). We removed highly correlated variables from the models (Spearman ρ = 0.7), retaining the variable within each functional category (e.g., air temperature) that had the highest correlation with winter flux. We further reduced the model by removing variables in reverse order of their relative influence, until further removal resulted in a 2% average increase in predictive deviance. We compared this model with one in which we included site level in situ data as explanatory variables. We used the geospatial model because it allowed us to upscale results and because the percent deviance (SI Table 6) and driving variables (SI Fig. 1) were similar between models.

We assessed BRT model performance using: 1. The correlation between predicted and observed values using the CV data (i.e., data withheld from model fitting), hereafter referred to as the CV correlation, and; 2. deviance explained by the model over the evaluation dataset (i.e., CV data), calculated as: % deviance = (CV null deviance - CV residual deviance)/CV null deviance *100. Further details of the BRT models can be found in the Supplementary Methods.

We obtained an estimate of model uncertainty by first obtaining the average internal root mean squared error (RMSE; 0.21 g C m−2 d−1) for the ensemble of boosted regression trees. We then made the assumption that this error applied equally to all grid cell areas within the domain. Scaling this error to the full domain (16.95 × 106 km2) and by the total number of days for the winter (October through April) period provided us with a winter flux error of 813 Tg C winter−1.

Spatial and temporal domain for mapping

We scaled the modeled flux data to the northern permafrost land area ≥ 49° N[41], which comprises 16.95 × 106 km2 of tundra and boreal lands (excludes glaciers, ice sheets, and barren lands; Fig. 1) with lake area removed. We defined the winter period as the months of October through April. Because the climate within this timeframe varies substantially across the permafrost zone, this month-based definition, while temporally consistent, may include some areas that are influenced by climate that would fall outside expected winter temperature ranges. Therefore, in a separate approach (presented in the Supplementary Method), we defined winter based on soil temperature, but we did not find substantial differences in regional flux budgets when using the two approaches (temperature-defined winter flux was ~ 5% higher, 1,743 Tg C, than when using the month-based winter period).

Spatial upscaling of fluxes

The BRT model was applied at a monthly time step from 2003 through 2017. For each month, the map predictions were applied to a raster stack of input predictors using the R ‘dismo’ package[46] for interface with the ‘gbm’ package and the ‘raster’ v2.6–7 predict function for geospatial model applications. A n.tree (# of trees) of 1,000 was selected for each model run. Output monthly mean estimates of daily CO2 flux (g CO2-C m−2 d−1) were generated for each 25-km grid cell. Total pan-arctic CO2 flux was obtained on a monthly basis by first calculating the terrestrial area for each grid cell by subtracting lake fractions (MODIS satellite product MOD44W) from each grid cell area. The fluxes were then scaled according to days per month and terrestrial area to obtain per grid cell totals.

We analyzed the pan-arctic flux data for annual temporal trends using the nonparametric Mann-Kendall test, which was run in the R ‘zyp’ package[47] with pre-whitening (Yue and Pilon method) to remove autocorrelation. We report Kendall’s correlation coefficient, τ, to describe the strength of the time-series and Theil-Sen slope to describe trends over time.

Comparison of BRT estimates with process-based models

We compared our regional winter flux estimates to: 1) outputs from five process-based terrestrial models estimated for the northern permafrost domain: National Center for Atmospheric Research (NCAR) Community Land Model (CLM) versions 4.5 and 5; Lund-Potsdam-Jena Dynamic Global Vegetation Model (LPJ-DGVM), Wald Schnee und Landscraft version (LPJ-wsl); CARbon DAta MOdel FraMework (CARDAMOM); and the NASA SMAP Level 4 Carbon (L4C) Version 3 product; 2) estimates for the northern permafrost domain derived from FluxCom, a global gridded machine-learning net ecosystem exchange (NEE) product; and 3) four process-based terrestrial models and eight atmospheric inversion models from the high latitude model intercomparison for the Regional Carbon Cycle Assessment and Processes (RECCAP) tundra and northern boreal domain[10]. See Supplementary Methods for further description of these models.

Projected CO2 flux

Inputs for the BRT model of future scenarios of winter CO2 flux were obtained from ensembles of Earth System Model (ESM) outputs from the Fifth Coupled Model Intercomparison Project (CMIP5) for RCP 4.5 and 8.5[2]. Inputs included: 1) Annual GPP; 2) mean annual summer LAI (July & August); 3) mean summer soil moisture (June, July, August); 4) mean monthly soil moisture; 5) mean monthly near-surface (2 m) air temperature; and 6) mean monthly soil temperature (layer 1) (SI Table 7). Ensemble mean RCP 4.5 and 8.5 predictor fields were bias-corrected using the delta, or perturbation method[48], based on historic ESM outputs and observed historical data and re-projected to EASE2 25 km grids.

In addition to the 0.21 g C m−2 d−1 error obtained based on the BRT model RMSE, we used the outcome from bootstrapped BRT model simulations to estimate additional, inherit prediction variability in the machine learning outcomes for current and future CO2 emissions[49] (see Supplementary Information).

For the CMIP5 RCP 4.5 and 8.5 simulations of respiration, we used an r1i1p1 ensemble mean from 15 models (see Supplementary Information).