Novel coupled permafrost–forest model (LAVESI–CryoGrid v1.0) revealing the interplay between permafrost, vegetation, and climate across eastern Siberia

Boreal forests of Siberia play a relevant role in the global carbon cycle. However, global warming threatens the existence of summergreen larch-dominated ecosystems, likely enabling a transition to evergreen tree taxa with deeper active layers. Complex permafrost–vegetation interactions make it uncertain whether these ecosystems could develop into a carbon source rather than continuing atmospheric carbon sequestration under global warming. Consequently, shedding light on the role of current and future active layer dynamics and the feedbacks with the apparent tree species is crucial to predict boreal forest transition dynamics and thus for aboveground forest biomass and carbon stock developments. Hence, we established a coupled model version amalgamating a one-dimensional permafrost multilayer forest land-surface model (CryoGrid) with LAVESI, an individual-based and spatially explicit forest model for larch species (Larix Mill.), extended for this study by including other relevant Siberian forest species and explicit terrain.

Following parameterization, we ran simulations with the coupled version to the near future to 2030 with a mild climate-warming scenario. We focus on three regions covering a gradient of summergreen forests in the east at Spasskaya Pad, mixed summergreen–evergreen forests close to Nyurba, and the warmest area at Lake Khamra in the southeast of Yakutia, Russia. Coupled simulations were run with the newly implemented boreal forest species and compared to runs allowing only one species at a time, as well as to simulations using just LAVESI. Results reveal that the coupled version corrects for overestimation of active layer thickness (ALT) and soil moisture, and large differences in established forests are simulated. We conclude that the coupled version can simulate the complex environment of eastern Siberia by reproducing vegetation patterns, making it an excellent tool to disentangle processes driving boreal forest dynamics.

1 Introduction

Boreal forests cover vast areas of the Northern Hemisphere with strong gradients in climatic conditions and environments. They established in the Northern Hemisphere after the Last Glacial Maximum, leaving only a thin stretch in the north bordering the Arctic Ocean of pristine tundra areas (Mamet et al., 2019; Bonan, 2008; MacDonald et al., 2008). Anthropogenic climate warming is leading to the relaxation of warmth deficit limits at the northern margins and hence invasion of the tundra at a rate that is still unclear (Berner et al., 2013; Rees et al., 2020). At the same time, large parts of boreal forests, especially in eastern Siberia, are exposed to increasing disturbances (such as fires and drought), potentially driving a forest transition from deciduous species to evergreen taxa (Bonan, 2008; Herzschuh, 2020; Herzschuh et al., 2016). Accordingly, wildfires lead to increased greenhouse gas emissions through burning biomass and by deepening of the seasonally thawed layer for decades. The forest transition furthermore reduces the albedo, leading to a net positive global warming feedback, which will likely not be offset by increased carbon sequestration of a denser understorey vegetation (Bonan, 2008). However, the involved forest dynamics and interactions with the atmosphere and soil need to be considered in sufficient detail to forecast more realistic projections and to better understand the consequences for the boreal permafrost ecosystems of Siberia (Kirpotin et al., 2021).

Forest modelling is typically done globally, including the carbon cycle and permafrost, but individual trees and all life history stages need to be considered for a precise simulation. Modern global models such as LPJ-GUESS (Zhang et al., 2013) include individual models. The global models are used to show that forests will change and advance north. However, migration lags are typically not represented, and only climate envelopes serve for the distribution of plant functional types (PFTs). Dispersal processes and complexities have recently been recognized (Snell, 2014; Snell and Cowling, 2015; Lehsten et al., 2019) but are not yet used as standard for simulations. Furthermore, most modelling schemes still start with established trees, which makes them more general and computationally effective for a global application but at the cost of losing important detail for ecosystem responses (see discussion in Kruse et al., 2016). Also, the use of representative grid cells on a large grid without considering landscape will cause deviations to an extent that is unclear as to whether the impact on results is large and significant or not. Individual-based models (IBMs) could help here as they have sufficient detail of represented species and ecosystems, but applications are therefore only possible on landscapes not continents (Grimm and Railsback, 2005; DeAngelis and Mooij, 2005). Nevertheless, IBMs are the best tools to understand a system and develop general responses that can then inform or guide global model development. Further, neither a radiative transfer scheme through a multilayer canopy nor detailed representation of permafrost is included in typical simulation approaches.

Here, we aim at creating a model system that can accurately assess detailed thermal and hydrological fluxes between permafrost and forest cover as recently developed by Stuenzi et al. (2021a, b). The new model will include a dynamic vegetation model, which has a full life cycle to allow intraspecies and interspecies interactions at all stages (seed–seedling–mature tree), leading to nonlinear behaviour of population dynamics as well as resolving a three-dimensional landscape that is available for Siberian treeline areas developed by Kruse et al. (2016, 2019a, b).

2 Methods

We further developed two models and start each description by using the overview part of the ODD protocol for describing individual-based models (Grimm et al., 2010). We describe the host model LAVESI in Sect. 2.1, a spatially explicit, individual-based model handling the full life cycle of tree species and interactions among individuals and their environment. The second model, CryoGrid, which is informed by the host model and delivers improved state variables back to LAVESI, is described in Sect. 2.2: a one-dimensional, numerical land-surface model that simulates the thermo-hydrological regime of permafrost ground by numerically solving the heat conduction equation.

2.1 The 2D vegetation model LAVESI

2.1.1 Model description

The Larix vegetation simulator (LAVESI) is an individual-based spatially explicit model that simulates larch stand dynamics (Kruse et al., 2016, 2018). Monthly temperatures of the coldest (January) and warmest (July) months as well as precipitation series can force this model. In addition, 6-hourly data on wind speed and direction are needed to simulate seed distribution and tree reproduction, growth, and death (Kruse et al., 2019a, 2018, 2016). Recently the model has been extended by including topography and landscape sensing of the individuals (Sect. 2.1.2), and further boreal forest tree species were introduced aside from the larch species the model was initially developed for (Sect. 2.1.3; see full changes in Appendix A).

Purpose

The novel coupled model LAVESI–CryoGrid v1.0 was set up to understand tree stand structure, migration, and population dynamics of boreal forests growing between the leading edge at the Siberian treeline ecotone and the southern limit in response to a changing climate and its feedbacks with permafrost soils.

Table 1 State variables of the model with dimensions and update rates ordered by the model's hierarchy levels. Based on the table provided in Kruse et al. (2016).

^a Asterisks mean updated once at establishment (trees), production (seeds), or initialization (abiotic environment). NDD₀ indicates the number of days exceeding 0 ^∘C, and AAT₁₀ is the sum of temperatures above 10 ^∘C.

Entities, state variables, and scales

The model consists of two hierarchical levels characterized by a set of variables (Table 1): (1) simulation areas characterized by the specific biotic and abiotic environment and (2) individual trees and seeds.

The individual simulation areas are variable and have a size of typically 510×510 m (for parameterization and simulation experiments) on which seeds and trees are exactly positioned by x and y coordinates. Using the basal diameter of individual trees, the plot is overlaid with a tree density grid with a resolution of 0.2×0.2 m.

Simulation runs proceed in yearly time steps. We performed simulations for years 1–2100, prolonged by Representative Concentration Pathway (RCP) prediction scenarios. Additionally, to reach stabilization of population dynamics and the forcing climate series, simulations were preceded by a stabilization period with a length of 1000 years (for parameterization and sensitivity analysis). All simulations start from bare ground introducing 5000 ha⁻¹ yr⁻¹ seeds in the first 50 years, and to allow for repopulation of simulation areas after extinction, 100 ha⁻¹ yr⁻¹ seeds are added every year to the simulation areas.

Process overview and scheduling

The simulation proceeds in yearly time steps from the beginning to the end of the input climate time series, which includes a stabilization period to ensure that emerging populations reach equilibrium with the environment. In each initialization phase of each simulation run, the weather data are processed and used to estimate maximum diameter growth (at basal and breast height) for each simulation year based on 10-year mean climate auxiliary variables (see details in Sect. 2.2.2 “Description of sub-models” in Kruse et al., 2016). Within the growth processes of the model, these variables are used to individually estimate the current diameter growth of trees constrained by their actual biotic (competition) and abiotic (landscape features: elevation, TWI, slope, soil moisture, active layer depth) environment (design concept: sensing). Stochasticity in the model was introduced by using random numbers generated with a pseudo-random number generator (mt19937_64, from the random library) to allow for different results between two or more consecutive runs of the model (design concept: stochasticity).

Within 1 simulation year, the following processes become consecutively invoked (see Fig. 2 in Kruse et al., 2016; detailed explanations for each process can be found in a corresponding section in 2.2.2 “Description of sub-models”). Update of environment. Interactions between neighbouring trees are local and indirect. Basal diameters of each individual tree are used to evaluate the competition strength. We use a yearly updated density map to pass information about competition for resources between trees (design concept: interaction). Further, a litter layer and the state variables of each grid cell are updated as well. Growth. The individual growth of basal diameter and, if a tree reached a height of 1.3 m, of breast height diameter is calculated from the maximum possible growth in the current year affected by the tree's density index and its abiotic environment. From the resulting diameters, the tree height is estimated differently for two height classes: smaller and greater than 1.3 m (design concept: collectives). Seed dispersal. Seeds in “cones” are dispersed from the parent trees at a set rate. The dispersal directions and distances are randomly determined from a ballistic flight influenced by wind speed and direction with decreasing probabilities for long distances and only to places lower than the release height. If dispersed seeds leave the extent of the simulated plot they are removed from the system, but optionally they could be introduced from the other side or only on the east–west margins, depending on the user's choice. Seed production. Trees produce seeds after the year at which they reached their stochastically estimated maturation height. The total amount depends on weather, competition, and tree size. Optionally, the pollen donor for the pollination of ovules of seeds produced can be selected by a wind-determined and distance-dependent probability distribution function using a von Mises distribution. Establishment. The seeds that lie on the ground germinate at a rate depending on current weather conditions and is constrained by the actual litter layer height. Mortality. Individual trees or seeds die, i.e. they become removed from the plot, at a specified mortality rate. For trees this is deduced from long-term mean weather values, a drought index, surrounding tree density, and tree age and size, plus a background mortality rate. Seeds, on the other hand, have the same constant mortality rate whether on trees and or the ground (design concept: emergence). Ageing. Finally, the age of seeds and trees increases once a year, and seeds are removed from the system when they reach a defined species age limit.

2.1.2 Addition of landscape sensing

Data from the digital elevation model (DEM) TanDEM-X 90 m were downloaded from the web service provided by the German Aerospace Center (DLR https://download.geoservice.dlr.de/TDM90/, last access: 31 July 2021; Krieger et al., 2007). Subsequently, the tiles were reprojected to the corresponding Universal Transverse Mercator (UTM) zone of the focus areas (Khamra N49, Nyurba N50, Spasskaya Pad N52). All tiles were merged for each subzone and resampled by linear interpolation from 90 to 30 m resolution using functions from the “raster” package in R (Hijmans, 2020). The results were imported in SAGA GIS version 2.3.2 (Conrad et al., 2015) and subjected to a basic terrain analysis tool using the standard parameters. The resulting rasters were water-masked using the cloud-based geospatial data analysis platform Google Earth Engine (GEE, Gorelick et al., 2017) to assess Sentinel-2 imagery between 1 May and 15 October 2018 with a cloud cover of less than 20 % and thresholds manually set for spectral band B12 (2190 nm) until all water was masked out by comparing them to true-colour satellite images. The DEM along with slope angle and terrain water index (TWI, moisture content) were cropped to 510×510 m (260 100 m²) areas for this study and exported as plain text files for import into LAVESI.

LAVESI reads these data provided in 30 m resolution and interpolates linearly from the closest four grid cells for each 20 cm grid tile of the environment grid. Based on empirical relationships of forest presence for combinations of slope angle and TWI established in a study by Shevtsova et al. (2021), an environment growth impact factor (Envirgrowth, 0–1) is calculated for each tile, and tree diameter growth at this position is reduced accordingly (Eq. 1, Appendix B):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\text{Envirgrowth}}_i& ={\displaystyle \frac{-\mathrm{0.045999}\cdot {\text{TWI}}_i+\mathrm{0.994066}}{\mathrm{2}}}\\ & +{\displaystyle \frac{\mathrm{0.85654}\cdot e^{-\mathrm{0.5}\cdot {\left({\text{Slope}}_i-\mathrm{8.78692}\right)}^{\mathrm{2}}/{\mathrm{6.90743}}^{\mathrm{2}}}}{\mathrm{2}}},\end{array}\end{array}$$

where TWI_i is the interpolated terrain water index of the 20×20 cm² environmental grid cell i, and Slope_i is the slope angle of the same grid cell i.

Seed dispersal has been improved. Seeds can now only be dispersed to places which are at the same elevation as or a lower elevation than the release height in the terrain.

Table 2 Species traits and model variable values in LAVESI either newly introduced for this study, adjusted from the initial version LAVESI v1.01 (*), or introduced in the predecessor LAVESI-WIND v1.0 (**). Values from * Abaimov et al. (1998), ** Sato et al. (2016), own analyses, educated guess, or parameter fitting.

2.1.3 Addition of species and estimating leaf area index (LAI)

Further species were added to the existing model presented in Kruse et al. (2016). To add a fast forward implementation of species in LAVESI, we modified the code so that the programme can be started with either one or all species in a mix simultaneously. The species are numbered (integer values), which are used internally to assess species-related variables (Table 2) when called for in the functions as necessary. Therefore, the code is independent from the species and allows adding species or functional types simply by adding a new line in the new specieslist.csv in the main folder of LAVESI.

For this study, we analysed field data from the Chukotka and central Yakutia 2018 expedition in the same way as we did for Chukotka (Biskaborn et al., 2019; Kruse et al., 2019b). In the area of central Yakutia, species belonging to the Pinaceae family form the forests. From these, two deciduous boreal forest tree species were sampled, Larix cajanderi Mayr. (LACA) and Larix gmelinii (Rupr.) Rupr., (LAGM), as were three evergreen species, Picea obovata Ledeb. (PIOB), Pinus sibirica Du Tour (PISI), and Pinus sylvestris L. (PISY) (Kruse et al., 2019b). While the two larch species are best adjusted to the harsh environment of northeastern Siberia and are able to grow on shallow active layers above permafrost, they differ mainly in their frost hardiness, and the species LACA can even endure colder temperatures in winter (Table 2). PIOB is a competitor for L. gmelinii growing at similar environmental conditions but preferring deeper thawed active layers of a minimum of 200 cm. On well-drained sites, PISY grows well and outcompetes the other species. In milder environments, LASI and PISI grow on similar sites as LAGM and PIOB.

Tree-ring width data were established from tree discs and cores collected from sites close to Lake Khamra and from the region Nyurba. The discs and tree cores were prepared by standard dendroecological processing steps: (1) sanding with progressively finer paper until tree rings are clearly visible, (2) making high-resolution images for a track with a binocular and attached camera, (3) detecting rings with CooRecorder (Cybis Elektronik & Data AB) and cross-dating, and (4) exporting individual tree-ring chronologies (more details in Kruse et al., 2020). Tree-ring width data per species were then imported to R using the dplR package (Bunn et al., 2020), and regression models were set up by fitting nonlinear functions using generalized least squares with the gnls function from the nlme package (Pinheiro et al., 2019). For each species, we extracted the median of the loess-smoothed (span: 1.5) yearly growth increase of individual trees and set up a generalized least squares regression using a nonlinear model. This was successful for LACA, LAGM, and PIOB, but not for PISI and PISY due to small sample sizes, for which current values of PIOB are used as a first estimate (Table 2). For each tree in the simulation, the maximum actual growth can be estimated with the following equation:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\text{TRW}}_{\mathrm{Species}j,\mathrm{year}t}& =\exp (\text{gdbasalconst}+\text{gdbasalfac}\cdot t\\ & +\text{gdbasalfacq}\cdot t^{\mathrm{2}}),\end{array}\end{array}$$

where TRW is the tree-ring width for species i at 1 year depending on the fitted parameters gdbasalconst, gdbasalfac, and gdbasalfacq.

Biomass data were prepared following the protocol of Shevtsova et al. (2020), and allometric relationships were established to empirically estimate the leaf area (LA) from total leaf biomass for each tree (Eq. 3), followed by a log–log linear regression forced to pass through the origin employing the basal diameter as an explanatory variable (Eq. 4). To estimate the LA for each tree, we used specific leaf area (SLA) parameters to translate from the dry weight of needles to leaf area (Eq. 3). For each species, the SLA was extracted from literature values: SLA_LAGM=120 cm² g⁻¹ (Xian-Kui et al., 2015), which was also used for the closely related sister species LACA, SLA_PIOB=50 cm² g⁻¹ (Konôpková et al., 2020), and SLA_PISY=50 cm² g⁻¹ (extracted from the most recent source Błasiak et al., 2021, although other values are reported; 34 cm² g⁻¹ in Reich et al., 1998, 40 cm² g⁻¹ in Mencuccini and Bonosi, 2001). For PISI no source for SLA values was found and we assume it is similar to PISY:

$$\begin{array}{}\text{(3)}& {\text{LA}}_{\mathrm{Tree}i,\mathrm{Species}j}={\text{BM}}_{\mathrm{dry}\mathrm{needle}}{\text{SLA}}_{\mathrm{Species}j}/\mathrm{100},\end{array}$$

where BM is the biomass of tree i in grams, and SLA is the specific leaf area for species j.

$$\begin{array}{}\text{(4)}& \log \left({\text{LA}}_{\mathrm{Tree}i,\mathrm{Species}j}\right)=a\cdot \log \left({\text{DB}}_{\mathrm{Tree}i}\right)+\mathrm{0}\end{array}$$

Here, a is the slope of the linear model fit and DB is the basal diameter of the tree i for species j.

During simulations runs with LAVESI, the LA for each individual tree is estimated based on the fitted linear regression model using the following equation.

$$\begin{array}{}\text{(5)}& \widehat{{\text{LA}}_{i,j}}=\exp (a\cdot \log \left({\text{DB}}_{\mathrm{Tree}i}\right))\end{array}$$

The leaf area index (LAI) of each CryoGrid 10×10 m grid cell in LAVESI is then the sum of leaf area values of present trees. When a tree crown area covers more than one cell, the value is distance-weighted on the closest grid cells. For each species, the crown radius is estimated from field data with a log–log linear regression, and the slope and y intercept are used in LAVESI: parameters crownradiusestslope and crownradiusestinterc, respectively (Table 2).

2.1.4 Addition of a dynamic litter layer and estimating active layer thickness (ALT)

A dynamic, growing litter layer with constant growth of 0.5 cm yr⁻¹ and stochastic disturbance effects was introduced in the Environmentupdate function of LAVESI. When the parameter litterlayer is switched on, each of the 20 cm grid cells have a chance that the litterlayerheight can be reduced. This is stochastically implemented and for each year there is a 10 % chance the litter layer is reduced by 10 %, a 9 % chance of a 25 % reduction, a 0.9 % chance of 50 %, a 0.09 % chance of 90 %, and a 0.01 % chance of a 99 % reduction. This leads to a litter layer of ∼15 cm in the areas of interest in simulation runs, as observed in the region of interest (Kruse et al., 2019b). With this functionality, locally acting insulation effects are included in the estimation of the actual ALT in one environment grid cell. The estimation of the maximum active layer thickness was already introduced in the original setup of LAVESI (Kruse et al., 2016) and still serves as a first estimate, thus making it possible to run a stand-alone simulation of LAVESI without coupling to the CryoGrid (see below).

2.1.5 Model validation

We compared results from simulations until the year 2015 with field inventories from the 2018 expedition and literature values, focusing on the following key regions: Lake Khamra (westernmost, warmest), Nyurba (intermediate, climate station), and Spasskaya Pad (easternmost for boreal forests of Yakutia) for which we used literature values for comparison. Values were in the range of expected results.

2.2 The 1D permafrost model CryoGrid

2.2.1 Model description

The model used to simulate the thermo-hydrological interactions between permafrost ground and the forest canopy is based on CryoGrid (originally described in Westermann et al., 2016). CryoGrid is a one-dimensional, numerical land-surface model that simulates the thermo-hydrological regime of permafrost ground by numerically solving the heat conduction equation. The CryoGrid model has recently been extended by a multilayer canopy module developed by Bonan et al. (2014) for use in boreal permafrost regions (see Stuenzi et al., 2021a, b, for model details). The multilayer canopy model provides a comprehensive parameterization of fluxes from the ground through the canopy layer up to a roughness sublayer. In combination with CryoGrid the canopy model replaces the standard surface energy balance scheme, while soil state variables are passed back to the forest module. Following Stuenzi et al. (2021b), a realistic canopy structure is simulated by allowing fractional compositions of deciduous and evergreen taxa within a simulated forest stand.

Purpose

The model CryoGrid-Vegetation (Stuenzi et al., 2021a) was set up to understand the heat and water exchange between the atmosphere, boreal forest, and permafrost. The coupled multilayer vegetation–permafrost model reproduces the energy transfer and thermal regime of typical boreal permafrost ecosystems at different study sites in boreal permafrost regions.

Entities, state variables, and scales

Model entities are multiple layers of atmosphere and vegetation (based on Bonan, 2008), as well as permafrost (based on CryoGrid, Westermann et al., 2016). The physically based, numerical land-surface model simulates the radiative heat and water transfer through the atmosphere, vegetation, and ground at a 1D scale.

The simulation proceeds at a 5 min time step. The numerical model simulates the aboveground and belowground temperature field based on temporally changing conditions at the ground surface and top of the canopy–atmosphere boundaries.

Process overview and scheduling

To simulate the thermo-hydrological regime of the permafrost ground CryoGrid solves the one-dimensional heat equation numerically including groundwater phase change. The canopy model was 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 model and provides a comprehensive parameterization of fluxes from the ground through the canopy up to the roughness sublayer. This allows the simulation of diverse forest canopy structures and their impact on the vertical moisture and energy transfer. The exchange of radiation, sensible heat, condensation, and evaporation at the different layers are simulated with a surface energy balance scheme based on atmospheric stability functions. In every time step top-of-the-canopy incoming radiation and precipitation are partitioned at each layer throughout the canopy. The change of internal energy in the subsurface domain (ground) over time is composed of fluxes across the upper (surface energy balance below the canopy) and lower (geothermal heat flux at 100 m depth, 0.05 W m⁻²) boundaries. The model simulates the evolution of the snow cover based on an extensive CROCUS-based snowpack scheme (Zweigel et al., 2021). Furthermore, rainfall and snowfall are intercepted throughout the canopy with only a fraction reaching the ground directly as throughfall. The remaining water and/or snow is added to the canopy layers as canopy water and/or snow, which either evaporates or reaches the ground as canopy drip or stem flow in the following time steps. The model is forced by standard meteorological variables, which can be obtained from automatic weather stations, reanalysis products, or climate models. The required forcing data include air temperature, precipitation (solid and liquid), wind speed, incoming shortwave and longwave radiation, humidity, and air pressure (Westermann et al., 2016).

Table 3 Study site parameters for CryoGrid.

2.2.2 Model validation and parameters

This entire model setup has previously been extensively validated for different study sites throughout our study region, including Nyurba (63.08^∘ N, 117.99^∘ E), Spasskaya Pad (62.14^∘ N, 129.37^∘ E), and Ilirney (67.40^∘ N, 168.37^∘ E) (Stuenzi et al., 2021a, b). Validation exercises were carried out based on measured and modelled ground surface temperature (GST), active layer thickness (ALT), soil moisture, Bowen ratio, and shortwave and longwave radiation below and above the canopy. Parameters defining the canopy, snow, and soil properties were set according to literature values, model documentation, and own measurements (see Stuenzi et al., 2021b, for constants and multilayer canopy parameter choices). Table 3 summarizes the parameter choices for the three different sites. Appendix C, Table C1 summarizes the commonly used CryoGrid parameters.

2.3 Coupling the models

The coupled model setup benefits from the detailed process implementation gained while developing the individual models and brings the 1D to a landscape simulation. Therefore, we can reproduce the energy transfer and thermal regime in permafrost ground as well as the radiation budget, nitrogen and photosynthetic profiles, canopy turbulence, and leaf fluxes, while at the same time predicting the expected establishment, die-off, and treeline movements of larch forests (Fig. 2). In our analyses, we focus on vegetation and permafrost dynamics and reveal the magnitudes of different feedback processes between permafrost, vegetation, and current and future climate in Siberia.

LAVESI serves as host model and can now be set to call individual CryoGrid instances in a given year. For this, the data in LAVESI are aggregated on a 10×10 m grid superimposed on the 20×20 cm grid. Key state variables are leaf area index (LAI), stem area index (SAI), fraction of deciduous species, litter layer height, organic layer height, albedo, and the soil moisture in percent (plant-available water, PAW), which are provided to CryoGrid. These values can either be sorted by LAI and exported for five quartiles (implemented but not used here) or from the three areas that are equal slices from left to right (used here, see Appendix B). When the output file is created, LAVESI can be set to either start CryoGrid directly via a system call or schedule the instance with a bash file for the workload manager slurm (Yoo et al., 2003). Based on the key state variables provided by LAVESI for each of the areas, CryoGrid starts three (or five) parallel simulations. Once the output has been written, LAVESI reads the file and produces for the three levels anomalies for available soil water and active layer thickness. With these anomalies, the 10×10 m CryoGrid grid in LAVESI is filled and, from this, the anomalies used to calculate the new values for each 20×20 cm environment grid cell. When the quartile mode is set, the state values are assigned to this grid calculated by linear interpolation of the LAI-sorted state values and anomalies are calculated as in the other mode.

The multilayer canopy model in CryoGrid requires a minimum LAI of 0.7 m² m⁻² and a minimum height of 1 m to successfully build the radiative transfer scheme from the atmosphere to the ground; therefore, forest cover below these values is ignored.

Figure 1 Scheme of focus areas covered on the Chukotka and central Yakutia 2018 expedition. Region II: boreal forest bi-stability summergreen–evergreen transition needed for the expansion of species in addition to Larix cajanderi and L. sibirica. Data in map: treeline (Walker et al., 2005); ESA CCI Land Cover forest classes (ESA, 2017).

Figure 2 Scheme of coupling CryoGrid and LAVESI as well as the involved processes.

2.4 Forcing data and landscape of focus areas

The meteorological forcing data required by the multilayer canopy-permafrost model (air temperature, air pressure, wind speed, relative humidity, solid and liquid precipitation, incoming longwave and shortwave radiation, and cloud cover) are obtained from ERA-5 (ECMWF Reanalysis, Hersbach et al., 2020). The data are extracted for the focus regions Nyurba (63.08^∘ N, 117.99^∘ E; covering sites EN18067, EN180-68, EN180-70), Spasskaya Pad (62.14^∘ N, 129.37^∘ E), and Lake Khamra (59.97^∘ N, 112.96^∘ E; covering EN18079–83; Fig. 1).

To provide a millennia-long time series for model spin-up of LAVESI these series were matched to historical climate data for the forcing retrieved from the 0.5^∘ × 0.5^∘ Climate Research Unit gridded Time Series (CRU TS version 3.23) monthly data (1901–2014) (Harris et al., 2020). By repeating the 20th century data in a loop, a 2100-year-long monthly climate series was established from 1 to 2100 CE for each focus region using the RCP2.6 prediction scenario.

2.5 Simulation experiments

We forced LAVESI simulations with the RCP2.6 climate scenario by calling CryoGrid first in 2015 and yearly in the following years, letting the simulation run until the year 2030. Simulation runs were started with the updated LAVESI version on an empty landscape with true topography starting at 1 CE to allow for spin-up and ending in 2100 CE. Into the empty landscape, seeds (5000 ha⁻¹ yr⁻¹) for initiating population establishment were introduced for the first 50 years. Subsequently, only 100 ha⁻¹ yr⁻¹ seeds were introduced to allow for re-establishment after a complete die-out of trees on the whole simulation area. Each simulation was rerun without calling CryoGrid to compare the differences when the improved active layer thickness and available soil water are used.

In addition to the simulation that uses equal proportions of seeds of each species introduced into the simulation area, we started individual simulations for each single species.

2.6 Statistical analyses

All statistical analyses in this study were performed in R 3.6.1 (R Core Team, 2019) and mostly using included standard functions, with the addition of functions from the package “lattice” (Sarkar, 2008) for plotting the data.

Figure 3 Thawing depth anomaly in the coupled simulation model for all focus areas (rows) and areas within the simulation areas (columns: I, II, and III) for year steps 2015–2030. LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris.

Figure 4 Soil moisture anomalies in the coupled simulation model for all focus areas (rows) and areas within the simulation areas (columns: I, II, and III) for year steps 2015–2030. LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris.

3 Results

3.1 Simulations results with LAVESI

A gradient of population densities (expressed in LAI values) forms, which negatively follows the TWI gradient on all sites (I: left, driest to III: right, wettest, Figs. 5–7, Appendix E, Fig. E1). Further, an increase in larch dominance towards nearly pure larch tree stands can be observed from Khamra (southwest, warmest) via Nyurba (intermediate) to Spasskaya Pad (northeast, coldest). Stand densities are highest in the mixed-species simulations at Nyurba (∼1.9 m² m⁻²), larger than at the warmers site Khamra (∼1.5 m² m⁻²), and lowest at Spasskaya (0.9 m² m⁻²). Regarding other species, PISY is present in mixed stands in small numbers and grows only in open stands in simulations with only this species, suggesting that this species prefers a certain environment (Appendix B and D). PIOB performs better in single monospecific runs, leading to larger LAI than in mixed stands, but also reaches dense populations in warm areas (Khamra) and has the smallest sizes at coldest sites (Spasskaya Pad, Figs. 5–7). LAGM grows under most conditions but not in the wettest areas (highest TWI values, Appendix D).

Figure 5 Leaf area index (LAI) values at Lake Khamra for the three areas within the simulation areas (I, II, III) on the same plot at which CryoGrid was called (upper row) and only LAVESI runs (lower row). LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris.

Figure 6 Leaf area index (LAI) values at Nyurba for three areas within the simulation areas (I, II, III) on the same plot at which CryoGrid was called (upper row) and only LAVESI runs (lower row). LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris.

Figure 7 Leaf area index (LAI) values at Spasskaya Pad for three areas within the simulation areas (I, II, III) on the same plot at which CryoGrid was called (upper row) and only LAVESI runs (lower row). LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris.

3.2 Comparing simulations with LAVESI and the coupled version

The values are very similar for the runs with all and individual species. In nearly all years, LAVESI's ALT values are higher by up to 20 cm (mean over all is 109.6±11.4 cm versus 96.1±10.2 cm, which is ∼14.1 %) at all focus regions (Fig. 3). The soil moisture anomaly fluctuates around 0 % at Lake Khamra, is lower in the coupled model for Nyurba by ∼10 %, and lower at Spasskaya Pad by ∼20 % than in simulations using only LAVESI (Fig. 4).

Smaller population sizes can be observed in all simulations, leading to a drop in LAI values when LAVESI is updated by CryoGrid (comparing lower to upper panels in Figs. 5–7). In CryoGrid coupled runs, species grow with less dense stands but still cover the same area. In the coupled runs, populations die out in some cases at the end of the simulation at year 2030 (Appendix D).

Table 4 Qualitative comparison of simulation results to expectations based on observations.

^a Kruse et al. (2019b); ^b Ohta et al. (2001), Sugimoto et al. (2002); ^c Mamet et al. (2019); ^d Liang et al. (2016).

4 Discussion

4.1 Simulation performance

Species preference matches observations and expectations (Table 4, Kuznetsova et al., 2010). Larches have a wide ecological niche and are widespread (Mamet et al., 2019). They are generalists and best adapted to the harsh Siberian environments that were predominantly wet but have now become drier with global warming (Churakova et al., 2021; Kharuk et al., 2021). Picea obovata grows best in the westernmost warm areas and reaches larger LAI and biomass than when growing in mixed stands competing with other species (Kharuk et al., 2007). This is as expected, and competition, as Wieczorek et al. (2017) show, seems to be a strong factor dampening the response of tree stands when climatic conditions improve. Further, the simulation of denser stands at the Khamra site contradicts the general observation that mixed-species stands are more productive and/or denser than the niches occupied (Liang et al., 2016), but depending on the stand structure it could be negative (Zeller et al., 2018) and is in line with the observation of Chen et al. (2003).

4.2 The coupled version LAVESI–CryoGrid

The simulations using site information for the three focus sites yield dense tree stands in LAVESI simulations but not in the coupled version. The coupled model results in smaller key soil parameters, active layer thickness (ALT), and plant-available water (PAW). PAW has a strong impact as trees grow poorly in conditions exceeding drought and waterlogging thresholds (e.g. Liang et al., 2014; Mamet et al., 2019; Lawrence et al., 2015; Barber et al., 2000). Drought leads to a higher mortality of trees, and in consequence, population simulations are driven close to extinction within the simulation duration of 15 years. However, when drought-adapted Pinus sylvestris occupies the niche, there is nearly no change and it could, in the end, colonize the simulation areas. This matches expectations and further implies that the model is reproducing the natural dynamics well (Table 4).

As fires become more intense and frequent under global warming, spruce or other species may become dominant rather than shade-intolerant larch species. In the currently naturally deciduous, larch-covered areas, evergreen taxa may invade and change the heat fluxes and energy balance with the threat of entering a positive feedback loop such as a deepening of the ALT (Bonan, 2008; Stocker et al., 2013; Stuenzi et al., 2021b).

In general, technical issues arise when coupling models and implementing I/O indirectly via output files. The two different time steps (years in LAVESI vs. 5 min in CryoGrid) and computational speeds lead to long computation times and high requirements of computational resources for the coupled version. To avoid any delay, a parallelization of CryoGrid simulations, as implemented here, is highly recommended, especially for dry study sites where a simulated year in CryoGrid can take up to 4 h. We find that simulations of homogeneous areas perform best, especially when the exchange is set up using three to five instances sorted by LAI. The constraint lies here and more instances would improve the representation of variants of deciduous- and/or evergreen-covered plots, but these improved LAVESI simulations come at the cost of computation time when not using the parallelized version as developed for this study.

5 Conclusions

The as-is application of LAVESI overestimates ALT values by around 14 %, and therefore we advise using the implemented correction from CryoGrid for forecasting forest dynamics in the proposed coupled version. The 3D simulations provide a way to understand permafrost distribution and interactions with vegetation. Further implementations, tracked online in the GitHub repository of LAVESI, include spatially explicit fire and trait variation and adaptation. This public sharing of the source code plus advances in both models allows easy exchange, development, and adaptation to further regions. This and the simple setup make the coupled model version easy to implement and thus offers wide applicability. However, fieldwork or literature values for present species and stand structure as well as soil moisture and temperature time series from remote areas are necessary to adjust parameters and adapt the model and species to local site conditions, which is an issue for the vast remote areas in Siberia. With increasing data from these remote areas, such as better satellite imagery coverage and resolution, the collection of more detailed field data (loggers recording soil temperatures and moisture in the active layer), and monitoring of permafrost dynamics and tree growth, the drivers of forest dynamics may be disentangled and thus improve the model.

Appendix A

Table A1 Complete ODD “overview” table for LAVESI. Improvements from left to right; normal text represents unchanged model parts, while bold font highlights novel developments.

Appendix B: Landscape defining the focus region's plot area

Figure B1 Elevation (DEM), slope angle, and terrain water index (TWI) define the environment growth impact (0 no growth possible; 1 good, no constraints) using an empirically fitted function for present forest growth at the area of interest Khamra.

Figure B2 Elevation (DEM), slope angle, and terrain water index (TWI) define the environment growth impact (0 no growth possible; 1 good, no constraints) using an empirically fitted function for present forest growth at the area of interest Nyurba.

Figure B3 Elevation (DEM), slope angle, and terrain water index (TWI) define the environment growth impact (0 no growth possible; 1 good, no constraints) using an empirically fitted function for present forest growth at the area of interest Spasskaya Pad.

Appendix C: CryoGrid model parameters and constants used

Table C1 Overview of the CryoGrid parameters used.

Appendix D: Spatial distribution of the leaf area index (LAI) for mixed-species and pure-species simulations at the focus regions in 5-year steps

Figure D1 Leaf area index (LAI) values of the CryoGrid grid aggregated at year 2015 at Lake Khamra. Upper row shows LAVESI–CryoGrid coupled; lower row shows LAVESI simulations. LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris.

Figure D2 Leaf area index (LAI) values of the CryoGrid grid aggregated at year 2020 at Lake Khamra. Upper row shows LAVESI–CryoGrid coupled; lower row shows LAVESI simulations. LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris.

Figure D3 Leaf area index (LAI) values of the CryoGrid grid aggregated at year 2025 at Lake Khamra. Upper row shows LAVESI–CryoGrid coupled; lower row shows LAVESI simulations. LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris.

Figure D4 Leaf area index (LAI) values of the CryoGrid grid aggregated at year 2030 at Lake Khamra. Upper row shows LAVESI–CryoGrid coupled; lower row shows LAVESI simulations. LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris.

Figure D5 Leaf area index (LAI) values of the CryoGrid grid aggregated at year 2015 at Nyurba. Upper row shows LAVESI–CryoGrid coupled; lower row shows LAVESI simulations. LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris. Simulations with PISI were not possible.

Figure D6 Leaf area index (LAI) values of the CryoGrid grid aggregated at year 2020 at Nyurba. Upper row shows LAVESI–CryoGrid coupled; lower row shows LAVESI simulations. LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris. Simulations with PISI were not possible.

Figure D7 Leaf area index (LAI) values of the CryoGrid grid aggregated at year 2025 at Nyurba. Upper row shows LAVESI–CryoGrid coupled; lower row shows LAVESI simulations. LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris. Simulations with PISI were not possible.

Figure D8 Leaf area index (LAI) values of the CryoGrid grid aggregated at year 2030 at Nyurba. Upper row shows LAVESI–CryoGrid coupled; lower row shows LAVESI simulations. LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris. Simulations with PISI and PISY coupled were not possible.

Figure D9 Leaf area index (LAI) values of the CryoGrid grid aggregated at year 2015 at Spasskaya Pad. Upper row shows LAVESI–CryoGrid coupled; lower row shows LAVESI simulations. LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris.

Figure D10 Leaf area index (LAI) values of the CryoGrid grid aggregated at year 2020 at Spasskaya Pad. Upper row shows LAVESI–CryoGrid coupled; lower row shows LAVESI simulations. LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris.

Figure D11 Leaf area index (LAI) values of the CryoGrid grid aggregated at year 2025 at Spasskaya Pad. Upper row shows LAVESI–CryoGrid coupled; lower row shows LAVESI simulations. LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris. Simulations with PISI coupled were not possible.

Figure D12 Leaf area index (LAI) values of the CryoGrid grid aggregated at year 2030 at Spasskaya Pad. Upper row shows LAVESI–CryoGrid coupled; lower row shows LAVESI simulations. LAGM: Larix gmelinii; PIOB: Picea obovata; PISI: Pinus sibirica; PISY: Pinus sylvestris. Simulations with PISI coupled were not possible.

Appendix E: Evolution of leaf area index (LAI) across the focus region simulation areas

Figure E1 Mean leaf area index aggregated east–west for each simulated focus area and time slice (2010–2030). Upper row shows LAVESI–CryoGrid coupled version; lower row shows only LAVESI.