Simulating interactions between topography, permafrost, and vegetation in Siberian larch forest

In eastern Siberia, topography controls the abundance of the larch forest via both drought and flooding stresses. For the reconstruction of these topographical effects, we modified a dynamic vegetation model to represent soil water relocation owing to within-grid heterogeneity of elevation, over-wet-kill of trees, and air temperature differences within-grid. After calibration, the model reasonably reconstructed the geographical distributions of observation-based-estimates of fundamental properties of plant productivity and thermo-hydrology. Thus, the model appropriately responded to environmental gradients in eastern Siberia. The modified model also partially reconstructed the topography control on tree abundance and thermo-hydrology status in eastern Siberia, although its geographical distribution was not always good. In the modified model, soil water redistribution increased the risk of over-wet-kill in lower elevation classes, whereas it reduced the risk of over-wet-kill for larch trees in higher elevation classes. We demonstrated that without considering the latter effect, forest collapse due to over-wet stress would happen throughout eastern Siberia under a forecasted climatic condition during the 21st century, which will deliver a much moister environment throughout eastern Siberia. Therefore, modeling the over-wet-kill of trees without considering topographical heterogeneity would result in the overestimation of forest collapse caused by the over-wet-kill of trees under an expected climate trend in eastern Siberia.


Introduction
In eastern Siberia, larches (Larix spp.) often exist in pure stands, constructing the world's largest coniferous forest (Schulze et al 1995), to which changes can significantly affect the earth's albedo and global carbon balance (Bonan et al 1995, Gibbard et al 2005. Siberian larch forest is located in an extreme environment for the existence of a forest ecosystem; hence, small differences in environmental conditions will alter its presence. There is a quantitative pattern that topographic properties control the abundance of larch forest via both drought and flooding stresses in this region; larch abundance appears to be constrained by drought stress in mountainous regions, but flooding stresses in plain areas (Sato and Kobayashi 2018).
The larch-dominated region in Siberia is underlaid by permafrost and plant productivity is partially under the control of active layer thickness (ALT); deeper ALT increases the water holding capacity of soil, thereby enhancing plant productivity (Beer et al 2007, Sato et al 2016, but see also Zhang et al 2011. Soil wetness inversely influences ALT. For example, the greater ice content in the soil of moister areas facilitates upward heat conduction from the permafrost, which results in a shallower ALT (Woo andXia 1996, Carey andWoo 1998). Such a thermo-hydrological system is expected to change drastically because the arctic region has warmed more than twice as fast as the global average (Cohen et al 2014) and near-surface permafrost and seasonally frozen ground areas are projected to decline substantially during the 21st century (Lawrence and Slater 2005).
Therefore, for the quantitative reconstruction and projection of geographical distribution and functions of the Siberian larch forest, topography and thermo-hydrological dynamics mediated heterogeneity of the environment should be taken into consideration. Previous studies have focused on modeling the influences of thermo-hydrological dynamics on larch forest using process-based models (Beer et al 2007, Zhang et al 2011, Sato et al 2016. However, topography mediated heterogeneities of the environment within simulation grids were often neglected. Kleinen et al (2012) integrated a dynamic vegetation model, LPJ-DGVM with a wetland module, which dynamically determines the fraction of inundated areas and probability density function (PDF) of the water-table-depth within each simulation grid. Then, this integrated model simulates peat accumulation, which can only happen in water-saturated areas. This wetland module is based on TOPMODEL, which is a physically based streamflow and water-tabledepth computational scheme. TOPMODEL calculates hydrological heterogeneity within watersheds; it inputs water-table-depth averaged over the watershed and a PDF of the combined topographic index (CTI) and then outputs the mean watershed base flow and a PDF of water-table-depth (Beven and Kirkby 1979, Stieglitz et al 1997, Gedney and Cox 2003. In this study, a dynamic vegetation model, SEIB-DGVM (Sato et al 2016), which works with a thermohydrology land-physics-model NOAH-LSM 2.7.1 (Ek et al 2003), was extended for the reconstruction of the topography-mediated-heterogeneity of this environment within simulation grid cells. This extension was achieved by the implementation of the TOPMODEL. After the calibration of the model for eastern Siberia, we conducted a model experiment to quantitatively understand how plants, thermo-hydrology, and topography work together. Then, we adapted the model to a forecasted climatic condition at the end of the 21st century and evaluated how these modifications alter the prediction of the abundance of the Siberian larch forest.
Specifically, we addressed the following questions. (1) Do the model modifications contribute to the reconstruction of topography control on thermohydrology and tree abundance in eastern Siberia? (2) Is the mechanism that is responsible for the reconstruction consistent with observation trend? (3) How do the model modifications change the projection of tree abundance in eastern Siberia at the end of the 21st century?

The original model
All simulations in this study were conducted by the integrated model of the SEIB-DGVM (Sato et al 2016), which reconstructs plant and carbon dynamics under specified climatic conditions, and the onedimensional (i.e. not considering the lateral flow of water and heat) land-surface-model NOAH-LSM to mechanistically link thermo-hydrology with plant dynamics. In the model, vegetation fraction controls incoming solar radiation and surface heat flux conductance on top of the soil layer. Soil moisture and its state (liquid or solid) influence both hydraulic conductance and heat conductance, and thereby feedback to soil moisture and its state. This integrated model roughly reconstructs the latitudinal gradients for permafrost presence, soil moisture, canopy leaf area index (LAI), and biomass over the entire larchdominated area in eastern Siberia. Sato et al (2016) provides a full description of this integration. Here, we only outline the integratedmodel and describe the modifications made for this study. The simulation unit of the model is spatially explicit; a 30 m × 30 m 'virtual forest' in which individual trees establish, compete, and die. The ground is vertically divided into 20 soil layers, each 0.1 m deep. Plant dynamics in the virtual forest work in conjunction with the carbon and water cycles. In the integrated model, the SEIB-DGVM provides the NOAH-LSM with daily updates of LAI and vegetation fraction, whereas the NOAH-LSM provides the SEIB-DGVM with daily updates of the vertical profiles of soil wetness and temperature.
The Siberian larch forest generally suffers from severe nitrogen deficiency (Schulze et al 1995), which limits leaf carbon gain (Vygodskaya et al 1997). However, the model does not consider nitrogen limitation on plant growth; hence, our model might overestimate the CO 2 fertilization effect.

Model modifications
The simulation area is in the Sakha Republic of Russia, which corresponds to the geographical extent to which permafrost temperature and ALT data (Beer et al 2013) are available. Each simulation grid, at the resolution of 0.5-degree, was divided into 56 × 56 pixels (1/112-degree resolution), which corresponds to calibration and validation datasets. Land pixels of each simulation grid were sorted in ascending order of elevation and almost equally divided into 10 groups (hereafter, referred to as 'elevation rank'). Elevation rank of each simulation grid is represented by a single model simulation, assuming that all pixels at the same elevation ranked in a simulation grid share identical environmental conditions. As all elevation ranks represent nearly the same number of pixels, we regard stated variables averaged over the elevation rank as grid representing values.
Therefore, within-grid heterogeneity was only controlled by elevation in the modified model. This assumption was rationalized by (1) considering that elevation is the most straightforward way to implement within-grid heterogeneity of air temperature and (2) elevation tightly correlates with slope and soil depth profiles, all of which control hydrology, and tightly correlates with the flooding record in eastern Siberia (Sato and Kobayashi 2018).
TOPMODEL, which dynamically determines the heterogeneity of water-table-depth within a watershed, was employed in this study. Instead of detailed topographic properties, TOPMODEL requires a CTI, which is defined to be ln (α i /tan (β i )), where α i is area draining through location i per unit contour length, and β i is the local slope at that location. α i and tan(β i ) can be regarded as, respectively, a dimensionless index for the area draining through location i, and a measure of how fast water can be transported downslope through that location. Regarding a simulation grid as a watershed, we calculated the mean CTI value using Marthews et al (2015) for each elevation rank of each simulation grid. As there are negative correlations between CTI and elevation throughout the simulation area (mean r 2 = −0.286, n = 2678), the elevation rank partially reflects the difference in CTI within each grid.
The specific procedure for the implementation of TOPMODEL into the SEIB-DGVM is explained in the supplemental information (SI) 1 (available online at stacks.iop.org/ERL/15/095006/mmedia). With the redistribution of soil water among elevation ranks, the air temperature of each elevation rank is adjusted by subtracting 0.6 • C for every 100 m higher elevation than the grid average. The bottom soil temperature at a depth of 2 m was taken from the observation of permafrost temperature by Beer et al (2013) and was assumed to be constant. This assumption is only due to technical reason for reducing computation demands. We should note that it will result in the underestimation of the rate at which permafrost thaws under a warming trend, and it ignores possibility of rapid landscape changes (e.g. thermokarsting and plateau subsidence) those can happen as a result of permafrost decay (Walvoord and Kurylyk 2016).
In the original model, trees can die due to the low-growth-efficiency (carbon starvation) and bioclimatic-limits. We additionally introduced death by over-wet stress; 20% of individual trees die when the monthly mean soil water within 50 cm of the active layer exceeds 74% of its soil-saturation-point. This evaluation of over-wet stress is conducted at the end of each month during the leafing period of larch. This formulation is a result of model calibration based on the observational data. Detail procedures are described in SI2.

Data
The model requires daily climatic variables for the following items: mean air temperature, daily range of air temperature, precipitation, downward shortwave and longwave radiations, wind velocity, and relative humidity. For the calibration and validation of the model, we employed historical climate data from the period between 1991 and 2000, which was mainly delivered from CRU-TS4.00 observation-based climatic data (0.5 Degree, 1901 ∼ 2015) (Harris et al 2020). During this period, number of meteorological station contributed to the CRU data is relatively stable in Russia (Harris et al 2014), and this period overlaps ranges those major validation datasets cover (larch LAI: 1998, GPP: 1998. As CRU data are provided at monthly time intervals, daily climatic variability within each month was supplemented using NCEP/NCAR daily climate data (Kalnay et al 1996). For daily wind data, NCEP/NCAR reanalysis was employed for simulations. For radiations, we made corrections based on the NASA/GEWEX SRB Release-3.0 monthly dataset (https://eosweb.larc.nasa.gov/project/srb/srb_table). SI3 explains the specific and detailed procedure of data conversion.
In the simulations for the 21st century projection, we used an MIROC-AOGCM output under the RCP8.5 scenario (Watanabe et al 2011). The spatial resolution of this AGCM output was 192 × 94 global grids and was linearly interpolated to a 0.5 • grid mesh, which corresponds to the spatial resolution for the simulations in this study. To provide consistency for the historical climate data, the 21st century data of MIROC were converted using the 20th century data of CRU; for each climatic item by month, averages from 1991 to 2000 data were subtracted from the 21st century projections of the MIROC and the averages from 1991 to 2000 of the CRU were added. Because the AGCM data are also provided as monthly time interval, daily climatic variability within each month was supplemented in the same way as the historical climate data. Geographical distributions of annual precipitation and annual mean air temperature of historical and future conditions are shown in figure S1.
The model requires elevation, CTI, and permafrost temperature as forcing data (figure 1, table S1). For the calibration and validation of the model, we employed tree LAI, gross primary production (GPP), living biomass, active layer thickness (ALT), and runoff fractions (figure 2, table S2). Elevation, CTI, permafrost temperature, and ALT were converted to 1/112 • using the assumption that these environmental variables are linearly changed between proximate grids. This resolution corresponds to the Larch-LAI data, which have the coarsest resolution among these datasets.

Calibration and validation
For the calibration and validation of the model, we used a simulated forest of 100 years, which can be regarded as the average forest age of a Siberian larch forest, because the average interval of stand-replacing fire of a Siberian larch forest is estimated to be approximately 200 years (Kharuk et al 2013).
To calibrate and validate the model, we conducted a 100 year simulation by repeatedly inputting climate data from the last 10 years of the 20th century (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000). Atmospheric CO 2 concentration is assumed to be 368 ppm, which was the observed global average in the year 2000. Each 100 year simulation started from bare ground and assumed no wildfires during this period. The results of the last 10 years of the simulation were averaged and used for validation and calibration.
For calibration, we reconstructed simulationgrid-averages of fundamental properties of plant productivity and thermo-hydrology: annual GPP, canopy LAI in July, ALT, and runoff fractions. Biomass was not used for the calibration, because its estimates in Siberia varies quite much among analysis (Houghton et al 2007). The calibration parameters and their assigned values are presented in table S3. The comparison of observation-based-estimates and simulated values (averages of simulation grids) is available in table S4. To validate the calibrated model, we compared geographical distributions of calibrated variables and living biomass between observationbased-estimates and simulated values. The geographical distribution of runoff fraction was not employed here, as it shows artificially binominal distributions (figure 2(i)). The calibrated model roughly reproduced the geographical distributions of the calibrated variables and biomass (figures 2 and S2), which demonstrated that the model reasonably responds to environmental gradients in eastern Siberia. Correlation coefficients (CCs) between observation-basedestimate and simulated values were 0.622 for tree LAI in July, 0.802 for annual GPP, 0.635 for biomass, and 0.584 for ALT.

Simulation and analysis
With the calibrated model, three different simulations were conducted under both current and future climatic conditions: (1) control run, (2) simulation without redistribution of soil water among elevation ranks (hereafter, no redistribution run), and (3) simulation without air temperature adjustment among elevation classes (hereafter, no temperature adjustment run). Except above modifications, all runs were conducted in accordance with the calibration and validation procedures. The 100 year run is not enough for bringing soil carbon pools to a quasi-equilibrium state, however, soil carbon pools do not feedback to the vegetation dynamics and thermo-hydrology in this model.
The comparisons of the three simulations are referred to as 'model experiment' in this study. To demonstrate how model modifications enable the reconstruction of within-grid heterogeneity, we calculated the CCs between the elevation rank and larch LAI and between the elevation rank and ALT for each half degree grid cell.
We also conducted a sensitivity test for the assumption of fire return interval on the simulation results: 50 and 200 year old forests were simulated and compared with a forest stand 100 years from bare ground (control run). Except for the simulation period, both runs were conducted in accordance with the control run.
Finally, we conducted a simulation using a fore-  in the year 2001 and the results in the last 10 years (2091-2100) are averaged for evaluation. During the simulation, atmospheric CO 2 concentration changes along with the RCP 8.5 scenario. In accordance with the simulations under historical climate condition, no redistribution run and no temperature adjustment run were conducted under the future climate condition.

Results
A clear geographical pattern existed for the CC of elevation rank and tree LAI in the observationbased data. A strong negative trend was observed for mountain regions, whereas a weak positive trend was observed for plain regions ( figure 3(a)). The modified model qualitatively and partially reconstructed this pattern (figure 3(c)), although its magnitude was generally higher than that for the observation-based estimate. The simulated negative CC in mountain regions should be primarily reconstructed by air temperature adjustment among elevation ranks because it nearly disappeared in the no temperature adjustment run (figure 3(g)), whereas it was maintained in the no redistribution run (figure 3(e)). The simulated positive CC in plain regions should be reconstructed by the redistribution of soil water among elevation ranks because the trend nearly disappeared in the no redistribution run (figure 3(e)), whereas it was maintained in the no temperature adjustment run (figure 3(g)).
A geographical pattern also existed for the CC of elevation rank and ALT in the observation-based data. A negative trend was observed for mountain regions, whereas a positive trend was observed for plain regions ( figure 3(b)). The modified model qualitatively and partially reconstructed the negative CC in the mountain region, whereas it did not reconstruct a consistent CC for the plain regions ( figure 3(d)). The positive trend in the mountain region should be reconstructed by air temperature adjustment among elevation ranks because it disappeared in the no temperature adjustment run (figure 3(h)), whereas it was maintained in the no redistribution run (figure 3(f)). The predominant positive trend in the CC occurred in the no temperature adjustment run (figure 3(h)), which indicated that soil water redistribution was responsible. , and soil temperature, respectively, simulated over example plots of plain area, of which the location is shown in figure 3. This region includes 180 simulation grids, which are averaged for each elevation rank. Rightmost row shows seasonal changes in LAI simulated over the same area. Panels on top, middle, and low columns, show simulation results averaged over elevation rank 10 (highest), 5 (near average), and 1 (lowest), respectively. Solid lines on the soil temperature show active-layer-thickness. Unit of water content 'fraction' indicates fraction of soil volume.
To evaluate how these different patterns occur between mountain and plain regions, seasonal changes in thermo-hydrological status were compared between mountain and plain example areas. The locations are shown in figure 3. Roughly consistent with figure 3(c), tree LAI was higher at higher elevation ranks in the plain area (figure 4), whereas the opposite trend occurred in the mountain area (figure 5). For both the plain and mountain regions, the content of liquid soil water was higher at lower elevation ranks, whereas its maximum value at within tree root depths (0-0.5 m) was higher for the plain region than for the mountain region. The soil layer within thaw depth was near-saturated with soil water at the lowest elevation rank of the plain area. Roughly consistent with figure 3(d), thaw depth was shallower at higher elevation ranks for the mountain region (figure 5), whereas no clear trend occurred with elevation ranks in plain regions (figure 4).
In the grid averages of the simulation result under the historical climate condition, soil water redistribution increased tree LAI, GPP, and biomass, whereas it decreased the runoff fraction ( figure S3). The decrease of runoff with increase of plant productivity generally occurs in the simulations of the SEIB-DGVM, because higher LAI intercepts more rain fall, and higher photosynthesis rate consumes more soil water. These changes were only conspicuous in the plain regions. Soil water redistribution also increased ALT, except for in the southern mountain region. However, temperature adjustment did not bring noteworthy change to the grid averages (figure S3).
Comparison of simulated forests of different ages showed that older forests have larger tree LAI, GPP, and biomass, while slightly shallower ALT and slightly smaller fraction of runoff (figures 2 and S4). In the CC between elevation rank and tree LAI, its positive trend in plane regions became weak and eventually became negative in older forest (figure S5), although this change was much less apparent compared to those of LAI, GPP, and Biomass. In the CC between elevation and ALT, no substantial changes occurred with forest age ( figure S5).
Under the forecasted climatic condition during 2091-2100, tree LAI becomes denser for the whole simulation areas in the control run (figures 6(a) and (b)). The no temperature adjustment run did not have a noticeable difference from the control run (data not shown). In the no water-redistribution run, an intense decreasing trend occurs for the whole simulation areas, except for high elevation and northern coast regions around longitude E140-160 ( figure 6(d)). Under the control run of future climate, the lowest elevation rank of both plain (figure S6) and mountain (figure S7) regions are nearly saturated with soil water status in its ALT, which will kill trees. Higher elevation ranks have much lower soil wetness. Under the no-water-redistribution run of future climate, ALT saturation of soil water occurs in both plain (figure S8) and mountain (figure S9) regions irrespective of elevation ranks, which is responsible for their low tree LAI.

Reconstruction of tree LAI differences among elevation class
The model reconstructs the topography control on tree abundance; negative correlations between elevation rank and larch LAI in mountain regions (figure 3(a)) were partially reconstructed by the control run ( figure 3(c)). The model experiment showed that this reconstruction was a result of air temperature adjustment with elevation, of which the intensity was higher for the mountain region, with its high within-grid deviation of elevation ( figure 1(b)), than in the plain region. In the simulations of eastern Siberia, colder environment generally decreased simulated larch LAI through a shorter growing season, lower photosynthesis efficiency, and shallower ALT; and accompanied by lower water holding capacity in the soil (Sato et al 2016). However, the model did not reconstruct the positive correlation in the southernmost mountain region around latitude N56-59 in the simulation. This area had an approximately 10 • C higher annual-mean-temperature (figure S1(a)) than the central mountain region and more than 500 mm of annual precipitation (figure S1(c)). This climatic combination would inhibit the above mechanisms for the downregulation of tree LAI.
The model experiment showed that within-grid soil water redistribution was responsible for reproducing positive correlations between elevation and larch LAI in plain areas. This positive correlation should be caused by the combination of soil water redistribution and over-wet-kill of trees, of which the condition was calibrated with tree mortality observed at the Spasskaya Pad experimental forest during August 2006 (SI2). Our simulation demonstrated that the same extent of this moist condition can occur throughout simulated plain areas under the climate condition at the end of the 21st century.

Reconstruction of ALT differences among elevation classes
The observation-based data showed the trend that ALT becomes shallower in higher elevation areas in mountainous region ( figure 3(b)). This pattern is intuitively understandable because ALT becomes shallower in colder environments and the within-grid temperature difference should be much higher in the mountain regions, which has much higher withingrid heterogeneity of elevation, than in the plain regions ( figure 1(b)). The simulation roughly reconstructs this trend and the model experiment showed that this reconstruction is caused by air temperature adjustment with elevation. The opposite trend was observed in the plain regions; ALT becomes shallower in lower elevation areas ( figure 3(b)). Moister environments in lower elevation classes would explain this pattern. Observational studies on arctic Canada have demonstrated that wetland soil has a shallower ALT than drier sites and this pattern is caused by the fact that wetland soil has a larger ice content, which requires more energy to thaw. Additionally, a larger ice content facilitates downward heat conduction to permafrost; heat conductance of ice is approximately four times larger than that of liquid state water (Woo andXia 1996, Carey andWoo 1998). The model experiment generated a trend consistent with the above explanation; soil water redistribution among elevation classes generally results in shallower ALT in lower elevation areas (figure 3(h)), where ice saturation in soil happens during winter ( figure 4). However, the model failed to reconstruct the geographical distribution of the negative correlation (figure 3(b) and (d)).

Sensitivity test for forest age
We evaluated simulated forests at 100 years old, which can be regarded as half of the average frequency of stand-replacing fire in Siberia (Kharuk et al 2013). However, much uncertainty remains in the estimate due to the complicated nature of wild fires in this region. For example, the average of overall fire frequency is 15 years, ranging from 4 to 43 years, and has not been observed to exceed 50 years (Katamura et al 2009). However, about 80% (Conard and Ivanova 1997) to 93% (Isaev et al 2002) of the area burned is thought to be surface fires, which only consumes litter and surface vegetation. In permafrost areas, wildfires deliver significant impacts on the forest conditions as well as hydrology (Yoshikawa et al 2003, Jafarov et al 2013. Indeed, forest age gave substantial influences on larch LAI, GPP, and biomass in our simulations (figures 4 and S4). In contrast, forest age provided only small or negligible influences on topographical controls on larch LAI and ALT within grid scale (i.e. CC between elevation rank and larch LAI, and CC between elevation rank and ALT) (figure S5). Although each simulation assumes that all forests had the same age, this high independency would demonstrate that our validation strategy, which employs these CCs, is robust against the assumption for forest age.

Further effort required for better simulation
The model cannot reconstruct the geographical distribution of the positive correlation between elevation and ALT (figures 3(b) and (d)). Additionally, the simulated CCs between elevation rank and larch LAI were much higher than those of the observationbased data (figures 3(a) and (c)). For a more accurate reconstruction of thermo-hydrology and vegetation distribution, further efforts are required to address issues below.
Part of the reason for this limited accuracy can be attributed to the assumption that elevation class represents all environmental heterogeneity within a grid. Another part of the reason for the inaccurate reconstruction can be attributed to the assumption that all forests had the same age in our simulation. This is because vegetation can alter thermo-hydrology dynamics through processes that involve canopy shading, snow interception, and litter accumulation, which function as effective thermal insulants between the soil and atmosphere (Jorgenson et al 2010).
Indeed, observational studies have shown that seasonal maximum of ALT is the highest immediately after stand replacing fire and gradually becomes shallower with the recovery of forest (Sato et al 2016). This is consistent with the fact that stand replacing fires in the Siberian larch forest consume most of the surface soil organic layer, of which the depth is typically 10-25 cm in mature forests (Sofronov and Volokitina 2010). Consistently, in permafrost ecosystems in Alaska, canopy removal following fire increases both ground heat flux (Chambers andChapin 2003, Liu et al 2005) and ALT (Viereck 1982).
Interactions between thermo-hydrology and vegetation can occur at finer geographical scales (Walvoord and Kurylyk 2016); hence, it is challenging to simulate its behavior at the geographical scale of eastern Siberia. For example, Hinkel (2003) noted the need for specific small-scale studies to understand the causes of the significant intra-site variation in ALT and to simulate moisture conditions found at the larger scale. Indeed, Wright et al (2009) showed that ALT is highly variable over a relatively short distance (0.25-1 m) and spatial variability is strongly correlated to soil moisture distribution, which is partially influenced by lateral flow converging to frost table depressions.
Finally, the projected influences of near-surfacepermafrost thaw on the existence of the larch forest differ among studies. For example, permafrost decay results in enhanced plant productivity for Siberian larch forest with a modified LPJ-DGVM (Beer et al 2007), because a deeper ALT would increase the water holding capacity. However, in the dynamic vegetation model of Zhang et al (2011), permafrost thaw results in the collapse of the Siberian larch forest because of the rapid drainage of soil water, which causes a large water deficit in summer. Therefore, more effort is required to improve thermo-hydrological representations for a better projection of the Siberian larch forest in a changing climate.

Simulation under a forecasted climate condition
In the control run, water redistribution should reduce the risk of over-wet-kill for larch trees in higher elevation classes, whereas it should increase this risk in lower elevation classes. In our simulations, the former effect should be higher than the latter one because water redistribution among elevation classes slightly increased the grid average of larch LAI in the plain areas under the climate of 1991-2000 (figures 6(a) and (c)) and significantly increased the grid average of larch LAI in nearly all simulated areas in the climate of 2091-2100 (figures 6(b) and (d)). This significant difference between the 20th and 21st centuries demonstrates that the 'refugia' effect in higher elevation patches becomes substantial under the trends of increasing annual mean temperature and precipitation, which are forecasted throughout eastern Siberia during the 21st century (figures S1(b) and (d)).
In our previous study using the original SEIB-DGVM, which ignores both over-wetting-kill and topographic heterogeneity within simulation grids, larch biomass increases under the forecasted climatic trends of the MIROC-based projection under the RCP 8.5 (Watanabe et al 2011) as a result of growing season extension and reduced water shortages (Sato et al 2016). This simulation also shows an increasing trend of soil water content during the growing season of larch; hence, it may result in the reduction of larch biomass if over-wet-kill is taken into consideration without the consideration of topographic heterogeneity. This study demonstrates the risk of overestimating over-wet-kill without the consideration of topographic heterogeneity within simulation grids.

Conclusion
The modified model partially reconstructs the topography control on tree abundance; a negative correlation between elevation and canopy LAI in the mountain region but a positive correlation in the plain regions. Additionally, the modified model partially reconstructs topography control on ALT; a negative correlation between elevation and ALT in the mountain region but a positive correlation in the plain regions. The model experiment shows that the positive correlations in the mountain region are the results of air temperature adjustment with elevation, whereas the negative correlations are the results of the combination of soil water relocation along elevation and over-wet-kill of trees. These mechanisms, which contribute to the reconstruction of observation trends, may be consistent with the actual mechanism in nature. However, the model has poor performance for the qualitative reconstruction of the geographical distribution of topographic controls on canopy LAI and ALT, which indicates the need for further effort to develop a more accurate representation of the interactions between climate and vegetation conditions on the thermo-hydrology dynamics.
We demonstrated, in a forecasted climatic condition of the 21st century, that the future may hold a much moister environment throughout eastern Siberia, which is enough to cause ubiquitous tree death owing to over-wet stress. Modeling the overwet-kill of trees without the consideration of topography heterogeneity may result in the overestimation of forest collapse caused by over-wet-kill of trees under the expected climate trend in eastern Siberia. Applicability of this method in other permafrost ecosystems would also be worth exploring.