Decadal increases in carbon uptake offset by respiratory losses across northern permafrost ecosystems

Tundra and boreal ecosystems encompass the northern circumpolar permafrost region and are experiencing rapid environmental change with important implications for the global carbon (C) budget. We analysed multi-decadal time series containing 302 annual estimates of carbon dioxide (CO2) flux across 70 permafrost and non-permafrost ecosystems, and 672 estimates of summer CO2 flux across 181 ecosystems. We find an increase in the annual CO2 sink across non-permafrost ecosystems but not permafrost ecosystems, despite similar increases in summer uptake. Thus, recent non-growing-season CO2 losses have substantially impacted the CO2 balance of permafrost ecosystems. Furthermore, analysis of interannual variability reveals warmer summers amplify the C cycle (increase productivity and respiration) at putatively nitrogen-limited sites and at sites less reliant on summer precipitation for water use. Our findings suggest that water and nutrient availability will be important predictors of the C-cycle response of these ecosystems to future warming.

Because the number of sites collecting data has increased since the beginning of our time series, we reran our models on a subset of the data to assess the robustness of our findings to temporal biases in data collection.We re-ran the models on the subset of observations from 2003-present, representing the most recent 20 years of data.This represents 85% and 81% of summer and annual NEE estimates, 87% and 74% of summer and annual GPP, and 85% and 81% of summer and annual Reco estimates, respectively.This analysis supported our core conclusion that annual CO2 uptake is increasing across non-permafrost ecosystems (P = 0.02), but not across permafrost ecosystems (P = 0.54), with estimates of both slopes within 1 g m-2 yr-1 of estimates based on the full dataset (Fig. 2 in main text).Long-term trends in summer NEE were also largely similar when using the shortened timeseries, although the nonpermafrost ecosystem C sink increased at a slightly faster rate (-4.1 g m-2 yr-1; P = 0.002 vs. -2.6 in full dataset), and at a slightly slower rate in permafrost ecosystems (-1.9 g m-2 yr-1; P = 0.13 vs. -3.0 in full dataset), with 95% confidence intervals in the slopes of both ecosystems largely overlapping between the subset and full datasets.Analysis of changes in the component fluxes across this shorter time period supported the idea of increased amplification rates of permafrost C cycling, with both GPP and Reco generally increasing.Conversely, in non-permafrost sites (where we detected no significant trends in GPP and Reco in the full dataset), this shortened window of observations led to no consistent trends in these component fluxes, with models often failing to converge.To assess the effect of broad-scale spatial patterns in data collection on our results, we subset our data by landmass, and re-ran the models for North America (including Greenland) and Eurasia.North America (largely Alaska) represents 82% of our annual observations of permafrost NEE (Table S2), and here, we found that the trend of increasing annual NEE (i.e.weakening land CO2 sink over time; Fig. 2b main text) became statistically significant.Similarly, evidence was stronger for increasing annual R eco across North American permafrost ecosystems (P < 0.05).Notably, the changes in both GPP and Reco across Eurasian permafrost ecosystems were roughly double those observed in North America, perhaps suggesting a greater amplification of annual C cycling in Eurasia, but the very limited sample size of annual permafrost fluxes here (Table S2) strongly limits interpretation of these results.While trends in permafrost sites may be driven by North American ecosystems, the trends we observed across nonpermafrost sites appeared to be stronger in Eurasia.The increased annual CO2 uptake (decreasing NEE) in non-permafrost ecosystems was statistically significant in Eurasia (P = 0.04), where we found evidence for increased annual GPP over time (P = 0.06), but no clear trend in annual Reco.In contrast, we found little evidence for consistent trends in NEE, GPP, or Reco in North American non-permafrost ecosystems.To further assess the uncertainty in the slope estimates of our timeseries models, we refit all models (with an identical random effects and autocorrelation structure) in a Bayesian framework using the brms package in R. We assigned the fixed effect parameters from the frequentist model outputs as prior distributions, with relatively uninformative priors for the random effects and covariance parameters (Student's T (ν=3, μ=0, σ=47) and LKJ η=1, respectively).Each timeseries model was run on 8 separate MCMC chains of 20,000 steps with a warmup of 12,000, resulting in 64,000 total samples in the posterior distribution.We then compared the slopes and 95% credible intervals from this distribution to the slopes and confidence intervals derived from the original frequentist models.We calculated P-values for the Bayesian slopes as the density value of posterior distribution at the null (0) divided by the density at the Maximum A Posteriori value for the slope (Mills 2019).Generally, the frequentist and Bayesian analyses agreed, with significant (or marginally significant) slopes based on one approach corresponding to significant slopes using the other approach.The only exception to this was the trend towards increasing annual GPP in non-permafrost systems, which was not statistically significant in our original analysis (slope = -6.7 g m-2 yr-1, P=0.25), but had a steeper slope that differed from zero based on the Bayesian analysis (slope = -10.2,P=0.004).
Because our analyses included a relatively small number of sites, the inclusion (or exclusion) of single sites may affect the conclusions of our models.We assessed the robustness of our findings using leaveone-out cross validation (LOOCV) based on cluster-wise exclusion using the cv package in R85.Briefly, for each timeseries model we iteratively removed one site at a time and recalculated the model to assess its predictive ability for the observations in the removed cluster86.We then compared the root mean squared error (RMSE, calculated using fixed effects) based on the predictions to the RMSE for the original model.This analysis suggested that our conclusions regarding NEE trends were likely robust to biases in site selection: The predictive RMSE based on LOOCV was lower than the model based RMSE for the models where we detected a significant trend (permafrost and non-permafrost summer NEE, annual non-permafrost NEE; Table S5).In contrast, GPP and Reco trends appeared to be more sensitive to removing sites.It was not possible to calculate the LOOCV RMSE for many of these models, because running the models on one or more of the subsets of data resulted in models failing to converge.In cases when GPP and Reco models converged for all subsets of data, the predictive RMSE was generally higher than that of the original model.The one exception to this was the model showing an increase in summer GPP in permafrost sites, where the predictive RMSE based on LOOCV was slightly smaller (1%, Table S5) than the original model.Collectively, these uncertainty analyses point to a greater degree of confidence in NEE timeseries than in GPP or Reco.This is perhaps not surprising, given that we have a larger sample size for NEE (Table S2).Additionally, is the variable measured NEE at these sites (using eddy covariance towers or chambers), while GPP and Reco are estimated by partitioning the NEE flux, which introduces additional uncertainty into these estimates.

Table S3 .
Summary of timeseries analyses results from models run on datasets subset from 2003present.

Table S4 .
Summary of annual timeseries analyses results with datasets subset by landmass (i.e.North America and Greenland vs. Eurasia).

Table S5 .
Comparison of slopes and 95% confidence intervals derived from frequentist models (reported in main text), with slopes and 95% credible intervals derived from Bayesian reanalysis.Model root mean square errors (RMSE) for frequentist models are shown in comparison to RMSEs recalculated using cluster-level leave-one-out cross validation (LOOCV).The differences in between model RMSE and crossvalidated RMSE are reported in the last column in units of percent change.