Introduction

The amplified warming of the Arctic-Boreal region strongly affects the thickness and duration of snow cover. Observations show that the start of the snow season is delayed, the snow insulation maximum is peaking sooner, and snowmelt is occurring earlier in the year1,2. While the average long-term trend across the globe is towards less snow, there are important regional differences in both snow depth and snow season length3,4. Notably, where temperatures remain below 0 °C, an increase in atmospheric moisture can still lead to thicker snow5. It remains challenging to predict such regional patterns in future snow conditions due to the high spatio-temporal variability in the warming rate and the timing, amount, and form of precipitation5,6,7. Still, model projections from the Coupled Model Intercomparison Project Phase 6 (CMIP6) suggest two contrasting snow trends that will dominate across the northern hemisphere in the future: (i) either an increase in extreme snow depth—exceeding historic 99.9th percentile snow events—e.g. in North-Eastern Siberia or (ii) a decrease in snow depth e.g. in Europe and Western Alaska8,9,10,11. These differences in snowfall depend on both short-term and local scale weather conditions and long-term regional-scale climate variability12. For example, precipitation increases in the Arctic have been linked to more evaporation due to sea-ice retreat, and northward moisture transport13,14,15. On top of those large-scale drivers, regional differences may occur due to shifts in the Arctic Oscillation16.

Subsequent changes in snow conditions have a large impact on both permafrost thaw and greenhouse gas exchange in the Arctic-Boreal region since snow is one of the key controlling factors over soil thermodynamics and biogeochemistry17. Increases in snow depth have caused soil warming and permafrost degradation in the discontinuous permafrost region, despite unchanged air temperatures1. On average, permafrost temperatures have increased by 0.4–0.6 °C decade−1—a trend that is likely to continue1,18. Moreover, wintertime carbon emissions from the permafrost region are larger than previously thought17, and this carbon loss could, in the future, be amplified even further since the cold season is warming faster (0.7 °C decade−1) than the snow-free season (0.4 °C decade−1)2. Soils in the permafrost region contain nearly half of the total global soil carbon pool (1400–1600 Pg C), and the release of only a fraction would substantially increase the concentration of greenhouse gases in the atmosphere19,20. Previously frozen carbon may become bioavailable through a deepening of the active layer depth (ALD), and, combined with increases in soil temperature, this stimulates heterotrophic respiration (Rh) rates. Through this series of events, changes in snow cover may lead to increased carbon loss from the Arctic-Boreal region.

Nonetheless, models often oversimplify cold season processes, such as snow- and permafrost dynamics, as well as sub-zero soil organic carbon (SOC) decomposition rates, despite their importance for the carbon and hydrological cycles17,21. Most importantly, changing snow cover has multiple impacts on soil temperature and moisture: earlier snowmelt and a declining areal extent of snow cover directly reduce surface albedo, causing a warming of the ground surface. In addition, changes in the snowpack’s properties, such as snow density, influence snow’s insulation capacity, leading to either a cooling or a warming of the soil.

There are still large uncertainties regarding the present and future state of the Arctic-Boreal carbon balance. In particular, whether the Arctic-Boreal region acts as a net sink or source of carbon - now and in the future22,23,24,25. Process-based models indicate a net carbon loss of 120  ± 80 Pg C by 2100 under high emission scenarios6,26,27, but there is also a potential for an increased uptake by vegetation6,27,28. However, it remains uncertain to what extent elevated plant carbon uptake (NPP) might offset the permafrost carbon feedback29. It is, therefore, imperative to better understand how snow-soil-vegetation interactions control this balance and to identify whether future responses diverge across the high latitudes.

To our knowledge, no prior study has examined how regional changes in snow depth and snow season length may affect the Arctic-Boreal carbon balance in the future, with CMIP6 forcing. In this study, we explore the impact of snow cover changes on soil biogeochemistry under different climate scenarios, using an arctic-specific version of the LPJ-GUESS Dynamic Global Vegetation Model, which recently incorporated a dynamic, multi-layer snow scheme21. We aim to (i) identify areas with contrasting changes in snow conditions (snow depth and snow cover duration) and (ii) quantify their influence on carbon pools and fluxes. We forced the model offline with daily atmospheric output from the MRI-ESM2-0 Earth System Model (ESM) from the CMIP6 archive (see Methods). We aggregated the outputs into four spatial groups defined by present-day (2000–2015) mean annual air temperature (Fig. 1), and analysed the model outputs by calculating the differences between historical (1995–2015) and future (2080–2100) conditions.

Fig. 1: Four sub-regions based on present-day (2000–2015) mean annual air temperature.
figure 1

The current permafrost extent estimated by ref. 66 is shown by the grey outlined area, which approximately corresponds to the extent of climatic groups 1 and 2. The spatial groups align well with the group-wise pattern of snow depth and near-surface soil temperature changes between the historical period and future scenarios (Fig. 3).

Full size image

Results

Snow cover duration—defined as the annual cumulative number of days with at least 5 cm of snow on the ground—decreased quasi-uniformly across all regions and climate scenarios (Fig. 2b, d and Fig. S1, see Supplementary Results 1). We found that the annual median snow depth decreased by 3–5 cm by 2100 under all future scenarios (Fig. 2a). However, the change in snow depth was not uniform across the study region (Fig. 2c). Annual mean snow depth increased in Siberia and Northern Canada, whereas snow depth decreased across Europe, Western Russia, southern Alaska, and southern Canada (Fig. 3). This spatial pattern of contrasting snow depth trends aligns well with our four climatic subgroups (compare Figs 1 and 3).

Fig. 2: Future changes in snow depth and snow cover duration.
figure 2

Median annual snow depth (a) and snow cover duration (b), with a 10-year running average. The black dashed line in a, b shows the simulated snow depth and snow cover duration using the CRU-NCEP climate forcing. Spatially grouped relative change (%) in annual median snow depth (c) and snow cover duration (d). Differences in c, d were computed by subtracting the historical period (1995–2015) from the future scenarios (2080–2100). Colours indicate the four climatic regions and the dotted horizontal line represents unchanged conditions between the present and future scenarios. The shaded boxes show the interquartile range, and the horizontal line is the median data value. The error bars show the minimum and maximum range of distributions.

Full size image
Fig. 3: Diverging spatial trends in snow depth across the Arctic region.
figure 3

Diverging trends in mean snow depth are shown as the difference between 2000 and 2015 and 2080–2100 under SSP1-2.6 (a), SSP3.7-0 (b), and SSP5-8.5 (c) scenarios.

Full size image

Median snow depth increased in the coldest region (group 1, Fig. 2c), while snow depth in the two warmest regions decreased. These changes in snow conditions lead to large differences between climate groups in the insulation capacity during the autumn and winter (Fig. S2k, n). The colder the region, the stronger soil temperatures increased—under all climate scenarios. Over 5–15 °C of warming occurred in the coldest region (group 1) during the winter (December-January-February, Fig. S2f), while simulated ALD more than doubled for the SSP3-7.0 and SSP5-8.5 scenarios (Fig. S3). Changes in ALD can be attributed to an increase in soil temperature throughout the year, but the strongest warming occurs during the winter (DJF, Fig. S2f–j).

Changes in snow and soil temperatures led to an increase in both NPP and Rh (Fig. 4a, b and Figs. S4a–j, S5). Absolute changes in NPP and Rh are small in autumn and winter due to low biological activity. A comparison to the Northern Circumpolar Soil Carbon Database version 2 (NCSCDv2) SOC dataset19 revealed that simulated soil carbon in the two colder regions (groups 1 and 2) is significantly underestimated by the model (Fig. S6), with a mean difference of 19.0 and 11 kg m−2 (−78 and −44%), for the respective groups. We also ascertained that simulated SOC is more underestimated for upland soils than peatland soils (Fig. S7). Modelled Rh is a function of soil carbon content, and consequently, our results underestimate the potential carbon release from the coldest regions. Therefore, we normalised Rh (divided by the SOC), to better show the response, given an equal amount of SOC. This shows that for the same amount of SOC, Rh increased the most in the coldest regions and under the warmest scenarios (Fig. 4b).

Fig. 4: Future changes in carbon fluxes per spatial group.
figure 4

Annual differences in simulated NPP (a), normalised Rh (b), and NEE (c). Changes were computed by subtracting the historical period from future scenarios. Colours show the four sub-regions. The shaded boxes show the interquartile range, and the horizontal line is the median data value. The error bars show the minimum and maximum range of distributions.

Full size image

Overall, the model output suggests that the ecosystem’s carbon uptake capacity (NEE) will increase in the future, especially for groups 2 and 3, and nearly no change for group 1 (Fig. 4c and Fig. S4k–o). While the simulated absolute changes in carbon fluxes were small, the coldest region experienced the largest relative changes in carbon fluxes. We observed a strong reduction in relative carbon residence time of more than 50% (Fig. 5). This reduction increased in the warmer scenarios, indicating a faster ecosystem response to changing climatic conditions.

Fig. 5: Annual relative changes (%) between the scenario and historical periods in relative carbon residence time.
figure 5

Colours represent the four sub-regions. The shaded boxes show the interquartile range, and the horizontal line is the median data value. The error bars show the minimum and maximum range of distributions.

Full size image

Discussion

Linking snow conditions and carbon flux changes

Despite non-homogeneous snow trends across the high latitudes, most studies exclusively report mean trends through time11,30. Spatially contrasting snow depth trends are seldom the primary focus of pan-Arctic studies3,10,31. However, the coldest regions (i.e., northern Siberia and northern Canada) are projected to experience increases in snow depth despite a concurrent decrease in snow cover duration. Our simulations show a decrease in snow cover of 14–44 days by 2100, with a more pronounced decrease for SPP3-7.0 and SSP5-8.5 (Fig. 2d). We did not detect spatial contrasts in snow cover duration changes between Eurasia and Northern Canada as reported by others3,6. Rather, the simulated changes were quite uniform across the region (SI Fig. S1 and Table S1).

Differences in how snow cover duration is computed may partially explain these differences. We used a threshold of 5 cm, like other studies31, but previous estimates have also been derived from snow water equivalent and snow-covered fraction in gridcells rather than snow depth10,32. Median changes in snow depth decrease by ~20% by 2100, which is in line with direct output from CMIP6 models30. The diverging snowfall trends above and below the −5 to −10 °C isotherm (groups 1 and 2) are in line with CMIP5 projections33, while increases in simulated snow depths in North-Eastern Russia and East Canada align with the largest increases shown by others31. The spatially divergent changes in snow depth affected the snowpack’s insulation capacity, which translated into spatially diverging impacts on soil temperatures (Fig. S2).

It is challenging to differentiate the impact of air temperature from snow depth changes since air temperature has both a direct (via heat transfer) and indirect (via rain/snow partitioning) effect on soil thermodynamics. However, the emergent changes in soil temperature were not uniform throughout the year. The strongest increase in snow insulation capacity occurred at the end of the cold season (MAM), with the largest increase simulated in the coldest region (Fig. S2l). During summer (JJA), projected changes in air temperature are rather uniform between the groups (Fig. S2c), but changes in soil temperature do show differences (Fig S2h), related to earlier snowmelt. Insulation capacity changes in autumn (Fig. S2n) can be attributed to a later start of the snow season or a complete lack of snow cover. The disconnect between air and soil temperatures, plus the absence of strong differences between the groups during summer, show that the strong soil warming observed in the coldest regions was primarily caused by changes in snow cover rather than air temperature.

The dominant influence of snow cover on soil temperature is strongly supported by snow fence experiments34,35, and a detailed evaluation of the snow scheme in the same arctic-specific version of LPJ-GUESS as in this study21. Based on that evaluation, we are confident that the model realistically replicates snow-soil-carbon flux interactions. This study shows that changes to snow conditions and soil temperatures influence carbon fluxes in a distinctively different way across the northern high latitudes. Perhaps the most compelling finding is that all process rates had the highest relative increase in the coldest group.

The largest absolute changes in carbon fluxes occurred in group 3 (where the mean annual air temperature is between −5 and 0 °C, Figs. 4 and S4), due to the non-linear temperature dependency of biological process rates and higher pre-existing biomass. Even though the largest air and ground temperature increase occurred in the two coldest groups, ground temperatures did not exceed freezing everywhere. This explains why the significant (>10 °C) warming in the coldest regions resulted in small absolute changes. Rh increased strongly in relative terms, but the influence on NEE is small due to the low soil carbon content and long relative carbon residence time. Interestingly, Rh in the coldest groups increases continuously with higher warming scenarios, whereas the increases in the two warmer groups level out (Figs. 4 and S4).

Previous research showed that observed winter respiration is much higher across the permafrost region than ecosystem models suggest17. This implies that the relative contribution of winter-time changes to the permafrost carbon feedback becomes more important with more warming. A key finding of this study is the strong decrease in relative carbon residence time, which is consistent in all applied climate scenarios (Fig. 5). Such decreases in carbon residence time were previously found in Alaska36, and carbon residence time is one of the main contributors to the uncertainty in vegetation response under future scenarios37. Interestingly, we show that observational estimates of mean soil carbon content do not differ strongly among the four climatic regions, while carbon residence time shows a large decrease. This suggests that permafrost carbon loss will be most strongly amplified in the coldest regions with the strongest projected increase in snow depth. This finding should prompt further research on the role of snow cover change on permafrost carbon loss.

Snow-soil-vegetation interactions

Snow events do not only affect soil thermodynamics and carbon emissions but also closely interact with vegetation dynamics, with certain species benefiting from environmental changes, while others are adversely affected38. Observational evidence indicates that the response of arctic vegetation to warming is more complex than the established thinking of a uniform greening trend39,40. The timing of snowmelt is an important factor that sets the stage for the snow-free season41. A shorter snow season yields earlier soil water availability and a faster warm-up of the soil, which could influence the seasonality of carbon fluxes. We show that shorter snow seasons lead to increased simulated future plant productivity (Fig. 4 and Fig. S5).

These findings are supported by previous studies linking earlier snowmelt to a longer growing season and increased total ecosystem productivity, although this does not necessarily lead to a greater net carbon uptake42. A synthesis of 11 sites across the Arctic showed that earlier snowmelt boosts ecosystem respiration, offsetting over 40% of growing season CO2 uptake43. Then again, earlier spring snowmelt has also been linked to increased boreal forest productivity44. Our results show low plant productivity in the two colder groups due to a bioclimatic limitation of vegetation establishment and growth. Further increases in NPP, due to longer growing seasons, may still offset increased Rh from warming permafrost soils.

Limitations

The represented soil-snow-vegetation interactions in the model affect both vegetation dynamics and the modelled carbon fluxes. Limitations regarding the applications include the low simulated plant productivity in the permafrost region and the inability to capture the size of old carbon pools (Figs. S6, S7). Simulated soil carbon pools are, in general, a more important determinant of carbon residence time than vegetation stocks in ecosystem models45. In this study, the calculated residence time is highly dependent on the SOC, and our estimates of these carbon pools are below observational estimates in the coldest regions. We observe a higher deviation from observed soil carbon estimates in upland soils (Fig. S7), which is a common problem across CMIP6 Earth system models46. Despite the uncertainties in the absolute carbon pools, the decreasing trend in carbon residence time across all seasons and spatial groups is robust and indicates that carbon stocks may become more responsive to climate change.

Most models are limited due to their generalised and coarse-scale parameterisations. This means that smaller, sub-grid-scale phenomena are not captured, even though arctic landscapes are highly heterogeneous40,47. Present-day ecosystem models do not represent divergent snow patterns at the sub-grid scale, apart from a few notable exceptions48. To achieve better estimates of the role of snow cover on high-latitude ecosystems, we need a stronger focus on local and sub-grid processes - such as thermokarst formation. The projected increases in snow depth may trigger abrupt thaw events, creating newly formed thermokarst features35,49,50. Thermokarst can lead to a large additional carbon loss from permafrost soils, but few models currently simulate these rapid changes51,52,53.

Our model results suggest that the Arctic-Boreal region will remain a carbon sink in the future, but this is rather uncertain for the coldest part of the Arctic due to an underestimation of SOC. Even though we observed an increase in plant productivity in the future, it remains challenging to assess whether increased carbon sequestration capacity may offset the potential for permafrost carbon loss18. We established that the lower simulated carbon loss is due to low SOC, while productivity is likely overestimated due to a lack of some crucial disturbances in the model, such as thermokarst, pest outbreaks, and vegetation damage due to frost droughts and rain-on-snow events. Of these, thermokarst alone may add 40 % on top of the carbon loss that can be expected from gradual thaw54. Therefore, we consider our results to be conservative, and it is likely that increasing snow depths, as a driver of permafrost carbon thaw, will cause more, not less, carbon loss in the future.

Conclusions

Our results show that non-uniform changes in future snow conditions can significantly influence Arctic-Boreal carbon fluxes through altered snow insulation capacity and snow season length. Carbon uptake strengthened during the growing season, while more carbon was released during the cold season—indicating an intensification of carbon cycling. These shifts in carbon fluxes can be partly attributed to climatic changes (warmer air temperature), but are also unequivocally influenced by snow conditions. We highlight the potential for regional amplification of the permafrost carbon feedback through future snow cover changes, with the coldest areas showing the largest relative changes in carbon fluxes.

One of the key findings of this study is the substantial decrease in carbon residence time in the future, which shows a strongly reduced storage capacity and a higher rate of carbon loss in permafrost-underlain areas. We show that accounting for the spatio-temporal variability in snow depth is essential, and a necessity to go beyond simple averaging of trends across the Arctic-Boreal region. Improving our understanding of the response of the permafrost region to warming is paramount to enable more robust projections of climate feedback originating from the changing Arctic.

Methods

Model description

The Lund-Potsdam-Jena General Ecosystem Simulator (LPJ-GUESS) is a process-based dynamic vegetation model55,56 that includes representations of vegetation establishment, competition for resources (light, water, nutrients), and mortality. The model can simulate key high-latitude processes, including soil freeze-thaw dynamics22,28,57. In this study, we use a customised Arctic version of LPJ-GUESS (version 4.1, subversion 11205), which includes a new advanced snow scheme21. The updated multi-layer snow scheme can capture internal snowpack dynamics and simulate a more realistic snow insulation capacity than the previous snow scheme. Snow cover may play an important role in shaping vegetation dynamics. LPJ-GUESS, as an individual-based vegetation model, can capture shifts in dominant vegetation types and eventual shrub expansion. In the model, snow conditions affect vegetation dynamics through the control of soil temperature and timing and amount of meltwater provided. Direct snow-vegetation interactions are not yet explicitly represented in the model and, therefore, are not the main focus of this study. Our study focuses on northern high latitudes (>60°N). The model is applied on a daily time step to simulate processes at a spatial resolution of 0.5 × 0.5°. The model uses 15 standard plant functional types (PFTs) to characterise prevalent vegetation types58, as well as a sub-set of PFTs specifically designed for arctic conditions (see Table S2 in Supplementary Methods 1 and S2.2 in ref. 28). For further details on the model’s structure, see refs. 55,56,57 and references therein.

Simulation protocol

We used the atmospheric output from three model runs by MRI-ESM-2.0 under the shared socioeconomic pathways (SSPs) and representative concentration pathways (RCPs) scenario matrix—SSP1-2.6, SSP3-7.0, and SSP5-8.5—to simulate a wide range of different future scenarios59. CMIP6 ensembles represent historical snow climatologies better, show a smaller winter bias, and have stronger trends in snow extent and snow mass than CMIP5 (see refs. 10,11). We applied a 500-year-long spin-up period, forced with a repeating 30-year de-trended historical static climatology to establish equilibrium vegetation conditions. To establish soil carbon stocks, an offline spin-up was used by saving monthly litter input, decomposition rates, C and N leaching, N deposition and vegetation uptake for years 140 to 220 of the spin-up—repeating these years offline for 40,000 years. Historical atmospheric greenhouse gas concentrations followed ref. 60. Statistically downscaled CMIP6 climate output was bias-corrected following the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) 3b protocol (using ISIMIP3BASD v2.5.061). Climatic forcing variables were surface air temperature (minimum, average and maximum), total precipitation, net incoming shortwave radiation, surface wind speed and relative humidity. Nitrogen deposition and CO2 for each scenario were derived from the CMIP6 Input datasets for Model Intercomparison Projects (input4MIPs62). Wetland extent was prescribed with the Boreal-Arctic Wetland and Lake Dataset (BAWLD; ref. 63).

Data analysis

We analysed changes in snow conditions, soil temperature and carbon fluxes in four distinct spatial groups (Fig. 1) to evaluate potential links between snow conditions and biogeochemical changes. Differences were calculated by subtracting the historical means for 1995–2015 from the mean for 2080–2100 under each climate scenario. We evaluated the annual and seasonal changes to quantify and compare biogeochemical changes. We evaluated the size of the simulated soil organic carbon (SOC) pools by comparing the total soil column estimates (150 cm) to the top 100 cm observational estimates from the Northern Circumpolar Soil Carbon Database version 2 (NCSCDv2.2) dataset64. We show the model-data comparison depending on the defined climatic groups (Fig. S6) and in upland and wetland classified gridcells (Fig. S7) in the Supplementary Information. Relative carbon residence time was calculated by dividing the amount of soil carbon by heterotrophic respiration, as shown in Eq. (1), where SoilC indicates the total soil carbon, Rh is heterotrophic respiration and τ is relative carbon residence time. Data regarding the absolute and relative changes in simulated carbon fluxes can be found in Table S3 in the Supplementary Information.

$$\tau =\frac{SoilC}{Rh}$$
(1)

Study limitations

The chosen ESM, MRI-ESM2-0, has a central-estimate climate sensitivity (ECS) of 3.1 °C65 fitting our application, but despite this robustness, our analysis might benefit from using more ESMs—preferably with different ECSs to provide a potential range of future snow conditions. However, we included in our analysis a wide range of climate scenarios with different amounts of warming, which all show the same direction of change for all four climatic regions. Therefore, we are confident that forcing LPJ-GUESS with the output from an ESM with a different ECS would show the same overall pattern. More information regarding the choice of climatic forcing and selected spatial groups can be found in the Supplementary Information (Supplementary Methods 2, Supplementary Results 2 and Figs. S8–S10). An in-depth comparison of forcing data is outside the scope of this project, but we found that LPJ-GUESS has less vegetation at high latitudes when forced with daily CMIP6 output, compared to the default monthly CRU-NCEP climatic forcing. As a result, the simulated vegetation- and soil-carbon pools were lower than observationally-constrained estimates64. Seasonal differences between CRU-NCEP and WFDE5—the dataset used to bias correct CMIP6 climatic forcing following the ISIMIP3b protocol61—is a potential cause of the deviations between the default LPJ-GUESS and CMIP6 forced simulations in the historical period.