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 setups are often static or are not able to capture important processes such as the vertical structure 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; Community Land Model) and is used to reproduce the energy transfer and thermal regime at a study site (63.18946^∘ N, 118.19596^∘ E) in mixed boreal forest in eastern Siberia. An extensive comparison between measured and modeled energy balance variables reveals a satisfactory model performance justifying its application to investigate the thermal regime; surface energy balance; and the vertical exchange of radiation, heat and water in this complex ecosystem. 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 magnitude in storage heat flux comparable to that simulated at the grassland site. Further, the end of season snow cover is 3 times greater at the forest site, and the onset of the snow-melting processes are delayed.

1 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 $-\mathrm{50}{}^{\circ }\mathrm{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 (IPCC, 2014; Scheffer et al., 2012). Moreover, during 2007–2016 continuous-zone permafrost temperatures increased by $\mathrm{0.39}(\pm \mathrm{0.15}){}^{\circ }\mathrm{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; 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 of the near-surface permafrost. Further, extensive ecosystem shifts such as a change in composition, density or the distribution of vegetation (Holtmeier and Broll, 2005; Pearson et al., 2013; Gauthier et al., 2015; Kruse et al., 2016; Ju and Masek, 2016) and resulting changes to the below- and within-canopy radiation fluxes (Chasmer et al., 2011) have already been reported. Changes to the vegetation – permafrost dynamics – can have a potentially high impact on the numerous feedback mechanisms between the two ecosystem components. Increased soil carbon release from thawing permafrost through the delivery of soil organic matter to the active carbon cycle (Schneider Von Deimling et al., 2012; Romanovsky et al., 2017) is modified by vegetation changes, which can compensate for carbon losses due to an increased CO₂ uptake (as observed at ice-rich permafrost sites in northwestern Canada and Alaska; Estop-Aragonés et al., 2018) or even further accelerate total carbon loss.

These vegetation–permafrost dynamics in eastern Siberia have been documented through exploratory and descriptive field studies showing a clear insulation effect of forests on soil temperatures (Chang et al., 2015). Further, the biogeophysical processes controlling the evolution of the ecosystem have been described by conceptual models (Beer et al., 2007; Zhang et al., 2011; Sato et al., 2016). Modeling schemes such as ORCHIDEE–CAN (Organising Carbon and Hydrology In Dynamic Ecosystems – CANopy; Chen et al., 2016), JULES (Joint UK Land Environment Simulator; Chadburn et al., 2015), Lund–Potsdam–Jena (LPJ DGVM – Dynamic Global Vegetation Model; Beer et al., 2007), NEST (Northern Ecosystem Soil Temperature; Zhang et al., 2003) or SiBCliM (Siberian BioClimatic Model; Tchebakova et al., 2009) have added a vegetation or canopy module, with defined exchange coefficients for the fluxes of mass and energy, to their soil modules. This is feasible for varying levels of complexity, and the models are capable of addressing a variety of different aspects such as forest establishment and mortality (Sato et al., 2016), unfrozen vs. frozen ground and fire disturbances (Zhang et al., 2011), or the evolution of the vegetation carbon density under diverse warming scenarios (Beer et al., 2007).

While all of these studies have significantly improved our understanding of essential mechanisms in boreal permafrost ecosystems, it is important to further understand how a dynamic evolution of the forest structure and canopy affects the thermal state and the snow regime of the ground, especially amid ongoing shifts in forest composition (Loranty et al., 2018). The existing model setups are often static or not able to capture important processes such as the vertical canopy structure or the leaf physiological properties which strongly control the energy transfer between the top-of-the-canopy atmosphere and the ground. To our knowledge, so far, none of the existing models is able to capture the important processes of the vertical canopy structure in combination with a physically based, highly advanced permafrost model. The novel, physically based model introduces 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. This allows us to quantify the surface energy balance dynamics below a complex forest canopy and its direct impact on the hydro-thermal regime of the permafrost ground below.

With a tailored version of a one-dimensional land surface model (CryoGrid – CG, Westermann et al., 2016) we perform and analyze numerical simulations and reproduce the energy transfer and surface energy balance in permafrost-underlain boreal 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 (Langer et al., 2016; Westermann et al., 2016; Nitzbon et al., 2019, 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 the atmosphere, forest canopy and ground. We take advantage of a detailed in situ data record from our primary study site as well as from a secondary, external study site. These data are 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 new canopy module on the surface energy balance of the underlying permafrost at a mixed boreal forest site in eastern Siberia.

2 Methods

2.1 Study area

Our primary study site is located southeast of Nyurba (63.18946^∘ N, 118.19596^∘ E) 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 being evergreen 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 for larch, respectively. These boreal forest environments experience 6 to 8 months of freezing temperatures reaching extremes of $-\mathrm{62}{}^{\circ }\mathrm{C}$ in winter and up to 35 ^∘C between May and September. The low annual average temperatures result in continuous permafrost and are therefore poorly drained, podzolized and nutrient-poor soils (Chapin et al., 2011). Annual precipitation showed an increasing trend from 1900 until 1990, mainly due to an increase in wintertime precipitation. Between 1995 and 2002, summertime precipitation decreased by 16.9 mm in August and 4.2 mm in July (see Table 1 in Hayasaka, 2011, for further details). The temperature trend from 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 72.08^∘ N, Larix sibirica Ledeb. from 60 to 90^∘ E, Larix gmelinii Rupr. between 90 and 120^∘ E, and Larix cajanderi Mayr. from 120 to 160^∘ E (see Fig. 1). Larch competes effectively with other tree taxa because of its deciduous leaf habit and dense bark. In more southern margins of eastern Siberia, such as our study area, larch is mixed with evergreen conifers (Siberian pine (Pinus sibirica and Pinus sylvestris), spruce (Picea obovata Ledeb.) and fir (Abies sibirica)) and hardwood (Betula pendula Roth, B. pubescens Ehrh. and 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). This study comprises a secondary site for complementary model validation which is situated 581 km east of the primary site (for further description see Appendix C).

Figure 1 (a) Overview image of the location of the automatic weather station (AWS) in the grassland in the lower-right corner and the location of the instrumented forest site in the upper-left corner (Brieger et al., 2019a). (b) Respective soil profiles with the depths of the organic-matter-dominated O horizon, the top layer of the mineral soil containing the decomposed organic layer (Ah horizon), and the subsoil mineral layer (B horizon) at the grassland site (left) and at the forest site (right). (c) Global evergreen needleleaf (dark green) and summergreen needleleaf (light green) boreal forest distribution and the boundary line between the discontinuous and continuous permafrost extent in brown. The primary study site is marked with a black cross; our secondary study site is marked with a grey circle. Data: © ESA CCI Landcover classes, aggregated to summergreen and evergreen needleleaf classes (after Herzschuh, 2019) and permafrost extent from Land Resources of Russia – Maps of Permafrost and Ground Ice (after Kotlyakov and Khromova, 2002).

Bonan et al. (2018)Westermann et al. (2016)Bonan (2019)Bonan et al. (2018)Westermann et al. (2016)Zweigel et al. (2020)

Table 1 Overview of the processes for which this study differs from the former CG parameterizations.

Figure 2 Schematic of the surface energy and water balance of the forest (a, c) and grassland (b, d) schemes for (a, b) the snow-free period and (c, d) the snow-covered period. The active layer (AL), the permafrost (PF), tree and grassland (dotted lines), and the snowpack (Snow) in the snow-covered period. In each of the four panels, incoming and outgoing longwave (L_in and L_out), incoming and outgoing solar (S_in and S_out), turbulent fluxes (Q_h and Q_e), the storage heat flux (Q_s), and precipitation (and interception) of rain- and snowfall are shown.

Figure 3 Surface energy balance for the snow-covered (28 October 2018–27 April 2019) and snow-free (10–27 October 2019 and 28 April–10 October 2019) periods at the ground surface of grassland and forest and at the top of the canopy of the forest (Forest TOC). Shown are the net radiation (Q_net) and sensible (Q_h), latent (Q_e) and storage heat flux (Q_s) for the model runs of the forest and grassland site as well as the measured values at the grassland site. The bars indicate mean values, while the whiskers show the corresponding standard deviations.

2.2 Meteorological and soil physical measurements

An automatic weather station (AWS, Campbell Scientific; for a detailed list of sensors, see Table A1) is installed at $\mathrm{110}\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ on a meadow next to the forest patch described above. The grassland is grazed by horses in summer, and deforestation occurred more then 50 years ago. The AWS records air temperature and relative humidity at two heights (1.1 and 2.5 m) above the ground, while wind speed and direction are measured at 3.2 m above ground. In addition, the station measures liquid precipitation, snow depth, and incoming and outgoing short- and longwave radiation and is equipped as a Bowen ratio station (B, see Appendix A). All meteorological variables were recorded with 10 min resolution and stored as 30 min averages. In order to install soil temperature and moisture sensors in the ground a soil pit was excavated in the immediate vicinity (2.5 m) of the AWS. The O horizon has a depth of 0.04 m; the A horizon at 0.1 m contains undecomposed roots, dead moss remains, dense rooting and some organic humus. The mineral soil is podzolized, sandy and dominated by quartz. The rooting depth is 0.18 m. Iron-rich bands were found at 0.4, 0.7 and 1.1 m. The active-layer thickness was 2.3 m in mid-August 2018 and in early August 2019. In this soil pit, soil temperature and moisture measurement profiles are installed from the top to the bottom of the active layer, consisting of eight temperature sensors (0.07, 0.26, 0.88, 1.33, 1.28, 1.58, 1.98 and 2.28 m) and four moisture probes (0.07, 0.26, 0.88 and 1.33 m). In addition the conductive ground heat flux in the topsoil layer is measured with a heat flux plate installed at 0.02 m depth. Further, we record the near-surface ground temperature with five stand-alone temperature loggers (iButtons, see Table A1) with a measurement interval of 3 h. These are installed in the upper 0.03 m of the organic soil at our forest site. The forest soil has a litter layer of 0.08 m and an organic-rich A horizon reaching a depth of 0.16 m. It is rich in organic and undecomposed material. Mineral soil is podzolized, and the rooting depth is 0.20 m. The average active-layer thickness between spatially distributed point measurements was 0.75 m in mid-August 2018 and 0.73 m in early August 2019. In a vegetation survey along a 150 m transect from the grassland into the forest, the tree height of every tree within a 2 m distance was estimated. Trees <2 m were measured with a measuring tape; trees >2 m were measured with a clinometer or visually estimated after repeated comparisons with clinometer measurements. Together, the instrumentation, with a variety of different loggers, records the spatial and temporal variances across the two sites which are representative for a large area of the mixed boreal forest domain in eastern Siberia (https://doi.org/10.1594/PANGAEA.914327, Langer et al., 2020; https://doi.org/10.1594/PANGAEA.919859, Stuenzi et al., 2020b). In addition, we use a secondary study site, located 581 km east of Nyurba, with an extensive measurement record of radiation and temperature data from below the forest canopy, allowing for a comprehensive model validation (see Appendix C).

2.3 Model description

2.3.1 Ground module

CryoGrid is a one-dimensional, numerical land surface model developed to simulate landscape processes related to permafrost such as surface subsidence, thermokarst and ice wedge degradation. The model version is originally described in Westermann et al. (2016) and has since been extended with different functionalities such as lake heat transfer (Langer et al., 2016), multi-tiling (Nitzbon et al., 2019, 2020) and an extensive snow scheme based on Crocus (Vionnet et al., 2012; Zweigel et al., 2020). The thermo-hydrological regime of the ground is simulated by numerically solving the one-dimensional heat equation with groundwater phase change, while groundwater flow is simulated with an explicit bucket scheme (Nitzbon et al., 2019). The exchange of sensible and latent heat, radiation, evaporation, and condensation at the ground surface are simulated with a surface energy balance scheme based on atmospheric stability functions. In addition, the model encompasses different options for simulating the evolution of the snow cover including the Crocus snowpack scheme. 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). The change of internal energy of the subsurface domain over time is controlled by fluxes across the upper and lower boundaries written as

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle \frac{\mathit{\delta }E}{\mathit{\delta }t}}=S_{\text{in}}-S_{\text{out}}+L_{\text{in}}-L_{\text{out}}-Q_{\mathrm{h}}-Q_{\mathrm{e}}-Q_{{\mathrm{h}}_{\text{precip}}},\end{array}$$

where the input to the uppermost grid cell is derived from the fluxes of shortwave radiation (S_in and S_out) and longwave (L_in and L_out) radiation at the same time regarding the latent (Q_e) and sensible (Q_h) heat added by precipitation (Q_precip) and storage heat flux (Q_s) between the atmosphere and the ground surface (Westermann et al., 2016).

For this study, we have adapted the multilayer canopy model developed by Bonan et al. (2014, 2018). The canopy model is 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.

2.3.2 Canopy module

The multilayer canopy model provides a comprehensive parameterization of fluxes from the ground, through the canopy up to the roughness sublayer. The implementation of a roughness sublayer allows for the representation of different forest canopy structures and their impact on the transfer of vertical heat and moisture. The concept is similar to the multilayer approach in ORCHIDEE-CAN (Chen et al., 2016; Ryder et al., 2016). In an iterative manner, photosynthesis, leaf water potential, stomatal conductance, leaf temperature and leaf fluxes are calculated. This improves model performance in terms of capturing the stomatal conductance and canopy physiology, nighttime friction velocity, and diurnal radiative temperature cycle and sensible heat flux (Bonan et al., 2014, 2018). The within-canopy wind profile is calculated using above- and within-canopy coupling with a roughness sublayer (RSL) parameterization (see Bonan et al., 2018, for further detail).

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 CLM plant functional type of evergreen needleleaf forest. Please note that all parameters defining the canopy are set as constant values so that vegetation is not dynamic and that changes in forest composition or actual tree growth are not considered in this study.

2.3.3 Snow module

The snow module employed in this study is based on Zweigel et al. (2020) (submitted, available upon request). The CryoGrid model is extended with 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 a snow scheme (Bonan et al., 2018). Following Vionnet et al. (2012) the microstructure of the snowpack 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 wind-speed-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 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 time step 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 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 snowpack itself (see Table A2). Variables exchanged based on the snow cover are the 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 Eq. (4).

2.3.4 Interactions between the modules

The vegetation module forms the upper boundary layer of the coupled vegetation–permafrost model and replaces the surface energy balance equation used for common CryoGrid representations. The top-of-the-canopy (TOC) surface energy balance is calculated by the vegetation module based on atmospheric forcing. The forest module numerically solves the energy balance of the ground surface below the canopy defined as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{S}_{{\text{in}}_{\text{canopy}}}-S_{{\text{out}}_{\text{ground}}}+L_{{\text{in}}_{\text{canopy}}}-L_{{\text{out}}_{\text{ground}}}-Q_{{\mathrm{h}}_{\text{ground}}}\\ \text{(2)}& {\displaystyle }& {\displaystyle }-Q_{{\mathrm{e}}_{\text{ground}}}-Q_{{\mathrm{s}}_{\text{ground}}}=\mathrm{0},\end{array}$$

where $L_{{\text{in}}_{\text{canopy}}}$ and $S_{{\text{in}}_{\text{canopy}}}$ are the incoming long- and shortwave radiation at the ground surface and the lower boundary fluxes of the multilayer canopy module, $S_{{\text{out}}_{\text{ground}}}$ is the outgoing shortwave flux from the ground, $L_{{\text{out}}_{\text{ground}}}$ is the outgoing longwave flux, $Q_{{\mathrm{h}}_{\text{ground}}}$ is the sensible heat flux, $Q_{{\mathrm{e}}_{\text{ground}}}$ is the latent heat flux, and $Q_{{\mathrm{s}}_{\text{ground}}}$ is the storage heat flux at the ground surface. The first six components of the subcanopy energy balance directly replace the respective components surface energy balance scheme of CryoGrid (see Eq. 1). $Q_{{\mathrm{s}}_{\text{ground}}}$ is calculated based on temperatures of the uppermost ground or snow layers that are passed from CryoGrid to the forest module. The storage heat flux is calculated as

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{Q}_{{\mathrm{s}}_{\text{ground}}}=k{\displaystyle \frac{T_{\mathrm{s}}-T_{\text{ground}}}{\mathrm{\Delta }z}},\end{array}$$

where k is the soil thermal conductivity, T_s is the soil surface temperature, T_ground is the actual ground temperature for the first layer below the surface and Δz is the layer thickness. Soil thermal conductivity is parameterized following Westermann et al. (2013, 2016) and is based on the parameterization in Cosenza et al. (2003). The thermal conductivity of the soil is calculated as the weighted power mean from the conductivities and volumetric fractions of the soil constituents water, ice, air, mineral and organic (Cosenza et al., 2003). 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 setup 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_{{\text{ground}}_{\mathrm{s}}}$) is described as follows

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle \frac{\mathit{\delta }{W}_{{\text{ground}}_{\mathrm{s}}}}{\mathit{\delta }t}}=f{P}_{\mathrm{R}}+D_{\mathrm{c}}-E_{\mathrm{c}}+D_{\mathrm{t}}-E_{\mathrm{t}}\end{array}$$

and consists of the direct throughfall (fP_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) as

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle \frac{\mathit{\delta }{W}_{\mathrm{c}}}{\mathit{\delta }t}}=(\mathrm{1}-f-f_{\mathrm{t}})P_{\mathrm{R}}-E_{\mathrm{c}}-D_{\mathrm{c}},\end{array}$$

where $\mathrm{1}-f-f_{\mathrm{t}}{P}_{\mathrm{R}}$ is the intercepted precipitation, and the retained trunk water (W_t) is represented as

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle \frac{\mathit{\delta }{W}_{\mathrm{t}}}{\mathit{\delta }t}}=f_{\mathrm{t}}{P}_{\mathrm{R}}-E_{\mathrm{t}}-D_{\mathrm{t}},\end{array}$$

respectively. Lateral water fluxes are neglected in this baseline, one-dimensional model setup.

2.4 Model setup and simulations

We ran 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, and initial water and 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 the surface, 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 uppermost 8 m section has a layer thickness of 0.05 m, followed by 0.1 m for the next 20, 0.5 m up to 50 and 1 m thereafter. The remaining CryoGrid parameters were adopted from previous studies using CryoGrid (see Table A2) (Langer et al., 2011a, b, 2016; Westermann et al., 2016; Nitzbon et al., 2019, 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 the 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.

For the canopy stratigraphy, we follow the parameterizations in Bonan et al. (2018) for the plant functional type evergreen 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 the 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 use with multilayer models (see 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 and calculated from below-canopy light measurements or by harvesting leaves and relating their mass to the canopy diameter. Ohta et al. (2001) have described the monitored deciduous–needleleaf forest site at the Spasskaya Pad research station (our secondary study site, see Appendix C), which has comparable climate conditions but is larch-dominated. The value of the tree plant area index (PAI), obtained from fish-eye imagery and confirmed by litter fall observations, varied between 3.71 m² m⁻² in the foliated season and 1.71 m² m⁻² in the leafless season. This value does not include the ground vegetation cover. Further, Chen et al. (2005) compared ground-based LAI measurements to MODIS values at an evergreen-dominated study area (57.3^∘ N, 91.6^∘ E) southwest of the region discussed here, around the city of Krasnoyarsk. The mixed forest consists of spruce, fir, pine and some occasional hardwood species (birch and aspen). They find LAI values between 2 and 7 m² m⁻². To assess the LAI we use data from the literature and experience from the repeated fieldwork at the described site. Following Kobayashi et al. (2010), who conducted an extensive study using satellite data, the average LAI for our forest type is set to 4 m² m⁻², and the stem area index (SAI) is set to 0.05 m² m⁻², resulting in a plant area index (PAI) of 4.05 m² m⁻² and nine vegetation layers for model simulations. The lower atmospheric boundary layer is simulated by 4 m of atmospheric layers.

We perform simulations over a 5-year period from August 2014 to August 2019. The model runs are initialized with a typical temperature profile of 0 m depth at 0 ^∘C, 2 m at 0 ^∘C, 10 m at $-\mathrm{9}{}^{\circ }\mathrm{C}$, 100 m at 5 ^∘C, 5000 m at 20 ^∘C. Spin-up period prior to the validation period is 4 years before we compare modeled and measured data. Test runs with a longer spin-up period of 10 years confirmed that only 4 years are sufficient when focusing on GST. The meteorological forcing data required by the model include air temperature, relative humidity, air pressure, wind speed, liquid and solid precipitation, incoming short- and longwave radiation, and cloud cover. ERA-Interim (ECMWF Reanalysis) data for the coordinate 63.18946^∘ N, 118.19596^∘ E were used to obtain forcing data for the total available period from 1979 to 2019 (Simmons et al., 2007).

3 Results

3.1 Model validation and in situ measurements

At our primary study site, the model is validated against ground surface temperature (GST) measurements of forested and non-forested study sites. The dataset used covers one complete annual cycle from 10 August 2018 to 10 August 2019. In addition, the model output is compared to radiation, snow depth, conductive heat flux, precipitation and temperature measurements of the AWS at the grassland site. The AWS was set up on 5 August 2018 and taken down on 26 August 2019. Data were recorded continuously, except for 40 d in late May and early June due to a power cut. The mean annual air temperature was $-\mathrm{7.3}{}^{\circ }\mathrm{C}$ with a maximum temperature of 33.1 ^∘C, a minimum of $-\mathrm{54.0}{}^{\circ }\mathrm{C}$ and an average relative humidity of 70.5 %. Precipitation is 129.8 mm yr⁻¹ (liquid). The maximum snow height at the grassland site is measured to be 0.5 m in February, and the ground was snow-covered for 181 d from 28 October 2018 to 27 April 2019 (values above 0.05 m snow height). A quality check of the radiation data revealed partly inconsistent incoming longwave radiation measurements for the time span of 1 November 2018–26 February 2019. During this period it is likely that the sensor is partially covered by snow, making it necessary to discard those measurements from the record. High-quality Q_net and L_in and L_out measurements, thus, only exist for the periods 28 to 30 October 2018 and 27 February 2018 to 27 April 2019. This data gap consequently also limits the period for which Bowen ratios are calculated and sensible heat fluxes (Q_h) and latent heat fluxes (Q_e) can be derived. The mean annual grassland albedo is 0.35 with an average of 0.30 during the snow-free and 0.48 during the snow-covered season. From December to February the albedo reaches its highest values with a mean of 0.7. Mean annual GST at 0.07 m depth is $-\mathrm{2.6}{}^{\circ }\mathrm{C}$ (range from 19.1 to $-\mathrm{24.9}{}^{\circ }\mathrm{C}$) with an average of $-\mathrm{11.4}{}^{\circ }\mathrm{C}$ in the snow-covered period and 8.0 ^∘C in the snow-free period. The average annual GST recorded in forested areas at a depth of 0.03 m is 1.9 ^∘C (range from 15.6 to $-\mathrm{23.4}{}^{\circ }\mathrm{C}$) with an average of $-\mathrm{9.3}{}^{\circ }\mathrm{C}$ in the snow-covered period and 5.6 ^∘C in the snow-free period.

We acknowledge that the target variable GST does not allow a detailed evaluation of the surface energy balance. Therefore, we further validate the model performance with additional measurements from an external study site (see Appendix C). Through the Arctic Data archive System (ADS) we have been provided with meteorological and radiation data from beneath and above the larch-dominated forest canopy at Spasskaya Pad for 2017–2018 (Maximov et al., 2019). These data are used for additional model validation and are added to the Appendix of our paper. Overall our analysis reveals a satisfactory agreement between modeled and measured components of the surface energy balance below the canopy. Thus, we argue that the performance of the model at the external study site justifies its application at the primary study site in Nyurba, where below-canopy fluxes were not acquired (see Appendix C).

3.1.1 Surface energy balance

In a first step we assess the surface energy balance by comparing the modeled net radiation (Q_net), sensible heat flux (Q_h), latent heat flux (Q_e) and storage heat flux (Q_s) at the forested site and the modeled and measured fluxes at the grassland site (see Fig. 3 and Appendix A). Turbulent fluxes at the forest ground are close to zero for both the snow-free and snow-covered periods. The TOC sensible heat flux is highest in the snow-free period (48.6 W m⁻²), resulting in the highest net radiation flux (83.5 W m⁻²). Forest ground net radiation flux is only a third (25.8 W m⁻²) of the TOC flux. The latent heat flux in the snow-free period is similar for the forest TOC (16.5 W m⁻²) and grassland sites (measured value of 22.1 W m⁻² and modeled value of 18.5 W m⁻²). During the snow-covered season, forest TOC and ground turbulent heat fluxes and net radiation are close to zero. Net radiation flux in the snow-covered period is smallest at the grassland site (modeled value of $-\mathrm{20.0}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ and measured value of 9.7 W m⁻²). The resulting storage flux is more than double at the forest ground (30.4 W m⁻²) for the snow-free period and slightly positive (0.5 W m⁻²) during the snow-covered period.

Figure 4 Modeled incoming solar and longwave radiation for the snow-covered (28 October 2018–27 April 2019) and snow-free (10–27 October 2019 and 28 April–10 October 2019) periods at the ground surface of forest (grey) and grassland (striped). Measured (black) incoming solar (for the same time periods) and longwave radiation (for 28–30 October 2018 and 27 February 2018–27 April 2019) are shown for the grassland site. The bars indicate mean values, while the whiskers show the corresponding standard deviations.

The average measured Bowen ratio (B) at the AWS is 1.04, with an average of 1.94 for the snow-free and 0.35 for the snow-covered periods. At the grassland site the model predicts a B of 1.8 for the snow-free period and −16.54 for the snow-covered period. This sums up to an annual average B of 1.09 for grassland. Modeled annual average B at the forest ground is more than double with 2.85 for the snow-free period and 2.01 for the snow-covered period. The top-of-the-canopy modeled annual average B is 2.99, with an average of −0.77 in the snow-covered period and 3.93 in the snow-free period.

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 44.2 W m⁻² higher at the forest ground. In the snow-covered period downward longwave radiation flux is 40.7 W m⁻² higher at the forest ground. This results in a surplus of energy of $+\mathrm{9.6}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ for the snow-covered and $+\mathrm{16.7}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ for the snow-free periods under the forest canopy compared to the grassland. The shortwave radiation reaching the forest ground is very small for both periods (9.9 W m⁻² for the snow-free and 2.8 W m⁻² for the snow-covered periods), showing that the canopy effectively intercepts (absorbs and reflects) most of the incoming shortwave radiation. Shortwave downward radiation at the grassland site is more than 19 times higher in the snow-free period and 18 times higher in the snow-covered period.

Figure 5 Average measured and modeled GST in the snow-covered period, average measured and modeled GST in the snow-free period, annual average measured and modeled GST, and the respective standard deviations in forest (a at 0.03 m depth) and grassland (b at 0.07 m depth) over a measurement period of 1 year (10 August 2018–10 August 2019). The bars indicate mean values, while the whiskers show the corresponding standard deviations. Unpaired t test between modeled forest and grassland GST shows a statistically significant temperature difference for the snow-covered and snow-free periods.

Figure 6 (a) Modeled (grey) and measured (black) average weekly GST in 0.03 m depth in forest and (b) in grassland in 0.07 m depth with standard deviations for modeled (light grey) and measured (dark grey) values. In addition, the measured duration of the snow-covered period is shaded at the grassland site. Please note that there is a data gap in the measured data at the grassland site in May–June (see Sect. 3.1 for further details).

3.1.2 Thermal regime of the ground near the surface

In a second step, we compare the annual average GST in the snow-free and snow-covered periods to understand the overall model performance and the relative temperature differences between the forest and grassland sites (see Fig. 5). We further discuss the annual cycle of the thermal development of the permafrost ground, the modeled and measured active-layer thickness, and the volumetric groundwater content at both of our study sites.

Measured and modeled average GST values are summarized in Fig. 5. The highest deviation between modeled and measured GST is found at the grassland site. Here, the model shows a cold bias of $-\mathrm{4.1}{}^{\circ }\mathrm{C}$ for the snow-covered period. Also, there is a cold bias of $-\mathrm{2.8}{}^{\circ }\mathrm{C}$ in the snow-free period in the forest. Overall, we find an average annual difference of 0.7 ^∘C between the two sites. This difference is 5.2 ^∘C for the snow-free season and 6.7 ^∘C for the snow-covered period, respectively.

For a more detailed understanding of the annual cycle of the thermal evolution of permafrost ground at our study sites, we compare the weekly averaged GST at the grassland and forest site (see Fig. 6). The more detailed analysis of the annual cycle reveals periods with distinct differences between the model simulations 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 show a much faster ground warming in spring. This difference is over 20 d at the forest site and 15 d 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 in measured values.

Figure 7 (a) Measured snow depth at AWS (black solid line) in the grassland, modeled snow depth in grassland (dashed light green) and modeled snow depth in forest (dashed dark green). (b) Active-layer thickness (ALT) dynamics, measured ALT at the grassland and forest sites (black, point measurements in 2018 and 2019 (grassland only)), modeled ALT in forest (dashed light green), and modeled ALT in grassland (dashed dark green) sites. Topsoil freezing starts at the beginning of October.

Figure 8 (a) Modeled longwave and solar radiation at the forest ground for an LAI of 1 m² m⁻² (Sparse, black) and LAI of 4 m² m⁻² (Nyurba, grey). (b) Modeled turbulent fluxes (Q_h and Q_e), net radiation (Q_net) and storage heat flux at the forest ground (Q_s).

Figure 9 Modeled (Wind min (grey) and Wind max (red)) and measured (black) average weekly GST in 0.03 m depth in forest with standard deviations for modeled (light grey and light red) and measured (dark grey) values. Red represents a simulation using the wind speed at the top of the canopy as the input value for snow compaction. Grey shows the standard model simulation using the wind speed at the canopy bottom for the snow compaction mechanisms.

To further investigate the temporal evolution of the permafrost ground, we compare modeled and measured active-layer thicknesses at both study sites. In the grassland, the modeled maximum active-layer thickness (ALT) is 2.35 m between 13 and 24 October 2018, with complete freezing occurring on 9 November and topsoil thawing starting on 3 May. The measured ALT in the grassland was 2.3 m in mid-August 2018 and early August 2019. The measured ALT in the forest was between 0.5 and 1.1 m in mid-August 2018. In the forest, the modeled maximum ALT is 2.05 m in October 2018 with freezing being completed on the 14 November. Topsoil thawing began on 23 June, 51 d later than in grassland. The modeled ALT in August 2018 was between 0.4 and 1.8 m and therefore overestimated by 0.3 m compared to the point measurements taken in August 2018 (see Fig. 7).

Moreover, the measured volumetric water content (VWC) in grassland reaches its maximum of 0.2 in August. The averaged measured VWC at the forest site in August 2018 was 0.3. The model can broadly reproduce this difference, but there is a model bias towards higher VWC for both sites. The modeled maximum water content in forest is 0.5 between the end of June and the beginning of July and is about 0.2 higher than in grassland, where the simulation shows a maximum VWC of 0.3 in August. The modeled winter ice content reaches a maximum value of 0.36 in grassland and 0.42 in forest.

4 Discussion

The model presented here is found to be capable of simulating the differences in the ground thermal regime between a forested and a non-forested site for permafrost underneath boreal forests. This can provide important insights into the range of spatial differences and possible temporal changes that can be expected following current and future landscape changes such as deforestation through fires, anthropogenic influences and afforestation in currently unforested grasslands or the densification of forested areas. The implemented scheme is able to simulate the physical processes that define the vertical exchange of radiation, heat, water and snow between permafrost and the 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 the annual cycle, which is 32 ^∘C higher than at the forest site. This vegetation dampening effect on soil temperature is well-described in the 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 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.7 ^∘C higher GST during the snow-covered period and 5.2 ^∘C lower GST during the snow-free period. Measurements reveal an average of 2 ^∘C higher GST in forest during the snow-covered period and 2.3 ^∘C lower GST during the snow-free period. Our model simulations show that the strong control on the thermal state of permafrost is a result of the combined effects of canopy shading, the suppression of turbulent heat fluxes, below-canopy longwave enhancement, increased soil moisture and distinct snow cover dynamics. These relevant processes controlling forest insulation will be discussed individually in the following subsections followed by a detailed discussion on model applicability and limitations.

4.1 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 the snow-free and snow-covered periods. The forest ground albedo therefore has little influence on the energy balance. This canopy shading effect makes the longwave radiation the largest source of radiative energy at the forest site for both time periods. A surplus of longwave radiation of 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 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 they are close to zero at the forest ground. In summary, 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 at the forest ground because of extremely low turbulent heat fluxes.

4.2 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 snowmelt, 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). Sugimoto et al. (2002) found that plants used rainwater during wet summers but meltwater 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 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 by holding the water table close to the surface, which improves plant water supply.

4.3 Insulating the litter and moss layer

The existence of a thick moss and organic layer on the forest ground can significantly lower ground temperatures due to the high insulation impact (Bonan and Shugart, 1989). The low bulk density and low thermal conductivity of the organic mat effectively insulate the mineral soil, causing lower soil temperatures and maintaining a high permafrost table. A thick moss–organic layer on the forest floor is an important structural component of the forest–permafrost relationship, controlling energy flow, nutrient cycling, water relations, and, through these, stand productivity and dynamics (Bonan and Shugart, 1989). At our forest site the moss coverage was found to range between 0 % and 40 % with thicknesses ranging from 0.005 to 0.02 m. This is a comparably thin moss layer, but it was taken into account in the ground setup for the forest site.

4.4 Snowpack dynamics

Snow cover dynamics on the other hand seem to be a highly important factor in regard to the thermal evolution of the underlying permafrost ground (Gouttevin et al., 2012). Snow cover is an essential ecosystem component, acting as a radiation shield, insulator and seasonal water reservoir. The forest canopy exerts a strong control on snow accumulation due to interception and reduced wind speeds (Price, 1988). To analyze characteristics in the snow cover evolution, we compare modeled and measured snow depths at the AWS (grassland site) with the modeled snow depth at the forest site. 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 snowpack in the forest accumulates to 1.2 m. The melting of snow starts around the same period but lasts longer up until the end of May. The snowpack at the forest site exhibits different characteristics than at the grassland site. The found differences are clearly induced by the canopy structure controlling snow interception within the canopy, mass unload from the branches to the ground, sublimation and snow compaction.

In addition, our simulations at the forest site show an increase in longwave radiation, a decrease in solar radiation and the oppression of turbulent fluxes, which additionally lead to the slower melting of snow, less snow compaction and therefore a higher snowpack. This is in agreement with earlier work on snowpack modeling in coniferous forests (Price, 1988). For example, Beer et al. (2007) note that vegetation effects such as solar radiation extinction and atmospheric turbulence have a far greater influence on snow cover dynamics in boreal forests of eastern Siberia than snow interception alone. In addition, Grippa et al. (2005) found that the leaf area index (LAI) and snow depth are highly connected.

To understand the high impact that coupling the vegetation has on the snow cover, we next study the surface energy balance simulated by our model for our specific forest study site and a hypothetical, sparse canopy with an LAI of 1 m² m⁻², while keeping all other parameters the same (see Fig. 8). In snow-free periods the mean 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). 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 longwave radiation remains the dominant energy component.

Snow depth analysis further reveals that for a sparse forest canopy the maximum snow depth reaches 0.23 m only, resulting in a maximum ALT of 0.88 m and an annual average δT of 0.8 ^∘C. This confirms that the thermal differences between forest and grassland sites are largely controlled by the impact of the canopy density on snow depth and density.

The snowpack at our primary site (mixed forest) reaches a maximum thickness of 1.2 m, which is in accordance with studies of boreal forest snow depths in other boreal regions such as in Canadian boreal regions, where, i.e., Kershaw and McCulloch (2007) found varying mean snowpack depths between 0.73 and 1.3 m in different forest types and depths of only 0.08 m in a tundra landscape. Further, Fortin et al. (2015) measured maximum snowpack heights between 0.7 and 0.9 m in a black-spruce-dominated forest–tundra ecotone. Similar values were also found in a study in mixed boreal forests in northeastern China (Chang et al., 2015) and in a more general large-scale approach for the circumpolar north (Zhang et al., 2018). This strong variability and heterogeneity in snow distribution has already been identified as a very important driver of the subsurface and hydrological regimes and runoff in unforested permafrost regions (Nitzbon et al., 2019).

Nevertheless, the strong delay between observed and modeled top ground thawing at the forest site (see Fig. 6) demands further investigation. The snow compaction currently used in the snow module is dependent on wind speed only (see Sect. 2.3.3). Due to the coupled forest canopy, modeled wind speed at ground level is strongly reduced to the minimum value of 0.1 m s⁻¹. Consequently, the snow compaction is remarkably low. To understand how much of the difference between modeled and measured spring GST can be explained by the underestimated snow compaction, we simulate an extreme case, using the above-canopy wind speed for the snow compaction processes (see Fig. 9). The use of the above-canopy wind speed and the resulting high snow compaction reduces the difference between modeled and measured GST in spring by about 50 %. This reduction arises from the lower insulation capacity of the thinner snowpack. However, the timing of topsoil thawing does not improve. This may be explained by the fact that the snow cover has, on the one hand, a lower depth but, on the other hand, a higher density, which results in the same snow water equivalent and an equally high amount of energy needed to melt the snow cover completely.

4.5 Applicability and model limitations

The presented model is largely able to reproduce recorded GSTs in forests. The detailed analysis of the annual cycle shows that the snowmelt period in spring is biased at the forest and grassland sites. In reality, the ground warms up faster than is modeled. In the forest this is most likely caused by a wrong representation of snowpack compaction and snowmelt. Our analysis reveals that an extreme case of snow compaction only partly reduces the difference between modeled and measured GST in spring. This points to more complex processes that control snowmelt in forest than are currently represented by the model. Thus, it would be highly desirable to obtain further field measurements in order to gain a better understanding of snowmelt processes in boreal forests.

An aspect not represented in the model is the moisture transport and migration in frozen ground including the forming of ice lenses and excess ground ice, which can have a high impact on the local micro-topography and the surface energy balance. Furthermore, lateral water flow and snow redistribution may be important processes to be investigated in the future, since they can strongly modify the ground thermal regime as well as the snowpack development.

Additionally, more detailed field studies and modeling exercises on the variation of canopy densities and structures should be carried out in order to obtain a better understanding of the impact of dynamic forest stand development on permafrost and vice versa. In combination, the above model limitations could explain a great part of the described GST and ALT differences between measurements and modeled simulations.

To simulate the needle tossing of deciduous larch, we have incorporated a leaf area index threshold between needle tossing and leaf-out (10 October–30 May) for simulations at our external validation site at Spasskaya Pad (see Appendix C). This tunes the model towards a more detailed representation of larch-dominated forests, which are particular to the secondary study site and large parts of eastern Siberia. The analysis reveals a satisfactory agreement between modeled and measured components of the surface energy balance below the predominantly deciduous forest canopy. The mixed forest cover at our primary study site only contains 7 % of deciduous larch trees. The LAI reduction implemented is therefore very small and had no noticeable effect on our results. This modification can be used to study 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 – Larix Vegetation Simulator, Kruse et al., 2016), 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.

5 Conclusions

This study presents a specific application of a novel, coupled multilayer forest–permafrost model which enables us to investigate the energy transfer and surface energy balance in permafrost-underlain boreal forest of eastern Siberia. By simulating interactions between the vegetation cover and permafrost, our modeling approach allows us to quantify and study the impact of the forest on the hydro-thermal regime of the permafrost ground below. An extensive comparison between measured and modeled energy balance variables (GST, Q_e, Q_h, Q_net, S_in and S_out) reveals a satisfactory model performance justifying its application to investigate the thermal regime and surface energy balance in this complex ecosystem. Despite overall good performance, the field measurements reveal model shortcomings during the snowmelt period. Based on this modeling exercise and field measurements, we investigate the thermal conditions of two landscape entities as they typically occur in the boreal zone. In regard to the forest insulation effect on permafrost and ongoing land cover transitions, this study delivers important insights into the range of spatial differences and possible temporal changes that can be expected following landscape changes such as deforestation through fires, anthropogenic influences, and afforestation in currently unforested grasslands or the densification of forested areas. The detailed vegetation model successfully calculates the canopy radiation and water budgets, leaf fluxes, and 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 over 30 ^∘C on the annual ground surface temperature amplitude of the permafrost. Further, forested permafrost maintains a higher soil water 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. In summary, we identify the following key points.

• i.

  The forest 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 forest canopy enhances the longwave radiation below the canopy by up to 20 %, similar to a greenhouse, which results in a comparable magnitude of storage heat flux for both the forest and the grassland sites.

• iv.

  Forested permafrost holds a higher groundwater content than the dry grassland site.

• v.

  Forest canopy shading leads to the slower melting of snow, less snow compaction and therefore a higher snowpack.

• vi.

  The differences in the thermal development of the forest and grassland sites are highly influenced by the depth, density and resulting insulation capacities of the snow cover, which are in turn controlled by the forest canopy density.

Appendix A: Bowen ratio and turbulent heat flux calculation

With the AWS equipped as a Bowen ratio station, B is calculated following Foken (2016) as

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\displaystyle }B={\displaystyle \frac{c_p}{L_{\mathrm{v}}}}\times {\displaystyle \frac{\mathrm{\Delta }T}{\mathrm{\Delta }q}},\end{array}$$

where the specific heat at constant pressure for moist air (c_p) is 1.006 $\mathrm{kJ}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ and the latent heat of vaporization of water (L_v) is 2260 kJ kg⁻¹; ΔT is the temperature difference between the two air temperature sensors at heights 0.115 and 0.252 m; and Δq is the difference in specific humidity calculated from measured relative humidity (ϕ), temperature and pressure.

Table A1 Sensors used for field measurements.

Kershaw and McCulloch (2007)Westermann et al. (2016)Boike et al. (2019)Westermann et al. (2016)Westermann et al. (2016)Langer et al. (2011a)Nitzbon et al. (2019)Boike et al. (2019)

Table A2 Overview of the CryoGrid parameters used.

Thereafter, latent heat flux (Q_e) is calculated as

$$\begin{array}{}\text{(A2)}& {\displaystyle }{\displaystyle }{Q}_{\mathrm{e}}={\displaystyle \frac{Q_{\mathrm{s}}-Q_{\mathrm{g}}}{\mathrm{1}+B}},\end{array}$$

and sensible heat flux (Q_h) is

$$\begin{array}{}\text{(A3)}& {\displaystyle }{\displaystyle }{Q}_{\mathrm{h}}=(Q_{\mathrm{s}}-Q_{\mathrm{g}}){\displaystyle \frac{B}{\mathrm{1}+B}},\end{array}$$

with the storage heat flux (Q_s) being

$$\begin{array}{}\text{(A4)}& {\displaystyle }{\displaystyle }{Q}_{\mathrm{s}}=L_{\text{in}}+S_{\text{in}}-L_{\text{out}}-S_{\text{out}},\end{array}$$

and Q_g being the convective ground heat flux.

Table A3 Ground setup for simulations. Depth in meters with all others in unitless volumetric fractions.

Table A4 Constants.

Bonan (2002)Bonan (2002)Bonan (2002)Bonan (2002)Bonan (2002)Bonan (2002)Bonan (2002)Bonan (2002)Bonan (2002)Bonan (2002)Bonan (2002)Bonan (2002)Bonan (2002)Bonan (2002)Bonan (2019)Bonan (2019)Bonan (2019)Bonan (2019)Bonan (2019)Bonan (2019)Bonan (2019)Oleson et al. (2013)Oleson et al. (2013)Bonan et al. (2018)Bonan (2019)Bonan (2019)Bonan (2019)

Table A5 Multilayer canopy parameters. PFT: plant functional type. VIS: visible. NIR: near-infrared. NET: needleleaf evergreen boreal.

Bonan (2019)Bonan (2019)Bonan (2019)Bonan (2019)

Table A6 Further ground parameters needed by the vegetation.

Appendix B: Direct and diffuse solar radiation components from cloud cover data

ERA-Interim cloud cover data (N) allow us to use a simple approach to differentiate the incoming shortwave radiation (S_in) into diffuse,

$$\begin{array}{}\text{(B1)}& {\displaystyle }{\displaystyle }{S}_{{\text{in}}_{\text{diffuse}}}=S_{\text{in}}\times \left(\mathrm{0.3}+\mathrm{0.7}\times (N/\mathrm{8}{)}^{\mathrm{2}}\right),\end{array}$$

and direct,

$$\begin{array}{}\text{(B2)}& {\displaystyle }{\displaystyle }{S}_{{\text{in}}_{\text{direct}}}=S_{\text{in}}-S_{{\text{in}}_{\text{diffuse}}},\end{array}$$

components, based on Younes and Muneer (2007).

Appendix C: External validation site “Spasskaya Pad”

Further validation of the model performance is performed for a well-studied research site at Spasskaya Pad at 62.14^∘ N, 129.37^∘ E. This additional validation site is located 581 km from our primary study site but allows for validating further model variables due to additional observational data. Spasskaya Pad is a continuous permafrost region, and the active-layer depth is about 1.2 m in larch-dominated forests. The main tree species is Dahurian larch (Larix gmelinii), and there is a stand density of 840 trees ha⁻¹. The understory vegetation (Vaccinium) is dense and 0.05 m high. In 1996 a 32 m observation tower was installed (Ohta et al., 2001) in larch-dominated forest. Through the Arctic Data archive System (ADS, http://ads.nipr.ac.jp/, last access: 3 September 2020) we have been provided with the most recent available meteorological and radiation data from beneath and above the larch-dominated forest canopy for the time period 2017–2018 (Maximov et al., 2019). Each variable used here is measured at exactly 5 min intervals, except radiation (1 min). Ventilated shelters cover air temperature and humidity sensors. Net all-wave radiation and the four components of radiation are measured every minute, and the data loggers record average, maximum and minimum values. Upward and downward longwave radiation is corrected using the sensed temperature at domes and sensor bodies. Ground temperature is measured at seven depths, and soil moisture is measured at five depths. A more detailed description of the sensors can be found in Table 1 in Ohta et al. (2001). We have set up and run a 6-year simulation for this study site (2013–2019), using ERA-Interim forcing data using the grid cell closest to the coordinate 62.14^∘ N, 129.37^∘ E. The measurement tower is situated in larch-dominated forest so that a simple leaf-off parameterization is implemented. Following Ohta et al. (2001) we define a partial leaf-off period from 10 October–30 May, resulting in a reduced winter LAI of 0.5 m² m⁻². The summer LAI is set to a constant value of 1.9 m² m⁻², and we use the measured average tree height of 18 m for setting up the canopy structure. In order to ensure consistent model validation with the primary study site we used identical soil parameters for the external study site. All soil parameters used are summarized in Table A3, while Table A4 summarizes the constants used. We make use of canopy parameters defined by the PFT deciduous needleleaf due to the dominance of deciduous larch. The subsurface (soil) stratigraphy extends to 100 m below the surface, where the geothermal heat flux is set to a standard value of 0.05 W m⁻² (Langer et al., 2011b). The ground is divided into separate layers in the model. The uppermost 8 m have a layer thickness of 0.05 m, followed by 0.1 m for the next 20, 0.5 m up to 50 and 1 m thereafter. All remaining model parameters were set to default values as defined in previous studies (see Table A2) (Langer et al., 2011a, b, 2016; Westermann et al., 2016; Nitzbon et al., 2019, 2020). Similar to the primary study site we use ground surface temperature (GST) as one of the target variables for model validation, measured and modeled at 0.2 m. In addition we use air temperature below the canopy, measured at the height of 1.2 m, net radiation (Q_net), latent (Q_e) and sensible (Q_h) heat flux, and incoming (S_in) and outgoing (S_out) shortwave radiation flux at the ground surface as additional target variables allowing for a comprehensive validation of the modeled heat and moisture exchange processes within and below the canopy.

We assess the surface energy balance by comparing the median weekly values of modeled and measured net radiation (Q_net), sensible heat flux (Q_h), latent heat flux (Q_e), and incoming (S_in) and outgoing (S_out) solar radiation at the forested site (see Fig. C1).

Modeled turbulent fluxes below the canopy are small during the snow-covered period, and measurement data are not available during this period. Modeled and measured sensible heat flux in the snow-free period differ by 0.1 W m⁻² only. Modeled latent heat flux is only a fourth of the measured value and therefore underestimated in our model. Modeled net radiation in the snow-free period (25.7 W m⁻²) is slightly above measured net radiation (19.4 W m⁻²). For the snow-covered period, median modeled net radiation is slightly below the measured median value. The incoming shortwave radiation measured and modeled for the forest site fit well with differences well below 10 W m⁻². The standard deviation of measured values is higher for all variables except snow-covered net radiation.

In a second step, we compare the modeled and measured annual median GST for the snow-free and snow-covered periods and air temperature below the canopy 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 (see Fig. C2).

The highest deviation between modeled and measured temperatures is found in the GST of the snow-free period. Here, the model shows a cold bias of $-\mathrm{2}{}^{\circ }\mathrm{C}$. For the snow-covered period the difference is 1.8 ^∘C. For the air temperature below the canopy the difference between modeled and measured in the snow-free period is 1.5 ^∘C; for the snow-covered period, the difference is again 1.8 ^∘C. This falls into the range of 1.5–2 ^∘C that is commonly used for validation purposes (Langer et al., 2013; Westermann et al., 2016).

Overall our analysis reveals a satisfactory agreement between modeled and measured components of the surface energy balance below the canopy. Thus, we argue that the performance of the model at the external study site justifies its application at the primary study site in Nyurba, where below-canopy fluxes were not acquired.

Figure C1 Modeled (grey) incoming and outgoing solar radiation (S_in and S_out) and turbulent fluxes (Q_h, Q_e and Q_net) for the snow-covered (28 October 2017–27 April 2018, above) and snow-free (10–27 October 2017 and 28 April–10 October 2018, below) periods at the ground surface of forest. The bars indicate median values, while the whiskers show the corresponding standard deviations.

Figure C2 (a) Modeled (grey) and measured (black) air temperature (^∘C) below the canopy and (b) ground surface temperature (^∘C) for both the snow-covered (28 October 2017–27 April 2018, above) and snow-free (10–27 October 2017 and 28 April–10 October 2018, below) periods at the forest site at Spasskaya Pad. The bars indicate median values, while the whiskers show the corresponding standard deviations.