Identifying important topics for model refinement in a widely used process-based model informed by correlative model analyses in a boreal forest

Models attempting to predict treeline shifts in changing climates must include the relevant ecological processes in sufficient detail. A previous correlative model study has pointed to nutrients, competition, and temperature as the most important factors shaping the treelines of Pinus sylvestris L., Picea abies (L.) H. Karst. and Betula pubescens Ehrh. in Finnish Lapland. Here, we applied a widely used process-based dynamic vegetation model (LPJ-GUESS) to (i) test its capability to simulate observed spatial and temporal patterns of the main tree species in Finnish Lapland, and (ii) to explore the model representation of important processes in order to guide further model development. A European parameterization of LPJ-GUESS overestimated especially P. abies biomass and the species’ northern range limit. We identified implemented processes to adjust (competition, disturbance) and crucial processes in boreal forests to include (nutrient limitation, forest management) which account for the model’s failure to (edaphically) restrict P. abies in Finnish Lapland and the resulting species imbalance. Key competitive mechanisms are shade and drought tolerance, nutrient limitation, fire resistance, and susceptibility to disturbances (storm, herbivory) which we discussed with respect to boreal ecology and promising model developments to provide a starting point for future model development.


Introduction
Arctic treelines are the (more or less) well-defined biome boundaries between dense forest and tundra which have shifted in the past and will continue to shift in the future.Today's spatial pattern of treelines is the combined result of historic developments and current processes (Payette et al. 2002;Seppä et al. 2002).Recent climate change has led to a rise in mean annual and especially winter temperatures in Finland between 1901 and 2000 (Jylhä et al. 2004).The effect of these climatic changes on tree growth is not straightforward.Rising temperature sums can improve growing conditions of established trees (Sveinbjörnsson et al. 2002;Danby 2007;Moen et al. 2008) and reduce seedling mortality due to consistently milder winter temperatures (Kullman 1997).However, non-climatic factors can limit a possible treeline advance, e.g.competition by shrubs (Nilsson et al. 1993;Weih and Karlsson 1999), herbivore pressure by reindeer (Anschlag et al. 2008;Aakala et al. 2014) and adverse soil conditions (Holtmeier et al. 2003;Sutinen et al. 2005).In the south of Finnish Lapland, the Tanaelv and Lapland Greenstone Belts form moist nutrient-rich soils dominated by Picea abies (L.) H. Karst.stands (Sutinen et al. 2005).They are bordered to the north (>68°N) by the Lapland Granulite Belt (Cagnard et al. 2011) which is at most 80 km wide.Its dry nutrient-poor soils dominated by Pinus sylvestris L. stands act as a dispersal barrier for P. abies (Sutinen et al. 2005).This could explain why P. abies is restricted to the south of Lapland while surviving in isolated outposts (natural and planted) far north of its current treeline (Oksanen 1995).
Numerous studies document elevational treeline advances due to climate warming in the 20th century in Finnish Lapland (Juntunen et al. 2002;Middleton et al. 2008;Aakala et al. 2014;Franke et al. 2015), and the Swedish Scandes (Kullman 1997;Kullman and Öberg 2009).However, non-climatic factors not affected by climate change prevent or at least slow down treeline advance in many places (Dalen and Hofgaard 2005;Holtmeier 2011).Concerning latitudinal treeline advances, there are far fewer studies available (see Harsch et al. (2009) for a global meta-analysis), and these studies report a slower-than-expected response of latitudinal treelines to climate change (e.g. in Russia (MacDonald et al. 2008), Norway (Hofgaard et al. 2013) and Canada (Masek 2001)).
Here, we followed up on results from a correlative model study that investigated processes determining the current treeline position of P. sylvestris, P. abies and Betula pubescens Ehrh. in Finnish Lapland (Schibalski et al. 2014).Despite the warming over the study period of 25 years, the underlying forest inventory data of 1978 and 2003 showed no clear latitudinal treeline advance.Site fertility, abundance of co-occurring species and growing degree days (GDD) were the most important predictors explaining the three tree species' occurrence (presence-absence) and basal area distribution (Schibalski et al. 2014).This study thus supported the critical role non-climatic factors can play in limiting tree line advance.The interesting questions of climatic vs. edaphic limitation, the role of competition, and the observed delay in the species response to recent climate change make Finnish Lapland a suitable case study for a second modeling approach.Thus, we applied the established and widely used process-based dynamic vegetation model LPJ-GUESS (Smith et al. 2001; http://iis4.nateko.lu.se/lpj-guess/) to predict the ranges and biomass of P. sylvestris, P. abies and B. pubescens in the same region and over the same time period as the correlative model study.
Our aims were (i) to test the general capability of LPJ-GUESS to simulate the spatial and temporal biomass pattern and ranges (treelines) of the three main tree species in Finnish Lapland, and (ii) to explore the representation of competition, climatic and edaphic factors in LPJ-GUESS, thereby revealing potential shortcomings and thus guiding further model development.To assess model performance (i), we compared the biomass patterns and range limits simulated by LPJ-GUESS with observed biomass (forest inventory data, 2011).To explore process representation in LPJ-GUESS (ii), we analyzed the parameterization of currently implemented processes with respect to the ecology of boreal forests and reviewed additional process implementations in other existing LPJ-GUESS versions.

General model description
LPJ-GUESS (Lund-Potsdam-Jena General Ecosystem Simulator) is a flexible biome-scale model for simulating vegetation biogeography and dynamics, as well as biogeochemical cycles at regional to global scales.It shares many ecophysiological process-representations with the widely used Lund-Potsdam-Jena Dynamic Global Vegetation Model, LPJ-DGVM (Smith et al. 2001;Sitch et al. 2003).But vegetation dynamics and vegetation structure are simulated at a greater level of detail, allowing the parameterization of individual species as opposed to broader plant functional types.Vegetation dynamics are simulated as the emergent outcome of growth and competition for light, space and soil resources among woody plant individuals and a herbaceous understorey based on their functional traits.Plant-physiological processes like photosynthesis and respiration, as well as the exchange of carbon and water between vegetation, soil, and atmosphere, are simulated on a daily basis.Vegetation growth, biomass allocation, establishment, and mortality are simulated once at the end of a simulation year.Tree mortality occurs (i) if the five-year average growth efficiency, calculated from annual net primary production, leaf carbon mass and specific leaf area (Smith et al. 2001), falls below a species-specific minimum growth efficiency (greff, in Supplementary file S1/ Table S1 (available at https://doi.org/10.14214/sf.6977)),(ii) as trees reach their maximum age, and additionally as a result of (iii) fire and (iv) a stochastic patch-destroying disturbance which recurs within an expected mean interval of, here, 200 years.This patch-destroying disturbance kills all trees in a patch and represents rare events such as pest calamities or windstorms (see Hickler et al. (2012) for more details).Wildfires are modelled based on temperature, fuel (litter) load and moisture (Thonicke et al. 2001) and affect trees according to their species-specific fire resistance (r fire , Suppl.file S1/Table S1).
Vegetation dynamics are simulated in each of a number (50 in this study) of replicate patches (0.1 ha) representing "random samples" of each simulated locality or grid cell.Each model grid cell is homogeneous in terms of soil texture, atmospheric CO 2 concentration and a set of climatic variables (daily temperature, precipitation, and radiation).Its size is determined by the spatial resolution of this input data (10 × 10 km in this study).Multiple patches are simulated to account for the variability within a landscape representative of the grid cell, as vegetation stands differ in their histories of disturbance (which are assumed to occur independently of neighbouring patches or grid cells) and stand development (succession).The output from individual patches is averaged to characterize the average vegetation per grid cell.
In this study, we used the version parameterized for major European tree species and plant functional types by Hickler et al. (2012), with an additional species-specific water supply function (Schurgers et al. 2011).Pre-defined bioclimatic limits determine the species which can establish and survive in a model grid cell.This fact makes LPJ-GUESS a "fitted" process-based model according to Dormann et al. (2012) -as opposed to "forward" process-based models without calibration.One of these bioclimatic limits is "minimum growing degree days for establishment" (GDD 5,min , Suppl.file S1/Table S1) which we can directly compare to statistical thresholds found in the correlative model study (Schibalski et al. 2014).

Species characterization in LPJ-GUESS
In LPJ-GUESS, the simulated (tree) species are described by leaf or needle functional traits, leaf area to sapwood cross-sectional area ratio, phenology, fire resistance, root distribution, bioclimatic limits for establishment and survival (minimum monthly winter temperature), as well as life history strategy (related to shade tolerance, see below).Establishment depends on minimum GDD (5 °C), maximum monthly winter temperature (representing the chilling requirement of northern species) as well as minimum plant-available water content of the upper soil layer during the growing season.The latter also influences the species-specific water supply function, with more water available for a given soil water content for species with a lower limit (Schurgers et al. 2011).All parameters are listed in Suppl.file S1/Table S1 (cf.Hickler et al. 2012).The simulations of this study were carried out in "cohort mode" in which cohorts of individuals recruited in the same patch in a given year are represented by a single average individual and are thus assumed to retain the same size and form as they grow.
In LPJ-GUESS, shade tolerance defines an important trade-off during succession: Shadeintolerant species like B. pubescens require more light for establishment (par min , Suppl.file S1/ Table S1) than shade-tolerant species.Shade-intolerant species also have higher maximum establishment rates (est max , Suppl.file S1/Table S1) under high-light conditions, but establishment rates rapidly decline as the canopy closes and less radiation reaches the forest floor (α, Suppl.file S1/ Table S1).They also suffer more from growth-efficiency mortality (greff, Suppl.file S1/Table S1) as the canopy closes and growth is diminished due to increasing competition for light.However, as a result of higher sapwood to heartwood conversion (turn sapwood , Suppl.file S1/Table S1), shadeintolerant species grow faster under high-light conditions.For a full description of the associated equations, see Hickler et al. (2012) and Smith et al. (2001).The associated parameters were fitted to yield realistic succession patterns in selected European forests, but not including sites from northern Scandinavia (Hickler et al. 2012).Silva Fennica vol. 51 no. 4

Environmental input data and setup of model runs
As environmental input data, LPJ-GUESS requires daily mean air temperature, precipitation sum, radiation, atmospheric CO 2 and soil texture.We used soil data from the 7th and 9th National Forest Inventory in 1978 (Kuusela and Salminen 1991) and 2003 (Tomppo et al. 2011), also used in Schibalski et al. (2014), to assign each plot one of the nine soil classes in LPJ-GUESS which differ in terms of water holding capacity and thermal diffusivity (Sitch et al. 2003; Table 4).In our study region, medium textures dominate (70%), but there are organic soils in the southern part of Finnish Lapland (21%, Suppl.file S1/Fig.S1).
For regional climate input, we used monthly mean and minimum temperature, precipitation and radiation in an interpolated 10 × 10 km grid from 1961 to 2003 (Venäläinen et al. 2005), and linear interpolation between monthly values to construct the daily inputs.In contrast, atmospheric CO 2 was given as annual averages, not further regionalized (Suppl.file S1/Fig.S2).From 1978 to 2003, mean monthly temperatures have increased significantly (p < 0.001, Wilcoxon rank sum test) in all months except June and December (Fig. 1; Suppl.file S1/Fig.S3).Similarly, growing  (1968)(1969)(1970)(1971)(1972)(1973)(1974)(1975)(1976)(1977)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002); the mean temperature of the coldest month was averaged over 17 years preceding the simulation year (1961-1977 and 1986-2002, respectively).Plus signs indicate median (intersection) and standard deviation (length of the arms).Solid signs/plain text mean that values were significantly higher in 2003 than 1978, broken signs/italic text means the opposite (p < 0.001, Wilcoxon rank sum test); colour in the upper two panels is used to clarify label assignment.For growing degree days and coldest month mean temperatures, individual grid cell values are shown in addition to their median and standard deviation.See Supplementary file S1/Fig.S3 for more detailed information on monthly mean temperature and precipitation sums.degree days have increased, but we found spatial differences across Finnish Lapland with decreases in some areas (Suppl.file S1/Fig.S4).Monthly precipitation sums have increased for all months but September over the 25 years (Fig. 1; Suppl.file S1/Fig.S3).
LPJ-GUESS grows vegetation from bare soil.To reach approximate equilibrium conditions, we let the model run for 1000 years before the actual simulation period .As input data for this spin-up, we recycled the oldest 30 years of historical climate data (with detrended temperatures).
To assess the model's capability to simulate the spatial and temporal biomass patterns and ranges (treelines) of the three main tree species in Finnish Lapland -aim (i) -, we ran the model with all three species together (called "multi-species" hereafter), thus including biotic interaction.We compared above-and belowground biomass [kg m -²] per species and grid cell with observed biomass data (described in section 2.2).In addition, we ran the model separately for each species alone, i.e. without the competition of the other two species (called "single-species").We were thus able to assess the influence of interspecific competition in LPJ-GUESS -aim (ii) -and gain insight into the species' performance independent of competing species.

Multi-Source National Forest Inventory data
We compared LPJ-GUESS biomass estimates from the multi-species run with biomass estimates from the Multi-Source National Forest Inventory (MS-NFI, Tomppo et al. 2008) for Finnish Lapland.These biomass estimates are a combination of field observations and satellite imagery from 2011, publicly available online (http://kartta.luke.fi/index-en.html).MS-NFI biomass estimates are provided as separate biomass components [10 kg ha -1 ]: living and dead branches, roots, stump, stem with bark and stem residual, as well as foliage or needles.To directly compare LPJ-GUESS results with the MS-NFI data, we added up the single biomass components for P. sylvestris, P. abies and broad-leaved trees (including B. pubescens).

Results
The total biomass, i.e. the biomass sum of all three species, was overestimated by LPJ-GUESS (Fig. 2a, observed and simulated biomass) in Finnish Lapland.However, the spatial trend of northwards decreasing biomass observed in the MS-NFI data was reproduced by LPJ-GUESS.
For P. sylvestris, the biomass range matched between the LPJ-GUESS output (multi-species run) and MS-NFI data (Fig. 2a), except for the far north (>69°N) where the LPJ-GUESS biomass predictions were too high.The spatial pattern of high and low biomass was not reproduced correctly as the simulated biomass increased towards the north, while the observed biomass actually decreases towards the treeline (Fig. 3b).Without the competition of the other two species (single-species run), the simulated biomass was much higher, obviously exceeding the observed values, but the spatial pattern of northwards decreasing biomasses was correctly captured (Fig. 2b).
For B. pubescens, we found a similar pattern: the range of biomass was similar between the LPJ-GUESS output (multi-species) and MS-NFI data (Fig. 2a), especially when taking into account that MS-NFI data comprised all deciduous species.In the far north, where no other deciduous species prevail, the MS-NFI estimate equalled B. pubescens biomass, and the match between LPJ-GUESS simulations and MS-NFI observations was good.In the south, LPJ-GUESS underestimated the MS-NFI data which includes other deciduous species coexisting with B. pubescens.Again, the correct spatial trend of northwards decreasing biomass in the single-species model run was effectively reversed when including competition (Fig. 3d).In Finnish Lapland, the range limit of B. pubescens is much less distinct than the two conifers' clear treelines, which was reflected by both observed and simulated biomass (Fig. 3d).
Finally, P. abies was greatly overestimated in both biomass range (Fig. 2a) and species range (LPJ-GUESS did not capture the distinct treeline at ~68.5°N, Fig. 3c).Although LPJ-GUESS simulated a decrease in biomass towards the north, the range limit of P. abies in the model was not reached and is obviously far north of the observed treeline (Fig. 3c).In contrast to the other two species, the multi-species and single-species model runs yielded virtually the same results for P. abies.(averaged over 1994-2003).To maximize visibility of spatial differences but retain comparability between observations and simulations, we used quantiles to define classes for each species and the total.This results in the irregular class spacing and reflects the different biomass distributions (cf.Fig. 2)."Natural" colours (white to black) cover the range of the observed values (i.e. the upper limit of the black class is always the maximum of the respective MS-NFI data); "artificial" colours (shades of magenta) cover the predictions that exceed the observed value range (model overestimation).
Simulated biomass increases from 1978 to 2003 were distributed throughout Lapland for P. abies and B. pubescens (Fig. 4); they were not associated with a treeline advance.In contrast, simulated biomass increases of P. sylvestris were concentrated in the far north of our study region (Fig. 4) with the greatest biomass increase (2.06 kg m -2 ) at 69.9°N.
We compared growing degree day thresholds in statistical response curves from the previous modelling study (Schibalski et al. 2014) with LPJ-GUESS parameter GDD 5,min (Suppl.file S1/ Table S1).For P. abies, the parameter value of 600 GDD fit the observed data very well -at least for 1978 (Fig. 5).For P. sylvestris, the parameter value of GDD 5,min is 500 GDD (Suppl.file S1/ Table S1), which was much lower than any threshold (1978 or 2003) in the observations (Fig. 5).

Total biomass overestimation
Total biomass was overestimated, which is in accordance with other studies applying LPJ-GUESS to Scandinavia.Smith et al. (2008) found LPJ-GUESS to overestimate conifer forest biomass, leaf area index and tree density in northern Fennoscandia unless the model was constrained by satellite data.As in similar dynamic vegetation models (Thornton et al. 2007;Sokolov et al. 2008;Zaehle et al. 2010), however, primary production, a key driver of the simulated biomass, in LPJ-GUESS is substantially lower in northern forests when accounting for nitrogen limitation compared to the unlimited model version (Smith et al. 2014).Including nitrogen cycling in LPJ-GUESS reduced the overestimation of gross primary production from 56% to 18% in boreal forests (Fleischer et al. 2015).This confirms the general assumption that forest growth in the region is heavily limited by soil nutrients, in particular, nitrogen (Vitousek and Howarth 1991;Lupi et al. 2013).Although nitrogen limitation has been implemented within the global LPJ-GUESS version based on broader plant functional types (Smith et al. 2014), these developments have not yet been combined with regional tree species parameterization (see also discussion of species-specific nutrient limitation in section 4.3.4).
Another process that could potentially reduce total biomass to the observed level is disturbance.In LPJ-GUESS, patch-destroying disturbances are a stochastic process, determined by the disturbance interval (parameter), which destroys all biomass in a patch.Our return interval for patch-destroying disturbance of 200 years, which was adopted as an average across Europe (Hickler et al. 2012), is probably too long.Increasing the disturbance frequency would effectively restrict biomass accumulation, especially that of the slower growing conifers.However, the susceptibility to disturbances in Finnish Lapland differs between species (see section 4.3.6 below).This should be accounted for to correct not only total biomass levels in the single-species model runs but also spatial patterns influenced by inter-specific competition in the multi-species model run.1968-1977 and 1993-2002; given as rug plot in the upper panels: top 1978, bottom 2003).Vertical lines mark the bioclimatic limit used in LPJ-GUESS for the respective species (GDD 5,min , Suppl.file S1/Table S1).Transparent bootstrapped confidence bands (0.95) were calculated following the procedure detailed in Coutts (2011) and Coutts and Yokomizo (2014), using the boot.cifunction in R (Canty and Ripley 2013).Note: The prevalence of Picea abies was very low in the 1978 dataset (7%) leading to the excessive bootstrapped confidence band (right).
The omission of forest management in our version of LPJ-GUESS (however, see Jönsson et al. 2015) surely contributes to the overestimation of the conifer biomass.Forest management in Finnish Lapland is limited mainly to the southern part (south of about 68°N) and the Lake Inari region.It consists of thinning and clear-cutting, as well as e.g.preparation of regeneration areas (Ylitalo 2013).Rotation times in Lapland range between 60 and 150 years compared to 40 to 100 years in the south of Finland (Äijälä et al. 2014).Thus, including region-specific management measures in the model could alleviate the biomass overestimation, especially where forests are used more intensively.

Treeline dynamics and growing degree days
In accordance with the findings of Schibalski et al. (2014), the spatial pattern of simulated biomass increase from 1978 to 2003 did not indicate a treeline advance of P. abies (B.pubescens lacks a clearly defined treeline in Finnish Lapland).For P. sylvestris, however, simulated biomass increases concentrated in the far north of our study region (Fig. 4).Indeed, field studies report improved growing conditions for P. sylvestris due to climate change and a potential for elevational treeline advance in Finland (Juntunen et al. 2002;Aakala et al. 2014;Franke et al. 2015), while others contest the idea of a rapid P. sylvestris treeline advance and stress other limiting factors (Holtmeier and Broll 2005;Kullman 2010;Holtmeier 2011).To our knowledge, however, there is only one study reporting a northward advance of P. sylvestris treeline due to recent climate change in Scandinavia (Hofgaard et al. 2013).The northward advance for P. sylvestris reported in Hofgaard et al. (2013) was very modest (10 m year -1 ) compared to B. pubescens treeline advance (340 m year -1 ).
Moreover, Hofgaard et al. (2013) stress the discrepancy between the empirical evidence of treeline advance and the often dramatic predictions of vegetation changes by DGVMs.Dynamic global vegetation models like LPJ-GUESS are suspected to overestimate the climate-sensitivity of treelines and thus climate change-driven treeline advance (Loehle and LeBlanc 1996;Hofgaard et al. 2009).While our application of LPJ-GUESS spanning merely 25 years did not predict a clear treeline advance (in line with empirical evidence), a long-term modelling study applying LPJ-GUESS simulated complete boreal forest coverage of Finnish Lapland by 2090 (Wolf, Callaghan et al. 2008).
In our application, LPJ-GUESS overestimated P. abies beyond its current treeline, although the parameter value of GDD 5,min matched the observed data (Fig. 5), indicating that temperature limitation is not what keeps P. abies from occupying the far north of Finnish Lapland.For P. sylvestris, correlative GDD thresholds were much higher than the GDD 5,min parameter value in LPJ-GUESS.Increasing GDD 5,min from 500 to 625 GDD (suggested by the response curves in Fig. 5) should efficiently restrict P. sylvestris in the north (Fig. 6).Statistical fine-tuning as suggested here has the potential to improve LPJ-GUESS parameterization for a particular time (or place, e.g.Pappas et al. 2015).However, as we can already see from the climate change over the 25 years, this correlation changes over time and parameters would need to be adjusted again to effectively restrain the species in the model.Here, "fitted" process-based models like LPJ-GUESS underlie the same equilibrium assumptions as do correlative models (Dormann et al. 2012).They are also subject to the same problems when these models are applied to ongoing climate change, and nonstationarity violates the underlying assumptions (or fitted parameters).Snell et al. (2014), who advocate using DGVMs such as LPJ-GUESS to simulate range shifts, are aware of this issue.They propose Bayesian methods for parameterization (van Oijen et al. 2005;Hartig et al. 2012) and point to "next-generation DGVMs" (Scheiter et al. 2013) which simulate plant individuals with potentially unique trait combinations.

Species imbalance
There was virtually no effect of the presence of the other two species on P. abies in the multi-species LPJ-GUESS run.The biomass values and spatial patterns of P. sylvestris and B. pubescens, however, were distorted by P. abies presence in the model.In the single-species run, P. sylvestris biomass decreased gradually towards the north, correctly indicating that the species slowly approached its range limit due to unfavourable growing conditions in the model, albeit further north compared to the observations (Fig. 2).In the multi-species run, however, highest P. sylvestris biomass was found in the north coinciding with the lowest P. abies (and B. pubescens) biomass.This suggests that P. abies was by far too competitive in the model and that competition plays a pivotal role in LPJ-GUESS.
The model's inability to correctly reproduce the occurrence pattern of especially P. abies in Northern Lapland is in accordance with recent studies applying LPJ-GUESS: in a Holocene vegetation reconstruction study, Miller et al. (2008) were not able to correctly model P. abies' occurrence in Finland and Sweden over time.In their simulations of the current treeline in Arctic Europe, Fang et al. (2013) found that LPJ-GUESS did capture the coniferous treeline, but failed to correctly predict species-specific treelines.P. abies occurred north of its observed treeline where it suppressed P. sylvestris as additional simulation experiments showed.Fang et al. (2013) attributed the competitive strength of P. abies to its shade tolerance.

Shade tolerance
Competition for light is crucial in closed-canopy forests which were predicted by our simulations (incorrectly in the far north).Shade tolerance-related parameters in LPJ-GUESS include minimum light requirement for establishment, maximum establishment rate and growth-efficiency-related mortality (Suppl.file S1/Table S1).Wramneby et al. (2008) demonstrated that LPJ-GUESS is highly sensitive to shade tolerance-related parameters and that unfortunately, the uncertainty of these parameters is very large.Pinus sylvestris is ranked intermediate shade-tolerant and is thus trumped by both competitor species, i.e.P. abies which tolerates shaded conditions as well as B. pubescens Fig. 6.Maps of growing degree days for 1978 (1968-1977) and 2003 (1993-2002) with treelines of Pinus sylvestris (solid) and Picea abies (broken).Treelines are defined as the marginal sites occupied by the respective species and had not changed between 1978 and 2003.
which benefits most efficiently from light abundance after disturbances.Shaded conditions are more probable in our case because total biomass was overestimated and thus shading must be greater than observed.This as well as the low disturbance frequency discussed in section 4.1 reduces the competitive advantage of B. pubescens.Consequently, P. abies effectively distorts P. sylvestris and B. pubescens biomass distribution in the multi-species run, resulting in LPJ-GUESS's failure to correctly simulate the species balance observed in Finnish Lapland.We thus concur with Wramneby et al. (2008) in that shade tolerance is a very important trait in LPJ-GUESS which can dominate over other physiological differences between species (Suppl.file S1/Table S1).
In the following, we discuss competitive advantages that P. sylvestris and B. pubescens might have over P. abies and why the two species apparently fail to play off their strengths in our simulations.Some processes are already implemented in the LPJ-GUESS version used for this study but need to be modified for our application, while others require further model development.

Drought tolerance
P. sylvestris outcompetes P. abies on dry, acidic, nutrient-poor sites along the Lapland Granulite Belt as known from field observations (Sutinen et al. 2005) and experiments (Ingestad 1979).Accordingly, dry conditions should favour P. sylvestris in LPJ-GUESS as it has 40% of its roots distributed in lower soil layers compared to 20% for P. abies and B. pubescens.Pinus sylvestris is thus able to take up more water at low soil moisture contents (see water uptake function, Suppl.file S1/Fig.S5) and requires less soil moisture for establishment (awc min , Suppl.file S1/Table S1).Soil texture in LPJ-GUESS influences the water holding capacity and thermal diffusivity of a soil, but Wolf et al. (2008) showed that the vegetation outcome is rather insensitive to different soil moisture and soil temperature representations (cf.Pappas et al. 2013).LPJ-GUESS is, however, very sensitive to changes in soil depth.As no soil depth data was available for the study region, a uniform soil depth of 1.5 m was assumed.However, shallower soils do exist, at least locally.This simplification thus weakens P. sylvestris' advantage due to drought tolerance and fails to edaphically restrict P. abies along the Lapland Granulite Belt.

Nutrient limitation
In order for nutrient limitation to affect the competitive strength of individual species in LPJ-GUESS, it would need to be implemented species-specifically.It is questionable, however, whether we have enough process understanding to include species-specific responses to soil nutrients in a process-based framework.
Mixed in with issues of soil fertility are also management decisions: on dry and nutrient poor sites, forest managers will favour the superior species P. sylvestris by actively thinning P. abies and B. pubescens, which are less successful on these sites anyway (Äijälä et al. 2014).This positive feedback increases the competitive advantage of P. sylvestris over P. abies and B. pubescens where forests are managed (but see discussion of forest management in section 4.1).

Fire
P. sylvestris is well adapted to survive moderate fires by its thick bark (heat insulation), high crown base (preventing crown scorching) and deep root system (Fernandes et al. 2008).This is collectively reflected in LPJ-GUESS by a four times greater probability to survive fires of P. sylvestris compared to P. abies and B. pubescens (r fire , Suppl.file S1/Table S1).Thus, frequent fires could favour P. sylvestris over P. abies in our simulations.However, since the beginning of the 20th century, anthropogenic fire suppression in Fennoscandia has greatly extended the interval of forest fires (Zackrisson 1977).Additionally, the mean fire interval increases from 20 years in the south of Finland (58°N, 1401Finland (58°N, -1998) ) to more than 500 years in the north (69°N, 1400(69°N, -2001(69°N, , Larjavaara et al. 2005)).Thus, forest fires in Finland today are infrequent, small (the mean burnt area for the whole of Finland was 537 ha in 1994-2003, Ylitalo 2013)), and no longer play an important role in forest ecology.In line with this, fires did not play an important role in our simulations, as LPJ-GUESS underestimates natural fire cycles in northern Scandinavia.

Other disturbances (wind damage, pest calamities, herbivory)
Other disturbances, however, do play an important role in our study region (Moen et al. 2008).In Lapland, a higher proportion of forest land is classified as damaged to some degree (58.5%) compared to the rest of the country (45.6%), mainly due to the direct and indirect effect of the harsher climate (Ylitalo 2013).Disturbances in Lapland are mainly due to abiotic factors (wind, snow, frost, drought, nutrient imbalance and fire), fungi and moose or reindeer damage (Ylitalo 2013).Importantly, however, the effect on species differs: P. sylvestris is less affected by abiotic disturbances, but more prone to insect damage than P. abies and B. pubescens (Nevalainen et al. 2010).Consequently, it is difficult to define an average return time for patch-destroying disturbances as currently implemented in LPJ-GUESS because not all species (i.e. the entire patch) are affected equally.Additionally, stand-replacing disturbances (e.g.storms or insect outbreaks) are complemented by small-scale disturbances on the level of individual trees in late-successional boreal forests (Aakala et al. 2009;Kuuluvainen et al. 2017): continuous tree falls of old individuals more susceptible to physical damage and decay maintain a mosaic of gaps even in the absence of larger disturbances (Caron et al. 2009).Here, the generalized process representation that encompasses a wide variety of potential disturbances fitting for different ecosystems in global applications is not detailed enough for our regional application.Required improvements are thus species-specific effects of disturbances, contagious spread among neighbouring grid cells where appropriate (e.g.insect outbreaks) and accounting for small-scale disturbances not affecting the entire patch.
Attempts to implement more detailed process representations of specific disturbances in LPJ-GUESS include storm damage (Lagergren et al. 2012) and pest outbreaks (Jönsson et al. 2012).In Lagergren et al. (2012), species-specific storm sensitivity ranked the three tree species in line with results from tree puling experiments (Peltola et al. 2000), i.e.P. abies (1.0) < P. sylvestris (0.5) < deciduous species (0.1).Wind damage increases the probability of pest calamities by providing brood trees for e.g.Ips typographus L., the spruce bark beetle (Komonen et al. 2011).Jönsson et al. (2012) coupled LPJ-GUESS with an I. typographus population model and thus successfully simulated observed outbreaks patterns across Sweden and the resulting damage of P. abies.A similar approach is needed for B. pubescens and Epirrita autumnata L., the autumnal moth.Mass outbreaks of E. autumnata in Lapland have caused the B. pubescens treeline to retreat (Lehtonen and Heikkinen 1995) and will likely be aggravated by climate change (Jepsen et al. 2008).
In Lapland, the effect of herbivores differs greatly between tree species, seasons and regions.Kozlov (2008) found foliar damage of B. pubescens by background insect herbivory to increase from the north (1-2% at 70°N) to the south (5-7% at 60°) of Fennoscandia.During summer, B. pubescens forests are intensely browsed by reindeer (Stark et al. 2007), while reindeer dig for lichens in winter, mechanically damaging P. sylvestris and P. abies seedlings (Helle and Moilanen 1993).High reindeer densities might even limit P. sylvestris recruitment to the extent of preventing treeline advance (Aakala et al. 2014).On the other hand, reindeer grazing reduces competition for P. sylvestris (which is not normally grazed itself), esp.Cladina lichens which negatively affect P. sylvestris mycorrhiza development (Brown and Mikola 1974).Thus, direct and indirect effects of herbivores differ among tree species, and net effects are far from unanimously discussed (Weisberg and Bugmann 2003) which complicates the inclusion of herbivory in LPJ-GUESS.Nonetheless, Zöckler et al. (2008) did include the effect of reindeer grazing in LPJ-GUESS simulations in a very simplified way which needs to be further refined to offset direct and indirect, positive and negative effects on individual tree species.

Scale issues with process-based vegetation models
As similar "fitted" process-based dynamic vegetation models, LPJ-GUESS has been parameterized at certain scales (globally by e.g.Smith et al. (2014), for Europe by Hickler et al. (2012)).Generally, an application on a smaller scale requires accounting for study region-specific processes and the ecology of the main tree species (Hickler et al. (2004) and Tang et al. (2012) for northeastern U.S.; Hickler et al. (2012) for Europe; Seiler et al. (2014) for Bolivia).Zhang et al. (2013) applied LPJ-GUESS to the whole Arctic at an accordingly coarse resolution and reported a good match between observed and predicted treelines, albeit of plant functional types rather than species.Furthermore, they assessed potential natural vegetation -a common LPJ-GUESS application (Wolf, Callaghan et al. 2008;Hickler et al. 2012;Zhang et al. 2013;Zhang et al. 2014) but difficult to validate with observations and recently critically discussed (Chiarucci et al. 2010;Loidi and Fernández-González 2012).In our study region, even the arctic version of LPJ-GUESS (Zhang et al. 2013) incorrectly predicted the whole of Finnish Lapland to be a boreal evergreen forest (while the northernmost part is only occupied by B. pubescens, Fig. 3d).The flexible model design of LPJ-GUESS makes regional adjustments possible, but the parameterization is in many cases challenging.One European parameterization, which reproduces European-wide potential natural vegetation types (Hickler et al. 2012), is clearly not applicable to the study area here, and we expect that the same is true for other smaller-scale regional applications.With our discussion of crucial processes in boreal forests we provide a starting point for future model development.
Total of 122 references.

Fig. 2 .
Fig. 2. Comparison of a) observed (MS-NFI data; red) and simulated biomass [kg m -2 ] (LPJ-GUESS, multi-species run; blue), and, b) multi-species (i.e. with competition) and single-species (without competition; sandy) LPJ-GUESS runs, by latitude bands (lines are means within 0.5° latitude bands).Symbols are transparent to visualize the distribution of values.Note the different range of biomass values for Betula pubescens.

Fig. 3 .
Fig. 3. Map comparison of (a) total and (b-d) species-specific biomass [kg m -2 ]: observed data (MS-NFI, 2011) and results from multi-species and single-species LPJ-GUESS simulations(averaged over 1994-2003).To maximize visibility of spatial differences but retain comparability between observations and simulations, we used quantiles to define classes for each species and the total.This results in the irregular class spacing and reflects the different biomass distributions (cf.Fig.2). "Natural" colours (white to black) cover the range of the observed values (i.e. the upper limit of the black class is always the maximum of the respective MS-NFI data); "artificial" colours (shades of magenta) cover the predictions that exceed the observed value range (model overestimation).

Fig. 5 .
Fig.5.Species-specific response curves of statistical occurrence models (cf.Schibalski et al. 2014, Fig. 7a) to growing degree days (averaged over the ten years preceding the inventory year, i.e.1968-1977 and 1993-2002;  given as rug plot in the upper panels: top 1978, bottom 2003).Vertical lines mark the bioclimatic limit used in LPJ-GUESS for the respective species (GDD 5,min , Suppl.file S1/TableS1).Transparent bootstrapped confidence bands (0.95) were calculated following the procedure detailed inCoutts (2011) andCoutts and Yokomizo (2014), using the boot.cifunction in R(Canty and Ripley 2013).Note: The prevalence of Picea abies was very low in the 1978 dataset (7%) leading to the excessive bootstrapped confidence band (right).
article id 6977 • Schibalski et al. • Identifying important topics for model refinement…