Accounting for Winter Warming Events in the Ecosystem Model LPJ-GUESS Evaluation and Outlook

Winter warming events (WWEs) are short‐lasting events of unusually warm weather, occasionally combined with rainfall, which can cause severe ecosystem impacts by altering ground temperatures and water fluxes. Despite their importance, how large‐scale ecosystem models perform in depicting the impacts of WWEs remain largely unknown. The frequency and intensity of WWEs will likely increase further in the future, making it necessary to understand their potential impacts on high‐latitude ecosystems. In this study, we evaluated the ability of the dynamic ecosystem model Lund‐Potsdam‐Jena General Ecosystem Simulator (LPJ‐GUESS) to represent the responses of subarctic ecosystems to future WWEs, and identified model gaps hindering more accurate estimates of these responses. In response to WWEs, the model simulated substantial ground cooling (up to 2°C in winter) due to reduced snow depth (insulation), with rain on snow (ROS) exerting a marginal influence on the ground temperature responses: these modeled responses are in apparent contradiction with the strong ground warming effect of ROS reported in most observational studies. The simulated ground cooling led to changes in biogeochemical fluxes that were substantial and often comparable in magnitude (but often opposite in direction) to those from altered winter climatologies. The mismatch between the modeled and the observed ground temperature responses to WWEs highlights LPJ‐GUESS's current limitations in realistically simulating some of the effects of WWEs. These limitations likely stem from the (a) absence of a surface energy balance, (b) lack of snow‐vegetation interactions, (c) daily time‐step, and (d) simplistic water retention scheme in LPJ‐GUESS.


Introduction
The Arctic is warming three times faster than the global average, with the strongest warming occurring in autumn and winter (Amap, 2021).Consequently, the occurrence and impacts of winter warming events (WWEs), that is, short-lasting extraordinarily warm spells, often accompanied by rainfall (rain on snow; ROS), are increasing rapidly and expanding geographically (e.g., Bartsch et al., 2010;Pan et al., 2018;Vikhamar-Schuler et al., 2016).These short-lived but intense events can cause significant environmental and societal impacts that can override the effects of long-term climatic trends (e.g., Phoenix & Bjerke, 2016).Climate models predict a further increase in the frequency and intensity of WWEs, making their study a research priority (Pascual et al., 2020).
WWEs can affect ground temperatures (GT) in multiple ways, mostly through: (a) direct heat transfer from the air, (b) latent heat release from refreezing melt and rainwater, and (c) changing the snowpack properties, such as depth and density, which influence the energy exchanges between the atmosphere, snow, and soil.By inducing the melting of the snowpack, WWEs can further influence GT by altering surface albedo and groundwater content, which can have impacts lasting to the growing seasons (Pascual & Johansson, 2022).These WWE-associated environmental changes can strongly alter microbial activity, greenhouse gas emissions (e.g., Natali et al., 2019), and permafrost and vegetation dynamics (Bruhwiler et al., 2021).Additionally, the unusually high air temperatures and the loss of snow insulation during and after WWEs can affect vegetation dynamics through changes in phenology and frost stress (e.g., Bokhorst et al., 2009Bokhorst et al., , 2010)).Together, these processes are likely to alter substantially the arctic carbon cycle.Given the important role of high-latitude ecosystems on the global carbon cycle (Hugelius et al., 2014;Tagesson et al., 2020), understanding the ecosystem impacts of WWEs in these regions is paramount for constraining the future changes in the global carbon budget and designing informed climate change mitigation and adaptation strategies.
Due to the co-ocurrence of multiple and interlinked processes, their complexity, and the lack of highly controlled manipulation studies, the processes and feedbacks related to WWEs are difficult to disentangle in observational data and challenging to implement in models.Recently, substantial efforts have been put in developing more realistic modeling schemes to account for some important high-latitude processes in future predictions, including climate-snow-soil interaction (Pongracz et al., 2021), vegetation cold hardening (Lambert et al., 2022), and rain water infiltration in the snow (Westermann et al., 2011).Yet, current large-scale ecosystem models might still suffer from these and other limitations in representing the impacts of stochastic climate extremes, including for example, coarse spatio-temporal resolutions, and lacking/simplistic hydraulic-schemes and snow-vegetation interactions (e.g., Tang et al., 2015), making it necessary to evaluate their performance and identify the most important gaps.
In this study, we aim to address this research gap by evaluating, for the first time, the ability of the dynamic ecosystem model LPJ-GUESS to represent the responses of subarctic ecosystems to future WWEs.Specifically, we focus on the modeled responses to predicted WWE scenarios of physical and biogeochemical variables during winter and the growing-season.Additionally, we aim to identify model gaps hindering more accurate representations of the ecosystem responses to WWEs.This work will likely provide key information about how to enhance the accuracy of model predictions of ecosystem changes in high-latitude regions and contribute to the ongoing efforts to mitigate the impact of climate change in the Arctic.

Study Sites
The Torneträsk area, in northern Sweden, is a topographically heterogeneous area that aligns along a strong westeast oceanic-continental climatic gradient, with precipitation and winter temperature decreasing eastwards due to the increasing continentality and the rain shadow effect caused by the Scandes Mountains.The area has experienced rapid climate warming (Callaghan et al., 2010) and a substantial increase in the frequency and intensity of WWEs (Pascual & Johansson, 2022).
Vegetation in the area varies following its climatic and altitudinal gradients.Birch forests occur below c. 600 and 800 m above sea level (m.a.s.l) (Van Bogaert, 2010).Tundra species dominate above the tree line, while in the lowlands, birch forests alternate with peat plateaus underlain by permafrost, and non-permafrost fens.
In this modeling study, we selected four sites representing these dominant ecosystem types in the Torneträsk area, including (a) a birch forest (∼370 m.a.s.l) located <10 km east of the Abisko Station (ANS) (Heliasz, 2012); (b) a tundra site (∼410 m.a.s.l), located <1 km to the southeast of the ANS (Michelsen et al., 2012 and references therein); (c) a peat plateau (∼380 m.a.s.l) known as Storflaket, located c. 6 km east of the ANS (Johansson et al., 2013), and (d) a fen (∼515 m.a.s.l) located c. 25 km west of the ANS, near the Katterjokk Station (Swedish Meteorological and Hydrological Institute (SMHI).The dominant vegetation species at these sites are found in Appendix A in Supporting Information S1.

Model Description
We used the Lund-Potsdam-Jena General Ecosystem Simulator (LPJ-GUESS) for our study.This process-based dynamic ecosystem model is widely used in regional and global scale studies (Smith et al., 2001(Smith et al., , 2014)).The model simulates vegetation dynamics (including vegetation establishment, mortality and competition, etc.), water, carbon and nitrogen cycles, and soil biogeochemistry.We used the latest version of LPJ-GUESS (version 4.1, Smith et al., 2014), which includes a recently developed dynamic, intermediate complexity snow scheme enabling the simulation of climate-snow-soil interaction (Pongracz et al., 2021).The model can simulate up to five snow layers, their physical and thermal properties, and their development throughout the cold season.Based on the individual snow layer properties (e.g., temperature, density, thermal conductivity), freeze-thaw processes in snow layers, and heat transport through the snowpack between the atmosphere and soil can be simulated, and ROS events can be accounted for.LPJ-GUESS also includes detailed representations of permafrost and wetland processes, including peatland hydrology, peatland-specific PFTs, and CH4 emissions (see Wania et al., 2009aWania et al., , 2009bWania et al., , 2010)).In the current model version, snow-vegetation interactions are not implemented, including for example, snow interception of vegetation, mechanical and physiological damage due to lack of snow insulation, and cold hardening.

Model Setup
Daily climate data provided by the ANS (ANS, 2020) and Katterjokk Station (SMHI), including air temperature, air temperature daily range, and precipitation, together with shortwave radiation (1913( -1984( , Sheffield et al., 2006;;1984-2018, ANS) , ANS) and annual CO 2 concentrations obtained from the Global Monitoring Laboratory (https://gml.noaa.gov/ccgg/trends/),were used to drive the model from 1913 to 2018.Given their vicinity and similar elevation (altitudinal range <100 m), the birch forest, peat plateau, and tundra sites were run with climate data from the ANS data set (1913-2018;ANS, 2020), whereas the fen site used data from Katterjokk Station (1973-2018;SMHI) and bias-corrected daily data  from the ANS.Soil property data was extracted from the WISE5min, V1.2 Soil Property Database (Batjes, 2012).More details about the setup and input data are found in Appendix B in Supporting Information S1.
Representative plant functional types (PFTs) were selected for each site (Table A1 in Supporting Information S1).We enabled high-latitude and wetland-specific plant functional types (PFTs) in the simulations to better capture site-specific vegetation conditions (see Wania et al., 2009 for more details).The PFT parameters at each site followed previous studies (e.g., Gustafson et al., 2021;Tang et al., 2015) (Table A2 in Supporting Information S1).

Model Calibration and Evaluation
Sobol sensitivity analyses (SA) (Saltelli, 2002;Saltelli & Annoni, 2010;Saltelli et al., 2008) were conducted to explore the influence of different snow-related parameters and parameter interactions on the estimated seasonal snow density, snow depth, snow temperature, and GT at the study sites (except for the tundra site due to lack of observational data for calibration and evaluation) (Appendix C in Supporting Information S1).A sampling of eight relevant parameters was conducted, using ranges based on literature values when available, and certain percentage variations from the original values (Table C1 in Supporting Information S1).Among the most influential parameters, we selected the parameter values that minimized the absolute differences between the measured and the modeled seasonal snow depth and GT at each site (2006-2012 at the peat plateau, 2001-2010 elsewhere) (Figures C4-C8 in Supporting Information S1).The model was subsequently evaluated with independent observational data (2011-2018) whenever possible (Appendix D in Supporting Information S1).

Model Simulations With Future WWEs
We conducted manipulation experiments to assess the model's ability to simulate the responses of different ecosystem processes to WWEs under various levels of frequency and intensity.We selected six different climate scenarios in the Coupled Model Intercomparison Project Phase 6 (CMIP6) (Eyring et al., 2016) from two general circulation models (GCMs) representing the high and low ends of climate sensitivities (CanESM5 and GFDL-ESM4) and three shared socioeconomic pathways (SSP119, SSP270, and SSP585) representing different levels of greenhouse gas projections.For each scenario (n = 6), we extracted daily meteorological data  for the gridcell near the Torneträsk area and bias-corrected it at a daily scale against the observed meteorological data using the period 1985-2014, following Hawkins et al. (2013).Detailed descriptions of the CMIP6 scenarios and the bias-correction method are found in Appendix E in Supporting Information S1.
The bias-corrected GCM's outputs were used to create three manipulation experiments in addition to the HIS-TORICAL runs (S0, daily observational data 1913-2018) (Table 1), In these experiments, the future monthly anomalies in the frequency and intensity of melt days (S1), ROS (S2), and both (S3), in the GCM's outputs (Table E1 in Supporting Information S1) were added to the 1985-2018 period within the HISTORICAL climate inputs, maintaining the long-term climate means as unchanged as possible.These anomalies were calculated based on four indices modified from Vikhamar-Schuler et al. ( 2016) (Table 1) and refer to the differences in melt and ROS days counts (frequency), and the positive degree days and rain-on-snow amounts (intensity), in the GCM's data for the periods 2071-2100 and 1985-2014.Anomalies were added as follows: if anomalies in melt days and positive degree days for a given month and scenario were for example, 2 days and 10°C respectively, we artificially increased the temperature of 2 non-melt days every such month to 0.1°C (those days with air temperature closest to 0°C), and added the remaining 9.8°C in equal amounts among the melt days occurring in such month, each year for 1985-2018.Likewise, if anomalies in ROS days and rain amounts for a given month and scenario were 2 days and 10 mm/month respectively, we generated 2 new ROS days every such month by artificially increasing the daily precipitation of dry melt days (if occurring, starting from the warmest) to 0.1 mm, and then adding 9.8 mm of rain equally distributed among the ROS days occurring in such month each year (for 1985-2018).If no melt days occurred in a given month, they were created beforehand as for scenario S1.To compare the effects of altered WWEs with those of altered long-term climate trends, we conducted an additional experiment (S4) in which we directly added future winter monthly anomalies (Table E2 in Supporting Information S1) to the historical daily air temperature and precipitation.

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 modeled seasonal snow depth, snow density, snow temperature, and GT (Figures C1-C3 in Supporting Information S1).
The parameter values yielding the lowest measured-modeled 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 (Figures C4-C8 in Supporting Information S1).
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.

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 D1 in Supporting Information S1), 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 D2 in Supporting Information S1) where the modeled 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 (Figures D2a, D2c, D2d in Supporting Information S1).
The model was able to simulate short-term fluctuations in snow depth and GT during WWEs at all sites, but not their magnitudes (Figure D3 in Supporting Information S1).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 modeled differences were largely linked to the simulated insulating capacity of the snowpack.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.

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 CO 2 heterotrophic respiration (R h ) was underestimated at the birch forest and overestimated at the peat plateau (Figure D4 in Supporting Information S1).
We speculate an overestimation of winter R h at the tundra -its similar winter R h values to the fen and peat plateau, and 6-fold larger values than the birch forest clash with observations across the Arctic (Natali et al., 2019)-likely due to the overestimated soil C pool at this site.The model underestimated the growing-season ecosystem respiration (R eco ) and gross primary production (GPP) in the tundra (Figures D5c and D5d in Supporting Information S1).At the birch forest, the model predicted a strong annual C sink in 2007-2010, but measurements indicated a neutral C balance (Figures D4a and D4b) due to the underestimated winter Rh (driven by the underestimated winter GT).At the peat plateau and fen sites, the model captured the annual fluctuations and magnitudes of CO 2 fluxes (Figures D4c-D4f).Across the four sites, the model was unable to capture the observed strong C source in September.
The differences in CH 4 fluxes between the peat plateau and fen sites were well captured by the model.However, the model underestimated by 55% the winter CH 4 fluxes at the peat plateau, likely due to the colder bias in modeled GT (Figures D6a and D6b in Supporting Information S1).

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 (Figures 1a-1d).
The modeled snow depth also decreased substantially under all WWE experiments.This might suggest that the major driver of the modeled 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 threshold 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 stronger 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, that is, 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 modeling results by Beer et al. (2018), who suggested that WWEs might cause GT cooling mainly by reducing the snowpack depth.However, there is increasing observational evidence that intense Journal of Geophysical Research: Biogeosciences 10.1029/2023JG007464 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 (Pascual & Johansson, 2022;Westermann et al., 2011).This long-lasting GT warming cannot be captured in LPJ-GUESS, likely 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 (Figures 1e-1h).The WWE impacts on non-winter groundwater content (0-50 cm depth) were marginal except for the fen site (decrease of up to 0.3 m 3 m 3 ) (Figure F1 in Supporting Information S1).
Groundwater content affects the thermal properties of soils, and its reduction likely contributed to the larger nonwinter GT cooling modeled at the fen.The reduced albedo following an earlier snow cover disappearance can contribute to faster GT warming in spring, but albedo is not represented in the current version of LPJ-GUESS.

Impacts on Biogeochemical Variables
Our results further showed large WWE-induced impacts on all biogeochemical variables, with magnitudes often comparable to impacts deriving from altered winter climatologies (Figure 2).GT and water availability are key drivers of the ecosystem C cycle, as they influence the start of the growing season, nutrient availability, vegetation dynamics, and soil R h .
The modelled impacts of the manipulation experiments on biogeochemical variables generally followed the same direction as those observed for GT and groundwater content.We noted substantial reductions in GPP (Figures 2q-  2t) under the WWE experiments S1 and S3 of up to 25% at the fen, 20% at the peat plateau, and 10% at the birch site, and smaller reductions under ROS events alone (S2).The tundra site showed weaker reductions in GPP than the other sites, ranging from a few percent to up to 5% in the milder WWE experiments (driven by SSP119), and increases of up to 10% under CanESM5 SSP585.In contrast, for S4, the model simulated weak changes in GPP in the mildest CMIP6 scenarios, but sizable increases in GPP in the CanESM5 SSP585: up to 75% in the peat plateau, c. 50% in the birch forest, c. 25% in the tundra, and c. 10% in the fen site.The impacts on R a mirrored those described for GPP (Figures 2m-2p).Noticeably, the relationship between GT, and both GPP and R a , is not linear because the processes involved in plant dynamics are multiple and complex.For example, observational and manipulation studies have reported considerable vegetation damage, delays in bud phenology, and reductions in vegetation greenness due to phenological and frost stress following WWEs (Bokhorst et al., 2009(Bokhorst et al., , 2010)).Domine et al. (2015) simulated substantial winter ground warming attributable to the snow interception of tallergrowing shrubs that reduces the horizontal redistribution of snow and therefore results in thicker (more insulating) snowpacks.Lambert et al. (2022) reported substantial improvements in the simulation of vegetation dynamics in high-latitudes when incorporating wintertime plant physiological changes (cold hardening) into modelling schemes.These impacts are not explicitly represented in the model, as it does not yet include snow-vegetation related processes such as vegetation interception of snow, cold hardening, and mechanical and physiological damage due to the loss of snow insulation following WWEs.
Increases in GT can stimulate microbial activity and accelerate litter and soil decomposition rates, resulting in increased R h (Natali et al., 2019).Accordingly, winter R h increased (decreased) by <5% at the fen, and up to 25% at the tundra and peat plateau (Figures 2a-2d), following the increases (decreases) in GT.Contrastingly, in response to the CanESM5 SSP585 scenario at the birch forest site, despite the overall lower winter GT, winter Rh increased by up to >200%.This may be explained because the modelled winter R h in the birch forest is very low (6 g C m 2 , compared to >30 g C m 2 at the other sites), and even small increases in R h during WWEs can result in substantial relative increases in winter emissions.A recent synthesis of in situ observations across the Arctic indicates that winter R h accounts for a substantial portion of the arctic's annual C budget (Natali et al., 2019).
Winter R h alone currently offsets ∼40% of the measured annual vegetation C uptake at our sites.Applying the winter R h response curve to GT by Natali et al. (2019) (Q10 = 2.9), WWE-induced impacts on winter GT of the magnitudes reported here could change total winter R h by up to 25%.Hence, realistically simulating the effects of WWEs on GT is important for improving estimates of greenhouse gas emissions in high latitudes.
The modelled CH 4 emissions at the peat plateau and fen sites decreased under all WWE experiments due to WWE-induced decreases in GT and groundwater content (Figure G1 in Supporting Information S1).At the fen site, emissions halved mostly due to lower water tables in the non-winter season.
Overall, the year-round ecosystem C exchange (i.e., NEE) of the birch forest, tundra, and peat plateau sites decreased (i.e., became a smaller C sink) considerably under most WWE experiments , generally between 20% and 50%, and up to 90% occasionally in the tundra site, due to larger reductions in vegetation C assimilation compared to the C losses from R h .Contrastingly, NEE at the fen site increased weakly (i.e., became a stronger C sink), mostly due to large reductions in R h .Noticeably, the changes in NEE caused by WWEs (S1-S3) are substantial and of similar magnitude to those caused by shifts in future winter climatologies (S4).This further indicates that despite their short duration, accounting for WWEs in LPJ-GUESS may alter significantly the modelled changes in high-latitude ecosystem C cycling.
Given the observed discrepancies in modelled versus measured GT responses to WWEs (Section 3.2.1), the modelled impacts on ecosystems C fluxes reported here should not be interpreted as a prediction of future impacts of WWEs, but rather as a sensitivity test of the current model's responses to future WWEs.Additionally, the inability of LPJ-GUESS to reproduce some observed ecosystem responses to WWEs highlight certain model limitations in predicting future changes in high-latitude ecosystems.

Model Limitations and Potential Future Developments
LPJ-GUESS has difficulties in capturing the physical changes in snow depth induced by WWEs.We identified the lack of surface energy balance in LPJ-GUESS as one of the main limitations to simulating WWE impacts on GT.Heat transfer between air-snow-soil layers is simulated using the Crank-Nicholson finite difference scheme (Crank & Nicolson, 1996).This indirect method of defining heat transfer may significantly affect the computed snow layer and ground temperatures, and the rate and magnitude of snowmelt events.Additionally, the model's daily time step is too coarse to capture the sub-daily freeze-thaw cycles and hydrological processes within the snowpack.The absence of a surface energy balance in LPJ-GUESS also implies an omission of the impacts related to changes in albedo including reductions in albedo due to shorter snow cover seasons and the exposure of vegetation above the snowpack following WWE-induced snowmelt.
In the model, heat transfer is affected by the changes in snow layer temperature and thermal properties, and the latent heat release upon freezing is not captured (Figure 3).The SA showed that the current model setup is not sensitive to changes in liquid water holding capacity in snow layers (Figures C1-C3 in Supporting Information S1), potentially due to the simplistic water retention scheme applied.Rainfall infiltration follows a bucket model approach, limited by each snow layer's maximum water-holding capacity.Precipitation from ROS events and melt water is quickly forwarded to the soil as runoff at the simulated time step (i.e., daily) when the maximum liquid holding capacity of a layer is reached.If the ground is frozen, the excess liquid water will not stay in the bottom snow layer, but drain out as surface runoff.This feature also prevents the formation of ice layers of high thermal conductivity within the snow layers, as many observational studies suggest (e.g., Langlois et al., 2017), which could influence the simulated GT.Observations suggest major ROS events may have durable impacts on GT (e.g., Westermann et al., 2011), but the current model setup cannot capture such persistent effects due to the simple water retention scheme and the lacking processes related to energy balance.
The LPJ-GUESS version used in this study has an intermediate complexity snow scheme, similar to many ecosystem models.These modules are developed and tuned to represent average conditions rather than capturing extreme and smaller-scale phenomena such as WWEs.There is a need for either further snow scheme development in LPJ-GUESS or evaluating a more complex, designated snow model to address the listed shortcomings and capture internal snow dynamics on a finer spatio-temporal scale.A recent extension of the model with detailed land surface processes and surface energy balance LPJ-GUESS LSMv1.0 (Martin Belda et al., 2022) could be used in future studies to assess whether the model-measurements mismatch is reduced when using the LPJ-GUESS with detailed energy balance and sub-daily processes.This land surface model version of LPJ-GUESS needs to be thoroughly evaluated and tested in high latitude environments.
Predictions of future changes in high-latitude ecosystems would also benefit substantially from incorportating processes related to snow-vegetation interactions, including snow interception by vegetation, and other lessunderstood processes such as vegetation cold hardening.These processes have been shown to exert significant influences on ground temperature, hydrology, and vegetation dynamics but are currently lacking in LPJ-GUESS and other ecosystem models (e.g., Lambert et al., 2022).Regional modeling studies are required to further understand how WWEs might affect the pan-Arctic carbon budget.

Conclusions
We applied future-simulated WWE experiments in LPJ-GUESS to evaluate the model's ability to represent the impacts of these events on four dominant ecosystems in subarctic Sweden.The modeled ecosystem responses to WWEs were substantial and of magnitudes often comparable to those of altered winter climatologies.These events induced reductions in GT, which altered numerous biogeochemical processes and resulted in substantial changes in ecosystem CO 2 uptake.However, the modeled decreases in GT following WWEs diverged from the in situ observations and the majority of the observation-based literature, which highlights certain limitations of LPJ-GUESS in capturing the effects of WWEs on ecosystem processes, adding uncertainties in future predictions of ecosystem changes.We identified the most important current model limitations to be (a) the lack of surface energy balance, (b) the absence of snow-vegetation interactions, (c) the model's daily time step, (d) and the simplistic water retention scheme applied in the LPJ-GUESS snow scheme.

Figure 3 .
Figure 3. Theoretical model of the winter warming events (WWE) impacts.Box colors indicate the overall direction of changes (red = increase, blue = decrease) in the variables based on the literature.Thick outlines and bold text indicate a disagreement between the literature and our simulations.Arrow colors show the direction of the change exerted by each process.Dashed lines refer to ROS-related processes.Gray boxes and lines refer to uncertain processes and responses.Green ticks and red crosses indicate WWE-related processes that are, or are not, implicitly implemented in the current model version.

Table 1
Description of the HISTORICAL and MANIPULATION RunsNote.T and P refer to air temperature and precipitation.The monthly anomalies were calculated for November to March.