Sensitivity of ecosystem-protected permafrost under changing boreal forest structures

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 vegetation-permafrost-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 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.

Zoom In Zoom Out Reset image size

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^\circ \times 0.5^\circ$. 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).

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_\mathrm{total}$, (gm⁻²)) (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).

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. Table E3 summarizes constants used. We perform model simulations over a time period of five years from August 2014 to August 2019. This equals a spin-up period of four years before comparing modeled and measured data. The meteorological forcing data (air temperature, air pressure, wind speed, relative humidity, solid and liquid precipitation, incoming long- and shortwave radiation, and cloud cover) are obtained from ERA-Interim (ECMWF Reanalysis) extracted for the three sites (N 63.08^∘, E 117.99^∘, N 62.14^∘, E 129.37^∘, and N 67.40^∘, E 168.37^∘) (Simmons et al 2007). We use ground surface temperature (GST, top 0.4 m of the soil column) as the major target variable for model validation (appendix D). We further analyze maximum yearly ALT and the available water for plants within this active layer (PAW).

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 ($\textrm{m}^{2}{\textrm{m}}^{-2}$)) between 0 and $7\,{\textrm{m}}^2{\textrm{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).

Zoom In Zoom Out Reset image size

The biome projection data from LPJ-GUESS model simulations data reveal that in the eastern sub-domain 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\,{\textrm{m}^2}{\textrm{m}}^{-2}$ and a maximum value of $2~\textrm{m}^{2}{\textrm{m}}^{-2}$, and $1.8~\textrm{m}^{2}{\textrm{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~\textrm{m}^{2}{\textrm{m}}^{-2}$ in the western region and doubles in the eastern region (appendix B). Currently, we find mean values around $3~\textrm{m}^{2}{\textrm{m}}^{-2}$ for the western study region and $1\,{\textrm{m}^2}{\textrm{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\,{\textrm{m}}^2{\textrm{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.

Zoom In Zoom Out Reset image size

The simulations clearly demonstrate that higher forest density leads to lower mean GST in the snow-free period. This trend is highly significant with $p~\lt$ 0.01 for Chukotka and Nyurba. The average snow-free 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 \lt$ 0.01). In the snow-covered period mean GST increase with larger LAI values at Chukotka and Nyurba ($p \lt$ 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.

Zoom In Zoom Out Reset image size

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\,{\textrm{m}^2}{\textrm{m}}^{-2}$. The lowest maximum ALT is simulated at the Chukotka site with a value of 0.2 m only for LAI $7\,{\textrm{m}^2}{\textrm{m}^{-2}}$ (figure 5). We find a significant trend ($p \lt$ 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~\textrm{m}^{2}{\textrm{m}}^{-2}$. The maximum ALT for the simulations without a forest cover is higher at all sites. The difference between LAI $1~\textrm{m}^{2}{\textrm{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~\textrm{m}^{2}{\textrm{m}}^{-2}$) is 0.22 m, while average ALT is 0.27 m for a low LAI ($1~\textrm{m}^{2}{\textrm{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 high-density forest.

Zoom In Zoom Out Reset image size

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 ($\lt0.2\,{\textrm{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.

Across the three study sites, we find a significant trend ($p \lt$ 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 \lt$ 0.01).

Zoom In Zoom Out Reset image size

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 \lt$ 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 \lt$ 0.01) towards higher available plant water with decreasing deciduous taxa at Chukotka and Spasskaya.

Zoom In Zoom Out Reset image size

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 permafrost-vegetation-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 co-authors. UH, ML and JB secured funding.

No competing interests are present.

B.1.1. Nyurba

B.1.2. Spasskaya Pad

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.

Zoom In Zoom Out Reset image size

The canopy model has been coupled to CryoGrid by replacing its standard surface energy balance scheme while soil state variables are passed back to the forest module. The vegetation module forms the upper boundary layer of the coupled model and replaces the surface energy balance equation used for common CryoGrid representations. The multilayer canopy model provides a comprehensive parameterization of fluxes from the ground, through the canopy up to the roughness sublayer, which allows the representation of different forest canopy structures and their impact on the vertical heat and moisture transfer.

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) (Westermann et al 2016) 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 liquid-air 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.

The subsurface stratigraphy extends to 100 m, where the geothermal heat flux is set to 0.05 W m⁻² (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, Westermann et al 2016, 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, Westermann et al 2016, Nitzbon et al 2019, 2020, 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).

Zoom In Zoom Out Reset image size

Here, we use root/shoot ratio ($R_\mathrm{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.

$$\begin{align} R_{\textrm{total}} = R_{\textrm{RS}} \times \mathrm{LAI} \times \frac{1}{\mathrm{SLA}_{\textrm{d/e}}} \times \frac{1}{F_{\textrm{carbon}}} \times \frac{1}{(1-F_{\textrm{water}})}, \end{align} \tag{ C.1 }$$

with the specific leaf area at the top of canopy SLA$_{d} = 0.008\,{\textrm{m}}^{2}\,{\textrm{g}}^{-1}{\textrm{C}}$ for deciduous and SLA$_{e} = 0.008\,{\textrm{m}}^{2}\,{\textrm{g}}^{-1}\,{\textrm{C}}$ for evergreen taxa, respectively. The carbon content of the dry biomass is $F_{\textrm{carbon}} = 0.5\,{\textrm{g}}{\textrm{C}}\,{\textrm{g}}^{-1}$ and the ratio of the fresh biomass that is water $F_{\textrm{water}} = 0.7\,{\textrm{g}}\,{\textrm{H}}_{2}{\textrm{O}}\,{\textrm{g}}^{-1}$ (Bonan et al 2018). Simulation results for LAI = $1-7~\textrm{m}^{2}{\textrm{m}}^{-2}$ and corresponding root biomass = 267–$1867~\textrm{g}~\textrm{biomass}~\textrm{m}^{-2}/\textrm{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.

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 ($\textrm{m}^{2}{\textrm{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 below-canopy 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~\textrm{m}^{2}{\textrm{m}}^{-2}$ and $7~\textrm{m}^{2}{\textrm{m}}^{-2}$. Stem area index is not varied here and set to $0.05~\textrm{m}^{2}{\textrm{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.