Thuya (Tetraclinis articulata (Vahl) Mast) is a Mediterranean forest species mainly occupying semiarid environments in North African countries, where it provides important ecological and economical services, such as biodiversity conservation, soil protection against erosion, fuelwood, timber for fencing, construction and handicraft, resins, etc. Despite the importance of the species, there is a severe lack of scientific knowledge as regards the management of these forests or modeling tools to support multifunctional forest management decision making. In the present work, we developed a stand-level integrated model for the management of Thuya forests in Tunisia. The model comprises a family of site index curves, built using the Generalized Algebraic Difference Approach (GADA) method, which provides predictions for stand growth, aboveground biomass, total and merchantable volumes, along with a non-linear system of stand level equations presented as stand density management diagrams (SDMD). The developed model has been used to define, characterize and compare four different management specific schedules for different site qualities and multifunctional objectives.
Site Index ModelStand Density Management DiagramSemiarid EnvironmentsStand-level ModelManagement ScheduleThinningModellingIntroduction
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 (Charco 1999). In Tunisia, the geographical distribution of the Thuya is limited to the areas close to the Cap Bon, Boukornine and Zaghouan (Rejeb et al. 1996). The estimated total area of Thuya in Tunisia is approximately 33 000 hectares, with more than 22 000 ha of homogeneous stands (DGF 1995). 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 (Tenbergen et al. 1995). 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 (Seigue 1985). 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 (Esteve-Selma et al. 2010).
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 (Charco 1999, El-Mouridi 2011), 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 (DREF 2002), the production and income per hectare being very low. These unfavorable conditions can be aggravated under scenarios of climate change (El-Mouridi et al. 2011).
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 (Nabil 1989, Schoenenberger 1995, Nicolás et al. 2004, Esteve-Selma et al. 2010). However, in recent years a number of papers focusing on growth and yield modeling have been published (Calama et al. 2012, Sghaier et al. 2013).
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 (Carmean 1975, Cieszewski & Bailey 2000), 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 (Newton & Weetman 1994). 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 (Barrio-Anta & Álvarez-González 2005, Castedo-Dorado et al. 2009):
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 sub-model 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 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 (Calama et al. 2012) and more recently, a distance independent diameter increment model has been constructed (Sghaier et al. 2013). 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 methodsStudy 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 (Ben Mansoura & Garchi 2001). 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 - DGF 1995). 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.
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 (Carmean 1972), and the modification proposed by Newberry (1991) for the topmost section of the tree, based on earlier studies (Dyer & Bailey 1987). 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. (2012) 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. (2011). 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 (Scott & Reinhardt 2001). 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. (2012). 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 (Ricardi et al. 2007), 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 (2000) and using the dummy approach to make the adjustment (Cieszewski et al. 2000).
Productivity is assumed to be dependent on an unobservable site-specific parameter X which can neither be measured nor defined (Cieszewski & Bailey 2000), 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 (Bailey & Clutter 1974) was assured by fitting them using the dummy variable approach (Cieszewski et al. 2000). 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 (2002) from the log-logistic model, which is equivalent to the Hossfeld model; model M3 was derived by Krumland & Eng (2005) from the Bertalanffy-Richards’ function while models M4-M5 were derived from the same base equation by Cieszewski (2004). Finally, model M6 was derived by Barrio-Anta et al. (2006) from the Lunqvist-Korf’s function. More information on GADA methodology and dummy approach fitting technique can be found in Gea-Izquierdo et al. (2008) or Diéguez-Aranda et al. (2005).
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):
AICd = n \log {\hat{ \sigma}}^2 + 2(p+1) - min \left [n \log {\hat{ \sigma}}^2 + 2(p+1) \right ]
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 (Huang et al. 2003).
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 (Gregoire et al. 1995), 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 (Zimmerman et al. 2001 - 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 (SAS Institute Inc 2004), 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 (Reineke 1933), relative spacing index (Wilson 1946) or the self-thinning rule (Yoda et al. 1963). 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 = N \left (\frac{25} {Dg} \right )^{-1.605}
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 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_{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 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_{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 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 Castedo-Dorado et al. 2009 for more details in fitting techniques). A simultaneous fit of the models was carried out by using again the SAS/ETS^{®} MODEL procedure (SAS Institute Inc 2004).
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 (El-Mouridi 2011).
ResultsSite 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).
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_{3} for M3 and b_{2} 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^{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 (Ruiz de la Torre & Ceballos 1979). 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_{1} and a_{3} as related to site productivity) was selected (eqn. 10):
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 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 (Diéguez-Aranda et al. 2005). 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 (Diéguez-Aranda et al. 2005). 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.
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_{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 200-5000 trees ha^{-1} for stand density, 2-12 m for dominant height, 5-120 m^{3} 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.
DiscussionFitted 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 (Cieszewski 2002, Diéguez-Aranda et al. 2005), 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 (Assman 1970, Clutter et al. 1983). 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 (Barrio-Anta & Álvarez-González 2005, Castedo-Dorado et al. 2009, López-Sánchez & Rodríguez Soalleiro 2009) 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 (El-Mouridi et al. 2011, Ricardi et al. 2007) 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:
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 (1985), 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.
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 (El-Mouridi 2011). 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 (Haddad et al. 2005). 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.
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^{3} ha^{-1}, and 30% total volume (15 m^{3} 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).
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 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.
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^{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 (DREF 2002), 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 (Bachoua & Voreux 1987, El-Mouridi 2011). 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 (DREF 2002) 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 (Drew & Flewelling 1979) 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. (2012) and the diameter increment model proposed by Sghaier et al. (2013), 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 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).
ReferencesAssman EThe principles of forest yield studies. Pergamon Press Ltd, Oxford, UK, pp. 506.1970Bachoua A, Voreux ChL’aménagement de la Tétraclinaie de l’Amsittène (Maroc) [The management of the Amsittene Tetraclins forest (Morocco)]. Forêt Privée 174: 57-68. [in French]1987Barrio-Anta M, Álvarez-González JGDevelopment of a stand density management diagram for even-aged pedunculate oak stands and its use in designing thinning schedules. Forestry 78: 209-216.2005Barrio-Anta M, Castedo Dorado F, Diéguez-Aranda U, Álvarez-González JG, Parresol BR, Rodríguez-Soalleiro RDevelopment of a basal area growth system for maritime pine in northwestern Spain using the generalized algebraic difference approach. Canadian Journal of Forest Research 36 (6): 1461-1474.2006Bailey RL, Clutter JLBase-age invariant polymorphic site curves. Forest Science 20: 155-159.1974Ben Mansoura A, Garchi SCaractérisation de la croissance et de la régénération du Thuya par une technique modifiée de mesure de distances [Characterization of the growth and the regeneration of the Thuya by a modified distances measurement technique]. Annales de l’INRGREF, Numéro spécial 2001, pp. 54-76. [in French]2001Calama R, Sánchez-González M, Garchi S, Ammari Y, Cañellas I, Sghaier TTowards the sustainable management of thuya (Tetraclinis articulata (Vahl.) Mast.) forests in Tunisia: models for main tree attributes. Forest Systems 21 (2): 210-217.2012Carmean WHSite index curves for upland oaks in the Central States. Forest Science 18: 109-120.1972Carmean WHForest site quality evaluation in the United States. Advances in Agronomy 27: 209-269.1975Castedo-Dorado F, Crecente-Campo F, Álvarez-Álvarez P, Barrio-Anta MDevelopment of a stand density management diagram for radiate pine stands including assessment of stand stability. Forestry 82 (1): 1-16.2009Charco JEl bosque mediterráneo en el norte de África: biodiversidad y lucha contra la desertificación [The Mediterranean forest in North Africa: Biodiversity and desertification]. Ed. Mundo Árabe e Islam - AECI, Madrid, Spain, pp. 370.1999Cieszewski CJ, Bailey RLGeneralized algebraic difference approach: theory based derivation of dynamic site equations with polymorphism and variable asymptotes. Forest Science 46: 116-126.2000Cieszewski CJComparing fixed and variable base-age site equations having single versus multiple asymptotes. Forest Science 48: 7-23.2002Cieszewski CJGADA derivation of dynamic site equations with polymorphism and variable asymptotes from Richards, Weibull, and other exponential functions. PMRC-TR 2004-5, University of Georgia, Athens, GA, USA, pp. 16.2004Cieszewski CJ, Harrison M, Martin SWPractical methods for estimating non-biased parameters in self-referencing growth and yield models. PMRC-TR 2000-7, University of Georgia, Athens, GA, USA, pp. 11.2000Clutter JL, Fortson JC, Pienaar LV, Brister GH, Bailey RLTimber management: a quantitative approach. John Wiley and Sons Inc., New York, USA, pp. 333.1983DGFRésultats du premier inventaire forestier national en Tunisie. [Results of the first national forest inventory in Tunisia]. Direction Générale des Forêts, Ministère de l’agriculture, Tunisie, pp. 88. [in French]1995Diéguez-Aranda U, Burkhart HE, Rodríguez-Soalleiro RModeling dominant height growth of radiata pine (Pinus radiata D. Don) plantations in north-western Spain. Forest Ecology and Management 215 (1-3): 271-284.2005DREFThuya: importance écologique et économique. [Thuya: Ecological and economic importance]. Terre et Vie 52: 4.2002Drew TJ, Flewelling JWStand density management: an alternative approach and its application to Douglas-fir plantations. Forest Science 25: 518-532.1979Dyer ME, Bailey RLA test of six methods for estimating true heights from stem analysis data. Forest Science 33: 3-13.1987El-Mouridi MCaractérisation mécanique de la loupe de thuya (Tetraclinis articulata (Vahl.) Mast) en vue de sa valorisation [Mechanical characterization of the magnifying glass of Thuya (Tetraclinis articulata (Vahl.) Mast) for its valorization]. PhD Thesis. Université Mohamed V-Agadal, Faculté des Sciences, Rabat, Morocco, pp. 119. [in French]2011El-Mouridi M, Laurent T, Famiri A, Kabouchi B, Alméras T, Calchéra G, El Abid A, Ziani M, Gril J, Hakam ACaractérisation physique du bois de la loupe de thuya (Tetraclinis articulata (Vahl) Masters) [Physical characterization of the magnifying glass wood of the Thuya (Tetraclinis articulata (Vahl) Masters)]. Physical Chemical News 59: 57-64. [in French]2011Esteve-Selma J, Martinez-Fernandez J, Hernández I, Montávez JP, Lopez JJ, Calvo JF, Robledano FEffects of climatic change on the distribution and conservation of Mediterranean forests: the case of Tetraclinis articulata in the Iberian Peninsula. Biodiversity and Conservation 19: 3809-3825.2010Gea-Izquierdo G, Cañellas I, Montero GSite index in agroforestry system: age-dependent and age-independent dynamic diameter growth models for Quercus ilex in Iberian open oak woodlands. Canadian Journal of Forest Research 38 (1): 101-113.2008Gregoire TG, Schabenberger O, Barrett JPLinear modelling of irregularly spaced, unbalanced, longitudinal data from permanent-plot measurements. Canadian Journal of Forest Research 25: 137-156.1995Haddad A, Lachenal D, Marechal M, Kaid-Harche M, Janin GCaractéristiques papetières de la p’te de bois de thuya de Berbérie (Algérie) (Tetraclinis articulata Vahl) obtenue par un procédé soude-anthraquinone [Paper characteristics of the wood pulp of Thuya (Algeria) (Tetraclinis articulata Vahl) obtained by a soda-anthraquinone process]. Annals of Forest Science 63: 493-498. [in French]2005Huang S, Yang Y, Wang YA critical look at procedures for validating growth and yield models. In: “Modelling forest systems” (Amaro A, Reed D, Soares P eds). CAB International, Wallingford, Oxfordshire, UK, pp. 271-293.2003Krumland B, Eng HSite Index systems for major young-growth forest and woodland species in northern California. California Forestry Report 4, California Department of Forestry and Fire Protection, Sacramento, CA, USA, pp. 219.2005Long JNA practical approach to density management. Forestry Chronicle 23: 23-26.1985López-Sánchez C, Rodríguez Soalleiro RA density management diagram including stand stability and crown fire risk for Pseudotsuga menziesii (Mirb.) Franco in Spain. Mountain Research and Development 29 (2): 169-176.2009Nabil MAEssai de synthèse sur la végétation et la phyto-écologie tunisiennes. I - Elément de botanique et de phyto-écologie [Synthesis assay on the Tunisian vegetation and phytoecology. I - Botany and phytoecology element] (vol. 4-6). Faculté des Sciences, Tunis, Tunisia, pp. 247. [in French]1989Newberry JDA note on carmean’s estimate of height from stem analysis data. Forest Science 37 (1): 368-369.1991Nicolás MJ, Esteve MA, Palazón JA, López Hernández JJModelo sobre las preferencias de hábitat a escala local de Tetraclinis articulata (Vahl) Masters en una población del límite septentrional de su área de distribución [Modelling the habitat preferences at local scale of Tetraclinis articulata (Vahl) Masters of a population in the Northern boundary of its distribution area]. Anales de Biología 26: 157-167. [in Spanish]2004Newton PF, Weetman GFStand density management diagram for managed black spruce stands. Forestry chronicle 70: 65-74.1994Reineke LHPerfecting a stand density index for even-aged forest. Journal of Agricultural Research 46: 627-638.1933Rejeb MN, Khaldi A, Khouja ML, Garchi S, Ben Mansoura A, Nouri MGuide pour le choix des espèces de reboisement: espèces forestières et pastorales [Guide for the reforestation species choice: forest and pastoral species]. INGREF, Ministère de l’agriculture, Tunisie, pp. 137. [in French]1996Ricardi CL, Prichard SJ, Sandberg DV, Ottmar RDQuantifying physical characteristics of wildland fuels using the Fuel Characteristic Classification System (FCCS). Canadian Journal of Forest Research 37: 2413-2420.2007Ruiz de la Torre J, Ceballos LÁrboles y arbustos de la España Peninsular [Trees and shrubs of peninsular Spain]. ETS Ingenieros de Montes, Madrid, Spain, pp. 400. [in Spanish]1979SAS Institute IncSAS/ETS® 9.1 User’s Guide. SAS Institute Inc, Cary, NC, USA, pp. 2416.2004Schoenenberger AGroupements végétaux des séries de végétation. La végétation forestière [Plant groupings of the series of vegetation. Forest vegetation]. In: “Essai de synthèse sur la végétation et la phyto-écologie tunisiennes (vol. II&III). Programme flore et végétation tunisiennes”, pp. 213-237. [in French]1995Scott JH, Reinhardt EDAssessing crown fire potential by linking models of surface and crown fire behavior. Researh Paper RMRS-RP-29, Rocky Mountain Research Station, USDA Forest Service, Fort Collins, CO, USA, pp. 59.2001Seigue ALa forêt circumméditerranéenne et ses problèmes [The circum-Mediterranean forest and its problems]. In: “Techniques agricoles et productions méditerranéennes (vol. 5)” (Coste R ed). Editions Maisonneuve et Larose, Paris, France, pp. 230-232. [in French]1985Sghaier T, Tomé M, Tomé J, Sánchez-González M, Cañellas I, Calama RDistance-independent individual tree diameter-increment model for Thuya (Tetraclinis articulata (VAHL.) MAST.) stands in Tunisia. Forest systems 22 (3): 433-441.2013Tenbergen B, Günster A, Schreiber KFHarvesting runoff: the minicatchment technique - an alternative to irrigated tree plantations in semiarid regions. Ambio 24: 72-76.1995Wilson FGNumerical expression of stocking in terms of height. Journal of Forestry 44: 758-761.1946Yoda K, Kira H, Ogawa H, Hozumi KSelf-thinning in overcrowded pure stands under cultivate and natural conditions. (Intraspecific competition among higher plants XI). Journal of the Institute of Polytechnics, Osaka City University, Series D 14: 107-129.1963Zimmerman DL, Núñez-Antón V, Gregoire TG, Schabenberger O, Hart JD, Kenward MG, Molenberghs G, Verbeke G, Pourahmadi M, Vieu PParametric modelling of growth curve data: an overview. Test 10 (1): 1-73.2001
Map showing the studied zone of Tetraclinis articulata in Tunisia.
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).
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.
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).
Summary statistics for the used data set. (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 (Reineke 1933).
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
Parameter estimates and goodness of fit statistics for site index models M1-M6. (Par): Parameter; (Est): Estimate; (SE): Standard Error; (H_{250}): dominant height in m for class 1 at 250 years as an estimate of maximum growth rather than the asymptote.
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.48E-14
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.04E-48
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.16E-3
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.857E-3
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.903E-3
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
Parameter estimates and goodness-of-fit statistics for SDMD system of nonlinear equations. (Par): Parameter; (Est): estimate; (SE): standard error; (RMSE): root mean square error; (p-value): probability associated with the hull hypothesis (bias = 0).
Eqn.
Responsevariable
Par
Est
SE
RMSE
Bias
R^{2}_{adj}
p-value
[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.98E-6
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.2E-5
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
Proposed silvicultural schedules for Thuya stands in NE Tunisia. Age (years); (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 15-10 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
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.