3D simulation of boreal forests: structure and dynamics in complex terrain and in a changing climate

Boreal forests contain approximately half of the terrestrial carbon (Gower et al 2001), which is mostly stored in soils and permafrost (Strauss et al 2013). However, changes in boreal forest vegetation structure, composition and productivity can affect how the below-ground carbon is exchanged with the atmosphere (Betts 2000, Bonan 2008, Hollinger et al 2010, Ma et al 2012, Groisman and Gutman 2013, Kuusinen et al 2014). For example, increased tree mortality from wildfires, insect outbreaks, and other disturbances, release carbon to the atmosphere. These disturbances and associated changes in species composition additionally can lead to albedo changes and more exposed soils, which may facilitate further release of carbon from soils and permafrost (Krankina et al 2005). Simultaneously, a warmer environment may facilitate increased productivity, especially in forests where precipitation and air temperature increase concurrently.

Modeling allows us to synthesize field data, to understand forest processes, and to extrapolate for conditions not currently observed, such as future climate. Any model is based on simplifications and assumptions, which carry with them a certain level of uncertainty. Our study employs a new spatially-explicit gap dynamics model, SIBBORK, which has been validated for forest composition and structure for southern- and middle-taiga ecotones (Brazhnik and Shugart 2015). SIBBORK has a high spatial resolution, and does not assume climate–vegetation links utilized in other models that have been used to assess potential changes in boreal vegetation (e.g., SiBCliM, Tchebakova et al 2009). SIBBORK does not include chilling requirements for growth and regeneration that could drive especially the southern boundary of species ranges, or frost tolerance that may limit the northward distribution. The species ranges in SIBBORK are reproduced based on interspecific competition in combination with species-specific tolerances to limitations on four environmental factors deemed of primary importance for forest vegetation processes (Pastor and Post 1986, Urban 1990): available light, soil moisture, soil fertility, and growing degree days. The model explicitly considers the fully-coupled three-dimensional light-environment including the effect of available light on tree growth, and the feedback with the amount of foliage affecting available light throughout the canopy. The model does address the multitude of other ways in which vegetation potentially affects the local microclimate. The current version of SIBBORK does not include permafrost or disturbances, whereas changes in both have been forecast for the near future in Siberia and may contribute to current treeline location and future migration. Nonetheless, SIBBORK contains several strengths that render it the appropriate model for this investigation. In contrast to many models that employ parameterizations at the plant functional type level (e.g., FORMIX3, Huth et al 1998), SIBBORK utilizes species-specific parameterizations, which include tolerances to environmental stressors and reflect differences in vascular and root structure. For example, larch is more drought-tolerant than spruce due to difference in root structures (Babushkina et al 2010), and pine is more drought-tolerant than broadleaf species due to differences in vascular structure (Bauweraerts et al 2014). The model outputs are at the individual tree level, in contrast to models that evaluate cohorts or size classes. Furthermore, SIBBORK does not require or assume landscape homogeneity, which allows for heterogeneous landscapes and associated gradients to be simulated. Vegetation responds to the environmental conditions on the landscape, and we are therefore able to resolve the differential responses of vegetation, which include differences in stand structure and composition, based on the simulated conditions on the landscape. Application of SIBBORK to the simulation of transition zones and the response of vegetation in these ecotones presents a new approach to evaluating treeline migration with changes in the thermal regime and simultaneously provides a rigorous test of this new model's capabilities.

There is uncertainty in how these changes in forest structure and composition will affect the role of Russian boreal forests in the global carbon cycle. To understand how the boundaries of the Siberian boreal forest may shift in the near future and to provide a rigorous test for a new gap dynamics model, we applied the new gap dynamics vegetation model SIBBORK to the investigation of potential changes at the northern and southern treeline transition zones in central Siberia, which has some of the fastest recorded warming rates (Gruza et al 2015). These shifts in temperature regimes are likely to drive change in forest structure in the near future, in as little as decades. This will likely have a significant effect on how carbon is cycled between the boreal ecosystem and the atmosphere.

A thorough description of this new gap dynamics model can be found in Brazhnik and Shugart (2015). Briefly, SIBBORK is an individual-based gap model, which tracks the establishment, growth, and mortality of individual trees on 10 m × 10 m plots, which are arranged in a grid along user-specified terrain. Tree location on the terrain is resolved at the plot scale for a maximum stand density. Since it retains plots and specifies environmental conditions at the plot-level resolution, SIBBORK is not a continuous space model sensu Grimm and Railsback (2005) for competition and other interactions. Tree establishment, growth and mortality are kept track of on an annual basis through the state variables of diameter at breast height (DBH) and crown base. The governing equations in SIBBORK address tree growth, the interactions among trees through shading, and the interactions between trees and their environment. The four environmental factors explicitly parameterized are growing degrees days above a base temperature of 5° (GDD₅), soil moisture as a function of precipitation, evapotranspiration and runoff, available light throughout the canopy as a function of foliage density, and soil fertility, which caps the annual gross primary productivity. The environmental conditions (incident radiation, air temperature, soil fertility) are specified at the plot-level, which significantly enhances the resolution of the conditions on the landscape and allows for simulation of vegetation response to changes in those conditions.

In a transition zone simulation, environmental conditions change along a transect or across a region. Artificially compressed gradients (i.e. radiation, temperature) can be imposed within the simulation domain. At the southern boreal treeline in central Siberia, the transition is along a steepened gradient due to the presence of complex terrain. On mountain slopes, transitions between ecotones occur along 100 s of meters, whereas similar transitions occur along 100 s of km on flat terrain (Holtmeier and Broll 2005, Soja et al 2007, Kharuk et al 2010a). Using the ASTER GDEM (METI and NASA 2011) resampled at 10 m × 10 m resolution, the incident radiation is computed using the Area Solar Calculator (Fu and Rich, 2002) in ArcGIS (ESRI v10.2 2014). An environmental lapse rate of 6.5 °C km⁻¹ is applied to extrapolate the temperature from the nearest weather station to the simulated landscape based on elevation (Polikarpov 1986). Downscaling of precipitation to the resolution of several meters is poorly understood, therefore, in the model, the entire simulated region receives the same monthly precipitation. Soil fertility, field capacity, and wilting point are specified for each plot based on soil type using data and soil nomenclature from Onuchin (1962) and Stolbovoi and Savin (2002).

A treeline is a transition between a region dominated by arboreal vegetation and one with a general absence of trees (Stevens and Fox 1991). These transitions occur at the northern or southern edge of tree dominance, or on a more local scale, such as along an elevational gradient or along the boundary of a body of water. The drivers of northern treeline are poorly understood, but likely include a combination of the following: the distribution of continuous permafrost (Kryuchkov 1973), the average position of the polar front (Pielke and Vidale 1995), the 10 °C–12 °C July isotherm (Anuchin 1986), the length of the growing season (Stevens and Fox 1991), incident radiation during the growing season (MacDonald et al 2000), winds (Holtmeier and Broll 2005), the balance between precipitation (P) and potential evapotranspiration (PET) (Hogg and Hurdle 1995), and other components of the local climate (Esper and Schweingruber 2004, Kharuk et al 2006, Sun et al 2011). In contrast, the southern treeline of the boreal forest experiences a large degree of climatic continentality, and has different potential drivers, which may include a combination of the local temperature/precipitation regimes (Chytry et al 2008), the 19 °C–20 °C July isotherm (Anuchin 1986), the P-PET balance (Kharuk et al 2013b), winds (Holtmeier and Broll 2010), fire regime (Ivanova et al 2010), and competition with shallow-rooted grasses, which may be better at utilizing the meager amounts of precipitation received during the growing season than deep-rooted, large-structured trees (Stevens and Fox 1991). Additionally, the microclimates associated with complex terrain, as well as resource and vegetation patchiness, can promote tree establishment and growth even when the surrounding conditions are not favorable to arboreal vegetation. Not all of these environmental factors are captured by SIBBORK, however, the simulation of transition zones, such as treelines, provides a rigorous test for the strength of SIBBORK in simulating three-dimensional terrain and associated gradients, and the resulting heterogeneous response of vegetation to the conditions on the landscape.

We selected the relatively data-rich northern treeline ecotone in the Karsnoyarsk Region of central Siberia for simulation with SIBBORK to better understand how the structure and dynamics of this vegetation may change with shifts in temperature and precipitation regimes. The topography at northern treeline is relatively flat (figure 2(a)), with poor drainage and low soil fertility on gleyzem soils (Stolbovoi and Savin 2002) underlain by continuous permafrost (Kotlyakov and Khromova 2002). Along the Yenisei River meridian, the treeline is observed between the World Meteorological Organization (WMO) weather stations of Dudinka (#23 074, 69.4° N, 86.18° E, 19 m above mean sea level or amsl) and Turukhansk (#23 472, 65.47° N, 87.56° E, 32 m amsl), 440 km to the south. Average monthly temperatures and precipitation totals were computed for these two stations using the meteorological record from 1936–1990 (NCDC 2005a, NCDC 2005b). The climates at the two stations set the conditions for the southern and northern edges of the simulated transect, with a linear interpolation of monthly temperatures in between to specify the thermal regime at the plot-level. The ArcGIS Area Solar Calculator (Fu and Rich 2002) was used to compute monthly radiation at the plot-level for Igarka (#23 274, 67.28° N, 86.32° E, 30 m amsl), located between Dudinka and Turukhansk, which was then validated against observations from the nearest World Radiation Data Centre (wrdc-mgo.nrel.gov) station in Olenek (68.30° N, 112.26° E). The northern treeline is discontinuous, dominated by larch (Larix sibirica), with a maximum productivity of up to 100 m³ ha⁻¹ (wood volume) found along the river and stream banks (Bartalev 2010, Bartalev et al 2011, Huttich et al 2014). Environmental conditions specified in SIBBORK simulation of northern treeline are listed in the top row of table 1.

Table 1.  Environmental conditions at northern and southern limits of the boreal forest in central Siberia.

LatN longE alt. (m)	Aspect	Slope	Plot size (m2)	Simulation domain (plots, ha)	Mean monthly precipitation (mm)	Dominant soil type(s)	Soil field capacity (cm)	Soil wilting point (cm)
69.4
86.4	None	Flat	100	12 × 181, 21.72	30–59	Gleyzems	62.5	31.25
19.0
51.4
94.2	North/south	sloped, ridge	100	12 × 581, 69.72	3–47	Cambisols, kastanozems	26.5	13.25
100–3000
Mean monthly temperature (°C)	Mean monthly precipitation (mm)	Mean monthly radiation (W m−2)	Mean monthly radiation (growing season) (W m−2)	Relative direct and diffuse solar radiation ratio (growing season)
−27.7 to +13.5	30–59	191	251	0.33/0.67
−30.1 to +19.9	3–47	285	303	0.5/0.5

Note: top row—topographic, climatological and edaphic conditions at Dudinka, near northern treeline.Bottom row—topographic, climatological and edaphic conditions near Kyzyl, in southern Alay–Sayan mountain region.

In contrast, the southern taiga treeline occurs in the complex terrain of the generally east–west-oriented Altay–Sayan mountain system (Ivanova et al 2010) and mountain ranges of northern Mongolia. Treelines are found for both the lower and higher elevational limits of forest. Although the north-western Altay–Sayan receives considerable precipitation on windward slopes (up to 1600 mm yr⁻¹; Chytry et al 2008), the south-central Altay–Sayan region, and specifically the Tannu-Ola, Sengilen, and Obruchev ranges are characterized by high continentality, with average July temperatures exceeding 20 °C in the valleys and basins, annual precipitation of 200–500 mm and PET of 150–300 mm yr⁻¹ (Samoilova 1973, Ivanova et al 2010). Both slopes are in a rain shadow, but the south-facing slopes receive more insolation, which drives higher PET demands and results in drier soils and vegetation dominated by grasses rather than trees.

The distribution of vegetation differs based on elevation and slope aspect, with coniferous forests dominating the north-facing slopes, and steppe and hemiboreal forests prevalent on south-facing slopes, called 'expositional forest-steppe' (Sedel'nikov 1979, Wilmking et al 2004, Dulamsuren et al 2005b, Kulagin et al 2006, Chytry et al 2008, Kharuk et al 2009, 2010a, 2010b, Timoshok et al 2014). Species composition differs among the different sub-regions within the Altay–Sayan Mountains. Although there is variability between individual mountain ranges, on average the south-facing slopes below 1500 m amsl are dominated by the steppe ecotone. An alpine belt, containing larch (Larix sibirica), Siberian fir (Abies sibirica), and Siberian cedar (Pinus sibirica, also called Siberian pine or Siberian Stone pine) stands is found 1500–2000 m, above which alpine tundra dominates. Along north-facing slopes, dark taiga, comprised of Siberian fir, Siberian spruce (Picea obovata), and Siberian cedar, extend up to 1750 m, with the subalpine belt found 1750–2100 m, and tundra above that (Gorchakovsky and Shiyatov 1978, Isachenko et al 1988, Istomov 2005, Chytry et al 2008, Kharuk et al 2013b), although in some locations the transition to tundra is found higher, at 2200–2400 m (Ovchinnikov and Vaganov 1999). The variability of vegetation belts along the Altay–Sayan mountain range is often on the order of 300–500 m and can be up to 1000 m (Monserud and Tchebakova 1996).

To simulate the southern treeline, an artificial digital elevation model for an idealized 3000 m tall 'mountain' with south- and north-facing slopes was generated. This DEM was used to compute the incident radiation in ArcGIS and the elevation-based temperature adjustment from the nearest WMO weather station located in Kyzyl (#36 096, 51.4° N, 94.23° E, 629 m amsl) using the 6.5 °C km⁻¹ lapse rate (Polikarpov et al 1986). GIS-generated monthly radiation was validated against observations from the nearest World Radiation Data Centre station in Irkutsk (51.16° N, 104.21° E). PET and GDD₅ were computed along the mountain slopes using plot-level monthly temperature and radiation input, and simulation-wide precipitation (200–500 mm yr⁻¹; Samoilova, 1973, NCDC 2005b). We did not adjust the precipitation received on the north- versus the south-facing slope of the idealized mountain in the simulation, and focused more on the average conditions representative of the Tannu Ola and Selingen ridges, which are most proximal to the southern treeline in the Krasnoyarsk region (figure 1).

Zoom In Zoom Out Reset image size

A soil fertility gradient was generated based on description of mountain soils in that region, with more fertile soils along the foothills, and shallower, less fertile soils toward the top of the ridge (Samoilova 1973, Chytry et al 2008, Chlachula and Sukhova 2011). Soil fertility in SIBBORK is specified as a cap on gross annual productivity (GPP). Simulated plots along the lower slopes were assigned a GPP cap of 6 t ha⁻¹ yr⁻¹, whereas GPP in plots at the top of the ridge was capped at 1 t ha⁻¹ yr⁻¹. The GPP limit was estimated based on the correlation between site index and forest productivity in forestry yield tables (Shvidenko et al 2006) and a map of biological productivity for Russia (Isachenko 1985). Soils in these mountains are mostly well draining cambisols, with kastanozems in the arid valleys where arboreal vegetation is generally absent (Samoilova 1973, Stolbovoi and Savin 2002, Ivanova et al 2010). Field capacity and wilting point were computed from a georeferenced database of soils (Stolbovoi and Savin 2002) for the average cambisols found on Sayan mountain slopes. Environmental conditions specified in SIBBORK simulation of southern treeline are specified in the bottom row of table 1.

To minimize computational expense, the spatial domain of the northern treeline simulation was limited to 2 km × 120 m (22 ha). A compressed thermal gradient was created via temperature adjustment to each row of plots within this swath of landscape. The generated July temperature gradient is shown in figure 2(a), and stretches between 7 and 16 °C. This temperature gradient corresponds to an approximately 700 km long transect along the Yenisei meridian, as determined from temperature interpolation between the Dudinka, Igarka, and Turukhansk WMO weather stations. Conversely, for the simulation of the southern treeline, an idealized 100–3000 m amsl mountain was created (6 km × 120 m, 70 ha), with half of the area on the north-facing and the other half on the south-facing slopes. Elevation-based adjustment to air temperature was applied, but the thermal gradient was not compressed, so the distances and elevation within the simulation represent actual distances and elevation gradients.

Zoom In Zoom Out Reset image size

For all simulations, the historical climate, which is described by monthly average temperatures and precipitation sums, was defined by the WMO pre-1990 record, as in Gruza et al (2015). Simulations were initialized from bare ground (no trees, but seeds available), with inseeding turned on for all species of interest to the timber industry in the Krasnoyarsk Region. The start year of the simulation was assigned based on the average age of the forest reported for each region, which is based on the mean fire return interval for the region of interest: ∼100 for northern treeline (Pautova 1976, Mitrofanov 1977, Deyeva 1985, 1987, Esper and Schweingruber 2004, Kharuk et al 2005, 2006), 60–80 years for southern Siberian mountains (Onuchin 1986, Houghton et al 2007, Buryak et al 2009, Ivanova et al 2010). Wildfire and other disturbances rarely occur at northern treeline, so the average age of the stands can be in excess of 100 years (Esper and Schweingruber 2004). In contrast, frequent disturbances, such as annual ground fires on the dry south-facing slopes (Ivanova et al 2010), significant variability in crown fire frequency over the 20th century, and an increased occurrence of insect outbreaks (Kondakov 2002) in the Altay–Sayan mountains maintain a younger forest, but the more moist north-facing and upper slopes have some old trees in excess of 400 years (Ovchinnikov and Vaganov 1999, Istomov 2005).

To approximate the average stand age in the simulation, estimated based on the fire return interval for each ecotone, the first year of the simulation was initialized from bare ground corresponding to the calendar year 1890 and 1910, in the simulation of the northern and southern treeline ecotones, respectively. In the year 1990, warming trends were applied as seasonal linear monotonic increases to the simulated monthly temperature based on the observed seasonal warming trends for each region (Buermann et al 2014, Gruza et al 2015). Both, historical and warmed air temperatures were adjusted at the plot-level for elevation relative to the nearest weather station. The historical climate parameters and the warming trends are specified in tables 1 and 2, respectively. Model outputs from 150 replicates were recorded at 5 year increments. The height and biovolume, both of which are function of DBH, were computed for each tree on each plot and averaged across the replicates. Raster maps of the most frequently dominating species (mode), as well as Lorey's height and biovolume, were generated at plot-level resolution. (Note: biovolume is a field-estimated parameter for the volume of the tree trunk, commonly reported in forestry yield tables and field surveys. This is in contrast to biomass, measurements of which require destructive sampling, and is often computed from field-estimated biovolume using bulk wood density, instead.) The precipitation regime was not altered alongside warming during these simulations, as no appreciable changes in growing season precipitation were observed in the Altay–Sayan mountains since the 1970 s (Kharuk et al 2013a, Gruza et al 2015), and although annual precipitation has somewhat increased at the northern treeline, this transition zone is not water limited and increase in soil moisture is not likely to significantly affect the arboreal vegetation growing there.

Based on very simple parameterization of climatological and environmental conditions, in which the plots receive an enhanced temperature adjustment based on the distance from the nearest weather station, creating a steepened temperature gradient across the northern treeline ecotone, equal receipt of precipitation on each plot, and GIS-based incident radiation for the latitude of interest, SIBBORK generates annual PET and GDD₅ at the plot level. In this simulated environment, SIBBORK plants saplings and grows trees based on species-specific tolerances and silvicultural information. In this case, SIBBORK planted and grew predominantly larches, the height, biovolume and density of which decreased to the north. North of a certain point along the air temperature gradient, no saplings were established. This constitutes the northern forest limit—the treeline—in the simulation (figures 2(b) and (c)).

The simulated treeline coincides with the 12 °C July isotherm and the 30 cm yr⁻¹ PET isoline in the simulation, similar to the conditions observed at the northern treeline (Kolosova 1982, Malysheva 1993). Trees become stunted north of 500 °C GDD₅. The established larch saplings are stressed due to insufficient heat during the growing season and succumb to stress-related mortality within the first few decades of life north of 400 °C GDD₅, and are not able to establish north of 300 °C GDD₅. This latter cutoff is in agreement with the Sayan Mountain model (Monseud and Tchebakova 1996), which also draws a division between tundra and taiga vegetation at 300 °C GDD₅. South of 500 °C GDD₅, canopy dominant larch trees achieve 20–24 m in height, with average biovolume of 102 ± 39 m³ ha⁻¹, and an admixture of spruce begins to appear. This compares well to observed biovolume of 75–100 m³ ha⁻¹ in larch stands between the 14 °C and 15 °C July isotherms (Bartalev, 2010). The likely reason the simulated biovolume appears a little high is because the small saplings are not excluded from the biovolume calculation. Bartalev (2010) does not report the minimum DBH used for biovolume calculation in his database.

The structure and composition of the simulated forest along these thermal transitions are shown in figures 2(b) and (c). If there are no disturbances and the forest is allowed to grow for several centuries, spruce plays a more dominant role south of 500 °C GDD₅ by year 400 (this corresponds to region 5 in figures 2(b) and (c)). To the south of that, mixed dark conifer taiga prevails. In the simulation of the northern treeline, trees essentially stop growing where they are supposed to. Species ranges are not specified in SIBBORK. Species stop growing whenever the environmental conditions are not favorable to seedling establishment, or where the trees are too stressed to obtain >10% of age- and species-specific optimal diameter growth and succumb to stress-related mortality, as is the case north of the treeline. To the south, species (e.g. larch) may be outcompeted by other species (e.g. spruce), for which the site conditions are more favorable. GDD₅ is not the only driving factor for the northern treeline and species distribution, but here provides a systematic way to describe the orientation of stands of different structures and composition along a (compressed) thermal gradient. The simulated species distribution closely resembles the observed species distribution between the tundra and Turukhansk, and simulates the location of the northern treeline very close to the observed treeline along the Yenisei River meridian, as shown in figures 2(c) and (d).

Within the first few decades of warming corresponding to year 1990–2020, the sparse, unstable, severely stunted larch saplings at the northern treeline are able to establish farther north, up to the new 400 °C GDD₅ isoline, with unstable saplings extending north to the new 300 °C GDD₅ isoline. Colonization to the north appears almost immediate, on the order of years to decades, as conditions change to allow regeneration. This may be an artifact of the sapling establishment subroutine, which plants trees at a DBH of 2.5 ± 0.25 cm. Siberian larch of this diameter can be 4.7–6.4 m tall and 16–40 years old, depending on site index, which makes it seem as though trees, not saplings, appear beyond the treeline mere years after the warming trend is applied. This rapid response could potentially be delayed through modification to the inseeding subroutine or if spatially-explicit dispersal processes were included in the model (Kharuk et al 2006). However, not enough data are available on small saplings with DBH less than 2.5 cm in this region for validation of a new sprouting subroutine. With regards to vegetation responses in the spatial and temporal domains, regeneration of larches has been observed >3 km from the nearest treeline (Kharuk et al 2006), and rapid treeline equilibration with environmental conditions is seen in some paleoecological as well as in more recent records (Clark 1998, Kullman and Kjallgren 2006), thus, the limit on seed dispersal may not need to be as stringent. As warming continues over the course of 150 years, however, the treeline in the simulation does not continue to advance northward. Instead, the stands south of the treeline accrue biomass, and a gentle decrease in maximum tree height is observed between what used to be the treeline under the historical climate, where full mature trees now grow, and the new treeline.

For the simulation of the southern boreal boundary, SIBBORK computes PET and GDD₅ on each plot based on temperature extrapolated from the nearest WMO station based on elevation and radiation computed for the north- and south-facing slopes of a given steepness. Slopes of both aspects receive the same amount of precipitation. SIBBORK appropriately simulates a drier environment on the south-facing slope and a 12 °C July isotherm around 2000 m amsl on the north-facing slope. Because of differential environmental conditions, different species are able to establish and become dominant at different elevations on the two opposite slopes. On the south-facing slope, no trees occur in the simulation up to 1400 m. Between 1400 and 2200 m, larches dominate, with stunted and stressed larches occurring close to the lower (1400–1500 m) and upper (2100–2200 m) elevational treelines. Larches at the lower tree line are stressed by insufficient soil moisture, whereas at the upper timberline they are stressed by insufficient heat. Above 2200 m no arboreal vegetation is simulated. On the north-facing slope, the vegetation is starkly different—dark conifers (fir, spruce, Siberian cedar) dominate the slope up to 1800 m. The Lorey's height of the mixed conifer forest reaches 22–24 m at lower elevations, and gradually decreases to 5 m at 1800 m. Above 1800 m, larches gain dominance, but these trees are significantly shorter (3.7 ± 0.7 m) than the stands observed on the south-facing slope at the same elevation. Above 2040 m the larch populations are unstable, severely stunted, and stressed by insufficient GDD₅. No arboreal growth occurs in the simulation above 2180 m on the north-facing slope. The simulated elevational distribution of boreal vegetation is remarkably similar to observed vegetation belts in the Altay–Sayan mountain system (table 3).

Within the first few decades of warming (this corresponds to years 1990–2030), the lower treeline on the south-facing slope retreats uphill up to 2000 m. Larches are able to establish sparsely at 1600–2000 m, but they are drought-stressed and do not survive more than a few decades. Meanwhile, on the north-facing slopes, the larch stands grow taller, and the admixture of spruce in the 1600–2000 m belt increases. As warming continues over the course of 150 years (this corresponds to years 2110–2150), the spruces almost completely replace the larches on the north-facing slope, while on the south-facing slope the lower treeline retreats to 2000 m. By the year 2100, the upper treeline advances to 3000 m, which is the highest elevation simulated.

Overall, the SIBBORK simulation of the environmental conditions, as well as vegetation structure and composition at the northern and southern boreal treeline ecotones matches field observations.

Larch dominates the northern treeline in central Siberia. The treeline is not a continuous wall of forest, but rather a gradual decrease in stand density from south to north, until just a few individual trees or clumps of trees are surrounded by tundra (Aksenov et al 2002, Holtmeier and Broll 2005). In the simulation, the northernmost stands are sparse and stunted, with tree height and biovolume increasing southward along the simulated transect. Just to the east of Dudinka, along the lower slopes of the Putorana Plateau, stands with average heights of 6.2 m and biovolume of 4–33 m³ ha⁻¹ have been observed, albeit for birch stands. In our simulation, northern treeline is comprised of larch stands with Lorey's height of 7.2 ± 2.6 m, and biovolume of 6–10 m³ ha⁻¹. Biovolume for larch-dominated treeline along the Yenisei River meridian is not available, and often only described as 'sparse larch' (Bartalev 2010).

The gradual decrease in stand density is appropriately simulated by SIBBORK, with few individual trees separate from the low density open canopy forest (region 2 in figure 2(b)), and biovolume increasing from 7 ± 1.3 m³ ha⁻¹ near the treeline to 45.9 ± 14.2 m³ ha⁻¹ at the 14 °C July isotherm. Spruce is not encountered until the transition from very poorly drained gleyzems soils to moderately well drained podzolic soil, 200–300 km south of the northern treeline along the Yenisei River (Stolbovoi and Savin 2002, Bartalev 2010, Bartalev et al 2011, Huttich et al 2014). Spruce begins to appear just north of the 16 °C July isotherm in the simulation, and dominates stands south of the 16.3 °C July isotherm that are more than 400 years old. Spruce in 200-year old simulated stands have a biovolume of 104 ± 78 m³ ha⁻¹ (not shown). In younger stands (<200-years old), spruce appears as an admixture just north of the 16 °C July isotherm, where at least 600 GDD₅ are accumulated annually. This agrees well with observations of 100–200-year disturbance regimes (predominantly fire) that maintain younger, larch-dominated stands (Pautova 1976, Mitrofanov 1977, Deyeva 1985, 1987, Esper and Schweingruber 2004, Kharuk et al 2005, Kharuk et al 2006), and observations of spruce stands with an average biovolume of 100–125 m³ ha⁻¹ at the northernmost range of spruce along the Yenisei River (Bartalev 2010), which is situated near the16 °C July isotherm.

Although Monserud and Tchebakova (1996) examined the thermal requirements for the tundra, subalpine taiga, montaine taiga, and dark-needled montane taiga in southern Siberia, the thermal limits for the ecotones are also likely to be valid for the northern treeline thermal gradient. They found that tundra generally exists in regions where the total thermal accumulation is less than 300 °C GDD₅, sparse and subalpine vegetation is found between 300 and 550 °C GDD₅, and dark taiga occurs in regions that acquire at least 750 °C GDD₅ during the growing season. The transitions between no arboreal vegetation and sparse larch stands occur in the simulation around the 300 °C GDD₅. The stands, dominated by larch, just as in the subalpine taiga ecotone, become thicker and taller between 300 and 500 °C GDD₅, at which point a small admixture of spruce (a dark needle-leaf conifer) begins to occur, with almost complete spruce and mixed dark taiga dominance south of the 16.3 °C July isotherm, which is approximately collocated with the 750 °C GDD₅ isoline. The model simulation appropriately reproduces the observed species distribution along a thermal gradient. Field observations from the northern treeline ecotone in Siberia reveal that individual trees can establish few (Kharuk et al 2006) to tens of kilometers (MacDonald et al 2000) from the nearest tree stand. These observations facilitate the interpretation of fast (decades) colonization of the tundra upon warming in our simulations, which do not include a spatially-explicit seed dispersal subroutine.

Northward expansion of sparse larch stands into the tundra, and the increasing density and crown closure of what used to be sparse stands near the northern treeline have been observed in various regions of northern Siberia, with rates of treeline advancement of 3–10 m yr⁻¹ near Khatanga (Kharuk et al 2006, Sun et al 2011), 1–4 m yr⁻¹ in the Polar Urals (Devi et al 2008), 4–14 m yr⁻¹ in the Polar Urals (Kharuk and Fedotova 2003), and comparable rates at numerous sites between those two locations (Esper and Schweingruber 2004), with likely increase in colonization rate as climate warms and as the new young treeline individuals reach reproductive maturity (Kharuk et al 2006) and are able to provide a more sheltered environment for establishment of new saplings. If the reproductive processes are temperature controlled, rapid advance of the northern treeline is likely (Grace et al 2002). Saplings that have established north of the 19th century treeline have experienced slow growth, often growing no more than 4 m tall in 50 years (Esper and Schweingruber 2004). SIBBORK also simulates slow growth of northernmost saplings, which reach 5.3 ± 1.0 m over the same time frame. Increasing crown closure has been observed near the historical northern treeline (Kharuk and Fedotova 2003, Kharuk et al 2006). SIBBORK appropriately simulates an increase in foliage density in trees located at the historical northern treeline following a period of warming.

The mountains of southern Siberia stretch for over 10 650 00 km² and include the Altay–Sayan ranges, as well as the Zabaikal'e region to the south and east of Lake Baikal, where forests and forest-steppe cover 80% of the territory (Samoylova 1973). Along the southern boundary of the boreal forest, ribbons of forest appear predominantly on north-facing and upper slopes in the Alay–Sayan mountain region. The forest belt, thus, has a lower treeline limit on south-facing slopes dictated by soil moisture conditions, and an upper timberline on slopes of all aspects driven by temperature (Samoilova 1973, Dulamsuren et al 2005a).

In SIBBORK, we are able to monitor which environmental factor exerts the greatest limitation on growth of each individual tree. In the simulation, the trees along the upper treeline are most limited by insufficient growing degree days, whereas the growth of trees along the lower treeline on the south-facing slope is most severely limited by insufficient plant available water and too great of a fraction of growing season spent in drought, where soil moisture is below the wilting point.

Larch dominates both the lower and upper treeline on south-facing slopes, with occasional admixtures of pine (Pinus sylvestris) and birch (Betula spp.) (Babushkina et al 2010). Stand density and tree height near and at the forest boundary are smaller than in the middle of the forest belt (Samoilova 1973), and our simulations appropriately reproduce this observed patchiness near the upper and lower treelines (figure 3). The larch stands at the lower treeline are stressed, and this boundary is likely to recede as temperatures warm without corresponding increases in precipitation (Dulamsuren et al 2009, 2011). A 300 m upward retreat of the lower treeline due to increased tree mortality was recently observed in the Altay mountains in just 6 years (Kharuk et al 2013b). Germino and colleagues found that higher sunlight, such as may be available at the upper or lower treelines on south-facing slopes, compared to north-facing slopes or flat terrain, in combination with warm temperatures and water stress, inhibit the growth of Picea and Abies saplings in montane regions (Germino et al 2002). These species are generally not found on south-facing slopes in the Altay–Sayan mountains, and in the Zabaikal'e region just to the east (figure 1). In our simulations, these species are appropriately absent from south-facing slopes, which are dominated by larch, with small patches of pine and birch (figure 3). Lorey's height for northern-aspect mid-slope forests in the simulation exhibit a range of 7–17 m (mean ± std: 9.9 ± 2.2 m), which compares well with observed heights in birch-fir-pine forests in the Khamar–Daban mountains of the Zabaikal'e region, just to the east of the Altay–Sayan complex (range: 7–14 m, mean ± std: 10.8 ± 3.6 m) in 40–60 year old stands (Onuchin 1986). Additionally, simulated forest biovolume on south- (range: 8–150 m³ ha⁻¹) and north-facing (range: 7–49 m³ ha⁻¹ for 80-year old stands; 100–181 m³ ha⁻¹ for 200-year old stands when plots with young saplings regenerating in recently-formed gaps with biomass of less than 1 kg/plot are excluded from calculation) slopes correspond well with observed biovolume in the Selingen range (range: N 100–175 m³ ha⁻¹; S sparse-150 m³ ha⁻¹; sparse denotes <100 m³ ha⁻¹, Bartalev 2010). Figure 3 represents the simulated species distribution along south- and north-facing slopes in the mountains of southern Siberia.

Zoom In Zoom Out Reset image size

In the early 1990s, forests of the Krasnoyarsk region were assessed to be a weak sink for atmospheric carbon (Krankina and Dixon 1994), however, with increasing air temperatures not accompanied by any significant changes in precipitation in the southern part of this region, forest vegetation is likely to be stressed from temperature-induced drought, especially in southern Siberia (Allen et al 2010), and the increased mortality is likely to shift the role of these forests into sources of atmospheric carbon. As warming forces forests to recede higher uphill (migration of the lower timberline), but the upper timberline advances at a slower rate, there is the likelihood of significantly decreasing forest cover in this region. It has been hypothesized that increasing PET demands as a result of increasing air temperatures accompanied by only meager, if any, increases in precipitation, will stress and decrease the productivity of the boreal forests in Siberia (Briffa et al 1998, Knorre et al 2006, Allen et al 2010), Europe (Dirnbock et al 2003), Canada (Hogg and Hurdle 1995, Ma et al 2012), and Alaska (Bonan and Shugart 1989, Barber et al 2000), and potentially result in biome shifts (Sitch et al 2008, Buermann et al 2014).

Regeneration and advancement of arboreal alpine vegetation in continental montane areas may be more sensitive to changes in the soil moisture than to the thermal regime (Holtmeier and Broll 2005), and as the climate in southern Siberia becomes warmer and drier, grasses may significantly outcompete tree saplings, especially on dry soils (Ivanova et al 2010). Others have forecast the rapid (decades to centuries) upslope migration of the lower and upper treelines (Tinner and Kaltenreider 2005), as well as the migration of the southern boreal forest boundary northward over the next few decades (Tchebakova et al 2009). The upward retreat of the lower treeline on the south-facing slopes in the Altay–Sayan mountains has been observed in the last decade, most pronounced on south- and southeast-facing slopes, where forest receded upslope by over 200 m in just 7 years (Kharuk et al 2013b). Similarly, increased drought-induced mortality of birch was observed in the forest-steppe region on south-facing slopes in Transbaikalia mountains to the east of Lake Baikal and only a few hundred kilometers east of the simulated region (Kharuk et al 2013a), effectively shifting the lower treeline upslope by 22–40 m within the first decade of the 21st century. Increased drought-induced mortality has been observed in aspen (Populus tremuloides) stands in the Rocky Mountains in the US (Anderegg et al 2013) and Canada (Hogg et al 2008) with similar rates of retreat.

With regards to the upper timberline, a shift of 160 m since the 1890s has been noted (Istomov 2005), which expands on the shift of 100 m observed in the Western Sayan Mountains between 1976 and 2000 (Kharuk et al 2010a), with the greatest regeneration above treeline observed on north- and northwest-facing slopes. Other field studies report similar rates of upward timberline migration in the latter half of the 20th century along the southern boundary of the boreal forest in Russia (63 ± 37 m in 42 years, Kharuk et al 2010b) and Scandinavia (165 ± 20 m in 50 years; Kullman 2000, 2001, 2004). In addition to the treeline shifts, there are also shifts of species ranges along the slopes, with SIBBORK appropriately reproducing the upward shift in Siberian cedar and fir populations, the beginnings of which have already been observed in Altay (Moiseev 2002). Shifts in optimum elevation of plant species, which focuses on the center of distribution range, rather than just the upper or lower montane limit, of 2.94 ± 1.09 m yr⁻¹ have been observed in the mountains of western Europe since 1971 (Lenoir et al 2008), showing that the observed warming with no appreciable change in precipitation regime affects the entire range of species, and not just the range boundaries. Furthermore, increased stand density and crown closure has been observed at the upper alpine treeline, sometimes in addition to, and sometimes in lieu of pronounced upslope timberline advance (Kullman 2001). This change in stand structure has also been appropriately simulated by SIBBORK.

Previous modeling studies have suggested that a significant loss of forest is likely to occur in mountain areas with a +2 °C change in average annual temperature in Europe (Dirnbock et al 2003) and North America (Chapin et al 2004), especially when the precipitation regime does not concurrently change or becomes drier. The +2 °C also appears as an important threshold for the substantial decrease in forest-covered area along the simulated mountain slopes of southern Siberia in our study. The area occupied by pine-larch forests in Tuva (southern Sayan) has decreased by 70% since 1990, due to increased tree mortality and fire activity, and decreased regeneration post disturbance (Buryak et al 2009, Ivanova et al 2010). Retreat of forests from the simulated southeastern Altay–Sayan region (Tuva basin) has been forecast based on observations of warming and drying conditions in the region (Tchebakova 2006). The replacement of forest-steppe with steppe at the lower treeline has already been observed following fire and permafrost retreat (Monserud and Tchebakova 1996). This ecotone shift is likely to persist with warming and an increase in drought mortality and fire frequency in the coming decades, and is simulated appropriately by SIBBORK, even though the simulation does not include explicit parameterizations of disturbances and permafrost.

Numerous environmental factors that may affect treeline location, such as permafrost, winds, snowpack, and seed dispersal, are not included in this version of SIBBORK, so it is entirely possible that the northern treeline, which appears to be a function of air temperature in the simulation, may have completely other drivers. The model essentially simulates the theoretical treeline, based on environmental conditions, and as those conditions change, the colonization of areas that become favorable or retreat from regions that become no longer favorable may occur significantly faster than observed in nature, as the model does not address the increased susceptibility of young saplings to stress, compared to mature individuals, and does not explicitly represent seed dispersal.

Crawford and colleagues suggested that the northern treeline in Siberia may actually retreat south with warming due to increased soil and permafrost melt resulting in paludification (Crawford et al 2003). Since permafrost is not included in the current version of SIBBORK, the possibility of this feedback is not represented in our simulations. Increases in stand density and canopy closure have been observed at the upper timberline in southern Siberia (Kharuk et al 2010a, b). In SIBBORK, however, the maximum stand density limit is static and specified in the model driver at the plot level. Although the simulation can contain less than the maximum number of stems per plot due to unfavorable conditions, the possibility of future environmental conditions facilitating denser stands may not be sufficiently accommodated by the static limit on maximum stand density. Wind is an important factor in maintaining not only the northern forest boundary, but also the upper timberline (Resler et al 2005, Kullman and Kjallgren 2006), and lack of parameterization for this factor may result in establishment of saplings in sheltered areas on the terrain not due to shelter from the wind, but rather due to decreased incident radiation and, therefore, PET demands, again, missing the exact combination of driving factors for the location of the treeline, but getting the gist of it with the limited environmental factors that are explicitly resolved in the simulation.

Our model does not consider the differences in allometric relationships observed between trees within a forest compared to those on the exposed treeline edge (Hogg and Hurdle 1995). It is assumed that the environmental conditions along the treeline edge will differ from those within the closed canopy, and the growth of treeline trees will be adjusted in the simulation based on the stress those individuals experience. Kharuk and his colleagues (Kharuk et al 2013a) have found that in areas with increased tree mortality from temperature-induced drought stress, vegetation underlain by patches of permafrost was not adversely affected by increased temperatures, because the thawing of the permafrost active layer provided sufficient soil moisture throughout the growing season. Our simulation can be improved to reflect this heterogeneous response using a spatially-explicit parameterization for permafrost and associated soil processes. Additionally, with more local precipitation data available from field sites in southern Siberia (Monserud and Tchebakova 1996, Kharuk et al 2009, 2010c, 2013b), it will be possible to adjust precipitation by elevation and slope aspect, which would improve the resolution of local environmental conditions and the associated vegetation response to changes in these conditions. Disturbances, such as wildfire and insect outbreaks, are not included in this version of the SIBBORK model. An increase in the frequency and extent of these disturbances has already been observed across the Eurasian boreal forest (Buermann et al 2014), and their incidence is expected to continue to increase. Presence of disturbance shifts the forest age structure to a younger forest, with greater dominance by pioneer species, such as birch and pine. Repeated fires inhibit regeneration of arboreal species, promote dominance of grass species (Kukavskaya et al 2013), and expedite the shift from a forest to a grassland or steppe biome (Kukavskaya et al 2014).

Our investigation extends and supports the forecasts by other modeling studies (Scholze et al 2006, Sitch et al 2008, Tchebakova et al 2009). More specifically, we extend the work of Tchebakova and her colleagues (Tchebakova et al 2009), which forecast changes to Siberian vegetation by 2020, and drastic changes in forest structure and the shift from the boreal forest to the steppe biome for the southern half of the Siberian boreal forest by 2080. Our model SIBBORK forecasts the same biome shift within the same timeframe, even though the model structure and input parameters for our SIBBORK model differs significantly from the static climate–vegetation equilibrium assumption of the Siberian BioClimatic Model SiBCliM (Tchebakova et al 2009). However, as SIBBORK has a much higher spatial resolution than SiBCliM, we are able to resolve that it is primarily the forests on south-facing slopes that will be affected by drought-induced mortality, whereas forests on north-facing slopes and on slopes that are windward to moisture-carrying airmasses (west-facing along Yenisei Kryazh, northwest-facing Altay), will change in structure and composition but will not be replaced by the steppe biome in that time frame. Montane ecosystems are considered particularly vulnerable to climate change (Fischlin et al 2007). Tchebakova and her colleagues also developed and applied the Sayan Mountain model to predict the potential vegetation in the mountainous regions of southern Siberia, with 0.17° × 0.17° spatial resolution. Our investigation extends their work through application of our model toward understanding vegetation shifts in the southern mountains of Siberia at a higher spatial resolution (10 m × 10 m × 1 m). Additionally, SIBBORK is able to simulate three-dimensional terrain and the distribution of vegetation along topography-associated environmental gradients, as well as changes in environmental conditions over time. Field observations have shown that the changes in forest structure and treeline location strongly depend on topographic features (Kharuk et al 2010b). SIBBORK offers a unique understanding of the heterogeneous response of vegetation in a compressed transition zone along an elevational gradient to spatially and temporally heterogeneous changes in environmental conditions.

The SIBBORK gap dynamics model keeps track of the establishment, growth, and mortality of millions of trees across tens of hectares of artificial or real terrain, and simulates heterogeneous response of vegetation to the conditions on the landscape. SIBBORK appropriately reproduces transition zone vegetation, with the structure and composition closely approximating the observed forests at the northern and southern boreal forest limits, including the upper and lower elevational treelines in the complex terrain of southern Siberian mountains. The ability to appropriately simulate vegetation distribution in central Siberia and at the transition zones, not only provides a rigorous test of the model's capabilities, but also helps us understand the potential response of these ecosystems with continued warming. Alarmingly, the loss of forest at the southern treeline in the simulation far exceeds the rate and amount of forest gained at the northern and upper elevational treelines, however, this is consistent with observations from Siberia, Scandinavia, and North America. A spatially-explicit gap dynamics model that is able to simulate forest structure and composition across tens of hectares provides a strong platform for understanding what near-future changes in forest vegetation and the associated carbon cycling are likely in different regions of the Siberian boreal forest.

Observed warming trends in the Siberian boreal forest exceed predicted trends for this time period (IPCC 2013, Gruza et al 2015). It is likely that not only will the forest vegetation be responding to larger climate changes than predicted, especially in mountainous southern Siberia (Tchebakova 2006), but also that the rate of response may be faster than anticipated. Boreal forests in Siberia have experienced significant decline, both in the wood stock (m³ ha⁻¹) and area covered since 1980s (Shvidenko and Nilsson 2002). The observed increase in wildfires and insect outbreaks in recent decades, which has resulted in the dieback of millions of hectares of forest in southern Siberia (Soja et al 2007) makes this investigation particularly pertinent. Furthermore, Hansen and colleagues have found that during the first decade of the 21st century, the ratio of forest loss to forest gain was 2.1 for boreal forests, with Russia experiencing greatest forest loss compared to the rest of the globe (Hansen et al 2013). This forest dieback has occurred with less than +2 °C increase in annual average temperature, demonstrating that the boreal forest ecosystems may be much more sensitive to warming than previously anticipated.

The authors are thankful to Mr Whittlesey for programming consultation and helpful comments on the earlier versions of this manuscript, and to Dr C Huttig for providing the dataset Hybrid Biomass and Land Cover Map for Central Siberia (Bartalev, 2010), portions of which were used in model validation and generation of figure 1 and figures 2(a) and (d) in this manuscript. We also express our gratitude to the two anonymous reviewers, whose comments and suggestions helped improve the manuscript.