GAPPARD : a computationally efficient method of approximating gap-scale disturbance in vegetation models

Models of vegetation dynamics that are designed for application at spatial scales larger than individual forest gaps suffer from several limitations. Typically, either a population average approximation is used that results in unrealistic tree allometry and forest stand structure, or models have a high computational demand because they need to simulate both a series of age-based cohorts and a number of replicate patches to account for stochastic gap-scale disturbances. The detail required by the latter method increases the number of calculations by two to three orders of magnitude compared to the less realistic population average approach. In an effort to increase the efficiency of dynamic vegetation models without sacrificing realism, we developed a new method for simulating stand-replacing disturbances that is both accurate and faster than approaches that use replicate patches. The GAPPARD (approximating GAP model results with a Probabilistic Approach to account for stand Replacing Disturbances) method works by postprocessing the output of deterministic, undisturbed simulations of a cohort-based vegetation model by deriving the distribution of patch ages at any point in time on the basis of a disturbance probability. With this distribution, the expected value of any output variable can be calculated from the output values of the deterministic undisturbed run at the time corresponding to the patch age. To account for temporal changes in model forcing (e.g., as a result of climate change), GAPPARD performs a series of deterministic simulations and interpolates between the results in the postprocessing step. We integrated the GAPPARD method in the vegetation model LPJ-GUESS, and evaluated it in a series of simulations along an altitudinal transect of an inner-Alpine valley. We obtained results very similar to the output of the original LPJ-GUESS model that uses 100 replicate patches, but simulation time was reduced by approximately the factor 10. Our new method is therefore highly suited for rapidly approximating LPJ-GUESS results, and provides the opportunity for future studies over large spatial domains, allows easier parameterization of tree species, faster identification of areas of interesting simulation results, and comparisons with large-scale datasets and results of other forest models.


Introduction
Forests are an important part of the Earth system, at present covering roughly 30 % of Earth's land surface, and are responsible for about half of the total terrestrial carbon (Fischlin and Midgley, 2007).Ongoing pressures on forest ecosystems including climate and land use change affect forest structure, composition and carbon storage, and changes in forests may in turn feed back to affect climate and ecosystem services (Fischlin and Midgley, 2007;Purves and Pacala, 2008).In order to assess the importance of forests in the Earth system and understand their sensitivity to ongoing environmental change, it is essential to have forest models that can be applied at continental to global scales.Large-scale forest models that can be used to address these questions are complex, as they should include a dynamic representation of forest demography, particularly with respect to forest disturbances and structure-related competition (Quillet et al., 2010;Bonan, 2008).The most widely used tools for assessing the role of plant cover in the Earth system are dynamic global vegetation models (DGVMs).All DGVMs simulate forest growth and include a representation of plant physiology and vegetation dynamics (Prentice et al., 2007), but the first generation DGVMs did not explicitly simulate forest structure, and showed important limitations in their ability to model competition and disturbances (Quillet et al., 2010).Recently, DGVMs have been developed that explicitly account for forest structural characteristics, improve the modeling of competition and small-scale disturbances, and, thus, lead to more realistic simulations of forest dynamics (e.g., Hickler et al., 2008;Sato et al., 2007;Fisher et al., 2010).Therefore, these new hybrids between original DGVMs and models that simulate forest structures, also called second generation DGVMs, have substantial advantages over the original DGVMs in terms of realism, but this typically comes at the cost of computational demand, which puts limits on the spatial domain or maximum resolution that can be simulated in a reasonable amount of time.
One commonly applied computationally time-consuming way of including dynamic forest structure into a DGVM is to apply the "gap model" approach (Shugart, 1984), in which forest dynamics are simulated on small patches that roughly represent the area of influence of one mature tree.Models relying on the gap model approach simulate the fate of individual trees, determined by growth and death processes and a stochastic establishment, leading to demographic stochasticity.Other stochastic elements can be climatic drivers and in particular stochastically appearing small-scale disturbances (disturbance stochasticity).A gap model simulates a number of replicate patches with the same external forcing (climate, soils) and aggregates these when providing grid-celllevel output.Due to the stochasticity, individuals and vegetation biomass on each replicate patch develop differently and simulations of many patches have to be averaged to yield the forest dynamics, requiring a lot of computational time.To obtain realistic results, Bugmann et al. (1996) recommended the use of 200 successive repetitions of simulations per stand for gap models.
An example of one second generation DGVM that applies a gap model approach is LPJ-GUESS (Smith et al., 2001;Hickler et al., 2004), which combines the plant physiological representations of the first generation LPJ-DGVM (Sitch et al., 2003) with the GUESS model of forest demographics (Smith et al., 2001).In LPJ-GUESS, most commonly 50 or 100 replicate patches are used (Koca et al., 2006;Miller et al., 2008;Hickler et al., 2008Hickler et al., , 2009;;Wramneby et al., 2008), but to save computational time, the number of patches is often smaller.Small-scale disturbances have a stronger effect on species composition, forest height, age structure and biomass than the demographic stochasticity (stochastic establishment and mortality).Demographic stochasticity varies only the numbers of individuals in cohorts (by drawing from a Poisson distribution in the establishment function and by imposing expected mortality rates as probabilities for stochastic death in the mortality function) and leads only to moderate deviations from the non-stochastic case.Small-scale disturbances, in contrast, have a strong effect on the simulated forest dynamics (Hickler et al., 2004;Gritti et al., 2006), because they are assumed to destroy all trees in a patch (Fig. 1; b 1 -b n ): after a disturbance all living biomass in that patch is removed to the litter (dead organic matter), and growth succession starts again from the bare ground.As a result, the total biomass (mean of all replicate patches) of a disturbance simulation is typically smaller than in an undisturbed run (Fig. 1a).Second generation DGVMs that use the gap model approach, and specifically LPJ-GUESS, are computationally expensive for two main reasons.They need to (1) simulate plant physiology in 5-50 age-based cohorts that represent the height structure of individual patches, and (2) simulate a high number of replicate patches to represent adequately the variability resulting from stochastic disturbance, establishment and mortality at the grid-cell level.Combined, these requirements increase the computational demand of a second generation DGVM by two to three orders of magnitude as compared to a first generation DGVM.Considering the great impact of small-scale disturbances in LPJ-GUESS and the large number of replicate patches required, we identified a need to develop a new approach to simulating small-scale disturbances that was both efficient and accurate.
LPJ-GUESS has been successfully applied to a wide number of studies over recent years, but because it is computationally expensive, global runs at high spatial resolution (e.g., on the 0.5 • grid commonly used by the LPJ-DGVM, Friedlingstein et al., 2006) are currently impractical without supercomputers.The LPJ-GUESS study that used the largest scale performed simulations for Europe with a 10 resolution, but modeled only 20 stochastic replicates for each grid cell to save simulation time (Hickler et al., 2012).LPJ-GUESS, like generally all DGVMs, has also limitations on a finer scale.Simulations on smaller areas with smaller cell sizes lead to the same numbers of cells simulated.Thus, computational resources are still stretched.Furthermore, applications on finer scales are limited because the parameterization of vegetation in LPJ-GUESS is not specifically adapted to local characteristics.In the first LPJ-GUESS version, tree species were classified into plant functional types (PFTs), since for research questions that address the global scale it is hardly possible to parameterize all species (Jarvis, 1995).Only recently first parameterizations of most common European tree species have been developed (Koca et al., 2006;Wolf et al., 2008a;Hickler et al., 2012) with the aim to reflect the species compositions within the relatively coarse grid cells (10-100 km resolution), which are characterized by mean climate.Additionally, this approach only focused on wide-spread species.These coarse parameterizations thus neglect the heterogeneity of both drivers and comparison data, and are unlikely to work also at higher resolutions (Hickler et al., 2012).Modeling of vegetation dynamics on finer grids results in a higher need for a more specific vegetation parameterization.LPJ-GUESS still does not incorporate the specific vegetation of every important region.For example, there are some area-wide LPJ-GUESS studies for northern Europe (Wolf et al., 2008b;Koca et al., 2006;Smith et al., 2008) in which vegetation was described as plant functional types, but forests were modeled realistically because there are only a few main species present in that region.The latest Europe-wide species parameterization does not include parameters for every and especially locally important tree species, and variations within species were not considered (Hickler et al., 2012).Another important gap of LPJ-GUESS and second generation DGVMs in general is that there is no parameterization specifically for the Alpine vegetation.For example, the two particularly for the Alpine region important tree species Larix decidua and Pinus cembra were still not parameterized for LPJ-GUESS.
In this paper we describe a new approach to simulating forest dynamics that dispenses with the need to simulate replicate patches in a second generation DGVM and is both as accurate as and substantially more computationally efficient than traditional models.We compare our new approach to standard LPJ-GUESS in a series of experiments along an environmental gradient in the Swiss Alps and demonstrate the quality of the new method.We adapt the tree species set of LPJ-GUESS to an Alpine environment and use Swiss National Forest Inventory data to adjust the vegetation parameters.We then suggest potential applications for this new, efficient model for addressing large-scale problems on the role of forests in the Earth system, and, with regard to the new parameterization, specifically in the Alps.

Material and methods
2.1 New approach to include small-scale disturbances

Basic assumptions
We maintain the idea that a forest consists of many patches, each of which is affected by disturbance independently, and that disturbances work on a yearly time step.We also maintain patch-destroying disturbances (i.e., living biomass state variables are set to 0 by a disturbance).
Consequently, at a given time T , the patches have different patch ages a that depend on the times T − a, when they were affected by a disturbance, and the distribution of patch ages P (a) at time T determines the forest state.A patch of age a has state variables (e.g., numbers per species and height class, sapwood, and heartwood mass, etc.) and output variables (e.g., biomass per height class).Here, we refer to both as output variables y(a) because each state variable can easily be treated as an output variable and because we apply the new method in a postprocessing way.
Our approach is based on the idea that a forest does not necessarily have to be represented by different replicate patches but can be calculated using a small number of undisturbed simulations starting from different time points, using the information of the patch age distribution.This includes a temporal upscaling of the information gained from such undisturbed, deterministic, and thus computationally efficient model runs.These runs include exactly the same model functions as the base models but do not simulate patchdestroying disturbances (if included in the base model), and switch off stochasticity.This, for example, means for LPJ-GUESS that the mortality and establishment functions are calculated in the deterministic mode: mean values to determine the number of dying and newly establishing individuals are used.We will refer to our new method as the GAPPARD method (approximating GAP model results with a Probabilistic Approach to account for stand Replacing Disturbances).

Constant drivers
If the drivers are constant (e.g., during the spin-up phase of a climate change simulation), we assume that each re-growth following a disturbance leads to the same development of output variables as the initial forest development (Fig. 1b), i.e., y(a) = y(x), with x the time since the start x 0 of the undisturbed simulation and a the time since the last disturbance in the disturbed simulation (i.e., patch age).In our simulations, we do not dynamically update several state variables that in reality would be indirectly affected by disturbance, including soil moisture and temperature, and the state of the snowpack.
To have age a at a given time T , a patch of a disturbed forest first must have encountered a disturbance and subsequently survived a years.Given the probability of a disturbance p dist , the probability that a patch afterwards survives a years without any disturbance is (1 − p dist ) a .Consequently, the probability that a forest patch has age a is (1) A special case is given by the patches surviving from the beginning (T = 0) to exactly time T ; they have not encountered any disturbance, but started from bare ground ("were killed") for sure (p dist = 1) at time 0: ( The expectation value Y of simulation result y is then given by To calculate the resulting expectation values, we first perform one simulation without disturbances (SWD) leading to y(x) for all time points x.Afterwards, Eq. ( 3) must be applied in a postprocessing step.
The method presented here is a modified version of the von Foerster equation (von Foerster, 1959), a general agestructured population dynamics approach, in which instead of the patch age distribution the individual age distribution is constantly changed during the simulation.

Changing drivers
When drivers change, disturbances occurring at different times have different impacts.For example, under low temperature conditions, succession after a disturbance will most probably be slower than in a warmer climate.In order to account for such transient drivers, we modify the standard method of running only one SWD (described above) by running several SWDs with different starting times s i ("nodes"), trajectories starting from si and s(i+1); dotted line: result after applying the GAPPARD method using both trajectories.b1), b2) weight of the trajectory on the calculation of a state variable at time point T with age a for all T − a that are between si and s(i+1), where T is a year after si to which the output variable is calculated and a is the patch age (time since last disturbance).A darker trajectory stands for a bigger influence.The solid lines characterize the weight of the trajectories for T − a.Here, the trajectory starting from s(i+1) has a higher weight because T − a is closer to s(i+1) than to si.

Fig. 2. Principle of GAPPARD applying changing drivers.
(a) dashed lines: developments of two subsequent trajectories starting from s i and s (i+1) ; dotted line: result after applying the GAP-PARD method using both trajectories.(b 1 ) and (b 2 ) weight of the trajectory on the calculation of a state variable at time point T with age a for all T − a that are between s i and s (i+1) , where T is a year after s i to which the output variable is calculated and a is the patch age (time since last disturbance).A darker trajectory stands for a bigger influence.The solid lines characterize the weight of the trajectories for T − a.Here, the trajectory starting from s (i+1) has a higher weight because T − a is closer to s (i+1) than to s i .i = 0 . . .n.This yields different trajectories of the output variables y i (s i + x), each starting at a different starting time s i with y i (s i ) = 0.
For each time point T of the output, E[y(T )] is determined similarly to Eq. (3).However, in this case instead of using one undisturbed y at one time point x = a in the summation, the two values that belong to two subsequent trajectories starting before and after the target time point T − a are used to describe the state of a patch with age a (Fig. 2).
The two output values y(s i + a) and y(s (i+1) + a) are then interpolated according to the distance of T − a to the nodes s i and s (i+1) , so that the trajectory with the node closer to the target time T − a has more weight than the other one.

Model application and evaluation
To evaluate the GAPPARD method, we applied it to LPJ-GUESS (LPJ-GUESS-G) and compared the results to the stochastic runs of LPJ-GUESS.We used only LPJ-GUESS to parameterize our tree species: if the parameters are valid for LPJ-GUESS, they must be also valid for LPJ-GUESS-G, because it is an upscaled version of LPJ-GUESS.

Location and climate data
We selected the Rhone Valley in the Swiss canton of Valais to test our new approach.The bottom of this valley is one of the driest regions in Switzerland, and the hillsides lead to steep gradients in environmental conditions.In the region, management generally did not affect the species compositions to such extents as in most Swiss regions, so that we were able to use recent forest data to parameterize our modeled tree species.We selected eight stands along a north-facing transect (Fig. 3, Table 2) that cover the vegetation zones where homogeneous forest areas exist (from ca. 150 m above the valley bottom to tree line).Each stand location was derived from climate data points of 100-meter grid in a way that the altitudinal distances between two stands would be approximately 200 m.Climate data are applied after the simulation year 1900.Up to 1900 we applied randomly selected values of the first 30 climate data years for the model spinup.For the 1901-1929 simulation period, we used CRU (Climatic Research Unit's monthly TS 1.2 dataset) data downscaled to a 100-meter grid (Mitchell et al., 2004).For the 1930-2006 simulation period, we used Swiss weather station data from the Federal Office of Meteorology and Climatology MeteoSwiss interpolated to a 100-meter grid by applying the Daymet method (Thornton et al., 1997).For the 2007-2100 simulation period, we used CRU climate data of the A1b climate scenario (Mitchell et al., 2004).Along with that scenario CO 2 levels reach 703 ppm in 2100 (IPCC, 2001, Annex II).
Based on the soil suitability map of Switzerland (Frei, 1976), we chose to use a low value of usable volumetric soil water holding capacity of 0.1 (fraction of soil layer depth) and a value for soil thermal diffusivity at 15 % water holding capacity of 0.8 mm 2 s −1 .These values correspond to the poorly developed soils on the slopes of the Rhone Valley.

Tree species parameterization
Using LPJ-GUESS, we optimized the parameterization of each tree species present in our study area to obtain the best possible fit to observed forest inventory data.We used the stochastic LPJ-GUESS for the parameter optimization because this model served as the reference for the subsequent model comparison, and applied the same optimized parameter set for simulations with the GAPPARD method.The tree species parameters we used are generally based on the existing LPJ-GUESS parameterizations for plant functional types (Hickler et al., 2004;Wolf et al., 2008a) and for species (Koca et al., 2006;Miller et al., 2008;Hickler et al., 2012).
In our experiments, we used the most abundant Swiss forest species, selected according to the species used in Lischke et al. (2006a), and then analyzed those of them that already had been parameterized for LPJ-GUESS.We excluded all LPJ-GUESS species not present in the Swiss Alps, and added parameterizations for three new species that are abundant in our study area: Larix decidua, Pinus cembra and Pinus mugo.For Larix decidua we generated an additional function to model its leaf phenology based on results presented by Migliavacca et al. (2008) (Appendix A2).
We further optimized the species-specific parameters used in our experiments so that model results would best match forest inventory data from the Swiss National Forest Inventories NFI1 (EAFV, 1988) andNFI3 (Brändli, 2009).We selected inventory data only from plots located south of the Rhone and within a 30 km distance of our simulation plots, and further stratified the inventory information into eight altitudinal classes analogous to the altitudes of the eight simulated stands.At each altitudinal class we calculated the mean and the standard deviation of the biomass of all living tree species and estimated carbon mass assuming that half of a tree's biomass is carbon (IPCC, 2003; see Appendix A1 for more details).
We used total forest carbon mass (above-and belowground) as the sole main metric for evaluating the model performance in light of the NFI data.Our aim was to optimize the model parameters so that for each altitudinal class the simulated total and main species carbon mass were similar to the NFI data.The parameters we optimized for all species included drought tolerance, minimum growing degree day sum, and maximum temperature for establishment (d tol, gdd5min e and tcmax e in further adjusted the allometric parameters for some species, in particular the steepness-influencing parameter in diameter to height relation (k allom2 in Table D3 of the Appendix).Further details on the parameterization, especially for the tree species newly added to LPJ-GUESS, are described in Appendix A2 and Tables D1-D3.
We applied one altitude-specific set of LPJ-GUESS average return intervals for generic, patch-destroying disturbances (RID, inverse of p dist in Sect.2.3.2).We assumed that the stands along the height gradient underlie different disturbances.In general, stands close to the valley bottom are more frequently disturbed by fire in the Valais region (Zumbrunnen et al., 2009).Furthermore, we assumed that more uphill stands are disturbed more frequently by storm events, rockfall and especially avalanches.On the one hand altitude is not a good explanatory characteristic for avalanche appearance, and Schneebeli and Meyer-Grass (1992) found that in spruce-and larch-dominated forests steepness favors the release of avalanches, but on the other hand a thin crown cover, big gap lengths and higher proportions of larches also increase the probability of avalanche release.The steepness is approximately the same for all our simulated stands except the most upper stand, which is on slightly flatter ground.At this highest elevation site however, shade-intolerant larch is very abundant; the trees are exposed to avalanches from higher altitudes, while the annual maximum snowpack is deeper, and trees grow with larger gaps surrounding them.First tests with the same RID values used for all stands led to less precise results, especially concerning the total modeled carbon mass (not shown).To create an altitude-specific disturbance distribution, we used an RID value of either 65 or 100 for each single stand, depending on which value suits best to predict the total carbon mass (RID in Table 2).
To evaluate the tree species parameterization, we used two indices: the Euclidean distance scaled to one and the percentage similarity coefficient (Bugmann, 1994).Details on both indices can be found in the Appendix.

Analyzed output variables
For our analysis, we used as investigated output variables (1) the total of the tree carbon mass of a species, and (2) this variable stratified by certain height classes (of 4 m height except the lowest being 2 m high).We examined the course of these variables from an LPJ-GUESS simulation including small-scale disturbances and 400 replicate patches, and compared the outcome to the results of LPJ-GUESS-G.

GAPPARD method versus stochastic LPJ-GUESS model runs
The crucial test of the GAPPARD method is its ability to reproduce the behavior of LPJ-GUESS disturbance runs in terms of total carbon mass, species composition and height structure.To apply the GAPPARD method, we used the RID values of the defined altitude-specific disturbance distribution (RID in Table 2) and took its inverse values for Eqs. ( 3) and ( 5).Using the new set of parameters and the altitudespecific disturbances, we simulated forest growth for all eight stands with a spinup time of 800 yr and a total simulation time of 1000 yr covering a simulation period from 1100 to 2100.
In addition to a full time SWD (including 800 yr spinup), we used four nodes from which we started SWDs that provide the input for the GAPPARD method to account for climate change : 1950, 2000, 2050 and 2080.All simulations ended with the year 2100.
Without disturbance, as forests become taller the model needs to calculate light interception in an increasing number of foliage layers (Prentice et al., 1993), which consumes computational resources.Therefore, we increased the depth of these layers from 2 to 5 m after 200 yr and to 10 m after 400 yr for the full time SWD.This change did not lead to a decline in result quality because without disturbances forests become homogenous in such a way that a less detailed light calculation does not have much influence.In addition, such old forests are rather rare so that such a simplification has an even smaller impact.
To examine the analysis, we first tracked the total carbon mass of the different species from the end of the spinup phase until 2080.Second, we mapped the carbon mass results of the different species along height classes for all stands and two time points of simulation: 1900 and 2080.We did not analyze simulation years after 2080 because they do not represent interpolated results between two nodes but a prolonging of 2080 conditions (Eq. 5, case a > y(s n )).Third, to quantify the quality of the results we calculated the root mean square error (RMSE) for each stand and each species, based on simulation results of a 10 yr resolution.The RMSE corresponds to the differences in carbon mass between two models (described in detail in Appendix C), and is calculated from the sums of all individuals of one species reaching a certain height class and an adjacent height class (tolerance to difference in height classes).For every species, each of these differences enters into the calculation of the RMSE as a fraction of the maximum possible difference appearing in that stand and the calculated simulation period.Hence, the maximum RMSE is one (completely different results).We calculated the RMSEs separately for two time periods: the spinup period and the climate change period.For both, we calculated the RMSE between LPJ-GUESS with 400 replicate patches and (1) LPJ-GUESS with 100, (2) LPJ-GUESS with 25 replicate patches and (3) LPJ-GUESS-G.
To test the model performance, we compared the simulation times of stochastic LPJ-GUESS model runs with 400, 100 and 25 replicate patches, and LPJ-GUESS-G.Listed times needed for the simulations using the GAPPARD method consider only initial runs.The computational time needed for the GAPPARD method is negligibly short.The simulations ran on one core of an AMD Opteron 2439 2.8 GHz processor.

Tree species parameterization
We successfully parameterized the LPJ-GUESS forest model using data from the Swiss National Forest Inventories (NFIs) of the Rhone Valley.We achieved a high accordance between the used NFI data and the simulated data of carbon mass for total living forest and for main species (Table 3).
Combined, the NFI1 and NFI3 data show an increase in forest biomass from the lowest site at 800 m to middle elevations of approximately 150 t ha −1 (from 100-120 to 250-270 t ha −1 ; from 5-6 to 12.5-13.5kgC m −2 ) followed further upslope by a decline in biomass to 130-170 t ha −1 (6.5-8.5 kgC m −2 ) at the upper alpine vegetation zone (see Tables D4 and D5 in the Appendix for more details).In general, we were able to simulate these altitudinal shifts, although we rather overestimated the NFI1 data and rather underestimated the NFI3 data.Consistent with the data, our simulations also show a small increase in forest carbon at all sites between the years in which the NFI was performed.We have reached percentage similarity coefficients higher than 0.93 and Euclidean distances scaled to 1 of smaller than 0.35 for the comparison of the total carbon mass over all altitudinal classes (T11 in Table 3).
With the new LPJ-GUESS parameters and functions, we were able to simulate the general pattern of a dominant Pinus sylvestris at stands closer to the valley bottom, dominant Picea abies and Larix decidua in forests at a mid-altitudinal elevation, dominant Pinus cembra and Larix decidua in forests at the upper alpine vegetation zone, and a continuous decrease of the proportion of broad-leaved tree species with altitudinal height (Figs.D1 and D2 in the Appendix).The newly parameterized LPJ-GUESS simulated the NFI distribution of most species with percentage similarity coefficient values of over 0.5 and Euclidean distances scaled to 1 of smaller than 1.0 (Table 3).There are only two major exceptions.First, the NFI data, especially at mid-altitudinal elevations, show high biomass of Larix decidua that LPJ-GUESS does not capture with the parameterizations, model functions, and initial modeling conditions we used.Concomitantly, the model generally predicts more Picea abies than is observed in the NFI data, so the sum of the carbon mass of   both species leads to a good coherence to the NFI data.Second, the modeled carbon mass of Quercus pubescens exceeds the NFI1 values in stand a, which is reflected in both used indices.Moreover, although there is basically no Tilia cordata and Carpinus betulus appearing in the NFI data, both species established as minor tree species in lower altitude stands in LPJ-GUESS.

Development of forest species composition and carbon mass with LPJ-GUESS
The climate change scenario we applied led to several changes in all simulated stands.Over the period from 1900 to 2000, the total carbon mass in most stands slightly increased.During the 21st century the carbon mass increased in all stands with highest increase at the beginning of the 21st century (black solid line in Fig. 4).At lower stands Quercus pubescens and Picea abies (stands a and b) or drought-tolerant broad-leaved species (stands b to e) profited from a decrease of Pinus sylvestris.In all stands, the proportion of broad-leaved species increased.Generally the increase of broad-leaved species was lowest at higher altitudes (< 0.1 kgC m −2 in stand h), and highest at low altitudes (approximately 1 kgC m −2 in stand a).Like Quercus pubescens and Pinus sylvestris, with climate change Picea abies established at higher altitudes.Although no Picea abies appeared in stands g and h in 1900, in 2080 this species made approximately one-third of the total carbon mass in stand g and one-fifth in stand h.
The three stochastic LPJ-GUESS model runs (using 400, 100 and 25 of replicate patches) in the long term showed similar results for tree carbon mass development (Fig. 4).However, it is clearly illustrated, that carbon mass can vary strongly for decades as well as for centuries (in the case of 25 replicate patches).A high number of replicate patches minimizes the chances of intensively altering output variables.However, as can be clearly seen in Fig. 4, even the results of the 400-replicate and 100-replicate simulations are quite different.Only for periods with extreme climatic situations, all of the three LPJ-GUESS runs were affected equally (e.g., in the beginning of the 1920).

Development of forest species composition and carbon mass with the GAPPARD method
Our comparison between LPJ-GUESS results using 400 replicate patches and LPJ-GUESS with GAPPARD (LPJ-GUESS-G) shows that the GAPPARD method successfully  D6 in the Appendix).The total carbon mass produced with LPJ-GUESS-G was in the range of LPJ-GUESS results for each stand and at any simulation time (Fig. 4).However, using the GAPPARD method smoothed the results over time so that changes of output variables only occur gradually with the simulation years, which can be seen in total carbon mass results and for single species.Also concerning species composition, LPJ-GUESS-G was successful in reproducing LPJ-GUESS results.Mean root mean square error (RMSE) values are in the range of the mean RMSE between LPJ-GUESS using 100 and LPJ-GUESS using 400 replicate patches.The only case where LPJ-GUESS-G produced a different species composition was in stand g.There, Picea abies established in the middle of the 20th century in the LPJ-GUESS stochastic simulation run, but arrived only about a century later in the LPJ-GUESS-G run.
Our comparison of simulation times shows that the GAPPARD approach performed approximately 11 times faster than LPJ-GUESS using 100 replicate patches (Table 4).The full time simulation without disturbance contributed to a great extent to the simulation time of the GAPPARD approach (results not shown) because there trees can get much older and taller.These taller trees claim unproportionally more simulation time than smaller trees mainly because of a more complex light calculation.
For the spinup period, the mean RMSE between LPJ-GUESS using 400 replicate patches and LPJ-GUESS-G was approximately 0.1 (calculated for species that produced at least 0.5 kgC m −2 in one of both models).For the simulation period, the mean RMSE between LPJ-GUESS and LPJ-GUESS-G was approximately 0.05.Reducing the numbers of replicate patches in a stochastic LPJ-GUESS run down to 25 resulted in a mean RMSE (between the 400-and the 25-replicate simulation runs) of approximately 0.09 for the spinup period and 0.07 for the simulation period, and was roughly four times faster than the 100-replicate simulation (Table 4).

Discussion
This study identified a novel, efficient method to run disturbance-driven models and also yielded an LPJ-GUESS parameterization and adaptation, and climate change simulations for a region with specific properties.
Table 4. Simulation times and mean RMSE.Sum of simulation times of all 8 simulated stands for 800 yr of spinup and simulation from 1900-2100 with LPJ-GUESS using replicate patches and using the new GAPPARD method.RMSE for model spinup of 800 yr (RMSE SPI) and for the simulation time 1900-2080 (RMSE SIM).The mean RMSE is calculated using all species that produced at least 0.5 kgC m −2 in one of both models at any point in time during the simulation.The reference to RMSE calculation is always LPJ-GUESS with 400 replicate patches.The 400-replicate run was done once, the others shown as means of 10 simulation runs that use the same random number of seeds for stochasticity (thereby producing the same output) to account for the different loads of the processor nodes.Simulations were run on one core of an AMD

Parameterization
For the first time the Alpine mountain forest species Larix decidua, Pinus cembra and Pinus mugo were parameterized for LPJ-GUESS.We were able to include both new Pinus species by using only existing model functions.One achievement of the parameterization is the newly added function to model the leaf phenology of Larix decidua (see Sect. 2.2.2 and Appendix A2).We used observed forest biomass from inventory data as the only variable upon which to optimize the LPJ-GUESS species-level parameters.Although this simple approach does not account for the properties influencing forest dynamics, such as tree age, height or width, it led to a considerable similarity between the newly parameterized LPJ-GUESS results and the NFI data.With the NFI data, we were able to use a relatively dense sample of plots close to the modeled stands, which makes the results highly reliable.Sources of error mainly concern the estimation of the biomass from the NFI data on the one hand, and the comparison with LPJ-GUESS allometry-based simulated values on the other hand.Using the LPJ-GUESS allometry function, trees with identical diameter always have the same height and therefore the same mass.As a consequence, it could be possible that the model results might match the carbon mass results of the NFI, but not the diameter at breast height.This may particularly be the case in high mountain settings where tree allometry is strongly influenced by meteorologically imposed constraints (e.g., wind), which we did not consider in our models directly.Such an inconsistency, which might vary for different species, could lead to an unrealistic simulation of the vertical structure of the forest, and thus to further deviations from observations.However, the good agreement between our simulations and the observed biomass data indicates that differences in allometry are compensated for in the model (e.g., through changes in stand density controlled by disturbance).Another source of error may be linked to the development of the forests in the used NFI plots.It is not clear in detail how strong the influence of the management in the specific NFI plots is.It is also not exactly known on which stands clear cuts or big disturbance events have taken place.We assume that the mean state of all used stands is representative of the overall disturbance regime, so that the constant probability of disturbance events we use can reproduce it.Despite all sources of error, our parameterization results show high consistency between NFI data and simulated values.Furthermore, the simulated carbon mass and forest species distribution in the different stands at the time of the NFI dating are plausible to a great extent.

Simulation results
Using the altitude-specific disturbance distribution and the new parameters, we were able to simulate total forest carbon mass and a species composition that largely reflects the NFI data.We were able to reproduce forest species composition and tree carbon mass without using a specific disturbance function like the LPJ-GUESS fire function (Thonicke et al., 2001).However, to account for feedback effects between fire and forest growth and the spreading of fire, and both especially in the context of climate change, an appropriate modeling of the fire function could be important for area-wide simulations of the modeled region.However, as our main concern was to present a new modeling technique, we accept approximations regarding missing modeling detail of disturbance types.The few discrepancies between the simulation results and the NFI data may be partly due to uncertainties in the interpretation of the data.It is very likely that the situation with the high carbon mass of Larix decidua and low one of Picea abies in mid-altitudinal stands is to a high degree a result of management practices in the Rhone Valley (Gimmi et al., 2010).In addition, including specific total disturbance years (removal of all living carbon mass of all replicate patches in one pre-defined year) did not lead to a significant increase in the biomass of early successional Larix decidua because Picea abies completely overgrew it after only 30 yr of succession (results not shown).Although this suggests that succession may occur too quickly in LPJ-GUESS, it can be assumed that in the absence of forest management the real forests analyzed here at altitudes between 1200 m and 1800 m would be strongly dominated by Picea abies (Frehner et al., 2005).For the species Carpinus betulus and Tilia cordata, the simulation results exceeded the NFI data, which possibly is enhanced by not modeling all existing broad-leaved summergreen species.It is generally known that Carpinus betulus does not exist in this part of the Rhone Valley (Welten and Sutter, 1982), which is also reflected by the NFI data.In contrast, Tilia cordata usually exists in this region (Welten and Sutter, 1982;Svejgaard Jensen, 2003), but for unknown reasons is found only on very few plots of the NFI data in the analyzed region.
The increase of tree carbon mass with climate change, also at the lowest site, can mainly be explained by the 20th century increase in atmospheric CO 2 concentrations (Fischlin and Midgley, 2007).Test simulations with a constant level of CO 2 led to a decrease in tree carbon mass in the 21st century in most plots (results not shown).This decrease can be traced back to a drier climate with more pronounced water stress so that growth rates might be reduced and respiration rates increased.With rising CO 2 levels this effect can be compensated, at least to a certain degree, because stomatal regulation maintains high production rates and prevents plant water loss.In our results, the increase of tree carbon mass is more pronounced in stands with more frequent disturbances (close to valley bottom or upper tree line).Two reasons may be responsible for that.First, species more adapted to the new climatic conditions could profit from opening gaps.Second, species adapted to previous climatic conditions (which survive more likely with less disturbances) could be affected negatively by the changing conditions opening space for new species.However, at the lowest stand (ca.800 m) the increase of carbon mass may be unrealistically large.Closer to the bottom of the Rhone Valley, extreme drought events causing forest diebacks are expected to occur more frequently in the near future (Rebetez and Dobbertin, 2004).But the modeled future carbon stock does not show the effect of such extreme events (compare with the dry year 1921 in Fig. 4) suggesting that the downscaled climate data that we used may underestimate the climate variability, and thus drought events, or that the CO 2 effect is too strong.In addition, the northfacing aspect of all the sites we simulated means that they may be relatively less sensitive to interannual climate variability than the inner-Alpine region as a whole.However, the carbon mass of LPJ-GUESS results does not continuously increase over time, mainly due to stochastic variances, so that long-term trends are much more significant than changes in carbon mass over decades.This is one reason why a quantitative analysis of the carbon mass development between NFI1 and NFI3 is not fully reliable.However, the limited sample size of the NFI data does not allow the strict quantitative analysis of changes in carbon stocks.
With climate change, the occurrence of Pinus sylvestris at stands closer to the Rhone Valley bottom will most probably be reduced and exchanged by broad-leaved droughtresistant species like Quercus pubescens (Rebetez and Dobbertin, 2004;Bigler et al., 2006).This is well reflected in our modeling results.But it is not entirely clear why in our simulations with LPJ-GUESS the carbon mass of Picea abies increases at the lowest stand, while the mass of Pinus sylvestris decreases, especially because the latter is more drought-resistant than Picea abies.We assume that the simulated decrease in Pinus sylvestris biomass is partly because, using LPJ plant physiological functions, it profits far less from increased CO 2 levels in comparison with broad-leaved summergreen tree species (Cheaib et al., 2012).Another reason for the decrease in Pinus sylvestris biomass might be that the lowest stand is still roughly 150 m above the valley bottom, so that conditions are still good enough for Picea abies.This is also reflected in the NFI data where the carbon mass of Picea abies at the lowest stand approximately doubled between NFI1 and NFI3 (Fig. 5, upper row).The general increase from the NFI1 to the NFI3 biomass confirms the increase of forest biomass, although this change also can be due to past forest management and the prevention of disturbance events.
It is unclear whether the modeled shift of species to higher altitudes as a consequence of climate change happened in a reasonable amount of time.In LPJ-GUESS, tree establishment of new species only depends on the environment and does not consider changes in and feedbacks to the seed pool.It is well known that modeling seed pools and the dispersal of seeds have a large potential to change simulation results (Lischke, 2005;Lischke et al., 2006b;Epstein et al., 2007;Neilson et al., 2005).Incorporation of seed dispersal and migration into LPJ-GUESS remains an open problem for future research.

GAPPARD method
With GAPPARD we utilized a modified version of the von Foerster equation.Several other approaches also used von Foerster types but during the simulations for each year and without using computationally efficient interpolation methods (Kohyama, 1993;Moorcroft et al., 2001;Falster et al., 2010).Despite the success of our method, there are some limitations.With the method presented here, it is currently not possible to include any spatial interactions between neighboring grid cells or patch-to-patch interactions.Therefore, seed dispersal or the spatial mass effect of LPJ-GUESS (establishment in a patch depends on carbon mass of other patches in a stand) cannot be applied yet.The stochastic mortality and establishment functions of LPJ-GUESS seem to have a much smaller effect on forest carbon mass and species composition than do stochastic small-scale disturbances.With a more significant impact on demographic stochasticity, LPJ-GUESS might lead to results that could not be reproduced with the GAPPARD method as adequately as we show here.In this case, the methods used in the forest models TreeMig (Lischke et al., 2006b) or TreeM-LPJ (Scherstjanoi et al., 2013) to model vertical and horizontal heterogeneity could provide possible solutions.Although it might not have an influence on the stands modeled here, one additional limitation of the GAPPARD method could be that the influence of climatic extreme events is not visible in the model output because of the linear interpolation between different initial undisturbed runs (see Sect. 2.1.2,Eq. 5) and because of the preset node positions.Thus, a central question for future applications would be the number and setup of starting points of simulations without disturbances (SWDs).Especially if a disturbance has happened before an extreme climatic event (e.g., extreme dry years) the following succession may show large differences to the mean calculated with our new method.Hence, nodes must be set such that longterm trends and short-term variabilities are depicted.Here, we use pre-defined starting points for the SWDs (nodes) chosen independently from climate data.However, despite of using this simplification concerning extreme climatic events we met our goals, also because unresolved spatial heterogeneity (e.g., microhabitats) implies that not all patches respond in the same way to extreme climate variability and extreme events.Still, this simplification might be the reason for the underestimation of Picea abies growth in stand (g).Our first SWD starts in 1950.An additional SWD starting in 1900 did not substantially change the model results (not shown).Hence, the combination of suitable establishment conditions in one or more specific years in the beginning of the 20th century in combination with disturbance events must be the reason why Picea abies establishes earlier in LPJ-GUESS.
Furthermore, it remains unclear whether a different interpolation method (e.g., spline interpolation) between two nodes could lead to even better results.In addition, errors may result from ignoring that some state variables are not set to 0 after a disturbance event in stochastic LPJ-GUESS simulations (e.g., amount of litter, soil water and snow layer) but with our method we assume they were.Another source of error can be explained by the effect of the smaller sample size in our method on the establishment of saplings compared to the gap model approach of LPJ-GUESS.Using replicate patches theoretically more combinations of different light conditions (e.g., very good light conditions after disturbances) and different climatic conditions can be attained.
It is important to note that at the moment the new method cannot serve to track output variables in the same way as the original LPJ-GUESS can.However, the method presented here is not designed to model carbon or nutrient cycles.Furthermore, the method is best applicable if state variables that are typically reset by disturbances do not influence the subsequent tree establishment.

Conclusions and outlook
With GAPPARD, we provided a new method of efficiently simulating the dynamics of tree biomass and forest species composition.It can be used for any output variable that can be produced with the deterministic run and that is reset by disturbance.GAPPARD can further be applied for any model that uses a gap model approach and that applies disturbances as stand-replacing events.Our simulations demonstrated that the GAPPARD method can also be used for simulations that consider the transient effects of changes in climate and atmospheric CO 2 concentrations.Moreover, the principle of the method could be applied to implement newly the effect of stand-replacing disturbances in any dynamic forest or vegetation model.
The GAPPARD method is particularly suitable for simulating a great number of stands in a fast way, and hence is applicable on larger scales.The results can be used to make first estimations about the development of output variables (e.g., species composition) or to identify hot spots of unusual or interesting simulation results, which then can be analyzed in more detail with the original models.
As a next step, we plan to apply the efficient method developed here in combination with the optimized specieslevel parameter set for Swiss tree taxa in Switzerland-and Europe-wide simulations on a 1 km grid.Furthermore, we plan to extend the method by implementing effects of demographic stochasticity, non-stand-replacing disturbances and spatial interactions.

Appendix A Parameterization details A1 Adaptation of NFI data
The selection of the NFI plots is based on the distance to the simulated stands, on the stand type, and topographical considerations.We used only NFI plots that were classified as accessible forest areas, and that are all located south of the Rhone, at most 30 km westward or eastward of the simulated stands, and at most 30 km southward of the southernmost stand.We also classified the chosen NFI plots according to their exposition, but the results were not sensitive to it (results not shown).
The NFI1 and NFI3 data of all plots are split into two parts.One part comprehends trees with diameter at breast height (DBH) higher than 12 cm (older trees), and the other trees with DBH lower than 12 cm (young trees).For the older trees, the biomass per area is estimated for each occurring species.The young trees in the NFI1 are classified into DBH classes of 0-4 cm, 4-8 cm or 8-12 cm, or are classified as 30-130 cm high.The young trees in the NFI3 are classified into height classes of 10-40 cm or 40-130 cm, or are classified as having a DBH of 0-12 cm.To estimate their biomass we used the mean values of the classes, applied it to the LPJ-GUESS allometry function and calculated the biomass considering wood density (Assmann, 1962).
For the parameter tuning, we utilized a set of simulation runs that all used 800 yr of spinup period and simulated forest developments from 1900 up to the years the NFI data were estimated.To smooth stochastic variations widely, we used a number of 400 replicate patches for the final parameter finetuning simulation runs.

A2 New plant physiological functions and parameters
We added three new tree species to LPJ-GUESS that have not been included before: Larix decidua, Pinus cembra and Pinus mugo.Hence, we had to parameterize them from scratch.Both Pinus species were applicable to existing functions of LPJ-GUESS.But first plausibility tests showed that these functions were not sufficient for Larix decidua, mainly due to the tree species' specific phenology.In LPJ-GUESS, the foliage of summergreen species is transferred to litter all at once on one simulation day (typically in autumn) when the maximum number of equivalent days with full leaf cover per growing season exceeds a certain value.For most species, this approximation has no significant negative influence because photosynthetic efficiency in general is reduced more suddenly.But especially for larches, leaf senescence can be a process lasting for months during which photosynthetic intensity is reduced stepwise.Based on Migliavacca et al. (2008), we included this physiological trait by defining a new phenology type for Larix decidua.It will be modeled like a summergreen species, but in autumn the phenological state of the larches will decrease with an S-shaped curve, depending on the number of days since the start of fall of leaves (sd) and the number of days with full leaf cover this year (md) due to the LPJ-GUESS leaf phenology function: We determined the other parameters of Larix decidua oriented on expert knowledge and literature about the species (Table D3).We defined it as a shade-intolerant species with a high ratio of leaf area to sapwood cross-sectional area (Oren et al., 1995).Although it is a boreal species that also grows under very cold conditions, it can establish under warmer conditions, too.Furthermore, saplings do not need much soil water for establishment.The parameters of the new Pinus species are mainly based on Pinus sylvestris parameters.However, both new Pinus species are more cold-resistant, have seeds that are less drought-resistant and their needles have a higher longevity.Moreover, Pinus mugo was defined as shade-intolerant.
An important issue was the parameterization of the drought tolerance (d tol in Table D3) of the three new species.Referring to Ellenberg (1986), Bugmann (1996) defined Larix decidua as a rather drought-intolerant species.However, he also listed other authors that had defined intermediate values for it.Eilmann and Rigling (2012) showed that Larix decidua is, in comparison to, for example, Pinus sylvestris, strongly affected by drought events.On the other hand, Lischke et al. (2006b) and Shuman et al. (2011) defined it as a very drought-tolerant species.Yan and Shugart (2005) used different larch species but also defined them as droughttolerant.Matras and Pâques (2008) noted that the response of Larix decidua to drought can vary strongly depending on stand conditions.Klimek et al. (2011) discussed the water consumption of Larix decidua seedlings compared to other conifers, and reported different observations ranging from same water uptake rates among all conifers to studies that show that Larix decidua consumes 10 times more than other conifers.If soil conditions allow it, this species can survive dry years better than other species that are more droughtadapted because its root system is very deep (Anfodillo et al., 1998;Valentini et al., 1994).Considering all this, we determined Larix decidua to have a high proportion of fine roots in the deeper soil layer and a moderately high drought tolerance (Table D3).Based on Bugmann (1996) and Lischke et al. (2006b), we also used a rather high drought tolerance for Pinus cembra and Pinus mugo.According to Valentini et al. (1994), Pinus cembra with its root system is also able to use groundwater from deeper layers, which is also in line with general knowledge.Thus, we determined a high proportion of fine roots in the deeper soil layer for both new Pinus species.
We modified several values of the latest existing LPJ-GUESS species parameterization (Hickler et al., 2012) (Table D3).Main changes address the soil water content needed for establishment, which we increased for Fagus sylvatica, Abies alba and Quercus robur, and decreased for Betula pendula and Picea abies.Furthermore, we introduced a new shade tolerance class, particularly for Picea abies (column "ns" in Table D1).With this, we contribute to Picea abies being less shade-tolerant than Fagus sylvatica or Abies alba but more shade-tolerant than intermediate shade-tolerant species (Bugmann, 1994;Roloff, 2010).
Based on the Switzerland-wide applied forest models FORCLIM (Bugmann et al., 1997) and TreeMig (Lischke et al., 2006b), we additionally adapted the allometry parameters.As a result, Betula Pendula, Pinus cembra and Pinus mugo have a higher stem diameter to tree height ratio than the other species.As another important issue, we changed the parameter of needed growing degree sum required for full leaf cover of Betula Pendula and Larix decidua to account for their comparatively fast budburst (Murray et al., 1989).
Another important change concerns the parameter of maximum 20 yr coldest month mean temperature for establishment (TMAX est ).This limit is not associated with plant physiological functions, but rather represents a surrogate for functions not implemented in LPJ-GUESS that are responsible for outcompeting cold adapted species under warmer climates.We removed this limit for all species but Pinus cembra.In a previous version, the implementation of TMAX est values of −1.5 (Miller et al., 2008;Hickler et al., 2012) or −2 (Koca et al., 2006) for Picea abies, and of −1 (Koca et al., 2006;Miller et al., 2008;Hickler et al., 2012) for Pinus sylvestris led to two discontinuities (results not shown).First, both species should have been able to establish in the lowest of the analyzed elevations (ca.800 m), which is reflected by the NFI data.But maintaining the limits led to a prevention of establishment of both because the climate is "too warm".Second, if the climate in a stand was near the limit of the TMAX est , a species did not grow during the spinup (1901-1930 climate) but established in the slightly colder 1940s.Although it became much warmer afterwards, the cold-adapted species, once established, did not become extinct (e.g., despite drought stress).This created the inverse picture of cold adapted species that grow better under warmer climates.A future application of this parameterization in regions with warmer winters might need further tests and possibly a reimplementation of TMAX est for some species.
Similarly, the drought tolerance parameter used covers only a part of plant physiological responses to drought.There is a high risk that its values will be defined varying from realistic values to also cover other plant-related effects that are not included in LPJ-GUESS (e.g., plant water storage, plant water conduction traits, certain stomata closure effects or absence caused by dispersal barriers).This complicated the parameterization of the drought tolerance.In accordance with the occurrence of species in the NFI data and Bugmann (1996) and Lischke et al. (2006b), Picea abies and Betula pendula are more drought-resistant than the drought tolerance values of Hickler et al. (2012) might reflect.Consequently, we increased the drought tolerance of Picea abies and Betula pendula (decreased d tol value, Table D3).For the same reason we decreased the drought tolerance of Abies alba, Quercus robur and Fagus sylvatica (increased d tol value, Table D3).We included C3 grass as a plant functional type into our modeling without changes to the existing LPJ-GUESS functions.We did not use a bole height to calculate the daily fraction of incoming photosynthetically active radiation (PAR).In other words, the leaf area of all species and the foliage layers to calculate the PAR are equally distributed vertically from ground to treetops.This could be a problem for the modeling of species that produce foliage high above the ground (e.g., Pinus sylvestris), and by that might have an advantage because they are less shaded.However, bole height is parameterized in LPJ-GUESS while in reality it shall be dependent on stand density and tree age, and therefore treated as a state variable.
Other broadleaved  2).See the lower right chart for descriptions.Note that for some species the two inner quartiles of the NFI data are located at 0 in all plots.
Fig. D1.Comparison of the NFI1 data with LPJ-GUESS runs.LPJ-GUESS results were produced using 400 replicate patches and altitudespecific disturbances (see RID in Table 2).See the lower right chart for descriptions.Note that for some species the two inner quartiles of the NFI data are located at 0 in all plots.Table D3.Specific tree parameters.b: boreal; t: temperate; st: shade-tolerant; ns: nearly shade-tolerant; ist: intermediate shade-tolerant; si: shade-intolerant; e: evergreen; s: summergreen; d: summergreen with decelerated senescence; phenramp: growing degree sum on 5 degree base required for full leaf cover; k latosa: ratio of leaf area to sapwood cross-sectional area; rootdist u and rootdist l: proportion of roots extending into upper and lower soil layer; chill b: changed chilling parameter (Sykes et al., 1996); d tol: drought tolerance, lower values show higher tolerance (minimum soil water content needed for establishment, averaged over the growing season and expressed as a fraction of available water holding capacity, and water uptake efficiency); gdd5min: minimum growing degree day sum on 5

Fig. 1 .Fig. 1 .
Fig. 1.Principle of GAPPARD.Development of any vegetation state variable (e.g.biomass of a species) y over time t.a) solid line: average development with disturbances; dashed line: development without disturbances.b1-bn) development of the state variable for patch 1 to n, stand-replacing disturbances appear with disturbance probability p. c) necessary information to calculate y with the GAPPARD method at time T. For years x1 to xT−1 the same development of y is applied.

Fig. 2 .
Fig. 2. Principle of GAPPARD applying changing drivers.a) dashed lines: developments of two subsequent

Fig. 3 . 42 Fig. 3 .
Fig.3.Terrain around and location of the modeled stands a-h.See Table2for detailed values of the stands' altitude.

Fig. 6 .
Fig. 6.Simulation results for stands (a), (d) and (h) for the simulation years 1900 and 2080.Carbon mass per tree species is plotted against tree height classes for (I) LPJ-GUESS using 400 replicate patches and (II) LPJ-GUESS-G (GAPPARD method used on LPJ-GUESS).Height class 1: trees of 2-6 m height; height class 2: trees of 6-10 m height, and so on.

G
Fig. D.1.Comparison of the NFI1 data with LPJ-GUESS runs.LPJ-GUESS results were produced using 400 replicate patches and altitude specific disturbances (see RID in Table2).See the lower right chart for

G
Fig. D.2.Comparison of the NFI3 data with LPJ-GUESS runs.Note that for some species the two inner quartiles of the NFI data are located at 0 in all plots.See Fig. D.1 for description.

Fig. D2 .
Fig. D2.Comparison of the NFI3 data with LPJ-GUESS runs.Note that for some species the two inner quartiles of the NFI data are located at 0 in all plots.See Fig.D1for description.

Fig. D3 .
Fig. D3.Single-species carbon mass with LPJ-GUESS using 400 replicate patches.Development of carbon mass until 2080.Bars indicate the NFI1 (1985) and NFI3 (2006) data.Black bar sections stand for broad-leaved species that were not modeled.
Fig. D5.Simulation results for stands (a)-(d) for the simulation years 1900, 2000 and 2080.Carbon mass per tree species is plotted against tree height classes for (I) LPJ-GUESS using 400 replicate patches, and (II) using the new GAPPARD method against LPJ-GUESS.Height class 1: trees 2-6 m height, height class 2: 6-10 m, and so on.

Table 2 .
Specific characteristics for stands a to h.Lat: latitude in Swiss coordinates (CH1903/lv03 projection; longitude is 638 300 m for all stands); Alt: altitude above sea level; Temp: mean annual temperature; Prec: sum of precipitation of main growing period (April-September); NFI-alt: altitudinal range set for NFI plots associated with the modeled stands; NFI1, NFI3: numbers of plots per altitudinal range; RID: return intervals for generic, patch-destroying disturbances.

Table D2 .
Climatic range parameters.The affiliations to species are shown in TableD3.

Table D5 .
Forest biomass of the NFI3 inventory and LPJ-GUESS results for the simulation year 2006 (all in t ha −1 ).See TableD4for descriptions.