Variability of the Surface Energy Balance in Permafrost Underlain Boreal Forest

Boreal forests in permafrost regions make up around one-third of the global forest cover and are an essential component of regional and global climate patterns. Further, climatic change can trigger extensive ecosystem shifts such as the partial disappearance of near surface permafrost or changes to the vegetation structure and composition. Therefore, our aim is to understand how the interactions between the vegetation, permafrost, and the atmosphere stabilize the forests and the underlying permafrost. Existing model set-ups are often static or are not able to capture important processes such as the vertical struc5 ture or the leaf physiological properties. There is a need for a physically based model with a robust radiative transfer scheme through the canopy. A one-dimensional land surface model (CryoGrid) is adapted for the application in vegetated areas by coupling a multilayer canopy model (CLM-ml v0) and is used to reproduce the energy transfer and thermal regime at a study site (N 63.18946◦, E 118.19596◦) in mixed boreal forest in Eastern Siberia. We have in-situ soil temperature and radiation measurements, to evaluate the model and demonstrate the capabilities of a coupled multilayer forest-permafrost model to in10 vestigate the vertical exchange of radiation, heat, and water. We find that the forests exert a strong control on the thermal state of permafrost through changing the radiation balance and snow cover phenology. The forest cover alters the surface energy balance by inhibiting over 90% of the solar radiation and suppressing turbulent heat fluxes. Additionally, our simulations reveal a surplus in longwave radiation trapped below the canopy, similar to a greenhouse, which leads to a comparable magnitude in storage heat flux to that simulated at the grassland site. Further, the end of season snow cover is three times greater at the forest 15 site and the onset of the snow melting processes are delayed. 1 https://doi.org/10.5194/bg-2020-201 Preprint. Discussion started: 18 June 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
Around 80% of the world's boreal forest occurs in the circumpolar permafrost zone (Helbig et al., 2016). Despite little human interference, and due to extreme climate conditions such as winter temperatures below −50 • C and very low precipitation, the biome is highly sensitive to climatic changes (ACIA, 2005;AMAP, 2011;IPCC, 2014) and thus prone to vegetation shifts. Boreal forest regions are expected to warm by 4 to 11 • C by 2100, coupled with a more modest precipitation increase 5 (IPCC, 2014;Scheffer et al., 2012). Moreover, during 2007-2016 continuous zone permafrost temperatures have increased by 0.39 (±0.15) • C (Biskaborn et al., 2019;IPCC, 2019). The northeastern part of the Eurasian continent is dominated by such vast boreal forest -the taiga. Due to its sheer size, the biome is not only sensitive to climatic changes, but also exerts a strong control on numerous climate feedback mechanisms through the altering of land-surface reflectivity, the emission of biogenic volatile organic compounds and greenhouse gases, and the transfer of water to the atmosphere (Bonan et al., 2018;10 Zhang et al., 2011). The forests are usually considered to efficiently insulate the underlying permafrost (Chang et al., 2015).
The canopy exerts shading by reflecting and absorbing most of the downward solar radiation, changes the surface albedo and decreases the soil moisture by intercepting precipitation and increasing evapotranspiration (Vitt et al., 2000). Additionally, the forest promotes the accumulation of an organic surface layer which further insulates the soil from the atmosphere (Bonan and Shugart, 1989). Changing climatic conditions can promote an increasing active layer depth or trigger the partial disappearance determine the energy transfer between the top of the canopy atmosphere and the ground. Therefore, there is a need for a physically based model with a robust radiative transfer scheme through the canopy for a detailed analysis of the vegetation's impact on the hydro-thermal regime of the permafrost ground below.
With a tailored version of a one-dimensional land surface model (CryoGrid, Westermann et al. (2016)) we perform and analyze numerical simulations and reproduce the energy transfer and surface energy balance in permafrost underlain boreal 5 forest of Eastern Siberia. CryoGrid has, so far, not included a vegetation scheme but has been used to successfully describe atmosphere-ground energy transfer and the ground thermal regime in barren and grass-covered areas Westermann et al., 2016;Nitzbon et al., 2019Nitzbon et al., , 2020. In our study, we have adapted a state-of-the-art multilayer vegetation model (CLM-ml v0, originally developed for the Community Land Model CLM by Bonan et al. (2018)). We tailor and implement this scheme to simulate the turnover of heat, water, and snow between atmosphere, forest canopy and ground. We take 10 advantage of a detailed in-situ data record from a study site, which is used to provide model parameters, as well as for model validation by comparing field measurements with simulation results. The main objectives of this study are 1. To demonstrate the capabilities of a coupled multilayer forest -permafrost model to simulate vertical exchange of radiation, heat, and water for boreal forests.
2. To investigate the impact of the detailed canopy module on the surface energy balance of the underlying permafrost at a 15 mixed boreal forest site in Eastern Siberia.

Study Area
The study site is located south east of Nyurba (N 63.18946 • , E 118.19596 • ) in a typical boreal forest zone intermixed with some grassland for horse grazing and shallow lakes. The forest is rather dense and mixed, with the dominant taxa evergreen 20 spruce (Picea obovata, 92%), deciduous larch (Larix gmelinii, 7.3%) and some hardwood birch (Betula pendula, <1%). The average tree height is 5.5 m for spruce and 12 m, respectively. These boreal forest environments experience 6 to 8 months of freezing temperatures reaching extremes of −62 • C in winter and up to 35 • C between May and September. The low annual average temperatures result in continuous permafrost and therefore poorly drained, podzolized and nutrient-poor soils (Chapin et al., 2011). The summer precipitation has decreased over the past 30 years (Hayasaka, 2011). The temperature trend from 25 1970 to 2010 for the Central Yakutian region is positive for spring, summer and fall and negative for winter (monthly surface temperature quantified using Climatic Research Unit (CRU) TS4.01 data (Harris et al., 2014;Stuenzi and Schaepman-Strub, 2020)). The treeline of Northern Siberia is dominated by the deciduous needleleaf tree genus Larix Mill. up to N 72.08 • . Larix sibirica Ledeb. from E 60 − 90 • , Larix gmelinii Rupr. between E 90 − 120 • and Larix cajanderi Mayr. from E 120 − 160 • (see Fig. 1). Larch compete effectively with other tree taxa because of its deciduous leaf habit and dense bark. In more southern 30 margins of Eastern Siberia, such as our study area, larch is mixed with evergreen conifers (Siberian pine (Pinus sibirica, Pinus sylvestris), spruce (Picea obovata Ledeb.), and fir (Abies sibirica)), and hardwood (Betula pendula Roth., B. pubescence Ehrh., Populus tremula L.) (Kharuk et al., 2019). Moreover, the ground vegetation is poor in diversity and dominated by mosses and lichens that form carpets. Larch has shallow roots, and grows on clay permafrost soils with an active layer of around 0.7 m and a maximum wetness of 20-40 %. Evergreen conifers and hardwood both prefer deeper active layers and a higher soil moisture availability (Ohta et al., 2001;Furyaev et al., 2001;Rogers et al., 2015). to simulate the evolution of the snow cover including the Crocus snowpack scheme as implemented by Zweigel et al. (2020).
The model is forced by standard meteorological variables which may be obtained from AWSs, reanalysis products, or climate models. The required forcing variables include air temperature, wind speed, humidity, incoming short-and longwave radiation, air pressure and precipitation (snow-and rainfall) . The change of internal energy of the subsurface domain over time is controlled by fluxes across the upper and lower boundaries written as where the input to the uppermost grid cell is derived from the fluxes of shortwave radiation (S in , S out ) and longwave (L in , L out ) radiation, at the same time regarding the latent (Q e ), sensible (Q h ), sensible heat added by precipitation (Q precip ) and storage heat flux (Q s ) between the atmosphere and the ground surface .
For this study, we have adapted the multilayer canopy model developed by Bonan et al. (2014Bonan et al. ( , 2018. The canopy model is 10 coupled to CryoGrid by replacing its standard surface energy balance scheme while soil state variables are passed back to the forest module. In the following, we describe the canopy module and its interaction with the existing CryoGrid soil and snow scheme. All differences towards former CryoGrid parameterizations are summarized in Table 1. The multilayer canopy model (Bonan et al., 2018) was developed based on the use with CLM soil properties. Following the notations summarized in Bonan (2019) we developed a CLM-independent multilayer canopy module which can be coupled to CryoGrid by integrating novel interactions and an adapted snow cover parameterization. In order to set necessary parameters of the canopy module we make use of values defined for the plant functional type deciduous needleleaf forest of CLM. Please note that all parameters defining the canopy are set as constant values so that vegetation is not dynamic and changes in forest 5 composition or actual tree growth are not considered in this study.

Snow module
The snow module employed in this study is based on Zweigel et al. (2020) (submitted, available upon request). Zweigel et al.
(2020) extend the CryoGrid model with a snow microphysics parameterizations based on the CROCUS snow scheme (Vionnet et al., 2012), as well as lateral snow redistribution. The CLM-ml (v0) multilayer canopy model has not yet been coupled to 10 a snow scheme (Bonan et al., 2018). Following Vionnet et al. (2012) the microstructure of the snow pack is characterized by grain size (gs, mm), sphericity (s, unitless, range 0-1), and dendricity (d, unitless, range 0-1). Fresh snow is added on top of the existing snow layers with temperature and windspeed dependent density and properties. After deposition the development of each layers microstructure occurs based on temperature gradients and liquid water content (Vionnet et al., 2012). Snow albedo for the surface layer and an absorption coefficient for each layer are calculated based on the snow properties. Solar radiation 15 is gradually absorbed throughout the snow layers and the remaining radiation is added to the lowest cell. Additionally, the two mechanical processes of mechanical settling due to overload pressure and wind compression increase snow density and compaction. During snowfall, new snow is added to the top layer in each timestep and mixed with the old snow based on the amount of ice. Once a cell exceeds the snow water equivalent of 0.01 m, which equals a snow layer thickness of 0.03 m, a new snow layer is built. Here, snow accumulates on the ground under the forest canopy. During the first snowfall, the surface 20 energy balance of the ground and snow is calculated for each respective cover fraction. After reaching a snow layer thickness of 0.03 m, the ground surface energy balance is calculated for the snow-pack itself (Table A2). Variables exchanged based on the snow cover are ground surface temperature, surface thermal conductivity and layer thickness of the layer directly under vegetation. Evaporation flux is subtracted from the snow surface. Top of the canopy wind speed is used to calculate the density of the falling snow. Additionally, snow interception is handled like liquid precipitation interception described in Equation 6.
where L incanopy and S incanopy are the incoming long-and shortwave radiation at the ground surface and the lower boundary fluxes of the multilayer canopy module, and S out ground is the outgoing shortwave flux from the ground, L out ground the outgoing longwave flux, Q h ground the sensible heat flux, Q e ground the latent heat flux, and Q s ground the storage heat flux at the ground surface. The first six components of the sub-canopy energy balance directly replace the respective components surface energy balance scheme of CryoGrid (Equation 2. Whereas Q s ground is calculated based on temperatures of the uppermost ground or 5 snow layers that are passed from CryoGrid to the forest module. The storage heat flux is calculated as where k is the soil thermal conductivity, T s the soil surface temperature, T ground the actual ground temperature for the first layer below the surface and ∆z the layer thickness. In Fig. 2 the energy fluxes expected for a forested and a grassland site in the snow-covered and snow-free periods each are illustrated schematically.  In the novel model set-up which allows for soil-vegetation interaction, the vegetation module receives ground state variables of the top 0.7 m of the soil layers. These state variables are soil layer temperature (T ground ) and soil layer moisture (W ground ), as well as the diagnostic variables soil layer conductivity (k ground ), and ice content (I ground ). The vegetation transpiration fluxes are subtracted from the ground soil layers within the rooting depth and evaporation fluxes from the ground surface.
Following the notation in Bonan et al. (2018) the rain and snow fraction reaching the ground (W grounds ) is described as follows and consist of the direct throughfall (f P R ), the canopy drip (D c ), the canopy evaporation (E c ), the stemflow (D t ) and the stem evaporation (E t ), which are based on the retained canopy water (W c ) where 1 − f − f t P R is the intercepted precipitation, and the retained trunk water (W t ) respectively.

Model setup and simulations
We run model simulations for forested and non-forested scenarios based on in-situ measurements recorded in 2018 and 2019.
The subsurface stratigraphies used in CryoGrid is described by the mineral and organic content, natural porosity, field capacity 15 and initial water/ice content. Some of these parameters could be measured at the forest and grassland sites and were used to set the initial soil profiles in the model. Table A3 summarizes all parameter choices for soil stratigraphies and Table A4 summarizes constants used. The subsurface stratigraphy extends to 100 m below surface where the geothermal heat flux is set parameters were adopted from previous studies using CryoGrid (Table A2) (Langer et al., 2011a(Langer et al., , b, 2016Westermann et al., 2016;Nitzbon et al., 2019Nitzbon et al., , 2020. We use ground surface temperature (GST) as the target variable for model validation. GST results from the surface energy balance at the interface between canopy, snow cover, and ground and provides an integrative measure of the different model components. In addition it is the most important variable determining the thermal state of permafrost. 25 For the canopy stratigraphy, we follow the parameterizations in Bonan et al. (2018) for the plant functional type deciduous needleleaf (see Table A5 and A6). This canopy stratigraphy can be described by two parameters: the leaf area index (LAI) measured at the bottom of the canopy defines the total leaf area. The leaf area density function on the other hand describes the foliage area per unit volume of canopy space which is the vertical distribution of leaf area. Leaf area density is measured by evaluating the amount of leaf area between two heights in the canopy separated by the distance. This function can be expressed by the beta distribution probability density function which provides a continuous representation of leaf area for the use with multilayer models (Bonan (2019) for further information). Here, we use the beta distribution parameters for needleleaf trees (p = 3.5, q = 2) which resembles a cone-like tree shape. LAI can be estimated from satellite data, calculated from below-canopy 3 Results

Model validation and in-situ measurements
The model is validated against ground surface temperature (GST) measurements of forested and non-forested study sites. The  More detailed insights into differences of available radiation at the ground surface are presented in Fig. 4. Here, the incoming short-and longwave radiation measured at the grassland site and modeled for the forest and grassland are shown. The longwave radiation dominates the incoming part of the radiation balance at both sites throughout the year. In the snow-free period downward longwave radiation flux is 64.0 W m −2 higher at the forest ground. In the snow-covered period downward longwave radiation flux is 39.5 W m −2 higher at the forest ground. This results in a surplus of energy under the forest canopy. The snow-covered periods), showing that the canopy successfully intercepts (absorbs and reflects) most of the incoming shortwave radiation. Snow-free shortwave down at the grassland site is more than 12 times higher in the snow-free period and 16 times higher in the snow-covered-period.

Thermal regime of the ground near the surface
In a second step, we compare the annual, snow-free and snow-covered period average GST to understand the overall model 5 performance and the relative temperature differences between the forest and grassland sites (Fig. 5).
14 https://doi.org/10.5194/bg-2020-201 Preprint. Discussion started: 18 June 2020 c Author(s) 2020. CC BY 4.0 License. Measured and modeled average GST values are summarized in Table 2. The highest deviation between modeled and measured GST is found at the grassland site. Here, the model shows a cold bias of −4.1 • C. Also, there is a warm bias of 0.9 • C in the snow-free period in the forest. Overall, we find an average annual difference of 0.9 • C between the two sites. For snow-free season this difference is 2.2 • C and snow-covered 2.1 • C, respectively. For a more detailed understanding of the annual cycle of the thermal evolution of permafrost ground at our study sites we 5 compare the weekly averaged GST at the grassland and forest site (see Fig. 6). The more detailed analysis of the annual cycle 15 https://doi.org/10.5194/bg-2020-201 Preprint. Discussion started: 18 June 2020 c Author(s) 2020. CC BY 4.0 License. reveals periods with distinct differences between the model simulation and the measured values. For both study sites the model produces a slight GST overestimation in summer and a prolonged thawing period in spring. The measured data shows a much faster ground warming in spring. This difference is over 20 days at the forest site and 15 days later at the grassland site. In addition, there is a cold bias by 5 • C in January at the grassland site. This bias is not seen at the forest site. Thawing starts later in model simulations than measured. To further understand the temporal evolution of the permafrost ground, we compare the modeled and measured active layer thickness at both study sites. In the grassland, the modeled maximum active layer thickness (ALT) is 2.35 m between 13. and

Snow cover assessment
To analyze differences in the snow cover evolution, we compare modeled and measured snow depth at the AWS and the modeled snow depth in forest. Snow depth modeled in grassland agrees well with the measured snow depth and reaches a maximum value of 0.26 m in late April (see Fig. 7). Towards spring, the snow pack in the forest accumulates to 1.2 m and snow melting starts around the same period but lasts much longer up until the end of May. To understand the high impact that coupling the vegetation has on the snow cover, we next look at the surface energy balance simulated by our model for our specific study site and a hypothetical, sparse canopy with LAI = 1 (see Fig. 8). Incoming solar radiation at the ground is 5 times higher in the sparse canopy simulation, which leads to a higher net radiation flux (Q net ) and higher snow melt rates. Turbulent fluxes (Q h and Q e ) are similar, which suggests that air circulation is blocked, even in a very sparse canopy. The high longwave radiation at the forest ground is persistent for the hypothetical sparse canopy as well, and

Discussion
The model presented here is found to be capable of simulating the ground thermal regime of permafrost underneath boreal forests. The implemented scheme is able to simulate the physical processes that define the vertical exchange of radiation, heat, water and snow between permafrost and canopy. Our simulations show that the forests exert a strong control on the thermal state of permafrost. At the grassland site, we find a much larger ground surface temperature (GST) amplitude of 60.35 • C over 5 the annual cycle, which is 19 • C higher than at the forest site. This vegetation dampening effect on soil temperature is welldescribed in literature (Oliver et al., 1987;Balisky and Burton, 1993;Chang et al., 2015). Earlier work by Bonan and Shugart (1989) found that forest soils generally thaw later and less deeply, and are cooler than in open areas. In the winter, forested soils are typically warmer relative to open areas. The tree cover can maintain stable permafrost under otherwise unstable thermal conditions (Bonan and Shugart, 1989). Our results are in agreement with these observations, but further demonstrate that the 10 impact of mixed boreal forest on the GST is strongest during the snow period and the summer peak with the warmest months.
Our model reveals an average of 6.5 • C higher GST during the snow-covered period and 1.5 • C lower GST during the snowfree period. Our model simulations unravel that the strong control on the thermal state of permafrost is a result of the combined effects of canopy shading, suppression of turbulent heat fluxes, below-canopy longwave enhancement, increased soil moisture 18 https://doi.org/10.5194/bg-2020-201 Preprint. Discussion started: 18 June 2020 c Author(s) 2020. CC BY 4.0 License. and distinct snow cover dynamics. These relevant processes controlling forest insulation will be discussed individually in the following subsections. Then, we will present a detailed discussion on the model applicability and limitations of the model.

Canopy shading and longwave enhancement
The surface energy balances simulated by the model are very different for grassland and forest. The forest canopy reflects and absorbs over 92% of incoming solar radiation for both snow-free and snow-covered periods. The forest ground albedo therefore 5 has little influence on the energy balance. This canopy shading effect makes the longwave radiation the only source of radiative energy at the forest site for both time periods. A surplus of longwave radiation by over 20% is largely trapped below the canopy due to extremely low turbulent heat fluxes, similar to a greenhouse. The increased longwave radiation results in a relatively strong storage heat flux. Despite shading, the storage heat flux at the forest site is similar in magnitude to that simulated at the grassland site. This explains the small difference in the modeled depths of the active layer at both sites. The heat flux plate 10 used for measuring the conductive heat flux at the grassland site is designed for use in mineral soil. Due to a certain amount of organic content in the upper soil layer under-or overestimated heat fluxes are possible (Ochsner et al., 2006). This could explain to some extent the difference between measurements and simulations. We further find, that during the snow-free period, the sensible and latent heat flux at the canopy top are high while being close to zero at the forest ground. In spring, we see a high negative latent heat flux at the forest ground which points at a condensation event. This might be due to a still existing 15 snow layer that suddenly gets exposed to warmer air. Further, in Fig. 8 we see that the greenhouse effect caused by longwave enhancement is similar for a canopy with a much lower LAI and turbulent heat fluxes are also small. In summery, we show that the canopy effectively absorbs and reflects the majority of incoming solar radiation, making canopy shading one of the main controlling mechanisms and that the canopy enhances the longwave radiation the forest ground, because of extremely low turbulent heat fluxes.

Soil moisture, canopy interception and evapotranspiration
According to the majority of studies, tree growth in permafrost areas is limited by summer air temperatures and available water from snow melt, water accumulated within the soil in the previous year and permafrost thaw water (Kharuk et al., 2015;Sidorova et al., 2007). The amount of precipitation in the eastern Siberian Taiga is characteristically small compared to other areas, therefore it is expected that permafrost plays an important role in the existence of these forests (Sugimoto et al., 2002). 25 Sugimoto et al. (2002) found that plants used rainwater during wet summers, but melt water from permafrost during drought summers. This indicates that permafrost provides the direct source of water for plants in drought summers and retains surplus water in the soil until the next summer. They conclude that if this system is disturbed by future warming, the forest stands might be seriously damaged in severe drought summers (Sugimoto et al., 2002). Our grassland site, which was supposedly forested until the 1950s, has dried up and has a much smaller organic layer and a maximum active layer thickness of 2.30 m. The 30 volumetric water content in forest soil is 10-20% higher despite the same amount of precipitation and a higher evaporative flux during the growing season. This points to the conclusion that permafrost plays an important role in regulating the hydrological conditions in this boreal forest area and holds a certain amount of ground water that is later available to plants.

Applicability and model limitations
The presented model accurately simulates the recorded GST. The detailed analysis of the annual cycle shows that the snow melt period in spring is biased at the forest and grassland sites. In reality, the ground warms up faster than modeled and in the forest this might be due to the high snow cover. Forest snow depth and density measurements would be highly desirable for future work. Further, snow redistribution is not included in our model, but could be an interesting parameter to look at.

5
Additionally, one aspect not represented in the model is the moisture transport and migration in frozen ground or the forming of ice lenses. These processes could lead to the observed, lower ALT at the forest site. Further, actual plant area index values based on collected field data could be incorporated in the future for a more detailed understanding of the exact canopy density, but were out of scope for this particular study. An existing, extensive tree survey data set could be used to calculate the actual plant area indices in the proximity of our study site. A second canopy aspect not considered here, is the reduced leaf area 10 index of deciduous larch stands outside of the comparably short growing season. Based on the ALT point measurements at the forest site it is clear that our model overestimates the ALT in the forest. One feasible explanation is the overestimation of the snow cover which leads to higher winter insulation. This could be due to a missing LAI constraint affecting deciduous taxa in the snow-covered period. The model used here is not particularly parameterized for the application in deciduous larch stands, and we do not have snow depth measurements in the forest. A complete removal of the leaf area mass would lead 15 to higher turbulent fluxes, higher solar radiation, hence higher melt rates and more snow compaction that would potentially counteract the greenhouse effect and snow dynamics described. Our forest site is partly covered by deciduous larch trees (7%), which cause a lower LAI in winter and therefore higher solar radiation before leaf-out in spring. This could change the surface energy balance at the forest ground. In our simulations under a hypothetical, sparse canopy, we find a much lower snow pack height, with a lower insulation effect against the extreme cold and an earlier onset of snow melting processes. In 20 combination, this causes a much lower ALT. To simulate the needle-tossing of deciduous larch, a leaf area index threshold could be incorporated for a certain time period between needle tossing and leaf-out in spring. This tunes the model towards a more detailed representation of highly sensitive larch-dominated forests, which are particular to large parts of Eastern Siberia.
This modification could reveal further taxa-specific interactions with permanently frozen ground. It would also be desirable to implement a spatially-explicit, dynamic vegetation model, such as the larch forest simulator (LAVESI, Kruse et al. (2016)) 25 to further analyze the dynamic vegetation distribution under the recognition of the found interactions. This would allow us to simulate the vegetation response to changes in permafrost temperature and hydrology dynamically over a large timescale, and across a wide range of boreal forest ecosystems in Eastern Siberia. In combination, the above model limitations could explain a great part of the GST and ALT differences between measurements and modeled simulations.

30
This study presents a specific application of a coupled multilayer forest-permafrost model to investigate the energy transfer and surface energy balance in permafrost underlain boreal forest of Eastern Siberia. The comparison of measured and modeled GST at a mixed forest and a grassland site, as well as the comparison of the modeled and measured radiation fluxes at the grassland site, justify the use of the physically-based modeling approach to investigate the thermal regime and surface energy balance in this complex ecosystem. The detailed vegetation model successfully calculates the canopy radiation and water budgets, leaf fluxes, as well as canopy turbulence and aerodynamic conductance. These canopy fluxes alter the below-canopy surface energy balance, the ground thermal conditions and the snow cover dynamics. We find a strong dampening effect of 19 • C on the annual ground surface temperature amplitude of the permafrost. Further, forested permafrost maintains a higher soil water 5 content by controlling water storage in the ground. The forest cover alters the surface energy balance by inhibiting most of the solar radiation and suppressing turbulent heat fluxes. Additionally, we reveal that the canopy leads to a surplus in longwave radiation trapped below the canopy, similar to a greenhouse. Therefore, and despite the canopy shading, the storage heat flux at the forest site is similar in magnitude to that simulated at the grassland site. Further, the end of season snow cover is much greater at the forest site and the onset of the snow melting processes are delayed. In summary, we identify the following key 10 points.
i. The canopy effectively absorbs and reflects over 90% of incoming solar radiation, making canopy shading one of the main controlling mechanisms.
ii. The vegetation cover suppresses the majority of the turbulent heat fluxes in the below-canopy space.
iii. The canopy enhances the longwave radiation below the canopy by up to 20%, similar to a greenhouse, which results in a 15 comparable magnitude of storage heat flux for both sites.
iv. Forested permafrost holds a higher ground water content than the dry grassland site.
v. Canopy shading leads to slower snow melting, less snow compaction, and therefore a higher snow pack.
vi. The differences in the thermal development of the forest and grassland sites are highly influenced by the depth, density and the resulting insulation capacities of the snow cover, which is in turn controlled by the canopy density.

20
Code and data availability. With the AWS equipped as a Bowen ratio station, B is calculated following (Foken, 2016) as where the specific heat at constant pressure for moist air (c p ) is 1.006 kJ kg -1 K -1 and the latent heat of vaporization of water (L v ) is 2260 kJ kg -1 . With ∆T being the temperature difference between the two air temperature sensors at heights 0.115 m and 5 0.252 m. ∆q is the difference in specific humidity calculated from measured relative humidity (φ), temperature and pressure.
Thereafter, latent heat flux (Q e ) is calculated as and sensible heat flux (Q h ) 10 with the storage heat flux (Q s ) and Q g as the convective ground heat flux.