Introduction 

Thuya (Tetraclinis articulata (Vahl.) Mast.) is a semiarid Mediterranean forest species which extends discontinuously from North Africa to south-western 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₂ 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 stand-level 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]):

1. 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;
2. stand state is defined by means of dominant height and current stocking density;
3. they do not require a specific sub-model to predict current stand density as a function of dominant height and age (usually the weakest equation in static yield tables);
4. they allow all the information and outcomes to be presented graphically;
5. 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 stand-level 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 sub-model and a system of compatible stand attribute equations presented as a SDMD. In this study, the site index sub-model 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 end-size diameter and aboveground biomass. In a previous study, tree-level 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 sub-model 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 stand-level models are not available. As a subsequent step, the site index sub-model 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 minimum-maximum 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 even-aged 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 clear-cutting, 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 10-12 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.

Fig. 1 - Map showing the studied zone of Tetraclinis articulata in Tunisia.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

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

Volume and biomass estimation

The stem curve taper function and the height-diameter 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 El-Mouridi 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 semi-ellipsoidal 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 north-eastern 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 site-specific parameter X which can neither be measured nor defined ([12]), but depends on soil conditions, and ecological and climatic factors. While GADA assures base-age invariance of the derived algebraic forms, the base-age 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 (M1-M6, see Tab. S1 in Appendix 1) with two site-specific 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 dominant-height growth modeling. Models M1-M2 were derived by Cieszewski ([13]) from the log-logistic model, which is equivalent to the Hossfeld model; model M3 was derived by Krumland & Eng ([29]) from the Bertalanffy-Richards’ function while models M4-M5 were derived from the same base equation by Cieszewski ([14]). Finally, model M6 was derived by Barrio-Anta et al. ([4]) from the Lunqvist-Korf’s function. More information on GADA methodology and dummy approach fitting technique can be found in Gea-Izquierdo et al. ([25]) or Diéguez-Aranda 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 σ² 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 continuous-time autoregressive error structure [CAR(x)], which expands the error terms in the following way ([48] - eqn. 3):

where e_ij is the j-th ordinary residue on the i-th 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 n-order autoregressive parameter to be estimated, and t_ij-t_i(j-n) is the time distance (years) separating the j-th from the j-th-n observations. To evaluate the presence of autocorrelation and the order of the CAR(x) to be used, graphs representing residuals versus lag-residuals 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 self-thinning 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₀) on the x-axis and number of trees per hectare (N) on the y-axis. The main idea behind SDMD is to express all the mathematical relations in the system as a function of N and H₀, 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₀), 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 re-expressed 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₀-a₂, b₀-b₃, c₀-c₃, d₀-d₂ are parameters to be estimated. eqn. 5, 6 and 9 constitute a simultaneous system of equations, where N and H₀ are exogenous variables and Dg, W_T and Vm_i are endogenous variables. Error components of the variables on the left-hand and right-hand side are correlated, so a simultaneous fit of the system of equations was carried out by using the full-information 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₀, 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 lag-residuals 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 second-order continuous-time autoregressive error structure CAR(2), the trends in residuals disappeared (Fig. 2 - third column).

Fig. 2 - Height-Lag1-Residuals and Height-Lag2-Residuals versus Height-Residuals for model M5 fitted without considering the autocorrelation parameters (first column), and using continuous-time autoregressive error structures of first and second order (second and third columns, respectively).

  Enlarge/Shrink   Download   Full Width  Open in Viewer

The parameter estimates for each model and their corresponding goodness-of-fit statistics are shown in Tab. 2. The two models M3 and M6 present a non significant parameter (b₃ for M3 and b₂ for M6). For the other four models, all the parameters (including the site-specific parameters for each tree) were highly significant, although model M4 showed worse performance in terms of RMSE, R²_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.

Thus, among the six models evaluated for height growth prediction and site classification of Tetraclinis articulata stands in north-eastern Tunisia, Model M5 (the GADA formulation derived from the Bertalanffy-Richards base function by considering a₁ and a₃ 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₀ and t₀ represent the predictor dominant height and age. Note that the site-specific 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 40-65 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 mid-interval curves for site indexes at 4 and 6 m.

Fig. 3 - Curves for site index of 3, 5 and 7 m at the reference age of 40 years overlaid on the trajectories of the observed heights over time for models M5.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

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.

Parameter estimates were used to develop SDMDs by superimposing the isolines for the different attributes on a bivariate axis defined by H₀ 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 200-5000 trees ha^-1 for stand density, 2-12 m for dominant height, 5-120 m³ ha^-1 for stand volume, 5-140 tons of dry matter ha^-1 for aboveground biomass and 0.05-0.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 self-thinning 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.

Fig. 4 - Stand Density Management Diagram SDMD for Thuya: (a) total volume and mean squared diameter; (b): aboveground biomass and mean squared diameter; (c): total volume and ratio of merchantable volume (timber volume with and end diameter over 15 cm).

  Enlarge/Shrink   Download   Full Width  Open in Viewer

  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 north-eastern 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 Bertalanffy-Richards’ 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 age-dominant height pair without compromising the consistency of predictions, and the predictions are unaffected by changes in the base-age. 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 stand-level 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 first-draft approach, in the absence of more accurate data, it can be considered adequate for estimating and fitting stand-level 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:

1. 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).
2. 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.
3. 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 with SDI varying between 250 and 500, avoiding higher stocking levels in order to prevent losses due to self-thinning.
4. 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 above-mentioned constraints, we propose to evaluate four different management schedules (Tab. 4) for mid-high (SI = 6 m) and mid-low (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.

1. 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 year-old stand with 1000 stems ha^-1, a dominant height of 10.5 m, standing volume of about 47 m³ ha^-1, and 30% total volume (15 m³ ha^-1) of timber with an end-size 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).
2. 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³ ha^-1 and aboveground biomass of 51 t ha^-1. Almost 50% of the extracted timber have an end-diameter 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.
3. 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 year-old stands are thinned every 10 years, 5-8 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³ ha^-1, of which almost 44% can be used for fencing and other uses requiring small pieces.
4. 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 40-50 trees ha^-1 up to an age of 120-150 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 self-thinning sub-model, since the zone of imminent mortality due to competition ([20]) cannot be represented and future size-density 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 stand-level attributes in Thuya stands in NE Tunisia. The models developed, along with the tree-level 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₂ fixation). Future efforts should focus on constructing dynamic climate sensitive models, and implementing them into decision-support 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 co-operation 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

  Supplementary Material

Fig. S1 - Bias and root mean square error (RMSE) in height estimation for site index models M1, M2 and M5 by 5-year age classes.
Fig. S2 - Relative error in height prediction (RE) related to choice of reference age for Model (M5) by 5-year age classes adjusted with CAR(2) and number of observation (n).
Fig. S3 - Site index predictions against age using model M5 and the stem analysis data.
Tab. S1 - Base models and GADA formulations considered.
Tab. S2 - Proposed stand density isolines for constructing SDMD.