A forest simulation approach using weighted Voronoi diagrams . An application to Mediterranean fir Abies pinsapo Boiss stands

Aim of study: a) To present a new version of the forest simulator Vorest, an individual-based spatially explicit model that uses weighted Voronoi diagrams to simulate the natural dynamics of forest stands with closed canopies. b) To apply the model to the current dynamics of a Grazalema pinsapo stand to identify the nature of its competition regime and the stagnation risks it is currently facing. Area of study: Sierra del Pinar de Grazalema (S Spain) Material and methods: Two large plots representative of Grazalema pinsapo stands were used to fit and validate the model (plus 6 accesory plots to increase the availability of mortality data). Two inventories were carried out in 1998 and 2007 producing tree size and location data. We developed a forest simulator based on three submodels: growth, competition and mortality. The model was fitted, evaluated and validated for Grazalema plots. The simulation outputs were used to infer the expected evolution of structural diversity of forest stands. Main results: Vorest has proved to be a good tool for simulating dynamics of natural closed stands. The application to Grazalema pinsapo stands has allowed assessing the nature of the main processes that are driving its development pathway. We have found that the prevailing size-asymmetric competition dominates the self-thinning process in small-sized trees. At the same time, there is an active tree-size differentiation process. Research highlights: • Vorest has proved to be a good tool for simulating natural stands with closed canopies. • The Grazalema pinsapo stand under consideration is currently undergoing a natural process of differentiation, avoiding longterm stagnation.


Introduction
Abies pinsapo Boiss belongs to the Mediterranean group of firs that are the most ancient and less evolutionarily developed European species of this genus.Together with the Algerian species Abies numidica and the Turkish Abies cilicica, Abies pinsapo seems to be one of the first species to differentiate from the hypothesised common tertiary ancestor of all current Mediterranean species of this genus, including the widely spread Abies alba (Linares, 2011;Xiang et al., 2009). 2 as a result of ancient human intervention in the Mediterranean basin (Linares, 2011), have caused its current extremely low genetic diversity (Scaltsoyiannes et al., 1999) and its current status as endangered species.The small size of Abies pinsapo populations presents a considerable extinction risk.There are only three small areas of pinsapo forest in Spain, with a total extension of 3750 ha, from which less than 50% (1758 ha) are actually forests with closed canopies (MA-GRAMA, 2015).Grazalema is one of these areas with nearly 500 ha of pure dense pinsapo stands that are the object of this study.The consequences of climate change involving an increase of temperature and a decrease of rainfall are likely to further reduce the chance of survival of these unique stands (Linares et al., 2009).Conservation management of these ecosystems requires a thorough understanding of its current dynamics that were heavily influenced by past human disturbances and subsequent conservation attempts (Linares & Carreira, 2009).It is assumed that the regeneration dynamics of pure Abies stands mostly rely on disturbance events opening up potential growing space (Morin, 1994).Abies dynamics involve long periods of biomass accumulation and increasing stocking (Matic, 2001) followed by eventual disturbance events (mainly pests and diseases) that provide sporadic opportunities for new cohorts of trees (Antos & Parish, 2002).Gap processes and patch dynamics have also been shown to play an important role in shadetolerant stands that have been undisturbed for a long time, leading to a slower and more continuous turnover of the canopy (McCarthy, 2001;Antos & Parish, 2002).The dynamics of stands in intermediate overstocked periods are hardly understood.Competitive interaction and resulting differentiation and mortality processes are crucial issues for the dynamics and the stand structure at these developmental stages (Gray & He, 2009).Competition and also growth in dense forest stands are highly dependent on resource partitioning among trees.Trees acquire growth resources gathering them from their surrounding space and therefore resource partitioning has much to do with space partitioning.
Voronoi diagrams divide space into cells according to the closest distance to a tree location.For every tree location, its Voronoi region contains the points of the stand area that are closer to it than to any other tree location.The original Voronoi diagrams, developed by Georgy Voronoi in the early 1900s, consider Euclidean space.Since then, many generalizations of the Voronoi diagrams have been developed considering different spaces, different points or different metrics.The book by Okabe et al. (2000) provides a good overview of the subject.Our model considers tree locations and the corresponding Voronoi cells are weighted by the size of each tree.We use different distance functions for computing distances between different trees.This way, we can simulate the different ability of neighboring trees to explore the surrounding space and can ensure an adequate partitioning of space among them.
For forest stands with closed canopies, the more sizesymmetric the competition the higher the risk of stagnation (Oliver & Larson, 1996), especially if drought episodes become more frequent as it has been hypothesized for pinsapo stands in Spain (Linares, 2011).Over the last decades Spanish pinsapo stands have undergone a dramatic change in their dynamics due to a significant shift in management, i.e. from intensive timber management to strict conservation.Currently, after more than fifty years of almost total suppression of traditional ways of forest use including firewood, charcoal and livestock grazing, they have reached quite a high level of biomass accumulation (Coca et al., 2001;Linares et al., 2009) and there is no recruitment of small trees.On the other hand there is abundant advance regeneration in the understory waiting for opportunities (Arista, 1995;Cuadros et al., 2003).
Size-symmetric competition in forest stands has been defined as a situation where each individual competes and obtains growth resources in proportion to its size whilst size-asymmetric competition is defined as a situation where larger individuals have a disproportionate competitive effect on smaller individuals leading to the mortality of smaller ones (Weiner, 1990;Hara, 1992;Lundqvist, 1994;Bauer et al., 2004).
Size-asymmetric competition is often associated with even-aged stands and size-symmetric competition with uneven-aged ones (Lundqvist, 1994).In even-aged stands, canopy closure results in the extinction of light below canopy level and therefore in the suppression and death of slightly smaller individuals that remain below this light threshold (Weiner & Damgaard, 2006;Lundqvist, 1994;Oliver & Larson, 1996).However, in uneven-aged forests, the canopy is not totally closed at any height because otherwise there would not be enough light for smaller tree cohorts to develop.As a consequence self-thinning is not limited to small trees, rather it is distributed over the entire size range (Lundqvist, 1994).The nature of competition has much to do with stand structure and mostly with light extinction due to canopy closure.Nevertheless, size-asymmetric competition requires that some previous size differentiation among trees has established before canopy closure (Oliver & Larson, 1996).Otherwise, the competition among equal sized competitor trees would lead to an equal partitioning of resources (Uchmansky, 2000).In a forest stand with a closed canopy, this can lead to a high risk of stagnation, because many trees will approach light starvation as a result of insufficient 3 Vorest for natural dense forests: Pinsapo stand application Sierra del Pinar de Grazalema, located in the south of Spain (36° 46' N;5° 23' W).
Pinsapo stands typically grow on weathered limestone soils and, due to the exposition of these mountains facing the coast, there are rainfall averages of 2,133 mm/year, being the wettest location in the Iberian Peninsula.Nevertheless Mediterranean climate prevails with a distinct summer drought period (accumulated average rainfall in July and August is less than 10 mm) and high rainfall variability between years of approximately 500 mm per year.
Grazalema is, as a whole, an irregular forest with wide size ranges (in terms of diameter at breast height, or dbh, of trees), although the current high overstory density is effectively preventing the recruitment of new cohorts of small trees.
The main monitoring layout includes two 2100m 2sized rectangular plots, representative of the Grazalema pinsapo stand, one of them, plot F, was used to fit and evaluate the model and the other one (plot V) to validate it.Additionally we analyzed mortality data from another six similar plots (M1 to M6) in the same stand.All plots were surveyed in 1998 and again in 2007 so that the data span a period of 9 years (Coca et al., 2001;Cuadros et al., 2003;Abellanas et al., 2005;Cuadros et al., 2005).Table 1 shows the main stand characteristics of experimental plots.Among the variables measured in both surveys, we studied dbh, tree locations and mortality.The usefulness of this stand for modelling lies in the absence of silvicultural interventions during the last decades, i.e. the only type of mortality that the stand has undergone is natural mortality.
The increase in mean and maximum stand dbh values over nine years reflects the rarity of in-growth and the prevalence of mortality by suppression over ageing or external factors.
growing space due to a lack of self-thinning (Oliver & Larson, 1996).This is the risk that threatens most artificially established dense even-aged pure stands.Silvicultural thinning is a way to avoid stagnation in such cases, just doing what nature is not able to do because of extreme homogeneity.
The rapid restocking that Grazalema pinsapo stands have undergone in the last decades and the resulting canopy closure achieved poses the question about the nature of the dynamics that they are currently experiencing and their needs for silvicultural treatments.Size differentiation and self-thinning processes may be most instrumental in maintaining the stand in a healthy good condition.An absence of these processes would require silvicultural interventions in order to avoid stagnation and decline.
The objective of this paper is both to present an improved version of the growth model Vorest (Abellanas et al., 2007;Abellanas et al., 2012) that allows to deterministically simulate the natural dynamics of forests stands with closed canopies and to apply it for diagnosing the nature of the dynamics that the studied Grazalema pinsapo stand is currently undergoing.Our aim is to identify the nature of the prevailing competition regime and the level of stagnation risk this Grazalema pinsapo stand it is potentially confronted with.Our findings inform what silvicultural management is required for conserving these unique stands.

Study site and experimental layout
The study site is a 460 ha natural stand of Abies pinsapo, which occupies the steep slopes (average gradient of 50%) of the northern side of the range Tree growth is considered in the model in terms of diameter increment (id).The model relates tree growth to its capacity to acquire new growing space and estimates it by the area available to access solar radiation, assuming that light is the most limited resource in this stand and therefore the one that dominates the competition among tree neighbors.
When the main forest canopy is completely closed, the growing space is shared between the neighboring trees and typically most trees occupy a smaller growing space than they would do if they grew in the open.This behaviour leads to a corresponding reduction in tree growth compared to an open-grown tree of the same size (dbh).
The area that an open grown tree of the same dbh can explore for light we refer to as potential growing space (PGS) (or zone of influence) of a tree and we estimate this area from the crown projection area of an open grown tree with the same dbh.Some studies on crown size and diameter growth support the notion that growth follows the acquisition of canopy space (Thorpe et al., 2010).
We define the available growing space (AGS) of a tree as the area in its surroundings not occupied by neighboring trees and therefore available to that tree to search for light.Finally, we define occupied growing space of a tree (OGS) as the portion of its available growing space (as a subset of its potential growing space) that the tree actually occupies in terms of crown biomass.
VOREST estimates the available growing space (AGS) of each tree by the area of the multiplicatively weighted Voronoi region assigned to it in the Section Competition model (Römisch, 1996;Sterba & Zingg, 2001) and computes the occupied growing space (OGS) by limiting the extent of the Voronoi region (AGS) by the range of its potential growing space (PGS) or zone of influence (radius in short) according to its size (dbh).
The Voronoi diagram subdivides space acording to the proximity to a given set of points (Okabe et al., 2000).Our points of interest are tree locations and the Voronoi regions correspond to the space closer to each of them rather than to another tree location.For taking into account size differences between trees, the model uses weighted Voronoi diagrams (Aurenhammer & Edelsbrunner, 1984).In multiplicatively weighted Voronoi diagrams, Euclidean distances are modified by a factor that depends on the tree to which a distance is measured.In our case, the larger the tree, the smaller the factor.In this way, the Voronoi partitioning simulates the fact that larger trees have access to more resources.As

VOREST model description
The model description follows a reduced version of the ODD (overview, design concepts and details) standard protocol for describing individual-and agent-based models proposed by Grimm et al. (2006Grimm et al. ( , 2010)).

Purpose
The purpose of the model is to simulate the current space-time development of natural forest stands with closed canopies and no significant, active recruitment.The main driving forces are growth, competition, and mortality.

Entities, state variables and scales
The model handles two hierarchical levels: The lower level consists of individual trees and the upper level of forest stands (plots).The basic entities are individual trees that are characterized by the attributes and state variables specified in Table 2.The time resolution is one year (one time step of the simulation represents one year), according to the growth rate of trees in temperate zones.The model time horizon is between one and a few decades.
The spatial resolution is one decimeter.The model landscape corresponds to forest plots or stands between 1,000 m 2 and a few hectares.Tree locations are measured in meters (to the nearest decimeter) and dbh of trees in centimeters (to the nearest half millimeter).

Vorest overview
VOREST is a simulation model implemented as a computer program in C++ that consists of three submodels, which correspond to the three basic processes The model simulates the growth of trees and their resulting size, survival or mortality and the space they occupy.A number of emergent properties, some highly relevant to forest management, are derived from the simulation, such as the spatial structure of the stand, both in terms of tree locations and dimensional aspects, as well as of population characteristics such as size distribution or stocking and the stand self-thinning dynamics.

c) Interaction
The main interaction between trees is competition.The partitioning of growing space between trees is simulated by linking resource availability with physical space.The area actually assigned to a tree can increase or decrease during the simulation, consistently varying its growing rate.

d) Stochasticity
Currently, there is no stochasticity in the model.All simulated processes are deterministic.

Initialization
The input data for the initialization are provided in a text file (ASCII).The file includes the coordinates, species and initial diameter at breast height (dbh) of every tree.The header rows of this file also include the initial number of trees and function parameters (growth, mortality and allometric functions).These data are used in each run of the simulation.The length of the simulation period is optional and can be specified by the user during the execution of the program.

Growth model
The growth model is based on the potential modifier method (Newnham, 1964;Botkin et al., 1972;Pretzsch, 2009;Pommerening et al., 2011).Potential diameter increment (dbh) curves have been derived from dominant trees of the stands.The growth model used to fit potential diameter increment functions is the first derivative of the Chapman-Richards function (Pretzsch et al., 2002;Pommerening et al., 2011;Pommerening & Särkkä, 2013): ( ) [1]   potential growing space (PGS) or zone of influence, (maximum spatial range that a tree can explore depending on its size), we consider the area that a tree of the same size would occupy when growing in the absence of competition.The weights used to construct the weighted Voronoi diagrams rely on the relative size of the trees.The zone of influence (PGS) corresponds with the crown projection area of an open-grown tree, which is estimated on the basis of its dbh using allometric models from the literature (Krajiceck et al., 1961, Hasenauer, 1997;Ek, 1974;Bella, 1967, Farr et al., 1989, Leech, 1984, Smith et al. 1992;Paine & Hann, 1982).We call relative area (RelArea) the ratio of the occupied growing space (OGS) of a tree and its potential growing space (PGS) and it is used as a surrogate for the growing capacity of a tree.

Vorest schedule
In each simulation step, the following processes are run in the following order: 1. Synchronous size (dbh) updating of all trees according to the defined growing function.2. Identification of natural neighbors of each tree and updating the weights of each tree.3. Computing the range of the potential growing space (PGS) of each tree according to its current size.4. Synchronous simulation of natural mortality.5. Calculation of the occupied growing space (OGS) for each tree.

Design concepts a) Basic principles
The basic principles underpinning the model can be summarized as follows: The growth of each tree within the stand is calculated by applying a reduction coefficient to its potential growth to incorporate the effect of competition (Newnham, 1964;Botkin et al., 1972;Pretzsch, 2009, Pommerening et al., 2011).The potential growth of a tree (in analogy to that of an open-grown tree) is estimated by the growth of dominant trees of the stand (Pretzsch et al., 2002;Pommerening et al., 2011).Competition is simulated by the reduction in the occupied growing space (OGS) compared to the potential one (PGS), an effect caused by its neighbors.
Natural mortality occurs after a period of reduced or zero growth, which is related to the size of the tree.6 trees as an isotropic process.Evidence of directional growth factors could not be confirmed in this study.The competition load of a tree i at time step t is incorporated in the model as the ratio of OGS and PGS: OGS i,t is the occupied growing space of tree i in year t, computed as the area of the weighted Voronoi region of the tree i restricted by the range of its zone of influence (radius) at time t and PGS i,t is the potential growing space of tree i in year t estimated as the crown projection area of an open grown tree of the same dbh.
If all neighbors of a tree are sufficiently far from a target tree (i.e. when the zones of influence (PGS) of every neighbor do not overlap with that of the target tree), the occupied growing space (OGS) of the tree equals the potential one (PGS) and thus relative area (RelArea) equals one, indicating an open-grown tree.Vorest uses relative area (RelArea) [3] to define the final growth function in equation [4]: id i,t is the diameter increment of tree i at time step t and ν is a specific local model parameter.

Mortality model
Natural mortality can be either the consequence of sustained suppression or of ageing.In any case it is the ultimate result of poor growing conditions (Bauer et al., 2004).Vorest implements a memory function for simulating the death of a tree depending on its relative increment of the previous five years (Berger et al., 2004;Murphy & Pommerening, 2010;Pommerening Where id i,t pot is the potential diameter increment of tree i at time t; DBH i,t is the diameter at breast height of tree i at time t and A, k, p are model parameters.
Actual increment is simulated as a fraction of this potential increment depending on the competitive status of the tree.

Competition model
Competition is related to the magnitude of the reduction in growing space that a tree undergoes due to the presence of adjacent trees that act as competitors.Vorest considers the "natural neighbors" as competitors of a subject tree as defined by the weighted Voronoi diagrams (see Figure 1), unless the radii of their potential growing space (PGS) do not reach that of the target tree.
As a surrogate for potential growing space (PGS), Vorest incorporates functions published in the literature (Hasenauer, 1997) to estimate the radius of the crown projections of open grown trees as a function of dbh for several species.The occupied growing space (OGS) of each tree in each time period is simulated by the weighted Voronoi region that corresponds to it (AGS), limited by the range of its potential growing space, PGS.The weights used to construct the Voronoi diagrams are dependent on the relative sizes of the subject tree and its neighbors as shown in equation [2].
where d 1,t and d 2,t are the arithmetic mean of the stem diameters of the nearest neighbors and the arithmetic mean stand diameter, respectively.The weight w i,t of a tree i at time t, is dependent on the relative size of the tree i (DBH i,t ) in relation to the mean dbh of its neighbors at time t (local component) and on the dbh of tree i relative to the mean stand dbh at time t (global component).c and g are scaling parameters satisfying c + g = 1, c >> g (we used c = 0.8 and g = 0.2).
The individual distance function assigned to each tree is continuously adapted at each time step to generate the weighted Voronoi diagrams.This individual distance function is derived from the Euclidean distance divided by the weight assigned to the tree (w i,t ) according to equation [2].Thus, there is a positive relationship between the extent of the area of each tree estimated as AGS and its relative size in relation to its neighbors.
At each time step the weighted Voronoi diagrams define two issues: Which are the nearest neighbors of each tree (Figure 1) and the relative load of competition among them.The use of individual distance functions allows simulating the competition signals sent out by pd i,t 5 is the relative dbh increment of tree i during the five year period previous to year t; DBH i,t is the dbh of tree i in year t and DBH i,t-5 denotes the dbh of tree i in year t -5.Tree i dies in year t if α and β are model parameters.
Figure 2 illustrates the simulation process of this version of Vorest.

Graphical output
Apart from various estimations, Vorest also provides a graphical simulation output, which shows the spatial distribution of living and dead trees and the growing space occupied and available in the stand.We addressed the spatial edge effect by periodic edge correction methods (Illian et al., 2008) constructing a simulated edge (referred to as plus-sampling, see Figure 3).

b) Competition model
To fit the competition model we needed to estimate initial values of RelArea and according to its definition in [3] also of PGS and OGS.To estimate crown radius we used the parameter values proposed by Hasenauer (1997) for fir (Table 3).
To fit the final growth model [4] an auxiliary iteration of Vorest was run in plot F with the sole purpose of computing the initial tree AGS (i.e. the weighted Voronoi regions) and PGS areas that allowed obtaining initial values of OGS, PGS and RelArea for all trees but without simulating any growth.With these initial values, the growth model was fitted using nonlinear regression and least-square methods as implemented in the R software.

c) Mortality model
As natural mortality is a sporadic event in forest stands and our time series only spanned nine years, we added mortality data of plots M1 to M6 to those of plot F in order to increase the availability of mortality data.
A threshold value of relative dbh increment of the previous five years (pd 5 ) typical of live trees of the species Abies pinsapo has been investigated (Berger et al., 2004;Pommerening et al., 2011;Pommerening & Särkkä, 2013).Figure 5 shows the density functions of pd 5 for dead and live trees.
We can see that pd 5 does not allow to sharply split dead and live trees, even if the average differences are statistically significant (mean pd 5 = 0.017 for dead trees and 0.033 for live trees; p = 0.0069).This is likely so because firs are very shade tolerant species that can survive a very long time in the shade of an overstory.The fate of these suppressed trees depends on the creation of canopy openings around them, which would allow them to resume active growth.In a previous study of the same plots (Cuadros et al., 2005), a random sample of 32 live trees (trees with more than 1.5 m height and less than 5 cm dbh) were    The predictive ability of the growth model has been evaluated by statistical estimates of bias (mean error, ME) [9] and precision (mean absolute error, MAE) [10]:

Application to a Grazalema Abies pinsapo stand
DBH i and DBH i  are the observed and the estimated dbh values of tree i, respectively, at the end of each simulation period.DBH is the mean dbh of the plot at the same time.N is the number of live trees and p is the number of parameters of the model.The growth model is also assessed graphically by means of plotting expected versus observed dbh.The same assessment procedure is used for the evaluation and validation.The data of the plot used to fit the model were considered for evaluation (plot F) and data from an external plot for validation (plot V).
To evaluate and validate the mortality model two similar statistics were used: Mean deviation (MD) [11] and mean absolute deviation (MAD) [12] for each simulated period (Nunes et al., 2011): cut and dated by tree-ring counting at the stem base.The average age in this stand is 45 years indicating a mean dbh increment of 1.1 mm/year.We estimated a maximum value of pd 5 of these suppressed live trees of 0.022.We analysed the relationship between pd 5 and dbh for live and dead trees.Figure 6 shows that most dead trees are below 20 cm dbh and that the lower the dbh the higher is the pd 5 value of dead trees.We realized that a fixed value of pd 5 as mortality threshold would not fit well to this condition and consequently defined a decreasing threshold survival function for pd 5 as in [6] producing lower pd 5 mortality values for trees with larger dbh values.
We have fitted α and β parameters of this function using again the statistical package quantreg of R based on a quantile of 0.05 of the data of live trees (Figure 6).

Model evaluation and validation
The model performance was analyzed comparing simulated versus actual data both for growth and mortality.
To evaluate and validate growth simulation, goodness of fit has been assessed with the statistics root mean squared error [7] and fit index [8]:

Model parameters
Table 3 shows the fitted values obtained for the model parameters.

Graphical output
Figure 3 provides a graphical output of the initial and the final nine-year simulation of plot F.

Model evaluation
The values of the evaluation statistics for the dbh estimations are shown in Table 4. Figure 7 shows graphically the goodness of fit for the period of time simulated (9 years) and the actual and simulated dbh distributions for plot F. For the mortality evaluation we obtained values of MD = -0.01 and MAD = 0.21.The number of surviving trees at the end of the simulation (nine years) was 80 while the actual number of survivors in 2007 in plot F was 79. y i is a dummy variable that takes the value 0 if tree i actually dies during the simulated period and takes the value 1 otherwise.p i is another dummy variable that accounts for the simulated death or survival of trees.Comparing the number of dead trees in the simulation with the actual mortality in the same period is another global assessment of the mortality model.

Simulation of pinsapo stand dynamics and structure
In order to answer the initial question about the current stagnation risk of Grazalema pinsapo stands, we have projected the development of both plots, F and V, with Vorest 60 years forward.The aim of this simulation was to analyse the current trend of forest structure rather than to project the stand very far into the future.We are particularly interested in size differentiation, since differentiation is a natural strategy that dense forest stands adopt to avoid stagnation (Oliver & Larson, 1996).
From the simulated stands we have derived dbh and basal area distributions for consecutive simulation periods and then graphically analyzed the expected stand structure evolution.The expected evolution of stand structural diversity has also been assessed computing the Shannon diversity index (H') [13] for 10cm-diameter classes and for the same simulation periods.
p i is the proportion of trees belonging to the ith 10 cmdiameter class and R is the number of classes.11 Vorest for natural dense forests: Pinsapo stand application progressive expansion of the frequencies along an increasing interval towards the upper classes.The initial distributions reflect the presence of two distinct cohorts with a mean dbh of about 40 cm, the oldest and less numerous one, and around 15 cm, the youngest and most frequent one.The latter most likely derived from the rapid restocking which was the result of putting the stands under strict protection during the last decades.There are some big trees with dbh values greater than 50 cm but their frequencies are quite low.
In the basal area distributions (Figure 9, bottom), we can see a trend towards an increasing relative contribution of progressively larger diameter classes to the whole basal area of the stand.The expected evolution of size differentiation has also been assessed with the Shannon index applied to 10cm-diameter classes (Figure 10).The simulation indicates a progressive increase in the Shannon index value involving a subsequent increase in size diversity with time.This tendency points towards an evident differentiation process in the size distribution.

Discussion
Vorest has proved a good tool to simulate current Grazalema pinsapo dynamics through an appropriate simulation of self-thinning of dominated trees and by continuously re-adjusting the growing capacity of both, previously dominant and newly promoted trees.The growing space allocation among trees using weighted Voronoi diagrams has proved to be a good approach to simulate the mainly size-

Model validation
The data of another plot (plot V), not used to fit any of the submodels, was included in the validation of the model.The plot consisted of 127 trees at the beginning of the simulation in 1998.Table 5 presents the values of growth validation statistics.Figure 8 shows the goodness of fit for plot V and the period of time simulated (9 years) along with the corresponding observed and simulated dbh distributions.
For the validation of the mortality model we obtained values of MD = -0.04 and MAD = 0.15.The number of surviving trees at the end of simulation period (9 years) was 108 while the actual number of survivors in 2007 in plot V was 103.

Forecasting stand structure in Grazalema pinsapo stands
Figure 9 shows the expected evolution of diameter distributions and relative basal area distributions by diameter classes (for plots F and V).
The observed trend of the diameter distributions is similar in two ways: There is a marked decrease in the initial peak located around 15 cm, and there is a 13 Vorest for natural dense forests: Pinsapo stand application The presence of one older cohort in the Grazalema pinsapo stand, has probably made the increasing size evenness possible, blurring the known effect towards bimodality that is the typical consequence of prevailing size-asymmetric competition (Weiner, 1986;Bauer et al., 2004;Eichhorn, 2010).On the other hand the simulation has suggested a likely decline in mortality rates for the next decades (not shown) that together with the reduction in basal area and number of trees per hectare could indicate a global decrease in competition load.
Pinsapo stands of Grazalema are currently undergoing a process of mortality by suppression caused by a relatively sudden restocking based on the radical change of management imposed a few decades ago by the establishment of strict protection rules.However, our simulation has also indicated a sufficient promotion of some of the youngest cohort trees that progressively colonise the upper diameter classes, thus changing stand biomass distribution towards a more uniform partitioning among the size classes and a simultaneous increase of the weight of intermediate and larger trees.This trend allows expecting an increasing stability of the stand in the near future even in the absence of exogenous disturbances.

Conclusions
Using the present version of Vorest we have found that weighted Voronoi diagrams are a good method of modelling resource partitioning among trees and provide correct feedback loops steering the growthcompetition processes in conditions of high stand density.
The simulation with Vorest has shown that there is enough asymmetry in competition relationships among pinsapo trees to promote an active self-thinning.But due to the high shade tolerance of firs and likely due to the presence of only few interspersed big trees of previous cohorts, these dynamics have also entered an increasing size differentiation process.
The slow suppression process is typical of shade tolerant species and in our case a comparatively low mortality-growth threshold has been confirmed.
Our simulation suggests a likely decline in competition pressure that probably can allow for new episodes of recruitment to occur, even if no catastrophic regeneration episodes take place.
Our study also revealed that there is currently a low risk of stagnation in the pinsapo stands of Grazalema, due to an increase in structural diversity that tends to prosper because of self-thinning and the promotion of younger cohort trees by old seed trees.asymmetric competition that dominates differentiation and self-thinning in the Grazalema pinsapo stands.
Nevertheless, the version of Vorest presented in this paper deals with the specific dynamic stage that the Grazalema pinsapo stand is currently undergoing based on its history and is not supposed to cover extensive simulation periods.As stated by Antos & Parish (2002), in this kind of forests, variation in dynamics and in the processes controlling changes can be large in the long term.
The model presented here has considered a stand development stage involving a closed canopy with a high level of competition that is currently blocking new recruitment.The aim of the simulated 60 years of stand development was to elucidate the level of stability we can expect in the stand in the near future without the intervention of silvicultural treatments or other kind of disturbances.
The long term dynamics of fir forests are currently much debated in many countries.It is generally assumed that the shade tolerance and the branching pattern provide the genus with a high competitive advantage.Abies dynamics involve long periods of biomass accumulation and increasing stocking (Matic, 2001).Both dominated stems and seedlings, grow comparatively slowly and can persist for very long periods until their conditions improve (Antos et al., 2000;Antos & Parish, 2002).Consistent with this statement, we have found very low mortality related growth-rate threshold values for pinsapo trees.
Seedling banks play a major role in the regeneration process (Morin & Laprise 1997;Stewart et al., 1991;Woods, 1984;Antos & Parish, 2002).The recruitment of advance regeneration can occur as a result of major disturbances affecting canopy cover as well as of less catastrophic events leading to patch dynamics and gap processes or to even almost subtle endogenous events affecting only one or a few trees (Antos & Parish, 2002).
We have shown that, in our case, high levels of stocking together with high frequencies in lower diameter classes do not cause stagnation and the dynamics lead to an increasing size diversity.In an old growth spruce-fir stand, Antos & Parish (2002) have also found that there is no need for major disturbances to maintain a long-term stability in this kind of forest.In their case study, it probably was due to the low density of the subalpine forest that there was no sign of stagnation features after a long time period without exogenous disturbing events.In our case and for a limited time scale, the observed trend is probably the result of the presence of a few larger trees acting as seed trees.
As other authors (e.g.Antos & Parish, 2002) have shown, the absence of disturbances not necessarily leads to a stagnation of stands with closed canopies, if size differentiation is encouraged.
It has frequently been demonstrated that tree growth in stands with closed canopies is largely a function of local neighborhood competition and that the nature of this competition is mainly size-asymmetric (Oheimb et al., 2011;Grabarnik & Särkkä, 2011).The nature of competition and not only its intensity has proved an important factor driving forest dynamics and the resulting structure (Toda et al., 2010).As Bauer et al. (2004) have shown, it seems that size-asymmetric competition is indeed a general outcome of tree neighborhood interactions.

Figure 1 .
Figure 1.Natural neighbors (j x ) of a tree i at time t according to the weighted Voronoi diagram approach.

Figure 3 .
Figure 3. Graphical output of Vorest.Left: Plot F before the start of the simulation, right: Plot F after nine years of simulation.Blue dots are live trees and red dots are dead trees.Green bound regions are occupied growing spaces (OGS) assigned to live trees: Darker green areas for plot trees and pale green areas for constructed edge trees.Interspersed grey regions represent unoccupied growing space.
plot F to fit the potential growth model.From measured dbh values (1998From measured dbh values ( , 2007)), mean annual increment values have been derived for each tree by linear interpolation.The potential growth function [1] was fitted to the largest values of dbh increment by nonlinear regression using the statistical package quantreg of R (R Development Core Team, 2011), which allows fitting a function to a specified quantile of a set of data points.In our case, it was adjusted to 0.99 quantile of data points (Figure4).The regression yielded estimations of model parameters A, k and p [1].

Figure 4 .
Figure 4. Annual diameter increments of plot F computed by linear interpolation from dbh values in 1998 and 2007.The envelope curve represents the potential dbh increment, which has been obtained by fitting the first derivative of Chapman-Richards function to percentile 99 of the growth data fitted by the quantreg R package.

Figure 5 .Figure 6 .
Figure 5. Density functions for relative dbh increments of the previous 5 years (pd 5 ) for live (green line) and dead (red line) trees.

Figure 7 .
Figure 7. Growth model evaluation (plot F) involving a nine-year simulation period.Observed vs expected dbh values of 2007 (left) and actual and simulated dbh distributions (right).

Table 1 .
Stand characteristics of the monitoring plots in the Sierra del Pinar de Grazalema.

Table 2 .
Entities, attributes and state variables of the Vorest model.

Table 3 .
Vorest model parameters for Abies pinsapo in Grazalema (numbers in square brackets refer to the equations to which the parameters apply.t is the percentile used to fit the potential growth function; see Figure4).SE -standard error.(*)Hasenauer,1997(valuesfor Abies alba).Vorest for natural dense forests: Pinsapo stand application

Table 4 .
Growth model evaluation (plot F).The statistics refer to dbh estimates.Definition of statistics explained in the text.