1. Ecology

Regional opportunities for tundra conservation in the next 1000 years

  1. Stefan Kruse  Is a corresponding author
  2. Ulrike Herzschuh
  1. Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research, Germany
  2. Institute of Environmental Sciences 6 and Geography, University of Potsdam, Germany
  3. Institute of 7 Biochemistry and Biology, University of Potsdam, Germany
https://doi.org/10.7554/eLife.75163
Share this article
Cite this article
  1. Stefan Kruse
  2. Ulrike Herzschuh
(2022)
Regional opportunities for tundra conservation in the next 1000 years
eLife 11:e75163.
https://doi.org/10.7554/eLife.75163
  2. Download BibTeX
  3. Download .RIS

Abstract

The biodiversity of tundra areas in northern high latitudes is threatened by invasion of forests under global warming. However, poorly understood nonlinear responses of the treeline ecotone mean the timing and extent of tundra losses are unclear, but policymakers need such information to optimize conservation efforts. Our individual-based model LAVESI, developed for the Siberian tundra-taiga ecotone, can help improve our understanding. Consequently, we simulated treeline migration trajectories until the end of the millennium, causing a loss of tundra area when advancing north. Our simulations reveal that the treeline follows climate warming with a severe, century-long time lag, which is overcompensated by infilling of stands in the long run even when temperatures cool again. Our simulations reveal that only under ambitious mitigation strategies (relative concentration pathway 2.6) will ∼30% of original tundra areas remain in the north but separated into two disjunct refugia.

Editor's evaluation

This study will be of interest to researchers studying climate change and Arctic tundra systems. The authors apply LAVESI, a machine-intensive and spatially explicit simulation of individual Siberian trees at the tundra-forest boundary, to call attention to the rapid reduction in the tundra biome as climate warming pushes forests toward the Arctic Ocean. This detailed modelling study predicts dramatic losses of tundra area by the middle of the millennium even under an ambitious climate mitigation scenario and highlights considerable risks of extinction.

https://doi.org/10.7554/eLife.75163.sa0

Introduction

Arctic regions have experienced strong warming over recent decades (Box et al., 2019). The warming is expected to continue drastically based on future changes simulated with prescribed ‘relative concentration pathway’ (RCP) scenarios (Meinshausen et al., 2011). Following climate simulations, mean July temperatures in the region could be, by the end of the 21st century, less than +2°C warmer with ambitious mitigation efforts (RCP 2.6), about +3.1°C (intermediate effort scenario RCP 4.5), or an extreme of +14°C warmer (high-emission high-warming RCP 8.5). In consequence, the treeline ecotone can expand and will squeeze Arctic regions from the south and have only limited refugial potential in the north (Hofgaard et al., 2012; MacDonald et al., 2007) as simulation studies predicted (Rupp et al., 2001). This turnover threatens the existence of various types of sensitive tundra with differing vegetation structure along latitudinal and elevational bioclimatic gradients located between the treeline in the south and the Arctic Ocean in the north (Walker et al., 2005; Pearson et al., 2013; Yu et al., 2009). The Siberian tundra areas are known for their high regional floristic diversity and species biodiversity (Pauli et al., 2012; Mod and Luoto, 2016; Arctic Climate Impact Assessment, 2004; Schmidt et al., 2017) and support special (indigenous) types of land use (e.g. reindeer herding, foraging). However, the expansion of forests in response to increasing temperature is not well understood and the uncertain timing of future tundra loss is challenging conservation efforts.

In Siberia, the climate velocity of isotherms within the tundra biome has been calculated as ∼290 m year−1 (Loarie et al., 2009). When assuming a direct equilibrial relationship linking the modern position of the treeline to its main limiting factor of July temperature (MacDonald et al., 2007), the treeline is expected to reach its maximum expansion within the first centuries of the millennium. Regionally, where the shoreline of the Arctic Ocean is close to the treeline, the 100–200-km-wide tundra corridor would be completely covered by forests (Figure 1). Moreover, even regions with wide tundra corridors in Siberia (∼600 km; Taimyr Peninsula and Chukotka) could be forested in the warmest scenarios. However, the potential migration rates are likely constrained by other factors, partly acting on a small scale (Corlett and Westcott, 2013). This and the fact that in Siberia the treeline is formed solely by larch species growing on permafrost while in North America this niche is covered by a mix of species (Mamet et al., 2019; Herzschuh and Jordan, 2019) make it difficult if not impossible to transfer knowledge from any other treeline region to inform Siberian tundra responses.

Figure 1
Download asset Open asset
Transects (blue lines) were placed starting in the treeline at field sites and extending to the shoreline of the Arctic Ocean (map); only tundra areas above the treeline (Walker et al., 2005) are considered.

The land area was grouped into four regions equidistantly separated between the transects (lower plot) and plots show climate forcing for 501–3000 CE based on relative concentration pathway (RCP) scenarios in these regions at the extremes and middle part of the transects; numbers give mean values for modern and future strongest warming under RCP 8.5. Map projection: Albers Equal Area.

Equilibrium-based estimations are not backed by empirical studies which are showing slower responses and either no treeline advance or rates of only a few metres per year in Siberia (Scherrer et al., 2020; Kharuk et al., 2018; Wieczorek et al., 2017b; Rees et al., 2020). Either the different response of regions can be related to site conditions (e.g. Rees et al., 2020; Mamet et al., 2019) or it is not well understood how complex processes such as seed dispersal or intra-species competition acting on fine and broad scales constrain the response to warming (Corlett and Westcott, 2013; Wieczorek et al., 2017b). Models that do not consider these processes lead to forecasts prone to overestimating treeline advance (compare Pearson et al., 2013; Snell and Cowling, 2015), but it is challenging to consider both small-scale processes and the broad-scale picture.

Models which contain processes to estimate competition between individuals for resources and seed dispersal in a spatially explicit environment and which handle each tree individually are called individual-based spatially explicit models (DeAngelis and Mooij, 2005). These models have become increasingly used lately to assist in formulating conservation management strategies (Zurell et al., 2021). For the northern Siberian treeline ecotone, the vegetation model LAVESI was developed, which incorporates a full tree (Larix) life cycle beginning with seed production and dispersal, germination, and establishment of seedlings, which then grow based on abiotic (temperature, precipitation, active-layer depth) and biotic (competition) conditions (Kruse et al., 2019a, Kruse et al., 2016). Although soil nutrients and mycorrhiza are not considered explicitly in the current version, it was shown that it can capture the nonlinear response of the treeline populations well at a regional level on the Taimyr Peninsula (Kruse et al., 2019a, Kruse et al., 2016; Wieczorek et al., 2017b).

To cope with the complex forest-tundra ecotone responses, which are not in equilibrium with climate (Hofgaard et al., 2012), we updated and made use of our individual-based spatially explicit model LAVESI. With our simulation study, we explore treeline migration and tundra area dynamics to provide important information as to which tundra areas can survive future warming and where preventive measures for protection would be necessary to conserve the unique Siberian tundra. This study is guided by two questions: ‘What are the treeline trajectories in Siberia for different climate mitigation scenarios (RCP)?’ and ‘Under which climate scenario can tundra retreat to refugia and recolonize areas after temperatures cool again, and where?’.

Results and discussion

Treeline migration trajectories under future climate changes

The simulations revealed that the treeline is not in an equilibrial relationship but its advance follows climate warming, potentially reaching 30 km decade−1 in the 21st century under warmest climate scenarios (Appendix 2—table 1). When relating and tracking the position of the treeline to the climate (=climate-analogue), we find a severe migration lag of hundreds of years until simulations reach the same position (Figure 2). This behaviour is longer than expected from the properties of treeline-forming larch species with generation times of decades and weak long-distance dispersal rates (e.g. Mamet et al., 2019; Kruse et al., 2019a; Scherrer et al., 2020; Clark, 1998). Other processes, not yet explicitly captured by the model, may constrain the migration process, such as the strong greening of the Arctic increasing interspecific competition (compare Pearson et al., 2013; Berner et al., 2020), the warming and abrupt thaw of permafrost (e.g. Stuenzi et al., 2021; Biskaborn et al., 2019), soil nutrients (e.g. Sullivan et al., 2015), and animal activity (e.g. Wielgolaski et al., 2017). Regardless, our simulated rates over the past few decades are generally in line with the few existing modern observations showing only a slow treeline advance (Mamet et al., 2019; Hansson et al., 2021; Rees et al., 2020; Shevtsova et al., 2020; Guo and Rees, 2020). Rates are especially high in Chukotka, however, in all forcing scenarios (Appendix 2—table 1): an observation that can be attributed to forest patches ahead of the current treeline, which serve as nuclei for tundra infilling (shown by, e.g. Snell and Cowling, 2015; Kharuk et al., 2018; Kruse et al., 2019a; Harsch and Bader, 2011). Generally, the simulated rates are high compared to those during post-glacial migration (e.g. Feurdean et al., 2013) but seem likely as future climate is unprecedented, especially in their rate of change.

Figure 2
Download asset Open asset
The trajectory of the simulated treeline position relative to its current position and the maximum at the shoreline versus its climate-analogue position for the four regions shows a migration lag of the treeline during the first centuries (each line segment represents 25 years and the length of the arrow head corresponds to the step length) until the simulated treeline is limited by climate.

Forests expand their area further and infilling proceeds when climate conditions cool and even overshoot in the long run with cooling back to 20th-century temperatures. The diagonal is where climate and the treeline are in equilibrium; below the diagonal, tree migration lags climate; above the diagonal is ‘overshooting’ and reaching locations actual climate would allow. For each relative concentration pathway (RCP) scenario, two are presented, one for the scenario as-is and the second for the cooling; scenario RCP 2.6* warms only at half the rate of RCP 2.6. See also Appendix 3—figure 1 for the year when the trajectory passes the equilibrium.

After the initial migration lag, we can see from our millennium-long simulations that the treeline reaches more advanced positions than the climate analogue (Figure 2, Appendix 2—table 1, Appendix 2—figure 1, Appendix 2—figure 2, Appendix 3—figure 1). This strong overshooting was unexpected especially in scenarios with cooling temperatures back to the 20th-century level. However, we designed our model to be pattern-oriented and is based on observations that, once established, trees can endure periods of unfavourable environments (compare Kruse et al., 2016). This leads to the survival of trees ahead of the treeline acting as nuclei for tundra colonization when conditions improve (Väliranta et al., 2011; Snell and Cowling, 2015). This is backed by observations of slower-than-expected retreat of elevational treelines because trees can persist until dramatic events cause the death of all well-established trees (Scherrer et al., 2020).

Tundra area dynamics caused by forest expansion

The area loss over time roughly follows a negative exponential function, which is steeper the warmer the climate scenarios are (Figure 3, Appendix 3—figure 1). The quickest loss, already by ∼2100 CE, is seen in the intermediate warming scenario RCP 4.5 (+3°C July temperatures) and the extreme RCP 8.5 (July temperatures +14°C). For these scenarios, tundra retreats to its minimum cover of 5.7% by the middle of the millennium (lower panels in Figure 3). Even under mitigation scenario RCP 2.6, the tundra area is reduced to 32.7% or in the very ambitious scenario 2.6* to only 54.9%. Other simulation studies show 50% loss during the 21st century (Callaghan et al., 2005; Kaplan et al., 2003; Pearson et al., 2013) but not the dynamics of the long run until the end of 3000 CE, as shown by our study.

Figure 3 with 1 supplement see all
Download asset Open asset
Forest position and tundra area at year 3000 CE for different climate mitigation scenarios and under potential cooling back to 20th-century temperatures after peak temperatures have been reached.

The area of tundra changes over time and can only partly recover when temperatures cool and forests recede (plots next of the maps show years 2000–3000 CE). Only tundra areas above the treeline (Walker et al., 2005) are considered. Map projection: Albers Equal Area.

Cooling back to the 20th century levels led to a slight tundra re-expansion of a few percent of the area by the end of the simulations. While re-expansion is small for RCP 2.6* and 2.6, and largest for RCP 4.5 on the Taimyr Peninsula where forests retreat from areas of ∼200,000km2 in the final 200 years of the simulation, under RCP 8.5 the effect of the cooling led to no tundra releases. In consequence, only if humans act ambitiously in limiting global warming to well below +2°C (RCP 2.6, 2.6*), and temperatures actually cool after the maximum, do simulations show that ∼10–30% of the tundra area could persist in the north and can re-expand after the mid-millennium maximum treeline extent. However, due to the overshooting effect, less than 10% of the area remains unforested. All in all, the expansion is in accordance with other studies and will likely put especially cold-climate tundra types at risk of extinction (e.g. Pearson et al., 2013; Walker et al., 2005).

Concluding opportunities for tundra conservation

Following our simulation results, tundra is reduced to ∼5.7% of its current area at maximum forest expansion under the most likely scenario RCP 4.5 after an initial migration lag. Only ambitious strategies to mitigate climate warming (e.g. RCP 2.6) could prevent such an extensive loss and retain ∼32.7% of the current tundra. In both cases, the northern refugia are located in Chukotka and the Taimyr Peninsula, currently covered by High Arctic and Arctic tundra types which have relatively high biodiversity (Arctic Climate Impact Assessment, 2004). A WWF report (Krever et al., 2009) finds that the current total cover of protected northern areas is insufficient to achieve the requisite 30% protection necessary for biodiversity conservation (Secretariat of the United Nations Convention on Biological Diversity, 2021). Especially cold-climate tundra types on Siberian islands and in Chukotka are underprotected (e.g. Pearson et al., 2013; Walker et al., 2005) and thus further protected areas need to be established.

Although the large refugia could enable tundra species survival and later recolonization (analogon forests, Harsch et al., 2009; Clark, 1998; Stralberg et al., 2020), the diversity in these ~2500 km distant fragments is threatened by the fundamental disadvantages of a small population (inbreeding/genetic drift, Ohsawa and Ide, 2008; Charlesworth and Charlesworth, 1987; Mona et al., 2014; May et al., 2013). A network of connected systems would be necessary with protection of potential refugia (e.g. Morelli et al., 2020) distributed across the ~4000-km-wide tundra area. Due to its individual-based nature, LAVESI can be a tool to assess population genetics and aid in the optimal placement of migration corridors and areas to protect (cf. Zurell et al., 2021).

Materials and methods

Model description

Request a detailed protocol

The model LAVESI is an individual-based and spatially explicit model that was developed to study the treeline ecotone growing on permafrost soils on the Taimyr Peninsula (Kruse et al., 2016; Wieczorek et al., 2017b; Epp et al., 2018) and since then it has been improved and modified for different purposes including for migration rate studies (Kruse et al., 2019a). To achieve the most realistic model, each life history stage of the larch individuals is handled explicitly, and the model’s processes adapted to observed patterns of surveyed tree stands for the dominant larch tree species, Larix gmelinii and Larix sibirica.

The simulation proceeds in yearly time steps from the beginning to the end of the input climate time series following a stabilization period to ensure that emerging populations reach equilibrium with the environment. Stochasticity in the model was introduced by using random numbers generated with a pseudo random number generator.

Within one simulation year, the following processes become consecutively invoked (Figure 4):

  • Update of environment: Interactions between neighbouring trees are local and indirect. Basal diameters of each individual tree are used to evaluate the competition strength. We use a yearly updated density map to pass information about competition for resources between trees. An active-layer depth is estimated based on an edaphic factor and the number of days exceeding 0°C following simplifications by Hinkel and Nicholas, 1995.

  • Growth: The maximum possible growth is estimated based on 10 years mean climate auxiliary variables (temperature of the coldest [January] and warmest [July] month, N⁢D⁢D0 ‘net degree days’, A⁢A⁢T1⁢0 ‘active air temperature’). The individual growth of basal diameter and, if a tree reached a height of 1.3 m, breast height diameter is calculated from the maximum possible growth in the current year affected by the tree’s density index. From the resulting diameters, the tree height is estimated.

  • Seed dispersal: Seeds in ‘cones’ are dispersed from the parent trees. The dispersal directions and distances are randomly determined from a ballistic flight influenced by wind speed and direction with decreasing probabilities for long distances, and, if dispersed seeds leave the extent of the simulated plot, they can be introduced from the other side only on the east–west margins.

  • Seed production: Trees produce seeds after the year at which they reached their stochastically estimated maturation height. The total amount depends on weather, competition, and tree size.

  • Establishment: The seeds that lie on the ground germinate at a rate depending on current weather conditions.

  • Mortality: Individual trees or seeds die, that is, they become removed from the plot, at a specified mortality rate. For trees this is deduced from long-term mean weather values, a drought index, surrounding tree density, tree age and size, plus a background mortality rate. Seeds, on the other hand, have the same constant mortality rate whether on trees or the ground.

  • Ageing: Finally, the age of seeds and trees increases once a year and seeds are removed from the system when they reach a defined species age limit.

Figure 4
Download asset Open asset
The model LAVESI invokes in each yearly time step processes to individually handle seeds and represent the full life cycle of simulated trees, leading to forest stand dynamics.

The climate environment drives establishment, growth, and mortality.

Many factors influence the success of tree establishment and growth, of which the following seem most important to understand treeline migration. Warmth, especially summer mean temperatures, cumulative temperatures, and the vegetation period length, are key climate variables (e.g. MacDonald et al., 2007) and are hence included in the model for estimating potential tree growth and tree mortality (Kruse et al., 2016). This is further influenced by the abiotic (permafrost, active-layer depth) and biotic (intraspecific competition) environment. For example, permafrost soils with a shallow active-layer depth limit the space available for rooting and nutrient stores (e.g. Sullivan et al., 2015) and frozen/cold soils prevent growth but water can support growth in dry summers (Sugimoto et al., 2002). The latter functioning could be lost in future hot environments, and in consequence it would decrease the treeline advance rate further and thereby increase the lag time. Accordingly, the active-layer depth is explicitly simulated in the model’s environment module, but soil nutrients are considered only implicitly. Individuals compete for space which serves as a surrogate for this. Permafrost degradation/abrupt thaw and warming (Biskaborn et al., 2019; Stuenzi et al., 2021) could lead to thermokarst processes that would locally cause swamping/waterlogging, leading to suppressed establishment/growth (Rees et al., 2020), which is not explicitly simulated and unlikely to limit tree invasion on a broader scale. Further limitations could arise from the lack of a suitable mycorrhiza symbiosis partner, but as shown by Hippel et al., 2021 even in northernmost areas on the Taimyr Peninsula mycorrhiza are present even before the invasion of trees in the Holocene. Additionally, dispersal rates of mycorrhiza with small spores are a magnitude faster than for taxa with larger seeds so that geographic distance is unlikely to be a limiting factor, similar to diatoms as shown by Stoof-Leichsenring et al., 2015. On the other hand, interspecific competition of recruits with present shrubs or grasses which are increasing more quickly – called the greening of the Arctic (Pearson et al., 2013; Berner et al., 2020) – is not explicitly simulated in LAVESI, although it is partly included in the germination probability and survival rate of individuals at the seedling stage. Further biotic interactions are implicitly included, such as herbivory by browsers (e.g. Wielgolaski et al., 2017), which reduces the survival and growth rate by damaging the shoots (compare krummholz study by Kruse et al., 2020) and can lead to a snow compaction cooling effect on soil temperatures (Beer et al., 2020), or pest outbreaks that are likely to increase in the future (e.g. larch silk moth, Fedotova et al., 2019). This would most likely lead to a further reduction of migration rates.

Model improvements

Request a detailed protocol

The source code was updated to LAVESI-WIND version 1.2 on the crutransects branch to be sufficiently light in terms of memory allocation (details in Appendix 1) and including landscape (developed in parallel for Chukotka mixed alpine-latitudinal treeline simulations in Shevtsova et al., 2021; Kruse et al., 2021). The model is freely available on GitHub (https://github.com/StefanKruse/LAVESI; Kruse et al., 2022). The source code of the version used in this article is permanently stored on Zenodo (https://doi.org/10.5281/zenodo.6344261).

Transect setup

Request a detailed protocol

At four representative locations of the Siberian treeline within 100–170° E, we placed linear transects, starting at a field site, and allowed open tree stands to establish up to the shoreline (Figure 1). From west to east: (1) Taimyr Peninsula, 573 km long, site 13-TY-09-VI, 72.15067∘ N (Wieczorek et al., 2017b, Wieczorek et al., 2017c), (2) Buor Khaya Peninsula, 137 km long, 14-OM-20-V4: 70.52671∘ N (Liu et al., 2020), (3) Kolyma River Basin, 146 km long, 12-KO-02/I: 68.38916∘ N (Wieczorek et al., 2017a), and (4) Chukotka, 626 km long, EN18022: 67.40102∘ N(Kruse and Stoof-Leichsenring, 2016; Kruse et al., 2019b). While the first three locations in the west cover predominantly flat terrain, the easternmost location in Chukotka has mountainous stretches.

Climate forcing

Request a detailed protocol

Temperature and precipitation time series were constructed using the following steps:

  1. Extract climate data in 10 km steps along transects beginning at the position of the field sites for each of the four focus regions from the grid cells of CRU TS 4.03 by distance-weighted interpolation (Harris et al., 2020).

  2. Establish a yearly δ18 O and ice layer thickness series for 501–1900 from the yearly resolution ice-core data from Severnaya Zemlya (north of the Taimyr Peninsula, available for years 934–1998; Opel et al., 2013), the only closest annual resolution archive allowing both temperature and precipitation reconstructions independent of plant-based proxies (pollen or tree rings). This was achieved by sampling 25-year-long blocks of yearly resolution data from the pre-industrial (<1800CE) and adding this variability to the AN86/87 core which dates back to 4703 BCE (Arkhipov et al., 2008).

  3. Adjust this locally to the 10 km steps climate data series per transect by building linear models of yearly temperature and precipitation for the overlapping years 1901–1998 and use this to estimate mean values for the period 501–1900.

  4. Sample 15-year-long blocks from CRU and adjust the mean temperature by the difference of the sampled block from the estimated mean values to achieve monthly data series.

  5. Prolong this series by coupled climate model long runs until 2300 CE with model output with MPI-ESM-LR prepared for CMIP5 for scenarios RCP 2.6, RCP 4.5, and RCP 8.5 and additionally a scenario that warms only at half the rate of RCP 2.6, which is here called RCP 2.6*.

  6. Extend these trends for 2100–2300 until 2500 CE and repeat in the following either the climate data of 1901–1978 in a loop, which in this study is called 20th-century cooling, or copy the period 2300–2500 until 3000 CE.

July temperatures (Figure 1), which have the strongest impact in the model due to a positive influence on tree diameter growth, increase in these scenarios by a median of +1.2°C (RCP 2.6*), +1.8°C (2.6), +2.1°C (4.5), and up to +5°C (8.5) compared to 1970–2000 until 2100, which partly cools until 2300 to +0.6°C (RCP 2.6*), +0.5°C (2.6), and reaches higher levels for the two warmest scenarios to +3.1°C (4.5) and up to +144°C (8.5).

Wind speed and direction data were obtained at 6-hourly resolution for 1979–2018 from EraInterim5 (dataset description in Dee et al., 2011) for each of the transects and supplied with the model for estimation of seed dispersal events. In a prior version, these data were only available for the Taimyr Peninsula and only for years until 2012.

Future climate changes will lead to an increase in extreme weather events (Masson-Delmotte et al., 2021), among others a higher ignition probability and thus wildfires are expected to which the Siberian boreal forests are sensitive (e.g. Shvidenko and Schepaschenko, 2014). Fire is not yet explicitly simulated in the model and in a hotter environment fires could lead to dramatic forest losses. Areas might be turned into steppe and not revert to the prior tundra. Increases in extreme rainfall events will lead to local flooding and waterlogging (Masson-Delmotte et al., 2021), a response not explicitly implemented in the model. The anticipated effect to the LAVESI predictions is that tundra colonization would most likely be slowed down in the lowlands but less so in the well-drained mountainous regions, which would further prolong the time-lagged response of the treeline migration.

Tree growth on transects

Request a detailed protocol

We implemented sensing of the environment for each individual tree via its y-position along the transect which is made up of 10 km spaced climate data. The maximum possible growth at a certain position (see initial publication Kruse et al., 2016) is calculated by interpolating the possible growth under the climate using the two closest climate variables.

The larch species L. gmelinii (Rupr.) Rupr. dominates the areas to the west of the Verkhoyansk Mountain Range ∼90–130° E and its sister species L. cajanderi Mill. grows in northeast Siberia Figure 1 (Abaimov, 2010). Hence, the already implemented species L. gmelinii was set to be present in the simulations and simulations tuned to fit observations. To initiate growth 100,000 seeds were introduced at the cold end with a distance based on a negative exponential kernel for the first 50 years. Similarly, extinction on a transect was prevented by permanently introducing 1000 seeds per 200 km length. Further seed introduction from the nonexplicitly simulated hinterland was estimated on a 500-m-long stretch at the cold site linking the production and release height to climate.

The model used for this study is solely forced by climate data and no vegetation feedback on permafrost soils and climate was implemented. Therefore, the expected decrease of the albedo leading to a positive feedback with climate when tundra is colonized by tree stands (e.g. Pearson et al., 2013) is neglected here. Considering this would likely further increase the rate of transition from tundra to forests.

Model tuning and validation

Request a detailed protocol

Simulation runs were forced with the prepared climate data for the transects for validation of the shape of the treeline as simulated by the model. The results of the year 2000 were compared to field data that was positioned along the temperature gradient on the simulated transects by using the respective CRU TS temperature data for each grid cell of each plot (example Figure 1). The field data collection included disturbed sites (e.g. in Chukotka), leading to very high stem counts, which are the number of trees exceeding 1.3 m in height per hectare. The model could be tuned by introducing a local adjustment based on elevation lapse rates of the relevant climate values used internally, which is necessary as the climate data may lead to over-/underestimations locally in a model that does not consider local topography: an issue that is especially relevant in the mountainous region of Chukotka. We made a moderate fit with the following amendments:

  • Taimyr Peninsula: TJan=−1.17∘C, TJul=−0.84∘C and PYear=−5.4 mm

  • Buor Khaya Peninsula: TJan=−0.23∘C, TJul=−2.48∘C, and PYear=−28.9 mm

  • Kolyma River Basin: TJan=+3.14∘C, TJul=+1.26∘C, and PYear=+85.7 mm

  • Chukotka: TJan=+4.46∘C, TJul=+4.30∘C, and PYear=+8.2 mm

From the individuals present along the simulated transects, the northernmost positions of the treeline are calculated. For the four focus regions these are +66, +88, +98, and +16 km beyond the treeline start location at the year 2000 (Appendix 2—figure 2), which are partly ahead of the observed positions based on visual inspection of satellite data ~+30, +30, +80, and +15 km beyond the field location selected as the starting point.

We assume that the fundamental drivers of tree growth remain stable in time. The space-for-time approach used here is based on using several sites from across the treeline in the regions of interest for model parameterization and validation that have experienced a range of past climate (Likens, 1989). However, climate warming will lead to unprecedented temperatures (Figure 1) that will affect the current stand structure and distribution. This may cause changes in the response of single considered processes or make other processes that are not currently important necessary (e.g. Blois et al., 2013; Likens, 1989). Keeping this in mind, the future simulations should be interpreted with caution, although a simulation study with LAVESI revealed realistically predicted migration rates and patterns for the recent decades of strong warming in the study region of Chukotka (currently in review).

Simulation setup and data processing

Request a detailed protocol

We ran three simulation repeats for each transect with the climate series prepared for this study as initial tests showed robust simulation results. To consider the available area from the start to the shoreline, we used 800-km-long and 300-km-long transects and 20-m-wide wrapping for the east and west boundary of each simulated stretch, thereby representing a slice of a continuous treeline.

The positions of three key variables of stand densities expressed in stems, which are trees >1.3 m tall, are extracted in 10-year steps: single-tree stands, defined as the northernmost position of forest islands ahead of the treeline with ≥1 stem ha−1; treeline stands, the northernmost position of a continuous forest cover with 1 stem ha−1≤density<100 stem ha−1; and forestline stands, the northernmost position of forest cover >100 stems ha−1 (see Figure 2 in Kruse et al., 2019a for a graphical representation). The position of the simulated treeline determined for the year 2000 is used as the baseline for the calculation of migration rates. The resulting values for each of the three key variables are used in 10-year steps to interpolate with a weighted average the expansion starting at the current treeline (e.g. Walker et al., 2005). This is similar to the biotic velocity of Ordonez and Williams, 2013. Opposed to this, we extrapolated the position of the modern observed treeline based on July temperatures (=climate-analogue) for the years 2000–3000 CE. Therefore, the current July temperature at the treeline position is determined and tracked along the transect at each time step of the climate forcing, which is similar to the climatic velocity of Loarie et al., 2009.

We used R version 3.6.1 (R Development Core Team, 2019) for data handling, statistical analyses, and graphical representation. The 30 arcsec WorldClim 1.4 data (Hijmans et al., 2005) along the transects and 2.5 min of a degree elevation for mapping was downloaded through the geospatial package raster version 3.0-12. (Hijmans, 2020). Additional geospatial packages sf version 1.4-2 and rgdal version 1.4-7 (Pebesma, 2018; Bivand et al., 2019) were used to reproject the treeline shape to Albers Equal Area projection and calculate weighted average buffers at each 10-km step along the treeline and merge the buffer polygons for each of the advancement values for the three key variables (single trees, treeline, forestline). The treeline constructed by the circumarctic vegetation map consortium (CAVM, Walker et al., 2005) is provided in segments that were simplified in QGIS version 3.10 (QGIS Development Team, 2020). For tundra area assessment, the land area polygons of Russia accessed via R package maps version 3.3.0 (code by Richard and Becker, 2018) were cleaned and the islands in the Arctic Ocean excluded. Linear models for temperature and precipitation reconstruction and climate-lapse rates were established with functions from the R base package (R Development Core Team, 2019).

Appendix 1

Model performance

During simulation runs, each Environment grid cell of 20 × 20 cm needs 10 bytes, an element of the Tree structure 64 bytes, and the Seed structure 32 bytes. In total, for very dense scenarios it needs for each square kilometre ~11 GB RAM. The parallelization was updated and instead of Standard Template Library (STL) list containers for tree and seed elements a vector structure was implemented to allow better support for OpenMP (https://www.openmp.org) parallelization for scaling on a high-performance cluster and the realization of large-scale transect simulations as needed for this study. The computation speed depends on the number of trees and reaches ∼3.5 s year−1 km−2 for very dense forests on eight cores that can be reduced when using more cores (Appendix 1—figure 1). However, due to overheads, scaling is not 1:1. Of the internal functions of LAVESI, the mortality is most computationally intensive as each seed and tree need to be handled individually.

Appendix 1—figure 1
Download asset Open asset
The computation speed increases with the number of available threads and reaches a plateau caused by overheads.

The processes of LAVESI scale differently and the process Environment is computationally the heaviest and thus relatively faster with more cores while processes such as Mortality show no clear benefit of more threads. Data were aggregated for 100 years at three different population dynamics stages: e, empty; m, mature old growth; d, densification.

Appendix 2

Migration rates

Appendix 2—table 1 provides a summary of simulated treeline advance and the corresponding position at the same year of the climate analogue. Appendix 2—figure 1 shows the position of the climate-analogue treeline position at the four transects, and Appendix 2—figure 2 gives an overview of the simulated maximum positions of the single-tree stands, treeline, and forestline at the same transects.

Appendix 2—table 1
Summary of simulated vs. climate envelope-based treeline advance.

Climate envelope solely based on July temperatures (following MacDonald et al., 2007). Values in bold reach the shoreline, but note that because of technical reasons the step sizes of the climate are in 10-km steps. Values are the result of three simulation repeats, and the standard deviation is stated next to the mean value.

TransectScenarioMaximum potential expansion based on climate and yearSimulated maximum expansion and yearTreeline migration rate 2000–2100 km decade−1Treeline migration rate 2000–2300 km decade−1
TaimyrRCP 2.6*210 km 2090192.2 ± 18.5 km 2733 ± 1620.8 ± 0.71.5 ± 0.1
RCP 2.6*c210 km 2090135.3 ± 1.4 km 2357 ± 583.4 ± 2.91.9 ± 0
RCP 2.6320 km 2090344.8 ± 2.9 km 2730 ± 8711.5 ± 2.39 ± 0
RCP 2.6c320 km 2090331.3 ± 1.4 km 2283 ± 298.6 ± 1.28.6 ± 0.1
RCP 4.5410 km 2297544.5 ± 0 km 2383 ± 623.4 ± 0.210.6 ± 0.1
RCP 4.5c400 km 2221405.3 ± 3.8 km 2567 ± 23124.7 ± 0.211.2 ± 0.1
RCP 8.5560 km 2082544.5 ± 0 km 2160 ± 029.1 ± 0.516.1 ± 0
RCP 8.5c560 km 2082544.5 ± 0 km 2140 ± 028.4 ± 0.416.2 ± 0
Buor KhayaRCP 2.6*130 km 204094.5 ± 0.9 km 2307 ± 750.1 ± 0.20.2 ± 0
RCP 2.6*c130 km 204091 ± 0.9 km 2497 ± 3180 ± 00 ± 0
RCP 2.6130 km 2020120.2 ± 1.4 km 2247 ± 121 ± 0.31.1 ± 0
RCP 2.6c130 km 2020122.3 ± 0.3 km 2633 ± 920.9 ± 0.61.2 ± 0.2
RCP 4.5140 km 2133124.5 ± 0 km 2303±60.6 ± 0.51.3 ± 0.1
RCP 4.5c140 km 2133124.5 ± 0 km 2287 ± 120.2 ± 0.11.2 ± 0.1
RCP 8.5140 km 2038124.5 ± 0 km 2120 ± 03.2 ± 01.2 ± 0
RCP 8.5c140 km 2038124.5 ± 0 km 2123 ± 123 ± 0.31.1 ± 0
KolymaRCP 2.6*130 km 2012120.8 ± 0.3 km 2330 ± 01.2 ± 0.10.7 ± 0
RCP 2.6*c130 km 2012119 ± 0.9 km 2217 ± 121.1 ± 0.20.7 ± 0
RCP 2.6150 km 2094134.3 ± 0.3 km 2230 ± 02.1 ± 0.11.2 ± 0
RCP 2.6c150 km 2094113.5 ± 17.3 km 2103 ± 921.3 ± 0.61 ± 0.1
RCP 4.5150 km 2284134.5 ± 0 km 2257 ± 122 ± 0.11.2 ± 0
RCP 4.5c140 km 2121134.5 ± 0 km 2307 ± 61.6 ± 0.51.2 ± 0
RCP 8.5150 km 2072134.5 ± 0 km 2143 ± 62.2 ± 0.21.2 ± 0
RCP 8.5c150 km 2072134.5 ± 0 km 2140 ± 02.2 ± 0.11.2 ± 0
ChukotkaRCP 2.6*300 km 2047285.2 ± 8.4 km 2173 ± 2925.1 ± 0.18.8 ± 0.1
RCP 2.6*c300 km 2047285.7 ± 3.2 km 2203 ± 2925.3 ± 0.29.2 ± 0.1
RCP 2.6350 km 2047363.5 ± 0 km 2423 ± 630.8 ± 111.3 ± 0.2
RCP 2.6c350 km 2047358.8 ± 6.4 km 2290 ± 8730.9 ± 0.211.3 ± 0.2
RCP 4.5380 km 2144604.5 ± 0 km 2747 ± 2927.2 ± 1.212 ± 0
RCP 4.5c380 km 2144377 ± 0 km 2297 ± 628 ± 0.211.9 ± 0.1
RCP 8.5610 km 2107604.5 ± 0 km 2160 ± 029.1 ± 0.319.6 ± 0
RCP 8.5c610 km 2107604.5 ± 0 km 2150 ± 029.2 ± 0.119.7 ± 0
Appendix 2—figure 1
Download asset Open asset
Climate-analogue treeline position along the four transects (columns) based on extrapolated July temperatures from the relative concentration pathway (RCP) scenarios (rows).

A general rapid expansion of the treeline (continuous cover from south ≥1 tree ha−1, green) can be seen. In each panel, the maximum treeline position and the corresponding year are provided next to migration rates in three periods in km decade−1. Blue colour at the top represents the shoreline.

Appendix 2—figure 2
Download asset Open asset
Simulated forest expansion dynamics along the four transects (columns) forced with climate data based on relative concentration pathway (RCP) scenarios (rows).

A general quick expansion of single-tree stands (northernmost position of forest islands ahead of the treeline with ≥1 stem ha−1, light green) and the treeline (northernmost position of a continuous forest cover with 1 stem ha−1≤density<100 stem ha−1, green) followed by the forestline (northernmost position of a forest cover >100 stems ha−1, dark green) is seen. Blue colour at the top represents the shoreline.

Appendix 3

Treeline migration trajectories

Appendix 3—figure 1 shows the treeline trajectories of the position based on simulations and the climate-analogue position at the four transects and the time when the simulated positions reach the first time the same position of the climate-analogue.

Appendix 3—figure 1
Download asset Open asset
Trajectories for all four regions.

Numbers are the first year when the simulated treeline position is equal to or farther north than the modern climate-analogue position. Colour of line segments ranges from yellow for year 2000 to blue in 300 CE.

Appendix 4

Tundra area dynamics

Appendix 4—figure 1 gives an overview of the tundra area development in 100-year time steps for each of the eight climate scenarios.

Year

Appendix 4—figure 1
Download asset Open asset
Forest expansion for 100-year time steps and all scenarios.

Modern tundra areas will become covered by at least open larch tundra forests under different climate scenarios and nearly reach the shoreline in warmest conditions.

Data availability

We publicly uploaded the model code on Github https://github.com/StefanKruse/LAVESI/tree/crutransects (copy archived at swh:1:rev:422dc634a9d47e36c2e5c3a78906901ab6f2cdc3), and stored the commit used in this publication on Zenodo https://doi.org/10.5281/zenodo.6344261 Simulation results are stored on Zenodo https://doi.org/10.5281/zenodo.6484111.

The following data sets were generated
    1. Kruse S
    (2022) Zenodo
    Branch crutransects of LAVESI-WIND v1.2.
    https://doi.org/10.5281/zenodo.6344261
    1. Kruse S
    2. Herzschuh U
    (2022) Zenodo
    Forest expansion for different warming scenarios simulated for 2010 to 3000 CE with LAVESI for Siberia.
    https://doi.org/10.5281/zenodo.6484111
The following previously published data sets were used

References

  1. Book
    1. Abaimov AP
    (2010) Geographical Distribution and Genetics of Siberian Larch Species
    In: Osawa A, Zyryanova OA, Matsuura Y, Kajimoto T, Wein RW, editors. Permafrost Ecosystems - Siberian Larch Forests of Ecological Studies. Springer Netherlands. pp. 41–58.
    https://doi.org/10.1007/978-1-4020-9693-8_3
  2. Book
    1. Arctic Climate Impact Assessment
    (2004)
    Impacts of a Warming Arctic-Arctic Climate Impact Assessment
    Cambridge University Press.
  7. Book
    1. Bivand R
    2. Keitt T
    3. Rowlingson B
    (2019)
    Rgdal: Bindings for the ’Geospatial
    Data Abstraction Library.
    1. Callaghan TV
    2. Björn LO
    3. Chernov Y
    4. Chapin III F
    5. Christensen T
    6. Huntley B
    7. Ims R
    8. Johansson M
    9. Jolly D
    10. Matveyeva N
    11. Panikov N
    12. Oechel W
    13. Shaver G
    (2005)
    Arctic tundra and polar desert ecosystems
    Arctic Climate Impact Assessment 1:243–352.
  32. Book
    1. Krever V
    2. Stishov M
    3. Onufrenya I
    (2009)
    National Protected Areas of the Russian Federation: Gap Analysis and Perspective Framework
    Moscow: WWF.
  33. Book
    1. Kruse S
    2. Stoof-Leichsenring KR
    (2016) Keperveem - Past and present vegetation dynamics at the most eastern extension of the Siberian boreal treeline
    In: Overduin PP, Blender F, Bolshiyanov DY, Grigoriev MN, Morgenstern A, Meyer H, editors. Russian-German Cooperation: Expeditions to Siberia in 2016 (berichte z edition). Bremerhaven, Germany, Alfred-Wegener-Institut: Alfred-Wegener-Institut, Helmholtz-Zentrum für Polar- und Meeresforschung. pp. 1–709.
    https://doi.org/10.2312/BzPM_0709_2017
  36. Book
    1. Kruse S
    2. Herzschuh U
    3. Stünzi S
    4. Vyse S
    5. Zakharov E
    (2019b)
    Sampling mixed-species boreal forests affected by disturbances and mountain lake mountain lake and alas lake coring in Central Yakutia
    In: Kruse S, Bolshiyanov D, Grigoriev MN, Morgenstern A, Pestryakova L, Tsibizov L, Udke A, editors. Russian-German Cooperation: Expeditions to Siberia in 2018. Berichte Zur Polar-Und Meeresforschung = Reports on Polar and Marine Research. Bremerhaven: Alfred Wegener Institute for Polar and Marine Research. pp. 1–734.
  45. Report
    1. Masson-Delmotte V
    2. Zhai P
    3. Chen Y
    4. Goldfarb L
    5. Gomis MI
    6. Matthews JBR
    7. Berger S
    8. Huang M
    9. Yelekçi O
    10. Yu R
    11. Zhou B
    12. Lonnoy E
    13. Maycock TK
    14. Waterfield T
    15. Leitzell K
    16. Caud N
    (2021)
    Working Group I Contribution to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change
    IPCC.
  57. Book
    1. QGIS Development Team
    (2020)
    QGIS Geographic Information System
    QGIS Association.
  58. Software
    1. R Development Core Team
    (2019) R: A Language and Environment for Statistical Computing
    R Foundation for Statistical Computing, Vienna, Austria.
    http://www.r-project.org
  66. Conference
    1. Shevtsova I
    2. Herzschuh U
    3. Heim B
    4. Pestryakova LA
    5. Zakharov ES
    6. Kruse S
    (2021)
    Future spatially explicit tree above-ground biomass trajectories revealed for a mountainous treeline ecotone using the individual-based model LAVESI
    Environmental Research Letters.

Decision letter

  1. Meredith C Schuman
    Senior and Reviewing Editor; University of Zurich, Switzerland
  2. Roman Dial
    Reviewer; Alaska Pacific University, United States

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Regional opportunities for Siberian tundra conservation in the next 1000 years" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Meredith Schuman as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Roman Dial (Reviewer #1).

Essential revisions:

1) Systematically discuss in the introduction or the appendix the main limiting factors of tree establishment and growth relevant for the study area, and mentioning those finally implemented in the model would add considerable value (i.e. limiting factors that prevent tree establishment and growth, permafrost degradation, soil nutrient development, biotic interactions (herbivory)).

This discussion would increase traceability of methods and assessment of relevance of results, but also further emphasize how much this study is an improvement over previous studies by including some of these processes largely neglected earlier. Some of this very relevant information is mentioned in the response letter, but only partly introduced in the revised manuscript.

2) Discussion of limitations of the modelling study is largely missing, including the following aspects:

– Is this vegetation model coupled with a climate model? If not, feedbacks of forest expansion with climate and permafrost are currently neglected. The model is tested along gradients in selected regions, but it remains uncertain if space-for-time approach will hold in the future and further north, esp. when large-scale feedbacks are included.

– What about disturbances and extreme weather conditions that might regionally impact treeline advance or tree survival? E.g. increasing tundra fire activity might strongly impact vegetation development. Also, droughts/flooding might lead to regional vegetation impacts, esp. at seedling stage. Extreme events and related disturbances are predicted to increase under climate change and a discussion on how they might impact predictions by the model is needed. Are these factors all only short–term and neglectable compared to the long–term perspective modelled? If yes, this should be mentioned.

3) The Short Report format is appropriate for this study but does entail space limitations. Please relegate some description of machine computations to the supplementary material, and use the space for a fuller description in the methods of how warming, wind, and water are included in model parameterization (rather than citing a previous paper behind a pay–wall). Related: Figure 3 is both a computational and a conceptual graphic. Many readers might prefer to understand how climate change is incorporated into LAVESI conceptually, at least as much if not more so than how much computational time is required to run it.

– Methods: "Climate forcing" and "Tree growth on transects" needs at least some connection beyond citation of a previous paper.

– Line 205+: "July temperatures (Figure 4), which have the strongest impact in the model" Through what mechanisms? It would help a field biologist or plant physiologist to understand and evaluate the model.

– Line 209+: "Wind speed and direction data" How are these used in the model? What role do they play? A naïve reader would assume they play a role in dispersal if from the south and a role in mortality if from the north.

– Lines 259–272: These could be reduced and simplified to make room for the more valuable information about how the model parameters and structure are structurally related to July temperature and winds in winter or summer and to precipitation. I realize that these topics are likely detailed in previous publications, but a brief review or sketch even if put in another appendix would be invaluable to evaluating and understanding this very interesting model.

4) Clean–up

– In general, figures and graphs need to be checked and improved; expand legend, indicate axes units consistently, some figures and their legend not matching in appendix, check references to figures and tables in text, etc.

– Figure 1. Is the upper left panel of Appendix 2 Figure 1 the same as this and so redundant? The sub–panel in Figure 1 is tantalizing but not clear how best to interpret. Mention what RCP2.6* is standing for. hard to interpret graph: 4 scenarios mentioned, but 7 lines criss–crossing? x– and y–axis units? very hard to separate 25 years segments, is entire time–range until year 3000 displayed? What are the units of the axes? Does 0 = current forest boundary and 1 = Arctic Ocean coast? add more detail to legend, so it is self–explanatory. Like the diagonal is where climate and treeline are in equilibrium; below the diagonal tree migration lags climate (as it always will when isotherms move faster than dispersal); and above the diagonal is "overshooting" when there is some sort of momentum in the movement of treeline? Also, maybe arrows rather than colors would more readily indicate the direction of time?

– Figure 2: annotate x–axis of graphs around maps with years. This figure is too small to be well readable.

– Figure 3 "Data were aggregated for 100 years at three different population dynamics stages: e=empty, m=mature old growth, d=densification." Relevance of this caption and figure is unclear (other than the conceptual graphic which is good but could be excellent with the addition of where temperature, the main driving variable, fits in) to the paper about forest migration and tundra disappearance?

– Figure 4 (map) coming first would help orient Nearctic researchers better to the frequent reference to regional geography. Hard to distinguish the different RCPs versus treeline, centre, shoreline due to similar line symbols (maybe just use colors for 3 different locations within transect or RCP scenarios instead of different W–E locations) (4 transect locations can be indicated in the graphs with their name, no need to color code this information). RCP legend: include dots (e.g. 2.6 instead of 26) and name RCP1.3 RCP2.6* for consistency with text.

– Check that Table 1 is provided: the link to "Table 1" goes to Appendix 1 Table 1 which Reviewer #1 struggled to read. Do we need so many decimal places? Here in Appendix 1 Table 1, what exactly are the {plus minus} referring to? Why two columns at the right margin of the table that appear to overlap in their years of coverage?

– Use consistent units (km/decade is likely more appropriate than m/year).

– L33/34: 'The Siberian tundra areas are known for their alpine species biodiversity' – unclear why the arctic–alpine species are mentioned specifically. The Siberian tundra area covered by this study is important as it shows a high regional floristic diversity, representing almost all arctic tundra types. E.g. lowland tundras with their vegetation are very important for migratory bird species, so it might be misleading to highlight the arctic–alpine species diversity only rather than emphasizing the regional floristic diversity.

– Lines 38–39: "the modelled climate envelope for tundra migrates north by ~ 290 m yr-1 (Loarie et al., 2009)." Perhaps a more accurate statement from that paper would read, "the climate velocity of isotherms within the tundra biome have been calculated as ~ 290 m yr-1" although in that paper the word "velocity" seems loosely applied and climate speed would seem more appropriate except when it's claimed a direction (like up in elevation or north in latitude).

– Lines 39–40: "When assuming a direct relationship linking the modern position of the treeline to its main limiting factor of July temperature…" Do you mean this "When assuming a direct equilibrial relationship linking the modern position of the treeline to its main limiting factor of July temperature…"? The paper would benefit from more consistent and integrated use of terms.

– Lines 46–49: "This and the fact that in Siberia the treeline is formed solely by larch

species growing on permafrost while in North America this niche is covered by a mix of species (Mamet et al., 2019; Herzschuh, 2020), make it impossible to transfer knowledge from any other treeline region to inform Siberian tundra responses." Suggested word change, because nowhere is a conifer treeline moving rapidly, right?

"This and the fact that in Siberia the treeline is formed solely by larch species growing on permafrost while in North America this niche is covered by a mix of species (Mamet et al., 2019; Herzschuh, 2020), make it difficult if not impossible to transfer knowledge from any other treeline region to inform Siberian tundra responses."

– Lines 50–52: "Equilibrium–based estimations are not backed by empirical studies which are showing slower responses and either no treeline advance or rates of only a few metres per year in Siberia (Scherrer et al., 2020; Kharuk et al., 2013; Wieczorek et al., 2017b; Harsch et al., 2009)." Suggest replacement of Hofgaard et al., 2012 with Rees et al. 2020 as most up to date review of treeline speeds. Rees, W.G., Hofgaard, A., Boudreau, S., Cairns, D.M., Harper, K., Mamet, S., Mathisen, I., Swirad, Z. and Tutubalina, O., 2020. Is subarctic forest advance able to keep pace with climate change? Global Change Biology, 26(7), pp.3965–3977.

– L54: why is the focus on intra–specific competition here? why not interspecific competition across taxa, including insect, mammal herbivory etc?

– Lines 80–81: Would the authors consider introducing their results with a more believable biotic velocity (if I am understanding "treeline advance follows climate warming" correctly) than 3 km y–1 (i.e. 30 km decade–1). That seems very unreasonable!

Also: check figure and table references throughout – e.g. should be Appendix Table 1, and on L112 (see below) Figure 2B cannot be identified.

– L112: see above – you refer to Figure 2b, but no Figure 2b.

– Line 251+: “single–tree stands” should these be “0 tree ha–1 < density {less than or equal to} 1 tree ha–1” because otherwise as written “[density] > 1 stem (tree > 1.3 m tall) per ha” would include the other two classes since both treeline and forest line have “densities > 1 stem per ha”. Mathematically, “forest cover not falling below 1 stem ha−1” is stated as “forest cover {greater than or equal to} 1 stem ha−1” and “forest cover not falling below 100 stems ha−1” is stated as “forest cover {greater than or equal to} 100 stem ha−1”. These confusing descriptions show up again in Appendix 1 Figure 2. Caption.

– Line 255+: Also unclear is this statement: "The determined treeline at year 2000 is used as a baseline for expansion and subtracted from the following years' values."

– Is "climate analogue" defined somewhere? Is the idea related to "climate velocity"?

– Is "simulated treeline position" related to the idea of "biotic velocity"?

– L 476: is appendix Figure 1 showing treeline position or as mentioned in legend forestline position?

– Appendix 1 Figure 1. Vs Appendix 1 Figure 2: It's not obvious why these are showing different patterns when the captions are so similar: Is Appendix 1 Figure 1. "climate velocity" sensu Loarie et al. (2008) And Appendix 1 Figure 2 "biotic velocity" sensu Ordonez, A. and Williams, J.W., 2013. Climatic and biotic velocities for woody taxa distributions over the last 16 000 years in eastern North America. Ecology Letters, 16(6), pp.773–781.

– Appendix 1 Figure 1:"Forest expansion dynamics along the four transects (columns) extrapolated based on July temperatures (climate–analogue) for climate data based on RCP scenarios (rows). What do the negative values represent in migration rates? And why are the rates averaged beginning with 2000 instead of averaged over the particular intervals? Intervals would seem more informative: 2000–2100, 2100–2300, 2000–2500. indicate unit of y–axis (or provide/explain in legend) and unit of migration rate. Treeline position in each panel – is this the northernmost treeline position identified and the corresponding year?

– Appendix 1 Figure 2: "Forest expansion dynamics along the four transects (columns) forced with climate data based on RCP scenarios (rows)." is green not treeline and dark green forestline? These definitions are quite hidden and come very late on line 251 while they are very important!

– Appendix 2 Figure 1. Needs some text labeling each figure by transect. "Numbers are the first year when the simulated treeline position is equal to or farther north than the modern climate–analogue position." I see no numbers other than the RCP scenario identifiers. No years indicated, no blue color, seems to be different graph. Should also be moved to Appendix 2 section.

– Appendix 2 Figure 2. The caption references Figure 2 but that seems an error. If it is not, please tell the reader what exactly they should look for in Figure 2. no red in figure mentioned in legend.

– Finally, some references to other individual based models like Rupp, T.S., Chapin, F.S. and Starfield, A.M., 2001. Modeling the influence of topographic barriers on treeline advance at the forest–tundra ecotone in northwestern Alaska. Climatic change, 48(2), pp.399–416.

https://doi.org/10.7554/eLife.75163.sa1

Author response

Essential revisions:

1) Systematically discuss in the introduction or the appendix the main limiting factors of tree establishment and growth relevant for the study area, and mentioning those finally implemented in the model would add considerable value (i.e. limiting factors that prevent tree establishment and growth, permafrost degradation, soil nutrient development, biotic interactions (herbivory)).

This discussion would increase traceability of methods and assessment of relevance of results, but also further emphasize how much this study is an improvement over previous studies by including some of these processes largely neglected earlier. Some of this very relevant information is mentioned in the response letter, but only partly introduced in the revised manuscript.

Following the recommendations of the reviewer, we have included new paragraphs with a more detailed presentation of the implemented abiotic/biotic limiting factors in the Methods section. This includes either a statement of how these are explicitly considered, or argumentation as to whether these are implicitly considered and thus already part of the model LAVESI.

“Many factors influence the success of tree establishment and growth, of which the following seem most important to understand treeline migration. Warmth, especially summer mean temperatures, cumulative temperatures, and the vegetation period length, are key climate variables (e.g. MacDonald et al., 2008) and are hence included in the model for estimating potential tree growth and tree mortality (Kruse et al., 2016). This is further influenced by the abiotic (permafrost, active-layer depth) and biotic (intraspecific competition) environment. For example, permafrost soils with a shallow active-layer depth limit the space available for rooting and nutrient stores (e.g. Sullivan et al., 2015) and frozen/cold soils prevent growth but water can support growth in dry summers (Sugimoto et al., 2002). The latter functioning could be lost in future hot environments and in consequence it would decrease the treeline advance rate further and thereby increase the lag time. Accordingly, the active-layer depth is explicitly simulated in the model’s environment module, but soil nutrients are considered only implicitly. Individuals compete for space which serves as a surrogate for this. Permafrost degradation/abrupt thaw and warming (Biskaborn et al., 2019; Stuenzi et al., 2021) could lead to thermokarst processes that would locally cause swamping/waterlogging, leading to suppressed establishment/growth (Rees et al., 2020), which is not explicitly simulated and unlikely to limit tree invasion on a broader scale. Further limitations could arise from the lack of a suitable mycorrhiza symbiosis partner, but as shown by Hippel et al. (2021) even in northernmost areas on the Taimyr Peninsula mycorrhiza are present even before the invasion of trees in the Holocene. Additionally, dispersal rates of mycorrhiza with small spores are a magnitude faster than for taxa with larger seeds so that geographic distance is unlikely to be a limiting factor, similar to diatoms as shown by Stoof-Leichsenring et al. (2015). On the other hand, interspecific competition of recruits with present shrubs or grasses which are increasing more quickly – called the greening of the Arctic (Pearson et al., 2013; Berner et al., 2020) – is not explicitly simulated in LAVESI, although it is partly included in the germination probability and survival rate of individuals at the seedling stage. Further biotic interactions are implicitly included, such as herbivory by browsers (e.g. Wielgolaski et al., 2017), which reduces the survival and growth rate by damaging the shoots (compare krummholz study by Kruse et al., 2020) and can lead to a snow compaction cooling effect on soil temperatures (Beer et al., 2020), or pest outbreaks that are likely to increase in the future (e.g. larch silk moth, Fedotova, 2019). This would most likely lead to a further reduction of migration rates.

2) Discussion of limitations of the modelling study is largely missing, including the following aspects:

– Is this vegetation model coupled with a climate model? If not, feedbacks of forest expansion with climate and permafrost are currently neglected. The model is tested along gradients in selected regions, but it remains uncertain if space-for-time approach will hold in the future and further north, esp. when large-scale feedbacks are included.

We did not couple LAVESI to a climate model but as stated used only climate forcing, thus allowing the model to be run in stand-alone mode with data from different sources. Nonetheless, as presented in Pearson et al. (2013) tundra conversion to forests will ultimately lead to further warming feedback and hence a faster rate of conversion could be expected. We added a paragraph at the end of the Methods subsection ‘Tree growth on transects’: “The model used for this study is solely forced by climate data and no vegetation feedback on permafrost soils and climate was implemented. Therefore, the expected decrease of the albedo leading to a positive feedback with climate when tundra is colonized by tree stands (e.g. Pearson et al., 2013) is neglected here. Considering this, would likely further increase the rate of transition from tundra to forests.”

Regarding the space-for-time approach, we included a clarification and short discussion at the end of the Model section ‘Model tuning and validation’ we added the text: “We assume that the fundamental drivers of tree growth remain stable in time. The space-for-time approach used here is based on using several sites from across the treeline in the regions of interest for model parameterization and validation that have experienced a range of past climate (Pickett, 1989). However, climate warming will lead to unprecedented temperatures (Figure 1) that will affect the current stand structure and distribution. This may cause changes in the response of single considered processes or make other processes that are not currently important necessary (e.g. Blois et al., 2013; Pickett, 1989). Keeping this in mind, the future simulations should be interpreted with caution, although a simulation study with LAVESI revealed realistically predicted migration rates and patterns for the recent decades of strong warming in the study region of Chukotka (currently in review).”

– What about disturbances and extreme weather conditions that might regionally impact treeline advance or tree survival? E.g. increasing tundra fire activity might strongly impact vegetation development. Also, droughts/flooding might lead to regional vegetation impacts, esp. at seedling stage. Extreme events and related disturbances are predicted to increase under climate change and a discussion on how they might impact predictions by the model is needed. Are these factors all only short–term and neglectable compared to the long–term perspective modelled? If yes, this should be mentioned.

The reviewer brings attention to an important point. In our model, currently (tundra) fire is only implicitly considered, however drought impact on growth and survival is considered explicitly. On the other hand, flooding and waterlogging, that may take place locally, are not. The impact on predictions could be that the tundra colonization is even slower, hence, further prolonging the time-lagged response of the treeline migration. In the long run, an increased number of extreme weather events could support a faster dieback and supports a faster retreat of the treeline, which seems to overshoot and not return as quickly to its climate-analogue position. We added a paragraph at the end of the Methods section ‘Climate forcing’: “Future climate changes will lead to an increase in extreme weather events (Masson-Delmotte et al., 2021), among others a higher ignition probability and thus wildfires are expected to which the Siberian boreal forests are sensitive (e.g. Shvidenko and Schepaschenko, 2013). Fire is not yet explicitly simulated in the model and in a hotter environment fires could lead to dramatic forest losses. Areas might be turned into steppe and not revert to the prior tundra. Increases in extreme rainfall events will lead to local flooding and waterlogging (Masson-Delmotte et al., 2021), a response not explicitly implemented in the model. The anticipated effect to the LAVESI predictions is that tundra colonization would most likely be slowed down in the lowlands but less so in the well-drained mountainous regions, which would further prolong the time-lagged response of the treeline migration.”

3) The Short Report format is appropriate for this study but does entail space limitations. Please relegate some description of machine computations to the supplementary material, and use the space for a fuller description in the methods of how warming, wind, and water are included in model parameterization (rather than citing a previous paper behind a pay–wall). Related: Figure 3 is both a computational and a conceptual graphic. Many readers might prefer to understand how climate change is incorporated into LAVESI conceptually, at least as much if not more so than how much computational time is required to run it.

– Methods: "Climate forcing" and "Tree growth on transects" needs at least some connection beyond citation of a previous paper.

We added further relevant details on the tree growth estimation based on climate forcing data in the Methods section’s first subsection ‘Model description’ – please see the above comment for a detailed text addition and below for minor changes at relevant positions in the text.

– Line 205+: "July temperatures (Figure 4), which have the strongest impact in the model" Through what mechanisms? It would help a field biologist or plant physiologist to understand and evaluate the model.

Here the impact on tree growth is meant. In short, warmer summer temperatures lead to a potential larger tree growth (growth is positively correlated with temperature and the vegetation period length), however, this is constrained by competition. This was added to the Methods section in the first subsection ‘Model description’.

For clarification we adjusted the text at the reviewer comment’s position to “July temperatures (Figure 1), which have the strongest impact in the model due to a positive influence on tree diameter growth, increase in these scenarios …”

– Line 209+: "Wind speed and direction data" How are these used in the model? What role do they play? A naïve reader would assume they play a role in dispersal if from the south and a role in mortality if from the north.

For clarifying the use of wind data we included further details in the Methods section in the first subsection ‘Model description’. See changes made above.

– Lines 259–272: These could be reduced and simplified to make room for the more valuable information about how the model parameters and structure are structurally related to July temperature and winds in winter or summer and to precipitation. I realize that these topics are likely detailed in previous publications, but a brief review or sketch even if put in another appendix would be invaluable to evaluating and understanding this very interesting model.

We thank the reviewer for the critical assessment of our manuscript. Following their suggestions, we removed the model performance plots from figure 3 and merged them, together with the text on lines 162-172, into the appendix and refer to it in the updated text of the manuscript.

Further, we included more details about how the forcing climate data link to the internal processes in the Methods section in the subsection ‘Model description’ after the first paragraph of a general description. Our changes also include a new conceptual graphic in figure 4 that shows the relevant processes of LAVESI and how climate warming will impact the individual processes in LAVESI.

The new text with a detailed model description is: “The simulation proceeds in yearly time steps from the beginning to the end of the input climate time-series following a stabilization period to ensure that emerging populations reach equilibrium with the environment. Stochasticity in the model was introduced by using random numbers generated with a pseudo random number generator. Within one simulation year, the following processes become consecutively invoked Figure 4: Update of environment: Interactions between neighbouring trees are local and indirect. Basal diameters of each individual tree are used to evaluate the competition strength. We use a yearly updated density map to pass information about competition for resources between trees. An active-layer depth is estimated based on an edaphic factor and the number of days exceeding 0 °C following simplifications by Hinkel and Nicholas (1995). Growth: The maximum possible growth is estimated based on 10-yearsmean climate auxiliary variables (temperature of the coldest (January) and warmest (July) month, NDD0 ‘net degree days’, AAT10 ‘active air temperature’). The individual growth of basal diameter and, if a tree reached a height of 1.3 m, of breast height diameter, is calculated from the maximum possible growth in the current year affected by the tree’s density index. From the resulting diameters, the tree height is estimated. Seed dispersal: Seeds in ‘cones’ are dispersed from the parent trees. The dispersal directions and distances are randomly determined from a ballistic flight influenced by wind speed and direction with decreasing probabilities for long distances and, if dispersed seeds leave the extent of the simulated plot, they can be introduced from the other side only on the east-west margins. Seed production: Trees produce seeds after the year at which they reached their stochastically estimated maturation height. The total amount depends on weather, competition, and tree size. Establishment: The seeds that lie on the ground germinate at a rate depending on current weather conditions. Mortality: Individual trees or seeds die, i.e. they become removed from the plot, at a specified mortality rate. For trees this is deduced from long-term mean weather values, a drought index, surrounding tree density, tree age and size, plus a background mortality rate. Seeds on the other hand have the same constant mortality rate whether on trees or the ground. Ageing: Finally, the age of seeds and trees increases once a year and seeds are removed from the system when they reach a defined species age limit.”

Figure 4 was newly added:

4) Clean–up

– In general, figures and graphs need to be checked and improved; expand legend, indicate axes units consistently, some figures and their legend not matching in appendix, check references to figures and tables in text, etc.

We have edited and improved the figures in the manuscript. Please find our specific responses to the comments below.

– Figure 1. Is the upper left panel of Appendix 2 Figure 1 the same as this and so redundant? The sub–panel in Figure 1 is tantalizing but not clear how best to interpret. Mention what RCP2.6* is standing for. hard to interpret graph: 4 scenarios mentioned, but 7 lines criss–crossing? x– and y–axis units? very hard to separate 25 years segments, is entire time–range until year 3000 displayed? What are the units of the axes? Does 0 = current forest boundary and 1 = Arctic Ocean coast? add more detail to legend, so it is self–explanatory. Like the diagonal is where climate and treeline are in equilibrium; below the diagonal tree migration lags climate (as it always will when isotherms move faster than dispersal); and above the diagonal is "overshooting" when there is some sort of momentum in the movement of treeline? Also, maybe arrows rather than colors would more readily indicate the direction of time?

The figure presented in figure 1 is indeed the upper left panel in Appendix 2 Figure 1. We moved all four plots to Figure 2 in the main text and reduced the subpanel. Further, as requested, we added more detail to the caption. The figures show 8 individual paths, for each RCP scenario and the cooling and normal forcing, which was added for clarification in the caption as well. Additionally, we tested using arrows to indicate the direction and segment length and decided that this improves the figure further, which shows the full 3000 year range.

The caption was changed to: “The trajectory of the simulated treeline position relative to its current position and the maximum at the shoreline versus its climate-analogue position for the four regions shows a migration lag of the treeline during the first centuries (each line segment represents 25 years and the length of the arrow head corresponds to the step length) until the simulated treeline is limited by climate. Forests expand their area further and infilling proceeds when climate conditions cool and even overshoot in the long run with cooling back to 20th century temperatures. The diagonal is where climate and the treeline are in equilibrium; below the diagonal, tree migration lags climate; above the diagonal is "overshooting" and reaching locations actual climate would allow. For each RCP scenario two are presented, one for the scenario as-is and the second for the cooling; Scenario RCP 2.6* warms only at half the rate of RCP 2.6. See also Appendix 3 Figure 1 for the year when the trajectory passes the equilibrium.”

– Figure 2: annotate x–axis of graphs around maps with years. This figure is too small to be well readable.

We improved the figure by adding x- and y-axis labels to all graphs of tundra area dynamics and by increasing the size of the figure.

– Figure 3 "Data were aggregated for 100 years at three different population dynamics stages: e=empty, m=mature old growth, d=densification." Relevance of this caption and figure is unclear (other than the conceptual graphic which is good but could be excellent with the addition of where temperature, the main driving variable, fits in) to the paper about forest migration and tundra disappearance?

We removed this figure and its corresponding text from the main text and brought it to the Appendix as explained above in response to Essential Revisions 3.

– Figure 4 (map) coming first would help orient Nearctic researchers better to the frequent reference to regional geography. Hard to distinguish the different RCPs versus treeline, centre, shoreline due to similar line symbols (maybe just use colors for 3 different locations within transect or RCP scenarios instead of different W–E locations) (4 transect locations can be indicated in the graphs with their name, no need to color code this information). RCP legend: include dots (e.g. 2.6 instead of 26) and name RCP1.3 RCP2.6* for consistency with text.

We moved Figure 4 to the beginning and refer to it in the introduction. Further, we adjusted the colour to the position on the transect and added in the title the region name; changed RCP legend names, and corrected the typo in the caption.

– Check that Table 1 is provided: the link to "Table 1" goes to Appendix 1 Table 1 which Reviewer #1 struggled to read. Do we need so many decimal places? Here in Appendix 1 Table 1, what exactly are the {plus minus} referring to? Why two columns at the right margin of the table that appear to overlap in their years of coverage?

– Use consistent units (km/decade is likely more appropriate than m/year).

Table 1 is actually located in the appendix. Unfortunately, we wrongly referred to this in the text, which is now corrected. We left the numbers with one decimal place if not referring to years, which we rounded. The plus/minus refers to the standard deviation of the results that stem from three simulation repeats. The two different periods over which the mean migration rate was calculated are to compare those to each other.

– L33/34: 'The Siberian tundra areas are known for their alpine species biodiversity' – unclear why the arctic–alpine species are mentioned specifically. The Siberian tundra area covered by this study is important as it shows a high regional floristic diversity, representing almost all arctic tundra types. E.g. lowland tundras with their vegetation are very important for migratory bird species, so it might be misleading to highlight the arctic–alpine species diversity only rather than emphasizing the regional floristic diversity.

We changed the text to mention also the high regional floristic diversity and added a citation to the corresponding text in the introduction.

The text is now: “The Siberian tundra areas are known for their high regional floristic diversity and species biodiversity (Pauli et al., 2012; Mod and Luoto, 2016; Arctic Climate Impact Assessment, 2004; Schmidt et al., 2017) and support special (indigenous) types of land use (e.g. reindeer herding, foraging).”

– Lines 38–39: "the modelled climate envelope for tundra migrates north by ~ 290 m yr−1 (Loarie et al., 2009)." Perhaps a more accurate statement from that paper would read, "the climate velocity of isotherms within the tundra biome have been calculated as ~ 290 m yr−1 " although in that paper the word "velocity" seems loosely applied and climate speed would seem more appropriate except when it's claimed a direction (like up in elevation or north in latitude).

“In Siberia, the climate velocity of isotherms within the tundra biome has been calculated as ~ 290 m yr−1 (Loarie et al., 2009).”

– Lines 39–40: "When assuming a direct relationship linking the modern position of the treeline to its main limiting factor of July temperature…" Do you mean this "When assuming a direct equilibrial relationship linking the modern position of the treeline to its main limiting factor of July temperature…"? The paper would benefit from more consistent and integrated use of terms.

Following your recommendation, we changed the sentence here. Furthermore, for a more consistent use of the terms in the paper, we updated the corresponding statements in the Results and Discussion sections first part in the beginning. Also, in the caption of Figure 2 (trajectories) we added that the diagonal means that the treeline position is in equilibrium with climate, see comment above.

Text here is now “When assuming a direct equilibrial relationship linking the modern position of the treeline to its main limiting factor of July temperature (MacDonald et al., 2008), the treeline is expected to reach its maximum expansion within the first centuries of the millennium.”

Text in Result and Discussion is now: “The simulations revealed that the treeline is not in an equilibrial relationship but its advance follows climate warming, potentially reaching 30 km decade−1 in the 21st century under warmest climate scenarios (Appendix 2 Table 1).”

– Lines 46–49: "This and the fact that in Siberia the treeline is formed solely by larch

species growing on permafrost while in North America this niche is covered by a mix of species (Mamet et al., 2019; Herzschuh, 2020), make it impossible to transfer knowledge from any other treeline region to inform Siberian tundra responses." Suggested word change, because nowhere is a conifer treeline moving rapidly, right?

"This and the fact that in Siberia the treeline is formed solely by larch species growing on permafrost while in North America this niche is covered by a mix of species (Mamet et al., 2019; Herzschuh, 2020), make it difficult if not impossible to transfer knowledge from any other treeline region to inform Siberian tundra responses."

We follow the statement of the reviewer that current observations show only slow migrations of treelines worldwide. However, we demonstrate that they potentially pick up speed, even attaining exponential rates. Accordingly, we adjusted the sentence as suggested.

– Lines 50–52: "Equilibrium–based estimations are not backed by empirical studies which are showing slower responses and either no treeline advance or rates of only a few metres per year in Siberia (Scherrer et al., 2020; Kharuk et al., 2013; Wieczorek et al., 2017b; Harsch et al., 2009)." Suggest replacement of Hofgaard et al., 2012 with Rees et al. 2020 as most up to date review of treeline speeds. Rees, W.G., Hofgaard, A., Boudreau, S., Cairns, D.M., Harper, K., Mamet, S., Mathisen, I., Swirad, Z. and Tutubalina, O., 2020. Is subarctic forest advance able to keep pace with climate change? Global Change Biology, 26(7), pp.3965–3977.

We understand that the reviewer meant Harsch 2009 here and in the following sentence. We replaced it by the suggested publication that was already cited elsewhere in the manuscript.

– L54: why is the focus on intra–specific competition here? why not interspecific competition across taxa, including insect, mammal herbivory etc?

We added to the Methods a section (see comment above) covering further interspecific competition and disturbances affecting tree taxa dynamics by e.g. pest outbreaks and herbivory.

Additionally, we changed the text in Discussion section/Line 86ff to:

“Other processes, not yet explicitly captured by the model, may constrain the migration process, such as the strong greening of the Arctic increasing interspecific competition (compare Pearson et al., 2013; Berner et al., 2020), the warming and abrupt thaw of permafrost (e.g. Stuenzi et al., 2021; Biskaborn et al., 2019), soil nutrients (e.g. Sullivan et al., 2015), and animal activity (e.g. Wielgolaski et al., 2017).”

– Lines 80–81: Would the authors consider introducing their results with a more believable biotic velocity (if I am understanding "treeline advance follows climate warming" correctly) than 3 km y–1 (i.e. 30 km decade–1). That seems very unreasonable!

It seems counterintuitive at first that rates of 3 km yr-1 are realistic knowing the debate about post-glacial migration (e.g. Feardeau 2013). However, rates of 500 m yr-1 have likely occurred. Knowing that small seeds have a good long-distance dispersal ability in this landscape as well as forming nuclei ahead of the treeline will enhance the speed of colonization (e.g. Kruse et al. 2019) such that these high rates may be realistic. In our study they were only obtained for the very hot scenarios, which have higher growth rates and the production of more seeds, both factors that positively influence the migration rate. We rephrased the first misleading sentence that could be read that rates even higher than 3 km yr-1 may be simulated and added a reference and a short debate about this in the Results and Discussion section.

The text there is now: “The simulations revealed that the treeline is not in an equilibrial relationship but its advance follows climate warming, potentially reaching 30 km decade−1 in the 21st century under warmest climate scenarios (Appendix 2 Table 1).”

In the end of the paragraph: “the simulated rates are high compared to those during post-glacial migration (e.g. Feurdean et al., 2013) but seem likely as future climate is unprecedented, especially in their rate of change.”

Also: check figure and table references throughout – e.g. should be Appendix Table 1, and on L112 (see below) Figure 2B cannot be identified.

– L112: see above – you refer to Figure 2b, but no Figure 2b.

We checked the reference to the table, which as stated above was misleading. Figure 2b was not linked in the document and is now updated to refer to Figure 3 lower panel.

– Line 251+: “single–tree stands” should these be “0 tree ha–1 < density {less than or equal to} 1 tree ha–1” because otherwise as written “[density] > 1 stem (tree > 1.3 m tall) per ha” would include the other two classes since both treeline and forest line have “densities > 1 stem per ha”. Mathematically, “forest cover not falling below 1 stem ha−1” is stated as “forest cover {greater than or equal to} 1 stem ha−1” and “forest cover not falling below 100 stems ha−1” is stated as “forest cover {greater than or equal to} 100 stem ha−1”. These confusing descriptions show up again in Appendix 1 Figure 2. Caption.

The definitions at that location were not correct as information was missing about continuous tree cover as a prerequisite for the variables treeline and forestline. Hence, we edited the text for clarification to: “The positions of three key variables of stand densities expressed in stems, which are trees > 1.3 m tall, are extracted in 10-year steps: single-tree stands, defined as the northernmost position of forest islands ahead of the treeline with ≥ 1 stem ha−1; treeline stands, the northernmost position of a continuous forest cover with 1 stem ha−1 ≤ density 100 stem ha−1; and forestline stands, the northernmost position of forest cover > 100 stems ha−1 (see Figure 2 in Kruse et al., 2019a, for a graphical representation).”

And we also corrected the definitions in the captions of Appendix 1 Figures 1 and 2.

– Line 255+: Also unclear is this statement: "The determined treeline at year 2000 is used as a baseline for expansion and subtracted from the following years' values."

We clarified the statement to: “The position of the simulated treeline determined for year 2000 is used as the baseline for the calculation of migration rates.”

– Is "climate analogue" defined somewhere? Is the idea related to "climate velocity"?

– Is "simulated treeline position" related to the idea of "biotic velocity"?

We use it in the results and discussion first paragraph “When relating and tracking the position of the treeline to the climate (=climate-analogue) …” but missed a proper definition in the Methods section. This was added to the subsection Simulation setup and data processing: “For comparison, we extrapolated the position of the modern observed treeline based on July temperatures (=climate-analogue) for years 2000-3000 CE. Therefore, the current July temperature at the treeline position is determined and tracked along the transect at each time step of the climate forcing.”

About the “simulated treeline position”; Yes, this key variable and its migration rate is the biotic velocity, similar to the climate velocity.

– L 476: is appendix Figure 1 showing treeline position or as mentioned in legend forestline position?

Corrected as stated above.

– Appendix 1 Figure 1. Vs Appendix 1 Figure 2: It's not obvious why these are showing different patterns when the captions are so similar: Is Appendix 1 Figure 1. "climate velocity" sensu Loarie et al. (2008) And Appendix 1 Figure 2 "biotic velocity" sensu Ordonez, A. and Williams, J.W., 2013. Climatic and biotic velocities for woody taxa distributions over the last 16 000 years in eastern North America. Ecology Letters, 16(6), pp.773–781.

It is not obvious from the caption that Appendix 1 (now Appendix 2) Figure 2 shows the simulated forest expansion dynamics, so we adjusted the caption for clarification to “Simulated forest expansion dynamics along the four transects (columns) forced with climate data based on RCP scenarios (rows). A general quick expansion of single-tree stands (northernmost position of forest islands ahead of the treeline with ≥ 1 stem ha−1, light green) and the treeline (northernmost position of a continuous forest cover with 1 stem ha−1 ≤ density < 100 stem ha−1, green) followed by the forestline (northernmost position of a forest cover > 100 stems ha−1, dark green) is seen. Blue colour at the top represents the shoreline.” And also the caption of Appendix 1 (now Appendix 2) Figure 1 was changed to “Climate-analogue treeline position along the four transects (columns) based on extrapolated July temperatures from the RCP scenarios (rows). A general rapid expansion of the treeline (continuous cover from south ≥ 1 tree ha−1, green) can be seen. In each panel the maximum treeline position and the corresponding year is provided next to migration rates in three periods in km decade−1. Blue colour at the top represents the shoreline.”

We agree with the reviewer that the migration rates we present are velocities. The velocities along the transects were estimated both for the climate-analogue treeline position (see above comments) and extracted from the simulations. Accordingly, we cite the studies suggested by the reviewer and include for clarity those terms adjusted to the text elsewhere, as mentioned in above comments and in the methods section to “The resulting values for each of the three key variables are used in 10-year steps to interpolate with a weighted average the expansion starting at the current treeline (e.g. Walker et al., 2005). This is similar to the biotic velocity of Ordonez and Williams (2013). Opposed to this, we extrapolated the position of the modern observed treeline based on July temperatures (=climate-analogue) for the years 2000-3000 CE. Therefore, the current July temperature at the treeline position is determined and tracked along the transect at each time step of the climate forcing, which is similar to the climatic velocity of Loarie et al. (2009).”

– Appendix 1 Figure 1:"Forest expansion dynamics along the four transects (columns) extrapolated based on July temperatures (climate–analogue) for climate data based on RCP scenarios (rows). What do the negative values represent in migration rates? And why are the rates averaged beginning with 2000 instead of averaged over the particular intervals? Intervals would seem more informative: 2000–2100, 2100–2300, 2000–2500. indicate unit of y–axis (or provide/explain in legend) and unit of migration rate. Treeline position in each panel – is this the northernmost treeline position identified and the corresponding year?

Negative values represent a receding treeline and stem from calculating the mean over yearly migration rates in the time periods and ending at a more southerly position than at the start.

We selected intervals beginning in 2000 until 2100, 2300 or 2500 for the purpose of comparing whether, over longer time periods, the first strong migration is offset and returns back to a more southerly position.

The y-axis unit is provided for all plots at the left margin “Treeline position [km]”.

The unit of the migration rate was added in the caption along with the statement that the treeline position is stated at the top in each panel “In each panel the maximum treeline position and the corresponding year is provided next to migration rates in three periods in km decade-1.”

– Appendix 1 Figure 2: "Forest expansion dynamics along the four transects (columns) forced with climate data based on RCP scenarios (rows)." is green not treeline and dark green forestline? These definitions are quite hidden and come very late on line 251 while they are very important!

As stated above, we have corrected the definitions in the Methods section and in the captions.

– Appendix 2 Figure 1. Needs some text labeling each figure by transect. "Numbers are the first year when the simulated treeline position is equal to or farther north than the modern climate–analogue position." I see no numbers other than the RCP scenario identifiers. No years indicated, no blue color, seems to be different graph. Should also be moved to Appendix 2 section.

– Appendix 2 Figure 2. The caption references Figure 2 but that seems an error. If it is not, please tell the reader what exactly they should look for in Figure 2. no red in figure mentioned in legend.

The figures and its legend were unfortunately misplaced and the caption for Appendix 2 Figure 1 referred to Appendix 2 Figure 2. We moved the panel plot of all four regions to the main part of the manuscript. In this new graph we added the name of the regions to each header of the individual plots.

– Finally, some references to other individual based models like Rupp, T.S., Chapin, F.S. and Starfield, A.M., 2001. Modeling the influence of topographic barriers on treeline advance at the forest–tundra ecotone in northwestern Alaska. Climatic change, 48(2), pp.399–416.

We added the citation in the first paragraph in the Introduction as it is an early simulation study similar to our setup, however it is not an individual based model, but is spatially-explicit in the sense that for seed dispersal 2x2 km grid cells are considered, and on each grid cell, processes capture the turnover of ecotypes, which is quite different from our approach of simulating each single tree from seedling stage onwards and from purely spatially-explicit processes such as seed dispersal and competition etc.

https://doi.org/10.7554/eLife.75163.sa2

Article and author information

Author details

  1. Stefan Kruse

    Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research, Potsdam, Germany
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Software, Validation, Visualization, Writing – original draft, Writing – review and editing
    For correspondence
    stefan.kruse@awi.de
    Competing interests
    The authors declare no competing interests
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1107-1958

  2. Ulrike Herzschuh

    1. Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research, Potsdam, Germany
    2. Institute of Environmental Sciences 6 and Geography, University of Potsdam, Potsdam-Golm, Germany
    3. Institute of 7 Biochemistry and Biology, University of Potsdam, Potsdam-Golm, Germany
    Contribution
    Conceptualization, Funding acquisition, Methodology, Project administration, Supervision, Visualization, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared

Funding

H2020 European Research Council (772852)

  • Stefan Kruse
  • Ulrike Herzschuh

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This study acknowledges support by the ERC consolidator grant (no. 772852) led by Ulrike Herzschuh. We like to thank Sven Willner for assisting in programming and OpenMP improvements and Cathy Jenks for proofreading and improving the article. Further, we thank two reviewers and the topical handling editor for their comments which improved the article.

Senior and Reviewing Editor

  1. Meredith C Schuman, University of Zurich, Switzerland

Reviewer

  1. Roman Dial, Alaska Pacific University, United States

Publication history

  1. Received: November 1, 2021
  2. Preprint posted: November 30, 2021 (view preprint)
  3. Accepted: April 26, 2022
  4. Version of Record published: May 24, 2022 (version 1)

Copyright

© 2022, Kruse and Herzschuh

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,194
    Page views
  • 353
    Downloads
  • 4
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Stefan Kruse
  2. Ulrike Herzschuh
(2022)
Regional opportunities for tundra conservation in the next 1000 years
eLife 11:e75163.
https://doi.org/10.7554/eLife.75163

Categories and tags

Further reading

    1. Ecology
    2. Evolutionary Biology

    Tradeoffs explain scaling, sex differences, and seasonal oscillations in the remarkable weapons of snapping shrimp (Alpheus spp.)

    Jason P Dinh, SN Patek
    Research Article Updated

    Evolutionary theory suggests that individuals should express costly traits at a magnitude that optimizes the trait bearer’s cost-benefit difference. Trait expression varies across a species because costs and benefits vary among individuals. For example, if large individuals pay lower costs than small individuals, then larger individuals should reach optimal cost-benefit differences at greater trait magnitudes. Using the cavitation-shooting weapons found in the big claws of male and female snapping shrimp, we test whether size- and sex-dependent expenditures explain scaling and sex differences in weapon size. We found that males and females from three snapping shrimp species (Alpheus heterochaelis, Alpheus angulosus, and Alpheus estuariensis) show patterns consistent with tradeoffs between weapon and abdomen size. For male A. heterochaelis, the species for which we had the greatest statistical power, smaller individuals showed steeper tradeoffs. Our extensive dataset in A. heterochaelis also included data about pairing, breeding season, and egg clutch size. Therefore, we could test for reproductive tradeoffs and benefits in this species. Female A. heterochaelis exhibited tradeoffs between weapon size and egg count, average egg volume, and total egg mass volume. For average egg volume, smaller females exhibited steeper tradeoffs. Furthermore, in males but not females, large weapons were positively correlated with the probability of being paired and the relative size of their pair mates. In conclusion, we identified size-dependent tradeoffs that could underlie reliable scaling of costly traits. Furthermore, weapons are especially beneficial to males and burdensome to females, which could explain why males have larger weapons than females.

    1. Ecology
    2. Evolutionary Biology

    Thermal phenotypic plasticity of pre- and post-copulatory male harm buffers sexual conflict in wild Drosophila melanogaster

    Claudia Londoño-Nieto, Roberto García-Roa ... Pau Carazo
    Research Article Updated

    Strong sexual selection frequently leads to sexual conflict and ensuing male harm, whereby males increase their reproductive success at the expense of harming females. Male harm is a widespread evolutionary phenomenon with a strong bearing on population viability. Thus, understanding how it unfolds in the wild is a current priority. Here, we sampled a wild Drosophila melanogaster population and studied male harm across the normal range of temperatures under which it reproduces optimally in nature by comparing female lifetime reproductive success and underlying male harm mechanisms under monogamy (i.e. low male competition/harm) vs. polyandry (i.e. high male competition/harm). While females had equal lifetime reproductive success across temperatures under monogamy, polyandry resulted in a maximum decrease of female fitness at 24°C (35%), reducing its impact at both 20°C (22%), and 28°C (10%). Furthermore, female fitness components and pre- (i.e. harassment) and post-copulatory (i.e. ejaculate toxicity) mechanisms of male harm were asymmetrically affected by temperature. At 20°C, male harassment of females was reduced, and polyandry accelerated female actuarial aging. In contrast, the effect of mating on female receptivity (a component of ejaculate toxicity) was affected at 28°C, where the mating costs for females decreased and polyandry mostly resulted in accelerated reproductive aging. We thus show that, across a natural thermal range, sexual conflict processes and their effects on female fitness components are plastic and complex. As a result, the net effect of male harm on overall population viability is likely to be lower than previously surmised. We discuss how such plasticity may affect selection, adaptation and, ultimately, evolutionary rescue under a warming climate.

    1. Ecology
    2. Evolutionary Biology

    Sensory conflict disrupts circadian rhythms in the sea anemone Nematostella vectensis

    Cory A Berger, Ann M Tarrant
    Research Article Updated

    Circadian clocks infer time of day by integrating information from cyclic environmental factors called zeitgebers, including light and temperature. Single zeitgebers entrain circadian rhythms, but few studies have addressed how multiple, simultaneous zeitgeber cycles interact to affect clock behavior. Misalignment between zeitgebers (‘sensory conflict’) can disrupt circadian rhythms, or alternatively clocks may privilege information from one zeitgeber over another. Here, we show that temperature cycles modulate circadian locomotor rhythms in Nematostella vectensis, a model system for cnidarian circadian biology. We conduct behavioral experiments across a comprehensive range of light and temperature cycles and find that Nematostella’s circadian behavior is disrupted by chronic misalignment between light and temperature, which involves disruption of the endogenous clock itself rather than a simple masking effect. Sensory conflict also disrupts the rhythmic transcriptome, with numerous genes losing rhythmic expression. However, many metabolic genes remained rhythmic and in-phase with temperature, and other genes even gained rhythmicity, implying that some rhythmic metabolic processes persist even when behavior is disrupted. Our results show that a cnidarian clock relies on information from light and temperature, rather than prioritizing one signal over the other. Although we identify limits to the clock’s ability to integrate conflicting sensory information, there is also a surprising robustness of behavioral and transcriptional rhythmicity.