Sensitivity of ecosystem-protected permafrost under changing boreal forest structures

Boreal forests efficiently insulate underlying permafrost. The magnitude of this insulation effect is dependent on forest density and composition. A change therein modifies the energy and water fluxes within and below the canopy. The direct influence of climatic change on forests and the indirect effect through a change in permafrost dynamics lead to extensive ecosystem shifts such as a change in composition or density, which will, in turn, affect permafrost persistence. We derive future scenarios of forest density and plant functional type composition by analyzing future projections provided by the dynamic global vegetation model (LPJ-GUESS) under global warming scenarios. We apply a detailed permafrost-multilayer canopy model to study the spatial impact-variability of simulated future scenarios of forest densities and compositions for study sites throughout eastern Siberia. Our results show that a change in forest density has a clear effect on the ground surface temperatures (GST) and the maximum active layer thickness (ALT) at all sites, but the direction depends on local climate conditions. At two sites, higher forest density leads to a significant decrease in GSTs in the snow-free period, while leading to an increase at the warmest site. Complete forest loss leads to a deepening of the ALT up to 0.33 m and higher GSTs of over 8 ∘C independently of local climatic conditions. Forest loss can induce both, active layer wetting up to four times or drying by 50%, depending on precipitation and soil type. Deciduous-dominated canopies reveal lower GSTs compared to evergreen stands, which will play an important factor in the spreading of evergreen taxa and permafrost persistence under warming conditions. Our study highlights that changing density and composition will significantly modify the thermal and hydrological state of the underlying permafrost. The induced soil changes will likely affect key forest functions such as the carbon pools and related feedback mechanisms such as swamping, droughts, fires, or forest loss.


Introduction
The boreal forest cover exerts a strong control on numerous climate feedback mechanisms (Bonan et al 2018. Globally, 80% of boreal forests are underlain by permafrost (Helbig et al 2016). The forest cover is considered to efficiently insulate the underlying, ecosystem-protected permafrost (Chang et al 2015) and therefore play an important role in the development of boreal regions and the stability of permafrost in a warming climate. Boreal regions are projected to warm between 4 • C and 11 • C by 2100, with a modest precipitation increase (Scheffer et al 2012, Meredith et al 2019. The change in air temperature and precipitation directly influences the vegetation cover development (Esper et al 2010, Kharuk et al 2015, Sato et al 2016, Ito et al 2020 and permafrost thaw (Meredith et al 2019), directly affecting soil water availability and root space limitation (Carpino et al 2018). The changing thermohydrological soil conditions may provoke changes in forest density and forest composition (Takahashi 2006, Kharuk et al 2013, Liu et al 2017, Kropp et al 2021 leading to extensive ecosystem shifts (Pearson et al 2013, Gauthier et al 2015, Kruse et al 2016. Forest composition and density exert a strong control on permafrost stability (Yi et al 2007, Chasmer et al 2011, Fisher et al 2016 and a direct feedback mechanism is expected to control the temporal ecosystem evolution (Bonan et al 1992, Carpino et al 2018, Loranty et al 2018. This feedback mechanism (figure A1) is, however, poorly understood and broad-scale vulnerability studies do not yet exist. The canopy exerts shading by reflecting and absorbing most downward solar radiation and by suppressing the majority of turbulent heat fluxes in the belowcanopy space (Chang et al 2015). Further, the canopy controls the surface albedo, which is much lower than in grasslands especially during snow-covered periods (Bonan and Shugart 1989). The canopy decreases soil moisture and leads to a reduced thermal conductivity through precipitation interception (Thomas and Rowntree 1992) and higher evapotranspiration (Vitt et al 2000). Additionally, the canopy slows snow melting in spring and reduces snow compaction because of the suppressed turbulent fluxes, which therefore leads to higher snowpacks under denser canopies (Stuenzi et al 2021). Finally, the vegetation cover promotes the accumulation of an organic surface layer Shugart 1989, Yi et al 2007) which further insulates the topsoil from the atmosphere. A change in the forest density modifies the withinand below-canopy energy and water fluxes (Chasmer et al 2011, Stuenzi et al 2021. The forest composition also has an impact on the ground surface energy and water balance. Most boreal forests are dominated by evergreen needleleaf taxa, but wide areas of the north-eastern Eurasian continent are dominated by deciduous needleleaf taxa. The needle-shedding of deciduous taxa impacts the within and below canopy fluxes (Tanaka et al 2008, Zhang et al 2011, Peng et al 2020, Stuenzi et al 2021, the litter and organic surface layers (Bonan and Shugart 1989) and the fire regime (Rogers et al 2015). Since evergreen and deciduous taxa can establish under similar climate conditions Schweingruber 2004, Kharuk et al 2009) the successful spread of evergreen taxa into currently larch dominated areas and vice-versa, mainly depends on the frequency of disturbance events, which have increased over the past decades (Shuman et  Detailed modeling studies are needed to incorporate the local, heterogeneous, and complex feedback mechanisms, caused by the vegetation type and its relationship with topsoil temperature, active layer thickness (ALT), and available plant water (Tchebakova et al 2009, Schuur and Mack 2018, Kropp et al 2021, Stuenzi et al 2021. It has been shown that the integration of ecosystem components such as permafrost is highly relevant for projections on biomass and vegetation cover (Ito et al 2020). Here, we fill this gap between vegetation cover model projections and the actual physical impact this vegetation cover change has on permafrost ground, and present a detailed coupled permafrost-multilayer canopy model, developed for use in permafrostunderlain boreal forest systems.
We analyze (a) the trends of the two ecosystem changes, boreal forest densification and plant functional type composition, based on biome LPJ-GUESS model projections for north-eastern Siberia. Based on the projected trends in ecosystem changes, we use CryoGrid to simulate the projected ranges of (b) forest densities and (c) plant functional type compositions for three different study sites throughout north-eastern Siberia to investigate the impact of canopy variability on the ground thermal and hydrological regime. We thereby study the effect of the projected trends over an extensive range of predominantly deciduous-dominated boreal forests, as well as over different climate characteristics within the polar climate regime. This study delivers important insights into the range of spatial differences and possible temporal changes to the permafrost condition that can be expected following landscape changes such as deforestation through fires or other anthropogenic influences, afforestation in currently unforested grasslands, or the climate-induced densification of forested areas.

Study region
The treeline of north-eastern Siberia is dominated by the deciduous needleleaf tree genus Larix Mill. (figure 1), even though in mixed forest stands, larch taxa are out-competed by evergreen taxa, which is thought to represent the late-successional stage (Kharuk et al 2007). Once established, larch forests are likely to stabilize through a complex vegetationpermafrost-climate feedback system. Mainly shallow active layer depths hinder the establishment of evergreen taxa (Herzschuh 2019) and in more southern regions of eastern Siberia, larch is mixed with evergreen conifers (pine, spruce, fir) and hardwoods (Kharuk et al 2019). The ground vegetation is generally dominated by mosses and lichens that form carpets. Larch has shallow roots and preferably grows on clay permafrost soils with a shallow ALT and   Kotlyakov and Khromova (2002)), and Copernicus Global Land Service, Leaf Area Index (LAI), (after Copernicus Global Land Operations (2021)). Adapted from Stuenzi et al (2021). CC BY 4.0. maximum wetness of 20%-40%. Evergreen conifers and hardwood both prefer deeper active layers and a higher soil moisture availability (Ohta et al 2001, Rogers et al 2015. To capture these spatial differences across boreal forests, we study the current forest composition and structure along a east-west transect represented by three different sites as specified in table 1, figure 1, and appendix B with figure B1.

Projected forest evolution
We use ESA CCI Land Cover satellite data to parameterize forest composition and Copernicus Global Land Service leaf area index (LAI) data to parameterize forest density under current climate conditions (figure 1). To understand the current plant functional type distribution and the projected changes we study projections of LAI and plant functional types simulated by the LPJ-GUESS model as part of ISIMIP2b (Frieler et al 2017). We analyze the forest change scenarios from 2006 until 2099 under three global warming scenarios (RCP 2.6, RCP 6.0, RCP 8.5). The LPJ-GUESS dynamic vegetation model combines an individual-and patch-based representation of forest dynamics with biogeochemical cycling from regional to global scales and is further described in Smith et al (2014). The model is forced with bias-adjusted climate data from the Hadgem2-es earth system model. The EWEMBI dataset (Lange 2019) served as the basis for the trend-preserving bias adjustment of the GCMs at a daily time step (as detailed in Frieler et al (2017)). The data selected cover a region from E 105 • -167 • and N 45 • -70 • at a spatial resolution of 0.5 • × 0.5 • . We have separated this area into two individual study regions: West (E 105.25 • -137.25 • ) and East (E 137.75 • -169.75 • ) (figure 1), because of the differences in temperature and current vegetation cover between these two regions. We analyze projected LAI for needleleaf evergreen and needleleaf deciduous plant functional types under the three warming scenarios at the transect sites. Note that these LAI values are averaged annual values over the entire study sites, and do not represent the full summer LAI in deciduous forests. Therefore, we also study the projected monthly LAI (only available for the combination of all PFTs) for August, which corresponds to the maximum LAI of deciduous taxa plus the LAI of all other PFTs, including needleleaf evergreen and hardwoods, under the available warming scenarios (RCP 6.0, RCP 8.5) for the period 2006-2099 (figure B2).

Coupled permafrost-vegetation model
The model setup is based on the permafrost model CryoGrid (originally described in Westermann et al (2016)), a one-dimensional, numerical land surface model that simulates the thermo-hydrological regime of permafrost ground by numerically solving the heat-conduction equation (Nitzbon et al 2019). The CryoGrid model was extended by a multilayer canopy module developed by Bonan et al (2014) for the use in permafrost regions (appendix C and Stuenzi et al (2021) for model details). Here, we add a parameterization for deciduous forest to simulate the leafless state of deciduous-dominated regions outside of the short vegetative period in summer. This is achieved by allowing for separate leaf area index controlled by static time windows defining leaf-on and leaf-off season (10 October-10 April) following literature values for east Siberia (Spasskaya Pad) (Ohta et al 2001). Further, a more realistic canopy structure is simulated by allowing fractional composition of deciduous and evergreen taxa within the simulated forest stand. In addition, we test a parameterization for coupling forest density (LAI) to fine root biomass (R total , (gm −2 )) (appendix D). Further, we have implemented a new relationship for phase partitioning of water in frozen soil (freeze curve) based on Painter and Karra (2014) (appendix C).

Model simulations and setup
We ran model simulations for a wide range of forest types and forest compositions at the three transect sites. Parameters defining the canopy, snow, and soil properties were set according to literature values, model documentation, and own measurements (appendix for details). Tables E2 and E4 summarize the ground and vegetation parameter choices for all three sites. For each study site, we conduct 70 simulations representing different forest types and forest compositions (figure 2). The range of different forest types considered are bench-marked based on the projected ISIMIP2b data described by canopy density (leaf area index, LAI (m 2 m −2 )) between 0 and 7 m 2 m −2 and fractions between deciduous needleleaf and evergreen needleleaf taxa (0%-90% deciduous) (figure 2). To test the statistical significance of the differences between the simulated mean GSTs for varying forest densities and compositions, we apply variance analysis (one-way ANOVA) with a significance level of 0.001. Data are controlled for normal distribution and homogeneous variance across all groups. Statistical analyzes were performed using R software (R Core Team 2016).

Forest evolution under climatic warming
The biome projection data from LPJ-GUESS model simulations data reveal that in the eastern subdomain of our study region an increase in evergreen taxa are projected for all warming scenarios (RCP 2.6, RCP 6.0, RCP 8.5), with a peak in the yearly mean value of 0.5 m 2 m −2 and a maximum value of 2 m 2 m −2 , and 1.8 m 2 m −2 around 2075 respectively, for the RCP 6.0 and RCP 8.5 scenarios, followed by a decrease towards the end of the century (figure 3). In the western sub-domain, the three global warming scenarios project an increase in deciduous taxa. The overall LAI for August under the RCP 8.5 warming scenario increases by 2 m 2 m −2 in the western region Left: photographs from the three study sites and the respective model set-ups for Nyurba (NYU), mixed forest with a LAI of 3 m 2 m −2 (mid-density), for Spasskaya (SPA), deciduous-dominated forest with a LAI of 3 m 2 m −2 (mid-density), and for Chukotka (CHK), deciduous-dominated with a LAI of 1 m 2 m −2 (low-density). Right: schematic illustration of the range of possible forest cover scenarios (Forest density: low-density to high density, and forest composition: evergreen, mixed (10%-50% of deciduous taxa) or deciduous dominated (60%-90% of deciduous taxa), and no forest cover) caused by either, climatic changes and/or disturbance events such as i.e. an extreme drought, a fire event, logging, or pest infestation. Each forest cover scenario is simulated at each of the three sites, NYU, SPA and CHK.  and doubles in the eastern region (appendix B). Currently, we find mean values around 3 m 2 m −2 for the western study region and 1 m 2 m −2 in the east. According to the annual data (figure 3), the mean LAI is currently dominated by evergreen taxa. In the western region, this dominance switches around 2050 under all three climate forcing scenarios. Under the strongest climatic warming scenario, deciduous LAI increases to a mean value of 2.4 m 2 m −2 , a value three times higher than the end of the century deciduous taxa projection under RCP 2.6. The projection data reveals an increase in needleleaf evergreen taxa at both sites for the coming decade, followed by a decrease in the western region under all climate scenarios. In the eastern region, the increase continues until 2060, where-after the LAI of needleleaf evergreen taxa stays constant under RCP 2.6 and decreases under both the RCP 6.0 and 8.5 scenarios. In the western region, deciduous taxa will continue increasing until the end of the century under all climate forcing scenarios. Based on these data, which are in agreement with other model projections for Eurasia (Shuman et al 2014, Meredith et al 2019, we can constrain the expected changes in plant functional type compositions and forest densities for the entire eastern Siberian permafrost underlain boreal forest region east of 105 • and north of 45 • . The overall forest density is projected to increase with warming temperatures under all warming scenarios and for both study regions.

Permafrost sensitivity under changing forest density
The simulations clearly demonstrate that higher forest density leads to lower mean GST in the snowfree period. This trend is highly significant with p < 0.01 for Chukotka and Nyurba. The average snowfree GST is 1 • C colder for the simulations with the densest canopy covers. For Spasskaya Pad (SPA), this trend is reversed, showing an increasing GST for denser canopies (p < 0.01). In the snow-covered period mean GST increase with larger LAI values at Chukotka and Nyurba (p < 0.01) (figure 4). Temperature values for the simulations without forest are higher at all sites and for both time periods except for the snow-covered period at the Nyurba site, where the simulation without forest cover is 1.3 • C colder. The maximum difference between a sparse forest cover and no forest cover is a temperature increase of 8.3 • C in the snow-covered period at Chukotka.
Our model simulations show that the projected forest development alone exerts a strong control on the thermal state of the permafrost, in addition to the expected effect of a warmer and dryer climate itself. At two study sites, higher forest density leads to a significant decrease in ground surface temperatures in the snow-free period, while leading to an increase at the warmest site, SPA.
The magnitude of the insulation effect on the annual GST change from no forest cover to a dense forest cover is −6.3 • C at Chukotka, −0.2 • C at Nyurba, and −2.5 • C at Spasskaya.
The impact of forest density on GST consequently alters the annual ALT dynamics. We find a decline in maximum ALT with increasing canopy density for two sites in our study region. Highest maximum ALT of 0.68 m is found at the Spasskaya site with a LAI of 0 m 2 m −2 . The lowest maximum ALT is simulated at the Chukotka site with a value of 0.2 m only for LAI 7 m 2 m −2 (figure 5). We find a significant trend (p < 0.01) of a decrease in ALT with an increasing canopy density in Chukotka but an insignificant trend in Nyurba. At the SPA site, our model predicts an increasing maximum ALT with an increasing canopy density from LAI 1 − 4 m 2 m −2 . The maximum ALT for the simulations without a forest cover is higher at all sites. The difference between LAI 1 m 2 m −2 and no forest cover is up to 0.33 m in Chukotka. The decrease in snow-free period insulation with higher forest density is strongest at the coldest site of Chukotka.
Here, the average maximum ALT of all simulations at highest forest density (7 m 2 m −2 ) is 0.22 m, while average ALT is 0.27 m for a low LAI (1 m 2 m −2 ). The maximum ALT under a dense forest canopy is thus found 0.05 m (−19%) lower than under a sparse canopy. At Nyurba we find an average maximum ALT value of 0.45 m for a sparse canopy as well as for a dense canopy. At SPA low-density forest results in a maximum ALT of 0.54 m, which is considerably lower than the mean value of 0.57 m (−5%) for highdensity forest.
In order to analyze the impact of forest density on soil hydrology, we investigate the total yearly available plant water within the active layer. We find a clear and significant trend at Chukotka and Nyurba, with a decrease in available plant water for higher forest densities. The Chukotka site reveals the available plant water to be three times higher for the simulation without forest cover. The soil moisture in the active layer steadily decreases with increasing forest density at Chukotka and Nyurba, whereas it remains constantly low for SPA. SPA is the driest site, both in terms of liquid and solid precipitation, which leads to a very low amount of available plant water together with a relatively shallow snow cover (< 0.2 m) during winter.
The available plant water found is up to four times higher for the non-forested simulation at the Chukotka site and up to two times higher at the Spasskaya site. This indicates that forest loss may trigger the development of wetter and potentially swampy soil conditions depending on precipitation, evaporation, and ALT. In contrast, forest cover loss leads to a reduction in available plant water (up to 50%) at Nyurba which is characterized by climate conditions similar to Spasskaya. These contrasting hydrological impacts were observed in the vicinity of the respective study sites of Spasskaya and Nyurba. The performed simulations, thus, reveal that boreal forest loss can amplify both the wetting and drying of sub-Arctic regions.

Permafrost sensitivity under changing forest composition
Across the three study sites, we find a significant trend (p < 0.01) in lower GST's with an increasing percentage of deciduous taxa in the snow-covered period (figure 6). A lower percentage of deciduous taxa leads to a significant increase in the mean wintertime GST at Chukotka and Nyurba. The forest enhanced insulation effect of evergreen canopies, compared to deciduous cover, reaches up to +2.7 • C during the snow-covered period at Chukotka. A cooling trend of lower percentages of deciduous taxa in the summer period is found at SPA and Chukotka (p < 0.01).
The magnitude of the insulation on the annual GST change from evergreen to deciduous forest cover is −2.3 • C at Chukotka, −0.3 • C at Nyurba, and −1.2 • C at Spasskaya.
Changes in deciduousness also affect maximum ALT and the available plant water (figure 7). At SPA, Chukotka and Nyurba we find statistically significant trends (p < 0.01) towards higher maximum ALTs with decreasing deciduous taxa. The difference in ALT between 90% and 0% deciduous taxa are +0.04 m (+15%) at Chukotka, +0.05 m (+11%) at Nyurba and +0.07 m (+11%) at SPA. We find statistically significant trends (p < 0.01) towards higher available plant water with decreasing deciduous taxa at Chukotka and Spasskaya.

Discussion
About 55% of the total global permafrost area is covered by boreal forest (Gruber 2012, Helbig et al 2016. The forest cover plays an important role in insulating and stabilizing the permafrost underneath. The magnitude of this is highly dependent on the forest density as well as on the forest composition and structure but this relationship has not yet been stud-  permafrost by covering a wide variety of forest densities and plant functional type compositions. We find forest density to significantly control the ground thermohydrological conditions at all sites, whereby trends strongly differ in magnitude and direction. The cooling trends of denser canopies at the wetter sites, and the warming trend at the driest site, are reflected in the ALT dynamics. At the coldest site, the maximum ALT under a dense canopy is 19% lower than under a sparse canopy. Forest loss leads to higher snow-free GSTs at all sites and higher snowcovered GST at Chukotka and Spasskaya, with a maximum temperature increase of +8.3 • C. In just five years the forest cover loss leads to a warming of the GSTs at the same order of magnitude as the projected temperature increase for boreal regions 4 • C-11 • C until 2100 (Meredith et al 2019). In the snow-covered period, a lower share of deciduous trees was found to lead to warmer GSTs at all three sites. This difference in insulation capacity between deciduous-and evergreen-dominated canopies is up to +2.7 • C at the Chukotka site and +1.5 • C at SPA. Deciduousness has a higher effect on the average GSTs in cold regions (Chukotka) and a significant effect on the snowcovered GSTs at all sites. We show that in addition to the previously described change in fire regime (Rogers et al 2015) and albedo decrease (Bonan and Shugart 1989), the lower insulation capacity of evergreen canopies will be an important factor in the spreading of evergreen taxa in eastern Siberia. The actual thermal and hydrological impact of the forest cover is therefore determined by the forest density and structure, highly dependent on the local climate and hydrological conditions, and therefore varies greatly between our study sites. We find that forest loss can amplify wetting as well as drying of the soil. The available plant water after forest cover loss is four times higher at the coldest site, two times higher at the warmest site, and 50% reduced at the driest site. Depending on precipitation and soil type, forest cover loss can induce both drying and wetting. Generally Very sandy soils explain the good draining conditions and the resulting drying trajectory at Nyurba, while the clay-containing soils at Spasskaya and Chukotka are drained less, and hence the forest cover change can lead to wetting.
In this study, we focus on the direct physical impact of forest change on the detailed thermal and hydrological conditions of permafrost ground underneath, rather than investigating the exact timing of these ecosystem changes because the simulations themselves are decoupled from projected climate forcing data. Because of the difference found in the forest cover's impact on the thermal regime of the permafrost ground, we argue, that specific, local and detailed land-surface models are needed to understand the complex dynamics in permafrost underlain boreal ecosystems. Further, higher detail in the simulated change to the thermal and hydrological conditions could be achieved by incorporating a change in the thickness or composition of the litter, moss, and organic layers over time, and by additionally simulating the plant functional type broadleaf, which can establish wherever sufficient precipitation is available (Kharuk et al 2009, Shuman et al 2011.
While knowledge about carbon sequestration through boreal forests is well-established, more and more studies have found that different processes can counteract the boreal forest's role as a carbon sink (Betts 2000, Bonan 2008). As such, a decreasing albedo due to afforestation has been found to lead to a positive climate forcing for certain regions (Bonan 2008, Stuenzi andSchaepman-Strub 2020). Further, forest loss can lead to reduced evapotranspiration and a resulting short-term positive forcing effect (Liu et al 2019), as well as to an increased surface albedo, mainly in the snow-covered-period, and hence, a strong cooling effect (Lyons et  We argue that the development of the forest cover does not only influence the future of the boreal forest's function as a carbon sink but also plays an important role in the stability of permafrost. We show that varying density and tree composition have a significant effect on the thermal and hydrological state of permafrost. The insulating effect of the forest cover depends on the local climatic conditions but significant impact was found at all sites. Finally, the structure and composition of forests are highly dependent on the local ecosystem resilience towards an increasing frequency and intensity of forest fires, rising air temperatures, and a decrease in precipitation (Shuman et al 2011, Mekonnen et al 2019. Especially, the favoring of different fire regimes between evergreen and deciduous taxa, as well as warmer and drier conditions, can lead to fast ecosystem shifts. Altered thermal conditions, soil drainage or higher soil wetness, enrichment in nutrients, and an increased active layer depth can all have a favoring effect on either evergreen needleleaf or deciduous hardwood expansion, lead to the complete loss of forest cover, or change the forest density (Takahashi 2006, Kharuk et al 2013, Liu et al 2017, Kropp et al 2021. Here, we show that these changes will cause a shift in the thermal and hydrological permafrost state, which potentially destabilizes tightly coupled ecosystem functions.

Conclusions
In this study, we can underlay the tightly coupled interplay between forest and permafrost development with a physically-based model and make predictions on the progression of ecosystem-protected permafrost under a variety of forest change scenarios. In summary, we identify the following key points: (a) A change in forest density clearly affects the ground surface temperatures at all sites. Temperature differences are highest at the coldest site and in the snow-free period. This is further reflected by a decrease in the maximum ALT of up to 0.05 m or 19% at the two colder sites. The direction of this trend highly depends on local climate conditions. (b) At all sites, simulations without a forest cover reveal higher maximum ALTs of up to 0.33 m and higher GSTs of more than 8 • C after only five years. The thermal impact of forest cover loss is on the same order of magnitude as the climate warming projected for the region until 2100. Complete forest loss is found to lead to a deepening of the ALT and a warming of GSTs at all sites, independent of local climatic conditions. (c) Depending on precipitation and soil type, forest cover loss can induce both drying and wetting. After forest cover loss, the available plant water is four times higher at the coldest site, two times higher at the warmest site, and 50% reduced at the driest site. (d) At all sites, deciduous dominated canopies reveal lower GSTs, especially during the snow-covered period. This difference in insulation capacity reaches up to +2.7 • C for pure evergreen stands and is likely an important factor controlling the spreading of evergreen taxa and controlling the resilience of ecosystem-protected permafrost.
In the light of increasing disturbances (such as fires and diseases) in boreal forests our conclusions have strong implications regarding permafrostvegetation-climate feedback mechanisms. Our simulations indicate a positive feedback between the successive establishment of evergreen taxa and active layer deepening which may accelerate further forest transformation and permafrost thaw. In addition, forest cover transformation will have a strong impact on the hydrological regime, which may further amplify climate-induced changes in near surface temperature and precipitation. In consequence, the feedback loop might be further amplified by increasing fire probability and disease vulnerability due to additional water stress. On the other hand, under wetter climate conditions, enhanced wetting can eventually lead to swamping and thermokarst causing forest dieback as observed in drunken forests.
SMS designed the study, developed and implemented the numerical model, carried out and analyzed the simulations, prepared the results figures and led the paper preparation. ML, SW, JB, UH and SK co-designed the study, and supported the interpretation of the results. SMS, ML, and SW implemented the code in the model and designed the model simulations. AG prepared the land cover projection data. SMS, UH and SK prepared and conducted the field work in 2018, SMS conducted the field work in 2019. SMS wrote the paper with contributions from all coauthors. UH, ML and JB secured funding.
No competing interests are present.

Data availability statement
The code is available on Zenodo  Figure A1. Interactions between the atmosphere, boreal forest cover, and permafrost. In red (energy) and blue (hydrology) are all the mechanisms changing due to climatic changes. Climatic change leads to a change in the forest density and structure, which leads to changes in all feedback processes between forest and permafrost. Additionally, climatic change leads to permafrost thawing and a change in the water availability, which also leads to forest density and structure changes.

Appendix B. Study sites and their climate B.1. Study sites B.1.1. Nyurba
The most western study site is located south east of Nyurba at N 63.08 • , E 117.99 • , and 117 m asl, in a continuous permafrost boreal forest zone intermixed with some grassland, agricultural usage, and shallow lakes. The soils are sandy, and nutrient-poor (Chapin et al 2011). The forest soil has a litter layer of 0.07 m and an A-horizon reaching a depth of 0.16 m. It is rich in organic and undecomposed material. Mineral soil is podzolized. The rooting depth is 0.20 m. The average ALT between spatially distributed point measurements was 0.75 m in mid-August 2018 and 0.73 m in early August 2019. The forest is rather dense and mixed, with evergreen spruce Picea obovata Ledeb. and deciduous larch Larix gmelinii Rupr. The average tree height is 8 m (6 m for spruce and 12 m for larch). This site has been used as the main validation site in Langer et al (2020), Stuenzi et al (2020Stuenzi et al ( , 2021.

B.1.2. Spasskaya Pad
The central study site is the well-described forest research site in SPA at N 62.14 • , E 129.37 • , at 237 masl (Ohta et al 2001, Maximov 2015. SPA is located in a continuous permafrost region, and the active-layer depth is around 1.2 m in larch-dominated forests. The soils are sandy loam, and nutrient-poor. The forest soil has a litter layer of 0.08 m and an A-horizon reaching a depth of 0.16 m, mineral soil is podzolized. The measured average tree height is 12 m. Understory vegetation (Vaccinium L.) is dense and 0.05 m high. This site has been used as an external validation site in Stuenzi et al (2021).

B.1.3. Chukotka
The most northern study area is located at Lake Ilirney in Chukotka at N 67.40 • , E 168.37 • , and 603 m asl. The treeline here is dominated by deciduous larch and underlain by continuous permafrost. The soil is clay dominated with a litter layer of undecomposed Betula roots, dead moss, and dense rooting (0.01 m). The organic horizon consist of organic black hummus with highly decomposed organic material, moss remains, and good rooting (−0.18 m). The thawed mineral sediment layer had a thickness of (0.37 m) with little roots, dark grey clay matrix (40%), and clasts (60%). The average measured tree height is 11 m. Figure B1. Monthly average temperature (red) and total monthly solid and liquid precipitation (blue) for the three study sites (based on ERA-Interim (ECMWF Reanalysis) data for the study site coordinates). The exchange of sensible heat, radiation, evaporation, and condensation at the ground surface are simulated with an surface energy balance scheme based on atmospheric stability functions. In addition, the model encompasses different options to simulate the evolution of the snow cover including the Crocus snowpack scheme (Vionnet et al 2012) as implemented by Zweigel et al (2021). 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)  and cloud cover (Stuenzi et al 2021). We implement an updated model for the phase partitioning among liquid water, water vapor and ice based on the paramterization in Painter and Karra (2014). The proposed relationship for phase partitioning of water in frozen soil shows an improved performance for unsaturated ground conditions by smoothing the thermodynamically derived relationship to eliminate jump discontinuity at freezing. The flow in freezing soil is solved by a modified nonisothermal Richards equation. This constitutive relationship is more applicable for soils with very low total water content, which is the case for some regions in south and eastern Siberia, or high gas content. Following experimental results in Painter and Karra (2014), the ratio of ice-liquid to liquidair surface tensions for noncolloidal soil, β = 2.2, and the smoothing parameter, ω = 0, with n and α following the parameterization in the van Genuchten-Mualem model (van Genuchten 1980). This leads to an improved model performance for very dry ground conditions at boreal study sites in eastern Siberia.

B.2. Monthly forest cover projection data
The subsurface stratigraphy extends to 100 m, where the geothermal heat flux is set to 0.05 W m −2 (Langer et al 2011b). The ground is divided into separate layers in the model, the top 8 m have a layer thickness of 0.05 m, followed by 0.1 m for the next 20 m, 0.5 m up to 50 m and 1 m thereafter. The remaining CryoGrid parameters were adopted from previous studies using CryoGrid (table E1) (Langer et al 2011a, Nitzbon et al 2019, Stuenzi et al 2021. The model runs are initialized with a typical temperature profile of 0 m depth: 0 • C, 2 m: 0 • C, 10 m: −9 • C, 100 m: 5 • C, 5000 m: 20 • C. The remaining CryoGrid parameters were adopted from previous studies using CryoGrid (Langer et al 2011a, 2011b, 2016, Nitzbon et al 2019, Stuenzi et al 2021. The subsurface stratigraphy is described by the mineral and organic content, natural porosity, field capacity and initial water/ice content. Some of these parameters could be measured at the forest sites and were used to set the initial soil profiles and current canopy cover (tree height, forest composition) in the model (AsiaFlux 2017, Langer et al 2020, Stuenzi et al 2020).

C.2. Fine root biomass
Here, we use root/shoot ratio (R RS ) of 0.32 as defined in Jackson et al (1996) to calculate the fine root biomass correspondent to each LAI value to evaluate the importance of constraining this parameter.
with the specific leaf area at the top of canopy SLA d = 0.008 m 2 g −1 C for deciduous and SLA e = 0.008 m 2 g −1 C for evergreen taxa, respectively. The carbon content of the dry biomass is F carbon = 0.5 gC g −1 and the ratio of the fresh biomass that is water F water = 0.7 g H 2 O g −1 (Bonan et al 2018). Simulation results for LAI = 1 − 7 m 2 m −2 and corresponding root biomass = 267-1867 g biomass m −2 /m 2 reveal no consistent trend across the three study sites ( figure C1). The change in ALT caused by an increasing root biomass is up to 0.05 m at Chukotka and SPA and therefore on the same order of magnitude as found for forest deciduousness. We argue that root biomass needs to be constrained further but this is outside of the scope of this study.

C.3. Canopy description
The canopy is described by the leaf area index, the stem area index, and the leaf density function. LAI describes the total leaf area, which can be measured by harvesting leaves and calculating the total mass to canopy diameter ratio or estimated from below-canopy light measurements. The most common form of LAI estimation is from satellite data and the variance in values is rather high. LAI is measured at the bottom of the canopy and defines the total one-sided leaf area or the total projected needle leaf area (m 2 m −2 ) of all leaves per unit ground area (Myneni et al 2002, Chen et al 2005. LAI can be estimated from satellite data, calculated from belowcanopy light measurements or by harvesting leaves and relating their mass to the the canopy diameter. To assess the LAI in our study region we use data from literature and satellite data. Following Kobayashi et al 2010) who conducted an extensive study using satellite data, the average LAI of the forest types in our study region vary between 1 m 2 m −2 and 7 m 2 m −2 . Stem area index is not varied here and set to 0.05 m 2 m −2 , following Bonan (2019) and Stuenzi et al (2021). The leaf area density function is also not varied here and 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. Further, the lower atmospheric boundary layer is simulated by 4 m of atmospheric layers.

Appendix D. Model validation and in-situ measurements
We compare modeled and measured annual, snowfree and snow-covered mean GST to understand the overall model performance regarding the thermal regime of the surface and the ground and the relative temperature differences between the model and measurements. 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. The model has previously been validated against GST, radiation, snow depth, conductive heat flux, precipitation and temperature measurements for Nyurba and Spasskaya (Stuenzi et al 2021). To validate the model for all study sites used here, the model is validated against GST measurements at all sites. The data sets used for Nyurba and Chukotka cover one complete annual cycle from 10 August 2018 to 10 August 2019 (iButton (DS1922L), Maxim Integrated, accuracy: 0.5 • C Figure D2. Weekly averaged modeled vs. measured GST data and the respective regression line with regression coefficient R for all three sites (measured and modeled data cover the time periods 1 September 2018-10 August 2019 for Nyurba (grey) and Chukotka (blue) and 1 September 2017 to 10 August 2018 for SPA (yellow)). These differences fall into the range of 1.5 • C − 2 • C that is commonly used for validation purposes (Langer et al 2013.
(−10 • C to 65 • C, ). For SPA we have soil temperature measurements acquired through the AsiaFlux Network and the most recent data available and covering one entire year was used (August 2017-August 2018) (Asia-Flux 2017). The average annual GST recorded at the warmest study site in Nyurba at a depth of 0.03 m is 2 • C with an average of −9.4 • C in the snow-covered period and 5.4 • C in the snow-free period. Chukotka is the coldest study site with average snow-covered GST of −11.3 • C and 2.8 • C for the snow-free period, 1.9 • C colder than Nyurba and 0.5 • C colder than SPA. Here, the average snow-covered GST is −8.6 • C (figure D1). The model can reproduce these GSTs at all sites with a slight cold bias for the snow-free periods in Nyurba (−1.4 • C) and Chukotka (−1 • C) and a warm bias for the snow-covered period in Nyurba (+1.5 • C). These differences fall into the range of 1.5 • C − 2 • C that is commonly used for validation purposes (Langer et al 2013.