Research Article


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


Begoña Abellanas

Departamento de Ingeniería Forestal. Universidad de Córdoba. Campus de Rabanales. Córdoba. Spain

Manuel Abellanas

Dpto. Matemática Aplicada. Facultad de Informática. Universidad Politécnica de Madrid. Campus de Montegancedo s/n. Boadilla del Monte. Spain

Arne Pommerening

Swedish University of Agricultural Sciences SLU. Faculty of Forest Sciences. Department of Forest Resource Management. Skogsmarksgränd. Umeå. Sweden

Dolores Lodares

Dpto. Matemática Aplicada. Facultad de Informática. Universidad Politécnica de Madrid. Campus de Montegancedo s/n. Boadilla del Monte. Spain

Simón Cuadros

Departamento de Ingeniería Forestal. Universidad de Córdoba. Campus de Rabanales. Córdoba. Spain



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:

  1. Vorest has proved to be a good tool for simulating natural stands with closed canopies.
  2. The Grazalema pinsapo stand under consideration is currently undergoing a natural process of differentiation, avoiding long-term stagnation.

Keywords: Vorest; stand dynamics; individual-based forest model; spatially explicit forest model; pinsapo.

Citation: Abellanas, B., Abellanas, M., Pommerening, A., Lodares, D., Cuadros, S. (2016). A forest simulation approach using weighted Voronoi diagrams. An application to Mediterranean fir Abies pinsapo Boiss stands. Forest Systems, Volume 25, Issue 2, e062.

Received: 11 May 2015. Accepted: 23 May 2016

Copyright © 2016 INIA. This is an open access article distributed under the terms of the Creative Commons Attribution-Non Commercial (by-nc) Spain 3.0 Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Funding: Spanish National Research Project from the Ministry of Science and Education: MTM2008-05043/MTM. Research contract from Andalusian Government (Junta de Andalucía. Consejería de Medio Ambiente).

Competing interests: The authors have declared that no competing interests exist.

Correspondence should be addressed to Begoña Abellanas:





Materials and methods






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). Some authors such as Alba-Sánchez et al., (2010) and Linares (2011) suggest an ancient segregation of the two Iberian Abies species: A. alba and A. pinsapo that could have remained isolated throughout the whole Quaternary. Abies pinsapo has occurred since then in two reduced areas on both sides of the Strait of Gibraltar, i.e. on European and African shores, between 1000 and 1800 m a.s.l. on the northern side, and between 1500 and 2150 m on the southern one. The genetic isolation and the limited area to which these populations have spread, aggravated by the intense habitat fragmentation 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 (MAGRAMA, 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 shade-tolerant 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 size-symmetric 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 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.

Materials and methodsTop

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 ­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 2100m2-sized 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.

Table 1. Stand characteristics of the monitoring plots in the Sierra del Pinar de Grazalema. N98, N07: Number of trees per ha in 1998 and 2007, respectively; G98, G07: Basal area (m2/ha) in 1998 and 2007, respectively; Dg98, Dg07: Mean quadratic stand diameter at breast height (cm) in 1998 and 2007, respectively; H98, H07: Mean stand height (m) in 1998 and 2007. St. dev. is standard deviation.

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.

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. (2006, 2010).


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.

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

The spatial resolution is one decimeter. The model landscape corresponds to forest plots or stands between 1,000 m2 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).

Process overview and scheduling

Vorest overview

Vorest is a simulation model implemented as a computer program in C++ that consists of three sub-models, which correspond to the three basic processes governing the main dynamics of the stands under consideration: growth, competition and mortality.

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 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.

b) Emergence

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.


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.

Sub models: Description

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):

Where is the potential diameter increment of tree i at time t; DBHi,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.

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

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 and are the arithmetic mean of the stem diameters of the nearest neighbors and the arithmetic mean stand diameter, respectively. The weight wi,t of a tree i at time t, is dependent on the relative size of the tree i (DBHi,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 (wi,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 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:

OGSi,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 PGSi,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]:

idi,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 et al., 2011, Pommerening & Särkkä, 2013) in relation to its size (dbh) according to equations [5] and [6].

is the relative dbh increment of tree i during the five year period previous to year t; DBHi,t is the dbh of tree i in year t and DBHi,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.

Figure 2. Vorest flowchart for natural stands with closed canopies.

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).

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.

Application to a Grazalema Abies pinsapo stand

Model fitting

a) Growth model

We used Grazalema plot F to fit the potential growth model. From measured dbh values (1998, 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 (Figure 4). The regression yielded estimations of model parameters A, k and p [1].

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.

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).

Table 3. Vorest model parameters for Abies pinsapo in Grazalema (numbers in square brackets refer to the equations to which the parameters apply. τ is the percentile used to fit the potential growth function; see Figure 4). SE – standard error.

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 (pd5) 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 pd5 for dead and live trees.

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

We can see that pd5 does not allow to sharply split dead and live trees, even if the average differences are statistically significant (mean pd5 = 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 (Cua­dros 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 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 pd5 of these suppressed live trees of 0.022.

We analysed the relationship between pd5 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 pd5 value of dead trees. We realized that a fixed value of pd5 as mortality threshold would not fit well to this condition and consequently defined a decreasing threshold survival function for pd5 as in [6] producing lower pd5 mortality values for trees with larger dbh values.

Figure 6. pd5 (relative dbh increment of the previous five years) versus the dbh values for live (green dots) and dead (red dots) trees. The red curve is a threshold mortality function fitted to the 0.05 quantile of live trees.

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]:

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]:

DBHi and 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):

yi is a dummy variable that takes the value 0 if tree i actually dies during the simulated period and takes the value 1 otherwise. pi 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.

pi is the proportion of trees belonging to the ith 10 cm-diameter class and R is the number of classes.


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.

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

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).

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.

Table 5. Growth model validation (plot V). Statistics refer to dbh estimations. 9-year simulation period.

Figure 8. Growth model validation (plot V, nine-year simulation period). Left: Observed vs expected dbh values of 2007. Right: Actual 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).

Figure 9. Simulation of dbh distributions (top) and basal area distributions by dbh classes (bottom), 60 year into the future. Simulation using the observed data from 2007 of plots F (left) and V (right).

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 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.

Figure 10. Current and projected development (60 years) of the Shannon index for 10cm- diameter classes. Plot F (left) and V (right).


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-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.

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.


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 growth-competition 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.

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.


Abellanas B, Coca M, Cuadros S, Oliet J, 2005. Análisis de la diversidad estructural del pinsapar puro en la sierra de grazalema. Influencia sobre la dinámica de la regeneración. Actas del IV Congreso Forestal Español. SECF, 4CFE05-004-T1: 115.
Abellanas B, Abellanas M, Vilas C, 2007. Vorest: Modelización de bosques mediante diagramas de Voronoi. Actas XII Encuentro de Geometría Computacional: 2249-2256. Valladolid, España.
Abellanas B, Abellanas M, Pommerening A, Lodares D, 2012. Vorest: un modelo de simulación basado en diagramas de Voronoi. Cuadernos de la Sociedad Española de Ciencias Forestales, 34: 27-37.
Alba-Sánchez F, López-Sáez JA, Benito-De Pando B, Linares JC, Nieto-Lugilde D, López-Merino L, 2010. Past and present potential distribution of the Iberian Abies species: a phytogeografic approach using fossil pollen data and species distribution model. Diversity Distrib 16: 214-218.
Antos JA, Parish R, Conley K, 2000. Age structure and growth of the tree-seedling bank in subalpine spruce-fir forests of south- central British Columbia. Am Midl Nat 143: 342–354.[0342:ASAGOT]2.0.CO;2.
Antos JA, Parish R, 2002. Structure and dynamics of a nearly steady-state subalpine forest in south-central British Columbia, Canada. Oecologia 130: 1, 126-135.
Arista M, 1995. The structure and dynamics of an Abies pinsapo forest in Southern Spain. Forest Ecol Manag 74: 81-89.
Aurenhammer F, Edelsbrunner H, 1984. An optimal algorithm for constructing the weighted Voronoi diagram in the plane. Pattern Recognition 17(2): 251-257.
Bauer S, Wyszomirski T, Berger U, Hildebrandt H, Grimm V, 2004. Asymmetric competition as a natural outcome of neighbor interactions among plants: results from the field-of_neighborhood modelling approach. Plant Ecol 170: 135–145.
Bella IE, 1967. Crown width/diameter relationships of open-growing jack pine on four site types in Manitoba, Canada. Dep. Forest & Rural Development Research Notes 23(1): 5-6.
Berger U, Hildenbrandt H, Grimm V, 2004. Age-related decline in forest production: modelling the effects of growth limitation, neighbourhood competition and self-thinning. J Ecol 92: 846-853.
Botkin DB, Janak JF, Wallis JR, 1972. Some ecological consequences of a computer model of forest growth. J Ecol 60: 849-872.
Coca M, Abellanas B, Cuadros S, Oliet JA, Miguel A, 2001. Caracterización estructural del pinsapar de la Sierra de Grazalema. Actas del III Congreso Forestal Español. Ed. Junta de Andalucía, TRAGSA, SECF: 257-264.
Cuadros S, Oliet JA, Abellanas B, Coca M, Padrón E, 2003. La regeneración en el pinsapar de la Sierra de Grazalema. II: Estructura y dinámica del regenerado consolidado en el pinsapar puro. Cuadernos de la Sociedad Española de Ciencias Forestales 15(2): 27-32.
Cuadros S, Ramírez A, Abellanas B, 2005. Epidometría basada en análisis de imagen y apoyada en un sig. Cuadernos de la Sociedad Española de Ciencias Forestales 19: 85-90.
Eichhorn MP, 2010. Spatial organisation of a bimodal forest stand. J Forest Res 15: 391–397.
Ek AR, 1974. Dimensional relationships of forest and open-grown trees in Wisconsin.Univ. Wisconsin Forest Research Notes 181. 7pp.
Farr WA, Demars DJ, Dealy JE, 1989. Height and crown width related to diameter for open-grown western hemlock and Sitka spruce. Can J Forest Res 19: 1203-1207.
Grabarnik P, Särkkä A, 2011. Modelling the spatial and space-time structure of foreststands: How to model asymmetric interaction between neighbouring trees. Procedia Environ Sci 7: 62–67.
Gray L, He F, 2009. Spatial point-pattern analysis for detecting density-dependent competition in a boreal chronosequence of Alberta. Forest Ecol Manag 259: 98-106.
Grimm V, Berger U, Bastiansen F, Eliassen S, Ginot V, Giske J, Goss-Custard J, Grand T, Heinz SK, Huse G, 2006. A standard protocol for describing individual-based and agent-based models. Ecol Modelling 198: 115-126.
Grimm V, Berger U, Deangelis DL, Polhill JG, Giske J, Railsback SF, 2010. The ODD protocol: A review and first update. Ecol Modelling 221: 2760-2768.
Hara T, 1992. Effects of the mode of competition on stationary size distribution in plant populations. Ann Bot 69: 509-513.
Hasenauer H, 1997. Dimensional relationships of open-grown trees in Austria. Forest Ecol Manag 96: 197-206.
Illian J, Penttinen A, Stoyan H, Stoyan D, 2008. Statistical Analysis and Modelling of Spatial Point Patterns. John Wiley &Sons, Ltd. West Sussex. 534 pp.
Krajiceck JE, Brinkman KA, Gingrich SF, 1961. Crown competition- a measure of density. Forest Sci 7: 35-42.
Leech JW, 1984. Estimating crown width from diameter at breast height for open-grown radiate pine trees in South Australia. Aust Forest Res 14: 333-337.
Linares JC, Camarero JJ, Carreira JA, 2009. Interacting effects of changes in climate and forest cover on mortality and growth on the southernmost European fir forests. Global Ecol Biogeogr 18: 485–497.
Linares JC, Carreira JA, 2009. Temperate-like stand dynamics in relict Mediterranean-fir (Abies pinsapo, Boiss.) forests from southern Spain. Ann Forest Sci 66: 610:1-10.
Linares JC, 2011. Biogeography and evolution of Abies (Pinaceae) in the Mediterranean Basin: the roles of long-term climatic change and glacial refugia. J Biogeogr 38: 619-630.
Lundqvist L, 1994. Growth and competition in partially cut sub-alpine Norway spruce forests in northern Sweden. Forest Ecol Manag, 65: 115-122.
MAGRAMA, 2015. 3º Inventario Forestal Nacional de España (IFN3). Ministerio de Agricultura, Alimentación y Medio Ambiente de España.
Matic S, 2001. Silver fir (Abies alba Mill) in Croatia. Ed: Akademija Sumarskih Znanosti – Hrvatske Sume, p.o. Zagreb. 895 pp.
McCarthy J, 2001. Gap dynamics of forest trees: a review with particular attention to boreal forests. Environ Rev 9, 1: 1-59.
Morin H, 1994. Dynamics of balsam fir forests in relation to spruce budworm outbreaks in the Boreal Zone of Quebec. Can J Forest Res 24, 4: 730-741.
Murphy ST, Pommerening A, 2010. Modelling the growth of Sitka spruce (Picea sitchensis (Bong)carr.) in Wales using Wenk’s model approach. Allg Forst J Ztg, 181(1/2): 35-43.
Newnham RM, 1964. The development of a stand model for Douglas-fir. PhD Thesis. Univ. British Columbia. Vancouver. Canada. 201 pp.
Nunes L, Tomé J, Tomé M, 2011. Prediction of annual tree growth and survival for thinned and unthinned even-aged maritime pine stands in Portugal from data with different time measurement intervals. Forest Ecol Manag 262: 1491-1499.
Oheimb G, Lang AC, Bruelheide H, David H, Forrester I, Wäsche I, Yu M, Härdtle W, 2011. Individual-tree radial growth in a subtropical broad-leaved forest: The role of local neighbourhood competition. Forest Ecol Manag 261: 499-507.
Okabe A, Boots B, Sugihara K, Chiu SN, 2000. Spatial Tessellations – Concepts and Applications of Voronoi Diagrams. 2nd edition. John Wiley, 671 pp.
Oliver CD, Larson BC, 1996. Forest Stand Dynamics. Ed: John Wiley and Sons, New York. 544 pp.
Paine DP, Hann DW, 1982. Maximum crown-width equations for southwestern Oregon tree species. School of Forestry Oregon State Univ. Forest Res Paper 46: 1-20.
Pommerening A, LeMay V, Stoyan D, 2011. Model-based analysis of the influence of ecological processes on forest point pattern formation- A case study. Ecol Modelling, 222: 666-678.
Pommerening A, Särkkä A, 2013. What mark variograms tell about spatial plant interactions. Ecol Modelling, 251: 64-72.
Pretzsch H, Biber P, Ďurský J, 2002. The single tree-based stand simulator SILVA: construction, application and evaluation. Forest Ecol Manag, 162(1): 3-21.
Pretzsch H, 2009. Forest dynamics, growth and yield. From measurement to model. Spriger-Verlag. Berlin. 664 pp.
R Development Core Team, 2011. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.
Römisch K, 1996. Durchmesserwachstum und ebene Bestandesstruktur am Beispiel der Kiefernversuchsfläche Markersbach. In: Hempel, G. (Ed.), 8. Tagung der Sektion Forstliche Biometrie und Informatik, pp. 84–103. Biotechnische Fakultät, Abteilung Forstwirtschaft, Ljubljana, Slovenia.
Scaltsoyiannes A, Tsaktsira M, Drouzas AD, 1999. Allozyme differentiation in the Mediterranean firs (Abies, Pinaceae).A first comparative study with phylogenetic implications. Plant Syst Evol, 216(3/4): 289-307.
Smith WR, Farrar RM, Murphy RA, Yeiser JL, Meldahl RS, Kush JS, 1992. Crown and basal area relationships of open-grown southern pines for modeling competition and growth. Can J Forest Res 22: 341-347.
Sterba H, Zingg A, 2001. Target diameter harvesting- a strategy to convert even-aged forests. Forest Ecol Manag 152(1-3): 95-105.
Stewart GH, Rose AB, Veblen TT, 1991. Forest development in canopy gaps in old-growth beech (Nothofagus) forests, New Zealand. J Veg Sci 2: 679–690.
Thorpe SC, Astrup R, Trowbridge A, Coates KD, 2010. Competition and tree crowns: A neighborhood analysis of three boreal tree species. Forest Ecol Manag 259: 1586–1596.
Toda M, Yokozawa M, Emori S, Hara T, 2010. More asymmetric tree competition brings about more evapotranspiration and less run off from the forest ecosystems: A simulation study. Ecol Modelling 221: 2887-2898.
Uchmansky J, 2000. Resource partitioning among competing individuals and population persistence: an individual-based model. Ecol Modelling 131:21–32.
Weiner J, 1986. How competition for light and nutrients affects size variability in Ipomoea tricolor populations. Ecology, 67: 1425–1427.
Weiner J, 1990. Asymmetric competition in plant populations. Trends Ecol Evol 5: 360-364.
Weiner J, Damgaard C, 2006. Size-asymmetric competition and size-asymmetric growth in a spatially explicit zone-of-influence model of plant competition. Ecol Res 1: 707–712.
Woods KD, 1984. Patterns of tree replacement: canopy effects on understory pattern in hemlock – northern hardwood forests. Vegetatio 56: 87–107.
Xiang QP, Xiang QY, Guo YY, Zhang XC, 2009. Phylogeny of Abies (Pinaceae) inferred from nrlTS sequence data. Taxon 58(1): 141-152.