76 Banko Janakari, Special Issue No. 4 Volume and taper equations are used for estimating timber volume and biomass of a tree. Despite their usefulness, precise and site specific equations are still lacking for commercially important tree species in Nepal. The study was carried out at Chandak Chatiya Mahila Community Forest in Bardia district and Lumbini Collaborative Forest of Saljhandi in Rupandehi district in western Terai of Nepal. A destructive sampling method was used and selected fifteen Sal trees (Shorea robusta Gaertn. f.) from Saljhandi (site 1) and eighteen trees from Bagnaha (site 2) randomly to calibrate an individual tree volume and a stem taper function. At first, a non-linear stem taper function was calibrated using stem diameters outside bark at different heights above ground as response variable and D (diameter at breast height), H (total height), h (height of interest) as predictors. Then, effect of crown characteristics on stem taper was evaluated. As stem HCB (height to crown base) was found to affect stem taper, its usefulness in existing stem volume equation was tested. Empirical relationships between V (stem volume) as a response variable and D, H, HCB and sites in Bardia and Rupandehi districts as predictors were established using a linear mixed modeling approach. Our result showed that, instead of H, use of HCB in stem volume equation increased model prediction accuracy and reduced prediction bias. Applicability of the suggested models for predicting individual S. robusta tree volume and stem taper is discussed. Key words: Crown characteristics, non-linear mixed model, Shorea robusta, taper, volume Modeling taper and volume of Sal (Shorea robusta Gaertn. f.) trees in the western Terai region of Nepal R. Silwal1*, S. K. Baral2 and B. B. K. Chhetri3 Tree volume and taper equations are vital for forest management, and they are lacking for commercial tree species (e.g. Sal-Shorea robusta Gaertn. f.) in Nepal. Currently, there is a growing interest in multiple-product timber harvesting. This requires precise stem taper and volume equations for improved prediction of volume at individual tree and stand levels. We also need to know what portion of a tree can be used for specific products, and need to identify the entire array of products that can be obtained from a specific stand. Despite their usefulness, volume and taper functions have been rarely studied in Nepal. Sharma and Pukkala (1990) have developed volume equation for twenty one species of Nepal including S. robusta. The volume equations compiled by Tamrakar (2000) and developed by DFRS (2006) are from small-sized trees. Elsewhere, a considerable amount of work on modeling tree volume and stem taper has been done (Clutter et al., 1983; Kozak et al., 1969; Max and Burkhart, 1976; Newnham, 1988). Although various methods have been proposed for developing taper equations (Bennett and Swindel, 1972; Demaerschalk, 1973; Demaerschalk and Kozak, 1977; Goulding and Murray, 1976; Kozak et.al., 1969; Kozak and Smith, 1966; Max and Burkhart, 1976), the information is either theoretical or limited primarily to softwood species (Martin, 1981). Specifically, in the Nepalese context, there are no available stem taper functions for S. robusta. Hence, this study aimed at (1) evaluating stem taper profile to identify tree characteristics that affect stem taper, and (2) use this information for improving existing stem volume equation for S. robusta. Materials and methods Study area The study was carried out at Chandak Chatiya Mahila Community Forest in Bardia district 1 National Reconstruction Authority, Singh Durbar, Kathmandu, Nepal. *E-mail: rmshsilwal@gmail.com 2 Northern Hardwoods Research Institute, University of Moncton, Campus d’ Edmundston, NB, Canada 3 Institute of Forestry, Pokhara Campus, Tribhuvan University, Pokhara, Nepal. 77 Banko Janakari, Special Issue No. 4 and Lumbini Collaborative Forest of Saljhandi in Rupandehi district in western Terai of Nepal (Fig. 1). Rupandehi is slightly warmer and wetter (higher average annual rainfall) than Bardia (Table 1). Fig. 1: Map of the study area: (a) physiographic map of Nepal indicating location of study districts (b) study sites showing GPS points of trees in administrative map of VDCs Climate, soil and vegetation Saljhandi, Rupandehi district The climate of the site is sub-tropical and sub- humid with regular monsoon between June and August. Frost occurs seldom and the annual average number of days with minus temperature is zero (Jackson, 1994). Mean total annual precipitation is 2452 mm of which more than 80% falls from June to September. Monthly mean minimum and maximum temperature are 17.8°C and 31.4°C, respectively with an absolute minimum of 4.3°C (Jackson, 1994). Saljhandi site is flat and fertile. The soil of this site is loamy, deep, well drained with adequate nutrients. According to the map of the Land Resource Mapping Project (LRMP), this area belongs to the class I, most suitable land for agriculture and forestry. Actual land use for this area is degraded tropical mixed hardwood forest. Soil physical and chemical properties are exceptionally good for forestry use (FRP, 1989). This site has Sal forest, which consists of more than 80% Sal trees. Other associated tree species are Terminalia belerica, Terminalia chebula, Terminalia alata, Anogeissus latifolia, Phyllanthus emblica, Semicarpus anacardium, Lagerstroemia parviflora, Syzygium cumini, Adina cordifolia, Mallotus philippinensis, Myrsine semiserrata and Cassia fistula (Ojha et al., 2008). Chandak Chatiya Mahila CF, Bagnaha, Bardia district This site has a sub-tropical monsoon climate with three distinct seasons in the annual cycle: hot season (March–June), Monsoon (July–October) and winter (October–February). About 90% of the precipitation occurs from July to September. The absolute maximum temperature (41°C) and minimum temperature (3.1°C) were recorded in May 1996 and January 1987, respectively. The recorded highest and lowest rainfalls were 2798 mm and 1592 mm in 1990 and 1992, respectively. Mean annual rainfall at Chisapani at the foot of the Chure hill is 2230 mm whereas it is 1560 mm at Gularia, in an agricultural area to the south of the study site. Most of the areas of Karnali and Bardia fall into Bhabar which is broad alluvial plain that slopes gently away from the base of the Churia to India. Bhabar deposits are composed of cobbles, boulders, and coarse sand layers amidst clay and silt (HMGN, 1971). The soils are well drained and relatively deep. The study area is predominantly underlain by sandy loams and followed by sands and gravels. (Dinerstein, 1979) More than 70% of the forest is covered by Sal trees. A vegetation study conducted by Dinerstein (1979) classified six major vegetation types. It was later modified by Jnawali and Wegge (1993) to seven major vegetation types. Major associated species are T. alata, Buchanania latifolia, Dalbergia sissoo Acacia catechu, S. cumini, M. phillippensis, Bombax ceiba, A. cordifolia, Casearia tomentosa, Mitragyna parviflora, Phragmatis karka and Arundo donax. Silwal et al. 78 Banko Janakari, Special Issue No. 4 Sampling and measurement Thirty-three trees were felled and measured from May to June 2013. The data were collected from trees growing in natural forests of two different locations in western Terai region of Nepal. Representative healthy and undamaged trees of different DBH classes were selected as sample trees. Fifteen trees from Saljhandi (site 1) and eighteen trees from Bagnaha (site 2) were selected randomly. DBH (1.3 m above ground level), and crown diameter of the sample trees were measured and a photograph was taken before felling the tree. Total height and height to the base of live crown (HCB) were measured separately. Stem profile data (diameter outside bark and height above ground) was obtained at eight points on the bole: stump (30 cm above ground) and approximately at 1/8, 1/4, 3/8, 4/8, 5/8, 6/8, and 7/8 of the total height (Martin, 1981). The stem section of the tree (total height minus stump height) was divided into sub-sections of maximum three meter length and measured the diameter at three points (Diameter at two ends and middle of the sub- section). The large branches were treated as poles and measured by using above process for volume estimation in the stem analysis. Average sample tree characteristics are presented in table 1. Modeling Approach Modeling stem taper At first, base model for stem taper was calibrated using Ormerod function (Ormerod, 1973), in which ‘tree’ and ‘site’ factors were considered as random effects component of a non-linear mixed effect model using nlme package (Pinheiro et al., 2017) available in R (RStudio Team, 2015). Auto-correlation due to repeated stem diameter measurement within a tree was modelled using corAR(x) function. Fig. 2: Relationship between tree level random effects and height to crown base (HCB) Observed heteroscedasticity was accounted for with a variance function (varpower function). An Akaike Information Criteria (AIC) was used to assess usefulness of the covariates in fixed as well as random effect components of the model. Site was dropped from the random effect component, as AIC did not justify its usefulness in the model. Tree level random effect (bi) showed that there was significant tree to tree variation on stem tapering. Therefore, we extracted tree level random effect for each tree, and correlated it with several tree variables e.g. height to crown base, crown diameter, crown projection area, crown volume to identify the important tree characteristics that influence stem taper. HBC was found to be positively correlated with tree level random effect (RE) parameters (Fig. 2). Addition of HCB in stem taper equation as a covariate reduced mean bias of the estimation (Table 3). ...........(T1) where, Table 1: Sample tree and site characteristics (standard deviation of the corresponding values is provided in parenthesis) Attributes Saljhandi, Rupandehi Bagnaha, Bardia Average Range Average Range Number of trees 15 NA 18 NA DBH (cm) 55.04 (17.61) 37.5-91 38.69 (27.45) 5.2-92.5 Height (m) 30.97 (3.27) 25-37 21.07 (9.44) 4.5-34.2 HCB (m) 11.1 (4.39) 3.32-17.9 4.4 (2.23) 1.05-10.4 Rainfall (mm) 2296 NA 2230 NA Temperature (˚C) 24 NA 19.5 NA NA - not applicable Silwal et al. Fig. 2: Reletionship between tree level random efffects and height to crown base (HCB) 79 Banko Janakari, Special Issue No. 4 dij = diameter of tree i at height j (cm) Di = diameter at breast height of trees i (cm) Hi = total height of tree i (m) h = height of interest from the ground (m) εij = error component of the model b0 and b1 = fixed effect parameters to be estimated. bi = tree level random effect parameter for tree i .......... (T2) where, HBCi = height to base of live crown of tree i (m) Modeling total stem volume Modeling stem volume was started with the calibration of Sharma and Pukkala (1990) equation using our data (equation V1). As tree height was not significant in the model, it was dropped from the model and then the resultant model was presented as equation (V2). Stem taper analysis showed that trees of lower HCB had more tapered stem, thus, we suspected that such trees might have smaller volume for the given diameter. Therefore, HCB was added to equation (V2) to test if the addition of HCB improves prediction ability of total stem volume or not. Eventually, equation (V3) was evaluated as final model. Since our data were collected from two different sites, we adopted mixed modeling approach to account site as random effect in the model. An AIC criterion was used to assess usefulness of the covariates in fixed as well as random effect components of the model. Effect of site as random effect was tested for intercept as well as for slope parameters. However, AIC and ANOVA test justify the usefulness of random effect only for slope parameter. Once the models were fitted, assumptions of the regression analysis were checked. The plot for the standardized residuals versus the fitted values showed that the final models did not violate any model assumptions (diagnostic plots not shown). Modeling was done in R (RStudio Team, 2015) using nlme package (Pinheiro et al., 2017). Ln(Vij) = (a + ai) + bLn(Di) + cLn(Hi) + εij ......(V1) Ln(Vij) = (a + ai) + bLn(Di) + εij ....................(V2) Ln(Vij) = (a + ai) + bLn(Di) + cLn(HCBi) + εij .....(V3) where, Vij = total stem volume (m3) of a tree i in site j a,b,c = fixed effect parameters to be estimated ai = site level random effect parameter to be estimated Results and discussion Stem taper Since we had a hierarchical dataset which was grouped as ‘multiple diameter measurements within a tree’ and ‘trees were clustered in a site’, mixed modeling approach was used. The model was parameterized several times keeping the fixed-effect and random effect specification constant. We first included correlation structure, then the variance function. According to AIC statistics, every additional covariance feature, i.e., the correlation structure, and the variance function significantly improved the likelihood ratio (Table 2). Then, we added new covariate in the final model. Stepwise model building process was shown in table 2. Table 2: AIC (Akaike Information Criteria) comparison of different model forms of stem taper equation (Note: The smaller the AIC, the better the model) Models AIC Likelihood ratio NLME 1326 NA NLME+ corAR 1323 4.97 NLME + corAR + VarPower 1103 222.4 Final model (NLME + corAR + VarPower + HCB) 1102 2.6 Sequential improvement is on distribution of residuals through different predicted values (Fig. 3). Figure 3 (A and B) shows that bias increased with increasing stem diameter. However, no such trend was found in figure 3 (C and D) besides a few extreme values and slightly increasing bias towards smaller diameter values. As smaller diameter section of a tree stem mostly belongs within the crown, this bias might have introduced by branches or stem swelling at the lower part of the branch junctions. Silwal et al. 80 Banko Janakari, Special Issue No. 4 Parameter estimates of two final models (with HCB and without HCB as covariate) are presented in table 3. As Eq. (T2) was less biased, it is advised to use this equation, when HCB information is available along with DBH and tree height (Table 3, Mean bias). We can see that how HCB influences stem taper in figure 4. All three lines are the predicted stem taper profiles of three trees that have almost the similar D (40–43 cm) but different total height and HCB. It clearly shows that trees with lower HCB generally have more tapered stem, which is similar to the results observed by Adu-Bredu et al. (2008), and MacFarlane and Weiskittel Fig. 3: Residuals of the taper model (A) NLME base model, (B) Model with autocorrelation function, (C) Model with autocorrelation and weight function, (D) Final model Table 3: Parameter estimates of stem taper equations. Pseudo R-squared=1-(residual sum of squares/total sum of squares), Mean bias=average (observed value-predicted value) Parameters Model with HCB (Eq. T2) Model without HCB (Eq. T1) Estimates SE P-value Estimates SE P-value b0 1.9682 0.01416 <0.01 1.9873 0.0084 <0.01 b2 0.0025 0.0015 0.09 b1 1.4798 0.0373 <0.01 1.481 0.0376 <0.01 RE parameters σi 2 0.00091 0.00104 σij 2 0.00037 0.00037 Pseudo R-squared 0.95 0.94 Mean bias -0.0859 0.1107 Silwal et al. 81 Banko Janakari, Special Issue No. 4 (2016). Even if HCB in our model found to be marginally insignificant (p=0.09), we still keep in the model as it reduced prediction bias. We also suggest future studies to prove its usefulness in predicting stem taper as our study was limited to a sample of 33 trees only. Fig.4: Predicted stem taper profiles for trees of given diameter, height and height to crown bases. D=diameter at breast height, H=total tree height, HCB=height to crown base Stem volume Modeling tree stem volume showed that tree height was not found to be significant at the local level (Table 4, Eq. V1). Instead of H, HCB was found significant in stem volume equation (Table 4, Eq. V3) and reduced mean bias (Table 4, Eq. V2 and Eq. V3) and increased pseudo R-squared (Table 4, Eq. V2: Pseudo R-squared=0.94 and Eq. V3: Pseudo R-squared=0.95). This must be due to the positive correlation between branchiness and main stem taper (Ver Planck and MacFarlane, 2014). Inclusion of site as random effect might have explained the variation on stem volume for different height trees at tree height for a given diameter tree is related to site factor (Feldpausch et al. 2011). We observed that lower HCB contributed on reductions in the main stem volume for a given diameter tree, but we were not able to explain the associated volume shifted away from the main stem (into branches). Therefore, future research should focus on this issue as it is necessary to assess branch wood volume, estimate total biomass and quantify carbon stock more precisely. Bias correction The models (V1), (V2) and (V3) predict stem volume on a natural logarithmic scale. We need to back transform it to the original scale. Since a linear back transformation of predicted values are associated with a log-transformation bias (Baskerville, 1972), a correction factor (CF) that accounts this log-transformation bias was presented in table 4. For the bias correction, the predicted values should be multiplied by the correction factor provided in table 4. Conclusion General volume equation in Nepal developed by Sharma and Pukkala (1990) needs calibration with local data before using it at local level. In addition, height to crown base was found to be an important variable that affects main stem volume of S. robusta for a given diameter tree at a site. Thus, we recommend to use HCB in stem volume equation for increasing prediction ability Table 4: Parameter estimates for stem volume equations (ai : Bardia =-0.0891 and Rupandehi = 0. 0891). Pseudo R-squared=1-(residual sum of squares/total sum of squares), Mean bias = average (observed value-predicted value), CF (correction factor for log transformation bias) = exp (σi 2 + σij 2 )/2 Parameters Equation (V1) Equation (V2) Equation (V3) Estimates SE P-value Estimates SE P-value Estimates SE P-value a -9.095 0.4737 0 -9.1664 0.33198 0 -9.3124 0.2727 0 b 2.5601 0.2188 0 2.51635 0.07189 0 2.4202 0.0810 0 c -0.073 0.3422 0.83 0.2694 0.1077 0.02 RE parameters σi 2 0.0847 0.07958 0.0220 σij 2 0.0775 0.07519 0.0681 Pseudo r-squared 0.94 0.94 0.95 Mean bias -0.0190 -0.0074 0.0023 CF 1.08 1.08 1.05 Silwal et al. 82 Banko Janakari, Special Issue No. 4 and reducing prediction bias. The derived stem volume equation is recommended to produce local volume table of S. robusta in Rupandehi and Bardia District. The stem taper equation is useful for forest managers to calculate stem volume up to any desirable merchantable limit (e.g. 10 cm or 20 cm top diameter). These models can be applied to similar stand condition (basal area, canopy cover, stand age) from where the study data were obtained. As this is the first model for S. robusta stem taper, it is recommended to test in other sites too. Care must be provided when using this model in other sites and predicting beyond the observed range of tree size. Acknowledgements This research was accomplished under an MSc research grant funded by the Forest Resource Assessment (FRA) Nepal. We are thankful to Mr. Sahas Man Shrestha, Director General and team of Department of Forest Research and Survey for technical inputs. We also thank Mr. Michael David Hawkes, Chief Technical Advisor and FRA Nepal team for financial and logistic support. We express our deep gratitude to District Forest Officers: Dipak Gyawali, Rupandehi, Dhananjaya Lamichhane, Bardia and their team for generous cooperation in the field. We appreciate the contribution of Mr. Shovit Koirala and Mr. Nav Raj Shahi in data collection during the field-work. We are grateful to the Lumbini Collaborative Forest User Group and Chandak Chatiya Mahila Community Forest User Group for providing support to carry out the fieldwork. We also thank an anonymous reviewer for helpful comments that improved the manuscript. References Adu-Bredu, S., Bi, A. F. T., Bouillet, J. P., Me, M. K., Kyei, S. Y. and Saint-André, L. 2008. An explicit stem profile model for forked and un- forked teak (Tectona grandis) trees in West Africa. Forest ecology and management 255 (7): 2189–2203. Baskerville, G. L. 1972. Use of logarithmic regression in the estimation of plant biomass. Canadian Journal of Forest Research 2 (1): 49–53. Bennett, F. A. and Swindel, B. F. 1972. Taper Curves for Planted Slash Pine. USDA Forest Service, Southeastern Forest Experiment Station, Asheville, NC. USA. Clutter, J. L., Fortson, J. C., Pienaar, L. V., Brister, G. H. and Bailey, R. L. 1983. Timber Management: A Quantitative Approach. John Wiley and Sons, USA. Demaerschalk, J. P. 1973. Integrated systems for the estimation of tree taper and volume. Canadian Journal of Forest Research 3 (1): 90–94. Demaerschalk, J. P. and Kozak, A. 1977. The whole-bole system: a conditioned dual- equation system for precise prediction of tree profiles. Canadian Journal of Forest Research 7 (3): 488–497. DFRS. 2006. Local Volume Tables for Major Tree Species in Dhaulagiri Area. Department of Forest Research and Survey (DFRS), District Forest Offices: Parbat, Baglung and Myagdi, Livelihoods and Forestry Programme, Dhaulagiri area, Nepal. Dinerstein, E. 1979. An ecological survey of the Royal Bardia Wildlife Reserve, Nepal. Part I: vegetation, modifying factors, and successional relationships. Biological Conservation 15 (2): 127–150. Feldpausch, T. R., Banin, L., Phillips, O. L., Baker, T. R., Lewis, S. L., Quesada, C. A., Affum-Baffoe, K., Arets, E. J. M. M., Berry, N. J., Bird, M., Brondizio, E. S., Camargo, P. de, Chave, J., Djagbletey, G., Domingues, T. F., Drescher, M., Fearnside, P. M., Franca, B. M., Fyllas, N. M., Lopez-Gonzalez, G., Hladik, A., Higuchi, N., Hunter, M. O., Iida, Y., Salim, K. A., Kassim, A. R., Keller, M., Kemp, J., King, D. A., Lovett, J. C., Marimon, B. S., Marimon-Junior, B. H., Lenza, E., Marshall, A. R., Metcalfe, D. J., Mitchard, E. T. A., Moran, E. F., Nelson, B. W., Nilus, R., Nogueira, E. M., Palace, M., Pati˜no., Peh, K. S. H., Raventos, M. T., Reitsma, J. M., Saiz, G., Schrodt, F., Sonke, B., Taedoumg, H. E., Tan, S., White, L., Woll, H. and Lloid, J. 2011. Height-diameter allometry of tropical forest trees. Biogeosciences 8: 1081–1106. FRP. 1989. Research Protocol for Natural Sal Regeneration Management, Jogikuti, Butwal. Forest Research Project, Katmandu. Unpublished. Goulding, C. J. and Murray, J. C. 1976. Polynomial taper equations that are compatible with tree Silwal et al. 83 Banko Janakari, Special Issue No. 4 volume equations. New Zealand Journal of Forest Science 5 (3): 313–322. HMGN. 1971. Soil Survey of Bardia Division, Forest Resources Survey Office, Department of Forest, His Majesty’s Government, Nepal, Kathmandu, Nepal. Jackson, J. K. 1994. Manual of Afforestation in Nepal. Second edition. Forest Research and Survey Center, Babarmahal, Kathmandu,Nepal Jnawali, S. R. and Wegge, P. 1993. Space and habitat use by a small re-introduced population of greater one-horned rhinoceros (Rhinoceros unicornis) in Royal Bardia National Park in Nepal. In Rhinoceros biology and conservation. Zoological Society of San Diego, San Diego, USA, 208–217. Kozak, A. and Smith, J. H. G. 1966. Critical analysis of multivariate techniques for estimating tree taper suggests that simpler methods are best. The Forestry Chronicle 42 (4): 458–463. Kozak, A., Munro, D. D. and Smith, J. H. G. 1969. Taper functions and their application in forest inventory. The Forestry Chronicle 45 (4): 278–283. MacFarlane, D. W. and Weiskittel, A. R. 2016. A new method for capturing stem taper variation for trees of diverse morphological types. Canadian Journal of Forest Research 46 (6): 804–815. Martin, A. J. 1981. Taper and Volume Equations for Selected Appalachian Hardwood Species. Broomall, PA: U. S. Department of Agriculture, Forest Service, Northeastern Forest Experiment Station, USA. Max, T. A. and Burkhart, H. E. 1976. Segmented polynomial regression applied to taper equations. Forest Science 22 (3): 283–289. Newnham, R. M. 1988. A Variable-form Taper Function. Petawawa National Forestry Institute, Canada. Ojha, S. K., Acharya, K. P., Acharya, B. and Regmi, R. 2008. Simple coppice management options for the Sal (Shorea robusta Gaertn. f.) forests in the Terai of Nepal. Banko Janakari 18 (1): 32–41. Ormerod, D. W. 1973. A simple bole model. The Forestry Chronicle 49 (3): 136–138. Pinheiro, J., Bates, D., DebRoy, S. and Sarkar, D., 2007. Linear and Non-linear Mixed Effects Models. R package version, 3. R Studio Team. 2015. RStudio: Integrated Development for R. RStudio, Inc., Boston, MA, viewed on 20/02/2017. http://www. rstudio.com/. Sharma E. R. and Pukala T. 1990. Volume Equation and Biomass Prediction of Forest Trees of Nepal, Publication 47. Forest Research and Statistics Division, Kathmandu, Nepal. Tamrakar, P. R. 2000. Biomass and Volume Tables with Species Description for Community Forest Management. Ministry of Forests and Soil Conservation, Natural Resource Management Sector Assistance Programme (NARMSAP), Tree Improvement and Silviculture Component, Kathmandu, Nepal. Ver Planck, N. R. and MacFarlane, D. W. 2014. Modeling vertical allocation of tree stem and branch volume for hardwoods. Forestry 87 (3): 459–469. Silwal et al.