Towards the sustainable management of thuya ( Tetraclinis articulata ( Vahl . ) Mast . ) forests in Tunisia : models for main tree attributes

The thuya (Tetraclinis articulata (Vahl.) Mast.) forests are one of the most important ecosystems in semiarid environments in north-western Africa, providing important economic profit and social services to local populations. However, lack of tools aiding sustainable management of these forests is detected. In the present work models for the main tree attributes as total height, crown diameter, height to crown base and stem form are developed for the species, using data from a net of plots installed in JbelLattrech region, in the NE Tunisia. Presented models allow characterizing the actual state and timber production of forests by using variables measured in typical forest inventories and conform a preliminary step for the future development of dynamic growth models.


Introduction
The thuya or araar (Tetraclinis articulata (Vahl.)Mast.) is a Mediterranean forest species from the family Cupressaceae which extends from north-western Africa to south-western Europe, covering a potential area of 2,500,000 ha, nowadays reduced to only 800,000 ha.
The main areas for the species are located in Morocco (600,000 ha), Algeria (160,000 ha) and Tunisia, where it occupies approximately 22,000 ha (Charco, 1999;DREF, 2002).In Tunisia, the geographical distribution of the thuya is limited to the northern-eastern region (fig.1), showing discontinuous forest patches from Cap Bon to Hammamet and Zaghouan areas (Rejeb et al., Models for Tetraclinis articulata tree-level attributes Considering the singularity of the species and the potential economic and social interest for the rural populations, it is necessary to define sustainable management schedules warranting wood and resin profit together with the preservation and improvement of the forests.The scientific literature on the species mainly focuses on phytosociological and autoecology studies (Nabil, 1989;Schoenenberger, 1995), with almost no reference to management practices (Bachoua and Voreux, 1987).The development of management schedules requires: i) tools for estimating actual yield of forests and ii) growth models which allow managers to predict forest evolution under different scenarios and management schedules.As a preliminary step for attaining both objectives, tree-level models relating common variables measured in forest inventories -as breast height diameter -with interesting tree attributes are demanded.In the case of thuya, the only available tree equations are the ancient one and two-entry volume equations developed for the forests in Moroccan Great Atlas (Achhal et al., 1985).In this context, the main objective of this work is to construct new models for predicting total tree height, crown diameter, height to crown base, stem taper curve and tree volume with flexible end-size products classification for thuya forests in NE Tunisia.
The thuya forests are sources of highly valuable services as soil erosion control, biodiversity conservation and CO 2 fixation, being also adequate for afforestation programs in arid and semiarid environments and highly-erosionated soils, especially under expected scenarios of climate change (Esteve-Selma et al., 2010).In addition, thuya timber is largely appreciated for construction, handcraft and cabinetmaking due to its hardness, veneer, durability and fragrance.Special interest has the root wood located in the enlarged coppice stump (lupia).A secondary traditional use has been resin tapping, from where the appreciated sandarac gum is obtained.
Timber overexploitation, inadequate resin tapping, grazing pressure and uncontrolled fires have resulted in highly degraded forests (Charco, 1999), only preserved from extinction due to the resprouting capacity of thuya.In the actual state of degradation lack of large size trees exists, so currently the main product from thuya forests is small woods for fuel and fencing (DREF, 2002), attaining very low productions per hectare.

Material
The present study was carried at JbelLattrech, a sandstone area highly representative for the species, covering about 4,000 hectares.The annual average rainfall is 393 mm, and average minimum-maximum temperatures for the coldest (January) and warmest month (August) between 8.4-30.5 °C respectively (Ben Mansoura and Garchi, 2001).The geology is dominated by Oligocene sandstone.Soils of the study zone are represented by soil of erosion on sandstone (calcaric fluvisols) and soil on dismantled calcareous crust (calcic cambisols).Within this area, a net of 50 permanent plots on pure even-aged stands of thuya was installed during the summer of 2009.Plots are circular, of variable radius (ranging between 5.3-16.3m), including 40 trees with breast height diameter > 5 cm, and were selected in order to cover the whole range of age, site quality and stand density detected.The stands are mainly coppices originated from clearcutting, though selective thinning is occasionally applied.Areas with signals (stumps) of recent tree harvesting (less than five years) were avoided.Every tree was positioned by using a tape and a compass, identifying trees coming from the same stump.Breast height diameter was measured with a caliper at each tree, while total tree height, height to crown base and crown diameter over two perpendicular directions were measured in a random sample of 10-12 trees per plot, using an hypsometer and a flexible tape.A single tree from the dominant stratum per plot was felled, and cut into logs at 0.5 m intervals, according with the recommendations by Prodan et al. (1997).Circumference length was measured with a flexible π-tape at the end sections of each log, and log volume was calculated using Huber's formula.Over bark total stem volume was obtained as the sum of the volumes of all the logs.Summary statistics for the sampled plots and trees are shown in Table 1.Height-diameter data, crown diameter, height to crown base data, and the stem profile for the data set are shown in Figure 2.

Height-diameter and crown attribute models
Twelve generalised height-diameter equations were selected as candidate functions to model the heightbreast height diameter relationship (Table 2).All functions tested are non-linear, constrained to pass through the point (0.0, 1.3) and include dominant height and diameter as covariates.Additionally, six biparametric nonlinear functions over breast height diameter were also evaluated for the crown diameter model (Table 2).All these functions were selected among those proposed  (2007).Height to crown base hcb was modeled as a function of total height using the nonlinear equation proposed by Dyer and Burkhart (1987), linearized by applying logarithmic transformation over the ratio between hcb and total height (see eqs. hcb1 and hcb2; Table 2).
The available fitting set consists of measurements taken from trees located within different plots.This hierarchical nested structure leads to lack of independence among observations coming from the same plot.
To alleviate this, candidate functions were fitted as multilevel nonlinear (dbh-h and crown diameter models) or linear (log(hcb/h)) mixed models, including random components at plot level.In the nonlinear models, first all the parameters were considered as mixed, and then reduced structures considering less random components were compared in terms of loglikelihood.In the case of linear mixed models, explanatory covariates were first selected using stepwise criteria, and then random plot components acting over intercept were included.Linear and nonlinear mixed models were carried using respectively MIXED and NLMIXED procedures in SAS/STAT® v.9.2 package.
Between models comparison was carried by using Akaike Information Criteria (AIC) together with efficiency (EF), root mean square (RMSE) and bias (E), evaluated over the conditional predicted value (including EBLUP's for the random parameters):

Stem curve and compatible volume function
The thuya stem form was described applying the compatible system of equations proposed by Fang et al. (2000) for slash and loblolly pine.This system includes a segmented taper function, and total tree and merchantable stem volume functions.The segmented taper function assumes that a tree stem can be divided into three parts, each of them with a different geometrical shape, separated by two inflection points occurring at unknown but estimable relative heights p 1 and p 2 , and its main expression is: Where I 1 = 1 if p1 ≤ q ≤ p2; 0 otherwise; I 2 = 1 if p2 ≤ q ≤ 1; 0 otherwise; and With d i : section diameter (cm) at height section h i ; dbh: breast height diameter (cm); h st : stump height (m); q = h i /h; h: total height (m); k=π/40000; a 0 , a 1 , a 2 , b 1 , b 2 , b 3 , p 1 , p 2 : parameters to estimate.From the previous taper model total stem volume V (m 3 ) is given by: And merchantable volume V i (m 3 ) up to a section diameter d i is directly given by: Parameters from the taper function [1] were independently estimated, and then used to develop volume functions [2] and [3] (method (ii) in Fang et al., 2000).Diameter section data are expected to be highly correlated, so a continuous second-order autoregressive structure was incorporated (Zimmerman and Nuñez-Antón, 2001).Dynamic fitting of model [1] considering autoregressive structure was carried using PROC MODEL on SAS/ETS ® v.9.2.package.Accuracy of the compatible system of equations was evaluated in terms of EF, RMSE, and E for section diameter, total volume and merchantable volume.In the latter, volume for the first two-meters log from the 14 sample trees with dbh > 15 cm was used to carry evaluations.

Height-diameter and crown attribute models
Table 3 includes comparison criteria for height-diameter and crown diameter models.No option considering two or more random parameters attained convergence.Among the converged functions for height-diameter, the three-parameter model h8 with parameter "a" as mixed attained the best performance: Where h = total height, ho: dominant height, Do: dominant diameter, dbh: breast height diameter, cw: crown diamenter, hcb = height to crown base, f(x) is a linear function of explanatory variables constrained to be < 0. Models for Tetraclinis articulata tree-level attributes -.
1 3 1 0 4760 1 Where u is a random plot term, normally distributed with mean zero and variance σ 2 (u)=0.0316;e is a residual error term, following a normal distribution with mean zero and variance σ 2 (e)=0.4599.
In the case of crown diameter model, all the function tested reached similar values, with model cw1 considering "b" as a mixed parameter attaining the best values in terms of EF and RMSE: Where v is a random plot term, normally distributed with mean zero and variance σ 2 (v)=0.0011;e is a residual error term, following a normal distribution with mean zero and variance σ 2 (e)=0.1475.
For height to crown base, selected linear model for logarithmic transformation over hcb/h included as explanatory covariates the inverse functions of the squared age, stand density and the squared breast height diameter, as well as the product among diameter-height ratio and the inverse of stand age: Where w is a random plot term, normally distributed with mean zero and variance σ 2 (w)=0.1125;e is a residual error term, following a normal distribution with mean zero and variance σ 2 (e)=0.2778.The expected marginal value for hcb after applying antilogarithmic transformation and bias correction (Flewelling and Pienaar, 1981)  Performance results for models [6] and [7] are also included in Table 3. for section diameter and tree volume, respectively (Table 3).In the case of merchantable volume [3], RMSE and E values are similar to those obtained for tree volume, though EF reaches 0.65, mainly indicating low variability in the values for the first two-meter log in the selected trees.

Discussion and conclussions
The developed tree-level models for thuya allow predicting variables whose measurement is largely expensive, as individual height, crown diameter or height to crown base, using explanatory covariates typically obtained in forest inventories, as breast height diameter.Other interesting attributes as crown volume (computed by assuming a geometrical shape for the crown) or competition indices can be easily derived.The formulation of these models as mixed, including a random plot component, allows to carry calibrations for new locations where a small sample of tree heights or crown size attributes are measured (see Calama and Montero, 2004, for details on model localization), improving the accuracy in the predictions.Moreover, the developed stem curve and compatible volume equations following the method proposed by Fang et al. (2000) permit a flexible classification of timber products, by defining the minimum end diameter of the logs for cabinet making, construction, fencing or fuel wood, what means an advantage over the traditional volume equations by Achhal et al. (1985).Additional benefits of the proposed modeling approach is that allow combining the stem curve and crown dimensions models with parameters as timber density for the species (790 kg/m 3 , Matouk, 2005) and crown bulk density to make estimates of aerial tree biomass and fixed CO 2 , which will help to determine the capacity of Tunisian thuya forests as CO 2 pools.Finally, these models can be seen as a preliminary step for the subsequent construction of growth and yield models which permit to predict the evolution of thuya stands under different scenarios and to define optimum scientific-based management schedules.

Figure 1 .
Figure 1.Thuya forests region in Tunisia and JbelLattrech study area.

Figure 2 .
Figure 2. Plot charts for total height and crown diameter as a function of breast height diameter (above); height to crown base as a function of total height and stem form profile relating relative diameter and relative height (below).

Table 1 .
Summary statistics for the data set used Where N: stand density, T: age, BA: basal area, dg: mean squared diameter, Do: dominant diameter, ho: dominant height, hm: mean stand height, Vol: total stand volume, SDI: Reineke's stand density index, dbh: breast height diameter; h: total tree height, cw: crown diameter, hcb: height to crown base, V: total tree volume.Models for Tetraclinis articulata tree-level attributes by Sánchez-González et al.

Table 2 .
Evaluated height-diameter and crown diameter models

Table 3 .
Fitting statistics and comparison criteria for the different models evaluated.Bolds indicate selected models Where Random: parameter of the models selected as random, n: number of observations, E: bias, RMSE: root mean square error, EF: modeling efficiency, AIC: Akaike's Information Criterion (see text for detailed description on these statistics).