STAND DENSITY MANAGEMENT DIAGRAMS OF Eucalyptus viminalis: PREDICTING STEM VOLUME, BIOMASS AND CANOPY COVER FOR DIFFERENT PRODUCTION PURPOSES

Stand density management diagrams (SDMD) provide a guide for forest density management taking into account stands attributes such as trees ́ diameter or volume. One of the most common species planted in Pampean plains of Argentina is Eucalyptus viminalis for multiple objectives: solid wood use or firewood in local markets, pulp for cellulose industry and to provide services for agriculture and cattle raising (windbreaks or cattle refuge). The objective of this study was to gather the available production information /inventory data and to develop a first SDMD for estimating standing volume, biomass and canopy cover of E. viminalis as a tool for forest managers aiming at different plantation purposes. Data to develop the SDMD were obtained from 161 plots, distributed along a climate and soil gradient. We also generated two predictive equations capable of estimating dominant height from the diameter of the trees as well as canopy cover from stand basal area. As an example of application, the SDMD was used to estimate the wood production of three alternative systems: a) an unmanaged plantation (simulating a common practice in the region), b) a mixed production system such as an agroforestry system, and c) a plantation that maximizes wood biomass or volume production. v.25 n.4 2019 STAND DENSITY MANAGEMENT DIAGRAMS OF EUCALYPTUS VIMINALIS: PREDICTING STEM VOLUME, BIOMASS AND CANOPY COVER FOR DIFFERENT PRODUCTION PURPOSES


INTRODUCTION
Several attributes of forests are related to the stocking of trees per unit land area (i.e. tree density), such as mortality or growth rate (current or mean annual increment), average size of the individuals (i.e. stem diameter) and other characteristics of the stand (i.e. canopy cover, carbon accumulation; Long and Vacchiano, 2014). These authors emphasize the relevance of estimating the boundary line that maximizes the combination of these attributes -density relationships.
In the case of even-aged stands, intraspecific competition increases following the increase in tree size reaching a self-thinning situation. In this case, the increase in tree dimension depends on the decrease of the number of trees per unit area (Long et al., 2004). This kind of attribute -density relationship allows for foresters to develop stand density management diagrams (SDMD). These are average stand-level models that had been developed for a large number of species growing as pure or mixed stands, increasing their applicability following several objectives (Newton, 1997). The basic purpose of a SDMD is to easily guide the management of forest density taking into account stands attributes such as wood volume, canopy cover or biomass production, and also, the amount of energy provided by biomass or the potential pulp production (i.e. Pérez et al., 2008;Pérez-Cruzado et al., 2011;Tang et al., 2016).
Although these models are widely used worldwide, their development and implementation require large amounts of information, making them poorly available or even non-existent in some regions, particularly for species that have not been traditionally used for intensive wood production. This is the case of Pampa region of Argentina, South America, where agriculture and cattle raising are the main primary production activities. In this region, eucalypt silviculture is not as technologically advanced compared to intensive agriculture or to eucalypt silviculture in other regions of the country (i.e. North Eastern region) with a welldeveloped forestry sector. In this regard, there are no SDMD for anyof the forestry species cultivated in the Pampa region. The most common species planted in this region is Eucalyptus viminalis Labill (SINAVIMO, 2018). This species is widely distributed in Southeastern Australia, its provenance region, and is planted also in several countries of Africa showing a broad capacity to grow under different ecological constraints. As in its origin, in Pampa plains of Argentina it is cultivated under different soil and climate conditions. In this region, E. viminalis is grown between 600 and up to 1000 mm of mean annual precipitation. The highest growth rate occurs at the center of the region, reaching mean stem volume productivity of 35 to 40 m 3. ha -1. y -1 , decreasing this variable to the West due to water deficit, and to the East, because of the increase of waterlogging (SINAVIMO, 2018) associated with fine textured soils. In average, annual productivity is around 25-30 m 3. ha -1. y -1 , reaching 250-300 m 3 . ha -1 of merchantable wood at the end of the rotation period (around 10-12 years; BA Buenos Aires Forestal 2010). Currently, the main destination of E. viminalis wood is for pulp use in a quite geographically distant cellulose industry due to its good kraft pulp quality (Neiva et al., 2014), and in a minor degree, for solid use in local markets of products of low industrialization value. However, the species has a large potential for the production of veneer sheets and plywood panels (Iwakiri et al., 2013) and for direct structural solid use, increasing its economic value. Since the the regional economy is based on agriculture and cattle raising, it is common that forest stands are small, mainly smaller than 1 ha, with main or complementary objectives of providing cattle refuge, acting as windbreaks or providing firewood for the local residents (i.e. Mónaco and Laclau, 2013). At the same time, the potential for the development of silvopastoral systems as well as other agroforestry systems is currently being highlighted and promoted by different institutions (e.g. INTA, IPCVA, AFoA ) due to the need of innovative production systems that are more environmentally friendly, economically diversified and with a higher capacity of climate change mitigation than current agroecosystems.
Due to its high growth potential and adaptability, several trials have been installed analyzing differences in survival and growth of E. viminalis provenances from Australia and local Argentinean seed lots in different quality sites (Mendoza, 1974;Cappa et al., 2010). In addition, from samples of trees growing in a similar range of ecological situations, Ferrere et al. (2014) have developed allometric models to estimate biomass of the leaves, branches and stems (with and without bark) of this species. In spite of this knowledge, there does not exist yet a compilation of this information in a SDMD type model to be available for both farmers and wood markets, allowing a rapid estimation of standing volume from inventory data (Schnell et al., 2012). Based on this, the objective of the present study was to gather the available production information and to develop a first stand density diagram for estimating standing volume, biomass and canopy cover of E. viminalis planted in the Pampas plains of Argentina.

Study area
The study was conducted in the Pampas plains, Buenos Aires Province, Argentina ( Figure 1). The climate of the region is humid temperate, with average annual temperatures of 14-15°C , and annual average precipitation from 1000 mm in the East to less than 600 mm . per . year -1 in the West, evenly distributed throughout the year (Paruelo et al., 2007;Barros et al., 2013). This causes a water balance (precipitation minus evapotranspiration) gradient from positive or near zero in the NE area, to negative ( -400 mm per year -1 ) in the SW area ( Figure 1). Regarding the soils, changes of depth and texture are observed, on an East-West direction, from relatively shallow loam clayey soils due to the presence of petrocalcic horizons (locally called "tosca") or rock, to sandy deep soils in the West (Paruelo et al., 2007).

SDMD construction
The data used to develop the SDMD were obtained from 161 plots, all established in Buenos Aires province ( Figure 1). Some of these data are unpublished (37 plots), but other come from the original database that supports the following publications: Caccia et al. (1991; n=14 plots), Ferrere et al. (2008;n= 22 plots), Ferrere et al. (2014; n=72 plots), Milione (2019;n=16 plots). Because some of these data were not digitalized before this study or the published articles or reports are written in Spanish, we include the data in the Appendix A. This database covers the main distribution of E. viminalis in Buenos Aires province and most of Pampa region: the northern plots were located near Paraná river (33°35'31.20"S -59°52'3.58"O) meanwhile the Southern and Eastern plots were placed near Tandil city (37°38'37.44"S -59° 4'27.10"O). The Western-most plots were located near Gral. Villegas city (35°03´57.9"S -62°52´02" O) ( Figure 1).
Quadratic mean diameter of each plot (Dq, cm), dominant height (H, m, the height of the trees with highest Dg of the plot; criteria adopted by the authors of the databases Ferrere and Milione; see Appendix A) and number of trees per hectare (N, trees . ha -1 ) were used as the basic variables to develop the SDMD. We use Dg and N data to develop SDMD of Caccia et al. (1991), Ferrere et al. (2008, both  The measured stands covered a broad range of tree density (from 127 to 3680 tree . ha -1 ) and sizes of the trees (mean Dg= 10 to 56 cm; mean H= 8 to 33 m; Table 1). Some of the stands were harvested several times, meanwhile other had never received any intervention. From this set of stands, we followed the methodology proposed by Scharf et al. (1998) andZhang et al. (2005) in order to objectively select those points situated in the upper boundary of Dg and N relationship for all available plots. For this purpose, three methods were applied. First, we manually chose those points considered as the limit of the relationship. As a second method, we divided the range of Dg into equal intervals (we chose 6 intervals) and selected the plot with the maximum N from each interval. Finally, as a third method, we applied a quantile regression (quantile 0.95). We selected 18 plots using the first method, 6 following the second method, and 7 for the quantile regression (third method). Some of these plots were selected by two or three methods applied. Finally, we adjusted the following equation to each dataset obtained from the applied methods. Where: N= tree density, trees .

FIGURE 1 Study region showing the sites of the different
Eucalyptus viminalis stands used to develop the SDMD (empty and filled dots). The climatic water balance (mm y -1 ) was estimated as the difference between mean annual precipitations minus Penman-Monteith potential evapotranspiration (based on Nosetto et al., 2008). Empty dots represent the plots used to estimate the relationship between basal area and canopy cover (database from Milione, 2019; Figure 4).

CERNE
GYENGE et al ha -1 , Dg= quadratic mean diameter, cm, b represent the exponent of Reineke´s equation (equation 2), a= y-intercept of the linear equation.
1.48v program (Wayne Rasband, National Institutes of Health, USA). We used at least 5 photos per plot. The relationship between forest canopy cover (COV) and stand basal area (G) and Dg was computed by means of a regression analysis. Stand basal area (G) showed the best adjustment following (Table 2). Where: f and g correspond to the y-intercept and the slope, respectively, of the linear equation 4. [1] The adjusted model 1 using three methodologies of data selection were compared using F tests with an α = 0.05 (Neter and Wasserman, 1974). Stand density index (SDI; Reineke, 1933) was estimated for each plot according to the following formula 2. Where: b is the slope of equation 1, 25 is the reference quadratic mean diameter, in cm.

[2]
A dominant height -quadratic mean diameter relationship (H -Dg) was developed using data from 101 plots from the database provided by Ferrere et al. (2014), Ferrere (unpublished data, n=85 plots), and Milione (2019, n=16 plots). This dataset was divided into a modelling dataset and a validation dataset (50% of the data on each). We adjusted the following model (i.e. Tang et al., 2016) taking into account that parameter c was adjusted to the data and not using a standard value of 1.3: For validation purposes, the goodness of fit and the prediction capacity of the models were evaluated through the coefficient of determination (R 2 ), the observed versus predicted plot (with the line y = ŷ marked), the mean absolute percent error (MA%E=100 [Σ‫|‬yi−ŷi‫|‬ ‫|‬yi‫|‬ -1 ] n-1 , where yi and ŷi are the measured and estimated values), and a paired t test (Mayer and Butler, 1993).
Forest canopy cover (COV, %) was estimated from the analysis of flat photographs taken in the 14 plots described in Milione (2019) by means of ImageJ Application of SDMD for developing thinning schedules and yield estimation As an example of use of the developed SDMD, we propose three different thinning schedules, based on tree density or canopy cover, each attempting to optimize different variables depending on different production objectives (Figure 2): A. Maximization of stem biomass production: thinning of 50% of the trees when the SDI reaches a value of 80% of the maximum SDI; B. Development of a mixed production system requiring solar radiation reaching the understory: managing tree density in order to maintain the canopy cover (COV) between 40 and 50%; and C. No intervention (control situation). We considered two types of intervention: in the model A, we "harvested" the trees with the lower Dg reserving large trees (thinning from below). For this reason, the Dg of the stand increased following the trajectory of H (equation 3) because this last parameter is insensitive to stand density (i.e. Valbuena et al., 2008;Tang et al., 2016). On the contrary, in the model B, the objective of the management was to maintain light reaching the understory with a regular spatial distribution of the trees (systematic thinning). In this case, the thinning implied a spatially regular intervention, which resulted in no change in the average Dg of the stand (Figure 2). Finally, we assumed no natural mortality occurred between thinning treatments so size density trajectory was drawn parallel to x-axis. A similar assumption was also adopted to develop SDMD for other species (Valbuena et al., 2008).

RESULTS
Non-significant differences were found between the linear models describing the maximum relationship between lnN and lnDg (eq [1]) using the three different data selection methods. Based on this, one model was adjusted using a combination of the data sets (n=24).  Table 2).
The parameters of the relationship between H and Dg (equation 3) were estimated with a high level of statistical significance (p < 0.0001) presenting a relatively good goodness of fit (R 2 = 0.796, Table 2). The goodness of fit between the predicted and observed volume values was also high (R 2 = 0.802). No bias was identified in the relationship between the observed H values (the validation dataset) and the estimated H using the equation 3 (Figure 3). MA%E was low (8.36%) and paired t test indicated no statistical differences between observed and estimated data.
Canopy cover increased following a nonlinear relationship from approximately 30% in those stands with basal areas lower than 20 m 2 . ha -1 to reach maximum values around 60%, when G was around 50 m 2. ha -1 , in spite of the different environmental conditions (Figure 4; Table 2, equation 4, R 2 = 0.802, p<0.001).

FIGURE 2 Stand density management diagram (SDMD) for
Eucalyptus viminalis in Pampas plains, Argentina, including three examples of different types of thinning schedules (model A: wood production optimization, thinning from below; B: mixed production system; systematic thinning; and C: unmanaged plantation, control) planned based on different production objectives. SDI = stand density index (SDI).

FIGURE 3
Relationship between observed and estimated H using equation 3.

FIGURE 4
Relationship between basal area (G, m 2 . ha -1 ) and stand canopy cover (COV, %) of Eucalyptus viminalis in the Pampas plains described in Milione (2019). Different symbols denote the different sites corresponding to the empty dots in Figure 1.

SDMD Construction
The first step in the construction of the SDMD was to fit the non-linear equations described in Table 2, relating stem volume without bark (VOL, m 3 . ha -1 ), stem biomass (BIOM, kg . tree -1 ) and canopy cover (COV, %) ( Table 2). Taking into account the preceding equations (

Modelled relationships between stand variables
The slope of the relationship between N and Dg was relatively high (less negative) compared with other species (slopes ranged between -1.63 and -2.00, Vospernik and Sterva 2015; or -1.47 to -2.24 in Condés et al., 2017), even compared to the value adopted as a universal constant (-1.605) estimated by Reineke (1933). Less steeped slopes of this relationship had been estimated for stands with a skewed diameter distribution (uneven-aged stands, Gül et al. 2005). In the present study, only stands of highest plantation density showed a skewed diameter distribution to the lowest diameter values (data not shown), probably influencing the observed slope value. In this regard, it has been proposed that fast-growing Eucalypts plantations present growth dominance after canopy closure, possible linked to their physiological plasticity (changes in resource use and resource use efficiency; Whitehead and Beadle, 2004;Fernández et al., 2011), which could enhance size differences between the trees. Added to this, most of the studied plots lacked management, favouring high levels of intraspecific competition for resources and the expression of growth dominance. In this regard, Binkley et al. (2003) pointed out that in E. saligna monoculture, the largest 20% of the trees accounted for about half of the current wood growth. This marked growth dominance could result in skewed size distribution; however more effort is needed analysing the size distribution of the studied species under different environmental conditions and silvicultural management.
The strong relationship between COV and G in spite of the environmental differences between sites, at least in the annual water balance (Figure 1), was similar with those described for E. grandis plantations in South Africa. In this case, after canopy closure, the leaf area index in all density plantation density treatments converged to a similar level, that showed a slightly change in response to the soil water availability imposed by climate conditions (du Toit, 2008). Based on this, the silvicultural management in a mixed production system, such as silvopastoril systems, could be easier than that developed with species reaching higher canopy cover due to their intrinsic higher shade tolerance. In this regard, the maximum leaf area index (leaf m 2 . per . ground . m 2 ) developed by Eucalyptus plantations under the best conditions of soil water availability and temperature reached 6 (Battaglia et al., 1998) against 12 in Douglas fir (e.g. Gyenge et al., 2009). With the exception of very young or very dense stands, the incomplete canopy cover of Eucalyptus plantations could be based on their intolerance to low irradiance and/or the fragility of their foliage (Schönau and Coetzee, 1989). In general, due to their leaf angles (60-80° from the horizontal plane) and highly clumped leaf area, solar interception is lower in Eucalyptus spp forests than in other forests of similar height such as coniferous forests (Whitehead and Beadle, 2004).
Estimated maximum stand density index (SDI) was 1306, which represents the maximum number of stems per hectare for a reference diameter of 25 cm (equation 2). This value was close to 1210, the value estimated for plantations of E. globulus (Reineke, 1933). We chose threshold values of 60 and 40 of SDI based on published information. In general, based on Long (1985), it is considered that 60% of maximum SDI is the lowest threshold of self-thinning, while a range around 25-35% of the maximum SDI represents the onset of competition and/or the full site occupancy. In the particular case of E. regnans, the lower limit to self-thinning was around 55%, whereas 38% of maximum SDI corresponded to the maximum limit in which dominant trees showed a growth response to thinning (Goodwin, 1990). For E. grandis in Brazil, those limits corresponded with 55 and 45% of the maximum SDI, respectively (Marangon et al., 2017). However, because we do not have information about the mortality rate in E. viminalis, during the simulations we did not consider natural mortality between two thinning interventions, as indicated above. Full site occupancy and onset of competition for the development of management recommendations of E. viminalis from the SDMDs, are indicated as black filled lines in Figure 5.
Maximum SDI could change depending on environmental conditions and/or the particular responses to them of each species (Condés et al., 2017). These authors found that meanwhile the intercept and slope (in absolute value) in Pinus sylvestris were higher in the less humid areas compared to wetter ones, both parameters decreased in the more xeric sites in Fagus sylvatica. So in a particular environmental gradient, different species could present patterns of increasing or decreasing the parameters of the relationship between number-size of the trees. It is possible that maximum SDI may change in the Pampas plains of Buenos Aires, basically due to the different soil texture (from Clay to Sandy soils -e.g. Milione, 2019) and also, due to the gradient in water balance (Figure 1) that limits the productivity of the stands (BA Buenos Aires Forestal, 2010, SINAVIMO, 2018. However, we do not have enough information to test the hypothesis that these environmental conditions lead to changes in this relationship, and only one general relationship has been applied to our whole data set. More effort is needed in order to calculate SDMD based on the particular environmental conditions of the sites.
Application of SDMD for developing thinning schedules and yield estimation As an example of the applicability of the SDMD, we propose three different approaches outlined previously (Figure 2).
From the models, the maximum stem volume and biomass at the end of the production cycle was produced in Model C (no intervention), followed by Model A (thinning of 50% of N when SDI reached its 80%) and B (managing the tree density in order to maintain COV between 40 and 50%, Table 3). In this exercise, we determined the end of the cycle when the trees reached 50 cm of Dg, but not considering the amount of years needed to reach that tree size. Model B needed 3 interventions (thinning) in order to maintain COV < 50% meanwhile Model A only two ( Figure 4, Table 3). Model A produced around 96% and 98% of the VOL and BIOM, respectively, produced by the Model C (Table 3). Model C provided the highest amount of trees with large Dg (Table 3). In the case of model B, that produced the lower amount of BIOM and VOL (69 and 70% of Model C), we also have to consider the additional herbaceous yield (agriculture or animal forage), providing an annual income to the system (i.e. Ferreira et al. 2012). On the other hand, it is expected that at higher relative densities (i.e. model C), individual tree leaf area, and thus tree growth, is much lower than that of open grown trees of the same species and age on the same site (Long et al., 2004). So, the time period that this stand needs to finish the cycle should be higher than in the other two models because the trees have lower mean annual increment compared to open stands. This was likely due to the fact that thinning is associated with increases in tree-and stand-growth efficiency (amount of wood produced by unit of leaf area) basically due to an increase of available resources to remaining trees, but also possibly due to an increase in light-use efficiency (Whitehead and Beadle, 2004). Finally, another alternative is to follow the model C but harvesting all the trees when these reach the maximum size (100% SDI). In this case, trees reach a Dg of around 28 cm, and the VOL and BIOM of the stand were around 604 m3ha-1 and 425 kg . ha -1 . This value was very close to those estimated for the stand with a higher SDI described in Milione (2019). The productivity of this alternative model is higher than the first intervention of the other proposed models (Table 3). For this reason, it is important to incorporate into any economic decision the rotation period in the unthinned alternative or model C, which is longer than in the other options.
As we mentioned before, the time dimension is not considered in these diagrams, which can be assessed using site index curves if this information is available (e.g. Luis and Fonseca, 2004). Given the same initial density, stands on more productive sites will reach earlier the onset of competition, canopy closure, full site occupancy, and self-thinning than in less productive sites (Long et al., 2004). Unfortunately, there are no information about site index estimation for E. viminalis in Argentina, and we are unable to estimate any site index because the uncertainty about the age of a high proportion of the stands used to build the SDMD. Considering the energy provided by E. viminalis wood (gross calorific value of 15 and 19 MJ kg -1 with 20-15 and 5-0% moisture, respectively; Pérez et al., 2008), the productive models evaluated estimate a range of 5 and 8 TJ ha -1 at the end of the production cycle. This is higher than the average estimated for other species of Eucalypts or those reported for poplar or willow (Pérez-Cruzado et al., 2011). These observations agreed with Pérez et al. (2008) who showed that the amount of energy produced by the genus Eucalypts is one of the greatest economic impact in the north of Spain. Schönau and Coetzee (1989) have proposed that it is usually not economically justified to grow eucalypt plantations for the production of saw timber or other highly valued products unless the mean annual increment at 10 years is more than 30 m3ha-1year-1, or more than 20 m3ha-1year-1 at 20 years. Our estimations, and those provided by official forest-related sources (BA Buenos Aires Forestal 2010; SINAVIMO 2018), indicate that the E. viminalis afforestation in Pampa plains, even with the low genetic improvement of the planted materials, can reach wood volumes higher than 300 m3ha-1 with average growth rates of 30 m3ha-1y-1, which gives an opportunity to provide both solid-and also pulpwood to the industry. However, more studies are needed focused on the response of the species, and of different provenances and genotypes within it, to the large environmental gradient present in the study region, as well on its response to silvicultural treatments such as pruning, thinning or fertilization.

CONCLUSIONS
The stand density management diagram (SDMD) developed for Eucalyptus viminalis afforestation in Pampa plains allows to estimate standing wood biomass and volume at the same time than canopy cover. It therefore provides a valuable tool to plan and evaluate alternative production systems based on the optimization of individual wood volume for solid wood industry, stand biomass for pulp or energy industries, or mixed production such as silvopastoral systems. The canopy cover component could also be used to plan open-canopy systems compatible with biodiversity conservation (e.g. Fiandino et al., 2018), lower water resources use (e.g. Licata et al., 2011) and minimization of potential negative impacts on soils, such as secondary salinization (e.g. Nosetto et al., 2008). All these environmental aspects, added to the consideration of the proportion of wood products that goes into long-lived products and to energy production, could allow planning more sustainable production systems in the framework of current global challenges, such as climate change.
Notwithstanding this, more effort is needed in order to estimate a site index equation that allows the SDMD to estimate productivity (e.g. biomass per year) and rotation period, a key variable in any economic evaluation of different alternatives.