Treeline dynamics in Siberia under changing climates as inferred from an individual-based model for Larix

Ecological Modelling 338 (2016) 101–121

Contents lists available at ScienceDirect

Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel

Treeline dynamics in Siberia under changing climates as inferred from an individual-based model for Larix Stefan Kruse a,∗ , Mareike Wieczorek a , Florian Jeltsch b,d,e , Ulrike Herzschuh a,c a

Department of Periglacial Research, Alfred-Wegener-Institute, Helmholtz Centre for Polar and Marine Research (AWI), D-14473 Potsdam, Germany Institute of Biochemistry and Biology, University of Potsdam, D-14476 Potsdam-Golm, Germany Institute of Earth and Environmental Science, University of Potsdam, D-14476 Potsdam-Golm, Germany d ZALF, Leibniz-Centre for Agricultural Landscape Research, Eberswalder Str. 84, D-15374 Müncheberg, Germany e Berlin-Brandenburg Institute of Advanced Biodiversity Research (BBIB), D-14195 Berlin, Germany b c

a r t i c l e

i n f o

Article history: Received 9 March 2016 Received in revised form 4 August 2016 Accepted 5 August 2016 Keywords: Forest change IBM ODD model description Larix gmelinii Permafrost ecosystem Time-lag effects

a b s t r a c t Siberian boreal forests are expected to expand northwards in the course of global warming. However, processes of the treeline ecotone transition, as well astiming and related climate feedbacks are still not understood. Here, we present ‘Larix Vegetation Simulator’ LAVESI, an individual-based spatially-explicit model that can simulate Larix gmelinii (Rupr.) Rupr. stand dynamics in an attempt to improve our understanding about past and future treeline movements under changing climates. The relevant processes (growth, seed production and dispersal, establishment and mortality) are incorporated and adjusted to observation data mainly gained from the literature. Results of a local sensitivity analysis support the robustness of the model’s parameterization by giving relatively small sensitivity values. We tested the model by simulating tree stands under modern climate across the whole Taymyr Peninsula, north-central Siberia (c. 64–80◦ N; 92–119◦ E). We find tree densities similar to observed forests in the northern to midtreeline areas, but densities are overestimated in the southern parts of the simulated region. Finally, from a temperature-forcing experiment, we detect that the responses of tree stands lag the hypothetical warming by several decades, until the end of 21st century. With our simulation experiments we demonstrate that the newly-developed model captures the dynamics of the Siberian latitudinal treeline. © 2016 Published by Elsevier B.V.

1. Introduction Taxa forming the boreal treeline (Holtmeier and Broll, 2005) are expected to expand their distribution area northwards under a future warmer climate (Pearson et al., 2013). This tundra-totaiga transition will be accompanied by an albedo decrease (Bonan, 2008), which in turn will raise regional temperatures. Such positive feedback can significantly contribute to global warming as revealed by Earth-system modelling studies (Stocker et al., 2013). However, the detailed tree population processes that ultimately determine the timing and magnitude of tree responses to a given climate signal are still not fully understood. For example, seed dispersal limitations and slow reproduction rates will lead to leading-edge vegetation-climate disequilibrium that could last for decades or centuries (Lenoir and Svenning, 2015). Furthermore, one should

∗ Corresponding author at: Alfred-Wegener-Institute, Helmholtz Centre for Polar and Marine Research (AWI), Research Unit Potsdam, Telegrafenberg A43, D-14473 Potsdam, Germany. E-mail address: [email protected] (S. Kruse). http://dx.doi.org/10.1016/j.ecolmodel.2016.08.003 0304-3800/© 2016 Published by Elsevier B.V.

keep in mind that a limited or delayed reaction to the temperature signal is also possible because increased temperatures could lead to decreased growth rates due to drought stress generally in more southerly treeline stands or on slopes (Lawrence et al., 2015; Barber et al., 2000). This can, however, be buffered by permafrost meltwater that is an important direct source of moisture for larches (Sugimoto et al., 2002), but if this meltwater accumulates in depressions it can hinder growth due to water stress (Lawrence et al., 2015; Kharuk et al., 2015). Additionally, fires play a role in larch population dynamics in the northern parts of central Siberia (Sofronov et al., 2000), with frequency decreasing from high fire return intervals of around 100 years in the south (61◦ N) to low intervals of 350 years in the northernmost forests (Kharuk et al., 2011). Hitherto, time-series of population dynamics from treeline areas of reasonable temporal length are lacking. Vegetation-plot studies along treeline transects can portray tree-stand characteristics and recruitment patterns under different climates and thus can provide hints about future changes. However, the complexity of the treeline response to climate can best be disentangled with

102

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

the help of numerical modelling. Among the various vegetation model types available, an individual-based spatially-explicit model could be most suitable for simulating vegetation dynamics at the treeline. Such models can capture the most relevant processes with respect to the establishment, structure, and expansion of tree stands, including individual tree ageing, competition among different life stages, seed dispersal, and seed-bank establishment. For example, Martínez et al. (2011) successfully applied an individualbased spatially-explicit model of alpine treelines to evaluate the importance of krummholz for treeline range expansion. Treelines are known as tundra-taiga ecotones (Holtmeier and Broll, 2005; MacDonald et al., 2008) that are often defined by tree heights of between 2 m and 8 m at the latitudinal upper and lower tree limit (Holtmeier, 2009). We use the term treeline here for a region spanning from the limits of dense northern taiga forests in favourable habitats to trees that reach a maximum height of 2 m due to deteriorating conditions. The vast Siberian treeline is unique with respect to its sole composition of Larix Mill., a deciduous needle-leaf tree genus (Abaimov, 2010; Franz, 1973b). Evidence from other treeline areas cannot, therefore, simply be applied to Siberia. There are three Larix species forming the treeline up to ∼72.5◦ N: Larix sibirica Ledeb. from roughly 60–90◦ E, Larix gmelinii (Rupr.) Rupr. towards the east where continuous permafrost prevails between 90 and 120◦ E, and Larix cajanderi Mayr. from 120 to 160◦ E. These species are well adapted to the harsh north Siberian permafrost environment experiencing cold winters up to −40 ◦ C and short summers with only 60–90 days exceeding the freezing point. The regional temperature regime causes a maximum summer active-layer depth ranging only between 0.3–0.6 m in some areas (Franz, 1973a). If the climate envelope shifts north, forests in Siberia will probably densify and expand northwards. Such a general relationship is corroborated by results from palaeoecological investigations of glacial/interglacial and Holocene treeline transitions (Andreev et al., 2011; MacDonald et al., 2008) and the few available field observations of modern tree-population processes at the treeline (Golubeva et al., 2013; Kharuk et al., 2006). Computer simulations that could help to understand the underlying processes at the Siberian treeline are rare. A rather detailed approach for northern boreal forests was recently published by Brazhnik and Shugart (2015, 2016): their forest gap model SIBBORK includes a 3-D light environment and allows for complex interactions between boreal forest species and their environment. Simulations suggest a 2 ◦ C increase in annual temperature will cause northwards range expansion of larch forests and steppe formation in the south (Brazhnik and Shugart, 2015). Zhang et al. (2013) applied the LPJ-GUESS in cohort mode to simulate high-latitude arctic vegetation. Their results suggest that the treeline will advance polewards in the future and be accompanied by reduced albedo in winter months due to strong densification. Similarly, the application of the ArcVeg model to the Yamal Peninsula (Yu et al., 2009) led to the conclusion that a significant increase in biomass will occur when climate warms by 2 ◦ C. Despite their detailed implementation of abiotic processes, such as radiation, water table level and available nutrients, which affect growth, and the computation of these on daily time steps, such models are optimized for a global application and as a consequence the life-cycle stages are highly generalized (Sato et al., 2007; see, for example, Sitch et al., 2003). These models might therefore miss important details of relevance to local stand dynamics, which could enhance or diminish the pressure of climate change via biotic interactions (e.g. seedling establishment under unfavourable abiotic conditions triggered by facilitation, Martínez et al., 2011). To study complex ecosystems where individual-based processes cause observed patterns, individual-based models are widely used in ecology (DeAngelis and Mooij, 2005; Grimm and Railsback, 2005; Jeltsch et al., 2008; Martínez et al., 2011; Peck,

2004; Seidl et al., 2012), but have not been available for larch forests growing on permafrost. As a consequence, we set up a rather simple model framework and successively implemented processes which generate patterns such as those observed at reference sites, until we created a final model that is able to capture the observed patterns. With our focus on understanding local stand interactions that can cause larger scale population dynamics, seeds and tree individuals are handled individually and have explicit coordinates. This model allows for spatially-explicit interactions among all life-history phases of larches and is therefore a handy tool to realistically simulate responses of population dynamics in a changing environment (past-to-present or present-to-future). Additionally, information can be inherited from one parent individual to the offspring, so a fruitful future study could be to implement different genetic markers and validate the model against results of population genetics that try to infer the phylogeography of a species (for example central and far-east Siberian species: Semerikov et al., 2007; Polezhaeva et al., 2010). Here we present an individual-based model (LAVESI = Larix Vegetation Simulator) to simulate the relevant processes of Larix gmelinii population dynamics at the arctic treeline in Siberia. First, we calibrated LAVESI and checked the influence of parameterization on model performance by conducting a sensitivity analysis. Second, we simulated tree stands forming the modern treeline area on the southern Taymyr Peninsula, north-central Siberia, emphasizing stand density in our analysis. And, finally, we investigated the response of population density and structure of tree stands to a variety of hypothetical warmer and colder climate scenarios. 2. Material and methods 2.1. Reference sites Field data were collected from two replicate plots at each of five sites. The five sites used in this study were visited during the vegetation periods of the summers of 2011 and 2013 (M. Wieczorek et al., submitted; Table 1, Fig. 1). Locations are situated along a southwest to northeast transect across the treeline and are representative of tree growth at the treeline, which passes through the Khatanga River flood basin at one of the most northerly locations for tree growth at around 73◦ N. At each site, larch individuals in two plots of at least 20 × 20 m were examined; variables such as height, basal diameter and diameter at breast height were measured, number of cones estimated and the position of trees recorded. 2.2. Description of the model LAVESI Following the approach of pattern-oriented modelling (Grimm and Railsback, 2005) focusing on forest stand densities of several reference sites, we started with a simple growth model and successively added functions to allow the model to simulate those observations. Based on literature about larches growing at the treeline in Siberia (e.g. Osawa and Kajimoto, 2010), we designed the basic framework to represent the life-cycle of different species of larches as completely as possible from seeds to mature trees. This version was adapted for the specific deciduous larch Larix gmelinii (Rupr.) Rupr. that grows in single species stands on prevailing permafrost and under harsh cold climates in northern central Siberia (Naurzbaev et al., 2002). Dendrochronological analyses show that tree growth in these areas mainly depends on temperature, especially summer temperature (Abaimov, 2010), and to some extent on precipitation (Sidorova et al., 2010). Thus, we coupled the simulations to monthly temperature and precipitation series as driving forces. Within one simulation step, which is one year, a variety of processes become consecutively invoked: growth, seed produc-

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

103

Table 1 Reference site specifications. Site

Latitude [◦ N]

Longitude [◦ E]

Mean annual temperature 1960–1990 [◦ C]

Mean annual precipitation 1960–1990 [mm]

Vegetation zone

Reference values for each plot [trees ha−1 ]

1 2 3 4 5

72.4 72.2 71.8 71.1 70.7

105.4 102.2 102.9 100.8 97.7

−13.95 −12.94 −13.11 −13.69 −13.37

254.9 265.6 261.6 293.9 350.9

Southern Tundra Forest Tundra Forest Tundra Forest Tundra Northern Taiga

200–425 725–1062 1150–1250 1163–1550 700–1100

Fig. 1. Location of reference sites (1–5) on the southern Taymyr Peninsula, north-central Siberia. The map shows the visited sites 1–5 in black triangles; the red dotted line represents the 12 ◦ C July isotherm for 1950–1990 and the blue dash-dotted line the 250 mm isohyet for the same period [Source: CRU TS 3.22 (Harris et al., 2014)]; the brown dashed line marks the recent northern tree limit position [Source: CAVM (Walker et al., 2005)]; grey areas >150 m a.s.l. [Source: ETOPO-5 (National Geophysical Data Center, 1988)]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

tion and dispersal, establishment and mortality (Fig. 2; see ‘2.2.2 Description of sub-models’). The model was implemented in C++ using standard template libraries and compilation. Model runs were performed using a standard computer. The program code is available from the authors on request.

2.2.1. The ODD-Protocol for LAVESI The model description follows the ODD protocol (named after the three sections: Overview, Design Concepts, and Details) for describing individual- and agent-based models (Grimm and Railsback, 2005; Grimm et al., 2006; using the updated version from Grimm et al., 2010), with the ‘design concepts’ stated in ‘2.2.1.3 Process overview and scheduling’.

2.2.1.1. Purpose. The model was set up to understand tree stand structure and the dynamics of Larix gmelinii (Rupr.) Rupr. popu-

lations growing in the Siberian treeline ecotone in response to a changing climate. 2.2.1.2. Entities, state variables and scales. The model consists of two hierarchical levels characterized by a set of variables (Table 2): (1) simulation areas characterized by the specific biotic and abiotic environment, and (2) individual trees and seeds. Each individual square simulation area covers 100 × 100 m on which seeds and trees are exactly positioned by x,y coordinates. Using the basal diameter of individual trees, the plot is overlaid with a tree density grid with a resolution of 0.2 × 0.2 m. Of the whole simulation area, the central 20 × 20 m represents the investigation plot ensuring a border of 40 m to the boundaries to minimize potential boundary effects. We tested a set of varying individual square simulation areas, from 25 m edge length to 500 m, and evaluated the strength of the boundary effect on the central 20 × 20 m plot within the simulated tree stand in comparison to a 100 × 100 m area size.

104

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

Fig. 2. Flow chart of the model’s processes for the simulation cycle of the LAVESI model which proceeds in yearly time steps after initialization. Plots are created, model parameters and variables set, and weather data is read and processed for each simulation year. Within each simulation cycle, sub-models are consecutively invoked and finally the biotic environment is updated prior to the next simulation cycle. Arrows show feedback between sub-models within one simulation year.

Table 2 State variables of the model with dimensions and update rates ordered by the model’s hierarchy levels. Hierarchy Level 1 (100 × 100 m) Biotic environment Tree density Abiotic environment Monthly mean temperatures Monthly precipitation sums Isotherms for January and July ‘net degree days’ NDD0 ‘active air temperature’ AAT10 Precipitation sum Maximum basal diameter growth and at breast height Level 2 (individuals) Trees Position (x,y coordinates) Year of establishment Diameter at basal and breast height Relative diameter growth at basal and breast height List of diameters at breast height Height Age Cones (yes/no) Height to bear cones Density index Seeds Position (x,y coordinates) Location (cone/ground) Age Long dispersed seed (yes/no)

Dimension

Update ratea

(0.2 × 0.2 m)−1

yearly

◦

C mm ◦ C – ◦ C mm cm

* * * * * * *

m AD cm – cm Cm years – Cm –

* * yearly yearly yearly yearly yearly yearly after maturing yearly

M – AD –

* * and after dispersal event yearly after dispersal event

a Asterisks mean updated once at establishment (trees), production (seeds) or initialization (abiotic environment); NDD0 number of days exceeding 0 ◦ C and AAT10 sum of temperatures above 10 ◦ C.

The decision for the simulation area extent was based on an analysis where we used the similarity between predefined index variables calculated for results of 30 simulation runs on different area sizes. The simulation-to-simulation similarity between the results was

converging with increasing area and reached a plateau with high similarity at 75 m edge length (not shown). The results therefore let us assume an optimal simulation area for the current version of the model of 100 m edge length to minimize boundary effects. Nev-

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

105

Fig. 3. For one of the simulated locations (site 2, a) the real climate data is shown next to the randomized weather used in the stabilization phase of the model; and (b) the corresponding simulated population variables ‘tree diameter growth’ (green) and ‘stemcount’ (brown) are shown for this climate. Solid lines are 11-year running mean values. The lower panels show trees (>1.3 m, in green) growing in the central 20 × 20 m plot and the tree density grid (in brown colours; mean value and standard deviation above the plot) at three stages in the course of the population development throughout the simulation period: (c) after 100 years first trees grow out of the seedling stage, (d) after 1000 years the population reaches its equilibrium under the forced climate conditions, and (e) the population at the final step of simulation, which is used for evaluating the model. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

ertheless, there might be other factors that could penetrate more than 100 m into a forest (Laurance et al., 2002), such as wind, and if attempting to incorporate such influences explicitly whilst simulating forest patches in an inhomogeneous environment, the edge impact analysis must be revisited. Simulation runs proceed in yearly time steps. We performed simulations for years 1919–2011, where robust climate series were available. Additionally, to reach stabilization of population dynamics and the forcing climate series, simulations were preceded by a stabilization period with a length of 1000 years (for parameterization, sensitivity analysis and Taymyr treeline application) or 5001 years (for temperature experiment). All simulations start from bare ground introducing 1000 seeds in the first 100 years and, to allow for repopulation of simulation areas after extinction, 10 seeds are added every year to the simulation areas.

2.2.1.3. Process overview and scheduling. The simulation proceeds in yearly time steps from the beginning to the end of the input climate time-series following a stabilization period of 1000 years

to ensure that emerging populations reach equilibrium with the environment. In each initialization phase of each simulation run, the weather data are processed and used to estimate maximum diameter growth (at basal and breast height) for each simulation year based on 10-years mean climate auxiliary variables (see details in ‘2.2.2 Description of sub-models’). Within the growth processes of the model, these variables are used to individually estimate the current diameter growth of trees constrained by their actual biotic environment (Design concept: Sensing). Stochasticity in the model was introduced by using random numbers generated with a pseudo random number generator (C++-function ‘rand’, using the function ‘srand’ for seeding and using a runtime value with the function call ‘time(0)’ to allow for different results between two or more consecutive runs of the model; Design Concept: Stochasticity). Within one simulation year, the following processes become consecutively invoked (Fig. 2; detailed explanations for each process can be found in a corresponding section in ‘2.2.2 Description of sub-models’, below):

106

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

Fig. 4. Sensitivity values (boxplots, left y-axis) and probabilities that a parameter change results in simulations which depart from the reference runs (lines, right y-axis) calculated for simulations at reference site 3 for two different development phases: exponential (left panels a & b at year 500) and equilibrium (right panels c & d at year 1083, i.e. 2011 AD). Values sorted from left-to-right by decreasing mean probability for 30 realizations with parameter value changes of +5% (a & c) and −5% (b & d).

2.2.1.3.1. Update of environment. Interactions between neighbouring trees are local and indirect. Basal diameters of each individual tree are used to evaluate the competition strength. We use a yearly updated density map to pass information about competition for resources between trees. (Design Concept: Interaction). 2.2.1.3.2. Growth. The individual growth of basal diameter and, if a tree reached a height of 1.3 m, of breast height diameter, is calculated from the maximum possible growth in the current year affected by the tree’s density index. From the resulting diameters, the tree height is estimated differently for the two height classes, smaller and greater than 1.3 m. (Design Concept: Collectives). 2.2.1.3.3. Seed dispersal. Seeds in ‘cones’ are dispersed from the parent trees, at a set rate. The dispersal directions are randomly determined with decreasing probabilities for long distances and, if dispersed seeds leave the extent of the simulated plot, they are removed from the system. 2.2.1.3.4. 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. 2.2.1.3.5. Establishment. The seeds that lie on the ground germinate at a rate depending on current weather conditions. 2.2.1.3.6. 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 and or the ground. (Design Concept: Emergence).

2.2.1.3.7. Ageing. Finally, the age of seeds and trees increases once a year and seeds are removed from the system when they reach a defined species age limit. 2.2.1.4. Initialization. The first step of the model is to read the model parameters and variables and the weather data. Thereafter, the determined amount of empty 100 × 100 m simulation areas is invoked and the density grid with 0.2 × 0.2 m tiles is laid over the whole area. Each density tile in this grid starts with an empty tree number field and a tree density value of one. In the simulation area, the following two life stages of Larix occur: (i) seeds, with their plots’ positions and exact location (i.e. in ‘cones’ on trees or on the ground), start with ages of zero years, and (ii) if a tree is established it obtains its plot’s position and exact location from the seed. Initial basal diameters equal the current maximum basal diameter growth rate, calculated from current weather at the tree’s position, and heights are calculated by height-diameter functions. Diameter at breast height, age and number of yearly produced seeds as well as tree density indices (which are defined in the sub-models section below), all start with a value of zero. Growth indices (see definition in ‘2.2.2. Description of sub-models’) at basal and breast height are set to one. The reproductive status is set to not bearing any cones and finally, the height of a tree when it can start to bear cones is defined. 2.2.1.5. Input data. Climate time-series drive the processes of the model, namely monthly mean values of temperature and monthly precipitation sums. In the sub-model ‘computation of weather data’ the model calculates ‘net degree days’ and ‘active air temperature’

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

– two relevant climate auxiliary variables from climate series, to which larch species are reported to be sensitive (Abaimov, 2010). The first is the number of days exceeding 0 ◦ C and the latter is the sum of days with temperatures higher than 10 ◦ C for each year. 2.2.1.6. Observation. Information from tree stands in simulated plots is post-processed and written to data files. Sampling is conducted on all individuals growing in a central 20 × 20 m plot within simulated areas as well as of population measures such as tree density. 2.2.2. Description of sub-models 2.2.2.1. Computation of weather data. In the initialization phase of each run, assembled monthly weather data are processed to estimate daily values which are the basis for calculating the sum of temperatures of days exceeding 10 ◦ C, called ‘active air temperature’ (AAT10 ), and the vegetation period length by the sum of days exceeding 0 ◦ C, called ‘net degree days’ (NDD0 ) for each simulation year. Using these values and the mean temperatures of January and July, long-term trend variables are derived to account for inertia in the system. To estimate quasi-daily temperatures from yearly July temperatures (TJuly ) and the day of the year (Nday ), a sinusoidal function was fitted to observation data.







T Nday , TJuly =



TJuly + sin





2



365 ∗ Nday − 100

− 25

(1)

To estimate the maximum growth rates for individuals growing in each plot, AAT10 , NDD0 , and precipitation sum are used to calculate indices that follow a set of rules which are developed from references concerning species limits and responses to these climate proxies and kept as general as possible. Because of the ability of L. gmelinii to grow at low temperatures (c. 250–300 ◦ C AAT10 ), and to form close stands if exceeding 600 ◦ C AAT10 (Abaimov, 2010), the index (AATI) of each year follows a positive linear relationship, with AATI = 1/600 ∗ AAT10 . The second climate variable derived from temperature data, NDD0 , has to exceed values of 60 days above 0 ◦ C to allow L. gmelinii to grow (Abaimov, 2010). Thus, the NDD0 index (NDDI) of each year also represents a linear relationship, resulting in negative values for NDD0 smaller than 60 days, with NDDI = 1/60 ∗ NDD0 − 1. A precipitation index (PRI) is calculated, based on a negative correlation between tree-ring width and precipitation (cf. Sidorova et al., 2010) if the year’s precipitation exceeds an estimated optimum value. Therefore, a negative linear relationship with a slope (fPR ) is assumed, which reaches a value of one at the optimum precipitation and increases the smaller the amount of precipitation is below its optimum, such   that for each year PRI = 1 − fPR ∗ PR − PRoptimum /PRoptimum . 2.2.2.2. Update of environment. Space, which in the presented model is a proxy for suitable area and resources, is limited, and individuals growing in a certain area compete for it with surrounding individuals. This is represented by a density map and updated once at the beginning of each year. Depending on its size (basal diameter as DIA), a tree can influence a defined area in its environment with a strength decreasing linearly with increasing distance from its stem. We implemented this kind of zone-of-influence approach around individuals in a way that has been widely applied in a variety of individual-based models (see review of models in Czaran and Bartha, 1992; or Martínez et al., 2011). The density influences (DE) for each tile in the vicinity of a tree follows DE (xdiff, ydiff ) =



DIAbasal xdiff 2 + ydiff 2 + 1

(2)

107

where xdiff, ydiff is the position of a density grid tile relative to the tree stem, determined by a maximum value for distance which is calculated for each tree by xydiffmax (tree) = scalingfactor ∗ (DIAbasal ) /100 where scalingfactor is a factor that can be used to assign the area of influence of a tree (a rough estimate of 10 leads to an area of approximately 16 m2 for trees with basal diameters of 80 cm, after investigations of Kajimoto et al. (2007)). The sum of the density-influence values of a tree (TDEneeded ) is used as a proxy for the needs that allow for optimal growth. Thus, the individual value is compared to the sum of all the DE values of the tiles in its zone-of-influence (TDEreal ) and an index (TDEI) is calculated for each tree in a year, which increases with tree density, by TDEI = 1 − (TDEneeded ) / (TDEreal ). This density index can optionally be scaled by using tree size (basal diameters DIA) to represent an increasing intensity of the demand for resources of bigger trees. The user can manipulate the subsequent response function by adjusting the exponent (initial value: 0) to achieve the updated tree density ˆ index TDEI.



ˆ (tree, year) = TDEI (tree, year) ∗ TDEI

1−

0.01 DIA (tree, year)

exponent (3)

2.2.2.3. Growth. Growth of all individuals is implemented as a yearly increment in diameter. While individuals below 1.3 m in height only increase their basal diameter, larger ones additionally possess a diameter at breast height. The baseline increment values were deduced from performances of tree stands in Khatanga at sites 2, 3, and 5 (see details in Section ‘2.1 Reference sites’) and are increased or decreased depending on the current weather and the biotic environment for each individual. Due to high variance in diameter growth of individuals exceeding a height of 1.3 m (Fig. A2), individual diameter measurements are separated at this threshold. For the resulting two height classes, we used the upper, i.e. third, quartile of realized growth rates as an estimate for maximum growth rates, which for individuals <1.3 m is the maximum basal diameter growth rate of 0.04 cm year−1 (n = 18) and for individuals ≥1.3 m the maximum diameter growth rate at breast height of 0.06 cm year−1 (n = 53). The latter was also assigned for growth at basal height of the same height class. For each simulation year the maximum yearly growth of diameters (Gmax ) at basal and breast height for all individuals is influenced by the current weather, namely precipitation by the index PRI, and temperature by the ‘active air temperature’ index (AATI) and ‘net degree days’ index (NDDI). Based upon the standard value for growth of diameter (Gstandard,diameter ), Gmax is = PRI ∗ Gstandard,diameter ∗ calculated for eachyear by Gmax,diameter   1/fAATNDD ∗ AATI + 1 − 1/fAATNDD ∗ NDDI where fAATNDD is a factor balancing the influence of the two climate variables AAT10 and NDD0 . The yearly maximum growth values are used to derive the current growth rates for each individual affected by biotic influences such as intraspecific competition, triggered by its density index ˆ (TDEI), which, for example, accounts for root competition (Sofronov et al., 2000). Thus, current diameter growth at basal or breast height (Grealized ) is calculated for each tree and year according to Grealized,diameter = Gmax,diameter ∗ (1 − TDEI). Finally a tree’s current height (H) is estimated by using the following height-diameter relationship parameterized by using the same reference sites as above: for height estimation of the simulated individuals <1.3 m a linear relationship between height and basal diameter (DIAbasal ) is assumed and a non-linear one between height and diameter at breast height (DIAbreast ) for taller trees:



H=

44.43163 ∗ DIAbasal , H < 1.3m (7.02 ∗



2

DIAbreast ) + 130, H ≥ 1.3m

(4)

108

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

2.2.2.4. Seed dispersal. Seeds produced from trees are not dispersed into the environment before spring of the following year (Abaimov, 2010). Their chance to be released from cones is set initially to 70%. Released seeds disperse randomly in all cardinal directions; all other seeds endure another year within the cones. Seeds that are dispersed beyond the borders of simulated areas are removed from the system. For estimation of dispersal distances, we randomly draw the dispersal distance for each seed from a Gaussian-type function that allows for dispersal into the surrounding area but not exactly at the position of the releasing tree (Levin et al., 2003; Nathan and Muller-Landau, 2000). This was combined with a fat-tailed dispersal function to allow for infrequent long-distance dispersal and also increased the probability for seed dispersal close to the dispersing tree (see Appendix A5). For the dispersal distance (DL) of nearly 95% of the seeds to be within a certain radius around a tree, a distancescaling factor (sdist) is adjusted. The final calculation weights the distances calculated for each function by 20% Gaussian type, and 80% fat-tailed (multiplier of 0.5 and 2 respectively): DLgaussianwithfat−tail (rn) = sdist ∗ 0.5 ∗



 

0.5 ∗ DLgaussian (rn) + 2 ∗ DLfat−tailed (rn)

 (5)

where DLgaussian (rn) is a Gaussian-type function, leading to seed dispersal with peak density away from the central standing tree but with  decreased chance for long-distance dispersal: 2 ∗ width2 ∗ −1 ∗ log (rn), where width is a scalDLgaussian (rn) = ing parameter defining the span of the dispersal kernel (= 40) and DLfat -tailed increases the amount of very short-distance dispersal and also emphasizes long-distance dispersal and is therefore called ‘fat-tailed’, DLfat−tailed (rn) = rn(−1∗(1+∝)) where ␣ is a parameter for scaling the fat end of the curve (= 0.5). 2.2.2.5. Seed production. Trees mature and start to produce seeds after reaching a height that was deduced from the investigated tree stands. Rules that determine the maturation of trees, defined as beginning to bear cones, were derived from information about trees growing at sites 2, 3, and 5 in Khatanga (see ‘2.1 Reference sites’ for details). Height is used as a proxy for age of maturation. In the first of the two plots at site 2, a small fraction of trees began to mature and bear cones below 2 m, while at other sites trees first began to mature when taller than 2 m (see Appendix A1). Mean fractions of cone-bearing trees within size classes 0–2, 2–3, 3–8 m and greater than 8 m were respectively 0.3, 31, 73 and 94%. Thus, the relationship between maturity and height can be described by a logistic function and was parameterized to roughly follow the fraction of trees bearing cones (F) with 50% at 3 m and 90% at 11 m and going asymptotically towards one: FCones (H) = 1.0 −

1.0 H − 1.0

(6)

The calculation was performed in the following manner: for a sequence of heights beginning at 2 m and up to 1000 m in centimetre steps, newly scaled fractions (from 0 to <1 in 200 steps, thus reaching from 0 to 199) were calculated. The resulting series of fractions for heights has duplicate entries for some height values and consists of 183 entries (see Appendix Table A1). To account for the observation of a minor fraction of individuals <2 m tall producing seeds (see Appendix A1), we implemented a 0.5% chance to start seed production randomly between 1 and 2 m. Once a tree reaches seed-production height, it produces seeds every year (S*) which are located on the trees in the ‘cones’. The amount is dependent on tree height (H), weather conditions and competition, and is mediated by the current relative basal diameter growth (DIAR). Triggered by a scaling factor (fs ), the last term of the following function weights

seed production by tree height, beginning with 50% for 1 m individuals and rising asymptotically towards 100%, reaching 75% at 2 m and 95% at 10 m. The final number of seeds is:



∗

S =



fS ∗ DIARbasal ∗

1.0 −

 H −1.0 50m

(7)

2.2.2.6. Establishment. Each seed on the ground has a probability to germinate (PG) that consists of a ‘background germination rate’ (bgg) improved with a certain strength (fwg) by a current weather quality estimate. This is represented by the ratio of maximum realizable basal diameter growth at the seed’s position (Gmax ) divided by the maximum basal diameter growth of the current year (Gstandard ):



PG (year) = bgg +

fwg ∗

Gmax,basal (year) Gstandard,basal



(8)

2.2.2.7. Mortality. The mortality rate of a tree (mtotal ) depends on abiotic and individual factors which are summed. Simulated trees are removed from the plot when they are dead. The ‘background mortality’ (mbackground ) represents random effects that are not incorporated into other factors. Observed age distribution is used for parameterizing the mortality function. Patterns showed a continuous or clumped establishment and could begin or cease in a certain year (Wieczorek et al., submitted). Where tree discs were available, the rings were counted at basal height (basal age BA) multiple times in different directions and the maximal value chosen. Ages are estimated from basal diameters (DB) for remaining individuals. If the DB was also not available (as was the case for many individuals <0.4 m), age was estimated from the height (H). Hence, at each site the maximum height of the group of individuals without information is detected, followed by a linear regression analysis of individuals that are taller than this (red lines in Appendix A3). The sampled variables often departed from a normal distribution and were therefore transformed to conduct a linear regression. However, they were also used for analysis if the transformation failed to establish a normal distribution. From the fitted functions for estimating the basal diameter, the minimum basal diameters (DBmin ) are deduced for the smallest heights (Hmin ) found, which are in turn used for the parameterization of the ˆ for following linear functions. To estimate the basal diameter (DB) heights (H) from zero to Hmin (blue lines in Appendix A3) we used the function: ˆ (tree) = DB

DBmin ∗ H (tree) Hmin

(9)

Thereafter, estimation functions for basal ages at the three sites including individuals with only measured DB are parameterized (red lines in Appendix A3). The smallest available age values estimated by these functions are then used for parameterization of the linear functions, and finally used to estimate the age of smaller individuals for DB beginning with zero to DB min (blue lines in Appendix A3): ˆ (tree) = TAmin ∗ DB (tree) TA DBmin

(10)

For individuals >DBmin , the function parameters are directly used to estimate tree age at basal height (TA) with link functions (I) of response or dependent variables (RV, DV) and the inverse form indicated by the exponent of −1 and a slope of the linear function (m) as well as a y-axis intersection (n): −1 TA (tree) = IDV (m ∗ IRV (DB (tree) + n))

(11)

Finally, as a rough estimate, a simple linear regression of the distribution of ages is applied to a range of individuals beginning

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

109

Table 3 Overview of model parameters and processes. 18 parameters that lack a robust estimate were adjusted in a parameterization procedure: simulations were performed using randomly drawn parameter values from a normal distribution (start values may be unrealistic and are used only to calculate the standard deviation needed for determining the range of the normal distribution and random values are drawn until lying within the lower or upper extremes) and finally calculating the median of those parameter estimates that were able to capture the observed field data best (see 2.3 Parameterization for details). Parameter

Growth Growth rates - basal diameter (H <1.3 m & ≥1.3 m) - diameter at breast height (H ≥ 1.3 m) Height-diameter function slopes (H < 1.3 m & ≥1.3 m) Precipitation optimum Precipitation index strength Precipitation growth influence exponent Scale of influence between the indices AAT10 and NDD0 Influencing factor for zone of trees Density mode Influence of relative growth on other processes Seed production, dispersal and establishment Dispersal mode Factor of seed productivity

Minimum age to begin to bear cones Seed background input for population initialization Years with seed background input Probability of seed release from cones Factor of dispersal distance Background germination rate Influence factor of weather on germination rate Mortality Background mortality rate Tree youth influence factor on tree mortality Span of tree youth mortality Maximum tree age Influence factor for trees older than the age limit on tree mortality Current tree growth influence factor on tree mortality Influence exponent on current tree growth mortality Weather influence factor on tree mortality Density influence factor on tree mortality Drought influence factor on tree mortality Seed mortality rate on trees (in cones) and at the ground

Value and dimension

Parameter adjustment

References

Start values (and lower and upper extremes in brackets)

Final value (median ± SD)

0.04 & 0.06 cm yr−1

–

–

data-based estimate

0.06 cm yr−1

–

–

data-based estimate

44.43163 cm cm−1 & 7.02 cm0.5 cm−0.5 250 mm 0.1 1.0

–

–

data-based estimate

0.0, 500 (0–500) −0.8, 1.0 (0.0–1.0) –

371.718 ± 14.622 0.095 ± 0.039 –

estimated estimated estimated

2.0

−96, 100 (1–100)

62.990 ± 14.799

estimated

10.0

–

–

Kajimoto et al. (2007)

‘zone of influence’ yes

– –

– –

– –

– Gaussian with fat-tail

–

–

6

−8.0, 20 (>0.0–20)

16.044 ± 1.356

15 years

–

–

functions follow Nathan and Muller-Landau (2000) literature-based estimate (Kruklis & Milyutin, 1977, cited in Abaimov, 2010) Abaimov (2010)

1000

–

–

–

100

–

–

–

0.7

0.0, 1,4 (>0.0–1.0)

0.639 ± 0.086

estimated

0.16

−1.68, 2.00 (>0.0–2.0)

0.277 ± 0.074

– 0.5

−1.0, 1.0 (0.0–1.0) 0.0, 1.0 (0.0–1.0)

0.401 ± 0.128 0.448 ± 0.113

literature-based estimate (Abaimov, 2010) estimated estimated

0.0005 yr−1 0.85

−0.099, 0.100 (0.0–0.1) 0.0, 1.70 (0.0–<1.0)

0.024 ± 0.003 0.793 ± 0.264

data-based estimate data-based estimate

0.5 609 years

−1.52, 2.00 (0.0–2.0) –

0.763 ± 0.068 –

10

0.0, 20 (0.0–20)

8.188 ± 2.965

data-based estimate Vaganov et al., 1999, cited in Abaimov (2010) estimated

0.1

−0.8, 1.0 (0.0–1.0)

0.561 ± 0.109

estimated

1.0

–

–

–

0.5

0.0, 1.0 (0.0–1.0)

0.425 ± 0.057

estimated

0.1

−0.8, 1.0 (0.0–1.0)

0.248 ± 0.098

estimated

0.1

0.0, 1.0 (0.0–1.0)

0.238 ± 0.073

estimated

0.8 yr−1

0.0, 1.6 (0.0–<1.0)

0.447 ± 0.121 & 0.558 ± 0.120

estimated (Abaimov, 2010)

AAT10 : sum of temperatures above 10 ◦ C ‘active air temperature’, NDD0 : number of days exceeding 0 ◦ C ‘net degree days’.

with an age of 25 years, resulting in a slope of c. 7.5·10−4 %. Hence, the background mortality rate was set to an approximate value of 5.0·10−4 % as a start value which was then increased during the parameterization process.

A weather-dependent term (mweather ) introduces the influence of the current weather. A piecewise function is computed using three long-term weather variables (Value) and an index is computed with the following rules: below a threshold for species survival—mean temperature of January (−40 ◦ C) and July (10 ◦ C)

110

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

and 60 ‘net degree days’—the index becomes one, or, if exceeding this threshold, the index values decrease with differing response strengths for these three variables. These weather indices are calculated by: Valuevariable − Valuethreshold Valuethreshold

IWvariable = 1 − frs ∗

(12)

where frs is the factor for response strength (January isotherm = 8, July isotherm = 2 or long-term mean ‘net degree days’ = 1), Value is either the value of the weather variable in a given year or its threshold value. These values are then used to calculate the additional weather mortality by: mweather = fweather ∗ (1.0−



1.0 + e

−



160−



1.0

variable

60∗IWvariable





∗ H −1

wexp

(13)

∗Gmax,basal

where fweather is the factor for total influence of weather on mortality rate, and wexp is an exponent for scaling the height influence on weather mortality. Furthermore‘ ‘drought mortality’ (mdrought ) due to drought stress is estimated from evaporation rates throughout the vegetation period approximating estimations of Dolman et al. (2004) who investigated larch tree stands in far eastern Siberia growing on continuous permafrost. The calculation of the precipitation deficit follows:



August

Precdiff = 1.5

mm ∗ NDD/ day

PREC

(14)

May

If current evaporation exceeds precipitation throughout the vegetation period, an index for drought mortality (IDM) is calculated which uses the ‘active air temperature’ of the year if it exceeds a given threshold:



IDM (year) =

1 ∗ 10



0.0 , Precdiff > 1andAAT <= 600◦ C



AAT −1 300◦ C

, Precdiff > 1andAAT > 600◦ C

(15)

This is lowered by the tree’s height due to the positive relationship between tree height and drought or fire resistance. Additionally, a factor (fdrought) triggers the total influence of drought stress on mortality rate, such that ultimately, for each tree a mortality rate due to drought stress can be estimated with the formula: √ mdrought = fdrought ∗ IDM ∗ (H −1 ) (16) An age-dependent factor (myouth ) is introduced that increases mortality of small individuals. A limit until which individuals are affected is adjusted by the exponent yexp. The maximum value is determined by the parameter fyouth and its course along the curve depends upon height without the current year’s performance (Hscaled ). Its calculation is: myouth = fyouth ∗ exp(−Hscaled )yexp

(17)

For very old trees exceeding a given age threshold of 609 years (Abaimov, 2010), a tenfold background mortality rate (mage ) is added. A ‘density-dependent mortality’ (mdensity ) accounts for increasing tree density that triggers negative impacts, such as reduced resource availability or higher probability of infection by a disease or pathogens in denser tree stands. The maximum influence of this factor is determined with a parameter (fdensity) mdensity = fdensity ∗ TDI

(18)

Finally, mortality is increased depending on the individual tree’s performance (mperformance ). We used the relative diameter growth (DIAR) at basal (for individuals <1.3 m) or breast height (for individuals ≥1.3 m). The strength of influence is adjusted by the exponent (pexp). The calculation of this factor follows: mperformance = 1 − (DIARdiameter )pexp

(19)

Seeds, which are either situated within cones on trees or lying on the ground, have a fixed rate of mortality (=80%). This estimate is based on germination energy values of seeds in accordance with investigations by Abaimov et al. (2000) who found that 9–20% of seeds produced in one year were still viable the following autumn. If a seed dies, it is removed from the system. 2.2.2.8. Ageing. At the end of each simulated year, the age of all individuals is increased by one. Seeds of L. gmelinii have been shown to be dormant for up to three years (Abaimov, 2010), so seeds in this simulation can germinate until reaching this age, thereafter they are removed from the system (Fig. A3). 2.3. Parameterization In the calibration of sub-models, we based parameter estimates on field investigations and from the literature (Table 3). Where a robust estimation was lacking, we developed a procedure to parameterize the model with field data to further increase the quality of these parameter estimates (Fig. A4). To obtain robust estimates for 18 uncertain parameters (with a given ‘range for parameterization’; see Table 3) we started eight simulations simultaneously at five locations (Table 1) with a maximum of 5000 repeats, in which the parameters were estimated from the performance of their predecessors. The model followed a three-step procedure: (1) randomly assigning the parameters from a normal distribution using the weighted average and standard deviation of simulation run weights; the ranges were the first two values with equal weights of 1, (2) performing one simulation with 1000 years stabilization and a subsequent 93-year simulation time (years 1919–2011), and (3) after the final year of the simulation period, resulting ‘basal areas’ and numbers of individuals in three height classes, ‘seedlings’, ‘saplings’ and ‘trees’, were compared to corresponding field data (Table 1) and a weight calculated as the mean of the absolute values of all four relative departing values. A weight of 100 was set as a maximum value and if no growth was possible at any simulated plot, the weight was divided by 1000. The list of parameters and their quality estimate/weight is therefore amended after each simulation repeat. If ten simulations were run, the starting range values were deleted. After a further ten repeats the current parameter set is added to the list by replacing the set with the lowest weight, to allow the model to converge faster at a ‘good parameter estimate’. We based our evaluation of the resulting parameter sets on the deviation of the number of trees compared to field observation. Only one of these simulations produced stable populations, i.e. at all locations trees were growing in the final year 2011. Here, 3.1% parameter sets had a departing value of less than 0.1. From the 16 simulations with a departing value <0.1 from the first 500 simulations the median values were calculated (Table 3) and further tested for reproducible results. In 100% of the simulations conducted with the resulting parameter set, trees could establish at the five locations. Of 100 repeats, 11% were with an absolute mean departing value of <0.1 (not shown). 2.4. Khatanga climate time-series Monthly mean temperature and precipitation sums are mandatory to run the model. For temperature and precipitation, we used

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

111

Table A1 Maturation heights, which are used to determine the beginning of seed production. Random number (RN, 0–182)

Maturation height (MH) [cm]

RN

MH

RN

MH

RN

MH

RN

MH

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

100–200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236

37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73

237 238 239 240 241 242 243 244 245 246 248 249 250 251 252 253 254 256 257 258 259 260 262 263 264 266 267 269 270 271 273 274 276 277 279 281 282

74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110

284 286 287 289 291 293 295 297 299 300 303 305 307 309 311 313 316 318 320 323 325 328 330 333 336 339 341 344 347 350 354 357 360 364 367 371 374

111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147

378 382 386 390 395 399 404 408 413 418 423 428 434 439 445 451 458 464 471 478 485 493 500 509 517 526 535 545 555 566 577 588 600 613 627 641 656

148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182

672 689 707 725 746 767 790 815 841 870 900 934 970 1010 1053 1100 1153 1212 1277 1350 1434 1529 1639 1767 1919 2100 2323 2600 2958 3434 4100 5100 6767 10100 20100

the Climate Research Unit data set CRU TS 2.1 available at a latitudinal and longitudinal resolution of 0.5 ◦ for the years 1901–2002 (Mitchell and Jones, 2005). We extracted the data for locations with field observations (Table 1). This was directly achieved for coordinates of sites 2, 3, and 5 at which open-forest occurred. These data series were prolonged by modern data for the period 2003–2011 from the weather station situated in Khatanga (data records beginning in 1920 with some gaps in 1934, 1944, 1945; 71.983◦ N and 102.467◦ E, WMO-ID: 20891), which is the nearest weather station to the investigation area with the longest available record (data accessed on 2nd May 2012 from Global Historical Climatology Network (GHCN)-Daily; Menne et al., 2012). For each local series, the mean of the CRU data for years 1961–2002 was compared to the station data and the difference added to achieve a local adaptation of the most recent period. Ultimately, time series running from 1901 to 2011 were available. Because of low variation in the gridded CRU data in this region before 1919 we used only years 1919–2011 of this series for the evaluation of the sub-models, sensitivity analysis, parameterization, and temperature experiment.

stabilization phase of 1000 years. This was followed by the empirical phase for 93 years (years 1919–2011) using weather time series for site 3. To evaluate the sensitivity of a ±5% changed parameter in the model we compared the results with one baseline simulation (30 runs) without any parameter changes. Using the results of the final year from these simulation runs, nine key variables were calculated for model evaluation. The general stand characteristics (1–5) are indicated by the number of seedlings (NSeedlings ), saplings (NSaplings ), and trees (NTrees ), as well −1 as the stem count (number  of trees  >1.3 m ha ) (stemcount). The mean age of trees >2 m AgeTree

was also calculated. The individ-

ual performance of trees by the cumulative seed (6) is calculated  production of all trees N Seeds . Additionally, variables with the purpose of tracing spatial patterns of individuals (7–9) were gained by applying Ripley’s K(t)-function (Ripley, 1977) to each simulation result separately for seedlings, saplings, and trees. Results were further processed to L-values for each radius r by L(r) = (K(r)/)0.5 . With these, we estimated indices for the three height classes (ripSeed, ripSapl, ripTree) via the following protocol:

2.5. Sensitivity analysis A local sensitivity analysis (Grimm and Railsback, 2005) was performed to highlight parameters with a strong influence on the model’s outcome as a way to check their parameterization (Cariboni et al., 2007). Thirty parameters were changed separately and simulations performed with ±5% of the reference values of each parameter (see Appendix Table A2). For each parameter, we carried out 30 simulations and started population growth by introducing 1000 seeds per year in the first 100 years of each simulation and a

• generate 39 random point patterns for each of the groups to estimate significance envelopes • calculate ratios between positive values and the upper confidence envelope and similarly the ratio between negative values and lower confidence envelope at each radius, up to a maximum of 5m • bin the observed values to one of three possible cases: where positive values fall below negative ratios they are counted as significantly uniformly distributed (=1), those that are equal as

112

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

Table A2 Characterization of 30 model parameters which were tested in the sensitivity analysis. For each parameter, the reference and varied values as well as their dimension are given. Parameter

Abbreviation

Reference value

Values −5%

+5%

Type

Units

Slope of height-basal diameter estimation function (H < 1.3 m) Factor of height-breast diameter estimation function (H ≥ 1.3 m) Factor of precipitation index strength Potential evapotranspiration Factor of dispersal distance Seed background input at start phase Permanent seed introduction Years with seed background input Basal diameter growth rate for: − H < 1.3 m − H ≥ 1.3 m Growth rate at breast height Maximum tree age Influencing zone of trees Background germination rate Influence factor of weather on germination rate Influence factor for trees exceeding age limit on tree mortality Density influence factor on tree mortality Tree youth influence factor on tree mortality Background mortality rate Span of tree youth mortality Drought influence on tree mortality Weather influence factor on tree mortality Precipitation optimum Seed mortality rate - on the trees - on the ground Probability of seed release from cones Factor of seed productivity Minimum age to begin to bear cones Factor balancing the influence of two climate variables AAT10 and NDD0 Current growth influence on tree mortality

sLinear,BasalDiameter sNon-Linear,BreastDiameter fPRI fEvapotranspiration Sdist NSeedInput NSeedinputpermanently Seed Input Years

44.43163 7.02 0.095415 1.5 0.277025 1000 10 100

42.2100485 6.669 0.09064425 1.425 0.26317375 950 9 95

46.6532115 7.371 0.10018575 1.575 0.29087625 1050 11 105

double double double double double int double int

cm cm−1 cm cm−1 – mm yr−1 – – – –

GMaximum,BreastHeight GMaximum,BreastHeight GMaximum,BreastHeight Maximum Age fHAI fBackgroundGermination fWeatherGermination fAgeMortality fDensityMortality fYouthMortality mBackground Yexp fDroughtMortality fWeatherMortality PROptimum

0.04 0.06 0.06 609 10 0.40063 0.447975 8.18785 0.248435 0.762855 0.02382 0.79295 0.237805 0.42539 371.717745

0.038 0.057 0.057 578 9.5 0.3805985 0.42557625 7.7784575 0.23601325 0.72471225 0.022629 0.7533025 0.22591475 0.4041205 353.1318578

0.042 0.063 0.063 640 10.5 0.4206615 0.47037375 8.5972425 0.26085675 0.80099775 0.025011 0.8325975 0.24969525 0.4466595 390.3036323

double double double int double double double double double double double double double double double

cm yr−1 cm yr−1 cm yr−1 years – – – – – – yr−1 – – – mm

PSeedMortality,Cones PSeedMortality,Ground PSeedRelease fS Maturing Age fAATNDD fGrowthMortality

0.44724 0.55803 0.63931 16.043965 15 62.89988 0.561275

0.424878 0.5301285 0.6073445 15.24176675 14 59.754886 0.53321125

0.469602 0.5859315 0.6712755 16.84616325 16 66.044874 0.58933875

double double double double int double double

yr−1 yr−1 – – years – –

AAT10 : sum of temperatures above 10 ◦ C ‘active air temperature’; NDD0 : number of days exceeding 0 ◦ C ‘net degree days’.

randomly or not specifically distributed (=2), and where negative ratios exceed the positive ones as significantly clumped (=3).

2014). For each of the 1815 locations we performed 10 simulation repeats with the same settings as used in the sensitivity analysis.

Finally, all resulting variables were computed to sensitivity indices that present the percentage of absolute change scaled to the change in parameter values of ±5% (S+ and S− ) calculated by:

2.6.2. Temperature experiment For this experiment, temperature series were constructed by specific increases and decreases of yearly temperatures. We used weather data over 61 years (1951–2011) from the Khatanga timeseries, which lies in the vicinity of investigated site 3 (71.8◦ N, 102.9◦ E). From this, we drew random years and generated seven 210-year-long temperature series in a simple setup by changing the monthly temperatures by −6, −2, −1, 0, +1, +2 and +6 ◦ C. Furthermore, to compare the effect of a direct change to incrementally changing temperatures over 210 simulation years, six additional time series with linearly increasing temperatures were generated which reach, after 100 years, the same final temperature changes as used for former directly-changed temperature series. We used the same randomly generated time series as the basis for all series throughout the experimental phases to ensure clear effect-separation from pure climate-fluctuation influence. Simulation time for the experiments was 210 years following a 5001-year-long stabilization phase with random weather from the time series generated for this experiment with a yearly mean value of c. −12.9 ± 1.4 ◦ C. Stabilized populations were then exposed for 30 times to normal, increased, or decreased temperature and the resulting population dynamics recorded.

S+/− =

V+/− −VREF VREF P+/− −PREF | P | REF

(20)

where V is the variable of interest derived from each simulation run and P is the parameter of interest, both plus (+) and minus (−) 5% of the estimated parameter, or with the reference value. 2.6. Model experiments With the parameterized model we ran simulations at single locations with climate series spanning the range of the Taymyr Peninsula and with modified temperature series for our temperature experiments. The results are available at https://doi.pangaea. de/10.1594/PANGAEA.863584. 2.6.1. Taymyr treeline application First, we used local climate data to force simulations at field sites and compare the results with observations. The settings were similar to those used for the sensitivity analysis: a 1000-year-long stabilization period with random weather followed by the simulation period from 1919 until 2011 (Fig. 3). We calculated summary statistics for 30 simulation repeats. In a second step, to test the model for a wide range of climates, we extracted weather data for the Taymyr Peninsula ranging from 63.75–79.75◦ N and 91.75–118.75◦ E at a 0.5◦ latitude and longitude resolution from the gridded data set CRU TS 3.22 (Harris et al.,

3. Results 3.1. Sensitivity analysis 3.1.1. Sensitivity of simulations at central treeline location Changing the parameters by ±5% only moderately influences the performance of simulated tree stands at site 3 in the equi-

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

librium phase: spatial patterns (P1-3; mean absolute values 0.3 ± 0.4% (S+) and 0.2 ± 0.4% (S−) variable−1 ) are generally less sensitive than stand characteristics (S1–S6; mean absolute values 0.9 ± 1.2% (S+) and 0.9 ± 1.2% (S−) variable−1 ) or individual performances (I1–I3; mean absolute values 1.6 ± 1.5% (S+) and 1.8 ± 1.5% (S−) variable−1 ). Thus, individual counts, i.e. number of seedlings, saplings, and trees, are as equally sensitive as individual performances to parameter changes but sensitivity values are especially low for the spatial pattern regarding variables with values <0.1 ± 0.1%. Overall, mean probabilities of variables to deviate from a reference value are low, i.e. for all parameters it is <65%, and in only eight out of 30 is it >50% (Fig. 4). In the development phase at year 500, when grouped into the three categories, mean sensitivity conserves the same order of importance while spatial patterns stay the same (S+ 0.28 ± 0.32 and S− 0.2 ± 0.3), and stand characteristics (1.5 ± 1.7 and 1.6 ± 1.7) and individual performances (3.1 ± 1.8 and 4.3 ± 1.8) increase. In 27 cases (5% of all tested combinations) the defined keyvariables strongly deviate (with more than 95% probability) from reference values. Parameters with the highest sensitivity values throughout this analysis are ‘youth mortality rate’ (fYouthMortality ; mean of absolute sensitivity value 2.5 ± 2.1%), ‘slope of heightbasal diameter estimation function of trees with height (H) <1.3 m’ (sLinear,BasalDiameter,dbrustho ; mean of absolute sensitivity value 2.3 ± 1.9%), ‘factor of height-breast diameter estimation function (H ≥ 1.3 m)’ (sNon-Linear,BreastDiameter ; mean of absolute sensitivity value 2.3 ± 1.9%) and ‘growth rate at base of trees with height <1.3 m’ (GMaximum,BasalHeight,H < 1.3 m ; mean of absolute sensitivity value 1.9 ± 1.6%) (in Fig. 4 all sensitivity values). The direction of the resulting change is contrary: while increases in ‘youth mortality rate’ lead to decreased population numbers, ‘growth rate at base of trees with height <1.3 m’ and the two height estimation parameters show a positive relationship. Based on the deviation probability, the next parameter is ‘background mortality’ (mBackground ) with a moderate mean sensitivity value and a negative relationship with population density. For this parameter, the impact is stronger if increased than if lowered. The order of single parameter variations changes slightly; the first two were ‘youth mortality rate’ and ‘background mortality’, followed by ‘slope of height-basal diameter estimation function of trees with height (H) <1.3 m’. A parameter that is used to start the population growth in the beginning by introducing seeds to the area becomes important, i.e. ‘seed input years’ (for the first 100 years of each simulation 1000 seeds are placed randomly into the area; sensitivity increased from 0.3 ± 0.4 to 1.9 ± 1.9). 3.1.2. Sensitivity at different population stages and treeline locations Sensitivity values are lower if populations reached their equilibrium at simulation year 1083 (which is the final year of simulation equal to AD 2011) compared to the same populations at their exponential stage at year 500 (Fig. 5). Variations of absolute sensitivity values are low at around 0.3% of the corresponding cumulative sensitivity sum. Whereas values are higher by 200 for most of the sites, at the northernmost/coldest site the sensitivity is higher at final equilibrium stages. At the same time, the number of trees ha−1 is higher (283 ± 103 compared to 273 ± 127 trees ha−1 ) which is not the case at the other locations, where c. 300 more trees emerged until the final year. Sensitivity values decrease with increasing population density, reaching a minimum value of c. 350 at populations of c. 1500 trees ha−1 . 3.2. Taymyr treeline application In the final year of simulations at AD 2011, the numbers of trees growing in the central 20 × 20 m plot are in the range of field data

113

Fig. 5. Cumulative sensitivity values of simulation results with ±5% varied parameter values at each location at two different population stages: exponential growth at simulation year 500 and equilibrium plateau at 1083. The variation of sensitivity values for all is c. 0.3%.

Fig. 6. Stand densities of trees from simulated populations at five local climates compared at the final simulation year (= AD 2011) to field survey data at the treeline on the Taymyr Peninsula ordered by decreasing latitude. Each simulation included 1000-year-long stabilization and an empirical phase for years 1919–2011.

for the four north to mid-treeline sites (1, 2, 3, and 4; Fig. 6). The simulation for the southern extreme site 5 overestimates stand densities by c. 800 trees ha−1 (one-sample t-test of the 30 simulated tree numbers to depart from the mean of the two reference sites for each location: for site 5 p-value « 0.001, site 3 p = 0.0016 and site 1 p = 0.038; df = 29). On the other hand, simulated tree densities at site 1 slightly underestimate stand densities by 55 trees ha−1 . Simulations with temperature and precipitations series from the Taymyr Peninsula result in populations of trees up to 2530 trees ha−1 (Fig. 7). In general, populations could emerge under a variety of temperature–precipitation combinations. They are denser under prevailing warmer and wetter climate which corresponds to the lowland area around 71◦ N and 100◦ E (Fig. 1), reaching more than 1500 trees ha−1 at the end of the simulation period. Towards colder temperatures and less precipitation, stand densities decrease to populations of <500 trees ha−1 , until dying out at a certain threshold, which is the case for most populations growing in weather that is characterized by July temperatures below the 12 ◦ C isotherm and which also slightly coincides with precipitation levels of <250 mm. Throughout this simulation experiment the highest stand densities emerge at combinations of high July temperatures (exceeding at least 10 ◦ C) and precipitation sums (minimum c. 250 mm year−1 ). In warmer climate with temperatures exceeding a yearly mean of roughly −13 ◦ C, simulated populations are denser than their cor-

114

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

Fig. 7. Simulated number of trees at year 2011 for all possible temperature and precipitation time series at locations on the Taymyr Peninsula derived from CRU TS 3.22 (Harris et al., 2014). Stand densities are interpolated over yearly mean temperatures and precipitation sums of years 1960–1990 and corresponding densities plotted in grey shades with contour lines: field data locations are plotted as dark grey triangles. The region spanning the simulated combinations of climate time series is surrounded by a dashed line. The bold dashed lines indicate July isotherms of the same period.

responding field observations (here climate data was derived from its geo-position and the corresponding grid cell in the CRU TS cell).

3.3. Temperature experiments Generally, forcing the model with higher temperatures leads to increased stand densities, such that increasing temperatures by +6 ◦ C allows the populations to nearly double their size from c. 1200 to 2000 trees ha−1 . However, cooling affects populations even more strongly so that a decrease of -1 ◦ C triggered a reduction to less than half of the reference population, i.e. 500 trees ha−1 . Additionally, abrupt and transient change scenarios of colder than −1 ◦ C could lead populations close to extinction in all repeated simulations (Fig. 8). If temperatures change abruptly by 1, 2, or 6 ◦ C, population densities (trees ha−1 ) depart after 3–65 years if colder, and 16–118 years if warmer, than the reference run from the 95% confidence interval of baseline runs (Fig. 8a). Furthermore, observed trends of populations growing in cooler temperatures show non-linear behaviour, i.e. a 2 ◦ C decrease leads either to a rapid dying out of 50% of the populations after 99–130 years or to a strong diminishing in stand densities to only 250 trees ha−1 with a continued negative trend after 100 years of the simulation time. On the other hand, populations forced with 6 ◦ C warmer temperatures increase linearly at first, reaching local maxima after around 50 years, fall back to local minima around year 100 and finally reach a quasistable state with c. 2000 trees ha−1 at nearly the end of 150 years simulation time. A more modest increase of 1 ◦ C leads to a smaller and delayed local maximum and to final populations with densities of 1700 trees ha−1 , which is about 140% of the reference value. For comparison, following a 1 ◦ C cooling only 40% of the reference stand densities remained.

In a second approach, temperatures were linearly incremented per year to achieve 1, 2, or 6 ◦ C change over 100 years. The general trends are similar to the corresponding abrupt change scenarios, however the offsets of observable changes in tree population dynamics occur >30 years later, i.e. under 6 ◦ C changes began after 36 years for cooling and 65 years for warming (Fig. 8b). Furthermore, final populations under increased temperatures could, after a time-lag of about 20 years, linearly increase after 100 years to 1900 trees ha−1 . In cooler scenarios, populations could survive until the end of simulation time even with -2 ◦ C forcing, but still with a downward trend in stand densities. Such populations have c. 200 trees ha−1 at the end of the simulation. 4. Discussion 4.1. Assessment of LAVESI sensitivity We developed the spatially-explicit model LAVESI, dealing with seeds and trees as fundamental units, for the purpose of exploring population dynamics of larch forests growing on permafrost soils at the latitudinal treeline in Siberia. Although, the model does not explicitly consider permafrost thaw depth, simulations suggest that this new framework is reliable for this task. Based on the results of local sensitivity analyses we conclude that simulations with LAVESI are robust to uncertainty in the input parameters and we thus consider that the model has a high credibility to capture the major processes involved in population dynamics. We find that ±5%-parameter changes in almost all parameter-variable combinations result in small changes to the model’s output. The impact is higher if populations are in an exponential/development phase compared to populations that reached equilibrium. As expected, denser populations (that emerged in warmer climates) are less sensitive to a change in parameters than

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

115

Fig. 8. Population dynamics of simulated tree stands forced for 210 years by abruptly changed temperatures (left, a), or incrementally changed temperatures (right, b), after a 5001-year-long stabilization phase with randomized weather. Temperatures are given to the right of the plots for each experiment. Grey shaded areas represent the 95% confidence interval around the mean value of 30 reference simulations shown with a black line.

sparse populations at the edge of the species range (for example open stands observed at Ary-Mas (Kharuk et al., 2006)). Under such conditions, the destiny of a population either to succeed or fail to establish could be triggered by small parameter changes of 5%. For example, increasing the youth mortality or background mortality rate by 5% leads to a stronger negative impact than when decreased by the same amount. Results can thus change by up to 10% and become twice as strong. Among the tested parameters, a few combinations result in significant changes to the model’s output at both population stages: exponential and at equilibrium. First, the changing of ‘youth mortality factor’ results in a significant increase in population growth. This shows that populations can expand if a higher fraction of individuals is able to survive high competition and mortality rates at an early life stage. Our finding confirms results from field experiments with Larix and Pinus (Barbeito et al., 2012) where high mortality rates characterize early life-history stages. Second, changing the ‘slope of height-basal diameter estimation function (H < 1.3 m)’, ‘factor of height-breast diameter estimation function (H ≥ 1.3 m)’, and ‘growth rate at base of trees (H < 1.3 m)’ result in a significant increase in population growth. To achieve a good estimate for tree growth, we used regionally available tree-height and diameter data as well as tree-ring data and estimated factors for the height-diameter functions and the baseline possible growth as the upper quartile of observed tree growth increments. However, this estimate of the overall baseline value is potentially biased because local site conditions and interactions affect individual growth. To improve this parameterization, in a later version of the model, maximum growth can be based on tree-ring data of ‘well developed’ trees as applied by Fyllas et al. (2010). Nevertheless, the actual diameter growth of a tree in the model is influenced by its abiotic and biotic environment, i.e. the current weather and competition, so simulations should reliably represent the growth response even under a variety of scenarios. The spatial structure of tree stands can generally be influenced by seed dispersal and interactions among emerging individuals (e.g. Rammig and Fahse, 2009). From the observation of L. gmelinii’s strong investment in its rooting system (Sofronov et al., 2000), it can be assumed that the below-ground competition for soil activelayer space (the seasonally thawing soil space) is very important for the spatial pattern of Larix stands. We find that slight parameter changes do not significantly affect the developed spatial structures

of Larix populations. However, the sensitivity analysis was conducted using location specific climate and thus only the sensitivity of modelled tree-stands under these specific climates is discussed. For example, under warmer climates, denser populations emerge so that the potential space for further tree establishment becomes reduced and the importance of the different model parameters that influence spatial interactions among individuals might change, such that final spatial variables might increase in sensitivity. The latter could only be observed under developing populations, and sensitivity values for spatial variables were higher for denser populations. At the final stage of simulations, the lowest sensitivity values are observed for denser populations. Plot size might also influence the simulations (Shugart and West, 1979; Smith and Urban, 1988; Bugmann, 2001). To minimize boundary effects at our simulated forest plots we tested varying plot sizes and selected an area of 100 × 100 m. With this, our model shows similar responses to changes in area as forest-gap models, as noted in a review by Bugmann (2001). He reported that if the area of gap models varies by ±50% between 100 and 1000 m, typically only small effects could be observed. The size in our approach is similar to those used in a forest model application by Pacala et al. (1996, 1993) and larger than those used in current global models with simulated individual-based representation of vegetation, such as SEIB-DGVM with 30 × 30 m (Sato et al., 2007) or LPJ-GUESS with 10 × 10 m2 plots (Smith et al., 2001). Additionally, only the central 20 × 20 m plot is investigated which is of a similar size to the field sampling plots and thus ensures direct comparability. 4.2. Larix stand simulation under the Taymyr Peninsula weather Forcing the model with temperature and precipitation series for locations where field observations are available yields simulated tree densities that match, in general, with those found in the field (Fig. 6). For example, northernmost site 1 at ∼72.4◦ N fits well, with 257 ± 140 trees ha−1 , to the reference value of 200–425 trees ha−1 , while farther south at ∼71.5◦ N simulations yield populations similar to those found at open-forest sites 2, 3, and 4. However, simulations for the southernmost site at ∼70.7◦ N yield populations twice as dense as those found at site 5. This overestimation at one location may originate from ecological processes that are not covered by the model which could cause population declines or setbacks to earlier successional phases such as dis-

116

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

ease and fire (Kharuk et al., 2011; Schulze et al., 2012). Fires are a frequent and important factor at the southern boundary of L. gmelinii-forests and beyond, where warm temperatures lead to fast biomass accumulation and pronounced drought phases so that a high fire frequency below every 100 years can occur (Sofronov et al., 2000; Kharuk et al., 2011). Fires are not explicitly handled in the current version and as a further development and to be applicable for future warmer scenarios one should extend the model. Nevertheless, the high number of old trees and few juveniles observed may hint that an overmatured stand was surveyed; one where regrowth is supressed by a factor that has not yet been captured. Light might become an important limiting factor for regrowth in dense stands and a shift from below-ground competition for activelayer space to above-ground competition for light may occur along the investigated north–south transect. In which case, the model’s processes should be reconsidered and revalidated should an application aim to analyse typical taiga tree stands (see for example the modelling approach of Pacala et al., 1996). However, the nonfit between modelled and field-based density inferences may also result from atypical site conditions with respect to substrate or local climate. The thaw depth of permafrost soils (ALT) during the vegetation period might limit or enhance tree development further (Kajimoto et al., 2003). In this version of the LAVESI, we did not consider this factor explicitly because we intended to build a minimum adequate model capturing the population dynamics. Due to the fact that ALT is strongly linked to air temperatures (Zhang et al., 2005; Stefan Formula, Hinkel and Nicholas, 1995) and as, during the parametrization process, the model was optimized for field sites with potentially different thawing regimes, we assume that the influence permafrost might have is already considered in the model. Nevertheless, explicit handling of permafrost dynamics could help to elucidate their influence on stand dynamics and would increase the credibility of simulations and should be tested in a future version of the model. Finally, the assumed climate, as inferred from interpolation of the climate series New et al. (1999), may overestimate local temperatures because it is interpolated for the Taymyr Peninsula from few climate stations and heavily based on one in Khatanga next to a river. This could mean that important details about temperature and precipitation differences some hundred kilometres farther inland are missing. However, the comparison indicates that, overall, the model reproduced the field-data well, particularly for the northernmost and intermediate open-forest stands for which the model was designed. The regional simulations represent current climates of the whole Taymyr Peninsula and emphasize that forest development is restricted by cold temperatures or dry conditions, or both. Growth in this treeline region is assumed to be influenced mostly by temperature (Abaimov, 2010). Similar to the observations of MacDonald et al. (2008), simulated populations emerged in forms of open stands to dense northern-taiga stands between the 10 and 12.5 ◦ C July isotherms. The simulation results indicate a slight interaction between the two climate variables and forest development. We find that even when it is warm enough for trees to grow in dense stands (annual average temperature exceeding −12 ◦ C), a trend towards sparser populations in drier climates is observed. This influence of precipitation is corroborated by results of carbon and oxygen stable-isotope analyses of Larix tree-rings from northern Siberia (Sidorova et al., 2009). Additionally it has been found that tree growth of larches in north-central and far-east Siberia can be suppressed by dry summers in the preceding year (Lloyd et al., 2010; Berner et al., 2013), but locations are also known that show no correlation between tree growth and precipitation (Kirdyanov et al., 2008). Nevertheless, it is quite likely that the simulations present a treeline ecotone close to equilibrium with climate while the real current vegetation is in a non-equilibrium state with recent strong warming of nearly +3 ◦ C in the last 50 years

(anomaly of 2008–2013 compared to the 1960–90 mean measured at the weather station ‘Khatanga’, 71.983◦ N and 102.467◦ E, WMOID: 20891). This assumption is corroborated by our temperature experiments (discussed in Section 4.3). Simulation results portray that populations need 50–100 years to adapt to a warmer climate even under linearly increasing temperature scenarios of up to +6 ◦ C in annual mean temperature. 4.3. Transient Larix response to hypothetical future temperature changes Due to strong feedback mechanisms at the Arctic coast, particularly in relation to sea-ice decline, a temperature increase by up to 6 ◦ C is not an unrealistic scenario for the late 21st century as revealed by GCM runs forced by the IPCC representative concentration pathway scenario (“RCP8.5” IPCC, 2014). To test the response rate of larches in North Siberia against temperature change we forced simulated populations at equilibrium with current climate with arbitrary temperature series. We focus on the increasing temperatures, but one should keep in mind that the elevated atmospheric carbon dioxide concentration will also most probably directly enhance tree growth (Wullschleger et al., 1995). This effect, however, is not incorporated into our model. Our temperature experiments reveal that tree recruitment positively follows increasing temperature, yielding increased tree density. This prediction matches the observation of Kharuk et al. (2006) who found a recruitment increase of larch tree stands in north-central Siberia during recent decades, which they relate to the observed regional temperature increase. Stand densification and range expansion into formerly treeless tundra have been observed at most northern latitudinal and elevational treelines (e.g. Esper and Schweingruber, 2004; Harsch et al., 2009; Kharuk et al., 2006; Mamet and Kershaw, 2012; Mathisen et al., 2014; Shiyatov and Mazepa, 2012; Shiyatov et al., 2005). Our simulations are corroborated by observations, i.e. populations linearly densify under warmer scenarios and even double their size after 100 years in simulations with temperature increases of 6 ◦ C. Increased precipitation during winter or decreased precipitation in the vegetation period, which may induce drought stress and affect tree growth (Sidorova et al., 2009), were not tested here but might reduce the observed response rate to temperature changes. Nevertheless, temperature is the main driver of tree growth in this area so our simulation results suggest that larches might migrate north, albeit potentially more slowly than simulated due to dispersal limitations (Kharuk et al., 2006). Our experiments show that after initial establishment of trees following a pronounced time-lag, populations can quickly densify under optimal climate conditions. This supports the results of a comparative run of five Dynamic Global Vegetation Models (Sitch et al., 2008), which yielded an increase in woody cover on the Taymyr Peninsula by around 50% compared to 1860 by the end of the 21st century. Such future tundra-to-taiga transition is likewise supported by an Arctic vegetation zone modelling study by Pearson et al. (2013), and a general increase in biomass due to climate warming by forest-gap modelling throughout Russian Eurasia by Shuman and Shugart (2009). But, short-term cold spells could potentially limit fast population adjustments because our simulations show far stronger decreasing stand densities with temperature reduction than limitations from temperature increases by similar absolute values. For example, reducing the temperature by 1 ◦ C leads to a population decline to around 40%, while a similar temperature rise leads to an increase to 140% of the reference values. This aligns with the suggestion of Körner (1998) that treeline populations suffer more from single cold winters than they can profit from single hot summers. However, the slowdown of forest advance by cold winters is not a likely scenario because an increase in Arctic winter temperatures has been observed over the last decades and is

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

predicted for the future (Harsch et al., 2009; Meyer et al., 2015). In our experiments, cold spells are an unlikely cause of the time-lag in the population growth response to increased temperatures. The adjustment to an abrupt 1 or 2◦ increase in temperature is observed 70 or 120 simulation years following the temperature change as an increase in stem density. This slow response is because trees have to reach 2 m height to be counted as trees. Also, a possible response that biomass increases more rapidly in warmer scenarios is not incorporated into this version of our model, and an estimation of the above- and below-ground biomass of an area (Kajimoto et al., 2006) from simulated tree dimensions would be needed for a possible future extension. This would help us to understand the interactions of changing temperatures, biomass accumulation and fire frequency. In this experiment we used modified time series with annually increased or decreased temperatures although the recent change in temperatures was seasonally different: whereas temperatures in summer and winter increased by c. +2 ◦ C, spring and autumn temperatures increased by c. +3 to +4 ◦ C (mean of 2008–2013 compared to 1960–1990 measured at weather station ‘Khatanga’). Our simulations suggest that under warmer temperatures, larch stands in northern Siberia will most probably densify and are likely to migrate north, replacing treeless tundra. Our experiments show this response may lag behind by decades, similar to the indications of Epstein et al. (2007). This slow response is revealed by using a new model which incorporates all life-cycle stages of larch and the interactions among individual trees of all sizes, from seedlings to mature trees. Finally, even if temperatures do not reach the high levels predicted for the most aggressive IPCC CO2 concentration pathways, densification and northward expansion of larch stands seem inevitable during the latter half of the 21st century.

117

throughout the treeline of Siberia. Furthermore, if information on stand connectivity due to seed or pollen dispersal among locations can be reliably assessed for the treeline region, for example from population genetic studies, the model could be parameterized for a landscape-wide simulation and finally to forecast the dynamics of the treeline ecotone and especially the position of its northern boundary in a more detailed and reliable manner. Acknowledgements We thank our colleagues from the joint Russian-German project, namely Luidmila Pestryakova, Anatoly Nikolaev and Alexey Kolmogorov for their great support and valuable discussions. We also thank two anonymous reviewers, for their recommendations which greatly improved the manuscript and Cathy Jenks for proofreading. The PhD positions of Stefan Kruse and Mareike Wieczorek were funded by the Initiating and Networking Fund of the Helmholtz Association and by the German Research Foundation, respectively. Appendix. Figs. A1–A4

4.4. Conclusions Individual-based spatially-explicit modelling of larch forests in northern Siberia reveals that tree stands respond with a pronounced time-lag to climate change. Should global warming continue at its current pace, it is very likely that tree stands will intensively densify in the coming decades. It is now feasible to simulate the complexities of the larch ecosystem more realistically because all life-cycle stages of larches are explicitly represented in our newly-designed model, from seeds to seedlings and finally to trees, on continuous x,y coordinates. The model performs well for north to mid-treeline populations but overestimates populations at the southern edge, so we need to revisit the processes and parameterization of the model to improve its performance in the southern treeline ecotone. Nevertheless, it helps to improve our overall understanding of the reaction of the treeline in Siberia towards past and present changing climates. To better understand how the treeline might develop in the future, the model allows us to test local processes that might drive or prohibit stand development in a variety of habitats/areas

Fig. A1. The fraction of mature trees separated into four height classes with breaks at 2, 3 and 8 m shown as bars for all plots and as points for replicated individual plots at sites 2, 3, and, 5 in the vicinity of Khatanga, Taymyr Peninsula. The estimated logistic function and graph is presented in red.

118

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

Fig. A2. Basal diameter growth and age for the three study sites 2, 3 and 5 (rows from top to bottom) by estimated basal diameter-height (left) and age-basal diameter relationship (right).

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

Fig. A3. Dispersal distances as based on dispersal functions for Gaussian (grey) and fat-tailed exponential (darkgrey) dispersal.

Fig. A4. Sketch of parameterization process for the LAVESI model to obtain estimates for parameter values which were lacking detailed reference information.

119

120

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121

References Abaimov, A.P., Erkalov, A.V., Prokushkin, S.G., Matsuura, Y., 2000. The conservation and quality of gmelin larch seeds in cryolithic zone of Central Siberia. In: Proceedings of the Eighth Symposium on the Joint Siberian Permafrost Studies Between Japan and Russia in 1999, National Institute for Environmental Studies, pp. 3–9. Abaimov, A.P., 2010. Geographical distribution and genetics of Siberian larch species. In: Osawa, A., Zyryanova, O.A., Matsuura, Y., Kajimoto, T., Wein, R.W. (Eds.), Permafrost Ecosystems – Siberian Larch Forests. Springer, Dordrecht, The Netherlands, pp. 41–58. Andreev, A.A., Schirrmeister, L., Tarasov, P.E., Ganopolski, A., Brovkin, V., Siegert, C., Wetterich, S., Hubberten, H.-W., 2011. Vegetation and climate history in the Laptev Sea region (Arctic Siberia) during Late Quaternary inferred from pollen records. Quat. Sci. Rev. 30, 2182–2199. Barbeito, I., Dawes, M.A., Rixen, C., Senn, J., Bebi, P., 2012. Factors driving mortality and growth at treeline: a 30-year experiment of 92 000 conifers. Ecology 93, 389–401. Barber, V.A., Juday, G.P., Finney, B.P., 2000. Reduced growth of Alaskan white spruce in the twentieth century from temperature-induced drought stress. Nature 405, 668–673. Berner, L.T., Beck, P.S.A., Bunn, A.G., Goetz, S.J., 2013. Plant response to climate change along the forest-tundra ecotone in northeastern Siberia. Global Change Biol. 19, 3449–3462. Bonan, G.B., 2008. Forests and climate change: forcings, feedbacks, and the climate benefits of forests. Science 320, 1444–1449. Brazhnik, K., Shugart, H.H., 2015. 3D simulation of boreal forests: structure and dynamics in complex terrain and in a changing climate. Environ. Res. Lett. 10, 105006. Brazhnik, K., Shugart, H.H., 2016. SIBBORK: a new spatially-explicit gap model for boreal forest. Ecol. Modell. 320, 182–196. Bugmann, H., 2001. A review of forest gap models. Clim. Change 51, 259–305. Cariboni, J., Gatelli, D., Liska, R., Saltelli, A., 2007. The role of sensitivity analysis in ecological modelling. Ecol. Modell. 203, 167–182. Czaran, T., Bartha, S.S., 1992. Spatiotemporal dynamic models of plant populations and communities. Trends Ecol. Evol. 7, 38–42. DeAngelis, D.L., Mooij, W.M., 2005. Individual-based modeling of ecological and evolutionary processes. Annu. Rev. Ecol. Evol. Syst. 36, 147–168. Dolman, A.J., Maximov, T.C., Moors, E.J., Maximov, A.P., Elbers, J.A., Kononov, A.V., Waterloo, M.J., van der Molen, M.K., 2004. Net ecosystem exchange of carbon dioxide and water of far eastern Siberian larch (Larix dahurica) on permafrost. Biogeosci. Discuss. 1, 275–309. Epstein, H.E., Yu, Q., Kaplan, J.O., Lischke, H., 2007. Simulating future changes in arctic and subarctic vegetation. Comput. Sci. Eng. 9, 12–23. Esper, J., Schweingruber, F.H., 2004. Large-scale treeline changes recorded in Siberia. Geophys. Res. Lett. 31, L06202. Franz, H.-J., 1973a. Die natürlichen zonen der sowjetunion. In: Physische Geographie Der Sowjetunion. VEB Hermann Haack, Geogr.-Kartogr. Anstalt Gotha/Leipzig, pp. 20–21. Franz, H.-J., 1973b. Das nordsibirische tiefland. In: Physische Geographie Der Sowjetunion. VEB Hermann Haack, Geogr.-Kartogr. Anstalt Gotha/Leipzig, pp. 282–285. Fyllas, N.M., Politi, P.I., Galanidis, A., Dimitrakopoulos, P.G., Arianoutsou, M., 2010. Simulating regeneration and vegetation dynamics in Mediterranean coniferous forests. Ecol. Modell. 221, 1494–1504. Golubeva, E., Hofgaard, A., Silenchuk, K., 2013. The morphometric structure of the Larix gmelinii recruitment at the northern limit of its range in the forest-tundra ecotone. Geogr. Environ. Sustain. 6, 86–92. Grimm, V., Railsback, S.F., 2005. Individual-based Modeling and Ecology. Princeton University Press, Princeton, NJ. Grimm, V., Berger, U., Bastiansen, F., et al., 2006. A standard protocol for describing individual-based and agent-based models. Ecol. Modell. 198, 115–126. Grimm, V., Berger, U., DeAngelis, D.L., Polhill, J.G., Giske, J., Railsback, S.F., 2010. The ODD protocol: a review and first update. Ecol. Modell. 221, 2760–2768. Harris, I., Jones, P.D., Osborn, T.J., Lister, D.H., 2014. Updated high-resolution grids of monthly climatic observations – the CRU TS3.10 Dataset. Int. J. Climatol. 34, 623–642. Harsch, M.A., Hulme, P.E., McGlone, M.S., Duncan, R.P., 2009. Are treelines advancing? A global meta-analysis of treeline response to climate warming. Ecol. Lett. 12, 1040–1049. Hinkel, K.M., Nicholas, J.R.J., 1995. Active layer thaw rate at a boreal forest site in Central Alaska, U.S.A. Arct. Alp. Res. 27, 72–80. Holtmeier, F.-K., Broll, G., 2005. Sensitivity and response of Northern Hemisphere altitudinal and polar treelines to environmental change at landscape and local scales. Global Ecol. Biogeogr. 14, 395–410. Holtmeier, F.K., 2009. Definitions, terminology. Mountain Timberlines: Ecology, Patchiness, and Dynamics, vol. 36. Springer Science & Business Media, pp. 11–28. IPCC, 2014. Climate Change 2014: Impacts, Adaptation, and Vulnerability. Part A: Global and Sectoral Aspects. Contribution of Working Group II to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge. Jeltsch, F., Moloney, K.A., Schurr, F.M., Köchy, M., Schwager, M., 2008. The state of plant population modelling in light of environmental change. Perspect. Plant Ecol. Evol. Syst. 9, 171–189.

Kajimoto, T., Matsuura, Y., Osawa, A., Prokushkin, A.S., Mark, A., Abaimov, A.P., 2003. Root system development of Larix gmelinii trees affected by micro-scale conditions of permafrost soils in central Siberia. In: Jun, A.B.E. (Ed.), Roots: The Dynamic Interface Between Plants and the Earth. Springer, The Netherlands, pp. 281–292. Kajimoto, T., Matsuura, Y., Osawa, A., Abaimov, A.P., Zyryanova, O.A., Isaev, A.P., Yefremov, D.P., Mori, S., Koike, T., 2006. Size–mass allometry and biomass allocation of two larch species growing on the continuous permafrost region in Siberia. For. Ecol. Manage. 222, 314–325. Kajimoto, T., Osawa, A., Matsuura, Y., Abaimov, A.P., Zyryanova, O.A., Kondo, K., Tokuchi, N., Hirobe, M., 2007. Individual-based measurement and analysis of root system development: case studies for Larix gmelinii trees growing on the permafrost region in Siberia. J. For. Res. 12, 103–112. Kharuk, V.I., Ranson, K.J., Im, S.T., Naurzbaev, M.M., 2006. Forest-tundra larch forests and climatic trends. Russ. J. Ecol. 37, 291–298. Kharuk, V.I., Ranson, K.J., Dvinskaya, M.L., Im, S.T., 2011. Wildfires in northern Siberian larch dominated communities. Environ. Res. Lett. 6, 045208. Kharuk, V.I., Ranson, K.J., Im, S.T., Il’ya, A.P., 2015. Climate-induced larch growth response within the central Siberian permafrost zone. Environ. Res. Lett. 10, 125009. Kirdyanov, A.V., Treydte, K.S., Nikolaev, A., Helle, G., Schleser, G.H., 2008. Climate signals in tree-ring width, density and ␦13 C from larches in Eastern Siberia (Russia). Chem. Geol. 252, 31–41. Körner, C., 1998. A re-assessment of high elevation treeline positions and their explanation. Oecologia 115, 445–459. Laurance, W.F., Lovejoy, T.E., Vasconcelos, H.L., Bruna, E.M., Didham, R.K., Stouffer, P.C., Gascon, C., Bierregaard, R.O., Laurance, S.G., Sampaio, E., 2002. Ecosystem decay of Amazonian forest fragments: a 22-year investigation. Conserv. Biol. 16, 605–618. Lawrence, D.M., Koven, C.D., Swenson, S.C., Riley, W.J., Slater, A.G., 2015. Permafrost thaw and resulting soil moisture changes regulate projected high-latitude CO2 and CH4 emissions. Environ. Res. Lett. 10, 94011. Lenoir, J., Svenning, J.-C., 2015. Climate-related range shifts – a global multidimensional synthesis and new research directions. Ecography 38, 15–28. Levin, S.A., Muller-Landau, H.C., Nathan, R., Chave, J., 2003. The ecology and evolution of seed dispersal: a theoretical perspective. Annu. Rev. Ecol. Evol. Syst. 34, 575–604. Lloyd, A.H., Bunn, A.G., Berner, L., 2010. A latitudinal gradient in tree growth response to climate warming in the Siberian taiga. Global Change Biol. 17, 1935–1945. MacDonald, G.M., Kremenetski, K.V., Beilman, D.W., 2008. Climate change and the northern Russian treeline zone. Philos. Trans. R Soc. Lond. B Biol. Sci. 363, 2285–2299. Mamet, S.D., Kershaw, G.P., 2012. Multi-scale analysis of environmental conditions and conifer seedling distribution across the treeline ecotone of Northern Manitoba, Canada. Ecosystems 16, 295–309. Martínez, I., Wiegand, T., Camarero, J.J., Batllori, E., Gutiérrez, E., 2011. Disentangling the formation of contrasting tree-line physiognomies combining model selection and Bayesian parameterization for simulation models. Am. Nat. 177, E136–E152. Mathisen, I.E., Mikheeva, A., Tutubalina, O.V., Aune, S., Hofgaard, A., 2014. Fifty years of tree line change in the Khibiny Mountains, Russia: advantages of combined remote sensing and dendroecological approaches. Appl. Veg. Sci. 17, 6–16. Menne, M.J., Durre, I., Vose, R.S., Gleason, B.E., Houston, T.G., 2012. An overview of the Global Historical Climatology Network-daily database. J. Atmos. Oceanic Technol. 29, 897–910. Meyer, H., Opel, T., Laepple, T., Dereviagin, A.Y., Hoffmann, K., Werner, M., 2015. Long-term winter warming trend in the Siberian Arctic during the mid- to late Holocene. Nat. Geosci. 8, 122–125. Mitchell, T.D., Jones, P.D., 2005. An improved method of constructing a database of monthly climate observations and associated high-resolution grids. Int. J. Climatol. 25, 693–712. National Geophysical Data Center, 1988. ETOPO-5 Bathymetry/topography Data, Data Announcement 88-MGG-02, Digital Relief of the Surface of the Earth. NOAA, National Geophysical Data Center, Boulder, Colorado. Nathan, R., Muller-Landau, H., 2000. Spatial patterns of seed dispersal, their determinants and consequences for recruitment. Trends Ecol. Evol. 15, 278–285. Naurzbaev, M.M., Vaganov, E.A., Sidorova, O.V., Schweingruber, F.H., 2002. Summer temperatures in eastern Taimyr inferred from a 2427-year late-Holocene tree-ring chronology and earlier floating series. Holocene 12, 727–736. New, M., Hulme, M., Jones, P., 1999. Representing twentieth-century space–time climate variability. Part I: development of a 1961–90 mean monthly terrestrial climatology. Am. Meteorol. Soc. 12, 829–856. Osawa, A., Kajimoto, T., 2010. Development of stand structure in larch forests. In: Osawa, A., Zyryanova, O.A., Matsuura, Y., Kajimoto, T., Wein, R.W. (Eds.), Permafrost Ecosystems – Siberian Larch Forests. Springer, Dordrecht, The Netherlands, pp. 123–148. Pacala, S.W., Canham, C.D., Silander Jr, J.A., 1993. Forest models defined by field measurements: I. The design of a northeastern forest simulator. Can. J. For. Res. 23, 1980–1988. Pacala, S.W., Canham, C.D., Saponara, J., John, A.S.J., Kobe, R.K., Ribbens, E., 1996. Forest models defined by field measurements: estimation, error analysis and dynamics. Ecol. Monogr. 66, 1–43.

S. Kruse et al. / Ecological Modelling 338 (2016) 101–121 Pearson, R.G., Phillips, S.J., Loranty, M.M., Beck, P.S.A., Damoulas, T., Knight, S.J., Goetz, S.J., 2013. Shifts in Arctic vegetation and associated feedbacks under climate change. Nat. Clim. Change 3, 673–677. Peck, S.L., 2004. Simulation as experiment: a philosophical reassessment for biological modeling. Trends Ecol. Evol. 19, 530–534. Polezhaeva, M.A., Lascoux, M., Semerikov, V.L., 2010. Cytoplasmic DNA variation and biogeography of Larix Mill. in northeast Asia. Mol. Ecol. 19, 1239–1252. Rammig, A., Fahse, L., 2009. Simulating forest succession after blowdown events: the crucial role of space for a realistic management. Ecol. Modell. 220, 3555–3564. Ripley, B.D., 1977. Modelling spatial patterns. J. R. Stat. Soc. 39, 172–212. Sato, H., Itoh, A., Kohyama, T., 2007. SEIB-DGVM: a new Dynamic Global Vegetation Model using a spatially explicit individual-based approach. Ecol. Modell. 200, 279–307. Schulze, E.-D., Wirth, C., Mollicone, D., von Lüpke, N., Ziegler, W., Achard, F., Mund, M., Prokushkin, A., Scherbina, S., 2012. Factors promoting larch dominance in central Siberia: fire versus growth performance and implications for carbon dynamics at the boundary of evergreen and deciduous conifers. Biogeosciences 9, 1405–1421. Seidl, R., Rammer, W., Scheller, R.M., Spies, T.A., 2012. An individual-based process model to simulate landscape-scale forest ecosystem dynamics. Ecol. Modell. 231, 87–100. Semerikov, V.L., Iroshnikov, A.I., Lascoux, M., 2007. Mitochondrial DNA variation pattern and postglacial history of the Siberian Larch (Larix sibirica Ledeb.). Russ. J. Ecol. 38, 147–154. Shiyatov, S.G., Mazepa, V.S., 2012. Climate-driven dynamics of the forest-tundra vegetation in the Polar Ural Mountains. Contemp. Prob. Ecol. 4, 758–768. Shiyatov, S.G., Terent’ev, M.M., Fomin, V.V., 2005. Spatiotemporal dynamics of forest-tundra communities in the Polar Urals. Russ. J. Ecol. 36, 69–75. Shugart, H.H., West, D.C., 1979. Notes: size and pattern of simulated forest stands. For. Sci. 25, 120–122. Shuman, J.K., Shugart, H.H., 2009. Evaluating the sensitivity of Eurasian forest biomass to climate change using a dynamic vegetation model. Environ. Res. Lett. 4, 045024. Sidorova, O.V., Siegwolf, R.T.W., Saurer, M., Naurzbaev, M.M., Shashkin, A.V., Vaganov, E.A., 2010. Spatial patterns of climatic changes in the Eurasian north reflected in Siberian larch tree-ring parameters and stable isotopes. Global Change Biol. 16, 1003–1018. Sidorova, O.V., Siegwolf, R.T.W., Saurer, M., Shashkin, A.V., Knorre, A.A., Prokushkin, A.S., Vaganov, E.A., Kirdyanov, A.V., 2009. Do centennial tree-ring and stable isotope trends of Larix gmelinii (Rupr.) Rupr. indicate increasing water shortage in the Siberian north? Oecologia 161, 825–835.

121

Sitch, S., Smith, B., Prentice, I.C., Arneth, A., Bondeau, A., Cramer, W., Kaplan, J.O., Levis, S., Lucht, W., Sykes, M.T., Thonicke, K., Venevsky, S., 2003. Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model. Global Change Biol. 9, 161–185. Sitch, S., Huntingford, C., Gedney, N., Levy, P.E., Lomas, M., Piao, S.L., Betts, R., Ciais, P., Cox, P., Friedlingstein, P., Jones, C.D., Prentice, I.C., Woodward, F.I., 2008. Evaluation of the terrestrial carbon cycle, future plant geography and climate-carbon cycle feedbacks using five Dynamic Global Vegetation Models (DGVMs). Global Change Biol. 14, 2015–2039. Smith, B., Prentice, I.C., Sykes, M.T., 2001. Representation of vegetation dynamics in the modelling of terrestrial ecosystems: comparing two contrasting approaches within European climate space. Global Ecol. Biogeogr. 10, 621–637. Smith, T.M., Urban, D.L., 1988. Scale and resolution of forest structural pattern. Vegetatio 74, 143–150. Sofronov, M.A., Volokitina, A.V., Kajimoto, T., Matsuura, Y., Uemura, S., 2000. Zonal peculiarities of forest vegetation controlled by fires in Northern Siberia. Eurasian J. For. Res. 1, 51–57. Stocker, B.D., Roth, R., Joos, F., Spahni, R., Steinacher, M., Zaehle, S., Bouwman, L., Prentice, I.C., 2013. Multiple greenhouse-gas feedbacks from the land biosphere under future climate change scenarios. Nat. Clim. Change 3, 666–672. Sugimoto, A., Yanagisawa, N., Naito, D., Fujita, N., Maximov, T.C., 2002. Importance of permafrost as a source of water for plants in east Siberian taiga. Ecol. Res. 17, 493–503. Walker, D.A., Raynolds, M.K., Daniëls, F.J.A., et al., 2005. The circumpolar Arctic vegetation map. J. Veg. Sci. 16, 267–282. Wullschleger, S.D., Post, W.M., King, A.W., 1995. On the potential for a CO2 fertilization effect in forests:estimates of the biotic growth factor based on 58 controlled-exposure studies. In: Woodwell G.M., Mackenzie F.T. (Eds.), Biotic Feedbacks in the Global Climatic System, pp. 85–107. Yu, Q., Epstein, H., Walker, D., 2009. Simulating the effects of soil organic nitrogen and grazing on Arctic tundra vegetation dynamics on the Yamal Peninsula, Russia. Environ. Res. Lett. 4, 045027. Zhang, T., Frauenfeld, O.W., Serreze, M.C., Etringer, A., Oelke, C., McCreight, J., Barry, R.G., Gilichinsky, D., Yang, D., Ye, H., et al., 2005. Spatial and temporal variability in active layer thickness over the Russian Arctic drainage basin. J. Geophys. Res.: Atmos. 110, D16101. Zhang, W., Miller, P.A., Smith, B., Wania, R., Koenigk, T., Döscher, R., 2013. Tundra shrubification and tree-line advance amplify Arctic climate warming: results from an individual-based dynamic vegetation model. Environ. Res. Lett. 8, 034023.