Model simulations of arctic biogeochemistry and permafrost extent are highly sensitive to the implemented snow scheme in LPJ-GUESS

The Arctic is warming rapidly, especially in winter, which is causing large-scale reductions in snow cover. Snow is one of the main controls on soil thermodynamics, and changes in its thickness and extent affect both permafrost thaw and soil biogeochemistry. Since soil respiration during the cold season potentially offsets carbon uptake during the growing season, it is essential to achieve a realistic simulation of the effect of snow cover on soil conditions to more accurately project the direction of arctic carbon–climate feedbacks under continued winter warming.

The Lund–Potsdam–Jena General Ecosystem Simulator (LPJ-GUESS) dynamic vegetation model has used – up until now – a single layer snow scheme, which underestimated the insulation effect of snow, leading to a cold bias in soil temperature. To address this shortcoming, we developed and integrated a dynamic, multi-layer snow scheme in LPJ-GUESS. The new snow scheme performs well in simulating the insulation of snow at hundreds of locations across Russia compared to observations. We show that improving this single physical factor enhanced simulations of permafrost extent compared to an advanced permafrost product, where the overestimation of permafrost cover decreased from 10 % to 5 % using the new snow scheme. Besides soil thermodynamics, the new snow scheme resulted in a doubled winter respiration and an overall higher vegetation carbon content.

This study highlights the importance of a correct representation of snow in ecosystem models to project biogeochemical processes that govern climate feedbacks. The new dynamic snow scheme is an essential improvement in the simulation of cold season processes, which reduces the uncertainty of model projections. These developments contribute to a more realistic simulation of arctic carbon–climate feedbacks.

1 Introduction

The Arctic is undergoing rapid warming, with some of the most pronounced changes occurring during the winter (Box et al., 2019; Natali et al., 2019). As a result, snow thickness, the extent of snow-covered area, and snow season length are decreasing, and this is expected to continue in the future (AMAP, 2017; IPCC, 2014). Snow is an important abiotic component of the Arctic system, since it provides an insulating cover for vegetation and soil. Snow insulation is recognized as the primary control over soil thermodynamics (Lawrence and Slater, 2010), and soil temperature is closely connected to physical (i.e. permafrost active layer depth) and biogeochemical (i.e. decomposition, greenhouse gas emission) processes (Peng et al., 2016). Observations show that snow cover changes have played a major role in a warming trend of permafrost soils of approximately 0.3 ^∘C per decade (Biskaborn et al., 2019; AMAP, 2017). This warming may lead to increased microbial activity, decomposition rates and bioavailability of previously frozen soil carbon. Since permafrost soils contain approximately 1600 Pg carbon, accounting for half of the global soil carbon storage (Hugelius et al., 2014), there is ample potential for these changes to lead to the release of the greenhouse gases CO₂ and methane. This has the potential to accelerate global warming (Schuur et al., 2015), which underlines the need for a better understanding of drivers and potential feedbacks to better predict the rate and magnitude of future carbon exchange.

Despite numerous field-based and modelling efforts to date, it is still uncertain whether the Arctic will act as a carbon source or sink in the future (McGuire et al., 2012; Virkkala et al., 2021). The predicted future carbon balance varies widely among models – between a loss of 641 Pg C and a gain of 167 Pg C under RCP8.5 (McGuire et al., 2018) – depending on the representation and level of detail of key processes such as soil temperature and vegetation dynamics (Schuur et al., 2015; McGuire et al., 2018). One of the key goals of model development is to decrease uncertainty of simulations by refining these processes. While extensive research has been carried out on the mechanics of the growing season, few studies have been directed at cold season processes. Recent studies suggest that the contribution of the non-growing season to the annual carbon budget may have been underestimated (Pirk et al., 2015; Mastepanov et al., 2013). A recent meta-analysis by Natali et al. (2019) found significant wintertime carbon loss and highlighted the large spread in model simulations of non-growing season greenhouse gas emissions. Models generally underestimated the observed winter flux emissions due to the inaccuracies in their simulation of cold season respiration. Natali et al. (2019) stress the need to revise the impact of environmental drivers and feedbacks in models. Collectively, these efforts demonstrate that the influence of cold season processes on the annual carbon balance is larger than previously suggested.

The ability of models to simulate physical and biological processes in the soil is limited by the complexity of their representation of snow. A recent snow-related model evaluation project analysed the performance of models with different complexity – focusing on variables such as snow-covered area and snow season length. This SNOWMIP found that a dynamic simulation of internal snowpack processes, such as density and temperature calculations, is critical to simulate snow thermal profiles (Krinner et al., 2018). In addition, more complex snow schemes perform better when simulating cold season processes (Vionnet et al., 2012; Slater et al., 2017; Wang et al., 2016). To balance computational efficiency and the need for detail, most ecosystem models use an intermediate-complexity multi-layer snow module (Vionnet et al., 2012; Krinner et al., 2018). Such schemes may not capture fine-scale internal snowpack processes such as the evolution of high-density wind slab layers, but they are complex enough to simulate key physical processes – compaction, freeze–thaw cycles, and liquid water retention – that influence the thermal dampening property of snow. Since LPJ-GUESS had a single-layer static snow representation, it was found to deviate from observational records of air–soil temperature relationships – simulating cooler winter conditions and performing poorly when compared to eight land surface models (Wang et al., 2016). This showed, combined with previous research, that the snow representation in LPJ-GUESS needed to be revised to better capture Arctic cold season conditions.

The primary aim of this study is to improve LPJ-GUESS's simulation of the insulating effect of snow. By developing and integrating a dynamic intermediate-complexity snow scheme, we also aim to improve the soil temperature and biogeochemistry simulation. To investigate the effect of the new snow scheme in LPJ-GUESS, we set out to quantify the impact on physical variables, i.e. the direct impact of snow insulation on soil temperature and permafrost conditions. To further evaluate the snow-related influence on biogeochemistry – such as changes in growing season length – we analyse a set of biogeochemical variables. Due to differences in soil temperature, we expect to see changes in ecosystem productivity, heterotrophic respiration, and soil carbon pools. Moreover, we analyse the changes to vegetation dynamics and composition. The updates to the model will allow us to assess other snow-related processes and feedbacks on a global scale, such as the impact on surface albedo and food access to herbivores, in the future.

2 Materials and methods

LPJ-GUESS is a process-based dynamic vegetation model, widely applied on regional and global scales (Smith et al., 2001, 2014). For this study, we used a customized arctic version of LPJ-GUESS 4.0 (subversion 9905). The model simulates soil freeze–thaw processes and is applicable to studies of processes at northern high latitudes (Miller and Smith, 2012). In this study we restrict simulations to the northern circumpolar region (above 60^∘ latitude) with a spatial resolution of 0.5^∘ × 0.5^∘. The CRUNCEP global reanalysis climate dataset version 7 was used as input for all of our model simulations (Viovy, 2016). We ran the model with a 500-year spin-up period to establish an equilibrium vegetation state and a 40 000-year offline spin-up period for soil conditions.

LPJ-GUESS simulates vegetation dynamics on individual and patch scales, taking into account growth, competition for resources, and disturbances. This feature makes it possible to assess how changes in environmental conditions affect vegetation distribution and composition. In this study, we applied 15 plant functional types (PFTs) that characterize arctic ecosystems (see Table S3 in the Supplement). Permafrost dynamics follow Wania et al. (2009) and are simulated using the physical characteristics of 15 soil layers, each 10 cm thick. Soil thermodynamics is governed by climate and snow conditions, and the thermal properties of each soil layer depend on the ice, water, air, mineral, and organic soil fractions. The layer-specific thermal properties define the rate of heat transfer through the soil column. For more details on the model structure, see Smith et al. (2001, 2014), Wania et al. (2009), and references therein. To assess the newly developed intermediate-complexity snow scheme's performance and influence, we conducted simulations with both the old Static and the new Dynamic snow schemes.

2.1 Static snow scheme

The Static snow scheme, which has been in use in LPJ-GUESS until now, treats snow as a single layer with constant values for thermodynamic parameters. Snowfall is simulated on any given day when precipitation (P, mm) occurs and air temperature (T_k, ^∘C) is at or below T_max (^∘C) – which is the temperature maximum at which precipitation occurs in snow form. T_k denotes the temperature of a layer k, in this case the air layer. Snow density (ρ_k, kg m⁻²) and snow thermal conductivity (K_k, $\mathrm{J}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$) are constant at 362 and 0.196, respectively. Snow heat capacity (C_k, $\mathrm{J}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{K}}^{-\mathrm{1}}$) is calculated by Eq. (1) (Fukusako, 1990).

$$\begin{array}{}\text{(1)}& C_k=\mathrm{1000}{\mathit{\rho }}_k(\mathrm{0.185}+\mathrm{0.00689}{T}_k)\end{array}$$

Compaction processes are not represented in the Static scheme. Snowmelt (melt, mm) is governed by air temperature and precipitation and follows a linear function as shown in Eq. (2) (Choudhury and DiGirolamo, 1998).

$$\begin{array}{}\text{(2)}& \text{melt}=(\mathrm{1.5}+\mathrm{0.007}\cdot P)(T_k-T_{max})\end{array}$$

The snowpack is homogeneous in its physical properties, and neither internal processes nor seasonal dynamics are simulated using the Static scheme. Using this approach, snow conditions are assumed to be uniform across the Arctic regardless of air temperature regime or seasonal snow dynamics. Due to the heterogeneity in Arctic surface and local climatic conditions, this scheme has a limited ability to represent the variability in high-latitude ecosystems – this is the main shortcoming of the Static snow scheme that this study sets out to improve.

Figure 1 Snowpack structure and physical properties. Where z shows layer depth, W water content, I ice content, ρ density, C heat capacity, K thermal conductivity, D thermal diffusivity, and T layer temperature. A detailed list of the used variables can be found in Table A1.

2.2 Dynamic snow scheme

The schematic structure of the multi-layer snow scheme is shown in Fig. 1. The occurrence of snowfall on any given day depends on air temperature and precipitation, using the same principle as for the Static scheme. Fresh snow density (ρ_fresh, kg m⁻³) is calculated by taking into account air temperature and wind speed, following Eq. (3), where a, b, and c are scaling parameters defined by Vionnet et al. (2012) (for parameter values see Table A1 in the Appendix).

$$\begin{array}{}\text{(3)}& {\mathit{\rho }}_{\text{fresh}}=a+b\cdot T_{max}+c{U}_{\mathrm{10}}^{\mathrm{0.5}}\end{array}$$

U₁₀ denotes the 10 m height wind speed (m s⁻¹), following the detailed snowpack model Crocus (Vionnet et al., 2012). To avoid unrealistically low snow density values that may occur in rare cases, the density minimum is set to 100 kg m⁻³.

The new snow scheme simulates internal snowpack dynamics with up to five snow layers, taking into consideration each layer's depth. Fresh snow either initiates a snowpack or is added to already existing snow layers. If the freshly fallen snow is added to the snowpack, the physical properties of the snow layer are updated. The number and thickness of snow layers are defined according to predefined thresholds: a new snow layer is initialized when an existing layer exceeds twice the prescribed threshold height (2 × 100 mm). If a single snow layer exists but does not reach the minimum height (set to 50 mm), the shallow snow layer properties are combined with the top soil layer. Thereafter, their properties (ice, air, and liquid water fractions and heat capacity) are scaled using weighted averages based on the layer's ice, water, and air fractions for the sake of computational stability. In the case where all five layers exceed the prescribed maximum threshold, the bottom layer accumulates snow in order to preserve and align vertical resolution near the surface of the snowpack. The snow layer density (ρ_k) and depth (z_k, m) relationship is described by Eq. (4), where I_k (kg m⁻²) defines the ice content of a layer (Lawrence et al., 2019).

$$\begin{array}{}\text{(4)}& z_k={\displaystyle \frac{I_k}{{\mathit{\rho }}_k}}\end{array}$$

The density of a snow layer changes through compaction, which is simulated by two processes: (1) mechanical compaction due to pressure from the overlying snow layers as shown in Eq. (5) (Best et al., 2011).

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial {\mathit{\rho }}_k}{\partial t}}={\displaystyle \frac{{\mathit{\rho }}_{k}g{M}_k}{{\mathit{\eta }}_k}}\exp \left({\displaystyle \frac{k_s}{T_{max}}}-{\displaystyle \frac{k_s}{T_k}}-{\displaystyle \frac{{\mathit{\rho }}_k}{{\mathit{\rho }}_{\mathrm{0}}}}\right)\end{array}$$

The increase in the snow layer's density (∂ρ_k) depends on the mass of overlying layers (M_k, kg). η_k (10⁶ Pa s) denotes the compactive viscosity factor, k_s is an empirical constant defined by Best et al. (2011) with a value of 4000 K, and ρ₀ is a reference density (50 kg m⁻²). Snow density may also change by (2) phase changes as a result of freeze–thaw processes within the layers. If a layer's snow or liquid water content changed during freeze–thaw events, its depth and density properties are recalculated, taking into account the snow and ice fractions of the layer as shown in Eq. (4).

In contrast to the Static formulation, phase changes within the snow layers depend on the layer's internal temperature, and this controls the melting process in the Dynamic snow scheme. This development enables the simulation of mid-winter melt events and ensures an improved representation of internal snowpack thermodynamics. Upon melt, each layer can retain a fraction of liquid water based on Eq. (6), where rw_min and rw_max are empirical constants and ρ_t is a reference density (Wang et al., 2013; Anderson, 1976).

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{W}_{\text{cap},max}=& I_k[{\mathrm{rw}}_{min}+({\mathrm{rw}}_{max}-{\mathrm{rw}}_{min})\\ & \cdot max\left(\mathrm{0},{\displaystyle \frac{({\mathit{\rho }}_{\mathrm{t}}-{\mathit{\rho }}_k)}{{\mathit{\rho }}_{\mathrm{t}}}}\right)]\end{array}\end{array}$$

If the liquid water content (W_k, mm) of a layer exceeds the maximum water holding capacity ($W_{\text{cap},max}$, mm), water is passed to the layer underneath following a simple bucket model. Rain-on-snow events (ROSs) are simulated if it rains while a snowpack is present. The energy of rainwater may induce phase changes in the snow layers. The overflow liquid water is forwarded to the underlying snow layers and lastly to the top soil layer to percolate to the soil or to be discharged as runoff.

Each layer is characterized thermodynamically by the following physical properties: density, temperature, thermal conductivity, heat capacity, and diffusivity (D_k, m² d⁻¹). Thermal conductivity is calculated using density as shown in Eq. (7) (Best et al., 2011), following a power function ($\mathrm{J}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{K}}^{-\mathrm{1}}$).

$$\begin{array}{}\text{(7)}& K_k=\mathrm{2.22}{\left({\displaystyle \frac{{\mathit{\rho }}_{{}_k}}{{\mathit{\rho }}_{\mathrm{0}}}}\right)}^{\mathrm{1.88}}\end{array}$$

Heat capacity is determined by taking into account snow layer density and temperature according to Eq. (1). The snow diffusivity is calculated by Eq. (8). Soil and snow layer temperatures are computed, taking into account each layer's thermal conductivity, heat capacity, and height, using the Crank–Nicolson finite difference method to solve Eq. (9) (Lawrence et al., 2019).

$$\begin{array}{}\text{(8)}& {\displaystyle }{D}_k={\displaystyle \frac{K{}_k}{C_k}}\text{(9)}& {\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(D\left(k\right){\displaystyle \frac{\partial T}{\partial z}}\right)\end{array}$$

Figure 2 Steps in the daily computational cycle for the Dynamic snow scheme. Blue arrows indicate workflow in case a snowpack is present on the ground, while orange arrows show steps when a snowpack is absent.

The computational cycle ends by rearranging the layers based on the depth thresholds, taking into account the potential liquid water content. First we re-calculate each snow layer's depth based on the amount of snow and liquid water using Eq. (4). We then re-arrange the layers by using the leaky-bucket approach, where the snow layers are filled up from the bottom layer (closest to the surface). If the threshold depth is reached, a new snow layer is initiated and the process continues until the total depth of the snowpack is distributed to the specific snow layers. The overflow meltwater is passed to the soil for percolation after this step. This cycle is repeated each day when there is a snow or rain-on-snow event. The daily snow cycle of the Dynamic snow scheme is depicted in Fig. 2.

Besides the changes in the representation of snow, the calculation of heterotrophic respiration below 0 ^∘C was changed following a recent data synthesis (Natali et al., 2019) to better represent arctic conditions. This adjustment was implemented for both the Static and Dynamic schemes. The minimum decomposition temperature was set to −20 ^∘C, and the Q10 value was set to 2.9. A comparison between the old and new functions is shown in the Supplement (Fig. S1 in the Supplement). This adjustment led to higher soil respiration in both schemes during the cold season compared to the old model set-up.

The implemented processes and physical representations are simpler than in dedicated, high-resolution snow models – such as Crocus and SNOWPACK (Lehning et al., 2002; Vionnet et al., 2012) – but reflect the model improvements identified as being most important in previous model inter-comparison studies (Krinner et al., 2018). These improvements enable us to simulate a more realistic range of snow conditions and soil thermal conditions across the Arctic.

Figure 3 Snowpack dynamics at the Zackenberg GeoBasis station. Density values for the layers are extrapolated – from three and five layers for the observational and modelled data, respectively. The colours of the snowpack indicate snow density.

2.3 Simulation set-up

The performance of LPJ-GUESS using both snow schemes was compared at the site and regional levels. Modelled properties were compared to observational datasets, when available. We quantified the correspondence between simulated and measured variables using statistical methods.

• i.

  Site-level comparison.

  To highlight the differences in snowpack dynamics using the two snow schemes, we compiled a detailed single-site model–data comparison of the internal snowpack structure (Fig. 3). Afterwards, we ran the model for five well-studied northern high-latitude sites in order to identify how well the two snow schemes can simulate snow and soil temperature at the site level. These sites are Abisko, Zackenberg, Bayelva, Kytalyk, and Samoylov – see site details in Table S1 in the Supplement. Measurements of snow depth and soil temperature were sorted and averaged on a daily basis – 10 years for the simulations and all available years for observations at each site. Model outputs were examined and compared to these time series to evaluate the snow schemes' ability to simulate snow depth and soil temperature seasonality adequately.

• ii.

  Regional simulations.

  We conducted simulations for a set of Russian sites (256 sites) which were part of the study by Wang et al. (2016), as a follow-up, and re-evaluate the snow insulation effect in LPJ-GUESS over a large region. First, snow depth and soil temperature data were sorted monthly for each site for the years 1980–2000. Site observations were provided by the All-Russian Research Institute of Hydrometeorological Information – World Data Centre (RIHMI-WDC; http://meteo.ru/, last access: 3 January 2019). Following this, averages were calculated for December, January, and February. The difference between soil (25 cm depth) and air temperature – henceforth ΔT – was used as a proxy to evaluate the strength of the model-simulated insulation effect. Snow depth, soil temperature, and ΔT series were grouped according to air temperature to evaluate the insulation capacity under different temperature regimes.

• iii.

  Pan-Arctic simulations.

  Finally, we conducted model simulations across the Arctic to assess the effect of changing the snow scheme on selected physical and biogeochemical variables and vegetation properties. When applicable, variables were averaged over December, January, and February to emphasize the effect on the winter season. Instead of the absolute results, we show the difference between the set of simulations, calculated as the difference between Dynamic and Static model outputs.

3 Results

3.1 Site-level simulations

Prior to the evaluation of the large-scale performance of the new Dynamic snow scheme, we conducted a single-site comparison to examine the validity of the results. These detailed snowpack observations from Zackenberg helped to determine whether the Dynamic scheme can simulate internal snowpack dynamics, snow depth, and snow density.

Figure 4 Seasonal cycle of (a) snow depth and (b) soil temperature at 25 cm depth for the studied sites, comparing model simulations and observations. Site statistics show the spread of monthly snow depth and soil temperature values for the respective sites – excluding the summer months (June–July–August).

We established the ability of the new snow scheme to simulate snow conditions by comparing a simulated snowpack with snow depth and density observations from Zackenberg (2013–2014 snow season). Figure 3 presents the observed and simulated snowpack by the Dynamic and Static schemes. This figure shows that the Dynamic scheme simulates comparable snow depth and that the simulated snow densities follow the observed snow density pattern through the snow season. Density values are compared qualitatively, since it is difficult to accurately align the observational and modelled layer densities. To be consequent, we used global climatic forcing data for all simulations in this study, including this site-scale comparison. This fact should be taken into account when interpreting the model–data comparison in this section – as some of the differences may be derived from the differences in climatic data.

There are lower densities early in the snow season, with fresh snow having low density, while density increases in late spring, during the melt season. The Static scheme with constant snow density simulated a somewhat higher-than-observed snow depth. Thermal properties in snow layers are derived from density, and this is especially important in the Dynamic snow scheme. Dynamically simulated density translates to more realistic thermal conductivity dynamics, which governs the rate of heat transfer through the snowpack. This feature is essential in simulating a reliable atmosphere–snow–soil heat transfer interaction. The difference in snow depth between the Static and Dynamic simulations is most visible at the end of the snow season, before the start of snowmelt – as indicated in Fig. 3, bottom panel.

Overall, the new scheme reproduces the snow dynamics over the cold season better than the Static scheme. Taken together, these results suggest that the Dynamic scheme is skilled in simulating the snowpack's internal structure and dynamics. Since the Static scheme has a constant snow density throughout the snow season, the Dynamic scheme is expected to better capture the seasonal behaviour of snow and soil conditions. The Zackenberg site comparison indicated that the Dynamic scheme successfully integrated these key processes affecting the density over the snow season. In this study, we used a global climate forcing dataset, which may explain some of the observed model–observation differences. The mismatch between snow observations and simulations is influenced by the use of the global model forcing dataset instead of site-specific temperature, precipitation, or snowfall time series.

Table 1 RMSE for soil temperature and snow depth for the applied snow schemes for the single site simulations.

We moved to a multi-site analysis to compare the Static and Dynamic snow schemes on five well-documented sites. To assess the performance of the two snow schemes, we composed seasonal cycles based on monthly averages of (a) snow depth and (b) soil temperature at 25 cm depth, shown in Fig. 4. The corresponding root-mean-squared error (RMSE) for each study site is shown in Table 1. Generally, the Dynamic scheme shows only minor improvements in the simulation of snow depth. Despite this, modelled soil temperatures are much closer to the observed values for all sites, especially during the winter months. This behaviour highlights that changing the internal snowpack dynamics with the Dynamic snow scheme had a significant effect on soil temperature, even when the simulated snow depth differed marginally. The changes in soil temperature are due to the differences in snow thermal properties, which significantly influenced the insulation capacity of the snowpack.

Figure 5 Comparison of the observed and modelled snow insulation effect at the Russian sites between observations and model simulations using the Dynamic and Static schemes. (a–c) Soil temperature and snow depth relationship. (d–f) Difference in air–soil temperature and snow depth relationship. Snow depth presented on the horizontal axis is classified in 5 cm depth bins. Colours indicate different air temperature regimes, and upper and lower bars show the 25th and 75th percentiles.

The implementation of the new snow scheme resulted in significant changes within the snowpack that are responsible for the improved soil temperature simulation. The Static scheme applies constant snow density and thermal conductivity, which defines the rate of heat transfer through the snowpack. In the Dynamic scheme, thermal conductivity is dependent on the dynamically updated density; therefore the new scheme can achieve a more realistic simulation of snow heat transfer dynamics throughout the snow season – depending on environmental conditions. The Dynamic scheme simulates snow thermal conductivity in a range from 0.04 to 0.5 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, which aligns well with literature estimates of 0.021–0.65 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$. This feature enables the simulation of a wide range of conditions across the Arctic, as opposed to the general conditions assumed by the Static scheme.

The statistical comparison (site statistics) shows that there is a smaller variance of modelled values of soil temperature using the Dynamic snow scheme, which indicates an improvement in comparison to the Static simulations' outputs. The RMSE (Table 1) also shows that the Dynamic scheme provides an improved fit of simulated soil temperature and snow depth at most sites. Overall, we conclude that, with the Dynamic scheme, the model is able to simulate snow and soil temperatures that correspond better with the observed ranges.

Table 2 Pan-arctic mean values for the studied variables for the Static and Dynamic simulations and their respective differences.

3.2 Russian site simulations

Following the Dynamic scheme's improved performance at the site level, we further evaluate the model's performance at the regional scale for the same sites as in the previous model intercomparison by Wang et al. (2016) that highlighted shortcomings in the snow scheme of LPJ-GUESS. Figure 5 shows the snow insulation effect over a set of Russian sites, using the two snow schemes, where the coloured bars show different temperature regimes. The figure is compiled from 20 winter season average values of near-surface soil temperature (25 cm depth) and snow depth per site. Due to the air-temperature-based classification, the number of samples per bin is not balanced, which led to an uneven number of values allocated to the different groups. The top row of Fig. 5 shows that the Dynamic snow scheme has better skill in simulating the relationship between soil temperature and snow depth than the Static scheme. It must be noted that there is a clear difference between the current Static scheme simulations and results reported by Wang et al. (2016), which is due to recent updates in the model, independent of the snow module, and the different climate forcing dataset used in this study.

It is apparent from the ΔT and snow depth relationship (Fig. 5, bottom row) that the Dynamic scheme reproduces the observed insulation effect well. Unlike the Static scheme, the new snow module can also simulate the different insulation behaviour depending on the air temperature regimes. The improved performance of the Dynamic scheme is confirmed by the root-mean-squared error (RMSE), shown in the Supplement (Table S2 in the Supplement). RMSE decreased significantly for both the soil temperature–snow depth and ΔT–snow depth relationships. This regional analysis confirmed that the new Dynamic snow scheme has an improved skill in simulating winter soil conditions.

3.3 Pan-Arctic simulations

To assess how the two snow schemes differ in simulating seasonal snow across the Arctic, we subtracted output variables from simulations with the Static module from those with the Dynamic module. We calculated average conditions for winter (December–January–February) and summer (June–July–August) for the period 1990–2015. The mean pan-arctic seasonal dynamics of snow depth, soil temperature, and upper soil water content are shown in the Table 2.

Figure 6 (a) Snow depth difference in winter months and (b) maximum ALD difference, calculated by subtracting the Static from Dynamic simulation outputs. Modelled permafrost extent is based on mean annual ground temperature (MAGT) and plotted against the permafrost cover estimate by Obu et al. (2019a) (TTOP model). Simulated absolute snow depth is shown in Fig. S2 in the Supplement and ALD in Fig. S3 in the Supplement.

3.3.1 Impacts on physical variables

Figure 6a shows the difference in simulated wintertime snow depth. The Dynamic scheme shows an overall lower snow depth across the Arctic with the most pronounced changes in coastal Norway and in western Siberia. On average, the snow depth for the Dynamic scheme is 6 cm lower due to the implementation of snow-related processes affecting snow density and consequently snow depth.

Figure 7 Near-surface soil temperature (25 cm depth) difference between the Dynamic and Static simulations, for winter (a) and summer (b) seasons. Differences are calculated by subtracting the Static from Dynamic simulation outputs. The absolute simulated soil temperatures using the two snow schemes are shown in the Supplement in Fig. S4.

The main aim of developing the new snow scheme was not only to enhance the simulation of snow depth but also to improve the simulation of snowpack properties that directly affect soil conditions. Therefore, we investigated how the internal structural changes in the representation of snow influenced soil temperatures. The soil temperature differences shown in Fig. 7 reveal that the new snow scheme influenced the winter season to a large degree, both within and especially outside of the permafrost region. Winter soil temperatures are higher with the Dynamic scheme, while it results in a cooler near-surface soil temperature during the summer. A closer look at the monthly soil temperature values in Fig. 11 showed that spring months are cooler for the Dynamic scheme but that the difference between the two schemes decreases towards the end of summer. This pattern is more pronounced in the permafrost-underlain regions. This shows that the Static snow scheme has too little insulation and results in soil temperatures that are too cold during the winter months, as we also show in the site simulations. Moreover, the Static scheme also does not insulate soils sufficiently during the springtime when air temperatures rise above 0 ^∘C, which allows the soil to warm up more quickly even in the presence of a snowpack.

The depth to which the top soil thaws during summer, and refreezes in winter, in permafrost areas is called the active layer depth (ALD). The difference in the seasonal maximum active layer depth for the model simulations is shown in Fig. 6b. Since the Dynamic scheme had warmer soil temperatures, the modelled permafrost extent is smaller than with the Static scheme. We compared our model simulations with a recent satellite-driven permafrost extent estimate by Obu et al. (2019a) – from here on referred to as the TTOP model. Modelled permafrost extent was defined by the area where the mean annual ground temperature (25 cm depth) was below zero. The Dynamic scheme's permafrost extent is much closer to the TTOP model's estimate, while the Static scheme simulates a much larger permafrost extent, as shown in Fig. 6c. The Dynamic scheme's computed areal permafrost cover, while improved compared to the Static scheme, still overestimates the TTOP model estimates by approximately 5 % (see Table 2).

Figure 8 Mean fractional water content of the upper soil column in May, June, and July, using the Static (a) and Dynamic (b) schemes and their difference (c). Differences are calculated by subtracting the Static from Dynamic simulation outputs.

Besides governing the physical state of permafrost, snow and soil temperature also have a large influence on the temporal and spatial patterns of soil water content. We show the simulated upper soil column water content in Fig. 8. This figure shows the mean fractional soil water content over May–June–July. Soil water content was calculated using the average seasonal liquid water content from the top 50 cm of the soil column. Soil water content is represented as the fraction of the available capacity between the wilting point and field capacity, and therefore frozen water is not included in these values. Figure 8c shows that there is a higher water availability within the permafrost region using the Static scheme. Water availability is a key driver of the start of the growing season, nutrient availability, and vegetation dynamics. The time-series analysis of upper soil water content highlights that the snowmelt rate is not significantly different between the schemes. Still, there is a large difference in soil temperature dynamics. The Static scheme's soil temperature increases more rapidly during the spring than the Dynamic scheme's soil temperature (see Fig. 11a and b). This results in an earlier onset of snowmelt and earlier increase in soil water availability and nitrogen mineralization. This affects productivity, which we assess in the coming sections. Although the difference in water content and nitrogen mineralization between the snow schemes converges towards zero as summer progresses, we show that the change in snow scheme had a lasting effect beyond the cold season.

Overall, the new snow scheme had a substantial effect on winter soil temperatures. As a result, summer conditions were also altered by the snow scheme updates. It is apparent that the largest changes in snow depth and temperature coincide. For instance, along the Norwegian coast and in central Siberia. Taken together, our results show that the Dynamic snow scheme improved the simulation of physical variables.

Figure 9 The difference in simulated soil (a) and vegetation carbon pools (b) between the two schemes. Differences are calculated by subtracting the Static from Dynamic simulation outputs.

Figure 10 The difference in normalized heterotrophic respiration for winter (a) and summer (b) between the two schemes. Differences between simulated Dynamic and Static winter (c) and summer (d) NPP. Differences between simulated Dynamic and Static winter (e) and summer (f) NEE. Differences are calculated by subtracting the Static from Dynamic simulation outputs. Simulations for the two schemes are shown in Figs. S5 (heterotrophic respiration), S6 (NPP), and S7 (NEE) in the Supplement.

3.3.2 Impacts on biogeochemical variables

Besides the impact on soil thermodynamics, we investigated how key biogeochemical components – such as productivity and carbon pools – were affected. The changes across seasons and permafrost conditions are summarized in Figs. 11 and 12.

Our simulated soil carbon pools (Fig. S6 in the Supplement) deviate from literature values (Hugelius et al., 2014) and are consistently lower across the Arctic. The main reason for this is the model's representation of soil organic matter processes. Soil carbon and nitrogen are represented by pools that exist in the top 50 cm of the soil column (Smith et al., 2014) and are thus only influenced by near-surface conditions. Moreover, peatlands are not explicitly represented. The differences in soil carbon between the schemes, as shown in Fig. 9a, coincide spatially with the highest differences in soil temperature. This suggests that the changes in soil temperature influence soil carbon in the model and therefore the rate of respiration from soils as well. Vegetation carbon pools (Fig. 9b) are higher in the non-permafrost region using the Dynamic snow scheme (see Table 2 for mean values). Since the evaluation of soil carbon is not the focus of this study, soil carbon outputs were used to normalize the heterotrophic respiration to be able to interpret the relative differences between schemes (Fig. 10a and b). To do so, we divided the heterotrophic respiration by the soil carbon estimates for the respective simulations using the two snow schemes. With the Dynamic scheme, summer soil respiration decreased across the Arctic. Winter respiration, on the other hand, increased, except for in eastern Siberia. These changes in soil respiration can be attributed to changes in soil temperature, as shown in Fig. 7.

Figure 11 Seasonal dynamics of snow depth, soil temperature (25 cm depth), and fractional soil water content within (a) and outside of the permafrost region (b). Seasonal dynamics of NEE and NPP within (c) and outside of the permafrost region (d).

The difference in net primary productivity (NPP) between simulations with the two snow schemes for both winter and summer is shown in Fig. 10c and d, where positive NPP means more carbon uptake by the vegetation. We note an impact of the different snow schemes on summer productivity, caused by the different soil thermodynamics and soil water availability during the spring and early summer period. This artefact is also visible in the simulated pan-Arctic NEE, in Fig. 10e and f, where negative NEE values indicate a stronger uptake of carbon by ecosystems. The positive difference in winter NEE (e) shows that there is a higher carbon release in the winter season for the Dynamic scheme in central Europe, western Siberia, and coastal Norway. The mean winter NEE of the Dynamic scheme more than doubled (Table 2). Compared to the Static scheme, which relates to both the change in soil respiration and NPP. Figure 11c and d show some interesting contrasts regarding the seasonal carbon fluxes. Permafrost-underlain regions (Fig. 10c) experience little difference in the simulated NEE. The Dynamic scheme simulates lower peak summer NPP. Winter NEE in the non-permafrost region is higher using the Dynamic scheme, indicating a larger carbon uptake by the vegetation. On the other hand, we can observe an increased sink capacity (more negative NEE) during the summer months.

3.3.3 Vegetation composition and distribution

Vegetation composition and distribution depend on the changes in physical and biochemical variables described in the previous section. Therefore, we investigated how the application of the two snow schemes affected vegetation distribution to determine if there are shifts in dominant plant functional types (PFTs) as a result of using different snow schemes. The dominant PFT for each simulated grid cell was determined by selecting the PFT with the highest maximum leaf area index (LAI) during the simulation years (1990–2015). Using the Dynamic snow scheme, roughly half of the sites are dominated by summergreen low shrubs and boreal needle-leaved evergreen trees (LSS with 25 % and BNE with 23 %; see Table S3). Prostrate dwarf shrubs, (SPDS), graminoid and forb tundra (GRT), and boreal needle-leaved summergreen trees (BNS) accounted for 20 %, 8 %, and 7 % dominance, respectively. For an easier comparison between the Static and Dynamic simulations, PFT classes were grouped into forest, open grass, shrubs, and no vegetation categories after determining the dominant PFT in each grid cell (classification based on Wolf et al., 2008). This classification showed that grassland classes dominate (56 %), followed by forest cover (36 %). Shrubs dominate at 29 % of the simulation sites. There is a negligible number of sites with mostly bare soil. When comparing the spatial pattern of dominant vegetation groups, we noted that there is only a marginal difference between the Static and Dynamic simulations (see Fig. S10 in the Supplement). Changes in group dominance between the Static and Dynamic simulations occurred at approximately 10 % of the sites; see Fig. S11 in the Supplement. The Sankey diagram shows the direction of change between the three groups. Myers-Smith et al. (2011) suggest that increased soil temperature leads to a shift to an increased forest height and shrub cover. Grass-to-shrub domination change (3.1 %) was the most prevailing change, which indeed points towards an increase in vegetation height. However, we observe a shrub-to-grass (1.95 %) shift in domination at the same time; therefore, we cannot conclude the main direction of changes.

Figure 12 Relationship between variables during the winter and summer seasons, within and outside of the permafrost region. The colour of the boxes indicates the direction and qualitative magnitude of changes in the variables, based on the relative difference (%) between Dynamic and Static schemes – shown in the boxes below the variables. Variables are shown as follows: GPP: gross primary production (g m⁻²); NPP: net primary productivity (g m⁻²); NEE: net ecosystem exchange (g m⁻²); R_a: autotrophic respiration (g m⁻²); R_h: heterotrophic respiration (g m⁻²); N_min: net nitrogen mineralization (kg N ha⁻¹); soil water content: fraction. The relative changes in soil water content (in %) are high in the winter due to the low fractional water content values as we only account for the liquid soil water and do not consider the amount of frozen water in the soil. These changes correspond to small absolute changes in fractional soil water content. ^∗ Note: temperature differences are shown as absolute differences for easier readability.

3.4 Cause-and-effect relationships

As many of the reviewed processes interact with each other in a complex, non-linear manner, change in one variable may not translate to a direct impact on another variable. To provide an overview of our findings regarding the physical and biogeochemical processes, we created a flow chart showing observed changes in modelled state variables and their connections. Figure 12 shows the difference between simulations using the two snow schemes – calculated by subtracting the Static from Dynamic results. Reddish box colours show that the Dynamic scheme had higher values, and blueish colours show that the Dynamic scheme simulated lower values than the Static scheme. The lightness and darkness of the colours indicate the magnitude of changes between the winter and summer seasons qualitatively.

Each box contains the computed difference, and Table 2 summarizes the mean changes in these key variables. Considering the spatial pattern across the Arctic, we conclude that the pattern of changes and differences between the Static and Dynamic simulations vary depending on the presence or absence of permafrost cover. For a more detailed evaluation, process relationships are therefore divided into permafrost and non-permafrost regions. Since snow depth only affects these variables indirectly, through insulation, it was not included in the feedback graph. The choice of snow scheme induced changes in near-surface temperature (T), which is the key governing factor over these variables. In general, higher soil temperatures during the winter season prompt a positive response in respiration, soil water content, and vegetation primary productivity. Soil temperature increased to a greater degree in the non-permafrost region during the winter season. The same increasing pattern is observed for heterotrophic (R_h) and autotrophic respiration (R_a), soil water content, and NEE. Nitrogen mineralization decreased in the wintertime, with a larger decrease outside of the permafrost region. In contrast, summer months' average soil temperature showed an overall decrease in the permafrost region. This is due to the different thermal soil and snow dynamics in the two applied snow schemes. We observed rapid heat loss using the Static scheme, resulting in insufficient insulation during the snow season. The same feature causes soil temperature to rise rapidly in the spring, when air temperature is already above zero but a snow cover remains present. Faster soil warming leads to increased soil water availability that affects productivity. Respiration and NEE are slightly reduced for both permafrost and non-permafrost regions in the Dynamic scheme during the summer. Differences noted for the summer are smaller for all variables than in the winter season.

4 Discussion

4.1 Snow scheme dynamics

The site-level analysis shows that the new Dynamic scheme is able to simulate snow height and density adequately due to the implementation of physical processes and a dynamic representation of snow properties. The integrated mechanistic compaction scheme and phase changes within snow layers make it possible to simulate heterogeneous snow density and thermal properties within the snowpack. This influences the simulated snow density directly by altering snowpack structure. Density regulates heat transport rate through snow layers by affecting thermal conductivity (Eq. 7): lower density results in a more insulating cover, whereas higher density and compacted layers are a better heat-transferring medium and exhibit lower insulation. In the Static scheme, snow density was assumed constant through the snow season and across all study sites. Such static snow representation is unsatisfactory when simulating arctic conditions (Krinner et al., 2018). The new snow scheme provides an improved framework for a mechanistic snow season simulation.

The single site simulations (Sect. 3.1) provide reasonably consistent evidence that the new snow scheme's implementation leads to significant changes in near-surface temperature simulation – especially at Abisko, Bayelva, and Samoylov. As shown by Chadburn et al. (2017), the site-wise model–data comparison is challenging since point measurements may not be representative of a larger area due to the complexity in topography and vegetation conditions. The model–observation fit may be improved by using site-specific climatic forcing instead of a global gridded dataset.

To avoid site-specific problems in the interpretation of simulations, we also evaluated the model at a regional scale. By comparing the results of the Russian site simulations (Sect. 3.2) with those of Wang et al. (2016), we conclude that the development of the representation of snow in LPJ-GUESS significantly improved air–soil temperature and snow depth–soil temperature relationships. The Dynamic snow scheme's insulation capacity followed a quasi asymptotic trend, increasing with snow depth and slightly levelling out after reaching the so-called effective depth at 30–40 cm (Slater et al., 2017). The insulation capacity was, in general, slightly lower than observations, with a notable underestimation when snow depth is below 20 cm. Nonetheless, these results are a vast improvement over the old Static scheme, as shown in Fig. 5.

The RMSE values (Table S2) also show that the Dynamic scheme better captured the observed the soil temperature and snow depth relationship than the Static scheme. RMSE was slightly higher for the coldest air temperature regime for both snow schemes. We note that the Static schemes' performance differs from what was shown in Wang et al. (2016). The reasons for these differences are developments of the model since then in components other than the snow scheme and also the different meteorological forcing used in this study. Our results indicate that the enhancement of snow-related processes improves the simulation of soil temperatures in LPJ-GUESS and that the model can be more reliably applied to assess the impact of environmental changes on the arctic carbon cycle.

4.2 Impact on physical and biogeochemical variables

The changes to the snow insulation capacity in the Dynamic scheme had a significant effect on permafrost conditions. Our pan-Arctic results showed that the Static scheme simulated near-surface soil temperatures that were too cold in winter and too warm in summer. Permafrost extent simulated with the Dynamic scheme agreed more closely with the permafrost estimate by Obu et al. (2019a), as shown in Fig. 6. Comparison of these findings with other studies where a new snow scheme was introduced reiterates that the model representation of snow strongly affects soil temperatures (Gouttevin et al., 2012; Wang et al., 2013). Reliable soil temperature simulations are essential to study the permafrost–climate feedback. Biskaborn et al. (2019) concluded that recent warming trends of permafrost soils are partly due to an increase in snow insulation, accelerating its degradation. Both field observations and modelling studies have identified this close link between snow and permafrost conditions (Johansson et al., 2013; McGuire et al., 2016; Lawrence and Slater, 2010). Identifying changes in permafrost-underlain areas is important because of the potential increase in organic matter decomposition and release of greenhouse gases. These aspects will be further evaluated in LPJ-GUESS with the new snow scheme.

We observed a general decrease in mean NPP during the winter and a marginal difference in the summer. Considering the presence of permafrost, however, we noted an increase in GPP and NPP for non-permafrost-underlain areas in summer. The significantly warmer winter soil conditions for the Dynamic scheme caused an increase in heterotrophic respiration – i.e due to faster litter decomposition rates and increased microbial activity. Accordingly, soil respiration increased during the winter in the non-permafrost region. During the summer, there is an overall minor decrease in soil respiration due to the lower soil temperatures simulated by the Dynamic scheme. The net effect of the above-discussed processes is an overall increase in carbon emissions during the winter and an increased uptake during the summer.

The impact of the new snow scheme on summer conditions was surprising. These differences were caused by the changes in spring snow and soil temperatures and soil water availability. During springtime, soils with the Static scheme warm more quickly, due to the lower insulation, which leads to an earlier thaw and increased soil water availability. The Dynamic scheme simulates a more realistic atmosphere–snow–soil heat transfer, leading to a slower temperature transition. The difference between the schemes diminishes towards the end of the summer. Overall, the simulated pan-arctic carbon fluxes are systematically lower than other published values (Efren et al., 2019; Rawlins et al., 2015; McGuire et al., 2012). Virkkala et al. (2021) estimate an annual NEE in the range of −46 to +10 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, while our simulations have a mean estimate around −35.5 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (see Table 2).

Besides the carbon fluxes, we also evaluated the simulated annual soil and vegetation carbon pools. Vegetation pools were more different when applying the Dynamic, while no clear differences were apparent for soil carbon (see Table 2). With the Dynamic snow scheme, the soil carbon pool is lower within the permafrost region and higher outside of the permafrost region. These results align well with the sensitivity study by Gouttevin et al. (2012), which highlighted that decreased soil carbon stocks can be attributed to a higher respiration rate and increased microbial decomposition rates.

The mean simulated soil carbon content was around 10 kg C m⁻², which is much lower than the 50–100 kg C m⁻² range suggested in literature (Hugelius et al., 2014; Hugelius, 2012). This difference is most likely due to the fact that we only simulated grid cells with upland soils, while peatlands were not represented. The inclusion of peatlands would have led to a larger amount of soil carbon, since these ecosystems are characterized by waterlogged soils in which decomposition is suppressed – although carbon can be released as methane. Also, all organic matter is considered to be in the top 0.5 m of the soil in the current version of LPJ-GUESS and is therefore only affected by the average soil temperature and moisture conditions down to 0.5 m, but not by conditions further down. These aspects will be taken into account in ongoing model development. Our analysis highlights that the observed differences between the Static and Dynamic schemes correlate well with the spatial pattern of near-surface soil temperature changes. This shows that the changes in soil temperature influence the soil carbon content in the model. The shortcomings in soil carbon simulation will be addressed and improved in the future, which will enable a more reliable carbon pool assessment.

4.3 Impact on vegetation dynamics

Satellite-based studies have identified an overall greening trend across the Arctic, in response to a warming from the 1980s until now. However, they also showed that this greening trend is not uniform and certain areas have actually experienced browning (loss of greenness) during this period (Berner et al., 2020; Myers-Smith et al., 2020). This may be partly due to damage to vegetation following extreme winter events (Phoenix and Bjerke, 2016). At the site level, a recent study by Niittynen et al. (2020) showed that winter thermal conditions are a strong control on vegetation patterns in arctic landscapes. Still, it is challenging to fully understand vegetation responses to warming solely from remotely sensed data or field observations, due to the scale dependency of interpreting trends in vegetation dynamics. Moreover, most field sites are highly concentrated in northern Scandinavia and Alaska, which leaves the full heterogeneity of the arctic and its ecosystems vastly under-sampled (Metcalfe et al., 2018). With ecosystem models, we can fill in spatial gaps, identify feedback loops, and assess potential future changes.

Following the assessment of the new snow scheme's impact on biogeochemical variables, we compiled the simulated vegetation conditions with the two snow schemes. We found that PFT domination changed marginally using the updated snow scheme in some grid cells. The main direction of change is a grass-to-shrub dominance shift in grid cells – shown in Fig. S11. The forest–shrub border did not shift much in most areas. However, the vegetation carbon pool was higher with the Dynamic scheme, which indicates that even though the changes are minor visually, they affected vegetation biomass. These comparisons show that changing the snow scheme in LPJ-GUESS affected vegetation distribution and composition, albeit on a small scale.

4.4 Outlook

It is well established that the Arctic is highly susceptible to climate change, and the ongoing warming has significant consequences for the arctic system – even if we implement the most strict mitigation measures (Bruhwiler et al., 2021; AMAP, 2017). One of these consequences is a change in snow conditions. In the near future, snow thickness will decrease, caused by air temperature and precipitation changes, inducing a decrease in snow-covered area in the region (AMAP, 2017; IPCC, 2014). Due to a later onset and earlier spring melt, the snow season is expected to shorten under a changing climate. Moreover, northern high latitudes are predicted to be rain-dominated in the future (Bintanja and Andry, 2017; Johansson et al., 2011). These changes will strongly influence soil thermodynamics, and the observed and projected changes will have a significant impact on arctic ecosystems (Bruhwiler et al., 2021). To be able to provide robust projections of the future, we need to account for a multitude of interlinked processes and feedbacks. Some of the current key areas are to assess the relative sensitivity of plant productivity to climate change, the development of decomposition rates, and their net effect on the carbon budget.

Besides an assessment of geophysical and biogeochemical processes, LPJ-GUESS can also be used to explore future vegetation trends to assess whether favourable growing conditions will induce further greening or whether new stressors will prompt local- to regional-scale browning (loss of productivity) (Myers-Smith et al., 2020). Studying future vegetation trends across the Arctic is important from a global perspective. A potential decrease in snow-covered area may significantly decrease surface albedo, which would enhance arctic warming. Consequently, changes in snow dynamics on a local scale influence carbon fluxes by altering soil thermal conditions and vegetation habitat. Evaluating snow–soil–vegetation feedbacks in future studies is therefore relevant to further investigate climate change impacts on the Arctic in global-scale land surface modelling.

5 Conclusions

This study shows that the representation of snow dynamics in a dynamic vegetation model significantly influences the simulated soil thermodynamics and related biogeochemical variables. We show, due to the improved snow insulation capacity, that the new Dynamic snow scheme simulates more realistic soil thermodynamics and permafrost extent than the old Static scheme. The improved simulation of permafrost cover can be attributed to significantly warmer winter soil temperatures, which compare well to observations across 256 locations in Russia. We further showed the importance of an accurate snow scheme for the simulation of biogeochemical processes. Our results show that the intermediate-complexity snow scheme had a significant impact on carbon fluxes. Heterotrophic respiration increased during the winter, which led to an increased carbon release during the cold season. We also identified differences in soil carbon content between the Static and Dynamic simulations. Although the modelled soil carbon content was lower than literature values, the spatial pattern of low and high soil carbon content aligns well with observations. A differentiation between the seasons and accounting for permafrost presence highlighted the differences between the two sets of simulations. Wintertime carbon emissions were higher using the Dynamic scheme, both within but especially outside of the permafrost region. The differences between the simulations were larger within permafrost-underlain areas for the physical variables. Besides spatial patterns, we explored seasonal differences, which showed that summertime conditions were also affected by the representation of snow. In contrast to warmer soils in winters, soils were cooler in summer using the Dynamic scheme – especially in permafrost-underlain areas – due to a delayed response to snowmelt. These differences between the old and new snow schemes underline the importance of further developing winter processes as they may significantly affect the annual carbon budget.

These findings contribute to our understanding of the impact of wintertime changes on the arctic carbon cycle. We show that an accurate, dynamic snow scheme is essential to investigate the full complexity of snow–soil–vegetation relationships. Models are valuable tools to aid our understanding of large-scale climate change impacts due to the sparse availability of observations in the Arctic. Addressing identified knowledge gaps in models is imperative to decrease the uncertainty around carbon balance estimates. Due to the large spread of observed and modelled seasonal and inter-annual cycles of carbon fluxes, it is not yet possible to determine with high certainty whether the Arctic will act as a carbon source or sink in the future (Fisher et al., 2014; McGuire et al., 2018). To decrease uncertainty in simulations, contemporary modelling efforts are directed, on the one hand, at model inclusion (account for key, still missing processes) and, on the other hand, at refining process formulations using observational data (McGuire et al., 2012; Fisher et al., 2014, 2018).

In this study, we aimed to improve the representation of the cold-season process using non-growing-season observations and findings. This enhances the versatility and applicability of LPJ-GUESS as a tool to address the remaining uncertainties regarding climate change impacts at northern high latitudes and its consequences on a global scale. With this model, we have the ability to investigate complex ecosystem interactions under changing environmental conditions at multiple scales, considering nitrogen cycling, permafrost processes (freeze–thaw cycles, hydrology), stochastic vegetation dynamics, and also potential land cover and land use changes. Realistic soil temperature simulations are the first step to improve the simulation of greenhouse gas emissions under different climate scenarios across the Arctic (Natali et al., 2019). Our results show that by improving a process that appears only relevant in winter, such as snow, we not only decrease the uncertainty regarding physical and biogeochemical parameters during the cold season, but we also improve simulations of soil conditions and the carbon cycle in the growing season. Further developments will aim at improving soil carbon content simulations and better assessing plant responses to future environmental conditions during the cold season. By accounting for snow–soil–vegetation interactions in all seasons of the year, we ensure more reliable projections of the future state of vegetation composition, permafrost stability, and greenhouse gas exchange in a rapidly warming Arctic.

Appendix A: Simulation details

Table A1 Table with used variables, their description, and units.