Introduction
Thuya (Tetraclinis articulata (Vahl.) Mast.) is a semiarid Mediterranean forest species which extends discontinuously from North Africa to southwestern Europe. The largest forest areas are concentrated in Algeria, Morocco and Tunisia, covering more than 1 000 000 ha in Maghreb ([11]). In Tunisia, the geographical distribution of the Thuya is limited to the areas close to the Cap Bon, Boukornine and Zaghouan ([37]). The estimated total area of Thuya in Tunisia is approximately 33 000 hectares, with more than 22 000 ha of homogeneous stands ([17]). The natural character and the limited geographic distribution area of the species make it of particular ecological interest in terms of biodiversity and conservation. In spite of its slow growth, the Thuya is used for afforestation in arid regions, where few species are able to grow ([45]). It is also used in landscaping, particularly in Spain, thanks to its pyramidal form which confers a decorative character similar to that of the fir tree ([43]). Thuya forests provide highly valuable services such as soil erosion control, biodiversity conservation and CO_{2} fixation, as well as suitability for afforestation programs in arid or semiarid environments and areas with severely eroded soils ([24]).
Thuya timber is widely appreciated for construction, handicraft and cabinetmaking due to its hardness, veneer, durability and fragrance. The root wood of the enlarged coppice stump (lupia) is of particular interest. A secondary traditional use of this species is resin tapping to obtain the appreciated sandarac gum. However, overexploitation of timber, inappropriate resin tapping, grazing pressure and uncontrolled fires have resulted in highly degraded forests ([11], [22]), which have only avoided disappearance because of the resprouting capacity of Thuya. Due to the current degraded state of these forests, there is a lack of large trees; consequently, the main production from Thuya forests at present comprises small pieces of wood for fuel and fencing ([19]), the production and income per hectare being very low. These unfavorable conditions can be aggravated under scenarios of climate change ([23]).
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 aimed not only towards profitable wood and resin production but also the preservation and improvement of the forests. Unfortunately, scientific literature on this species is scarce, and mainly focused on phytosociological and autoecological studies aimed at determining the potential distribution area of the species and habitat preferences on a local scale ([32], [41], [34], [24]). However, in recent years a number of papers focusing on growth and yield modeling have been published ([7], [44]).
The development of management schedules requires: (i) tools for estimating current yield of forests; and (ii) growth models which allow managers to predict forest evolution under different scenarios and management schedules. Site index models and stand density management diagrams (SDMD) are among the preferred modeling tools at stand level. Site index, defined as the height of dominant trees at a given age, is the most commonly used method for indirectly estimating the site quality (mostly a combination of soil fertility and climate) in forest management ([9], [12]), and is generally considered the main component for developing stand level models and yield tables. A stand density management diagram (SDMD) is a set of static standlevel functions and models which are represented in a graphic way to show the relationships among stand yield, mean tree size, stand density and other stand attributes (e.g., ratio of merchantable timber volume to total timber volume, stand stability, etc.) at all stages of stand development ([35]). From a static perspective, SDMDs constitute an adequate alternative when no data from permanent plots are available, since they can be used in the same way as classical growth and yield tables  providing yield outcomes for management schedules at given densities  though SDMDs show several advantages ([3], [10]):
 they do not require a site index classification, so stand age is not required, and they are easy to apply by using traditional data recorded in management inventories;
 stand state is defined by means of dominant height and current stocking density;
 they do not require a specific submodel to predict current stand density as a function of dominant height and age (usually the weakest equation in static yield tables);
 they allow all the information and outcomes to be presented graphically;
 they are flexible enough to allow the simulation of initial spacing and different thinning schedules for stand density evolution in order to achieve different management objectives.
The main objective of the present work is to construct a standlevel growth and yield model to predict the development of the Tetraclinis articulata stands within the NE Tunisia region. The proposed model comprises a site index submodel and a system of compatible stand attribute equations presented as a SDMD. In this study, the site index submodel has been constructed using the Generalized Algebraic Differences Approach (GADA) model, which allows a flexible classification of sites in terms of stand age and dominant height. The nonlinear system of equations includes stand level attributes such as mean squared diameter, total stand volume, merchantable volume up to a flexible endsize diameter and aboveground biomass. In a previous study, treelevel attributes for the species such as total height, crown length, diameter and stem form were modeled ([7]) and more recently, a distance independent diameter increment model has been constructed ([44]). The aforementioned studies provide preliminary steps in the development of the present work. The site index submodel can be used to obtain dynamic projections of dominant height, thus it can be coupled with a SDMD to simulate the evolution of a stand in those cases where dynamic standlevel models are not available. As a subsequent step, the site index submodel and the developed SDMD were jointly used to characterize, compare and evaluate different silvicultural schedules for the multifunctional management of Thuya forests in Tunisia.
Materials and methods
Study area and experimental design
This study was carried out at Jbel Lattrech, a sandstone area highly representative for the species and covering about 4000 ha (Fig. 1). The annual average rainfall is 393 mm, and average minimummaximum temperatures for the coldest month (January) and the warmest (August) range between 8.4 and 30.5 °C, respectively ([6]). The geology is dominated by Oligocene sandstone. Soils in the study zone are represented by soil from the erosion of sandstone (calcaric fluvisols) and soil from dismantled calcareous crust (calcic cambisols  [17]). Within this area, a network of 50 temporary plots in pure evenaged stands of Thuya was established during the summer of 2009. The plots were circular, of variable radius (ranging between 5.3 and 16.3 m), in order to include 40 trees with breast height diameter > 5 cm in the plot. Plots were selected in order to cover the whole range of stand development, site quality and stand densities detected in the study area. Apart from those criteria, plots were uniformly distributed within the territory, and separated from each other more than 500 m in order to prevent spatial correlation. The stands were mainly coppices originating from clearcutting, although selective thinning was occasionally applied. Areas with signs of recent harvesting (with tree stumps < 5 years old) were avoided. Each tree was positioned using a tape and a compass, and trees growing from the same stump were identified. Breast height diameter was measured in two perpendicular directions using a caliper on each tree, while total tree height, height to crown base and crown diameter in two perpendicular directions were measured on a random sample of 1012 trees per plot, using a hypsometer and a flexible tape. Measurements were taken once just after plot installation, and a new inventory is expected to be made every 10 years.
One dominant tree per plot was felled and discs were extracted at stump height (0.30 m) and at 0.5 m intervals along the stem. The number of rings was counted at each cross sectioned point and then converted into stump age. The age at a given level of the stem corresponds to the difference between the current age and the number of rings counted up to this level. As the crosssection lengths did not coincide with periodic height growth, we adjusted the height/age data from the stem analysis to account for this bias by using the Carmean method ([8]), and the modification proposed by Newberry ([33]) for the topmost section of the tree, based on earlier studies ([21]). These corrections remove the bias when we assume that the height of the section is the maximum height attained at a given age.
Summary statistics for the main variables of interest in the sampled plots are shown in Tab. 1.
N
): stand density; (t
): age; (SI
): site index; (BA
): basal area; (Dg
): mean squared diameter; (D
_{0}): dominant diameter; (H
_{0}): dominant height; (Hm
): mean height; (V
): total standing volume; (V
_{10}): merchantable volume up to an end diameter of 10 cm; (V
_{15}): merchantable volume up to an end diameter of 15 cm; (W
_{T}): total aboveground biomass; (SDI
): stand density index ([36]).
Parameter  Mean  STD  Min  Max 

N (stems ha^{1}) 
1860  865  479  4533 
t (years) 
57  10  36  77 
SI (m) 
5.06  1.08  2.83  7.78 
BA (m^{2}/ha) 
12.90  6.39  2.32  27.74 
Dg (cm) 
9.41  1.55  6.32  12.95 
D _{0} (cm) 
12.76  2.32  7.85  17.43 
H _{0} (m) 
6.52  1.38  4.20  11.00 
Hm (m) 
5.46  1.09  3.61  9.17 
V (m^{3}/ha) 
36.14  22.13  5.16  111.61 
V _{10} (m^{3}/ha) 
16.76  14.34  0.195  73.46 
V _{15} (m^{3}/ha) 
3.43  4.96  0  29.33 
W _{T} (t /ha) 
45.05  26.09  6.22  128.98 
SDI

384  183  75  762 
Volume and biomass estimation
The stem curve taper function and the heightdiameter relationships developed by Calama et al. ([7]) were jointly used to estimate stem volume for all the trees within the plot. Total standing volume was thus computed as the sum of single tree volumes over the whole plot. In a similar way, for each plot a stem taper curve was also used to compute stand merchantable volume up to different end size diameters ranging from 1 to 23 cm. The ratio (R
_{i}) between stand merchantable volume up to a flexible end diameter d
_{i} and total stand volume was also computed. Stem dry biomass was derived from single tree stem volume using the basic wood density for the species (590 kg m^{3}) estimated by ElMouridi et al. ([23]). Crown biomass was estimated at tree level as the product of crown volume and crown bulk density, defined as the mass of available canopy fuel per unit of crown volume ([42]). Crown volume (V
_{crown}) was computed for each tree assuming that the Thuya crown has a semiellipsoidal shape, thus crown volume equals (eqn. 1):
where d
_{crown} and h
_{crown} are the crown diameter and the height to the crown base, respectively. These values were observed in sample trees and estimated for the rest of trees within the plot by using the models developed by Calama et al. ([7]). As regards crown bulk density, since no specific data for Tetraclinis articulata were found in the literature, we used the average value found for the species of the family Cupressaceae ([38]), which is 0.06 lb ft^{3}, or in metric units, 0.9612 kg m^{3}. Tree aboveground biomass was then computed as the sum of stem biomass and crown biomass. Total aboveground biomass at plot level was then computed as the sum of the individual tree biomass over all the trees within the plot (Tab. 1).
Site index
In this study we propose to construct a site index model for Tetraclinis articulata stands in northeastern Tunisia using the methodology based on the generalized algebraic differences approach (GADA) suggested by Cieszewski & Bailey ([12]) and using the dummy approach to make the adjustment ([15]).
Productivity is assumed to be dependent on an unobservable sitespecific parameter X
which can neither be measured nor defined ([12]), but depends on soil conditions, and ecological and climatic factors. While GADA assures baseage invariance of the derived algebraic forms, the baseage invariance of the model parameter estimates ([5]) was assured by fitting them using the dummy variable approach ([15]). In the present work, six dynamic GADA equations (M1M6, see Tab. S1 in Appendix 1) with two sitespecific parameters were tested. These GADA models were derived from three different growth functions (each one bearing three parameters) commonly used for the development of site index curves and dominantheight growth modeling. Models M1M2 were derived by Cieszewski ([13]) from the loglogistic model, which is equivalent to the Hossfeld model; model M3 was derived by Krumland & Eng ([29]) from the BertalanffyRichards’ function while models M4M5 were derived from the same base equation by Cieszewski ([14]). Finally, model M6 was derived by BarrioAnta et al. ([4]) from the LunqvistKorf’s function. More information on GADA methodology and dummy approach fitting technique can be found in GeaIzquierdo et al. ([25]) or DiéguezAranda et al. ([18]).
Site index model selection was based on statistical criteria derived from classical residual analysis: root mean square error (RMSE), mean residual (Bias) and adjusted coefficient of determination. Models were also compared in terms of Akaike’s information criteria differences (AICd), a measure of the relative information lost when a given model is used to represent the process that generates the data (eqn. 2):
where n
is the total number of observations, p
is the number of parameters in the model and σ
^{2} is the estimator of the error variance of the model.
Another important step in evaluating the models was the graphical analysis of the residuals and examination of the appearance of the fitted curves overlaid on the trajectories of the dependent variables for each plot. Visual or graphical inspection is an essential step in selecting the most appropriate model because curve profiles may differ drastically, even though fit statistics and residuals are similar ([28]).
The fact that data originate from stem analysis create an autocorrelation among observations within the same tree (correlation between the residuals within the same tree). Parameter estimates are not biased by this autocorrelation, but estimators of standard errors are ([26]), which invalidates the standard hypothesis testing. To account for this possible autocorrelation, we modeled the error terms using a continuoustime autoregressive error structure [CAR(x)], which expands the error terms in the following way ([48]  eqn. 3):
where e
_{ij} is the jth ordinary residue on the ith tree (i.e., the difference between the observed and the estimated heights of the tree i at age measurement j), with d
_{n} = 1 for j > n
and it is zero for j ≤ n, ρ
_{n} is the norder autoregressive parameter to be estimated, and t
_{ij}t
_{i(jn)} is the time distance (years) separating the jth from the j
thn
observations. To evaluate the presence of autocorrelation and the order of the CAR(x) to be used, graphs representing residuals versus lagresiduals from previous observations within each tree were visually examined. The dummy variables method and the CAR(x) error structure were programmed by the use of the SAS/ETS^{®} MODEL procedure ([40]), which allows for dynamic updating of the residuals.
Construction of the Stand Density Management Diagram
A key issue in constructing a SDMD is to characterize the growing stock relating average tree size to density, by means of density indexes such as the stand density index ([36]), relative spacing index ([46]) or the selfthinning rule ([47]). In the present work we propose to construct a SDMD using Reineke’s stand density index SDI
, which relates number of trees per ha (N
) to mean squared diameter (Dg
 eqn. 4):
SDI
is a useful index to define management schedules and thinning rules. Stand development in a SDMD is shown by representing dominant height (H
_{0}) on the xaxis and number of trees per hectare (N
) on the yaxis. The main idea behind SDMD is to express all the mathematical relations in the system as a function of N
and H
_{0}, and therefore to define isolines for different levels of stocking density (SDI
) as well as the rest of the stand attributes of interest by solving for N
. Thus all the functions and models included in the system can be solved in terms of these two variables, in order to plot isolines for the different values of the attribute variables over the diagram. Another key concern in constructing a SDMD is to determine the main basic functions for defining the stand state. In our SDMD we propose to fit the following system of nonlinear equations (eqn. 5  eqn. 8):
The eqn. 5 relates mean squared diameter (Dg
) as a function of stand density (N
) and dominant height (H
_{0}), while eqn. 6 and eqn. 7 define the aboveground biomass (W
_{T}) and the stand volume (V
) as an allometric relationship between the average tree size (Dg
), the dominant height of the stand and the number of trees per hectare. eqn. 8 represents the merchantable volume ratio up to a given diameter D
_{i} (R
_{i}) as a function of Dg
. eqn. 7 and eqn. 8 can be reexpressed in terms of total merchantable volume up to a section diameter D
_{i} (Vm
_{i}), defined as the product among V
and R
_{i }(eqn. 9):
where a
_{0}a
_{2}, b
_{0}b
_{3}, c
_{0}c
_{3}, d
_{0}d
_{2} are parameters to be estimated. eqn. 5, 6 and 9 constitute a simultaneous system of equations, where N
and H
_{0} are exogenous variables and Dg, W
_{T} and Vm
_{i} are endogenous variables. Error components of the variables on the lefthand and righthand side are correlated, so a simultaneous fit of the system of equations was carried out by using the fullinformation maximum likelihood technique. Since the number of observations in the database was different for mean squared diameter and aboveground biomass (one observation per plot) and merchantable volume (several observations per plot, one for each 1 cm top diameter threshold up to 23 cm), a special structure was proposed that assigns to each merchantable volume from the same plot the values of the rest of stand variables of the plot, and the inverse of the number of observations per plot was used as weighing factor (see [10] for more details in fitting techniques). A simultaneous fit of the models was carried out by using again the SAS/ETS^{®} MODEL procedure ([40]).
Once the nonlinear system of simultaneous equations is solved, the next step includes the definition of isolines for the different levels of SDI, Dg
and W
_{T} (Tab. S2 in Appendix 1). As Vm
_{i} is not easily solved for N
and H
_{0}, we split eqn. 9 to define isolines for V
and R
_{i}. To draw the isolines for R
_{i} we selected an end diameter of 15 cm, since logs over this diameter can be used for handicraft purposes ([22]).
Results
Site index curves
To determine the order of the CAR(x) to be used to correct the autocorrelation, we first fitted all the six studied models using nonlinear least squares without expanding the error terms to account for autocorrelation. It appeared that a trend in residuals as a function of lagresiduals within the same tree is evident in all models, as expected because of the longitudinal nature of the data used for model fitting. Fig. 2 (first column) provides an example with M5. After fitting the models by a secondorder continuoustime autoregressive error structure CAR(2), the trends in residuals disappeared (Fig. 2  third column).
The parameter estimates for each model and their corresponding goodnessoffit statistics are shown in Tab. 2. The two models M3 and M6 present a non significant parameter (b
_{3} for M3 and b
_{2} for M6). For the other four models, all the parameters (including the sitespecific parameters for each tree) were highly significant, although model M4 showed worse performance in terms of RMSE, R^{2}_{adj} and AICd, thus further comparison is focused on models M1, M2 and M5. Graphs showing bias and root mean squared error in height estimations by age classes were analyzed (Fig. S1 in Appendix 1). In the case of bias, model M1 performed slightly better, both at juvenile and at old ages, although the values were higher for the remaining ages. For RMSE, models M1 and M5 displayed a similar behavior, showing lower values than those obtained by M2, especially for the older ages. An additional key criterion for comparison is based on model behavior at ages close to the maximum longevity of the species, which is over 250 years ([39]). At this overmature age, differences between models can be more than 2 m (Tab. 2), with model M5 presenting a horizontal asymptote of 19 m, which seems to best describe the growth of the species under the Tunisian conditions.
Model  Par  Est  SE  t  p > t  RMSE  Bias  R ^{2} _{adj}  AICd  H_{250} 

M1 
b
_{1}

32.24E5  428178  7.51  <0.0001  0.2238  0.031  0.9863  4.29  21.3 
b
_{2}

4.467E8  5.48E14  8.15E21  <0.0001  
b
_{3}

0.912  0.031  29.87  <0.0001  
ρ
_{1}

0.947  0.013  74.02  <0.0001  
ρ
_{2}

0.907  0.013  69.41  <0.0001  
M2 
b
_{1}

45.089  11.709  3.85  <0.0001  0.2257  0.033  0.9860  13.02  21.4 
b
_{2}

4.002E25  6.04E48  6.63E72  <0.0001  
b
_{3}

0.868  0.027  32.71  <0.0001  
ρ
_{1}

0.947  0.013  73.14  <0.0001  
ρ
_{2}

0.911  0.013  71.81  <0.0001  
M3 
b
_{1}

8.16E3  0.002  4.67  <0.0001  0.2229  0.030  0.9864  0  19.8 
b
_{2}

1.082  0.114  9.48  <0.0001  
b
_{3}

0.527  0.318  1.66  0.0982  
ρ
_{1}

0.948  0.013  74.40  <0.0001  
ρ
_{2}

0.908  0.013  69.68  <0.0001  
M4 
b
_{2}

5.857E3  0.002  3.34  0.0009  0.2303  0.034  0.9855  33.86  20.1 
b
_{3}

2.598  0.096  27.08  <0.0001  
ρ
_{1}

0.948  0.013  74.14  <0.0001  
ρ
_{2}

0.911  0.013  72.79  <0.0001  
M5 
b
_{2}

7.903E3  0.002  4.52  <0.0001  0.2256  0.033  0.9860  12.09  19.0 
b
_{3}

0.539  0.015  35.14  <0.0001  
ρ
_{1}

0.947  0.013  74.06  <0.0001  
ρ
_{2}

0.907  0.013  70.07  <0.0001  
M6 
b
_{2}

199.063  115.70  1.72  0.0862  0.2263  0.032  0.9860  14.95  25.7 
b
_{3}

0.065  0.023  2.90  0.0040  
ρ
_{1}

0.947  0.013  73.39  <0.0001  
ρ
_{2}

0.910  0.013  71.85  <0.0001 
Thus, among the six models evaluated for height growth prediction and site classification of Tetraclinis articulata stands in northeastern Tunisia, Model M5 (the GADA formulation derived from the BertalanffyRichards base function by considering a
_{1} and a
_{3} as related to site productivity) was selected (eqn. 10):
with (eqn. 11):
and (eqn. 12):
where Y
is the predicted dominant height (meters) at age t
(years), and Y
_{0 }and t
_{0} represent the predictor dominant height and age. Note that the sitespecific parameters are discarded in a way similar to that used to discard the autocorrelation coefficients, because the general use of the model involves making predictions using observed height and its associated measurement age in new individuals ([18]). Dynamic equations derived using the GADA allow direct prediction of dominant height to be made at any stand age from any other reference height and age. Similarly, site index at a reference age can be estimated from current dominant height. The best reference age to define the site index is determined by plotting relative prediction errors against prediction age ([18]). Fig. S2 in Appendix 1 suggested the selection of site index from ages ranging 4065 years (lowest relative errors), while a decrease in the number of observations after 40 years leads us to propose a reference age of 40 years for Tetraclinis articulata stands in NE Tunisia. Site index SI
was then defined as the dominant height reached at 40 years of age. Since site index is a fixed stand attribute which should be stable over time, a plot of the site index predictions against total age and the stem analysis data was also developed (Fig. S3 in Appendix 1), revealing the consistency of site index predictions over time except for the youngest ages, under 10 years.
Based on all the previous findings, a family of polymorphic curves for site index 3, 5 and 7 m at reference age of 40 years was constructed and overlaid on the trajectories of the observed heights over time for the selected model (Fig. 3). Based on this, we propose to define three site classes, limited by the midinterval curves for site indexes at 4 and 6 m.
Stand Density Management Diagram
Parameter estimates for the simultaneous fit of the nonlinear system of eqn. 5, eqn. 6 and eqn. 9 are shown in Tab. 3. All the parameters in the model are significant at <0.0001 level. Models for W
_{T} and V
_{m} explained the main part of the observed variability, while the model for Dg
only explained 0.51 of observed variability. However, the three models lead to unbiased predictions (we cannot reject the null hypothesis of bias equaling zero) and low values of RMSE.
Eqn.  Response variable 
Par  Est  SE  RMSE  Bias  R ^{2} _{adj}  pvalue 

[5] 
Dg

a
_{0}

5.9277  0.3564  1.0665  0.0177  0.5137  <0.0001 
a
_{1}

0.5567  0.0227  
a
_{2}

0.0777  0.0063  
[6] 
W
_{T}

b
_{0}

0.00003  3.98E6  2.6699  0.0122  0.9896  0.7391 
b
_{1}

2.2495  0.0522  
b
_{2}

0.4747  0.0277  
b
_{3}

1.0711  0.0072  
[9] 
V

c
_{0}

0.00002  0.2E5  1.8680  0.0559  0.9916  0.3673 
c
_{1}

2.3262  0.0415  
c
_{2}

0.5416  0.0222  
c
_{3}

1.0814  0.0091  
R
_{i}

d
_{0}

0.5266  0.0615  
d
_{1}

3.1928  0.0298  
d
_{2}

3.0749  0.0620 
Parameter estimates were used to develop SDMDs by superimposing the isolines for the different attributes on a bivariate axis defined by H
_{0} and N
. The range of values for axis and isolines is close to the range of observed values in the fitting dataset. These ranges were 2005000 trees ha^{1} for stand density, 212 m for dominant height, 5120 m^{3} ha^{1} for stand volume, 5140 tons of dry matter ha^{1} for aboveground biomass and 0.050.40 for R
_{i}. SDI
isoline of 700 would define a line close to the maximum biological stand density (the identified maximum SDI
value was close to 750), while values over 600 would indicate densities where mortality due to selfthinning is expected. Thus, the bivariate space above SDI
isoline 700 is not considered for defining management schedules. The above equations defining the SDMD can be represented in different ways: Fig. 4a shows the classical management diagram including the total volume, the mean squared diameter and the stand density index evolution; Fig. 4b represents the aboveground biomass development; while Fig. 4c represents the total volume and ratios of merchantable volume up to an end diameter of 15 cm.
Discussion
Fitted models
The site index model developed in this work to describe the dominant height growth of Tetraclinis articulata stands at Jbel Lattrech in northeastern Tunisia represents a substantial advance as regards modeling the growth of such species, because no previous site index model existed. The model selection procedure indicated that the generalized algebraic differences approach (GADA) derived from the BertalanffyRichards’ base equation using the dummy variable technique to make the adjustment, resulted in the best compromise between biological and statistical aspects, producing the most adequate site index curves. It has all the desirable properties of a dominant height function ([13], [18]), providing time invariant predictions, thus ensuring that projections over the same period of time are equivalent, regardless of the length or number of projection intervals. This dynamic equation can compute predictions directly from any agedominant height pair without compromising the consistency of predictions, and the predictions are unaffected by changes in the baseage. The site equation is mathematically sound, always computes consistent predictions, and has the property of being continuous and reciprocal.
The site index curves allow the definition of 3 site qualities (Fig. 3), defined by mean dominant heights of 7, 5 and 3 m at 40 years of age, with asymptotic values of about 19, 14 and 9 m, respectively. These asymptotic values seem to slightly exceed the real height values that Tetraclinis articulata can reach under natural conditions in Tunisia. However, this has not great importance since age, when asymptote is reached, is far from natural longevity of the species. According to the data from the sampled plots, 16% of the studied stands have a dominant height which does not exceed 4 m at 40 years of age, 14% present a dominant height which exceeds 6 m at the same age and, finally, 70% have a dominant height between 4 m and 6 m. In this regard, the proposed curves constitute per se an interesting tool for classifying Thuya stands in the region, providing a basis for developing management plans.
The site index model provides the basis for a family of standlevel models, including mean squared diameter, total and merchantable volume, and aboveground biomass, which together constitute a useful tool for the management of the stands by facilitating the construction of stand density management diagrams and yield tables ([1], [16]). In our case, three nonlinear models for predicting mean squared diameter, merchantable volume up to a flexible end diameter (compatible with total volume) and total aboveground biomass were selected from the existing literature ([3], [10], [31]) and were simultaneously fitted to obtain accurate and uncorrelated estimates of the parameters. While the mean squared diameter and total volume models are based on real observed data, total aboveground biomass is based on estimates of stem and crown biomass, derived from applying timber basic density and crown bulk density obtained from the literature ([23], [38]) to the observed real data for stem and crown volume. Although this method is a firstdraft approach, in the absence of more accurate data, it can be considered adequate for estimating and fitting standlevel biomass models.
Proposed management schedules for Thuya forests in NE Tunisia
Defining management schedules within the framework of site index models and a stand density management diagram requires different stand management factors:
 The final state of the stand according to the main objectives of the forest management. This final state can be defined in terms of stocking volume, average tree size or stand attributes (e.g., mean squared diameter, dominant height).
 Rotation age derived from the specific objectives of the management and the expected products. In our modeling context, it can be derived from the objective dominant height and the site index by using the site presented index model.
 Upper and lower stocking density limits which define thinning regimes. According to the proposal by Long ([30]), adequate thinning regimes would vary between 30% and 60% of
SDI
_{max}. Based on this, in our case we propose thinning regimes withSDI
varying between 250 and 500, avoiding higher stocking levels in order to prevent losses due to selfthinning.  The main objectives of forest management: at the best sites, it is possible to carry out intensive silviculture oriented towards obtaining timber dimensions by the rotation age which would be suitable for handicraft and cabinetmaking. Such uses require timber pieces with an end diameter of at least 15 cm ([22]). Other less intensive thinning regimes for either high or low quality sites could be proposed where the objective of the management is to produce wood suitable for fencing, fuel or pulp ([27]). In that case, we must consider that in low quality sites even extending rotations does not result in valuable timber pieces, thus shorter rotation lengths are proposed.
Based on the abovementioned constraints, we propose to evaluate four different management schedules (Tab. 4) for midhigh (SI
= 6 m) and midlow (SI
= 4 m) site quality stands. In each case, we assume a young stand (20 years old), with an initial stand density of 4000 stems ha^{1}, derived from coppicing. Initial dominant height of the stand, according to the site index model, is 2.26 m and 3.49 m for site indexes 4 and 6, respectively.
H
_{o}): dominant height (m); (N
): stems ha^{1}; (Dg
): mean squared diameter (cm); (AB
): basal area (m^{2} ha^{1}); (Vol
): standing volume (m^{3} ha^{1}); (W
_{T}): standing aboveground biomass (t ha^{1}); (R
_{i15}R
_{i10}): rate of merchantable volume up to an end diameter of 1510 cm; (N_ext
): stems ha^{1} extracted; (V_ext
): extracted volume (m^{3} ha^{1}); (W_ext
): extracted aboveground biomass (t ha^{1}); (Cum_V
): cumulative (standing + extracted) volume (m^{3} ha^{1}); (Cum_V
): cumulative (standing + extracted) biomass (t ha^{1}).
Schedule  Age 
H
_{o}

N

Dg

AB

Vol

W
_{T}

R
_{i15}

R
_{i10}

N_ext

V_ext

W_ext

Cum_V

Cum_W


1.1  20  3.5  4000  6.2  12.2  22  28  0.00  0.05  0  0  0  22  28 
30  4.8  4000  7.5  17.5  40  49  0.00  0.18  1500  14  17  40  49  
30  4.8  2500  7.7  11.8  26  32  0.00  0.22  0  0  0  40  49  
40  6.0  2500  8.8  15.0  39  47  0.02  0.35  0  0  0  52  63  
50  7.1  2500  9.6  18.0  52  62  0.06  0.46  1000  19  23  66  78  
50  7.1  1500  10.0  11.7  33  39  0.08  0.50  0  0  0  66  78  
60  8.0  1500  10.7  13.5  42  49  0.13  0.57  0  0  0  75  88  
70  8.9  1500  11.3  15.2  51  58  0.18  0.63  500  15  18  84  98  
70  8.9  1000  11.7  10.8  35  41  0.21  0.65  0  0  0  84  98  
80  9.7  1000  12.3  11.9  41  47  0.26  0.69  0  0  0  90  104  
90  10.4  1000  12.8  12.9  47  53  0.31  0.72  1000  47  53  95  110  
1.2  20  3.5  4000  6.2  12.2  22  28  0.00  0.05  0  0  0  22  28 
30  4.8  4000  7.5  17.5  40  49  0.00  0.18  2000  18  23  40  49  
30  4.8  2000  7.9  9.8  21  26  0.01  0.24  0  0  0  40  49  
40  6.0  2000  8.9  12.5  32  38  0.03  0.37  0  0  0  50  61  
50  7.1  2000  9.8  14.9  43  51  0.07  0.47  2000  43  51  61  73  
2.1  20  2.3  4000  4.9  7.5  10  13  0.00  0.00  0  0  0  10  13 
30  3.2  4000  5.9  11.0  18  24  0.00  0.03  1000  4  5  18  24  
30  3.2  3000  6.1  8.6  14  18  0.00  0.04  0  0  0  18  24  
40  4.0  3000  6.9  11.2  22  27  0.00  0.11  1000  7  8  26  33  
40  4.0  2000  7.1  7.9  15  19  0.00  0.14  0  0  0  26  33  
50  4.8  2000  7.8  9.6  21  25  0.00  0.23  500  5  6  31  39  
50  4.8  1500  8.0  7.5  16  20  0.01  0.25  0  0  0  31  39  
60  5.4  1500  8.6  8.8  20  25  0.02  0.34  500  6  8  36  44  
60  5.4  1000  8.9  6.2  14  17  0.03  0.37  0  0  0  36  44  
70  6.1  1000  9.5  7.0  17  21  0.05  0.44  1000  17  21  39  48  
2.2  20  2.3  4000  4.9  7.5  10  13  0.00  0.00  0  0  0  10  13 
30  3.2  4000  5.9  11.0  18  24  0.00  0.03  0  0  0  18  24  
30  3.2  4000  5.9  11.0  18  24  0.00  0.03  0  0  0  18  24  
40  4.0  4000  6.7  14.2  28  35  0.00  0.10  0  0  0  28  35  
50  4.8  4000  7.4  17.2  39  47  0.00  0.18  4000  39  47  39  47 
 Schedule 1.1 is proposed for the highest fertility sites, and is based on an intensive thinning regime and large rotation age in order to obtain high quality timber for use in handicraft. The target stand is a 90 yearold stand with 1000 stems ha^{1}, a dominant height of 10.5 m, standing volume of about 47 m^{3} ha^{1}, and 30% total volume (15 m^{3} ha^{1}) of timber with an endsize diameter exceeding 15 cm (for end use in handicraft). Three thinnings are proposed at the ages of 30, 50 and 70 years (dominant heights 5, 7 and 9 m).
 Schedule 1.2 is also proposed for high fertility sites, aiming to provide fuelwood and biomass in a short rotation regime (50 years). The target stand density is 2000 stems ha^{1} with a dominant height of 7 m, a standing volume of 43 m^{3} ha^{1} and aboveground biomass of 51 t ha^{1}. Almost 50% of the extracted timber have an enddiameter over 10 cm, and can be used for fencing. The target stand is attained by applying a single thinning at age of 30 years, when 23 t ha^{1} of dry biomass are extracted for fuelwood.
 Schedule 2.1 is proposed for low quality sites in order to obtain a continuous supply of fuelwood and smaller diameter pieces for fencing. The rotation age is extended to 70 years, but since 30 yearold stands are thinned every 10 years, 58 t ha^{1} of dry biomass are obtained for fuelwood. At rotation age, 1000 stems ha^{1} are felled, with a standing volume of 17 m^{3} ha^{1}, of which almost 44% can be used for fencing and other uses requiring small pieces.
 Schedule 2.2 is based on a pure simple coppice over a 50 year rotation period without thinning, which can be considered the business as usual practice nowadays. Extracted biomass at rotation age is approximately 47 t ha^{1}, its final use being fuelwood.
Apart from thinning regime and rotation age, other relevant decisions include the method of regeneration felling, the proposed silvicultural system and other tending activities. Given the resprouting capacity of the species, we propose a simple coppice method based on vegetative reproduction from resprouts from the stump. The simple coppice method is commonly used today in Tunisia and north Africa ([19]), although the use of coppice with standard methods preserving 4050 trees ha^{1} up to an age of 120150 years has been proposed for Thuya forests ([2], [22]). Due to the huge competition for underground resources (water and nutrients) in the arid environments where Thuya grows, thinnings from below are proposed by selecting the most vigorous shoot(s) within the stump. In order to prevent stump resprouting after thinning and favor allocation of resources to the selected shoots, goat grazing must be promoted ([19]) during this period. Conversely, in order to ensure vegetative regeneration, grazing activity must be avoided after coppice felling during the first 20 years of stand development. Finally, an additional constraint is the role of these stands in preventing soil erosion. On very steep slopes, the coppicing should be carried out in two steps, in order to ensure the maintenance of sufficient plant cover to prevent erosive processes.
It is worth to be noticed that the use of the SDMDs developed in this work is limited by the absence of a selfthinning submodel, since the zone of imminent mortality due to competition ([20]) cannot be represented and future sizedensity lines are not explicitly modeled. Another point to be considered is that when using these diagrams, the data range used for fitting the models should be respected; its extrapolation beyond this range must be assessed with caution. Finally, it should be taken into account that the developed diagrams have been constructed using data from temporary plots, thus temporal projections should be carefully interpreted.
Conclusions
We present a family of site index curves and a set of nonlinear equations for standlevel attributes in Thuya stands in NE Tunisia. The models developed, along with the treelevel models proposed for the species by Calama et al. ([7]) and the diameter increment model proposed by Sghaier et al. ([44]), constitute the first set of basic, flexible tools for the classification and management of these ecologically and economically important forests, which in many arid areas represent the final serial stage of tree forests before shrublands. In this regard, the constructed tools support the forest manager’s decision for making different management objectives (handicraft timber, fencing material, fuelwood, CO_{2} fixation). Future efforts should focus on constructing dynamic climate sensitive models, and implementing them into decisionsupport systems which allow the definition of optimal management schedules under current and changing climate scenarios.
Acknowledgements
This study was carried out within the framework of the bilateral cooperation between Tunisia (INRGREF) and Spain (INIA) (project number: A/017275/08), and was financed by the Ministry for Scientific Research and Technology of Tunisia (MSRT) and the Spanish Agency of the International Cooperation and Development (AECID).
References
::Google Scholar::
::Google Scholar::
::CrossRef::Google Scholar::
::CrossRef::Google Scholar::
::Online::Google Scholar::
::Google Scholar::
::CrossRef::Google Scholar::
::Online::Google Scholar::
::CrossRef::Google Scholar::
::CrossRef::Google Scholar::
::Google Scholar::
::Online::Google Scholar::
::Online::Google Scholar::
::Google Scholar::
::Online::Google Scholar::
::Online::Google Scholar::
::Google Scholar::
::CrossRef::Google Scholar::
::Google Scholar::
::Online::Google Scholar::
::Online::Google Scholar::
::Online::Google Scholar::
::Online::Google Scholar::
::CrossRef::Google Scholar::
::CrossRef::Google Scholar::
::CrossRef::Google Scholar::
::CrossRef::Google Scholar::
::Online::Google Scholar::
::Google Scholar::
::CrossRef::Google Scholar::
::CrossRef::Google Scholar::
::Google Scholar::
::Online::Google Scholar::
::Google Scholar::
::CrossRef::Google Scholar::
::Online::Google Scholar::
::Google Scholar::
::CrossRef::Google Scholar::
::Google Scholar::
::Google Scholar::
::Google Scholar::
::Online::Google Scholar::
::Google Scholar::
::CrossRef::Google Scholar::
::Google Scholar::
::Online::Google Scholar::
::Google Scholar::
::CrossRef::Google Scholar::