Note. T and P refer to air temperature and precipitation. The monthly anomalies were calculated for November to March based on periods of 2071-2100 and 1985-2014.
3 Results and discussion
3.1 Sensitivity analyses and model calibration
The Sobol indices showed that the dens_min (minimum snow density) and compaction_rate (the rate at which snow is compacted over time), and their interactions, mainly influence the modelled seasonal snow depth, snow density, snow temperature, and GT (Figure D1-3).
The parameter values yielding the lowest measured-modelled differences in seasonal snow depth and GT are 50, 60, and 90 for dens_min, and 0.5, 0.8, and 0.9, for compaction_rate, at the birch forest, peat plateau, and fen sites, respectively (Figure D4-8).
For the tundra site, we applied parameter values of 90 and 0.8 for minimum snow density and compaction rate, as these showed the best agreement with the available growing-season observational data.
3.2 Evaluation of physical and biogeochemical variables in the historical period

3.2.1 Evaluation of physical variables

The seasonal patterns of snow cover were captured in the model at all sites. However, snowpack depth was considerably underestimated at the birch forest and the low-elevated fen site, but overestimated in low vegetation settings such as the peat plateau (Figure E1), likely due to the model lacking lateral transport and trapping/deposition of wind-blown snow. This mismatch affected the snow insulating capacity and caused substantial cold and warm biases in winter GT in the birch forest and peat plateau, respectively, but not in the fen (Figure E2) where the modelled snowpack depth exceeds the depth of the maximum insulating effect (40 cm; Zhang, 2005). The growing season GT was well captured at the birch forest and fen sites, but was underestimated in the peat plateau (Figure E2a,c,d).
The model was able to simulate short-term fluctuations in snow depth and GT during WWEs at all sites, but not their magnitudes (Figure E3). At the birch forest and peat plateau, larger than measured fluctuations in GT were modeled during WWEs and these differences were much smaller in the late winters. The modelled differences were largely linked to the simulated insulating capacity of snowpacks. Noticeably, stronger modeled reductions in snow depth occurred during early- and mid-winter WWEs, when the energy needed to warm-up and subsequently melt the thin and fresh (less dense) snowpacks is smaller. The modeled and measured GTs at the fen site remained around 0 °C throughout the winter due to the strong insulation of the thick snowpack, but the model captured the observed snowpack responses to WWEs.
3.2.2 Evaluation of biogeochemical variables
The modeled maximum LAI of 1.6 and 1.5 fell within the observed ranges at the birch forest (Heliasz, 2012) and tundra sites (Simin et al., 2021). Following the biases in GT, winter CO2heterotrophic respiration (Rh) was underestimated at the birch forest and overestimated at the peat plateau (Figure E4). We speculate an overestimation of winter Rh at the tundra –it presents similar winter Rh values than the fen and peat plateau, and 6-fold larger than the birch forest– due to the overestimated soil C pool at this site. The model underestimated the growing-season ecosystem respiration (Reco) and gross primary production (GPP) in the tundra (Figure E5c,d). At the birch forest, the model predicted a strong annual C sink in 2007-2010, but measurements indicated a fairly neutral C balance (Figure E4a,b). At the peat plateau and fen sites, the model captured the annual fluctuations and magnitudes of CO2 fluxes (Figure E4c-f). Across four sites the model was unable to capture the observed strong C source in September.
The differences in CH4 fluxes between the peat plateau and fen sites were well captured by the model. However, the model underestimated by 55% the winter CH4 fluxes at the peat plateau, likely due to the colder bias in modeled GT (Figure E6a,b).
3.3 Effects of enhanced WWEs
3.3.1 Impacts on physical variables
The WWE experiments (S1-S3) caused an overall reduction in winter GT (up to 2 °C) in all ecosystem types, except for the runs driven by WWEs from CanESM5 SSP585 at the tundra and the peat plateau (Figure 1a-d). The modeled snow depth also decreased substantially under all WWE experiments. This might suggest that the major driver of the modelled WWE effects on winter GT is snow insulation, which decreases as a combination of the snowpack depth reduction, and the increased thermal diffusivity (caused by higher thermal conductivity and lower heat capacity after freeze-thaw processes), facilitating the heat exchange between atmosphere and soil (Figure 3). However, a “tipping point” may be crossed above a certain WWE magnitude, when the strong cooling effects of the reduced snow insulating capacity are counteracted by the stronger warming effects of longer and more extreme WWEs. This point is reached faster under shallow snowpacks where snow insulation is weak and its further reduction has a smaller effect on GT compared to the overall warmer conditions.
In contrast to WWEs, altered future winter climatologies (S4) increased the modeled winter GT at all sites, from less than 1 °C to as much as 4 °C at the birch forest, tundra, and peat plateau in the warmest scenario, i.e., CanESM5 SSP585, despite snowpack reductions of >80%. The modeled GT warming effect diminishes under thick snowpacks, due to the snow insulation-related process described above: at the fen site (thickest snowpack) we can only see a GT warming under CanESM5 SSP585.
Our results agree with the modelling results by Beer et al., (2018), who suggested that WWEs may cause GT cooling mainly by reducing the snowpack depth. However, there is increasing observational evidence that intense ROS events cause substantial and long-lasting GT warming in winter (e.g., Hansen et al., 2014) through the latent heat released from refreezing of infiltrated water at the bottom of thick snowpacks (Westermann et al., 2011; Pascual et al., 2022). This long-lasting GT warming cannot be captured in LPJ-GUESS due to the lack of essential processes describing the energy and water exchanges between the atmosphere, snowpack, and soil (Figure 3).
In non-winter, the impacts of the manipulation experiments on GT followed those seen for winter GT, although smaller in magnitude (Figure 1e-h). The WWE impacts on non-winter GWC (0-50 cm depth) were marginal except for the fen site (decrease of up to 0.3 m-3m-3)(Figure F1). GWC affects the thermal properties of soils, and its reduction likely contributed to the larger non-winter GT cooling modeled at the fen. The reduced albedo following an earlier snow cover disappearance can contribute to faster GT warming in spring, but this process is not represented in the current version of LPJ-GUESS.