Introduction 

Thuya (Tetraclinis articulata [Vahl.] Mast.) is a semiarid Mediterranean forest species which extends discontinuously from North Africa to south-western Europe. The natural character and the limited geographic distribution area of the species make it of particular ecological interest in terms of biodiversity conservation. Thuya forests provide highly valuable services such as soil erosion control, biodiversity conservation and CO₂ fixation, and are suitable for afforestation programs in arid or semiarid environments, where few species are able to grow ([48]), and areas with severely eroded soils ([19]).

In addition, these stands are of potential economic and social interest for the rural populations, not only in Tunisia where they cover 33000 ha ([16]), but also in Algeria and Morocco, covering more than 1000000 ha in Maghreb ([12]). Thuya timber is widely used in construction, handicraft and cabinetmaking due to its hardness, veneer, durability and fragrance. A secondary traditional use of this species is resin production to obtain sandarac gum. However, overexploitation for timber, irrational resin tapping, grazing pressure and uncontrolled fires have resulted in highly degraded forests ([12], [18]). The extinction of this species has only been avoided by its resprouting capacity. The current degraded state has left these forests devoid of large trees, so the main production presently comprises small pieces of wood for fuel and fencing ([17]). However, the profitability of these forests could be improved by exploiting other uses of the wood such as in industry, production of decorative objects, etc.

Improving the profitability of wood and resin production, along with the need to preserve and improve Tunisian thuya stands, which play a key ecological role in arid or semiarid environments, necessitates the development of models for quantitatively assessing the production of wood at different forest sites and under different management schedules to ensure a sustained yield. A first step towards this objective was the development of tree-level models relating common variables measured during forest inventories, such as breast height diameter, with relevant tree attributes like total tree height, crown diameter, height to crown base, stem taper curve and tree volume ([8]). The diameter-growth dynamics of the trees was studied by developing a distance-independent individual tree diameter increment model ([46]). A stand-level growth and yield model has recently been developed consisting of a site index sub-model and a system of compatible stand attribute equations presented in the form of a stand density management diagram ([47]). The next step would be to characterize the stand diameter structure.

Fonseca et al. ([21]) state that detailed information on the variability of tree diameters in natural and altered stands is key to sustainable forest management and is needed to assess carbon sequestration in forest ecosystems. Tree diameter distributions play an important role in stand management. For example, data regarding the number of trees per diameter class in a stand where a given silvicultural treatment has been applied are of particular interest because the diameter size will determine the industrial use of the wood and thus the price of the different products ([24]). Diameter distribution data at stand level also provides a more scientific basis for forest managers to choose which trees should be harvested ([45]).

Different types of parametric density functions have been used to describe tree diameter distribution in forest stands, including the Normal ([37]), Log-normal ([5], [37]), Gamma ([39]), Beta ([13], [42]), Johnson’s SB ([27], [42], [21]) and the Weibull distribution ([2], [42]). Mateus & Tomé ([33]) used Johnson’s distribution for modeling the diameter distribution of first-rotation eucalyptus plantations in Portugal. Gorgoso et al. ([25]) compared the accuracy of the Weibull, Johnson’s SB and beta distributions for modeling even-aged stands of Pinus pinaster Ait., Pinus radiata D. Don and Pinus sylvestris L. in northwest Spain, and Binoti et al. ([4]) evaluated the effectiveness of fatigue life, Frechet, Gamma, Generalized Gamma, Generalized Logistic, Log-logistic, Nakagami, Beta, Burr, Dagum, Weibull and Hyperbolic distributions for describing diameter distribution Tectona grandis L. stands subjected to thinning treatments at different ages in Brazil. Sghaier & Palm ([45]) used Pearson’s type I function for modeling diameter distribution of Pinus halepensis Mill. stands in Tunisia.

The objective of this study is to model the diameter distribution of the main Tetraclinis articulata stands in Tunisia. To achieve this objective, we compared the Normal distribution, given its ease of use, and the two-parameter Weibull distribution, since it has been shown to be the most suitable in many cases ([7], [34]).

  Material and ethods 

Study area and data collected

The natural Tetraclinis articulata stands at Jbel Lattrech forest in north-eastern Tunisia cover an area of approximately 4.000 hectares and are essentially mono specific and homogeneous. The climate is semi-arid, with an annual average precipitation of 393 mm, of which 95 % of rain falls between September and May. The average minimum temperatures for the coldest month (January) and maximum for the hottest month (August) are 8.4 and 30.5 °C, respectively ([3]). The study zone is located on Oligocene sandstones where there are groupings of Tetraclinis articulata - Lavandula stoechas, dominated by Cistus monspelliensis L., Genista aspalathoides Lam., Erica arborea L., Brachypodium ramosum Roemer et Schultes, Avena bromoides Gouan and Lupinis angustifolius L. ([36]).

The data used to develop the diameter distribution model were collected in 2009 from 50 temporary circular sample plots of variable radius and established such that they cover a wide range of age, density and site quality combinations. The size of the plots ranged from 88 to 835 m², depending on the stand density, so as to include at least 40 measured trees per plot. The diameter at breast height of all trees in each plot (down to a minimum diameter of 5 cm) was measured to the nearest 0.1 cm. Additionally, the heights of a random sample of 11 trees per plot (including the 4 largest trees) were measured to estimate mean height. Dominant height was calculated from the four largest trees per plot, and age was estimated from stem analysis on fallen trees by counting the number of rings at 0.30 m (see [8] for more details on data collection). Tab. 1 shows the main stand variables measured in the 50 plots.

Diameter distribution unctions

We tested two distribution functions: the Normal function, known for its ease of use ([29]), and the Weibull function, most frequently used because of its flexibility, simplicity and existence of a closed form solution for its CDF ([10], [35], [41], [24]).

Normal distribution

The equation for the normal probability density function is (eqn. 1):

where x is the random variable, and m and σ are its arithmetic mean and standard deviation, respectively.

The cumulative probability distribution of the normal function is (eqn. 2):

The two parameters which define the normal distribution (arithmetic mean and the standard deviation) can be estimated by using the following expressions ([15] - eqn. 3, eqn. 4):

where n indicates the number of trees per plot and x_i (cm) the diameter at breast height of each tree.

Weibull distribution

The variety of shapes that can be taken by the Weibull distribution makes it particularly suitable for representing tree diameter distribution. The equation for the Weibull probability density function is as follows (eqn. 5):

for x ≥ a, a ≥ 0, b > 0, c > 0 or 0 if not. The cumulative probability distribution of the Weibull function is (eqn. 6):

for x > a. This function is defined by three parameters: a is the location parameter, b is the scale parameter and c is the shape parameter.

Parameter a, commonly termed the location parameter, identifies the lower bound of the diameter distribution. For fixed values of b and c, changes in the parameter a simply shift the entire distribution along the x-axis. Parameter b is the scale parameter, and point x = a + b corresponds approximately to the 63rd percentile of the distribution. Parameter c defines the shape of the Weibull distribution. When c < 1, the distribution has a reverse-J shape. In the special case where c = 1, the Weibull distribution is reduced to the exponential distribution. When c > 1, the distribution is mound-shaped, approximately equivalent to the Normal distribution for c = 3.6. When c is between 1 and 3.6 the Weibull distribution is positively skewed, and is negatively skewed when c is greater than 3.6 ([2]).

The most frequently used methods to estimate the parameters of the Weibull distribution, are the maximum likelihood method ([28]) and the method of moments ([6]). Parameter a can be considered the smallest possible diameter in the stand and thus should be between zero and the minimum observed value in some cases ([2]). In some studies, the parameter a is arbitrarily fixed at 0.5 d_min ([30]) or at zero ([24], [44]), thus reducing the function to a two-parameter Weibull distribution which is easier to model and provides similar results to those of the three-parameter Weibull function, at least for even-aged, single-species stands ([34], [32]).

The two-parameter Weibull function (with the location parameter a = 0) was considered in this study as follows (eqn. 7):

where x is the random variable, b the scale parameter and c the shape parameter that controls the skewness.

Maximum likelihood estimator (MLE)

The maximum likelihood method is a commonly used procedure for the Weibull distribution in forestry because it has certain desirable properties ([30]). Estimation of the parameters using maximum likelihood has been found to produce consistently better goodness-of-fit statistics compared to other methods, but it also puts the greatest demands on the computational resources ([11]). If we consider the Weibull probability density function given in eqn. 7, then the likelihood function (L) will be (eqn. 8):

Taking the logarithms from eqn. 8, differentiating with respect to b and c respectively, and satisfying the following equations ([14], [38], [20] - eqn. 9, eqn. 10):

where n equals the number of sample observations in a Weibull distribution and x_i the diameter of each tree. The value of c must be obtained from eqn. 9 by using standard iterative procedures and then it is used in eqn. 10 to obtain b.

Method of moments (MOM)

The method of moments is another technique commonly used for parameter estimation. In the Weibull distribution, the k moment readily follows from eqn. 7 and is ([30] - eqn. 11):

where Γ (Gamma function) is Γ(s) = ∫ x^s-1 e^-x dx (s>0). Then from eqn. 11, we can find the first and the second moment as follows (eqn. 12, eqn. 13):

which gives (eqn. 14):

where σ² is the variance of tree diameters in a plot, and m₁, m₂ are the arithmetic and quadratic mean diameter in a plot, respectively.

When σ² is divided by the square of m₁, the expression for obtaining c is (eqn. 15):

In order to estimate b and c, we need to calculate the arithmetic mean diameter d and the variance σ² of the observed distribution and obtain the estimator of c in eqn. 15. eqn. 15 was resolved by an iterative procedure. When the value of the location parameter (a) is zero, the scale parameter (b) can then be calculated directly using the following equation ([24], [44] - eqn. 16):

where d is the arithmetic mean diameter.

Method comparison

To determine whether the adjusted theoretical distribution provides a good representation of what is observed, we used the chi-square test (χ²). This test requires some extreme classes to be grouped when there are an insufficient number of observations ([15]). For this reason, the following goodness-of-fit statistics were computed for each method: the mean bias, which reflects the deviation of the model with respect to observed values and the root mean-square error (RMSE), which is a measure of the precision of the estimates.

For the calculation of these statistics, the observed value Y_i and the theoretical value Y_i were calculated for each diameter class per plot (2 cm diameter class). In addition, the error index ([43]) was calculated (eqn. 17):

where Y_i and Y_i are the observed and predicted number of trees in diameter class k for the i-th plot. The sum includes all diameter classes in the i-th plot. In this error index, we can multiply the absolute difference between the observed and predicted number of trees per diameter class by the volume, using an equation developed by Sghaier et al. ([47]): v = 2.182 · 10^-2 · d^2.89, R² = 0.917. The diameter used in this equation corresponds to the middle diameter class. In this case (eqn. 18):

This index can also be expressed relatively as a percentage by dividing it by the total volume of the observed distribution in the plot ([29] - eqn. 19).

The latter error index was systematically calculated for each plot in this study.

Parameter prediction: distribution parameters with stand variables

In order to predict the parameters of the distribution functions directly from the stand variables using simple regression models, correlation analyses were carried out by estimating Pearson’s correlation coefficients between the estimated parameters corresponding to each plot and different stand variables. The stand variables chosen were: age (A), density (N), dominant height (H_d), mean height (H_m), quadratic mean diameter (d_g), arithmetic mean diameter (d), 25% percentile (P₂₅), median (P₅₀), 75% percentile (P₇₅) of the stand, and natural logarithm transformations for the quotient of the 25% percentile, the 50% percentile, the 75% percentile and the mean diameter by the quadratic mean diameter (LP₂₅, LP₅₀, LP₇₅, and Ld, respectively).

To predict the fitting parameters from the yield table using the Normal and Weibull distribution, the mean diameter and the variance must be known (first and second order moments of the distribution, respectively). The variance can be directly obtained from the arithmetic and the quadratic mean diameter ([22]) by the following expression (eqn. 20):

The mean diameter can be modeled with eqn. 21, which ensures the prediction of dÌis lower than the quadratic mean diameter (d_g - eqn. 21):

where X is a vector of independent stand variables in a fixed instant and β is the vector of parameters to be estimated ([24], [44]).

To select the model, the bias, the mean square error and the adjusted determination coefficient were computed.

  Results 

Tab. 2 shows summary descriptive statistics for the estimated parameters.

In order to adjust the models to predict the parameters of the theoretical distributions using the stand variables, Pearson’s linear correlations for estimated parameters were calculated (Tab. 3).

The highest Pearson correlation between the parameters of Normal and Weibull distribution and the stand variables was obtained with d_g and its natural logarithm transformations. The values obtained for Pearson’s linear correlation between the Ld parameter and the two estimated parameters of the Weibull distribution using the Maximum Likelihood and moments approaches were -0.463 and -0.459 for parameter b and 0.884 and 0.908 for parameter c. In the case of the shape parameter (c), the highest values obtained for Pearson’s linear correlation were for the variable Ld (Tab. 3). The values obtained for the coefficients of determination were 91% and 96% for models that relate the parameter c to the variable Ld (Tab. S1 in Appendix 1).

As regards the variance of the Normal distribution and the scale parameter (b) for the Weibull distribution, the results also show good linear relationships between these parameters and the measured stand variables. It appears that there is a logarithmic relation between the variance and quadratic mean diameter (d_g). The coefficient of determination was 78% (Tab. S1 in Appendix 1).

Note that the mean diameter is needed to estimate the shape parameter (c) in both parameter prediction methods of the Weibull distribution as well as in the Normal distribution. Since the yield tables generally only give the quadratic mean diameter (d_g), an equation was employed to predict the mean diameter from the quadratic mean diameter (d_g), the mean height (H_m) and the age (A) of the stand (Tab. S2 in Appendix 1).

According to the fitted prediction equation of the mean diameter (Tab. S2 in Appendix 1), the mean quadratic diameter (d_g), the mean height (H_m) and the age of the stand (A) must be known. Of these three parameters, only the mean height (H_m) cannot be directly obtained from yield tables and was therefore estimated by the fitting equation (eqn. 22) using the dominant height (H_d) disaggregated from yield tables as an independent variable (eqn. 22):

with R² = 0.989 and CV% = 11.03. The expected theoretical values for each plot are calculated using the three studied distributions with the two parameter prediction methods (estimated and prediction parameters). The observed, estimated and predicted values were grouped for each plot in 2 cm diameter classes. The results of the different comparison criteria are shown in Tab. 4.

In the parameter estimation, the Normal distribution is the most suitable for all the comparison criteria evaluated, with the exception of bias. Furthermore, the two parameter estimation approaches of the Weibull function gave similar results, although the maximum likelihood approach (MLE) was found to be slightly better.

In contrast to the parameter prediction in the first model construction approach (parameter estimation), the results obtained for parameter prediction revealed that the Weibull distribution was better for all the comparison criteria used. In this case, the maximum likelihood method also provided the best results, with the exception of the bias, for which the method of moments gave the lowest value.

In addition to comparing the distributions according to the Reynolds’ index error (EI %), which is calculated for each plot, an analysis of variance (ANOVA) was performed for each model construction method. The results obtained reveal that there is no significant difference between the three distributions in each of the two model construction methods (α = 0.75 for parameter estimation and α = 0.57 for parameter prediction, respectively).

The Weibull distribution, using the method of moments, occupied the same intermediate position with the exception of the bias, for which it gave the lowest value in both cases (estimated and prediction parameters). To analyze these results in greater depth, some additional graphs were produced. The first two, (Fig. 1 and Fig. 2), present the values of Bias and RMSE in each diameter class. Fig. 3 shows the mean value of the error index (EI %) in each density class.

According to the values for bias (Fig. 1), all the distributions underestimated the frequency of trees in the smallest and largest diameter classes, and overestimated them for the 10 cm diameter class under both parameter construction methods. With the exception of the 8 cm diameter class, for which the Weibull maximum likelihood method gives the highest value, the distributions compared present a similar level of precision, though the Normal distribution is slightly better, especially for parameter prediction.

As regards the mean square error (Fig. 2), although there are slight differences in the smallest diameter classes (lower or equal to 10 cm), the distributions gave similar results for diameter classes of 12 cm and above, being the Weibull distributions slightly superior in the case of parameter prediction.

For the Reynolds error index (EI %), Fig. 3 shows the variation of the mean value of this criterion according to stand density class (trees ha^-1). For both model construction approaches, the Weibull distribution, using the maximum likelihood method to estimate the parameters, appears to give the most stable results of the three studied methods. In addition to its stability across all the density classes, this method also provides the most accurate results for the parameter prediction approach (Fig. 3).

According to the distribution of trees per diameter class in each of the 50 measured plots, three types of distribution shape were observed (symmetric distribution, right dissymmetric distribution and reverse-J distribution). In order to better appreciate the capacity of the adjusted models to predict the number of trees per diameter class at plot level, one representative plot per distribution type was selected. Graphs for estimated vs. observed values were constructed for a symmetric distribution (Fig. S1 in Appendix 1), a right skewed distribution (Fig. S2 in Appendix 1) and a reverse-J shaped distribution (Fig. S3 in Appendix 1) for the parameter estimation and parameter prediction methods. The three figures reveal that all the studied models behave similarly, though the Weibull function appears to be marginally better than the others, in particular for the parameter prediction approach.

The use of the model for predicting the distribution of trees by diameter class

The diameter distribution models can be used independently by measuring a number of stand variables, or together with the yield table, which provides information that allows to predict the parameters of the distribution employed. The example which is provided here uses information from the yield table elaborated as part of this study for modeling the growth and production associated with the management and conservation of Tetraclinis articulata forests in Tunisia ([47]). Tab. 5 presents the predicted Weibull PDF (probability density function) parameters for different ages of the best site quality (H_d = 7 m at 50 years of age) in the case of low stand density (Reineke stand density index: SDI < 400).

Fig. 4 shows the curves associated with these distributions and Tab. 6 represents the distribution of trees (in number and as a percentage) per diameter class for various stand ages.

Based on the tree-level model attributes presented in Calama et al. ([8]), which includes height-diameter relationship, a crown attributes model and a stem curve equation, it is possible to obtain information about the distribution of end-use timber volumes, or height distribution by diameter class.

  Discussion 

Our Pearson’s correlations analysis results are very similar to the findings of Gorgoso et al. ([24]) for Betula alba stands in northwest Spain, using the two-parameter Weibull function, except for the parameter Ld, which expresses the logarithm of the mean diameter divided by the quadratic mean diameter. In addition, the good linear correlations which are obtained between the scale parameter (b) of the Weibull function and the stand characteristics are analogous to those observed by Campos & Leite ([9]), Leite et al. ([31]) and Binoti et al. ([4]). Linear models that relate b to the logarithm of the quadratic mean diameter (d_g) show a very high determination coefficient. Similar linear models for this relationship, with values of the adjusted determination coefficient close to 99%, were obtained by Alvarez González ([1]) for Pinus pinaster in Galicia and García-Guëmes et al. ([26]) for Pinus pinea stands in Valladolid (Spain).

The results obtained in this study reveal that the Normal distribution, though less flexible, provides results which are as satisfactory as those of the Weibull distribution. Lejeune ([29]) reported similar findings for Picea abies in Belgium. However, the Weibull distribution using the maximum likelihood method to estimate parameters appears to provide a slightly superior method in terms of RMSE, MAE and EI%, especially with the parameter prediction approach.

A number of authors have reported that the two-parameter Weibull function provides the most simple and accurate approach for modelling diameter distributions, including Gorgoso ([23]), and [24]) for birch stands (Betula alba L.) in north-western Spain, Maltamo et al. ([34]) for Pinus sylvestris and Picea abies (L.) Karst. stands in Finland, Alvarez González ([1]) for Pinus pinaster stands, and Condés ([14]) for different species in Spain. Furthermore, the two-parameter Weibull function has been used in a large number of related studies due to its flexibility and high correlation of its parameters with stand characteristics ([40], [4]). Therefore, given its flexibility and the very strong correlation between the PDF parameters and stand characteristics, the two-parameter Weibull distribution, using the maximum likelihood method, appears to be the most suitable for describing diameter distributions in this study.

Gorgoso et al. ([25]) and Zhang et al. ([49]) in their studies in north west Spain and the north-eastern part of North America, respectively, obtained better results with the 3-parameter Weibull distribution, using the maximum likelihood method rather than the moments approach.

  Conclusion 

In this study, a model has been developed for predicting the distribution of trees by diameter class using stand variables from Tetraclinis articulata forests in north-eastern Tunisia. The developed model, used in conjunction with the previously developed models for tree level attributes ([8]), the distance-independent individual tree diameter-increment model ([46]) and the stand based growth and yield model presented in Sghaier et al. ([47]), constitutes a useful tool for the sustainable management of Tetraclinis articulata forests in north-east Tunisia.

  Acknowledgements 

This study (Project Number: A/017275/08) was carried out within the framework of the bilateral co-operation between Tunisia (INRGREF) and Spain (INIA), and was financed by the Ministry for Scientific Research and Technology of Tunisia (MSRT) and the Spanish Agency of International Cooperation and Development (AECID), for which they are greatly acknowledged. We also thank Adam Collins for the English review of the manuscript.

  References

  Supplementary Material

Tab. S1 - Regression equations retained for the parameter prediction using the stand variables (50 plots).
Tab. S2 - Parameter estimations and statistics used for modeling the mean diameter (dÌ) of the stand.
Fig. S1 - Observed vs. estimated values for the parameter estimation and parameter prediction approaches according to each test function for the plot number 36 with a symmetrical distribution of tree diameters.
Fig. S2 - Observed versus estimated values for the parameter estimation and parameter prediction method according to each test function for the plot number 50 with a right dissymmetrical distribution of tree diameters (positively skewed).
Fig. S3 - Observed versus estimated values for the parameter estimation and parameter prediction approaches according to each test function for the plot number 44 with a reverse-J distribution of tree diameters.