Joint simulation of stand dynamics and landscape evolution using a tree-level model for mixed uneven-aged forests

Abstract

• Context

A century and more after a major reforestation program, large areas in the French Southern Alps have moved to a landscape mosaic of old pine plantations and new, heterogeneous and uneven-aged, mixed stands. These conditions are challenging foresters in silvicultural practices and management choices.

• Aims

The aims of this study are to understand, analyze, model, and simulate the ongoing phenomena, and to propose a decision-making tool.

• Methods

An individual-based forest dynamics model considering recruitment, growth, and mortality (as related to the spatial arrangement of stands and species, to site conditions and competition) and a simulation system including spatial sampling are designed and calibrated to allow simulation of both silviculture treatments at the stand level and management strategies at the forest or landscape level.

• Results

By keeping track of the trees while simulating at the forest level, they offer an alternative to upscaling strategies and a suitable tool for prediction of stand and forest characteristics in situations influenced by strong driving forces such as colonization and forest maturation.

• Conclusion

This approach is a straightforward means for adjusting forest management to trends such as expansion of shade-tolerant species; as spatial and temporal variation in site conditions are accounted for, it is also a promising way towards predicting their warming-induced upward shift.

1 Introduction

Modeling and simulation are commonly used to predict the evolution of a forest stand according to different silvicultural options or to estimate the possible evolution of larger areas such as a whole forest, a hillslope, a watershed or any part of a forest landscape subject to dynamics such as succession, colonization and migration, or decline induced by climate change (Pretzsch et al. 2008; Mendoza and Vanclay 2008; Porté and Bartelink 2002). In multi-layered forests with regeneration appearing under a canopy, or in any context with trees of very different sizes and ages, the average-tree concept at stand scale has no real sense. Thus, when considering such a complex size and age structure at the forest or landscape scale, the challenge is to allow simulation over large areas while using an individual-based model able to account for local structure either natural or human-induced. The «gap model» approach (Bugmann 2001) is the most commonly used modeling strategy for this purpose. But, despite the increase of computer power, it is unworkable, and also probably useless, to simulate the recruitment and growth of every individual, either a seedling or a mature tree, all over a very extensive area. Applying upscaling methods to go from tree and stand level to forest or landscape level has been attempted in some cases. Scaling issues in the frame of forest modeling are discussed and illustrated in detail by Urban et al. (1999), Bugmann et al. (2000), and Urban (2005). These approaches use more or less complex modeling techniques often leading to a meta-model derived from data simulated by a fine-scale model. They require much care in order to avoid incorrect simplifying assumptions that can lead to a departure from the fine-scale dynamics and behavior, and from the predictive ability of the tree-level model (Bugmann et al. 2000; Fleming et al. 2002). Thus, while upscaling in landscape ecology and forestry appears to be a specialized field of study in itself, its outcomes do not maintain a direct link between tree level and forest level. Therefore, understanding any coarser-resolution evolution through looking at tree evolution is not possible unless getting back to the original tree-level model (Urban et al. 1999). Moreover, specifying silvicultural treatments involving fine management of regeneration in mixed, uneven-aged stands and simulating their effect on a wide scale might be troublesome, if not impossible. Finally, Urban et al. (1999) suggest that the way the tree-level model will be simplified towards a large-scale meta-model depends on the application and the upscaling procedure will therefore lead to several application-specific meta-models (e.g., scope mainly on the successional dynamics, or on timber management).

As an alternative to upscaling, we present a combined modeling and simulation approach relying on an abstraction of the forest as a grid of cells like in gap models, and on using a flexible spatial sampling procedure. To our knowledge, such a spatial sampling system for simulation has not been reported in the literature. This new approach was elaborated in the context of Mediterranean mountains in the South-East of France where restoration and forest maturation resulted in a complex mosaic of various forest stands. Besides describing the structure and components of a new tree-level forest dynamics model designed for mixed uneven-aged forests in this context, this paper focuses on the approach developed to simulate the evolution of large areas (forest, hillslope) while keeping the tree-level information, thus allowing fully compatible simulations of silviculture treatments at the stand level and of management or planning strategies at the forest or landscape level. The paper also addresses some examples of potential application.

2 Materials and methods

2.1 Species, study area

Extended areas in the French Southern Alps were the scene of large scale reforestation at the end of the nineteenth century and during the first half of the twentieth century: soil erosion was controlled and forested areas devastated by overexploitation and overgrazing were successfully restored.

For several decades now, this restoration has moved into a second, dynamic phase, as European beech (Fagus sylvatica) and silver fir (Abies alba), two shade-tolerant species which had remained within residual stands, are recolonizing the planted pine (Pinus nigra, Pinus sylvestris, Pinus uncinata) stands (Dreyfus 2008). Such an evolution is also observed in other countries (e.g., silver fir regeneration under pine stands in Italy, Tonon et al. 2005, in Spain, Oliva and Colinas 2007). The outcome of this “maturation process” is a complex landscape mosaic (old pine plantations and new uneven-aged mixed stands) that opens management questions as well as interesting prospects for the associated biodiversity and for the sustainability of the forest itself as climate change has become a new threat (Cailleret and Davi 2010).

2.2 Data

Three datasets were used for modeling growth and mortality, each of them concerning either “seedlings” (height <1.30 m) and “saplings” (diameter at breast height, or dbh, <7.5 cm), or “trees” (dbh ≥ 7.5 cm):

• D1: for P. nigra, a network of thinning experiments and permanent plots, including advanced regeneration stage and old mature stands, surveyed by INRA for several decades within the Mediterranean range of this species in Southern France, provided a complete dataset of individual size and increments, and of mortality rates. It covers various site conditions, initial stand structures and thinning characteristics.

• D2: as a similar network is lacking for the other species studied, we used temporary plots of the French National Forest Inventory (IFN) over the whole study area (11 French departments under the influence of both mountain and Mediterranean climates). For every tree above 7.5 cm in dbh, the dataset provides a 5-year radius increment (increment core), as well as present diameter, bark thickness and height. The age is measured on two or 3 three dominant trees per plot. Diameter at stump height is given both for living trees and for trees that apparently died within the last 5 years. We selected the plots where one of the five species studied (including P. nigra) is predominant; plots with evidence of cutting within the last 5 years were rejected.

• D3: a survey design based on nine intensive sites (0.49 to 1.04 ha in area, each partitioned by a grid of 100 m² cells) established and monitored by INRA URFM, between 950 and 1,400 m a.s.l. on north and south facing sides of the mount Ventoux (south-east of France, 44° 10′ 28″ N, 5° 16′ 44″ E, WGS84), a typical situation demonstrating the phenomena of interest (colonization, recent forest decline). Each site represents a step in the maturation process from a pure pine plantation to a mixed stand with pine and beech/fir together in the adult dominant layer, through steps with mixed regeneration establishing and growing under the pine canopy. All overstorey adult trees and saplings were counted and measured (dbh) in each cell, while seedlings were counted in two height classes : <30 cm, 30-130 cm. Total height was also measured on a sample of 30 adult trees in each site. In each cell, at least five sample seedlings or saplings were described (species, height, dbh) and measured again 5 years later. Age at 30 cm from the ground was measured on a sample by coring adult trees and tall saplings, by counting the branch whorls on seedlings and other saplings of coniferous species (silver fir and the three pine species are monocyclic species with well marked annual whorls of branches), and, for beech, by counting rings on sacrificed individuals (Table 1).

  Table 1 Characteristics of growth and mortality survey datasets (D1, D2, D3)

  Two other datasets were used for calibrating the recruitment submodel:

• D4: a complete inventory (trees measured in dbh, seedlings/saplings counted by height class) was carried out by the French Forest National Service (ONF) on 214 temporary plots (basal area: 0.0 to 100.8 m²/ha, number of trees: 0 to 3,080/ha, of seedlings/saplings: 0 to 13,060/ha) each one composed of four or five 100 m² cells (total area = 12.0 ha). These temporary plots were established on the nodes of a grid covering the same elevation range as the INRA survey sites (D3) on both sides of the mount Ventoux.

• D5: GIS maps of stand composition and silvicultural regime, and of site conditions, were also made available by ONF for the managed forests on the mount Ventoux. The stands map indicates the proportion of the three main species and dominant age of the main species for each of 1,756 management units (11,608 ha of forest stands).

2.3 Outlines of the modeling method

The model has something in common with gap models in that it is individual-based, with space abstracted as a grid of cells (in this case, 100 m² cells), and tree position within the cell is not considered by the model (but it can be set in simulation for visualization purposes). The dynamics of a stand results from the recruitment, growth and mortality processes for individual trees, according to species, site conditions and competitive interactions with neighbors. For each process, the model includes a submodel based on ecological considerations (like decomposing growth in a potential growth component as a function of the species, of tree age and of site index, and reduction components related to local inter-tree competition and to individual tree status, see, e.g., Arney 1985) and on empirical approaches (fitting process comparing several shapes for each relation or component). The appropriate predictor variables are chosen according to a basic understanding of how trees and forest stands behave and how factors like site, age and competition affect growth and mortality. Among various expressions for each factor, the best variable was then selected according to the results of test statistics.

2.4 Submodels

2.4.1 Recruitment submodel

Recruitment of new individuals in a forest stand is the result of the spatial pattern of seed trees, of reproduction and seed production, dispersal by various agents, germination and establishment of young seedlings in more or less favorable conditions (ground cover and vegetation, predation, diseases). Seed production and dispersal are not considered explicitly in this recruitment submodel, but datasets D4 and D5 provided the opportunity to build a basic recruitment submodel which form is shared by the five species. For each species, the number of seedlings counted on D4 plots appeared to be related to the distance (obtained through D5 maps and a GIS) of the nearest mature stand (i.e., mean dbh of trees of the same species is over 15 cm and their age over 40 years; same criteria for the five species) which can be considered as the preeminent seed source. Considering the assumptions that fecundity is linearly proportional to the basal area of seed trees (Clark et al. 1999; Nanos et al. 2010), and that the contribution of a tree to seed dispersal is proportional to its fecundity (Sagnard et al. 2007), the data were standardized to a reference basal area value (arbitrarily set to 40 m²/ha) for the species of interest. For this basal area of reference, the recruitment equation is as follows (see also Fig. 1):

$$ Log({\text{smallsdgs}} + 1.1) = {\text{r1}} + {\text{r2}} \cdot {\left( {{\text{neardist}} + 15} \right)^{{ - {\text{r3}}}}} $$

(1a)

with a slightly different form for F. sylvatica:

$$ Log\,({\text{smallsdgs}} + 1.1) = {\text{r1}} \cdot {e^{{ - {\text{r2}} \cdot {{\left( {{\text{neardist}} + 15} \right)}^{\text{r3}}}}}} $$

(1b)

where smallsdgs is the number of seedlings with height <0.3 m on a 100 m² cell (between about 0.1 and 0.3 m as seedlings of the current year are excluded) and neardist (m) the distance to the nearest mature stand of the same species; r1–r3 are parameters to be estimated. The small seedlings counted on a plot originate in fact from several annual dispersal and germination events: an annual number of new recruited seedlings is then estimated according to mean age of small seedlings for each species (mean values obtained from D3, data not shown). The terms 1.1 and 15 in Eqs. 1a and b are fixed values intended to control the model behavior when there is no seedling or when the nearest stand is around the plot (neardist = 0). This relationship implies that the other stands and the possibly interspersed seed trees (between the target plot and the nearest stand) are not considered; therefore, they contribute to a high residual variability which is modeled as a stochastic component (not presented) intended for random drawing during simulation runs.

Fig. 1

Number of small seedlings (height <0.3 m) recruited annually on a 100 m² cell as a function of the distance to the nearest source stand of the same species (see Eqs. 1a and b). Each curve is drawn to the maximum distance observed (ca. 500 m for Pinus nigra, 1,500 m for Fagus and 2,000 m for Abies)

As recruitment involves sources and targets distributed over the whole forest, it is referred to below as “flows” or “recruitment flows”.

2.4.2 Diameter growth submodel for trees

For adult trees of the pine species (D1 and D2), diameter increment is expressed as a function of a dominant height increment (growth potential term), a local stand density and a tree competition status (in two reduction terms):

$$ \begin{array}{*{20}{c}} {{\text{DINC5}} = 0.1 \cdot {\text{a1}} + \frac{\text{a2}}{{1 - {\text{a3}} \cdot Log\left( {\text{HdomINC5}} \right)}}} \\ { \cdot \left( {1 - {\text{a4}} \cdot {e^{{ - {\text{a5}} \cdot {\text{NB}}{{\text{A}}^{{ - {\text{a6}}}}}}}}} \right)} \\ { \cdot \left( {1 + {\text{a7}} \cdot {\text{ICS}}} \right)} \\ { \cdot min\left[ {1.025,\left( {1 + {\text{a8}} \cdot \left( {{\text{Hdom50}} - 15} \right)} \right)} \right]} \\ \end{array} $$

(2a)

where: DINC5 (cm) is the 5-year diameter increment.

HdomINC5 (m) is an individual tree potential height increment computed using a site index curve (for the species considered) and individual tree age at the beginning (t) and at the end (t + 5 years) of the increment period. It is an estimate of the growth of the tree, at present age and in these site conditions, if it would have grown as a dominant tree (i.e., in the upper part of the canopy of an even-aged stand) since its birth; this value is taken as a reference.

NBA (m²/ha), the neighboring basal area, is the basal area of the plot (in datasets D1 or D2, a few hundreds of m²).

ICS (cm/m), a standardized individual competition status for each adult tree, is obtained by dividing tree diameter by the dominant height on the plot. Arney (1985) proposed this ratio as a taper index standardized among all trees in the stand as well as among stands of various age and site. It is positively related to live-crown ratio which is an another good indicator, although not so easily observable and predictable, of competition experienced by a tree and of its ability to withstand present competition.

Hdom50 (m) is the local site index, i.e. the estimated stand dominant height at age 50.

a1–a8 are parameters to be estimated.

For adult trees of Abies and Fagus (D2 dataset), two shade-tolerant species, the site index can hardly be estimated because the present tallest trees may have survived and grown for years under other trees that now have disappeared. Therefore, Hdom50 and individual HdomINC5 as defined for Eq. 2a are not considered reliable enough. The tree mean height growth (since its birth) was then used as an individual tree growth potential; it is depending on site conditions and tree age, but also on competition experienced by the tree along its life.

The form of the equation is slightly different for the two species:

A. alba:

$$ \begin{array}{*{20}{c}} {{\text{DINC5 }} = 100 \cdot {\text{b1}} \cdot 1 / \left( {1 - {\text{b2}} \cdot Log\left( {H/{\text{Age}}} \right)} \right) \cdot {\text{Ag}}{{\text{e}}^{{ - b3}}}} \\ { \cdot 1/\left( {1 - {\text{b4}} \cdot {e^{{ - {\text{b5}} \cdot {\text{NBA}}}}}} \right)} \\ { \cdot \left( {1 - {e^{{ - {\text{b6}} \cdot {\text{ICS}}}}}} \right)} \\ \end{array} $$

(2b)

F. sylvatica:

$$ \begin{array}{*{20}{c}} {{\text{DINC5}} = 100 \cdot {\text{c1}} \cdot {{\left( {1 - {e^{{ - {\text{c2}} \cdot H/Age}}}} \right)}^{\text{c3}}} \cdot {\text{Ag}}{{\text{e}}^{\text{c4}}}} \\ { \cdot 1/\left( {1 + {e^{{{\text{c5}} \cdot {\text{NBA}}}}}} \right)} \\ { \cdot \left( {1 - {e^{{ - {\text{c6}} \cdot {\text{ICS}}}}}} \right)} \\ \end{array} $$

(2c)

where H (m) is tree height and Age (years) is tree age. The mean past height growth (H/Age) indicates a potential for further growth whereas the two others terms indicate reduction by current competition. b1–b6 and c1–c6 are parameters to be estimated.

2.4.3 Diameter growth submodel for seedlings and saplings

Diameter increment is expressed as a function of a cumulated growth lag, relative tree height and three competition indices. The two first terms reflect a growth potential and the other terms indicate reduction by competition:

$$ \begin{array}{*{20}{c}} {{\text{DINC5}} = 5 \cdot 0.1 \cdot {\text{d1}} \cdot {e^{{ - {e^{{{\text{d}}2 - {\text{d}}3 \cdot {\text{Dlag}}}}}}}}} \\ { \cdot {e^{{ - \left( {{\text{d4}} - {{H} \left/ {{H\max }} \right.}} \right)}}}^2} \\ { \cdot 1/\left[ {1 + {\text{d5}} \cdot (({\text{SDLp}} + 1)/100)} \right]} \\ { \cdot 1/\left[ {1 + {\text{d6}} \cdot {{(({\text{SDLfa}} + 1)/100)}^{\text{d7}}}} \right]} \\ { \cdot 1/\left[ {1 + {\text{d8}} \cdot {{(({\text{NBA}} + 0.01)/60)}^{\text{d9}}}} \right]} \\ \end{array} $$

(3)

Dlag (between 0 and 1) is the DBH/Dthmax ratio with Dthmax (cm) the dbh a tree would have in a virtual, “ideal”, situation when its height is the height (hdomi) of a dominant tree without competition, at the same age and in the same site conditions. ${\text{Dthmax}} = {100} \cdot {\text{hdomi / HoDBHmin}}$, where HoDBHmin (cm/cm) is the observed minimum value (species-specific) of H/DBH ratio for small trees. Although Dthmax is a theoretical value, Dlag allows to compare saplings for their deviation from this reference.

Hmax (m) is the maximum tree height within the same 100 m² cell.

SDLp (cm/100 m²) is a local index computed as the sum of the diameters of all pines (adult, saplings, seedlings, of P. nigra, P. sylvestris, and P. uncinata) with dbh greater or equal to the dbh of the subject, within the same 100 m² cell only.

SDLfa (cm/100 m²) is the same as SDLp, but considering larger individuals of F. sylvatica and A. alba instead of pines.

NBA (m²/ha) is the neighboring basal area, on a group of nine cells: the one bearing the tree and the eight neighboring cells.

d1–d9 are parameters to be estimated.

Since the seedlings were studied in more or less dense stands, not in large open areas more suitable for light demanding species, dataset D3 contains only few diameter increment data for pines, except for P. nigra. For P. sylvestris and P. uncinata a more simple relationship was therefore used:

$$ \begin{array}{*{20}{c}} {{\text{DINC5}} = {\text{(Dthmax5 - Dthmax)}}} \\ {{\text{.(1 - exp( - 0}}{.01}{\text{.RelSpcg))}}} \\ \end{array} $$

where Dthmax is computed at the beginning of the increment period (i.e. at Age t), and Dthmax5 is computed at the end of the increment period (i.e. at Age t + 5 years).

RelSpcg = 10,746/(Hmean·Ndbh^0.5) is a relative spacing index similar to Hart–Becking index, with Ndbh, the number(/ha) of individuals taller than 1.3 m in the same cell, and Hmean, the arithmetic mean height of these individuals.

2.4.4 Height–diameter submodel for trees

This relationship is based on the height-diameter curve equation proposed by Mønness (Omule and MacDonald 1991). It is conditioned to go through the point defined by the dominant diameter and dominant height of the stand. We related the s parameter (i.e. shape of the curve) to RelSpcgMod, another index derived from Hart–Becking spacing index:

$$ H = {\text{Hdom}} \cdot \frac{{\left( {1 - {e^{{ - s \cdot {\text{DBH}}}}}} \right)}}{{\left( {1 - {e^{{ - s \cdot {\text{Ddom}}}}}} \right)}} $$

(4)

where H (m) and DBH (cm) are tree height and dbh, and

$$ s = {\text{s0}} + {\text{s1}} \cdot {\text{RelSpcgMod}} $$

with

$$ {\text{RelSpcgMod}} = {\text{10746/(Hdom}} \cdot {\text{Ndb}}{{\text{h}}^{{{ 0}{.25}}}}) $$

Hdom (m) and Ddom (cm) are the dominant height and dominant diameter of the plot, and Ndbh, the number (/ha) of stems with height above 1.3 m, within the plot; s0 and s1 are parameters to be estimated.

2.4.5 Height growth submodel for seedlings and saplings

Height growth is a crucial process by which small individuals can come into higher storeys. This relationship predicts the ratio of tree height increment HINC5 to the potential growth HdomINC5 of a similar (same species, age and site conditions) virtual tree having grown in dominant position since its birth. Comparing seedlings (dataset D3) under dense or partly opened canopies, we observed that the cumulated lag in terms of height was very unfavorable to subsequent height growth. Two variables indicate either “one-sided”, asymmetric (SDLfa) or “two-sided”, symmetric (NBA) competition, in two growth reduction terms. The equation is as follows (see also Fig. 2):

$$ \begin{array}{*{20}{c}} {HINC5/{\text{HdomINC5}} = } \hfill \\ {{1/[1} + {\text{t1}} \cdot {\text{Hla}}{{\text{g}}^{\text{t2}}} \cdot {\text{((1 - H/Hmax}}{{)}^{\text{t3}}} + {0}{.01)]}} \hfill \\ { \cdot {1/[1} + {\text{t4}} \cdot {{{\text{( (SDLfa}} + {1)/100)}}^{\text{t5}}}{]}} \hfill \\ { \cdot {1/[1} + {\text{t6}} \cdot {{{\text{( (NBA}} + {0}{.1)/50)}}^{\text{t7}}}{]}} \hfill \\ \end{array} $$

(5)

where Hlag is the lag of present tree height as a fraction of hdomi: Hlag = (hdomi − hi)/hdomi.

Fig. 2

Negative effect of local asymmetric competition by Abies and/or Fagus larger individuals (canopy trees and taller saplings i.e. with dbh greater or equal to the dbh of the subject; SDLfa is the sum of the diameters of these individuals, within the same 100 m² cell) on the height growth of a sapling of either Abies, Fagus or Pinus nigra. The sapling height increment is expressed in relative value (as divided by its potential height growth, HdomINC5, see Eq. 5). Two basal area levels (NBA, neighboring basal area for the cell bearing the tree and the eight neighboring cells) are considered : 15 m²/ha (thin lines) and 40 m²/ha (thick lines); Pinus nigra regeneration is sensitive also to this factor. In the case presented, sapling height (H) is fixed to 3.5 m with a lag of 30% (Hlag = 0.30); maximum tree height (Hmax) on the same cell is 15 m

t1–t7 are parameters to be estimated.

2.4.6 Mortality submodel

The mortality submodel deals with competition-induced mortality. The equations were calibrated using mortality rates observed on five-year periods (D2 for adult trees, D3 for seedlings/saplings). The probability of mortality, probMort5, was fitted as a function of local stand density NBA and of relative tree size.

For trees:

$$ {\text{probMort5}} = ({\text{m1}} + {\text{m2}} \cdot {\text{NBA}}) \cdot {e^{{({\text{m3}} + {\text{m4}} \cdot {\text{NBA}}) \cdot {\text{DBH}}/{\text{Ddom}}}}} $$

(6)

where DBH is tree dbh and Ddom is the dominant diameter of the plot.

For seedlings/saplings (see also Fig. 3):

Fig. 3

Effect of individual competition status (local relative height H/Hmax) and local basal area (two levels of NBA: 40 m²/ha, thick lines, 15 m²/ha, thin lines) on the probability of mortality for a seedling/sapling (see Eq. 7)

$$ {\text{probMort5}} = 5 \cdot \left[ {1 - {e^{{ - {\text{mo}} \cdot {\text{NB}}{{\text{A}}^2} \cdot [{\text{l/min}}(l,H/H\max ) - 1]}}}} \right] $$

(7)

where H is tree height and Hmax is the maximum height of any individual within the same 100 m² cell.

m0–m4 are parameters to be estimated.

2.5 Simulation methodology

The model is embedded in some modules of the Capsis simulation platform (http://www.inra.fr/capsis, Dufour-Kowalski et al. 2011). The length of the simulation period can reach 100 years or more. The simulation step is set to 5 years in accordance with individual increment and mortality rate data used for submodels calibration.

A cohort approach (Bugmann 1996) is used: all trees of a given species established in given year and cell are assumed to be identical. In the simulator, a cohort of N identical trees is represented by only one tree instance, with a “number” attribute equal to N, and its evolution is computed only one time, not N times. Mortality within a cohort is simulated by reducing N according to the probability predicted by the submodel. As this probability is a decimal value between 0 and 1, the number of individuals computed for removal from the cohort is usually not an integer value. Therefore, a random drawing is done on the decimal part: for example, if M = 4.17 trees should be removed in the cohort, N will be reduced by 5 if a uniform random drawing in [0, 1] yields a value less than 0.17, or by 4 otherwise. The same way is used for getting the initial number of a new cohort as the value predicted by the recruitment submodel is usually not an integer value.

2.5.1 Input data and model initialization

In order to run a simulation, information about the stands and the forest, and about the site conditions, is required. Within the frame of forest management, tree and stand information is usually available in the form of forest compartment maps and of tables of usual characteristics: stand composition (species) and structure (high forest, coppice, even-aged, or uneven structure for age and size), stems/ha, top height, dominant age, approximate basal area, sometimes with subdivision according to a few layers (canopy, understorey, regeneration). Except when simulating the evolution of a single stand with a complete inventory (in dbh) and some height measurements, a preliminary step aimed at completing tree-level and cell-level information must be performed before launching the simulation itself. The age of trees, saplings and seedlings is estimated through statistical size-age relationships or distributions fitted in the study region.

As mentioned above, the growth submodels take site conditions into account through the conventional forestry concept of site index (dominant height reached at a given reference age—50 years in this study—for a pure, even-aged stand, depending on species and site) which value is used in these submodels for computing a potential or a reference for diameter and height growth of any tree. Despite much research effort is devoted to directly linking site factors to tree and forest growth, site index still remains the most available reference for the species and site conditions in the study area. Depending on the forest and the species, site index values (one for each species) used in our model are obtained (and assigned to each 100 m² cell through reading the corresponding maps) in two alternative ways: (1) either through traditional site type maps drawn up by foresters using site classification keys (i.e. each type is identified by topography, soil, and/or indicator plants) and completed by correspondence tables linking the site types to site index levels for each species of interest, (2) or directly by means of specific relationships predicting site index from environmental information given by a DEM, by raster maps interpolated from spatial climate databases and soil/geologic maps, allowing for building site index maps through a GIS. The spatial variability of site conditions is distributed in patches in the first case (all cells in the same patch have the same site index) and to pixels in the second case (with the resolution of the environmental information layers) allowing for more progressive variation between cells. In the “traditional site maps” case, the climate is usually taken into account indirectly through aspect and altitudinal limits. In the “site index raster maps” case, the climate is taken into account either directly if it belongs to the set of layers used to estimate site index, or indirectly through altitude and aspect maps. In the study area, most, if not all, approaches dealing with site conditions consider elevation and aspect, two factors linked to climate (especially temperature) through standard relationships.

2.5.2 Subsequent stages in each simulation step

At the beginning of the simulation step, all cells are inspected in order to check which ones are currently possible sources for dispersal: in fact, trees may have reached an age and a size which make them become seed trees and thus make their cell become a seed(ling) source; the power of this source can increase or decrease as its basal area in seed trees increases due to growth or is reduced by a thinning or by mortality. For each species, a flow between each potential target cell and its present nearest source cell is then computed; in order to obtain the number of seedlings in each new cohort generated, the value predicted by the recruitment equations (Eqs. 1a and b) for a source having 40 m²/ha in basal area (reference value, see 2.4.1) is adjusted according to the actual basal area of seed trees in the source cell. Growth and mortality are then simulated.

When simulating for a whole forest composed of many management units, thinnings cannot be designed in interactive mode for each of them. Some management rules, depending on the type of stand and on its development stage, are selected or defined by the user at the initialization and are applied to each unit (i.e. to the cells and trees representing this unit) at the end of each simulation step. A thinning (or regeneration cutting, etc.) is triggered and simulated only when the management unit has reached some threshold(s) selected or defined by the user; otherwise, no cutting is simulated for this unit until the next step and check. Damages to seedlings involved by cuttings (felling, skidding) are simulated in a simplified way according to our observations (unpublished data).

Finally, every kind of cell, stand and forest attributes are updated according to new tree dimensions and numbers.

2.5.3 Spatial sampling schemes

In a context of colonization or migration, predicting the evolution of a large area cannot be done by simulating the evolution of each stand independently; dispersal and regeneration processes require that surrounding stands in a wide range (i.e., the forest all around) are taken into account, as well as the evolution of each of these stands. A major feature in our approach is the ability to simulate the evolution of a whole forest by considering the evolution of each stand thanks to a limited sample of cells, as an alternative to a complete grid laying over the whole forest. Each stand can be abstracted according to several options: (a) a single group of nine cells located near its center, (b) several groups (of nine cells) evenly distributed over a network, more or less dense since the distance between groups can be varied at will (c) strips of cells (three-cell wide) crossing the stand in one (or two orthogonal) direction(s), (d) a complete grid of cells. The option can be different for each stand in the forest and is mentioned in the input file (text file or ESRI shapefile spatial data format). This simulation feature is aimed at decreasing time and memory requirement. Moreover, it allows for devoting more simulation effort to some particularly interesting parts of the forest: colonization front, zone of forest decline and regression, stands with high economic or ecological value.

In order to evaluate the spatial sampling approach and to compare different sampling schemes, a set of simulations were run on a rectangular (2 km × 1 km) pure stand, with uniform site conditions, which is progressively colonized by an invasive species starting from an adjacent stand of this second species (0.05 km × 1 km, Fig. 4). The colonization front is 1-km wide with a maximum range of 2 km. Three spatial sampling schemes are compared to the complete grid (“CG”): nine-cell groups every 100 m (“100S” sampling), or every 200 m (“200S” sampling) and orthogonal transects (transect sampling scheme “TRS” sampling; two transects, three-cell wide; the longest is along the colonization direction). Twenty simulations were run for each of these four options. In each simulation run, the variability between cells at the same distance (i.e., the shortest distance to the source stand) is due to random effects both in initial stand creation (slightly different at each run) and in simulation of the different processes in the model. For each repetition, nine-cell groups characteristics were averaged at distance intervals of 200 m (where all of the four schemes have cells) from the source stand. The remaining variability of these mean values between runs is analyzed using Tukey’s HSD test, a single-step multiple comparison procedure and statistical test, used in conjunction with a one-way analysis of variance (“aov”, then “TukeyHSD” functions in R software, R Development Core Team 2011) with the spatial option (three sampling schemes + complete grid) as the single factor (four levels, with 20 observations at each level). This test is performed to find which means are significantly different from one another, independently for each distance (every 200 m, from 200 to 1,800 m, starting from the source stand). Main attention was paid to the comparison between the complete grid (as the reference) and each of the three sampling schemes.

Fig. 4

Simulation of silver fir colonization in the frame of an evaluation trial (see also Table 2) with four spatial options. Colonization progress after 50 and 90 years is shown for one of the twenty simulation runs for each spatial option. The source, on the left (0.05 km × 1 km), is a silver fir stand (initial characteristics: basal area 23 m²/ha, 1980 stems/ha, age 40 years, dominant height 11 m); the colonized stand (2 km × 1 km) is a Scots pine stand (trees not shown, initial characteristics: basal area 10 m²/ha, 300 stems/ha, age 40 years, dominant height 11 m). Three spatial sampling schemes are compared to the complete grid (CG; 20,000 cells; cell limits are not shown): nine-cell groups every 100 m (100S sampling; 1,711 cells), or every 200 m (200S sampling; 435 cells) and orthogonal transects (TRS sampling; two transects, three-cell wide; the longest is along the colonization direction; 888 cells)

3 Results

3.1 Forest dynamics model

Parameter estimates and ordinary least squares regression statistics (NLIN procedure of the SAS software package, SAS Institute, Cary, NC) for each submodel and species are given in the Appendix. Three results are illustrated briefly hereafter in order to show differences among species that can strongly influence the evolution of the stands and of the whole forest.

Seedlings of fir or beech have been found in plots farther than 1,000 m from any fir or beech stand, respectively, while pine seedlings are much less abundant and almost lacking out of pine stands. This situation shows up in the calibrated recruitment submodel (Fig. 1, Eqs. 1a and b).

Figure 2 illustrates the strong negative effect of deep shade species (Abies and Fagus) on the height growth of their own regeneration and even more on pine seedlings (Eq. 5).

In addition, Abies and Fagus appeared to be less sensitive than P. nigra to the lag in height (Hlag) cumulated since tree birth (not shown).

Figure 3 shows the negative effect of competition pressure (neighboring basal area, NBA) and the positive effect of individual seedling/sapling status (local relative height) on survival probability. Competition-induced mortality rate (Eq. 7) is quite different among species, Abies and Fagus being far less sensitive than the two pine species shown.

3.2 Simulating colonization and migration

Figure 4 deals with a simulation for one of the evaluation trials (cf. Section 2.5.3) and Table 2 provides statistics about the twenty simulations that were run for each spatial option in this trial. The source is a silver fir stand; the colonized stand is a Scots pine stand with a rather clear canopy. Four variables are computed for each cell in the pine stand to indicate silver fir colonization progress: basal area of silver fir (BAFir), number (N3) and basal area (BA3) of fir saplings with height over 3 m, basal area of mature silver fir (BAM, age over 40 years and dbh over 15 cm). Table 2 shows the values at the end of the simulation period i.e. after 90 years, and for an intermediate stage (50 years) when some of the new silver fir trees in the pine stand become seed trees themselves (BAM > 0) and the initial fir stand is no longer the only possible source for new regeneration. Each of these four variables follows the same general pattern of sharp decrease with distance from the fir stand. The values indicate that spatial and temporal patterns are quite similar among all spatial options. The range and the standard deviation in the results of the 20 simulation runs is higher for the three sampling schemes than for the complete grid: the reason is that these sampling options correspond only to 9% (100S), 2.25% (200S) and 4.4% (TRS) of the complete area in term of number of cells (or of stand area covered by these cells), and the spatial differences between two runs can therefore be higher than for the complete grid option. For TRS, the variability of the values at each distance is especially high (as indicated by min, max and standard deviation) because, for each simulation run, there is only one group of nine cells at each distance (except at 1,000 m, which is the position of the orthogonal transect). According to the statistical tests, there is no significant difference for 100S sampling compared to the complete grid option, whatever the variable of interest, and this spatial sampling scheme is therefore able to represent the invasion pattern and speed in the same way as the complete grid. Some significant biases are observed for the 200S at 200 m from the fir stand; and for the TRS at 200 m and some other distances.

Table 2 Silver fir colonization under a Scots pine canopy after 50 and 90 years: comparison of 3 spatial sampling schemes with the complete cell grid

In terms of computing time and memory used, for this situation (2 km², 18 simulation steps) while the complete grid option requires 40 s and 837 MB of memory (for 70,000 cohorts at the beginning and almost 200,000 cohorts at the end), the requirement are only 5 s and 97 MB (from ca. 11,000 to 16,000 cohorts at the end), 3 s and 52 MB (from ca. 15,000 to 27,000 cohorts), 4 s and 68 MB (from ca. 13,000 to 19,000 cohorts) for the 200S, 100S, and TRS sampling schemes, respectively. These values are obtained with a generic but rather powerful personal computer (64-bit multi-core processor, 2.20 GHz, with multithreading, and 8 GB of RAM). Memory requirement is in fact the most limiting factor and the computing time may increase dramatically when the memory used is close to the maximum amount of RAM available. For all the pine, beech and silver fir stands of mount Ventoux (that is for ca. 65 km²), 200S scheme (17,730 square 100 m²-cells, 186,000 cohorts at the beginning) requires 75 s and 3.1 GB for 18 simulation steps (without forest management), 100S scheme (56,360 cells, 593,000 cohorts at the beginning) requires nearly 300 s and 6.7 GB, while for the complete grid (ca. 650,000 cells), even the initial step cannot be generated with 8 GB of RAM.

This situation was chosen in order to have a colonization strong and fast enough to permit the detection of differences between spatial options (complete grid and sampling schemes). Several tests (not presented) in other situations (e.g., fir under a more dense pine stand or a beech canopy, or with beech as the invasive species) demonstrated lower overall colonization speed with similar results in terms of spatial sampling scheme effect.

Figure 5 presents simulations starting from a real situation and dealing with natural colonization in Scots pine stands by A. alba present in an adjacent mature stand. Some spatial sampling options are shown at the same time. The diagrams in Fig. 5a show the progression of the fir under the pine canopy. Though pine adult trees are present in the fir stand, pine regeneration is not able to grow under the fir canopy; it is not able either to develop markedly within the pine stand. Figure 5b maps show that spatial sampling can represent the colonization phenomenon both for spread and for local variability. Every characteristic of the trees appearing on these maps or diagrams can be used to compute any desired statistics (mean values, dimensional or spatial variability and structure, within each species and/or between species) regarding trees, cells, stands or the whole area. These values are useful for reporting their simulated evolution and are obtained directly by averaging cell-groups characteristics in a way that mimics a usual systematic statistical forest inventory; they could be obtained as well by using first any kind of spatial interpolation to account for possible non-linear distributions in space (such as a regeneration gradient within a stand). It should be noticed that the way these summary characteristics are computed and reported has no influence on the subsequent simulation steps as only characteristics of the (instantiated) trees and cells are used for simulating dispersal and regeneration, and growth or mortality.

Fig. 5

Simulation of fir colonization in pine stands. a Diagrams of a 400 m long transect showing fir (in dark) regeneration spread and height growth within a pine (in bright) stand neighboring a fir stand, at time t, t + 30 years, t + 60 years. b On these maps (at t + 60 years), pine trees are hidden in order to show more clearly the progression of the fir in the pine stands, covered either with a complete grid of cells (cell width is 10 m), with nine-cell groups (60 m between groups in both directions), or transects of cells (stand on the left) and nine-cell groups (40 m between groups, stand on the right). The white rectangle represents the part shown in Fig. 5a

In a way to illustrate another promising application of the model, Fig. 6 presents a simulation where site conditions are evolving. The expected warming in the study region has been estimated by the MEDCIE (2008), namely according to IPCC scenario A1B. The corresponding shift in elevation range (often used as a proxy for simulating the effect of climate change on temperature, e.g. Heegaard and Vandvik 2004) was computed for the whole scene and for each cell by using the present equivalence between altitude (50-meter DEM, from IGN, www.ign.fr) and temperature (assessed through interpolating within the 1 km grid of AURELHY meteorological database, Bénichou and Le Breton 1987). In the course of the simulation, the growth potential for each species on each cell is thus gradually modified according to the relationships estimating site index from site conditions, including altitude virtually decreasing with the simulated warming. The present altitude range of each species in the study region, including a correction according to local aspect (Ladier 2004), is used to add a simplified climate-induced component to the mortality submodel: annual mortality rate is increased to an arbitrary value of 0.05 (i.e., within the range of drought-correlated mortality rates indicated in Allen et al. (2010), for some Abies species within the Mediterranean area) only for the cells shifted out of the species present altitude range due to temperature increase.

Fig. 6

Simulated silver fir evolution from 2010 to 2050 and to 2090 on the northern slope of mount Ventoux (south-east of France). On these 3D-diagrams, Abies alba stand basal area is indicated by the shade level (white: no fir); other species are not shown. The larger white areas near the mountain crest represent stone fields where forest is unlikely to progress; they are not taken into account in the simulation. The whole scene covers ca. 8 km × 5 km

In this simulation, silver fir is increasing in basal area till 2050 while expanding towards both the bottom and the top of the hill, reaching a maximum extent. Then, its basal area declines and silver fir disappears at the lowermost part of its 2050 range, where it is located, from then on, out of its specific altitude range; meanwhile, it still slightly expands towards the top with the highest levels of basal area in 2090 situated at the uppermost part, as its range and best site index values are brought upward by the simulated warming.

3.3 Simulating the effects of management

In order to demonstrate and compare the role of specific growth traits on the behavior of pine and beech, the simulated management options shown in Fig. 7 were designed in a way to be completely neutral regarding species: none of them is directly favored by the simulated selective thinnings. The proportion of each species is computed as the percentage of number of stems greater than half the local top height, that is only individuals expected to be in the overstorey or likely to reach it later. In the highest part of the forest where beech was already established in favorable site conditions, its progression is strong whatever the management. In the south-west part of the forest, the altitude is less favorable to beech and its progression is to some extent hindered by a low basal area (and low canopy cover) management allowing to indirectly favor the development of pines.

Fig. 7

Simulated evolution of beech from year 2000 to 2090 in the eastern part of south-facing slope on mount Ventoux; the whole scene covers ca. 8 km × 6 km with altitude ranging from 1,000 m to 1,600 m a.s.l. in the north-west of the forest. On these maps, the amount of Fagus sylvatica is indicated by the shade level (white: no beech; pine species not shown). The simulated management aims to maintain stand basal area either between 30 and 48 m²/ha (High Basal Area option) or between 13 and 20 m²/ha (Low Basal Area option) with more thinnings. The square frame indicates the region where Low BA option favors the pine

In real conditions, the difference between pine and beech evolutions would probably be emphasized if foresters cut more pines in areas where beech regeneration is highly vigorous while removing some beech at lower elevation where its competition is not so strong and the pine (namely P. nigra, in this altitude range) is interesting both for wood production and because it is more likely to survive in a warmer or drier climate.

4 Discussion

4.1 Key model and simulator features

The importance of individual size is reflected by the results of growth and mortality analyses and by the structure of the corresponding submodels: height growth of seedlings or saplings is strongly reduced by asymmetric competition (Fig. 2) which, by design, tends to decrease with tree size, and mortality risk is highly depending on tree relative height (Fig. 3). Hence, a tree-level model is crucial for representing the structure of a forest mixing individuals which sizes can easily differ by a factor of 100, and for predicting its evolution. Except for a planted overstorey, tree layers (understorey, regeneration layer) often have indistinct vertical and horizontal limits and together constitute a composite arrangement of individuals (with different species in many cases).

Although the cohort approach brings a strong improvement, simulation using individual-based models remains time-consuming and memory-demanding. While several modeling approaches involving complex analyses and simplified assumptions have been designed in a view to overcome this drawback (Bugmann 2001), none of them, however, seemed fully satisfactory, especially in keeping the richness offered by the tree level. The spatial sampling proposed in this study appears to be an efficient alternative for reducing simulation cost on a large area (or for a high number of patches) while keeping the possibility of describing more precisely (with a complete grid of cells) some stands immersed in a forest landscape and interacting with the other stands through dispersal/recruitment flows. Unlike stratified samples simply combining the possible ranges of environmental conditions and stand types (Urban et al. 1999), this sampling is truly spatial and allows taking into account the consequences of dispersal processes. TreeMig (Lischke et al. 2006) is another compromise between structure and efficiency with detailed modeling of seed dispersal and recruitment, and spatially interacting cells (with a standard cell size of 1 km) allowing for simulation on scales as large as a region or subcontinent. TreeMig is based on a distribution approach (trees distributed in height classes in each cell) and therefore might be less convenient for forest management simulation.

The cell size (100 m²) used in our model is in the lowest part of the range of cell sizes in gap models (Bugmann 2001); it allows simulating small gaps as well as bigger ones whether these gaps result from cuttings or from natural disturbance. The simulation should not use cells of other sizes, nor groups of less than nine cells, because the fitted growth relationships use indices for one and for nine 100 m² cells. To use larger cells, these relationships should be fitted again, but they would surely not be able to represent the effects demonstrated at the 100 m² resolution as a coarser resolution would not render the very short (within one 100 m² cell) and short range (group of nine 100 m² cells) competition effects that have been put into light by the model presented here.

Two kinds of interaction between 100 m² cells are considered in the model: a close interaction as tree growth is depending on some competition indices taking into account the near neighborhood, that is the eight cells around the subject cell (preliminary analyses, not presented, have shown that enlarging neighborhood does not bring any significant improvement), and an interaction with cells on a much wider range through “recruitment flows”.

The spatial sampling strategy together with the cohort approach and the overall simplicity of the model makes it possible to simulate forest dynamics on large areas with an individual-based model without requiring any complex upscaling model (if it is in any way scientifically achievable) and then avoiding some difficulties (Fleming et al. 2002) and information losses that might be linked to this upscaling. Rietkerk et al. (2002) state that “processes and patterns can be better understood and described, if based on cross-scale observations or modeling”. We would add that knowing the characteristics of the basic fine-scale components (trees) when examining the system (tree populations and communities) at any scale (cell, group of cells, patch or stand, whole forest or landscape) is a most favorable situation for analyzing its functioning and behavior.

Finally, it is worth mentioning that the structure of the forest dynamics model and of the Capsis simulation module allow to replace any relationship by a more efficient or more suitable version without significant change in the simulator, provided that the new version complies with this parsimonious framework, namely the type of entities modeled, the way space, site conditions and forest are abstracted and the kind of processes involved.

4.2 Recruitment

The recruitment submodel fitted relies on observation of seedlings distribution and is able to operate at the forest level in a context where colonization involves source stands which are not in the immediate vicinity of target stands. It is homogeneous with respect to species, ensuring consistency and avoiding the risk to get differences in species behavior due to model content rather than to real ecological traits or reasons. Such a submodel overcomes one of the weaknesses of most traditional gap models where seeds of all species are assumed uniformly available, which can result in overestimation of migration rates and diversity (Price et al. 2001). For A. alba and F. sylvatica, this submodel could be improved in the short term through recent advances (e.g., Sagnard et al. 2007) and continued attempts to model some parts of the complex chain of reproduction and dispersal processes leading to recruitment. As for the three pine species, alternative models for dispersal/recruitment of these species within closed stands are still lacking as far as we know (recruitment and dispersal models for P. nigra and P. sylvestris investigated by Debain et al. (2007) relate to open areas).

4.3 Growth and competition

In gap models, some resources (light, sum of degree-day) are treated explicitly, as well as competition for light, allocation of growth, and even more detailed ecophysiological processes (photosynthesis, respiration) in some recent “physiology-based gap models” (Bugmann 2001). But physiological approaches require intensive data which are still lacking for a range of complex mixed stand structures in the study area. Also, while height and height growth are crucial in multilayered stands, allocating photosynthates rather to height growth or to diameter growth in varying conditions is still challenging (Le Roux et al. 2001). One major departure of our model from the gap model concept concerns the growth submodel which is directly tackling stem dimensions variations (bypassing the physiological processes involved) and focusing thoroughly on the assessment of the different competition components responsible (together with site, age, and species) for these variations. For this purpose, we used a set of indices readily available as they only require number and stem size of trees. By the way, these indices are quite relevant in the frame of forest management and easy to link to operational silvicultural guidelines. A similar approach is used in the landscape model FORMOSAIC (Liu and Ashton 1998).

Our approach proved to account for some important specific traits or behaviors. Data analysis and model calibration suggested that pine species regeneration is more sensitive to competition within the range of situations studied (dataset D3), i.e., canopies either complete or recently opened by regeneration cuttings, which is consistent with the well-known low tolerance of these species to shade. On the other hand, pioneer species like pines are more likely to expand in open areas. In fact, Fig. 2 indicates that P. nigra seedlings grow as fast or faster than beech or fir seedlings when basal area (NBA) is low and not composed by species casting deep shade (low values of SDLfa). Until now, the model is not intended for predicting the evolution of open areas which were therefore neutralized when simulating on a large area including such situations.

Among the indices we used, the sum of diameters of larger trees (SDL), split in two parts accounting respectively for deep shading (SDLfa) and for light shading species (SDLp), brought interesting results about competitive interactions. Like basal area of larger trees (Wykoff 1990), the sum of diameters of larger trees estimates asymmetric competition. This index was used by Dhôte (1991) in a quite different situation: even-aged beech stands. In our case study, SDL proved to be more efficient than the basal area of larger trees in explaining growth of understorey trees and seedlings (preliminary results, not presented). When using diameter, small trees or saplings of the understorey, closer in height and possibly inducing direct shading, get a weight which is not as low (as when using basal area) compared to high diameter trees of the canopy, which crowns, though larger, are farther and might let some indirect or lateral light reach the regeneration layer. SDL was split by separating pine species (SDLp) with lighter crown from shade tolerant Fagus and Abies (SDLfa) supposed to cast deeper shade. A similar separation was done by Yang et al. (2003) for mortality modeling.

We found that SDLp does not significantly affect the height growth of seedlings/saplings of any of the studied species, which is probably due to a rather high level of light transmission through pine crowns. But both SDLp and SDLfa have a significant effect on diameter growth, which is probably less sensitive to light resource but more sensitive to below-ground resources, namely water: in fact, these resources are used by pine competitors as well as by other species.

None of the indices used in this study depend on the accurate location of trees within cells. Tree and seedling locations are expensive to obtain, and we rather collected less intensive data on a wider range of situations (stand structure, composition, age and density, site conditions). But the major reason is that the benefit of distance-dependent and of some light indices is still questionable (e.g. Stadt et al. 2007) and can be outweighed by the cost of the spatial or crown-level information required. Moreover, the variability of individual growth is partly due to other biotic or abiotic factors which are not measured either.

4.4 Regular mortality

Tree mortality is a complex phenomenon as it can be induced slowly by competition, or rapidly by diseases, pests, climatic events (e.g., Crecente-Campo et al. 2009). In this study, the regular, non-catastrophic mortality recorded in the datasets is modeled considering that it is competition-induced, although the finishing stroke on an individual severely weakened by competition can be a dry year, a wind blow, a snow break or a biotic damage, with interactions that cannot be disentangled. This continuous mortality is modeled both for seedlings/saplings and trees through an approach linking mortality rate to local stand density (local competition intensity) and local relative size (individual competition status). Other common approaches link mortality to current growth (e.g., Kunstler et al. 2005; Wunder et al. 2007). Unless having a high frequency in sample monitoring, growth before death can hardly be investigated on seedlings when they quickly disappear after dying, which is often the case in mountainous conditions because of snow, water run-off, etc.; the 5-year interval between subsequent observations in this study is too long to obtain such data.

Another question still open in the frame of non-catastrophic mortality concerns tree longevity and a possible maximum age as a survival limit. In this study, this was addressed by attributing an increased probability of death to individuals with age reaching the maximum age observed in the National Forest Inventory (dataset D2; except for P. nigra, which has been planted in this area only from the end of the nineteenth century onwards), without trying to relate this limit to site conditions.

4.5 Present and future site conditions

Site conditions (altitude, aspect and climate, topography, geologic features and soil) influence basic processes at tree-level (growth, mortality) and thus partly determine forest composition and structure (e.g., Abrudan and Mather 1999). The model presented can use a rating of site conditions according to altitude, aspect and type of bedrock, drawn up for the five species studied, and on their present altitudinal range (Ladier 2004) in addition to site index curves available for even-aged stands of these species. Site index curves are based on the very ancient observation (as indicated by Batho and Garcia 2006) that the height reached by trees at a given age is an indicator of site conditions. This basic assumption was refined by considering the dominant height of the stand; the stand should be pure and even-aged for ensuring that the trees measured for evaluating this index have always been in dominant position. In our model, site index curves are used to estimate a potential height or a potential height growth for any tree, including seedlings and overtopped trees, in pure as well as in mixed stands. In mixed, uneven-aged stands, only a minority of trees can actually reach this value, but this potential is a useful reference to characterize a tree of a certain species, at a given age and in specified site conditions.

On the other hand, the climatic component of site conditions has been changing significantly in the recent past and cannot be considered as steady over the life span of a tree. Cailleret and Davi (2010) indicate that the optimal radial growth potential of the two colonizing species A. alba and F. sylvatica was shifted upward by the recent global warming in our study area. Concerning the simulation of climate-induced mortality, we assumed that the species will decline quickly when they get out of their altitude range as defined by their present limits, which seems reasonable for fir and beech on the mount Ventoux although the possible range of a species can be broader than its realized range (e.g., Keane et al. 2001). This allowed for a simulation of the outlines of silver fir (and beech, not shown) migration as the result of decline at lower altitude and of progression towards higher altitude through the dispersal/recruitment and the growth predicted by the model, within the limits of their moving altitude range.

Ongoing research on forest decline (e.g., Rachedi 2008, for silver fir in our study area) will help clarify and better estimate mortality thresholds by using some physiological approaches (as in one kind of the hybrid modeling approaches discussed by Fontes et al. 2010). As the soil water holding capacity can partly compensate for rain deficit, large scale probabilistic maps of soil characteristics (Emberger 2010; Piedallu et al. 2010) will likely increase simulation reliability. Indeed, a main key point to determining the consequences of climate change on tree growth is the ability to predict the evolution of site conditions (not simply using a shift in elevation as a proxy for climate warming) whatever the modeling approach (site index approach, or process-based approach) used subsequently for simulating tree growth and mortality, and forest evolution. Also, as adaptation and migration often occur together (Davis and Shaw 2001), the advances in the knowledge of local genetic adaptation will have to be taken into account when simulating the future range of the species.

4.6 Application to forest management

This model allows simulation of silvicultural practices at tree and stand scale, and of landscape changes at the stand and forest scale. Although presented separately here (Figs. 6 and 7), a migration induced by climate change and any management option can be simulated concurrently in order to elaborate silvicultural strategies aiming at reducing the risk of decline or at promoting migration of some species. The simulator still needs some development to allow choosing the “target species” according to site conditions and to stand location in relation to species range as moved by climate warming. Such a decision-making tool is likely to help foresters design anticipating management options (Coates and Burton 1997).

Some components of the model might still be improved while keeping the overall framework and a continuous evaluation process should be conducted in order to measure both its present efficiency and the gain from the improved submodels. Evaluation cannot rely only on any independent datasets (seldom available, anyway) but should also combine qualitative appraisal of model behavior, examination of process formulation, and sensitivity analyses (Soares et al. 1995; Bugmann 2001).

Initialization is a major issue in simulation at the forest or landscape scale. The quality of the prediction depends on the information input (stands, trees, site conditions), namely on the proportion of the real, directly usable, data and of the virtual, pre-simulated, ones. The pre-simulation aims at completing the data structure and most often involves some stochastic procedures. As a result, the “initial forest” described for starting the simulation itself, while being quite similar in overall to the real situation, can be different from it at any point of the scene; one should bear this in mind when considering simulation results, namely in the frame of forest management, or for comparison with real data. While the lack of spatially continuous data layers on large areas has long been considered as a limitation to landscape modeling (Baker 1989), new technologies (e.g., airborne LiDAR) are likely to partly fill the gap. However, field measurement will probably remain the only means to get reliable regeneration information (abundance, species, individual sizes); taking advantage of location devices, GIS and large scale forest structure description, new inventory designs and statistic sampling for regeneration are needed in a view to be workable on a whole landscape.

5 Conclusions

The association of simulation and modeling is especially fruitful when the aim and scope of simulations are clearly defined and taken into account from the model conception. Thus, the model and the simulation approach presented in this paper contribute to significant advances in representing the mechanics of a rather complex system, at stand and forest scales, allowing to extend and refine our understanding of how some patterns emerging at large scale are shaped by a few ecological processes, depending on site conditions diversity and evolution, and on competition and tree characteristics at fine scale. Though simple, the model appeared to properly represent the species behavior as far as necessary to reproduce stand and landscape evolution and patterns. This is in agreement with the assertion that predictive ability does not require model complexity (Didion et al. 2009).

While avoiding useless simulation in large parts of the whole area, thus reducing time and memory requirement, the spatial sampling approach used for simulating forest evolution allows retaining tree/seedling level. In fact, tree-level must be addressed both for ecological reasons related to knowledge of ecosystem functioning and for technical motivation as the individual tree remains the basic object handled by any silvicultural intervention.

Finally, while providing a framework for further improvement of the model and of the associated knowledge, this approach allows to simulate at reasonable cost the evolution of a forest as influenced by some leading natural trends and by conventional or innovative management options.

Appendix: Parameter estimates and regression statistics for each submodel and species

Hereafter, “R ²” denotes the coefficient of determination, “mean” is the mean of dependent variable, “CVR” is the residual coefficient of variation (standard deviation of the residuals divided by the mean of the dependent variable), in percent, and “N” is the sample size.

Recruitment submodel (see also Fig. 1)

$$ Log(smallsdgs + 1.1) = r1 + r2 \cdot {\left( {neardist + 15} \right)^{{ - r3}}} $$

(1a)

Abies alba

(R ² = 0.67, mean = 0.55/year, CVR = 105%, N = 59)

Parameter	Value	Standard error
r1	−0.51	0.56
r2	11.65	3.69
r3	0.40	0.15

Pinus nigra

(R ² = 0.13, mean = 0.39/year, CVR = 155%, N = 179)

Parameter	Value	Standard error
r1	−0.046	0.129
r2	5	fixed
r3	0.75	0.09

Pinus sylvestris

(R ² = 0.06, mean = 0.17/year, CVR = 137%, N = 214)

Parameter	Value	Standard error
r1	0.10	0.06
r2	1.04	2.13
r3	0.70	0.84

Pinus uncinata

(R ² = 0.05, mean = 0.21/year, CVR = 162%, N = 188)

Parameter	Value	Standard error
r1	0.10	fixed
r2	1.02	0.88
r3	0.56	0.28

$$ Log(smallsdgs + 1.1) = r1 \cdot {e^{{ - r2 \cdot {{\left( {neardist + 15} \right)}^{{r3}}}}}} $$

(1b)

Fagus sylvatica

(R ² = 0.27, mean = 1.02/year, CVR = 90%, N = 214)

Parameter	Value	Standard error
r1	4.17	4.28
r2	0.30	0.63
r3	0.34	0.29

Diameter Growth submodel for Trees

Pinus nigra, P. sylvestris, P. uncinata

$$ \begin{array}{*{20}{c}} {DINC5 = 0.1 \cdot a1 + \frac{{a2}}{{1 - a3 \cdot Log\left( {HdomINC5} \right)}}} \\ { \cdot \left( {1 - a4 \cdot {e^{{ - a5 \cdot NB{A^{{ - a6}}}}}}} \right)} \\ { \cdot \left( {1 + a7 \cdot ICS} \right)} \\ { \cdot min\left[ {1.025,\left( {1 + a8 \cdot \left( {Hdom50 - 15} \right)} \right)} \right]} \\ \end{array} $$

(2a)

(R ² = 0.88, mean = 0.36 cm/year, CVR = 22%, N = 654)

Parameter	Value	Standard error
a1	−45.98	11.68
a2	64.10	12.46
a3	0.0212	0.0093
a4	0.401	0.060
a5	19.72	6.42
a6	1.23	0.13
a7	0.232	0.054
a8	0.015	0.0033

Equation 2a was fitted only for P. nigra. As P. sylvestris and P. uncinata have many traits in common with P. nigra, we used a correspondence (set up using D2; not presented) between diameter increment of each three species, as depending on HdomINC5, NBA, and ICS.

Abies alba

$$ \begin{array}{*{20}{c}} {DINC5 = 100 \cdot b1 \cdot 1 /\left( {1 - b2 \cdot Log\left( {H/Age} \right)} \right) \cdot Ag{e^{{ - b3}}}} \\ { \cdot 1/\left( {1 - b4 \cdot {e^{{ - b5 \cdot NBA}}}} \right)} \\ { \cdot \left( {1 - {e^{{ - b6 \cdot ICS}}}} \right)} \\ \end{array} $$

(2b)

N.B. : In this case, ICS is computed using dbh under bark, and DINC5 is the increment under bark. The switch between under bark and over bark dbh values is achieved using bark thickness relationships (established using D2, not presented). NBA is over bark.

(R ² = 0.61, mean = 0.38 cm/year, CVR = 32%, N = 644)

Parameter	Value	Standard error
b1	0.122	0.061
b2	4.09	2.56
b3	0.151	0.053
b4	0.786	0.052
b5	0.014	0.007
b6	3.55	0.55

Fagus sylvatica

$$ \begin{array}{*{20}{c}} {DINC5 = 100 \cdot c1 \cdot {{\left( {1 - {e^{{ - c2 \cdot H/Age}}}} \right)}^{{c3}}} \cdot Ag{e^{{c4}}}} \\ { \cdot 1/\left( {1 + {e^{{c5 \cdot NBA}}}} \right)} \\ { \cdot \left( {1 - {e^{{ - c6 \cdot ICS}}}} \right)} \\ \end{array} $$

(2c)

N.B.: underbark and overbark values: see Eq. 2b

(R ² = 0.69, mean = 0.25 cm/year, CVR = 29%, N = 686)

Parameter	Value	Standard error
c1	1.28	0.18
c2	0.1	fixed
c3	0.86	0.05
c4	0.056	0.053
c5	0.020	0.002
c6	0.84	0.07

Diameter Growth submodel Seedlings and saplings

$$ \begin{array}{*{20}{c}} {DINC5 = 5 \cdot 0.1 \cdot d1 \cdot {e^{{ - {e^{{d2 - d3 \cdot Dlag}}}}}}} \\ { \cdot {e^{{ - \left( {d4 - {{H} \left/ {{H\max }} \right.}} \right)}}}^2} \\ { \cdot 1/[1 + d5 \cdot ((SDLp + 1)/100)]} \\ { \cdot 1/[1 + d6 \cdot {{((SDLfa + 1)/100)}^{{d7}}}]} \\ { \cdot 1/[1 + d8 \cdot {{((NBA + 0.01)/60)}^{{d9}}}]} \\ \end{array} $$

(3)

When parameter d4 is non-significant (“NS”), the modifier ${e^{{ - \left( {d4 - {{H} \left/ {{H\max }} \right.}} \right)}}}^2$is not included (i.e. it is replaced by 1).

Abies alba

(R ² = 0.83, mean = 0.27 cm/year, CVR = 40%, N = 864)

Parameter	Value	Standard error
d1	12.99	0.80
d2	0.1	fixed
d3	1.24	0.19
d4	-	NS
d5	−0.071	0.024
d6	0.72	0.08
d7	2.08	0.15
d8	1.38	0.16
d9	2	fixed

N.B.: d5 < 0 : i.e. SDLp is favorable to diameter growth of Abies alba seedlings/saplings.

Fagus sylvatica

(R ² = 0.51, mean = 0.25 cm/year, CVR = 47%, N = 1,065)

Parameter	Value	Standard error
d1	6.34	0.58
d2	0.25	0.07
d3	5.01	1.10
d4	-	NS
d5	0.048	0.029
d6	0.51	0.067
d7	1.59	0.20
d8	0.46	0.12
d9	2	fixed

Pinus nigra

(R ² = 0.45, mean = 0.46 cm/year, CVR = 35%, N = 448)

Parameter	Value	Standard error
d1	8.57	0.62
d2	−0.056	0.333
d3	14.28	4.97
d4	0.19	0.05
d5	0.40	0.07
d6	1.47	0.53
d7	1.65	0.44
d8	13.59	6.65
d9	8	fixed

D3 contains too few growth samples or site/stand conditions for Pinus sylvestris or P. uncinata; for simulation, we used the relationship fitted for Pinus nigra.

Height–Diameter submodel for Trees

$$ H = Hdom \cdot \frac{{(1 - {e^{{ - s \cdot DBH}}})}}{{(1 - {e^{{ - s \cdot Ddom}}})}} $$

(4)

With $s = s0 + s1 \cdot {\rm Re} lSpcgMod$ (see § 2.4.4)

Abies alba

(R ² = 0.97, mean = 17.34 m, CVR = 12%, N = 5,689)

Parameter	Value	Standard error
s1	0.000257	0.000024
s0	0.0130	0.00214

Fagus sylvatica

(R ² = 0.92, mean = 12.55 m, CVR = 14%, N = 16,081)

s1	0.000535	0.000016
s0	0.0300	0.00161

Pinus nigra

(R ² = 0.97, mean = 12.54 m, CVR = 11%, N = 9,846)

s1	0.000131	0.000015
s0	0.0659	0.00175

Pinus sylvestris

(R ² = 0.94, mean = 10.28 m, CVR = 14%, N = 24,397)

s1	0.000342	9.175E-6
s0	0.0276	0.00129

Pinus uncinata

(R ² = 0.93, mean = 9.76 m, CVR = 15%, N = 3,696)

s1	0.000089	0.000021
s0	0.0505	0.00304

Height Growth submodel for Seedlings and saplings (see also Fig. 2)

$$ \begin{array}{*{20}{c}} {HINC5/{\text{HdomINC5}} = } \hfill \\ {{1/[1} + {\text{t1}} \cdot {\text{Hla}}{{\text{g}}^{\text{t2}}} \cdot {\text{((1 - H / Hmax}}{{)}^{\text{t3}}} + {0}{.01)]}} \hfill \\ { \cdot {1/[1} + {\text{t4}} \cdot {{{\text{( (SDLfa}} + {1) /100)}}^{\text{t5}}}{]}} \hfill \\ { \cdot {1/[1} + {\text{t6}} \cdot {\text{(NBA}} + {0}{.1) /50}{{)}^{\text{t7}}}{]}} \hfill \\ \end{array} $$

(5)

When parameter t3 is non-significant (“NS”), the term ${\text{ ((1 - H / Hmax)}}{{ }^{\text{t3}}}{ } + { 0}{.01)}$ is not included (i.e. it is replaced by 1).

Abies alba

(R ² = 0.67, mean = 0.38 m/m, CVR = 65%, N = 1,173)

Parameter	Value	Standard error
t1	34.1	3.5
t2	8	fixed
t3	-	NS
t4	0.29	0.06
t5	2.18	0.30
t6	0.12	0.05
t7	0.50	0.39

Fagus sylvatica

(R ² = 0.30, mean = 0.75 m/m, CVR = 56%, N = 1,161)

Parameter	Value	Standard error
t1	6.98	1.39
t2	2	fixed
t3	5.32	1.23
t4	0.15	0.05
t5	1.67	0.35
t6	0.009	0.053
t7	2	fixed

Pinus nigra, Pinus sylvestris and P. uncinata

(R ² = 0.65, mean = 0.62 m/m, CVR = 33%, N = 760)

Parameter	Value	Standard error
t1	6.69	1.11
t2	2.43	0.30
t3	14.81	1.57
t4	1.39	0.46
t5	2.39	0.53
t6	5.06	2.56
t7	10.48	2.80

For Pinus sylvestris or P. uncinata, D3 contains too few individuals with height growth measurements and the corresponding range of site/stand conditions is quite narrow; for simulation, we used the relationship fitted for P. nigra.

Mortality submodels

Trees

$$ {\text{probMort5}} = {\text{(m1}} + {\text{m2}} \cdot N{\text{BA)}} \cdot {{\text{e}}^{{{\text{(m3}} + {\text{m4}} \cdot N{\text{BA)}} \cdot {\text{DBH/Ddom}}}}} $$

(6)

DBH and Ddom are under bark values. N is the number of combinations of classes of NBA (8 classes ranging between 10 and 60 m²/ha) by classes of DBH/Ddom (7 ranging between 0.3 and 1.1). The mean is the mortality rate for 5 years.

Abies alba

(R ² = 0.32, mean = 0.0186, CVR = 93%, N = 48)

Parameter	Value	Standard error
m1	−0.0088	0.035
m2	0.0022	0.0012
m3	−2.09	1.25
m4	0.0021	0.0178

Fagus sylvatica

(R ² = 0.40, mean = 0.0037, CVR = 72%, N = 56)

m1	−0.0014	0.0016
m2	0.00026	0.00008
m3	−0.33	0.70
m4	−0.011	0.015

Pinus nigra

(R ² = 0.55, mean = 0.0060, CVR = 65%, N = 49)

m1	0.0028	0.0089
m2	0.00086	0.00033
m3	−2.48	1.08
m4	−0.0012	0.0257

Pinus sylvestris

(R ² = 0.76, mean = 0.0186, CVR = 22%, N = 42)

m1	0.0066	0.0064
m2	0.0013	0.0003
m3	−0.48	0.32
m4	−0.028	0.009

Pinus uncinata

(R ² = 0.14, mean = 0.0195, CVR = 54%, N = 53)

m1	0.00051	0.00523
m2	0.00082	0.00026
m3	0.92	0.53
m4	−0.037	0.014

Seedlings/Saplings (see also Fig. 3)

$$ probMort5 = 5 \cdot \left[ {1 - {e^{{ - m0 \cdot NB{A^2} \cdot \left[ {1/min(l,H/H\max ) - 1} \right]}}}} \right] $$

(7)

N is the number of combinations of classes of NBA (from 10 to 60 with classes of width 10 m²/ha) by classes of H/Hmax (from 0.1 to 1.0 with classes of width 0.1). The mean is the mortality rate for 5 years.

Abies alba

(R ² = 0.85, mean = 0. 0175, CVR = 115%, N = 27)

Parameter	Value	Standard error
m0	2.803E-7	8.569E-8

Fagus sylvatica

(R ² = 0.30, mean = 0.0095, CVR = 210%, N = 29)

m0

4.607E-7

1.147E-7

Pinus nigra

(R ² = 0.56, mean = 0. 0395, CVR = 122%, N = 25)

m0

1.408E-6

3.186E-7

Pinus sylvestris

(R ² = 0.94, mean = 0. 0295, CVR = 51%, N = 11)

m0

1.007E-6

3.027E-7

Dataset D3 contains too few samples for Pinus uncinata; for simulation, we used the relationship fitted for Pinus nigra.