The significant role of snow in shaping alpine treeline responses in modelled boreal forests

Treelines across the Northern Hemisphere are shifting upward and northward in response to global warming, particularly in boreal forests, where climate change progresses more rapidly at high elevations and latitudes. These shifts intensify competition for resources, threaten endemic alpine species, and disrupt established ecological relationships, leading to biodiversity loss. However, significant heterogeneity and regional variation exist in how treelines respond to environmental changes, with many underlying drivers and constraints still poorly understood.

This study aims to enhance understanding of treeline dynamics at three alpine study sites – located in Canada, Russia, and Alaska – and to improve vegetation model predictions under changing climatic conditions. We evaluated the relative impact of key factors influencing treeline migration velocity and examined the effects of varying snow regimes on treeline migration within the alpine treeline ecotone. To achieve this, we incorporated a novel snow module into the vegetation model LAVESI (Larix Vegetation Simulator), enabling the integration of precipitation outside the growing season, snow accumulation, and snowmelt processes. This module allows for explicit modelling of the positive and negative impacts of snow depth on tree growth and treeline migration, while accounting for stochastically occurring extreme events and capturing full weather variability.

Our findings reveal site-specific responses to factors driving treeline shifts and forest expansion, with localised conditions playing a critical role in shaping migration dynamics. The Canadian and the Russian sites demonstrate clear insights into primary migration drivers, while the high variability at the Alaskan site indicates more complex local dynamics and greater predictive uncertainty. The study highlights the significant role of snow in modulating migration potential, as snow accumulation creates favourable conditions for seedling germination and growth while also posing risks of increased mortality from snow loads or avalanches. These results underscore the importance of incorporating snow-related processes into vegetation models to improve the accuracy of predictions for boreal forest dynamics.

Overall, this study provides valuable insights into tree migration processes, highlighting the varied predictability of treeline responses across regions. These findings carry significant implications for refining vegetation models and guiding conservation strategies to sustain alpine tundra resilience in the face of accelerating climate change.

1 Introduction

1.1 Importance of treeline research

Treelines across the Northern Hemisphere are shifting upward and northward in response to global warming (Harsch et al., 2009; Holtmeier and Broll, 2010; Lu et al., 2021). This trend is particularly evident in alpine treelines within boreal forests, as climate change is progressing more rapidly at high elevations and latitudes than elsewhere (IPCC, 2022). Several studies report elevational treeline advancement in the mountainous regions of Eurasia, such as the Altai Mountains (Cazzolla Gatti et al., 2019) and Polar Urals (Devi et al., 2008, 2020), as well as in North America, including the Kenai Mountains (Dial et al., 2007) and Yukon Territory (Danby and Hik, 2007).

The upward expansion of trees into higher elevations intensifies competition for light, water, nutrients, and habitat, threatening endemic alpine species that depend on open tundra or meadow ecosystems. This increased competition can limit the growth and survival of native alpine plants adapted to specific environmental conditions now altered by tree encroachment (Greenwood and Jump, 2014). As a result, treeline advancement may drive specialised flora and fauna to extinction if they cannot adapt or migrate quickly enough (He et al., 2023), leading to biodiversity loss and disruption of ecological relationships (Hansson et al., 2021). Current studies reveal substantial regional variation how both elevational and northern treelines respond to environmental changes (Holtmeier and Broll, 2010). Despite ongoing research, the specific drivers and, especially, the constraints shaping alpine treeline dynamics in boreal forests under climate change are not yet fully understood (He et al., 2023).

1.2 Time lag in treeline response to climate warming

Numerous favourable and unfavourable factors influence tree establishment within the treeline ecotone and, consequently, treeline advance. First and foremost, climatic factors determine the position of the treeline. During his expeditions in the late 18th and early 19th centuries, Alexander von Humboldt observed a pronounced correlation between latitudinal and elevational treelines and corresponding temperatures. Humboldt integrated this revelation into his pioneering concept of isotherms, lines connecting elevations with equivalent temperatures. This framework enabled him to systematically connect mountains across the globe based on their treelines, revolutionising our understanding of climate patterns and their spatial distribution on a global scale. Körner (1998) later formulated the “growth limitation hypothesis” (GLH), which explains cool air temperatures during the growing season as the primary factor determining the upper limits of treelines. Contemporary observations supporting the GLH show that the global treeline aligns with a Humboldtian isotherm, where mean growing season air temperatures (GSAT) are around 6.7 °C (Körner and Paulsen, 2004) and 6.4 °C (Körner, 2012; Paulsen and Körner, 2014). Treelines warmer than this mean GSAT 6–7 °C isotherm (Millar et al., 2020) indicate that factors beyond cold temperatures also significantly influence the treeline position (Maher et al., 2021). Multiple studies highlight a notable time lag between climatic shifts and the corresponding latitudinal or elevational treeline response (Holtmeier, 1985; Hofgaard and Wilmann, 2002; Lloyd, 2005; Zheng et al., 2024). The fact that only a little more than half of global treelines have advanced in response to anthropogenic warming (Harsch et al., 2009) further underlines the need to integrate factors beyond GSAT into vegetation models.

Vegetation models are crucial for predicting tree responses to changing climate. Numerous simulation studies have employed dynamic global vegetation models (DGVMs), such as SEIB-DGVM (Spatially Explicit Individual Based DGVM) (Smith et al., 2001) or LPJ-GUESS (Lund-Potsdam-Jena General Ecosystem Simulator) (Sato et al., 2007). This study employs LAVESI (Larix Vegetation Simulator) (Kruse et al., 2016), an individual-based, spatially explicit model developed to simulate latitudinal migration of larch treelines in northern Siberia.

Studies using LAVESI, which incorporates many variables influencing treeline migration, also indicate that vegetation models often overestimate the transition from treeless tundra to forest (Kruse et al., 2016). While temperature is widely recognised as the primary driver of tree growth in this region, dispersal constraints may significantly hinder the northward migration of treelines. These findings underscore the pronounced time lag in tree population responses to climate change. However, despite the inclusion of many relevant factors in LAVESI, essential fine-scale processes critical to accurately simulating vegetation dynamics may still be underrepresented.

1.3 Factors influencing alpine treeline migration

The establishment and upward expansion of trees are governed by a complex interplay of abiotic climatic influences, species-specific traits, abiotic non-climatic constraints, and biotic interactions (Fig. 1).

Figure 1 Overview graphic for the key factors influencing alpine treeline migration.

1.3.1 Abiotic climatic factors

Treeline migration is strongly linked to rising summer temperatures (Almqvist et al., 1998; Holtmeier and Broll, 2007; Kullman, 2007; Devi et al., 2008, 2020; Grigoriev et al., 2022) and prolonged growing seasons, where an increase in warm days enhances photosynthetic performance (Holtmeier and Broll, 2007). In addition to GSAT, winter temperatures also play a critical role, as mild winters reduce seedling mortality, whereas early and late frosts hinder establishment (Kullman, 2007; Holtmeier and Broll, 2007). Frequent freeze-thaw cycles pose a greater risk than prolonged low temperatures (Gross et al., 1991). Elevated atmospheric CO₂ levels, however, do not significantly enhance tree growth, as carbon availability is not a limiting factor (Körner, 2003; Wieser et al., 2019).

Soil moisture is a crucial determinant, as increased summer precipitation enhances seedling emergence and survival (Devi et al., 2020), whereas drought imposes a significant constraint (Bailey et al., 2021). Otherwise, both paludification and waterlogging also impede treeline advance (Crawford, 2003; Holtmeier and Broll, 2007). Winter precipitation also has a major influence on the elevational treeline. Low snow fall, especially in wind-exposed sites, can expose seedlings and saplings to climatic stress and herbivory, hindering their survival (Gross et al., 1991; Holtmeier, 2003, 2005a, b; Kullman, 2005; Holtmeier and Broll, 2007). In contrast, snow cover provides essential protection, promoting forest expansion in snow-rich and consequently, sheltered sites (Holtmeier and Broll, 2007; Kirdyanov et al., 2012; Hagedorn et al., 2014). Increased winter snow depth enhances seedling establishment, aiding treeline advance (Germino and Smith, 2000; Smith et al., 2003; Holtmeier, 2005a; Grigoriev et al., 2022). However, persistent snow cover can shorten the growing season and delay tree establishment (Holtmeier and Broll, 2007), although the delayed onset may coincide with warmer temperatures and increased moisture availability (Walker et al., 1999). Excessive snowfall can also impede tree growth through destructive snow movements such as avalanches, snow slides, and snow creep, damaging crowns and branches (Holtmeier and Broll, 2007; Autio, 2006; Holtmeier, 2005a). Furthermore, wet snow exacerbates damage and promotes the spread of snow fungi, increasing seedling and sapling mortality (Holtmeier, 2005a, b; Holtmeier and Broll, 2007; Barbeito et al., 2013).

Wind can have very different effects on the advance of the treeline. Low to moderate wind speeds promote the dispersal of seeds and pollen, but high wind speeds cause physiological stress and frost-drying (Holtmeier, 1985; Holtmeier and Broll, 2007; Kruse et al., 2018). Seed accumulation occurs primarily in wind-protected shelters (Kullman, 2005; Resler et al., 2005; Anschlag, 2006; Holtmeier and Broll, 2007; Beloiu et al., 2022).

1.3.2 Abiotic non-climatic factors

Topography, geomorphological processes, and local site conditions, including microclimate and soil characteristics such as temperature, moisture, and mineralisation, can exert a greater influence on treeline advancement than rising temperatures (Holtmeier and Broll, 2005, 2007). In steep-sided high mountain valleys, treeline elevation is largely shaped by regional orography, with factors like avalanches, mass wasting, unstable slope debris, and fragmented or absent soil cover preventing forests from reaching their thermal limit (Holtmeier and Broll, 2005, 2007). Permafrost thawing further constrains treeline advance by destabilising steep slopes, and intensifying erosion in potential forest habitats (Lloyd, 2005; Holtmeier and Broll, 2007). Wildfires can generate favourable seedbeds, but also disrupt treelines through altered soil temperatures and snow regimes (Brehaut and Brown, 2022).

1.3.3 Species-dependent factors

As long as there are no disturbances, the height of the treeline is climatically determined (Körner, 2020); though recruitment-related factors influence the local appearance and demography of the treeline (Wieser et al., 2019). The success of development from seed to seedling to sapling is contingent upon various factors including the availability of mature, germinable, and viable seeds, seed dispersal and dispersal method, the presence of suitable seedbeds, and the rates of establishment, growth, and survival of seedlings and saplings (Holtmeier, 1995; Juntunen and Neuvonen, 2006; Holtmeier and Broll, 2007, 2019). In cold climates, short vegetation thrives in the relatively warm microenvironment near the ground, while trees reach maturity only if they withstand the surrounding air temperature (Wilson et al., 1987; Grace, 1989). Warm summers enhance seed viability, seedling establishment, and growth to a size that ensures winter survival (Karlsson and Weih, 2001; Holtmeier and Broll, 2007). In arid continental climates, soil moisture is the dominant limiting factor for the development of tree seedlings and saplings (Holtmeier and Broll, 2005, 2007). Trait variability fosters diverse variants and adaptability to environmental changes. Inheritance drives adaptation, promoting successful traits, enhancing population spread and resilience (Holtmeier and Broll, 2007; Gloy et al., 2023). Overaged trees exhibit reduced reproductive capacity, limiting treeline advance (Holtmeier and Broll, 2007). Global and regional studies often overemphasise broad drivers such as temperature, overlooking other crucial variables (Holtmeier and Broll, 2019).

1.3.4 Other biotic factors

During early treeline expansion, open tree stands accumulate more snow than treeless areas (Holtmeier, 2005a). As tree density increases, a warmer microclimate develops, mitigating summer frost damage by shielding seedlings and saplings from intense solar radiation after cold nights and strong winds (Örlander, 1993; Germino and Smith, 2000; Germino et al., 2002; Smith et al., 2003; Johnson et al., 2004; Holtmeier 2005a; Holtmeier and Broll, 2007; Sigdel et al., 2020). However, higher tree density intensifies both intra- and inter-specific competition, including above- and below-ground resource competition, which may ultimately limit further treeline migration (Kruse et al., 2016; Zheng et al., 2024). In addition to competition from trees, non-arboreal vegetation, such as grassland and dwarf shrub vegetation, further restricts seedling access to nutrients, moisture, and light (Holtmeier, 1995; Löffler et al., 2004; Dullinger et al., 2003, 2004; Holtmeier and Broll, 2007; Liang et al., 2016). In contrast, bare mineral soils favour seedling establishment by reducing competition (Hobbie and Chapin, 1998; Löffler et al., 2004). Therefore, moderate grazing may actually facilitate tree growth by reducing competing vegetation (Holtmeier, 1995; Löffler et al., 2004; Holtmeier and Broll, 2007).

Beyond plant competition, treeline dynamics are also shaped by natural and anthropogenic disturbances, including mass outbreaks of leaf-eating insects like the autumnal moth and grazing by reindeer and other pasturing animals, both of which affect seedling survival and forest expansion (Holtmeier et al., 2004; Holtmeier and Broll, 2007; Wieser et al., 2019; Hansson et al., 2021; Brehaut et al., 2022).

To enhance our understanding of treeline dynamics and to improve the predictive capabilities of vegetation models under climate change, this study examines the relative impact of various factors on alpine treeline migration. Specifically, we extended the LAVESI vegetation model by integrating a novel snow module that enables us to assess the influence of varying snow depths on treeline advancement.

This leads us to investigate two key research questions:

1. What are the most significant factors influencing the migration rates of alpine treelines, and how do these factors compare in terms of their respective impacts?

2. Does the rate of upward migration of alpine treelines within the alpine treeline ecotone vary in response to changes in snow regime?

2 Methods

2.1 Model description

In our study, we used the individual-based and spatially explicit model LAVESI (Larix Vegetation Simulator) (Kruse et al., 2016). LAVESI, initially developed for Larix gmelinii in northeastern Siberia, was designed with the aim of simulating the life cycle of larch species as completely as possible. Following pattern-oriented modelling, a simple model was iteratively refined to more accurately reproduce observed patterns (Grimm and Railsback, 2005; Grimm et al., 2005). Within one simulation step, which is one year, the relevant processes are incorporated consecutively as submodules. The previous version of LAVESI already incorporates many of the forementioned factors influencing treeline migration rates. At the beginning of each annual cycle, an “Environment update” takes place. On the one hand, climate variables and a weather index, namely monthly means of temperature and monthly precipitation sums, are processed to estimate daily values. On the other hand, competition is assessed using an annually updated map of tree population density. LAVESI then calculates an annual cycle of seed production, seed dispersal, establishment and competition-dependent tree growth. Finally, ageing and mortality complete a LAVESI simulation step. As part of the validation process, model outputs are compared with actual forest characteristics. While treelines are fitted, this is not done in a spatially explicit manner. LAVESI has been extended to include dispersal processes related to wind speed and direction (Kruse et al., 2018), landscape topography and further boreal forest species (Kruse et al., 2022), adaptability to environmental changes by implementing trait variation and inheritance (Gloy et al., 2023), and forest fire (Glückler et al., 2024). For North America, disturbance effects caused by insects were also included in the model.

2.2 Description of impact factors

In the initialisation phase of each simulation step, which is one year, the monthly weather data are processed by LAVESI (Kruse et al., 2016). The length of the growing season is calculated from the sum of days with a temperature greater than 0 °C, called “net degree days” (NDD0), for each simulation year. The summer temperature is based on the sum of the temperatures of the days above 10 °C, the so-called “active air temperature” (AAT10). The average January temperature is used as an indicator for the winter temperature. Summer precipitation is the mean precipitation for the months of May to August; winter precipitation – e.g. temperature-dependent snowfall and accumulation – is part of the new snow module. A wind module, where seed and pollen dispersal are controlled by wind speed and direction, has been implemented in LAVESI-Wind (Kruse et al., 2018). Increased tree mortality due to wind exposure is implemented in the snow module of this new version of LAVESI. In LAVESI, the mortality rate is determined by a combination of abiotic and individual factors; once mortality occurs, the tree is removed from the simulated plot (Kruse et al., 2016). Species-dependent factors of seed production, age of maturity, number of seeds introduced, seedling establishment, and seedling mortality are already included as sub-modules in the initial model (Kruse et al., 2016). Wind-dependent seed dispersal distance, which is positively correlated with the height of the releasing tree, as well as long-distance dispersal, are part of the wind module in LAVESI-Wind (Kruse et al., 2018). The impact of overageing is considered as a specific mortality rate, which is included as a sub-module in the initial version of the model (Kruse et al., 2016). As part of the implementation of topography and landscape sensing of individuals, the slope and the topographic wetness index (TWI) derived from a Digital Elevation Model (DEM) were first implemented in Kruse et al. (2022). The availability of the seedbed is derived from the litter layer height, the dynamics of which were first added to LAVESI in Kruse et al. (2022). The estimation of the maximum active layer thickness, which is strongly coupled to air temperature, was introduced in the original setup of LAVESI (Zhang et al., 2005; Kruse et al., 2016). Tree density dependent mortality due to competition was also included in the original version of LAVESI (Kruse et al., 2016), whereas tree density dependent facilitation is part of the new snow module. Competition with non-arboreal vegetation such as grassland and dwarf shrub vegetation is only indirectly implemented in the seedbed availability factor.

2.3 New model development: The snow module

In this study, we examine the driving forces and, in particular, the constraints that influence the responses of alpine treelines. The current version of the LAVESI model successfully simulates the dynamics of the Siberian and North American latitudinal treelines. Since many factors influencing treeline migration are common to both latitudinal and elevational treelines, we adapted LAVESI to account for the specific characteristics of alpine treelines. As part of this adaptation, various effects of snow on treeline migration are incorporated into the model to better capture the unique dynamics of alpine environments.

In recent decades, there have been a number of efforts to incorporate snow dynamics into individual-based gap models. In the SEIB-DGVM (Spatially Explicit Individual-Based DGVM) model, Sato et al. (2007) introduce a method for distinguishing between rainfall and snowfall based on an empirical function of daily mean air temperature. This approach allows for the accumulation of snowfall in a designated snow pool, which subsequently melts in response to soil temperature variations. UVAFME (University of Virginia Forest Model Enhanced) (Wang et al., 2017) has been significantly updated from its original version FAREAST (forest gap model to simulate dynamics and patterns of eastern Eurasian forests) (Xiaodong and Shugart, 2005). However, neither individual-based gap models account for snow accumulation or melt processes. The revised version of UVAFME by Foster et al. (2017) addresses this gap by implementing a simplified snowmelt submodel using the degree-day method proposed by Singh et al. (2000). In their framework, if daily air temperatures drop below 5 °C, precipitation is classified as snowfall and added to the snow cover; conversely, temperatures above 0 °C trigger the melting process for that day.

In previous versions of the individual-based, spatially explicit vegetation model LAVESI (Kruse et al., 2016, 2018, 2022; Gloy et al., 2023; Glückler et al., 2024), precipitation was limited to the growing season months of May to August, neglecting any contributions from precipitation outside this period. Consequently, processes related to snow accumulation and melt were not incorporated into the model. As part of the adaptation of LAVESI for elevational treelines, a range of snow-influenced processes were incorporated into the model, addressing both positive and negative ecological effects across varying snow conditions. Snow depth is calculated from precipitation using a snow-to-water ratio, which can vary significantly depending on factors such as the amount of ice in the snow-producing clouds, the structure of snowflakes, and ambient temperature. As a general guideline, we apply an average snow-to-water ratio of 10:1 (National Weather Service, 2024). Following the approach of Sato et al. (2007), which builds on the precipitation partitioning method of Ito and Oikawa (2002), processes of snow accumulation and snowmelt were integrated into the LAVESI model. Ito and Oikawa (2002) applied a gradual transition in which 50 % of total precipitation is considered snow at 2 °C. We initially adopted this approach, but subsequently refined the temperature threshold during model tuning. By comparing simulated and observed snow accumulation patterns, the threshold was iteratively adjusted to improve model accuracy. Unlike models such as SEIB-DGVM and UVAFME, which primarily focus on temperature and moisture dynamics without fully integrating snow effects, the revised version of LAVESI explicitly considers snow depths variability and its impact on treeline dynamics. Snow-related processes are governed by the interplay between snow depth, timing of snowmelt, tree height, and stand structure, all of which influence tree establishment, growth, and mortality (Fig. 2).

Figure 2 Flowchart illustrating snow-related processes influencing alpine treeline migration, as implemented in the snow module of the vegetation model LAVESI (Larix Vegetation Simulator).

Tree mortality is influenced by the interaction between tree height and snow depth. If the snow depth is lower than the tree's height, we assume that the mortality rate of the tree increases by about 50 % due to climatic stress. This value is a simplified estimate. We further assume that this effect decreases linearly with increasing tree height and that it disappears once a tree reaches approximately 5 m in height.

In contrast, if a seedling or young sapling is assumed to be completely covered by snow, we estimate that it benefits from insulation, leading to about a 30 % increase in growth for that year. This growth bonus is again an assumed approximation.

We also estimate that snow cover positively influences germination. Specifically, we assume a linear increase in germination probability with snow depth, up to a maximum improvement of +30 % at 1 m snow depth. Additionally, we include a simplified function for the timing of snowmelt: if snow melts before day 110 of the year, we estimate an increase in germination probability; if melting occurs later, germination probability is assumed to decrease. This relationship is defined in the model as a linear adjustment relative to the threshold day (DOY_threshold = 110).

Finally, we incorporate extreme snowfall and avalanche events as stochastic disturbances. In our assumptions, if daily snowfall exceeds the local 99th percentile, trees affected by this event are assigned a 50 % increase in mortality. This figure is again to be understood as an estimated penalty representing increased stress and damage under extreme conditions.

Avalanche risk is determined from a maximum snow depth during a winter season. A dynamic threshold is derived from the top 0.1 % of snow depth values across the modelled climate time series. When an avalanche occurs, it induces mortality primary at the leading edge of tree stands, with damage decreasing downslope due to the protective effect of surrounding vegetation. Trees affected by an avalanche exhibit reduced seed production for up to 10 years, depending on the severity of the impact.

In the updated version of the vegetation model, we incorporated new processes to account for how tree density modulates both competitive and facilitative interactions under varying snow conditions (Zheng et al., 2024). Tree mortality due to competition – i.e. the negative effect of increasing tree density on individual survival – was already implemented in previous model versions using competitive strength evaluations (Kruse et al., 2016). As tree density increases, competition for light, water, and nutrients intensifies, resulting in higher mortality rates.

However, in this study, we also account for facilitative effects, which may mitigate mortality under certain conditions. Specifically, increased tree density can improve the microclimate and physical environment – for instance by providing wind shelter or promoting snow accumulation – which in turn may reduce exposure-related mortality. Facilitation strength is modelled as a function on tree density: it increases linearly with tree density up to a value of 0.5 (dimensionless units), where it reaches its maximum. Within this range, the facilitative effect linearly reduces tree density-dependent mortality. Beyond a tree density of 0.5, the facilitative effect gradually diminishes and vanishes completely when tree density reaches 1.0 At densities exceeding 1.0, only competitive effects remain active in the model, and facilitation no longer occurs. This implementation allows the model to represent a realistic balance between competition and facilitation, particularly under snow-affected conditions. It captures how sparse stands may suffer from exposure (e.g. to wind and snow abrasion), while moderately dense stands benefit from protective interactions, and overly dense stands experience intensified competition.

Enhanced winter snow fall also affects soil moisture availability during the growing season. This feedback is incorporated into the general weather processing during the initialisation phase, modifying the species-specific weather index. As a result, the moisture-related climate-growth response of individual tree species is adjusted accordingly.

The LAVESI-Mountain Treelines 1.0 model code, along with the input data used in this study, is available on GitHub (https://github.com/StefanKruse/LAVESI/tree/circumboreal_mountain, Kruse et al., 2025, last access: 12 September 2025). The code, input data, and simulation output will be permanently archived on Zenodo (Haupt et al., 2025).

Table 1 Characteristics of the study sites.

Table 2 Adjustments to LAVESI simulation settings across study sites.

Figure 3 Map of the study sites: Fort McPherson West, Canada (CA); Lake Ilirney, Russia (RU); and Road to Central, Alaska, USA (AK).

Figure 4 Photographs of the study sites: (a) Fort McPherson West, Canada (CA); (b) Lake Ilirney, Russia (RU); and (c) Road to Central, Alaska, USA (AK).

2.4 Study sites

The model had to be parameterised for the three study sites in Fort McPherson West, Canada (CA), at Lake Ilirney, Russia (RU), and at the Road to Central, Alaska, USA (AK) (see Table 1; Figs. 3–4), for which the respective adjustments were made in the input files parameter.txt and declaration.h (see Table 2). Satellite images of the study sites, along with a digital elevation model (DEM), slope, and topographic wetness index (TWI) for the simulation area, are provided in Fig. 5.

Figure 5 Satellite images of the study sites: (a) Fort McPherson, Canada (CA), (b) Lake Ilirney, Russia (RU), and (c) Road to Central, Alaska/ USA (AK), including the simulation area (red box). Digital elevation model (DEM), slope, and topographic wetness index (TWI) for the simulation area. Esri, i-cubed, USDA, USGS, AEX, Geoeye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community.

2.5 Model inputs and simulation scenarios

Monthly precipitation and temperature data for the historical period (1400–1950 CE) were obtained from the Max Planck Institute Earth System Model (Transient MPI-ESM Glac1D_P3) (Kapsch et al., 2022). This past climate data was bias-corrected to align with modern data from the CRU-TS dataset (Harris et al., 2020; available for 1900–2020 CE), ensuring consistent means during the overlap period (1900–1950 CE) to avoid discontinuities. The corrected dataset was extended to 2020 CE using CRU-TS. The first 620 years are both to simulate the current conditions of the populations and as the spin up-phase to allow the forests to establish from seeds introduced over the entire area. Climate projections for 2020–2100 CE were derived from the CMIP6 MPI-ESM-2-LR model (Eyring et al., 2016; specific version: r1i1p1f1) and similarly bias-corrected to CRU-TS, using the overlap period (2015–2020 CE) to ensure consistency.

Monthly wind data for 1400–1950 CE was obtained from the transient MPI-ESM Glac1D_P3 model. Since the vegetation model LAVESI required 6-hourly wind data, a downscaling approach was applied. Modern wind data from the ERA5 dataset (Hersbach et al., 2023), covering 1940–2020 CE, was used to create a 6-hourly resolution dataset for the overlapping period and extend the dataset to 2020 CE. ERA5's detailed temporal resolution facilitated realistic downscaling and interpolation. Future wind data from the CMIP6 MPI-ESM-2-LR model (2020–2100 CE), provided at daily resolution, was interpolated to 6-hourly means to meet LAVESI's requirements.

The environmental input data for LAVESI was derived from the Ensemble Digital Terrain Model (EDTM) with a 30 m resolution (Ho et al., 2023). Key information, including elevation, slope, and the topographic wetness index (TWI), was extracted from the DEM of the simulation area.

Simulation runs are performed using a transect representing a slice of the alpine treeline with a size of 1 km × 0.1 km.

To evaluate and compare the influence of various factors on the migration rate of the alpine treeline and to identify key drivers and assess how variations in these factors influence treeline dynamics, a sensitivity analysis was performed. Most of the factors examined in this study have already been implemented in the original LAVESI model (Kruse et al., 2016) or its subsequent versions. These factors are detailed extensively in the corresponding publications. A comprehensive list of species traits and model parameter values used in this study is provided in the Supplement. The only newly introduced component in this study was the snow module, which is specifically designed to simulate the impact of snow cover dynamics on the treeline. A description of this snow module is provided earlier in this manuscript. For the sensitivity analysis following Kruse et al. (2018), simulations were run in which the following 17 factors were increased by 10 % and decreased by 10 % compared to an average reference value: growing season length, summer temperature, winter temperature, summer precipitation, wind speed, seed production, maturation age, seed intro number permanent, dispersal distance, long distance dispersal, seedling establishment, seedling mortality, overageing mortality, slope, twi, seed bed availability, and active layer (see Fig. 1).

In addition, further simulations were run with the forest fire module (Glückler et al., 2024), the insect disturbance module, or the module for trait variation and inheritance (Gloy et al., 2023). The results of each of the simulations were compared with the results of the reference run.

Finally, simulations were conducted using the new snow module, with a 10 % increase and decrease applied to the following factors within the module: winter precipitation, wind exposure, and facilitation. The results were then compared to those obtained from the reference snow module run.

We compared the effect of all these settings on treeline migration rate, where the treeline is defined as the northernmost position where forest density does not fall below a threshold of one hundred trees per hectare, where a tree must be >1.3 m tall.

In total, 45 different simulation scenarios were computed (for a structured overview, see Table A1 in the Appendix). These include those evaluating model sensitivity to predicting variables, as described before, as well as the simulations with and without the forest fire module, the insect disturbance module, the module for trait variation and inheritance, and the new snow module including the evaluation of the model sensitivity to predicting variables within the snow module.

Each of the simulation scenarios was repeated 40 times.

2.6 Statistical analyses of simulation

The output files generated by the LAVESI model provide information on stem count and tree density for each 5×5 m grid cell within the simulation area. The statistical analysis was performed using the statistical tool R (R Core Team, 2023).

Stem counts for all tree species were aggregated across spatial positions for every fifth year between 2020 and 2095. The aggregated data are visualised as line graphs using the “ggplot2” R package (Wickham, 2016) to illustrate changes over time across different scenarios. Additionally, spatial distributions of tree densities are depicted with separate raster plots created for the years 2020, 2055, and 2095 CE. Treelines were determined using threshold-based criteria of one hundred trees per hectare.

The migration rate of the alpine treeline was recorded for each simulation by calculating the change in treeline position between 2020 and 2095 for each scenario. Using the diff function in R, we measured incremental changes between consecutive treeline positions, capturing the position change between each measured interval. These differences were then averaged to obtain the mean migration rate, which represents the average advancement of the treeline over time. The significance of the migration rates was assessed using a paired t test, conducted with the t.test function in R.

The relative change in migration rate was calculated for both the increased and decreased scenarios to quantify the sensitivity of the treeline to each factor. All resulting migration rates were computed to sensitivity values that present the percentage of change in migration rate scaled to the 10 % change in the factor. The sensitivity values are calculated from the following Eq. (1):

$$\begin{array}{}\text{(1)}& S={\displaystyle \frac{\mathrm{\Delta }Y\times \mathrm{100}}{\mathrm{\Delta }P}},\mathit{}\end{array}$$

where ΔY represents the change in the migration rate derived from each simulation run and ΔP is the 10 % change in the parameter. The significance of the sensitivity values was evaluated using a one-sample t test, performed with the t.test function in R (R Core Team, 2023).

We conducted a redundancy analysis (RDA) using the rda function from the “vegan” R package (Oksanen et al., 2024) to assess the relationships between environmental variables and treeline migration rates across the three study sites (CA, RU, and AK). The response variable consisted of treeline migration rates from 45 different simulation scenarios. In these scenarios, key factors influencing treeline migration were manipulated by switching them on or off, or by increasing or decreasing their magnitude compared to an average reference value. Each scenario was simulated across 40 runs per study site. Environmental predictor variables were derived from weather data for each study site, spanning the period from 1990 to 2020. The climate variables, including mean annual temperature, mean January temperature, mean July temperature, mean summer (JJA) temperature, accumulated annual temperature (AAT), degree days above threshold (DDT), total annual precipitation, and summer (JJA) precipitation, were calculated and averaged for each study site. Stepwise forward selection, using the ordistep function from the “vegan” R package, was performed to identify the most significant climate predictors of treeline migration rates. The process began with a null model and progressively added variables based on the Akaike Information Criterion (AIC) and F-statistics. Significant variables (p<0.05) were retained in the final model, with significance evaluated using 999 permutations. For the RDA plot, site (migration rates), species (scenario), and environmental (weather data) scores were extracted and visualised using the “ggplot2” R package.

3 Results

3.1 Comparative analysis of factors influencing alpine treeline migration

To quantify and compare the influence of key predicting variables on simulated treeline migration rates, a sensitivity analysis was conducted across all three study sites, based on ±10 % variations in individual variables. On average, a 10 % increase in predicting variables results in a +2.98 % change in migration rate (SD = 8.76) at site CA, +1.89 % (SD = 7.14) at site RU, and +5.22 % (SD = 32.45) at site AK. A 10 % decrease leads to changes of −2.30 % (SD = 8.43) at site CA, −0.97 % (SD = 6.39) at site RU, and −7.21 % (SD = 39.51) at site AK. The minimum and maximum significant sensitivities observed are +0.43 % and +4.50 % for site CA and +0.50 % and +2.90 % for site RU. For site AK, no significant values were detected (Table 3).

Table 3 Sensitivity analysis results for treeline migration rates. Mean sensitivity values (in % %⁻¹) and standard deviations for 17 predicting variables, each adjusted by ±10 % from an average reference value, resulting in 34 distinct scenarios across three different study sites.

* Highly significant with p<0.05, and the remaining values non-significantly different from the reference run. Values are the mean of 40 simulations.

Site AK exhibits the highest mean sensitivity values; however, these results are characterised by substantial variability, preventing any significant findings. In contrast, site RU generally shows lower sensitivity magnitudes across all scenarios, although some scenarios do reach significance. Site CA displays several scenarios with significant positive or negative sensitivities, indicating a stronger influence on migration dynamics compared to the other sites.

The direction of sensitivity across all study sites is consistent for the majority of predicting variables, with positive and negative sensitivities aligning similarly between sites. There are four exceptions where the direction of influence differs between study sites, but these deviations are minimal and not statistically significant.

Notably, four scenarios demonstrate significant sensitivity at both sites CA and RU (p<0.05). Increased seedbed availability has the most pronounced positive influence on migration, with mean sensitivities of 4.32 % at site CA and 2.51 % at site RU (seedbedavailability_plus10%). In contrast, reduced maturation age strongly hinders migration potential, as reflected by negative mean sensitivities of −3.84 % at site CA and −3.10 % at site RU (maturationage_minus10%). Favourable conditions for seedling establishment also plays an important role, promoting migration with mean sensitivities of 3.39 % at site CA and 2.83 % at site RU (seedlingestablishment_plus10%). Lastly, higher wind speeds facilitate migration through enhanced seed dispersal, as indicated by mean sensitivities of 2.99 % at site CA and 2.61 % at site RU (windspeed_plus10%).

3.2 Impact of snow-related processes on alpine treeline migration rates

The comparison of the results from additional simulations – which included the forest fire module (Glückler et al., 2024), the insect disturbance module, the trait variation and inheritance module (Gloy et al., 2023), and the snow module implemented in this study – with those of the reference run revealed that, across all three study sites, only the simulations incorporating the influence of snow-related factors resulted in a significant change in the treeline migration rate (Fig. 6).

Figure 6 Box plots showing migration rates from 40 simulation runs at three study sites: (a) Fort McPherson, Canada (CA), (b) Lake Ilirney, Russia (RU), and (c) Road to Central, Alaska, USA (AK). Each panel compares the reference scenario with four scenarios in which specific environmental modules were activated. Note that for the Russian study site, the insect disturbance module was not implemented in LAVESI.

The sensitivity of treeline migration rates to ±10 % changes in snow-related predicting variables, relative to the snow module reference run, was analysed across all three study sites. On average, a 10 % increase in predicting variables leads to a +4.35 % change in migration rate (SD = 12.44) at site CA, +1.30 % (SD = 6.87) at site RU, and +4.15 % (SD = 13.28) at site AK. Conversely, a 10 % decrease in predicting variables results in a −6.49 % change (SD = 14.47) at site CA, −3.51 % (SD = 7.86) at site RU, and −2.66 % (SD = 12.31) at site AK. The minimum and maximum significant sensitivities observed are +3.91 % and +6.53 % for site CA, +2.99 % and +4.15 % for site RU, and +4.79 % and +5.30 % for site AK (Table 4).

Table 4 Sensitivity analysis results for treeline migration rates within snow computation simulations. Mean sensitivity values (in % %⁻¹) and standard deviations for three predicting variables, each adjusted by ± 10 % from an average reference value, resulting in six distinct snow computation scenarios across three different study sites.

* Highly significant with p<0.05, and the remaining values non-significantly different from the reference run. Values are the mean of 40 simulations.

At site CA, several scenarios demonstrate strong positive or negative sensitivities, with nearly all results reaching statistical significance (p<0.05). In contrast, site RU exhibits more moderate sensitivities, although some scenarios still achieve statistical significance. Despite relatively high mean sensitivity values at site AK, substantial variability in the results causes fewer statistically significant outcomes.

The sensitivity analysis for snow computation reveals consistent directional influences across the three study sites, with one notable exception: windexposure_plus10%. While sites CA and AK exhibit positive sensitivity values, site RU displays a negative sensitivity value. The sensitivity values at sites RU and AK are statistically significant (p<0.05).

The scenario winterprecipitation_plus10% has statistically significant impacts across all sites, suggesting that increased snow accumulation during winter enhances migration potential. Both facilitation_plus10% and windexposure_minus10% show significant sensitivity at sites CA and RU (p<0.05). Increased facilitation positively influences migration potential, with the strongest effects observed at site CA, while reduced wind exposure negatively impacts migration at both sites. Although site AK displays similar patterns, the variability in the data prevented statistical significance.

Across all three study sites, simulations without the snow module (Snow computation OFF) produce different tree density values over time compared to those incorporating snow-related processes (Snow computation ON). For all sites, the inclusion of snow processes results in an increase in tree density by the years 2055 and 2095 (Fig. 7).

Figure 7 Tree density over time. Tree density in 2020, 2055, and 2095 CE at study site Fort McPherson West, Canada (CA) under (a) the reference scenario (run 34) and (b) the scenario with activated snow computation (run 2). Tree density in 2020, 2055, and 2095 CE at study site Lake Ilirney, Russia (RU) under (c) the reference scenario (run 13) and (d) the scenario with activated snow computation (run 3). Tree density in 2020, 2055, and 2095 CE at study site Road to Central, Alaska, USA (AK) under (e) the reference scenario (run 15) and (f) the scenario with activated snow computation (run 3). (Run with the median migration rate calculated from 40 runs for each simulation scenario). Each plot represents data from the simulation run with the median migration rate derived from 40 simulation runs per scenario. Tree density values (dimensionless units) were log-transformed to enhance data visualisation.

Figure 8 Redundancy Analysis (RDA) biplot showing (a) site scores and (b) species scores along the first two principal axes. Plot (b) provides a close-up of the centre region (grey rectangle) of plot a to highlight tightly grouped species. Site scores (migration rates) are represented by points, colour-coded according to their study site. Species scores (simulation scenarios) are depicted as arrows, colour-coded to indicate whether snow computation is activated (blue: Snow computation ON) or not (red: Snow computation OFF). Black arrows in both plots represent the direction and strength of environmental variables influencing migration rates.

3.3 Relationships between environmental variables and treeline migration rates

The redundancy analysis (RDA) identified two significant constrained axes, collectively explaining 31.32 % of the variance in treeline migration rates across simulation scenarios and study sites (Fig. 8). The site scores reveal distinct clusters of study sites along the RDA1 axis, reflecting geographical differences. Sites CA and RU group closely together on one side of the axis, while site AK is distinctly separate on the opposite side (Fig. 8a). Along the RDA2 axis, the species scores indicate that scenarios incorporating snow-related processes (Snow computation ON) are clearly distinguished from those without snow effects (Snow computation OFF) (Fig. 8b).

The stepwise forward selection in RDA identified total annual precipitation and mean annual temperature as the most significant climate predictors of treeline migration rates across the simulation scenarios and study sites. Total annual precipitation is the first variable included in the model (AIC = 903.76, F=43.80, p=0.001), followed by mean annual temperature (AIC = 898.55, F=7.25, p=0.001). Total annual precipitation shows a positive correlation with treeline migration rates at sites CA and RU, but a negative correlation at site AK. Conversely, mean annual temperature exhibits an opposite effect: positively correlating with migration rates at site AK and negatively correlating at sites CA and RU.

4 Discussion

4.1 Site-specific responses to factors influencing alpine treeline migration

4.1.1 Spatial variability in migration rates

The RDA analysis reveals a clear spatial grouping of locations based on migration rates, with sites CA and RU clustering together, while site AK is distinctly separate. Notably, different tree species dominate at the respective sites (see Fig. B1 in the Appendix). Based on species-specific traits, one would expect Picea glauca to exhibit a faster migration rate in Alaska than Picea mariana in Canada based on faster and taller growth. However, our findings contradict this expectation, indicating that local climatic conditions likely exert a stronger influence on migration dynamics than intrinsic species traits alone. This highlights the importance of unique environmental conditions at each site on migration dynamics. Such geographical and ecological variability in treeline responses aligns with previous findings, reinforcing site-specific factors as the central factor to shaping treeline migration (Harsch et al., 2009; Hansson et al., 2021, 2023).

4.1.2 Sensitivity analysis: site-specific patterns

The sensitivity analysis further uncovered marked differences in responses to ±10 % changes in predicting variables across the three sites (CA, RU, AK). Site AK exhibited high mean sensitivity values but with substantial variability, reflecting the influence of complex and less predictable interactions between ecological and climatic factors. At this site, no single factor emerged as statistically significant, which may be due to stronger interactions among multiple drivers that obscure the effects of individual parameters. Another possible explanation is a greater influence of stochastic processes, further contributing to the uncertainty observed. Although studies explicitly linking stochastic processes to treeline dynamics are currently lacking, Shoemaker et al. (2020) and Lindo et al. (2023) emphasise their broader ecological relevance, particularly with regard to biodiversity patterns and ecosystem functioning. Consequently, migration dynamics at site AK are inherently more uncertain compared to the more consistent and significant sensitivities observed at sites CA and RU. These findings underscore the need to account for localised environmental conditions when predicting treeline responses to climate change.

4.1.3 Consistency and anomalies in sensitivities

Despite these site-specific differences, the overall consistency in the direction of sensitivities across all sites enhances confidence in the key drivers of treeline migration identified by this study.

Table 5 Expected and observed direction of influence of predicting variables on migration rate.

While many observed sensitivities align with ecological expectations, some findings are unexpected and warrant further investigation to clarify these interactions (see Table 5). For instance, the positive relationship between maturation age and migration rate is counterintuitive, as earlier maturation typically shortens generation time and reduces exposure to juvenile mortality before the first reproductive event (Stearns, 2000). However, later maturation offers advantages such as larger tree size at maturity and higher fecundity (Stearns, 2000). Moreover, taller, older trees likely enhance migration potential through improved pollen and seed dispersal distances, consistent with the dispersal model used (Kruse et al., 2018). This extended dispersal likely reduces local competition among seedlings, enhancing migration potential, even when maturation age increases.

Another unexpected finding is the association of higher seedling and overageing mortality with increased migration rates. One explanation could be that, while higher mortality intuitively reduces population growth, it may simultaneously reduce competition and increase resource availability for surviving seedlings, thereby enhancing their establishment and growth (Kruse et al., 2016). This highlights the complex balance between population dynamics and resource competition in determining migration outcomes, as also demonstrated by Zheng et al. (2024).

The positive sensitivity of migration rates to slope angles also stands out, as steeper slopes are typically considered physical barriers to migration. This may be attributed to the creation of diverse microhabitats in complex terrain, which reduce competition and provide unique niches for tree establishment, as observed by Soliveres et al. (2011).

4.1.4 Important factors for treeline migration

The significant sensitivity observed at sites CA and RU underscores the importance of certain environmental factors in promoting treeline migration. For example, reducing maturation age negatively impacts migration potential, likely because of the reduced tree height at maturity, which limits pollen and seed dispersal distances as discussed earlier (Stearns, 2000; Kruse et al., 2018).

Seedbed availability emerges as a crucial factor for migration potential, with sensitivity analysis confirming that increases in suitable seedbeds significantly enhance the migration rate, as expected. Suitable seedbeds provide favourable conditions for tree establishment, especially bare mineral soils, with their reduced competition from dwarf shrubs or grasses, which offer ideal conditions for wind-mediated tree seeds (Löffler et al., 2004; Hobbie and Chapin, 1998; Holtmeier and Broll, 2007; Liang et al., 2016).

Seedling establishment strongly supports treeline advancement, as higher rates of successful seedling establishment drive treeline progression and maintain population stability (Holtmeier and Broll, 2007), though delays in processes such as seed dispersal, establishment, or growth may slow migration rates, as noted by Holtmeier and Broll (2019).

Wind speed is another important factor influencing treeline migration, as moderate wind velocities enhance seed and pollen dispersal, aligning with findings from Holtmeier and Broll (2007) and Kruse et al. (2018).

These findings align closely with existing literature, highlighting the complex interplay between factors such as maturation age, seedbed availability, seedling establishment, and wind speed in shaping treeline migration processes. The observed patterns are consistent with broader ecological principles related to tree migration under changing environmental conditions. However, the site-specific variations in sensitivity emphasise the critical need to account for localised conditions in predictive simulation models.

Overall, the simulations analysing various factors affecting alpine treeline migration reveal statistically significant influences on migration rates at sites CA and RU, offering valuable insights into the key factors driving migration potential under current conditions. In contrast, the variability observed at site AK underscores the need for further investigation into local factors influencing migration sensitivity. These results collectively reinforce the importance of integrating both broad-scale ecological drivers and site-specific conditions to enhance our understanding of treeline dynamics.

4.2 Climate drivers of alpine treeline migration across study sites

Our findings align with previous research showing that climate-induced treeline shifts are shaped by a complex interplay of warming temperatures and precipitation variability, which jointly influence seedling establishment and growth (Kullman, 2005, 2007; Holtmeier and Broll, 2005, 2007; Körner, 2012). Here, total annual precipitation and mean annual temperature emerge as the most influential drivers of treeline migration accounting for over 30 % of the observed variance in migration rates across the different scenarios. Traditionally, growing season temperature has been regarded as a critical factor driving treeline shifts and forest expansion (Dial et al., 2007; Devi et al., 2008, 2020). However, our findings underscore the importance of considering year-round climate conditions in shaping treeline dynamics, a perspective consistent with studies emphasising the importance of mean annual temperature in influencing treeline migration (Cazzolla Gatti et al., 2019; Danby and Hik, 2007).

In the Polar Urals, winter precipitation has more than doubled over the 20th century (Devi et al., 2020), while in the Kenai Mountains, Alaska, mean annual temperatures have risen by 2 °C since the 1950s, with winter temperatures increasing by as much as 4 °C (Dial et al., 2007). These trends underscore the necessity of considering winter climate, in addition to growing season conditions, when studying treeline dynamics.

Although precipitation outside the growing season is often simply added to its annual sum, primarily in relation to soil water balance, our results highlight the multifaceted role of snow, which extends beyond its contribution to water availability (see Sect. 4.3).

Site-specific responses to climate drivers

The migration responses to the climate drivers temperature and precipitation vary between the sites (see Table 1). Site AK is characterised by relatively higher mean annual temperatures and lower precipitation compared to sites CA and RU.

At site AK, warmer mean annual temperatures positively influence migration. This effect is likely due to an extended growing season and improved soil warmth, which enhance seedling germination and establishment. These findings align with evidence suggesting that warmer conditions promote tree germination, establishment, and survival during the growing season (Almqvist et al., 1998; Holtmeier and Broll, 2007; Kullman, 2007; Devi et al., 2008; 2020; Grigoriev et al., 2022). In contrast, at the cooler, wetter sites CA and RU, higher temperatures are associated with reduced migration rates. This negative effect may be attributed to drought stress, that may override temperature benefits at upper treeline boundaries, as warming-induced soil moisture limitations could hinder tree growth and prevent treeline advance in a warmer world (Bailey et al., 2021; Gruber et al., 2022).

The influence of precipitation also varies across the sites. At sites CA and RU, higher precipitation is positively correlated with accelerated migration rates. This is likely due to improved soil moisture, which supports seedling establishment and early growth in these regions (Holtmeier and Broll, 2007). These findings are consistent with studies showing that treeline migration and forest expansion are linked to increased moisture availability, particularly from higher winter precipitation (Hagedorn et al., 2014; Grigoriev et al., 2022). Conversely, at site AK, higher precipitation is associated with slower migration rates, possibly due to negative factors such as waterlogging or paludification implicitly modelled by an optimum precipitation and adverse effects on potential growth in the model (Kruse et al., 2016). These findings highlight the potential for moisture-related limitations on tree expansion under certain site conditions (Crawford, 2003; Holtmeier and Broll, 2007).

Additionally, the different tree species present across sites likely play a role in shaping these site-specific responses. Species-specific stand structure dynamics influence how various species respond to environmental conditions, further contributing to the observed variation in migration rates (Chhetri et al., 2020). For instance, the deciduous species in Siberia (RU) contrast with the evergreen species in North America (CA and AK), which may affect their climatic sensitivities. Traits such as seed production, growth rate, and the ability to recover from disturbances can modulate species' responses.

4.3 The role of snow in alpine treeline dynamics

4.3.1 Significant impact of snow-related processes on alpine treeline migration

Our study further demonstrates that snow-related processes, incorporated into the LAVESI snow module, consistently reveal significant changes in migration rates across all three study sites. This finding aligns with the work of Huang et al. (2023), who show that seasonal snow-cover patterns – particularly spring snowline elevation and snow-cover duration – affect treeline elevation in the eastern Himalaya. Similarly, Niittynen and Luoto (2018) demonstrate that incorporating snow persistence data into species distribution models enhances predictive accuracy and refines distribution mapping in high-elevation and high-latitude ecosystems. Additionally, Autio (2006) predicts that warmer temperatures could lead to increased snowfall, potentially resulting in heavier snow loads that may damage or even destroy tree crowns and branches. Collectively, these findings underscore the critical role of snow as an environmental predictor shaping treeline dynamics in alpine and subarctic ecosystems.

In contrast, recent studies evaluating environmental modules with the model LAVESI, such as explicit forest fire (Glückler et al., 2024), biotic (insect) disturbances, and trait variation and inheritance (Gloy et al., 2023), have highlighted the significance of various processes in influencing tree migration. However, our simulations do not reveal notable impacts from these factors at alpine treeline ecotones. This discrepancy may be attributed to the relatively short duration of our simulations, which might not have allowed enough time for these processes to fully manifest. For the study sites of our study, findings emphasise the immediate and dominant influence of snow-related processes on treeline dynamics, compared to longer-term processes like fire or trait evolution.

4.3.2 The dual role of snow and model enhancements

The separation between scenarios with and without snow effects along the second RDA axis further underscores the substantial impact of snow-related processes on migration rates. Snow acts as both a facilitator and a hindrance, depending on environmental context. Incorporating snow effects is thus critical for accurately assessing treeline migration in snow-dominated regions (e.g. Niittynen and Luoto, 2018; Huang et al., 2023).

Compared to earlier implementations of snow effects in vegetation models, our new snow module represents a significant advancement. Instead of merely adding precipitation outside the growing season to the annual total, the module captures the complex and multifaceted role of snow, extending beyond its contribution to water availability. Hansson et al. (2023) note that average precipitation changes often show weak correlations with treeline migration, likely due to model limitations that neglect variables like moisture stress and extreme precipitation events. Similarly, Kharuk et al. (2022) observe that under harsh conditions, factors such as wind and soil moisture often outweigh the role of snow. By incorporating previously overlooked factors such as wind exposure and microtopographical features, the enhanced module provides a more nuanced representation of treeline dynamics, capturing the complex interactions between snow and other environmental drivers. These processes – modelled on findings from other studies and existing knowledge – are not explicitly validated with data in this study but are tested to better estimate their influences.

4.3.3 Sensitivity analysis and site-specific variations

The sensitivity analysis conducted with the snow module reveals significant site-specific variations in migration rates in response to ±10 % changes in predicting variables across the three study sites (CA, RU, AK). The higher sensitivity values observed at site CA likely reflect more extreme or variable environmental conditions, whereas the relatively lower sensitivity at site RU suggests a more stable environment where migration processes are less reactive to changes. Site AK exhibits high sensitivity values, but the lack of statistically significant results in some scenarios suggests complex interactions and local factors not fully captured by the model.

The contrasting sensitivity observed in the scenario with a 10 % increase in wind exposure highlights the complex role of wind in tree migration. At sites CA and AK, positive sensitivity values indicate that wind may enhance migration by aiding seed and pollen dispersal, which are crucial drivers of migration potential in these regions (Holtmeier and Broll, 2007; Kruse et al., 2018). Wind-driven dispersal facilitates seed movement over long distances, likely explaining the positive association between wind exposure and migration at these sites (Holtmeier and Broll, 2007). However, the negative sensitivity observed at site RU suggests that wind exposure here may exert physiological stress that inhibits treeline migration. High wind velocities, for example, can damage trees through frost-drying and limit growth, a pattern documented in other wind-exposed treeline areas (Holtmeier, 1985). Moreover, seeds carried by strong winds tend to accumulate in sheltered areas, complicating migration in exposed environments (Anschlag, 2006). Under such conditions, topographic factors such as wind exposure can outweigh potential dispersal benefits and prevent the advance of the treeline despite significant climate warming, as observed in other wind-exposed environments (Kullman, 2005; Beloiu et al., 2022). Additionally, wind erosion of exposed soils can reduce available seedbeds beyond the treeline, further inhibiting migration (Holtmeier et al., 2004; Anschlag, 2006). This duality underscores the complex role of wind in shaping treeline dynamics across varying environments.

Positive sensitivity to increased winter precipitation across all three study sites aligns with established ecological principles emphasising the role of moisture availability (Walker et al., 1999) and protective cover during critical early growth stages (Gross et al., 1991; Holtmeier, 2005a). Snow cover can provide crucial shelter for seedlings and saplings, protecting them from frost damage and herbivory during winter months, thereby supporting seedling survival (Juntunen and Neuvonen, 2006; Holtmeier and Broll, 2007). When snow melts, the resulting increase in soil moisture enhances seedling establishment during the growing season, especially in regions where winter precipitation forms a significant portion of the annual water supply (Germino and Smith, 2000; Smith et al., 2003; Hagedorn et al., 2014). However, excessive snow accumulation can also have negative effects, such as snow creep or slides that may threaten seedling survival (Holtmeier, 2003, 2005a, b; Kullman, 2005). Moreover, wetter snow conditions may promote snow fungi infections, resulting in increased mortality among seedlings of evergreen conifers (Holtmeier, 2005a, b; Holtmeier and Broll, 2007; Barbeito et al., 2013). This dual role of snow as both a protective factor and a potential hazard could explain why increased winter precipitation generally supports migration but may pose risks under extreme conditions (Holtmeier and Broll, 2007). Warmer temperatures, combined with increased snowfall, as observed in some regions, can promote forest expansion but also lead to snow-related damage (Autio, 2006). These findings underscore how snow-rich environments can facilitate forest advancement by providing shelter and moisture during critical growth stages (Hagedorn et al., 2014; Kirdyanov et al., 2012; Devi et al., 2020; Grigoriev et al., 2022).

Enhanced facilitation by surrounding trees, as opposed to competition, has a strong positive effect on migration at sites CA and RU. This likely arises from positive interactions within the growing tree population, where clusters of trees can moderate harsh conditions by creating warmer microclimates, dampening wind, and accumulating snowpack, thereby insulating young trees from extreme cold and herbivory (Germino and Smith, 2000; Germino et al., 2002; Holtmeier and Broll, 2007; Bader et al., 2021; Zheng et al., 2024). Conversely, reduced wind exposure consistently shows a negative impact on migration at sites CA and RU. This could be due to the important role wind plays in dispersing seeds across the landscape, as less wind exposure may reduce the ability of seeds to travel efficiently, thereby hindering migration (Holtmeier and Broll, 2007). Although reduced wind exposure may alleviate physiological stress on established trees, it can also delay or even prevent treeline advance by reducing seed dispersal in wind-exposed environments where natural shelter is limited (Holtmeier, 2005a; Kullman, 2005).

The sensitivity analysis highlights the complex interactions between snow accumulation, wind exposure, and facilitation processes in driving treeline migration. These processes are further shaped by site-specific factors such as topography and climatic variability, which can either enhance or inhibit migration under changing environmental conditions. These findings highlight the necessity of a comprehensive approach when assessing treeline migration in snow-dominated regions.

Our enhanced model represents a major advancement in simulating alpine treeline processes. By integrating snow effects on seedling survival and soil conditions, alongside factors like wind exposure, facilitation by neighbouring trees, and microtopography, the model provides a more comprehensive and accurate framework for understanding treeline migration in high latitudes and elevations. This holistic approach is crucial for improving predictions and informing conservation strategies in sensitive ecosystems. To refine our understanding, we challenged the model by incorporating various processes derived from general models and previous studies to assess their effects. Although explicit validation was beyond the scope of this study, this approach identifies key areas for future investigation to improve the accuracy of treeline dynamics modelling.

4.4 Constraints on treeline advance

The complexity of treeline migration is highlighted in our results, which show that site-specific climatic factors have the strongest influence on migration rates, though the relative importance of these factors varies by location. This aligns with the widely accepted view that the elevational position of treelines is climate-dependent in the absence of significant disturbance regimes (Körner, 2020; Maher et al., 2021). However, our study underscores that treeline dynamics are not driven by climate alone but also by additional environmental variables, such as local site conditions and the impacts of disturbances, both natural and anthropogenic. These non-climatic influences complicate the relationship between climate warming and treeline advance, as certain detrimental factors cannot be offset by even the positive effects of a warmer environment (Holtmeier and Broll, 2005, 2007, 2019; Hansson et al., 2021).

A key finding of our research is the importance of snow-related processes, which have been incorporated into the LAVESI snow module, in affecting tree migration rates within the treeline ecotone. Wieser et al. (2019) suggest that, as mean surface temperatures rise, the influence of extreme weather events (e.g., early or late frosts or summer soil dryness) may become more significant than gradual temperature increases. These extreme events, coupled with increasing snow cover due to climate humidification, may act as a constraint on treeline advance in alpine regions.

Moreover, research on forest dynamics in high-elevation ecosystems highlights the need to consider both long-term warming trends and climate variability, as these drive episodic dieback and tree invasions (Bugmann and Pfister, 2000; Millar et al., 2004). While changes in climate averages are important, episodic events such as extreme cold episodes or changes in moisture availability can profoundly impact treeline dynamics. These episodic effects have often been underrepresented in vegetation models.

This points to a potential explanation for the time lag between climate warming and observed treeline migration, as these short-term events and longer ecological processes complicate the response of tree populations to warming. The updated LAVESI version addresses this gap by explicitly accounting for the influence of varying snow depths on tree growth and treeline migration. By incorporating stochastically occurring extreme events, the model captures the complete spectrum of weather variability.

Another significant contributor to the time lag in treeline migration is the slow pace of population processes, such as seed dispersal, establishment, and growth. These slow processes can result in centuries-long delays between individual tree development and broader forest establishment (Lloyd, 2005; Holtmeier and Broll, 2007, 2019). Favourable conditions for seedling establishment may not align with good seed years, further delaying the upward shift of the treeline (Gruber et al., 2022). Despite these delays, over time, slow processes at the advancing forest edge will eventually bring tree populations into alignment with the climatic treeline and its fundamental niche (Körner, 2021). The persistence of relict populations from past warm periods, such as northernmost dwarf spruces in the Canadian tundra, further supports this idea of delayed treeline responses (Nichols, 1976; Holtmeier, 1985). Similarly, the legacy of past environmental conditions strongly shapes the distribution of plant communities in both the treeline ecotone and lower alpine regions (Hofgaard and Wilmann, 2002). The integration of snow processes into the LAVESI model enhances the accuracy of simulations related to seed dispersal, establishment, and growth, ultimately improving predictions of treeline advancement. However, the localised legacy effects of past populations vary significantly at small spatial scales. Addressing these variations requires further region-specific research to refine vegetation models and enhance their predictive capabilities.

Dial et al. (2022) propose that part of the delay in treeline advance may be due to conifers being on the brink of a stochastic, climate-driven invasion of tundra after centuries of stability. However, many treelines worldwide show slower upward shifts compared to the rate of densification, with tree populations rapidly increasing under favourable conditions (Kruse et al., 2016; Wang et al., 2016; Shi et al., 2022). The LAVESI model is capable of simulating not only treeline migration but also tree density within the specified simulation area. Additionally, shifts in reproductive methods – such as a move toward seedling-based regeneration in species like Engelmann spruce (Picea engelmanii) and subalpine fir (Abies lasiocarpa) – could facilitate greater dispersal distances and potentially accelerate treeline advancement in the future (Holtmeier and Broll, 2007). An upward shift of the alpine treeline in cold climate zones may also be limited by the fact that while short vegetation thrives in the relatively warmer microenvironment near the ground, trees can only reach full size if they can withstand the surrounding ambient temperatures (Wilson et al., 1987; Grace, 1989). Incorporating such characteristics that are partly species-specific and partly population-specific would enhance the accuracy of vegetation models. Moreover, integrating insights into the ecological processes that shape spatial patterns within treeline ecotones – particularly those related to tree cover, clustering, and stature – could further improve the prediction of treeline dynamics under changing climatic conditions (Bader et al., 2021). Our study sites are characterised by conifer-dominated treelines, comprising predominantly evergreen species in North America and deciduous conifers in Siberia. In ecotones with mixed forest composition, warming-induced acceleration of successional shifts from broadleaved to coniferous dominance warrants particular attention (Sigdel et al., 2024). Our findings underscore the complexity of treeline dynamics, suggesting that both climate factors, including snow-related processes, and slow ecological processes, such as seedling establishment, growth, survival, reproduction, and forest development, play crucial roles in shaping treeline responses to climate change. The gaps highlight the need for further refinement of vegetation models to better explain the time lag between climate warming and treeline advance.

5 Conclusions

In conclusion, this study highlights the distinct site-specific responses to factors influencing treeline shift and forest expansion, emphasising the critical role of localised conditions in shaping migration dynamics. While the findings reveal commonalities across sites, they also underscore unique interactions between environmental factors and local contexts. Notably, the responses observed at the Canadian and the Russian sites provide clearer insights into the primary drivers of migration, whereas the high variability at the Alaskan site suggests more complex or less predictable local dynamics. This variability highlights the need for further investigation to reduce uncertainties in model predictions.

A key finding of this study is the significant role of snow in modulating migration potential, as snow accumulation creates favourable moisture conditions and protects seedlings from extreme cold, particularly in snow-dominated environments. These results emphasise the importance of incorporating snow-related processes into vegetation models to improve predictions of boreal forest dynamics under changing environmental conditions.

In boreal regions, altered snow regimes and increasingly frequent climatic extremes may affect treeline advance, with cascading effects on forest structure, biodiversity, surface albedo, and the long-term carbon sequestration potential of these ecosystems. As warming accelerates, these dynamics are likely to intensify, highlighting the need for refined, region-specific models to anticipate ecological feedbacks.

Overall, this study enhances our understanding of tree migration processes and underscores the varied predictability of migration responses across different sites. These insights carry valuable implications for enhancing model accuracy and guiding conservation strategies aimed at sustaining alpine tundra resilience in the face of rapid environmental change.

Appendix A

Table A1 Overview of simulation scenarios. Summary of the simulation scenarios and their respective multipliers for predicting variables used in the sensitivity analysis.

Appendix B

Figure B1 Stem count [trees ha⁻¹] over time. Stem count from 2020 to 2095 CE at the study site Fort McPherson West, Canada (CA), under (a) the reference scenario (run 34) and (b) the scenario with snow computation enabled (run 2). Stem count from 2020 to 2095 CE at the study site Lake Ilirney, Russia (RU), under (c) the reference scenario (run 13) and (d) the scenario with snow computation enabled (run 3). Stem count from 2020 to 2095 CE at the study site Road to Central, Alaska, USA (AK), under (e) the reference scenario (run 15) and (f) the scenario with snow computation enabled (run 3). (Runs based on the median migration rate calculated from 40 simulations per scenario). Each plot represents data from the simulation run with the median migration rate derived from 40 simulation runs per scenario. Stem count values are shown for the tree species occurring at the respective study site.