A HYBRID MODEL TO ESTIMATE NATURAL RECRUITMENT AND GROWTH IN STANDS FOLLOWING MOUNTAIN PINE BEETLE DISTURBANCE by Derek Felix Sattler B.Sc., University of Alberta, 1998 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Forestry) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) January, 2009 © Derek Felix Sattler 2009 ii Abstract A method of linking SORTIE-ND and PrognosisBC was developed for the purpose of predicting natural regeneration and forecasting future stand conditions in mountain pine beetle (Dendroctonus ponderosae Hopkins - MPB) attacked stands in the Interior Douglas-fir (IDF) and Sub-Boreal Spruce (SBS) biogeoclimatic ecosystem zones of central and southeastern British Columbia. PrognosisBC, a spatially-implicit growth model, lacked a submodel suitable for predicting natural regeneration in unsalvaged MPB-disturbed stands. To fill this gap, estimates of regeneration (trees < 7.5 cm diameter at breast height - DBH) were supplied to PrognosisBC using the light-mediated forest dynamics model SORTIE-ND and the linked model was used to forecast future stand conditions. In order to improve results, a density-dependent system of crown allometry equations to predict crown depth and crown radius was developed and then added to SORTIE-ND. The equations used stand-level measures of stems ha-1, basal area (m2 ha-1), and the basal area of trees taller than the target tree to explicitly account of the effects of crowding on the crown axes. Additionally, crown radius and crown depth were used as dependent regressors. The equations were fit using a nonlinear three-stage least squares estimator and generally provided good estimates of crown depth and crown radius for lodgepole pine (Pinus contorta var. latifolia), hybrid spruce (Picea engelmannii x glauca (Moench) Voss), Douglas-fir (Pseudotsuga menziesii var. glauca (Beissn.) Franco) and trembling aspen (Populus tremuloides Michx.). Tests of the hybrid model with the improved system of crown allometry equations were performed using reconstructed plot data collected from natural stands disturbed by MPB 25-years ago. The hybrid model provided good estimates (small mean bias and low root mean square error) for the basal area of advance regeneration (2 < DBH < 7.5 cm) for lodgepole pine (Pinus contorta var. latifolia). The best estimates were achieved when trees <7.5 cm DBH were transferred from SORTIE-ND to PrognosisBC 15-years after MPB-disturbance. For trees < 2 m in height, poor estimates of stems ha-1 where obtained. Despite the shortcomings with respect to trees < 2 m tall, the results suggest that linking SORTIE-ND and PrognosisBC is an effective method of building a hybrid model capable of being used in MPB-disturbed forests. However, full parameterization of the SORTIE-ND model is likely needed to obtain accurate estimates for all sizes of natural regeneration. iii Table of Contents Abstract ....................................................................................................................................................................... ii Table of Contents ....................................................................................................................................................... iii List of Tables ............................................................................................................................................................... v List of Figures ............................................................................................................................................................ vi Acknowledgements ................................................................................................................................................... vii Co-authorship Statement .......................................................................................................................................... viii 1. Predicting natural regeneration following mountain pine beetle disturbance ................ 1 1.1 Introduction ........................................................................................................................................................... 1 1.2 Objectives .............................................................................................................................................................. 4 1.3 Thesis structure ...................................................................................................................................................... 4 1.4 Literature cited....................................................................................................................................................... 5 2. A system of nonlinear simultaneous equations for crown height and radius ................. 9 2.1 Introduction ........................................................................................................................................................... 9 2.2 Methods ............................................................................................................................................................... 12 2.2.1 Study area ...................................................................................................................................................... 12 2.2.2 Sampling approach and data description ....................................................................................................... 14 2.2.3 Development of Crown Models .................................................................................................................... 15 2.2.4 Parameter estimation using simultaneous regression .................................................................................... 17 2.2.5 Accuracy of Fitted Crown Models ................................................................................................................ 21 2.2.6 Model Evaluation and Validation .................................................................................................................. 21 2.3 Results ................................................................................................................................................................. 22 2.3.1 System of Crown Allometry Equations ......................................................................................................... 22 2.3.2 Model Evaluation and Validation .................................................................................................................. 31 2.4 Discussion ........................................................................................................................................................... 32 2.5 Conclusions ......................................................................................................................................................... 36 2.6 Literature cited..................................................................................................................................................... 37 3. A hybrid model to estimate natural recruitment and growth in stands following mountain pine beetle disturbance ...................................................................................... 43 3.1 Introduction ......................................................................................................................................................... 43 3.2 Methods ............................................................................................................................................................... 45 3.2.1 Study area ...................................................................................................................................................... 45 3.2.2 Data description and preparation ................................................................................................................... 47 3.2.2.1 Reconstruction of tree-lists .................................................................................................................. 48 3.2.3 The SORTIE-ND Model ............................................................................................................................... 48 3.2.3.1 SORTIE-ND growth and mortality ..................................................................................................... 49 3.2.3.2 SORTIE-ND natural regeneration ....................................................................................................... 50 3.2.4 Parameterization of SORTIE-ND Behaviours .............................................................................................. 50 3.2.4 The PrognosisBC Model ................................................................................................................................. 53 3.2.4.1 PrognosisBC Growth and Mortality ...................................................................................................... 53 3.4.2 PrognosisBC Parameters ................................................................................................................................. 54 3.2.5 Hybrid Model Structure ................................................................................................................................ 55 3.2.6 Simulations ................................................................................................................................................... 56 3.2.7 Evaluation ..................................................................................................................................................... 60 3.3 Results ................................................................................................................................................................. 62 3.3.1 Stand Simulation Mean Biases ...................................................................................................................... 62 3.2 Ecological expectations ....................................................................................................................................... 67 3.2.1 Density of seedlings and GLI ........................................................................................................................ 67 3.2.2 Growth rates of regenerated trees .................................................................................................................. 69 3.4 Discussion ........................................................................................................................................................... 71 3.5 Conclusions ......................................................................................................................................................... 75 iv 3.6 Literature cited..................................................................................................................................................... 75 4. Concluding Remarks ..................................................................................................... 81 4.1 General discussion and conclusions..................................................................................................................... 81 4.2 Contributions and criticism .................................................................................................................................. 81 4.2.1 Density-dependent crown models ................................................................................................................. 81 4.2.2 A hybrid model constructed by linking SORTIE-ND and PrognosisBC ........................................................ 82 4.3 Management implications .................................................................................................................................... 83 4.4 Future directions .................................................................................................................................................. 84 4.5 Literature cited..................................................................................................................................................... 85 v List of Tables Table 2.1 Summary statistics of trees in the model dataset by biogeoclimatic ecological zone. (Std.= standard deviation). ........................................................................................................................................................................................ 15 Table 2.2 Pearson’s correlations between the two crown axes and the predictor variables explored for the crown models. ........................................................................................................................................................................................ 16 Table 2.3 Estimated coefficients and their standard errors (in parentheses) for first-stage equations. ............................ 24 Table 2.4 Estimated coefficients and their standard errors (in parentheses) from N2SLS and N3SLS fits of the crown depth model in system of equations. ............................................................................................................................... 25 Table 2.5 Estimated coefficients and their standard errors (in parentheses) from N2SLS and N3SLS fits of the crown radius model in system of equations. .............................................................................................................................. 26 Table 2.6 Mean bias (observed - predicted), Root Mean Squared Error (RMSE) and Pseudo- R2 results using N3SLS parameter estimates. ........................................................................................................................................................ 27 Table 2.7 Mean bias (observed - predicted) and RMSE (in parentheses) in metres by density class (stem ha-1) and by basal area (m2 ha-1) class for N3SLS ............................................................................................................................... 28 Table 3.1 Mean stem density (stems ha-1) and basal area (m2 ha-1) by tree size class in 1981 (estimated through reconstruction) for the 13 stands used in the hybrid modelling simulations. .................................................................. 58 Table 3.2 Mean stem density (stems ha-1) and basal area (m2 ha-1) by tree size class in 2006 for the 13 stands used in the hybrid modelling simulations. ......................................................................................................................................... 59 Table 3.3 Mean bias and root mean square error (in parentheses) of basal area and stem density classes for four transfer simulations and independently run SORTIE-ND simulations......................................................................................... 63 Table 3.4 Mean bias and root mean square error (in parentheses) of tree height classes for four transfer simulations and independently run SORTIE-ND simulations .................................................................................................................. 66 vi List of Figures Figure 2.1 Location of stands sampled for the Williams Lake dataset (LeMay et al. 2007) ........................................... 13 Figure 2. 2 Relationship between observed and predicted aspen crown depth (a) and crown radius (b), and Douglas-fir crown depth (c) and crown radius (d). ............................................................................................................................ 29 Figure 2.3 Relationship between observed and predicted lodgepole pine crown depth (a) and crown radius (b), and hybrid spruce crown depth (c) and crown radius (d). ...................................................................................................... 30 Figure 3.1 Location of stands sampled for the Williams Lake dataset (LeMay et al. 2007) ........................................... 46 Figure 3.2 Schematic of the hybrid model approach. Overstory (O/S) and understory (U/S) tree-lists of the same population are imported into SORTIE-ND and PrognosisBC at Time 1, which is shortly after disturbance. At Time 2, PrognosisBC simulations are paused while SORTIE-ND is stopped. Understory trees are exported from SORTIE-ND, formatted and imported into PrognosisBC, replacing trees within the same size class. PrognosisBC is restarted with an updated tree-list and projected to points in time within the mid-term planning horizon (Time 3) .................................. 55 Figure 3.3 Average periodic growth (cm / year) from the 13 simulated stands for pine (a and b, respectively) and spruce (c and d, respectively) seedlings (trees < 2 m tall) versus predicted global light index when simulations were run to completion in SORTIE-ND. Average periodic growth rates from each stand are plotted in 5 year intervals over a 25 year simulation. ............................................................................................................................................................... 68 Figure 3.4 Predicted periodic growth (cm/yr) for pine between 5.5 and 7.5cm DBH in SORTIE-ND and PrognosisBC using the 15 year tree-list transfer. The gap in the periodic growth trends marks when understory trees were transferred from SORTIE-ND to PrognosisBC. Each trend line is the predicted growth from a single stand. Symbols on the trend lines identify individual stands, which are run in both SORTIE-ND and PrognosisBC ................................................... 70 vii Acknowledgements The completion of this thesis would not have been possible without the help I received along the way. In particular, I would like to gratefully acknowledge my supervisor, Dr. Valerie LeMay, who set aside countless hours of her time to help me through the many challenges that this thesis presented. Her insight, encouragement and patience were invaluable. I am certain that the lessons I have learned from working with her will keep me in good standing in my future endeavors. Thank you. I am also greatly in dept to my committee members, Dr. David Coates, Dr. Bruce Larson and Dr. Peter Marshall. Amidst busy schedules, they found the time to help when I needed it most. I am very thankful for the funding that was provided by the British Columbia Forest Science Program. Without this program and the willingness to invest in my research, this thesis could not have been completed. I am extremely grateful to Shelagh Hayes. Her support and encouragement during my time at UBC was unwavering. Finally, I am grateful to the wonderful friends I have made in the Faculty of Forestry. They kept my spirits up and gave me a helping hand when it was needed. viii Co-authorship Statement I, Derek Sattler, was the principal investigator and contributor for all components of this thesis. This includes the identification and design of the research program, research, data analysis, and write-up of all manuscripts in this thesis. Additional contribution to the identification and design of the research program and help in the preparation of the manuscripts was received from Dr. Valerie LeMay. Dr. LeMay’s contribution to the identification and design or the research program included the initiation of a pilot project that helped identify SORTIE-ND as a growth model that could be linked to PrognosisBC for the purpose of estimating natural recruitment. In the preparation of this thesis, Dr. LeMay served as the principle editor, which included re-writing minor sections in the two research-based manuscripts. 1 1. Predicting natural regeneration following mountain pine beetle disturbance 1.1 Introduction Lodgepole pine trees (Pinus contorta var. latifolia) in British Columbia’s Interior Forest Region are experiencing a Mountain Pine Beetle (Dendroctonus ponderosae Hopkins, MPB) outbreak that has reached epidemic levels (Eng et al. 2006). The amount of timber affected by MPB is beyond the industrial capacity to extract and process. Furthermore, forest managers are recognizing the importance of retaining some affected areas on the landscape to support the diversity and health of future forests. Thus, a large portion of MPB-disturbed stands will go unsalvaged (Hawkes et al. 2006). It is anticipated that much of the Province’s mid- to long-term timber supply will originate from these stands (Patriquin et al. 2008). This has placed considerable importance on the need to develop quantitative growth and yield models that capture the natural dynamics of MPB-disturbed stands. Critical to these models is the estimation of natural regeneration following MPB-disturbance. Growth and yield models are an important component in the development and implementation of sustainable forest management practices (Boisvenue et al. 2004). They allow for new silvicultural systems to be tested prior to implementation and can be used to provide supplementary data in addition to costly experimental trials. Models are particularly useful for testing long-term management strategies (Coates et al. 2003), such as the management of current stands for mid- and long-term timber supply (Hyytiainen et al. 2006). Although models use a variety of different approaches to predict forest growth, most are constructed by linking together submodels that correspond to a specific biological process (Robinson and Ek 2003). To be useful, a model must include the full range of identifiable process within the system being modeled (Ek et al. 1997). In unmanaged forest systems, a process that poses a major challenge to modelling is natural regeneration (Miina and Saksa 2006; Sagnard et al. 2007). Largely, this is due to the high variability of regeneration, which is influenced by a long list of biotic and abiotic factors (Maguire and Forman 1983). These include differences in reproductive strategies, shade tolerance, interspecific competition within and across vertical stratum and substrate preference (LePage et al. 2000). Regeneration submodels may use empirical or biologically-based methods to predict 2 recruitment (Vanclay 1994). Empirically based regeneration submodels generally require the use of long-term plot data containing extensive regeneration measurements. For example, empirical equations that use regression methods to predict regeneration from plot variables need historical plot data to calibrate equation parameters. For empirically-based regeneration imputation methods, such as most similar neighbour (MSN) analysis, historical plot data are matched, through shared plot variables, to units with missing regeneration measurements, providing them with an estimate of recruitment (Hassani et al. 2004).However, regeneration submodels of this type run into problems when: 1) long-term plot data for the target population are scarce; and 2) current biological conditions within the target population no longer match the conditions of long- term plot data. The growth model commonly applied to forests in southeastern British Columbia (BC), PrognosisBC, uses an MSN approach to predict natural regeneration (Hassani et al. 2004; LeMay et al. 2006). Within unmanaged stands, the high variability of natural regeneration recruitment and growth has been difficult to characterize using MSN models (Martin et al. 2002). For MSN methods to improve, a substantial amount of plot data are required since forests in this area are composed of heterogeneous stands with variable biotic and abiotic conditions. However, due to the cost of collecting data, many forest conditions are under-sampled. Thus, the ability of the submodel to accurately predict regeneration is limited (Moeur and Stage 1995, Hassani et al. 2004). There remains considerable interest in finding an alternative system of estimating natural regeneration that is compatible with the PrognosisBC model (Zumrawi et al. 2005). Forest disturbance regimes within southeastern BC consist of frequent small-scale gap disturbance events such as windthrow and less frequent large-scale disturbances such as fire and forest insect outbreaks (Temesgen et al. 2003). A model suitable to these complex forest systems would be an asset to silvicultural managers. To this end, the large and small tree growth submodels in PrognosisBC have been designed for use in managed or unmanaged stands of varying species composition and age classes (Boisvenue et al. 2004; Zumrawi et al. 2003). The BC Ministry of Forests and Range (2002) has been calibrating these submodels for several Biogeoclimatic Ecosystem Classification (BEC) subzone variants of the southeastern portion of the Province (Zumrawi et al. 2003). To strengthen PrognosisBC, a natural regeneration submodel capable of estimating recruitment in complex stands of southeastern BC is needed. Models that simulate a wide range of biological processes carry many advantages in terms of 3 estimating natural regeneration. The processes influencing natural regeneration are many. Thus, modelling as many of these processes as possible, should in theory, improve estimates (Vanclay 1994). This also provides greater control over the simulated environment, which encourages the model to be used for hypothesis generation, extension into new populations and testing of new management regimes (Robinson and Monserud 2003). Conversely, many empirical models only use predictor variables whose influence on regeneration is well understood. This approach offers little additional insight into other biological factors. Lastly, regeneration submodels that incorporate a diverse range of forest processes will be more effective in long-term projections, as the growing conditions are dynamically updated over the course of the projection period. It has been suggested that regeneration submodels of this type can be used where empirical methods have fallen short (Mäkelä 2003). A pilot study by LeMay et al. (2006) examined methods of estimating natural regeneration within forests disturbed by MPB. Promising results were obtained from SORTIE-ND, a spatially explicit forest dynamics model that incorporates biological theory into mechanistic functions that predict light, growth, mortality and recruitment (Astrup 2006). Forests of the study area were characterized as being heterogeneous in species composition and structure. A major driver of stand development is variable levels of disturbance caused by MPB. Stands disturbed by MPB experience a form of natural thinning. The increased light availability in the understory can stimulate natural regeneration or the release of advanced regeneration. Coates and Hall (2005) noted that survival of understory trees is largely dependent on a combination of shade tolerance and proportion of overstory mortality. Hawkes et al. (2004) found that shade tolerant species in the understory, such as spruce, showed increased growth in the years following MPB attack. Shade intolerant species in the understory, such as lodgepole pine, also experienced increased growth when they were in the advance regeneration stratum. Based on the results in LeMay et al. (2006) and the known relationships in forest structure that exists in these stands, crown size of overstory trees was identified as an important factor influencing understory light levels, which in turn controlled natural regeneration. Furthermore, it was hypothesized that the radius of overstory tree crowns would be largely influenced by limited growing space and direct physical contact by neighbouring trees. Crown architecture theory suggests that trees in open stands will have longer crowns along the vertical axis and larger crown radii (Oliver and Larson 1996; Temesgen et al. 2005). This trend is reversed in dense stands, where the effects of competition, either through direct physical contact or indirect means (e.g., competition for light), cause crown depth to shorten and crown radius to decrease. To improve 4 estimates of natural regeneration from SORTIE-ND, these biological principles needed to be incorporated into the model. From a biological point of view, the regeneration submodel in SORTIE-ND seemed well suited to these forests conditions. Moreover, LeMay et al. (2006) noted that with additional field data, many processes in SORTIE-ND could be parameterized for the study area. This raised the question: could the SORTIE-ND regeneration submodel be used to provide PrognosisBC with estimates of natural regeneration for stands disturbed by MPB? Furthermore, would this linked system provide defensible estimates of growth over a mid-term planning horizon of 25 years? Since most forest models are a conglomerate of submodels, it is possible to substitute one submodel for another. Robinson and Ek (2003) used this approach to develop Forest 5, a hybrid model of forest growth and stand dynamics for the Great Lakes region. Using existing models, they were able to extract individual submodels capable of simulating a discrete biological process, and link them together to create a functional growth model capable of simulating single and multiple cohort stands. There were several advantages to using this approach. First, using existing submodels allowed for quick assembly. Second, the approach afforded flexibility in terms of selecting submodels suitable for the target population and the type of data available. Lastly, the interchangeability allowed for easy comparisons between competing modules. Thus, empirical submodels could be contrasted against biologically- or processed-based submodels. Because of the modularization of forest processes in PrognosisBC and SORTIE-ND, it seemed reasonable to apply this approach. 1.2 Objectives There are two main objectives for this thesis. The first objective was to develop models for predicting overstory tree crown size in unmanaged stands affected by MPB based on quantifiable measures of tree size and stand density. This was identified as a possible issue in using SORTIE- ND to forecast regeneration in a preliminary study. The second objective was to develop and test a hybrid model capable of predicting natural regeneration in MPB-disturbed stands in the Cariboo Forest Region of British Columbia. This second objective focused on the integration of two models, the empirically based model PrognosisBC, and the mechanistic model SORTIE-ND. 1.3 Thesis structure This thesis is structured in a manuscript-based format. This chapter provides the background 5 information necessary to understand why the main questions of the thesis were asked. Chapter 2 deals with the development of models used to predict crown size in stands disturbed by MPB. The focus of Chapter 3 is on describing the hybrid model that was developed by linking SORTIE- ND to PrognosisBC after the new crown models were added to the SORTIE-ND model. Chapters 2 and 3 were written with the intention of being submitted for publication; hence, there are some minor redundancies in the methods since analyses were performed on the same sample data. Furthermore, there is extra attention to detail when describing specific methods. This is particularly true in Chapter 2, where the methods used to fit the crown models are thoroughly described. In preparation for submission to a peer-reviewed journal, the level of detail used to describe specific methodologies will be reduced. The final chapter provides conclusions and discusses the relevancy of the findings to the field of biometrics and growth and yield modelling. Based on the outcomes, future directions of study are proposed. 1.4 Literature cited Astrup, R. 2006. Modeling growth of understory aspen and spruce in western boreal Canada. Ph.D. Dissertation. University of British Columbia, Vancouver, BC. 156 p. Boisvenue, C., Temesgen, H., and Marshall, P.M. 2004. Selecting a small tree height growth model for mixed-species stands in the southern interior of British Columbia, Canada. Forest Ecology and Management 202: 301-312. British Columbia Ministry of Forests and Range. 2002. Ground sampling procedures: IDFdk4 and IDFxm sampling to calibrate PrognosisBC. Research branch, Cariboo forest region. 24pp. Coates, K.D. and Hall, E.C. 2005. Implications of alternative silviculture strategies in Mountain Pine Beetle damaged stands. Technical report for Forest Science Program Project Y051161. Bulkley Valley Centre for Natural Resource Research and Management, Smithers, BC. 31p. Coates, K.D., Canham, C.D., Beaudet, M., Sachs, D.L., and Messier, C. 2003. Use of a spatially explicit individual-tree model (SORTIE/BC) to explore implications of patchiness in structurally complex forests. For. Ecol. and Manage. 186: 297-310. 6 Ek, A.R., Robinson, A.P., Radtke, P.J., and Walters, D.K. 1997. Development and testing of regeneration imputation models for forests in Minnesota. For. Ecol. and Manage. 94: 129-140. Eng, M., Fall, A., Hughes, J., Shore, T., Riel, B., Walton. A., and Hall, P. 2006. Update of the infestation projection based on the 2005 provincial aerial overview of forest health and revisions to “the model” (BCMPB.v3). British Columbia Ministry of Forests and Range, Victoria, British Columbia. Hassani, B., LeMay, V.M., Marshall, P.L., Temesgen, H., and Zumrawi, A.A. 2004. Regeneration imputation models for complex stands of Southeastern British Columbia. For. Chron. 80:271-278. Hawkes, B., Taylor, S., Stockdale, C., Shore, T.L., Beukema, S., and Robinson, D. 2006. Predicting mountain pine beetle impacts on lodgepole pine stands and wood debris characteristics in a mixed severity fire regime using PrognosisBC and the fire and fuels extension. pp. 123-135. In: Lagene, L., J.Zelnik, S. Cadwallader and B. Hughes, editors. Mixed Severity Fire Regimes: Ecology and Management. November 17-19, 2004, Spokane, Washington. Washington State University Coop Extension Service/The Association of Fire Ecology, Pullman, Washington, Vol. AFE MISC03. Hawkes, B., Taylor, S., Stockdale, C., Shore, T.L., Alfaro, R.I., Campbell, R., and Vera, P. 2004. Impacts of mountain pine beetle on stand dynamics in British Columbia. Pages 175-195 in T.L. Shore, J.E. Brooks, and J.E. Stone, editors. Mountain Pine Beetle Symposium: Challenges and Solutions. Kelowna, British Columbia October 30-31, 2003, Natural Resources Canada, Canadian Forest Service, Pacific Forestry Centre, Information Report BC-X-399, Victoria, BC. 298p. Hyytiainen, K., Ilomaki, S., Mäkelä, A., and Kinnunen, K. 2006.Economic analysis of stand establishment for Scots pine. Canadian Journal of Forest Research 36: 1179-1189. LeMay, V.M., Lee, T., Scott, R.E. Sattler, D.F, Robinson, D., Zumrawi, A-A. and Marshall, P. 2006. Modeling Natural Regeneration Following Mountain Pine Beetle Attacks in the Southern and Central Interior of British Columbia: Results for Year 1. Internal report for Nat. Res. Can., MPB Standard Contribution Agreement, PO # 8.35. 70 pp. http://www.forestry.ubc.ca/prognosis/extension.html accessed November, 2008. 7 LePage, P.T., Canham, C.D., Coates, D., and Bartemucci, P. 2000. Seed abundance versus substrate limitation of seedling recruitment in northern temperate forests of British Columbia. Can. J. For. Res. 30: 415-427. Maguire, D.A., and Forman, R.T.T. 1983. Herb cover effects on tree seedling patterns in a mature hemlock-hardwood forest. Ecol. 64: 1367-1380. Mäkelä, A. 2003. Process-based modelling of tree and stand growth: towards a hierarchical treatment of multiscale processes. Can. J. of For. Res. 33: 398-409. Martin, R.A., Hassani, B., LeMay, V., Marshall, P., Temesgen, H., and Lencar, C. 2002. Development of regeneration imputation models for the IDF dk1, dk2, and dk3 subzone variants of the Kamloops and Cariboo forest regions. Report prepared for: Forest Renewal BC, PAR02002, Loc. 724187. http://www.forestry.ubc.ca/prognosis/documents/impute_kamcariboo.pdf accessed, November, 2008. Miina, J. and Saksa, T. 2006. Predicting regeneration establishment in Norway spruce plantations using a multivariate multilevel model. New Forests 32: 265-283. Moeur, M. and A.R. Stage. 1995. Most similar neighbour: an improved sampling inference procedure for natural resource planning. Forest Science 41, 337-359. Oliver, C.D. and Larson, B.C. 1996. Forest Stand Dynamics. John Wiley and Sons, Inc. ISBN 0- 471-13833-9, New York, 520 p. Patriquin, M.N., Lantz, V.A., Stedman, R.C., and White, W.A. 2008. Working together: a reciprocal wood flow arrangement to mitigate the economic impacts of natural disturbance. Forestry 81: 227-242. Robinson, A.P. and Ek, A.R. 2003. Description and validation of a hybrid model of forest growth and stand dynamics for the Great Lakes region. Ecol. Modelling 170: 73-104. Robinson, A.P. and Monserud, R.A. 2003. Criteria for comparing the adaptability of forest growth models. For. Ecol. and Manage. 172: 53-67 8 Sagnard, F., Christian, P., Dreyfus, P., Jordano, P., and Fady, B. 2007. Modelling seed dispersal to predict seedling recruitment: Recolonization dynamics in a plantation forest. Ecol. Modelling 203: 464-474. Temesgen, H., LeMay, V. and Mitchell, S.J. 2005. Tree crown ratio models for multi-species and multi-layered stands of southeastern British Columbia. For. Chron. 81: 133-141. Temesgen, H., LeMay, V.M., Froese, K.L., and Marshall, P.L. 2003. Imputing tree-lists from aerial attributes for complex stands of south-eastern British Columbia. For. Ecol. and Manage. 177: 277-285. Vanclay, K. J. 1994. Modelling forest growth and yield: Applications to mixed tropical forests. CABI, New York. pp. 103-133. Zumrawi, A.A., LeMay, V., Temesgen, H., and Snowdon, B. 2003. Expanding PrognosisBC to the Chilcotin variant of the interior Douglas-fir dry cool subzone: Large tree radial growth models. Report prepared for: Forestry Innovation Investment, Project No: FIIRFP RES2002/03. Zumrawi, A.A., LeMay, V., Marshall, P.M., Hassani, B.T. and Robinson, D. 2005. Implementing a PrognosisBC regeneration sub-model for complex stands of Southeastern and Central British Columbia. Collaborative Project between Research Branch, Ministry of Forests. Forest Resources Management Department, University of British Columbia, and ESSA Technologies, Inc. Vancouver, British Columbia, Canada. Report submitted to: Forest Science Program, Project No: Y051355, 157 p. http://www.for.gov.bc.ca/hre/pubs/azumrawi.htm accessed, June 2007. 9 2. A system of nonlinear simultaneous equations for crown height and radius1 2.1 Introduction An extensive list of inferences can be drawn from the size and shape of a tree’s crown. The list of inferences includes tree vigor (Assmann, 1970; Valentine et al. 1994), basal area increment (Hasenauer and Monserud 1996; Monserud and Sterba 1996), wood quality (Kellomaki et al. 1999; Valentine and Mäkelä 2005), and incident radiation (Canham et al. 1999). Because of their utility, crown models are often used in forest growth models. Several approaches have been developed to obtain precise estimates of crown size. Examples range between highly empirical- (Hasenauer and Monserud 1996; Soares and Tome 2001; Bechtold 2003) to biologically-based models (Biging and Gill 1997; Mäkelä 1997). Predictor variables usually include easily recorded measurements of individual tree size, which reflect the structural relationship between the size of the crown and the size of the tree bole (Baldwin and Peterson 1997; Gilmore 2001). Occasionally, variables that reflect the effect of growing space on crown size are also included (Rouvinen and Kuuluvainen 1997; Rudnicki et al. 2004; Meng et al. 2007). A simple estimate of crown size can be achieved by predicting two axes: the length of the crown along the vertical axis, known as crown depth, and the radius of the crown, measured at the base of the crown or at the point of maximum crown radius (Pacala et al. 1993, 1996). This simple method of describing crown size appears to be appropriate for a wide range of species and growing conditions (e.g., even - and uneven-aged stand conditions). A benefit of this method is that reliable estimates can often be achieved using simple models, requiring predictor variables that are inexpensive to collect. Furthermore, crown depth may be modeled as the ratio of live crown length to total tree height, which is a dimensionless measure (Hasenauer and Monserud 1996). Thus, the estimate of the vertical axis can be transformed between a dimensional and non- dimensional scale in order to suit spatial or non-spatial growth models. Of interest in this study is the role of crowns in the spatially explicit forest growth and dynamics model, SORTIE-ND (Coates et al. 2003; Astrup 2006). SORTIE-ND stems from the original SORTIE model, which was developed for use in deciduous stands of northeastern United 1 A version of this chapter will be submitted for publication. Sattler, D.F. and LeMay, V.M. Estimation of crown radius and crown depth following mountain pine beetle infestation. 10 States (Pacala et al. 1993, 1996; Canham et al. 1994). Modifications to the SORTIE-ND model have made it practical to use as a silvicultural decision making tool. One of the key components of SORTIE-ND is a light submodel which predicts incident radiation at any given point within a stand based on the sky brightness distribution for a given latitude. Within a simulated treed plot, light attenuation occurs as rays pass through individual tree crowns. The amount of shade cast by each tree is a function of its crown size and species specific light transmission coefficients (Pacala et al. 1993, 1996). Tree growth and survival is largely a function of predicted light levels, measured as the percent of full sun received at a point, reported as a global light index (GLI). Thus, estimates of crown size are important to several processes in the SORTIE-ND model. Within the model, the assumed crown shape is a vertical cylinder. Cylinder dimensions are based on estimates of crown depth and crown radius obtained using simple allometric equations. The choice of equations to use includes exponential or Chapman-Richards functions where estimated crown depth and crown radius are expressed as: 1 1ˆ bDBHaRadiusnCrow ×= (2.1) 2 2ˆ bHeightaDepthnCrow ×= (2.2) [ ] 131ˆ 31 CDBHbeaiRadiusnCrow ×−−+= (2.3) [ ] 241ˆ 42 CHeightbeaiDepthnCrow ×−−+= (2.4) where ,,,,,,,,,,, 12143214321 Ciibbbbaaaa and 2C are parameters to be estimated, DBH is the diameter outside bark at 1.3 m above ground in cm and height is the total tree height in m. The simplicity of the equations avoids the pitfalls of model over-specification. However, the models do not explicitly account for the effect of crowding on crown size (Hynynen 1995; Hasenauer and Monserud 1996; Pinno et al. 2001). Therefore, there is a tendency for the equations to overestimate crown axes in dense stands and underestimate axes in open stands (Astrup 2006). This, in turn, will have an affect on estimates of tree growth and regeneration. To be better represent stands with varying levels of crowding, a crown model that explicitly accounts for the effects of density is needed. In addition to the effects of density, there is a strong, biologically-based correlation between crown depth and crown radius. Each axis is influenced by a similar combination of internal 11 physiological and external processes (Oliver and Larson 1996). The hierarchy of branch control that extends from the main stem to first-order, second-order and all subsequent branch orders is the physiological process that governs crowns form. The ability of trees to maintain their characteristic crown shape or allow for shifts in the depth or radius of the crown in response to varying light conditions is largely a function of the level of control exerted by the terminal bud over the lateral branches (Oliver and Larson 1996). The effects of side shade are known to influence crown depth depending on shade tolerance and latitude (Canham et al. 1994; Canham et al. 1999). The same factors also influence radial crown growth; lateral branches are required to photosynthesize at a rate capable of sustaining their own growth (Sprugel and Hinckley 1988). As a result of these shared physiological and mechanical influences, an inherently strong relationship exists between the depth and radius of the crown. Because of the strong reciprocal relationship that exists between crown radius and crown depth, good estimates of these crown axes can be obtained by using the opposite crown axis as a predictor variable. An equation that predicts one crown axis by using the other axis as a predictor variable is an example of a simultaneous equation; the dependent variable clearly has an effect on the predictor variable, and vise versa. If both crown axes are fit in this manner, the two equations form a system of simultaneous equations. Fitting simultaneous equations or systems of simultaneous equations using ordinary least squares (OLS) results in biased estimates of the coefficients. The source of the bias is a correlation between the error term and the predictor variable (violation of one of the assumptions of OLS) (Kutner et al. 2005; Studenmund 2006). Bias of this type is termed simultaneity bias (Zellner and Theil 1962; Gallant 1975; Judge et al. 1985). To remove simultaneity bias, two- and three-stage least squares (2SLS and 3SLS) regression methods are often used. For nonlinear models, these are referred to as N2SLS and N3SLS. In forestry, simultaneous regression methods are often used to take advantage of the biologically- based reciprocal relationships that exit between growth-related variables (e.g., Borders and Bailey 1986; Hasenauer et al. 1997; Huang and Titus 1999). Because of the reciprocal relationship that exists between crown depth and crown radius, two- or three-stage regression methods are appropriate when estimates of both crown axes are needed, as is the case in the SORTIE-ND model. The main objective of this study was to develop equations to estimate crown depth and crown radius. Specific criteria were used to guide model development. First, based on the known effects 12 of crowding on crown size, the models should be able to account for changes in density. To satisfy this criterion, relationships between spatially implicit measures of density and crown size were explored. The density related variables were restricted to easily obtained measures, including the basal area of trees taller than the target tree (BALHT) (m2 ha-1), stem density (stems ha-1), Curtis` relative density (RD), basal area (m2 ha-1), and quadratic mean diameter (i.e., DBH of tree with average basal area; QDBH (cm)). Second, the equations should take advantage of the reciprocal relationship between the crown axes; hence simultaneous regression methods were used. Specifically, nonlinear models were fit using N2SLS and N3SLS. Each crown axis was included in the equation used to predict the opposite axis. Lastly, the equation forms should enable the models to provide biologically tractable estimates of crown depth and crown radius. Accordingly, various equation forms were explored. The chosen crown equations will serve as alternatives to the current allometric equations in the SORTIE-ND growth model. The focus was on obtaining estimates for lodgepole pine (Pinus contorta Dougl. ex Loud. var. latifolia Engelm), hybrid spruce (Picea engelmannii x glauca (Moench) Voss), Douglas-fir (Pseudotsuga menziesii var. glauca (Beissn.) Franco) and trembling aspen (Populus tremuloides Michx.). Data for this study were collected from stands located in the central interior of British Columbia (BC) where the most recent infestation of mountain pine beetle (Dendroctonus ponderosae Hopkins, MPB) was 25-years ago. The data were used in a previous study described by LeMay et al. (2007). The impact of historical MPB disturbance in these stands has led to the creation of highly variable stand densities where ongoing changes to the amount of light and available growing space has resulted in changes in crown shape. A secondary objective of this study was to gain a greater understanding of how tree size and stand-level density variables influence crown depth and radius in MPB-disturbed stands. 2.2 Methods 2.2.1 Study area Data for this study included stands located within 250 km of Williams Lake, BC (52°08’18.59”N 122°08’31.07”W). Many of the sampled stands were located on the Chilcotin Plateau, situated within the former Cariboo Forest Region (Figure 2.1). 13 Figure 2.1 Location of stands sampled for the Williams Lake dataset (LeMay et al. 2007) This area has a long history of disturbance events caused by MPB (Stockdale et al. 2004; Aukema et al. 2006), although evidence suggests that stand replacing and ground fires have also played a major role in influencing stand structure (Hawkes et al. 2004). Historical records gathered from aerial surveys (Wood and Unger 1996), and tree ring analyses (Hawkes et al. 2004) indicate that the most recent significant disturbance caused by MPB, aside from the current outbreak, occurred during a period beginning in the late 1970s stretching until 1984/85. The combination of disturbance events has resulted in the creation of uneven-aged, mixed-species stands throughout much of the Chilcotin plateau (Heath and Alfaro 1990). The topography of the study area is characterized by gentle slopes and an elevation between 900 and 1500 m. Lodgepole pine is either present or the dominant species in most of the upland stands in this area. Overall, species composition in stands from the Williams Lake dataset included lodgepole pine (about 80%), Douglas-fir (7%), and hybrid spruce (5%). Trembling aspen and subalpine-fir were present in some stands, but in small numbers overall. The sampled stands were located in three major biogeoclimatic ecosystem classification (BEC) zones (Meidinger and Pojar 1991): Sub-Boreal Spruce (SBS), Sub-Boreal Pine Spruce (SBPS), and 14 Interior Douglas-fir (IDF). Due to climatic similarities and for data simplification purposes, stands located within the SBPS zone were pooled with stands located in the SBS zone and are categorized as SBS stands. 2.2.2 Sampling approach and data description Sampling of the Williams Lake dataset was carried out in the summer of 2006. A brief description of the sampling approach is presented here, with greater details found in LeMay et al. (2007). Stands that had been attacked by the mountain pine beetle roughly 25 years ago were targeted for sampling. Circular plots of 11.28 m radius were systematically located in each of the 53 selected stands, beginning with a random starting point. For each plot and for live trees greater than 7.5 cm DBH, the DBH (cm), total tree height (m), and height to live crown (m) were measured. Maximum crown diameter at crown base (m) and crown diameter 90 degrees to this were measured and averaged for two randomly selected live and healthy trees of every species tallied in the plot. The ratio of live crown to total tree height was measured for these selected trees also, and later converted to crown depth (m). Based on the sizes of adult trees, structural composition of the stands, and information provided through inventory maps, it is estimated that the average age of adult trees in the sampled stands ranged between 70 and 90 years. Several of the stands that were sampled for this study showed evidence of the recent MPB outbreak. Many of the lodgepole pine trees were in stages of green (i.e., attacked within the last year) or red attack (i.e., attacked two to three years ago). Measurements of crown size were not collected for these trees. However, they were included in the estimates of live overstory tree density. Of the 53 stands sampled, 16 stands were located in the IDF zone and 37 stands located in the SBS zone. For the 16 stands in the IDF, the mean stem density of for trees > 7.5cm DBH was 677 (166; standard deviation given in parentheses) stems ha-1 with a mean basal area of 17.44 (5.88) m2 ha-1. For the stands in the SBPS/SBS zone, the mean stems ha-1 was 803 (322) and the mean basal area was 19.18 (8.13) m2 ha-1. Stands were generally composed of mixed species; however, almost all included lodgepole pine. For each stand, additional stand density measures were calculated, including Curtis’ RD, QDBH (cm), BALHT (m2 ha-1), stems ha-1, basal area (m2 ha-1). The model dataset included 650 trees with crown measures, 140 from the IDF zone and 302 from the SBS zone (Table 2.1). 15 Table 2.1 Summary statistics of trees in the model dataset by biogeoclimatic ecological zone. (Std.= standard deviation). Interior Douglas-fir Zone DBH (cm) Height (m) Crown radius (m) Crown depth (m) Species Number of trees Mean Std Mean Std Mean Std Mean Std Aspen 18 14.51 4.83 13.82 4.13 1.8 0.45 5.42 2.17 Douglas-fir 39 18.95 12.7 13.99 6.7 1.83 0.58 8.7 3.95 Lodgepole Pine 102 13.74 5.22 12.69 3.71 1.19 0.46 5.03 2.13 Hybrid spruce 8 12.56 3.87 9.32 20.6 1.60 0.43 8.39 1.85 Sub-Boreal Spruce Zone DBH (cm) Height (m) Crown radius (m) Crown depth (m) Species Number of trees Mean Std Mean Std Mean Std Mean Std Aspen 22 16.51 6.05 16.97 5.31 1.87 0.46 5.81 2.46 Douglas-fir 17 18.42 9.18 19.31 7.68 1.69 0.52 8.31 5.10 Lodgepole Pine 350 14.55 5.28 12.78 4.41 1.12 0.39 5.44 2.33 Hybrid spruce 70 16.21 7.54 13.85 6.02 1.50 0.51 10.01 4.19 2.2.3 Development of Crown Models Three criteria were used to guide development of a set of crown models that could serve as alternatives to the standard crown equations used within the SORTIE-ND model. First, the models should yield logical predictions for crown depth and crown width. For instance, estimates of crown depth should never be greater than total tree height. This criterion is often achieved through thoughtful selection of the functional form the dependent and independent variables or the form the equation itself (e.g., linear, nonlinear, exponential, logistic). The second criterion was that the equations should make use of the strong relationship that exists between crown depth and crown radius to improve predictions of the respective crown axes. Simultaneous regression techniques that avoided uncoupling this relationship were preferred. Lastly, the model should minimize the prediction bias across a range of stand densities. Given these criteria, the following nonlinear model form was selected to estimate crown depth: ( ))(1ˆ Xfe HeightDepthnCrow + = (2.5) where DepthnCrow ˆ is the estimated crown depth as the response variable, Height is defined as 16 for Eq. 2.1 to 2.4, and )(Xf is a linear function of other tree and stand-level variables. This form follows a logistic with the height of the tree whose crown is being estimated operating as the upper asymptote. This ensures that under all circumstances, the estimated crown depth will be less than the height of the tree. A power regression model of the following form was selected for crown radius: )(ˆ XfRadiusnCrow = (2.6) where RadiusnCrow ˆ is the estimated crown radius, as defined in Eq. 2.1 to 2.3, and )(Xf is a nonlinear function of tree and stand-level variables. Several tools were used to guide the selection of predictor variables. First, relationships among all pairs of variables were examined using simple scatter plots and Pearson’s simple correlations (Table 2.2). In keeping with second criterion, crown depth and crown radius were included in the list of predictor variables for the opposite axis. Table 2.2 Pearson’s correlations between the two crown axes and the predictor variables explored for the crown models. Crown Depth Species DBH Height Crown Radius Basal area Stems ha-1 BALHT QMD RD Aspen 0.249 0.446 0.404 -0.200 -0.177 -0.465 -0.142 -0.002 Douglas-fir 0.801 0.721 0.677 0.03 -0.331 -0.608 0.251 -0.066 Lodgepole Pine 0.364 0.224 0.399 0.006 -0.137 -0.105 0.176 -0.029 Hybrid spruce 0.792 0.665 0.595 0.153 -0.059 -0.279 0.227 0.121 Crown Radius Species DBH Height Crown Depth Basal area Stems ha-1 BALHT QMD RD Aspen 0.563 0.310 0.404 0.019 0.132 -0.125 0.119 0.170 Douglas-fir 0.795 0.649 0.677 0.009 -0.335 -0.596 0.248 -0.085 Lodgepole Pine 0.647 0.315 0.399 0.018 -0.158 -0.195 0.189 -0.029 Hybrid spruce 0.648 0.322 0.595 -0.022 -0.273 -0.208 0.203 -0.079 17 Following this, Eq`s 2.5 and 2.6 were transformed to linear equivalents using logarithms, thus enabling the use of variable selection tools. Maximum R2 improvement was used to obtain the best two-variable model, three-variable model, and so on. Additionally, variance inflation factors (VIF) were used to assess multicollinearity (Kutner et al. 2005). Based on the relationships among variables, maximum R2 variable selection, and biological theory, a final set of variables was chosen for Eq’s 2.5 and 2.6. Parameter estimates for Eq’s 2.5 and 2.6 were obtained for aspen, Douglas-fir, lodgepole pine and hybrid spruce using simultaneous regression methods. The model dataset used to fit the equations included observations pooled from both BEC zones. 2.2.4 Parameter estimation using simultaneous regression A system of simultaneous linear equations can be described as (using the notation of LeMay, 1990): γ β ε= + +i i i i i iY y X (2.7) where iY is an n by 1 matrix of n observations for the ith endogenous (i.e., dependent) variable on the left-hand side of an equation in the system of g equations, iy is an n by g-1 matrix of the observations for the g endogenous variables, excluding the ith endogenous variable, iγ is a g - 1 by 1 matrix of coefficients associated with the endogenous variables on the right-hand side of the equation, iX is an n by k +1 matrix of the observations for all of the exogenous (i.e., predictor) variables of the system plus a column of 1’s for the intercept, i β is a k +1 by 1 matrix of coefficients associated with the exogenous variables and the intercept, and i ε is an n by 1 matrix of the errors associated with the ith endogenous variable of the system. For each equation of the system, it is assumed that i ε has zero mean (E ( i ε )= 0), observations are independent, and variances are equal over all observations, resulting in the following variance-covariance matrix for the errors: nii ii ii ii ii Iσ σ σ σ = =Ω L MOM L L 00 0 00 00 18 where iiσ denotes the variances, which are all the same for all n observations. With a system of equations, error terms may be correlated across equations. As a result, the error-covariance matrix of the system of two equations becomes: n nn nn I II II ⊗ = = ΩΩ ΩΩ =Σ 2212 1211 2212 1211 2212 1211 σσ σσ σσ σσ where 11σ and 22σ are the variances for the first and second equations of the system, respectively, 12σ is the cross-equation covariance of the error terms, and ⊗ is the Kronecker product. As a result, efficiencies can be gained by fitting the equation as a system. For the crown allometry models (two equations), the simultaneous system of linear equations given in Eq. 2.7 can be modified to a system of nonlinear equations generally described as (Amemiya 1977): i i i i iY = f (y ,X ,θ ) (2.8) where if is a nonlinear function of endogenous and exogenous variables on the right-hand side of the equation, and i θ is vector of parameters to be estimated. The problem of simultaneous equation bias exists in the nonlinear form as in the linear form, and the use of cross-equations correlations can improve efficiency. Since the models chosen to provide predictions of crown depth and crown radius are nonlinear with respect to the parameters, N2SLS and N3SLS (Amemiya 1974; 1977) methods were used to estimate the parameters of the system of models. Parameter estimation is achieved through two stages in the case of N2SLS and three stages in N3SLS. In both methods, the first stage of parameter estimation involves the use of instrumental variables to estimate the endogenous variable on the right-hand side of all equations of the system (Hillier 2006). The estimated endogenous variables from the first stage replace their counterparts that occur on the right-hand side of equations in the system, thereby removing simultaneity bias. 19 For the crown allometry system of equations, estimated crown radius would replace measured crown radius and be used a predictor variable for the crown depth equation. Similarly, estimated crown depth from the first stage would replace crown depth where it occurs as a predictor variable for the crown radius equation. In essence, the estimation of the endogenous variable through the use of the instrumental variables removes the simultaneity bias caused by correlations between the error terms and endogenous variables on the right hand side of equations in the system. General guidelines for instrumental variables include selecting: 1) variables that are independent of the equation errors, such as the exogenous variables of the system of equations, or variables not included in the system but part of the model dataset; 2) low order polynomials of these aforementioned variables; and 3) derivatives with respect to the parameters, assuming that the derivatives are independent of the errors. At a minimum, the number of instrumental variables should equal the number of exogenous variables used in any one equation of the system. Commonly, OLS is used in the first stage to fit a linear function to estimate each endogenous variable from the instrumental variables (labeled as Z) as follows: -1 ˆ ( )β ′ ′=k k k k iZ Z Z y (2.9) kki Zy βˆˆ = (2.10) where kZ is the n x k matrix of instrumental variables; and kβˆ a k x 1 matrix of coefficients. In the second stage of N2SLS and N3SLS, the system of equations is fitted, where any endogenous variable on the right-hand side is replaced by estimates from the first stage equation, resulting in: -1( , , ) ( ) ( , , )θ θ′ ′ ′ ⊗ nf y x I Z Z Z Z f y x (2.11) where nI is an n by n identity matrix, and θ is the second stage estimator that minimizes the objective function, the sum of squared errors for the system of equations. In the second stage fit, errors across equations of the system are considered independent, resulting in the use of the identity matrix. 20 To improve efficiency by accounting for cross-equations correlations of error terms, the third stage of N3SLS replaces the identity matrix of Eq. 2.11 with a consistent estimate of Σ , defined here as the S matrix. The S matrix is constructed from the residuals following the second stage to obtain estimated variances of errors for each equation and covariances of errors across equations. The third and final stage of N3SLS uses generalized least squares to findθ that minimizes the objective function: -1 -1( , , ) ( ) ( , , )θ θ′ ′ ′ ⊗ f y x S Z Z Z Z f y x (2.12) where S is an estimate of Σ . It can be shown that the N3SLS estimator is more efficient than the N2SLS estimator, since the information in any cross-equation correlations is used in obtaining the N3SLS fit. The gain in efficiency decreases as the error variance-covariance matrix, Σ , approaches the identity matrix (i.e., Eq. 2.12 approaches Eq. 2.11). N2SLS and N3SLS carry the advantage that no assumption as to the probability distribution is used in the least squares fits. Therefore, the system of crown allometry models was fit using N2SLS and N3SLS using SAS and the PROC MODEL (SAS Institute Inc. 2003) procedure. A Gauss search algorithm was used to find the estimates of θ that minimized the N2SLS and N3SLS objective function. To seed the search algorithm, a set of starting parameters was provided for parameters in the simultaneous system of crown allometry models. These starting parameters were obtained by fitting single equation linear approximations of the nonlinear regression models using the exogenous variables, but not the dependent regressors on the right- hand side of the equations. To test whether the first stage of the N2SLS and N3SLS fits removed simultaneity biases, a Hausman specification test was used (Wu, 1973). This test is based on an m-statistic, where m is the total number of parameters estimated in the system of equations, and is used to test the null hypothesis of no measurement error in the variables on the right-hand side of the equations. The m-statistic follows a 2χ distribution under the null hypothesis; the 2χ was compared to a critical value from a 2χ table with m degrees of freedom. To test if there was any gain in efficiency in using N3SLS versus N2SLS, estimated standard errors of the coefficients were compared. Smaller standard errors following N3SLS estimation would suggest cross-equation correlations among equations following N2SLS. 21 2.2.5 Accuracy of Fitted Crown Models Accuracies of the crown allometry models were assessed by calculating fit statistics for the N2SLS and N3SLS estimators for each equation in the system. Estimated errors (residuals) for crown depth and radius were summarized to obtain mean bias and root mean square error (RMSE) defined as: Mean Bias ∑ = − = n i ii n YY 1 )ˆ( (2.13) and, RMSE ∑ = − = n i jj n YY 1 2)ˆ( (2.14) where jY is the actual value (i.e., measured crown depth or crown radius) for measurements j to n, jYˆ is the predicted value from the fitted equation; and n is the number of trees. A pseudo R 2 statistic, appropriate for nonlinear models, was calculated as: Pseudo – R2 ( ) ( )∑ ∑ = = − − −= n i j n i jj YY YY 1 2 1 2 ˆ 1 (2.15) as an additional diagnostic (Schabenberger and Pierce, 2002). The mean bias and RMSE values were calculated for all trees of the species, and then separately for each of four density classes with a class width of 500 stems ha-1. This was also repeated for basal area, using 10 m2 ha-1 intervals. To examine the impact of the stand density variables in the equations, partial-F statistics were calculated by first fitting the system of equations represented by Eq. (2.5) and Eq. (2.6) using only the tree-size variables (DBH, height, crown depth, and crown radius). Density variables were then added and results were compared to the fit using tree-level variables only. 2.2.6 Model Evaluation and Validation Since Eq’s 2.5 and 2.6 are nonlinear, the estimated coefficients are asymptotically unbiased and consistent only if sample sizes are large (Kutner et al. 2005). To assess these properties for 22 the sample size used in this study, bootstrap resampling was used to re-estimate the coefficients and their standard errors for the system of equations. The models were refit 1000 times using the same regression methods employed in the original fit. Large differences between the coefficients obtained in the original fit of the system of equations and those obtained through bootstrapping suggest that a larger sample size is needed. For a detailed description of general bootstrap sampling methods see papers by Efron (1987) and DiCiccio and Efron (1996). To assess the predictive capabilities of the chosen models, the predicted sum of square statistic (PRESS) (Quan 1988) commonly used to validate single equation models was extended to the system of equations used in this study. To calculate the PRESS statistic for the system of equations, the ith observation from the model dataset was deleted, then, the remaining observations were used to refit the models using a N3SLS estimator. Refitting of the system of equations in this manner was repeated n times, where n is sample size of the model dataset. Each time the coefficients for the system of equations were re-estimated, a predicted value, )(ˆ jjY , was obtained for the jth observation that was deleted, where )( jj is used to denote that predicted value is for the jth observation that was deleted. The PRESS statistic was calculated from the sum of squared prediction errors over all n observations: ( )∑ = −= n i jjj YYPRESS 1 2 )( ˆ (2.16) where jY is the jth deleted observation. The PRESS statistic and was compared to the error sum of squares (SSE) obtained in the original fit of the system of equations. PRESS values close to SSE support to the internal validity of the model and suggest that the means square error (MSE) is a reasonably good indicator of the model’s predictive abilities (Soares and Tome 2001). 2.3 Results 2.3.1 System of Crown Allometry Equations For the system of equations, DBH (cm), total tree height (m) and the ratio of height to DBH (H/D) (i.e., slenderness coefficient) were selected as predictor variables to account for the effects of tree size on crown radius and crown depth. By definition, the simultaneous system of equations also included crown radius and crown depth as predictor variables; crown radius used on the right-hand-side of the equation in the crown depth model and crown depth used on the right- 23 hand-side of the equation in the crown radius model. For the study area, the effects of density on crown size were best described through a combination of stem ha-1, basal area (m2 ha-1) and BALHT. BALHT was preferred over the more commonly used basal area of trees larger by diameter, as it was regarded to be a better surrogate for the effects of shade cast by neighbouring trees (Vanclay et al. 1995; Schwinning and Weiner 1998; Coomes and Allen 2007). From the list of variables selected for use in the system of equations, DBH, DBH2, 1/DBH, height, height2, H/D, BALHT, stems ha-1 and basal area (m2 ha-1) were used as instrumental variables in the first stage equations for N2SLS and N3SLS. These first stage equations were: DBH HeightTPHBABALHT DBH 1HeightDBHHeightDBHdâCr 1111 1 2 1 2 1111 ×+×+×+×+ ×+×+×+×+×+= jihg fedcba (2.17) and DBH HeightTPHBABALHT DBH 1HeightDBHHeightDBHthˆC 2222 2 2 2 2 2222 ×+×+×+×+ ×+×+×+×+×+= jihg fedcba (2.18) where dâCr is the first-stage estimated crown radius, and thˆC is the first-stage estimated crown depth, and 1a to 1j and 2a to 2j are sets of species-specific parameters for the first-stage equations. Estimates of first-stage coefficients are given in Table 2.3. 24 Table 2.3 Estimated coefficients and their standard errors (in parentheses) for first-stage equations. Variable Crown Depth Crown Radius Aspen Douglas-fir Pine Spruce Aspen Douglas-fir Pine Spruce Intercept 27.11 (22.5) 7.57 (9.1) 3.60 (5.8) 19.5 (14.4) -2.72 (2.86) 0.95 (0.729) 1.07 (0.394) 0.652 (0.95) DBH -2.47 (1.8) 1.21 (0.5) 0.02 (0.4) -2.59 (1.19) 0.21 (0.23) -0.047 (0.05) 0.019 (0.029) 0.005 (0.08) DBH2 0.05 (0.1) -0.01 (0.01) 0.005 (0.01) 0.039 (0.017) -0.003 (0.01) 0.001 (0.001) 0.001 (0.001) 0.001 (0.001) 1/DBH -102.1 (96.4) -84.68 (61.5) -4.93 (28.10) 25.5 (68.99) 13.3 (12.2) 3.71 (4.95) -4.98 (1.91) 1.09 (4.51) Height 1.31 (1.14) -2.14 (1.26) 0.254 (0.38) 4.687 (1.54) 0.019 (0.15) 0.129 (0.10) -0.071 (0.03) 0.077 (0.101) Height2 -0.014 (0.02) 0.038 (0.02) -0.012 (0.01) -0.086 (0.03) -0.0002 (0.003) -0.002 (0.002) 0.001 (0.001) -0.001 (0.002) H/D -0.710 (6.7) 13.25 (9.69) 0.509 (2.67) -29.62 (11.13) 0.197 (0.86) -1.19 (0.78) 0.339 (0.18) -0.831 (0.73) BALHT 0.056 (0.03) 0.014 (0.03) -0.001 (0.01) -0.006 (0.019) 0.007 (0.01) -0.002 (0.003) -0.0001 (<0.000) 0.001 (0.001) Basal area -0.341 (0.15) -0.124 (0.09) -0.002 (0.04) -0.03 (0.080) -0.033 (0.02) -0.0014 (0.01) 0.001 (0.003) 0.004 (0.01) Stems ha-1 -0.003 (0.01) 0.002 (0.003) -0.001 (<0.00) -0.001 (0.002) 0.001 (0.001) 0.0001 (<0.000) -0.0001 (<0.000) -0.001 (0.001) For the system of equations fit in the second- and third-stage, DBH, height, dâCr , thˆC , BALHT, stems ha-1 and basal area (m2 ha-1) were used. The system of equations was therefore: 3 3 3 3 3 3 3 3 3 3 CrownDepth 1 exp( ( )) ˆ( ) = DBH Height Crad BALHT BA TPH Height f x f x a b c d e f g ε= + + + × + × + × + × + × + × (2.19) and 4 4 4 4 4 4 4 ˆCrownRadius DBH Height Cht BALHT BA TPHrb c d e f ga ε= × × × × × × + (2.20) where 3a to 3g and 4a to 4g are sets of species-specific parameters, and 3ε and 4ε are error 25 terms. A comparison of N2SLS and N3SLS estimated coefficients and approximate standard errors is shown in Tables 2.4 and 2.5. The coefficients for BALHT and dâCr showed notable differences for estimates of aspen crown depth. For BALHT, the coefficient switched from a small negative value in N2SLS to a small positive value following N3SLS. For dâCr , the coefficient dropped from -0.407 to -1.014. For estimates of Douglas-fir crown radius, the coefficient for height switched from a negative value following N2SLS to a positive value following N3SLS. For lodgepole pine, the coefficients obtained through N2SLS and N3SLS remained fairly similar with the exception of DBH. For this variable, the coefficient changed from -0.19 following N2SLS to 0.20 with N3SLS. Table 2.4 Estimated coefficients and their standard errors (in parentheses) from N2SLS and N3SLS fits of the crown depth model in system of equations. Variable Second Stage Estimates Third Stage Estimates Aspen Douglas-fir Pine Spruce Aspen Douglas- fir Pine Spruce Intercept -0.382 (0.78) -5.359 (3.57) -0.989 (0.66) -2.546 (1.15) 0.041 (0.74) -4.940 (2.41) -0.625 (0.65) -2.491 (1.09) DBH 0.060 (0.05) -0.106 (0.10) -0.019 (0.07) -0.064 (0.11) 0.092 (0.05) -0.115 (0.07) 0.020 (0.07) -0.074 (0.10) Height -0.072 (0.05) 0.056 (0.08) 0.138 (0.04) 0.172 (0.09) -0.064 (0.05) 0.075 (0.05) 0.115 (0.04) 0.163 (0.08) Crown Radius -0.407 (0.92) 1.99 (2.02) -0.485 (0.89) -0.561 (1.06) -1.014 (0.84) 1.945 (1.30) -1.003 (0.92) -0.459 (0.94) BALHT -0.004 (0.01) 0.010 (0.02) 0.001 (0.001) 0.006 (0.01) 0.001 (0.01) 0.006 (0.01) 0.001 (0.001) 0.004 (0.01) Basal area 0.038 (0.06) 0.056 (0.05) -0.009 (0.01) 0.006 (0.03) 0.012 (0.06) 0.053 (0.03) -0.007 (0.01) 0.015 (0.03) Stems ha-1 0.001 (0.001) 0.001 (0.001) 0.001 (0.001) 0.001 (0.001) 0.002 (0.001) <0.000 (0.001) 0.001 (0.001) 0.001 (0.001) 26 Table 2.5 Estimated coefficients and their standard errors (in parentheses) from N2SLS and N3SLS fits of the crown radius model in system of equations. Variable Second Stage Estimates Third Stage Estimates Aspen Douglas-fir Pine Spruce Aspen Douglas- fir Pine Spruce Intercept 0.134 (0.17) 0.458 (0.69) 0.442 (0.15) 0.438 (0.30) 0.069 (0.09) 0.873 (1.15) 0.487 (0.16) 0.415 (0.28) DBH 0.58 (0.15) 0.733 (0.39) 0.932 (0.1) 0.811 (0.21) 0.603 (0.14) 0.942 (0.25) 0.948 (0.09) 0.760 (0.19) Height -0.94 (0.65) -0.115 (0.23) -0.483 (0.07) -1.00 (0.23) -1.125 (0.62) 0.065 (0.16) -0.499 (0.07) -1.057 (0.23) Crown Depth 0.49 (0.28) -0.227 (0.38) 0.145 (0.15) 0.739 (0.35) 0.650 (0.28) -0.663 (0.20) 0.109 (0.15) 0.797 (0.34) BALHT -0.032 (0.11) -0.010 (0.01) -0.002 (0.02) 0.024 (0.02) -0.028 (0.10) -0.006 (0.01) -0.007 (0.02) 0.017 (0.01) Basal area 0.163 (0.30) -0.165 (0.23) 0.029 (0.08) 0.201 (0.11) 0.195 (0.28) -0.389 (0.18) 0.044 (0.07) 0.236 (0.10) Stems ha-1 0.359 (0.27) 0.096 (0.20) -0.093 (0.05) -0.120 (0.11) 0.472 (0.28) 0.075 (0.18) -0.103 (0.05) -0.099 (0.10) The Hausman specification test indicated that the calculated m-statistic was less than the critical value from a 2χ distribution with 14 degrees of freedom and α = 0.05 for all four species,. Thus, the null hypothesis of no measurement error in the variables on the right-hand side of the equations was not rejected. Both N2SLS and N3SLS estimators benefitted from the use of the instrumental variables in the first-stage through the removal of simultaneity bias. For nearly all of the estimated coefficients, there was a gain in efficiency when the N3SLS was used, as evidenced by smaller standard errors. Due to the gain in efficiency obtained through the third- stage of parameter estimation, the N3SLS estimator was preferred over the N2SLS estimator. An evaluation of fit statistics (Mean bias, RMSE and Pseudo- R2) using the N3SLS estimates suggested that fitting the models for crown depth and crown radius in a system of equations provided good estimates of the two crown axes over a range of stand densities. This was most evident in the estimated mean biases for crown depth and crown radius, which were low for all four species (Table 2.6). There was a general tendency for predictions of crown depth and crown radius to be slightly lower than values measured in the field on average. The one exception was the estimate of crown depth for aspen, which was larger than mean value recoded in the field on average. The largest mean bias was observed for estimates of crown depth on hybrid spruce. The estimated standard deviations (RMSE) were generally small for both crown radius and crown depth, and suggested that the fitted system of equations provided accurate estimates. For 27 lodgepole pine and hybrid spruce, estimates of crown depth differed from actual values by roughly 2 m, while estimates of crown radius differed from actual values by nearly 0.3 m. The largest estimates of RMSE were observed for Douglas-fir, where estimates of crown depth differed from actual values by 3.4 m on average, while estimates of crown radius differed by 0.57 m from actual values on average. Table 2.6 Mean bias (observed - predicted), Root Mean Squared Error (RMSE) and Pseudo- R2 results using N3SLS parameter estimates. Dependent variable Bias (m) RMSE (m) Pseudo- R2 Aspen Crown Radius 0.001 0.329 0.419 Crown Depth -0.004 1.401 0.624 Douglas-fir Crown Radius 0.001 0.570 0.228 Crown Depth 0.014 3.393 0.284 Lodgepole pine Crown Radius 0.001 0.285 0.519 Crown Depth 0.001 2.095 0.156 Hybrid spruce Crown Radius 0.003 0.331 0.532 Crown Depth 0.029 2.054 0.716 When the fit statistics were calculated by basal area (m2 ha-1) and stems ha-1 classes, a clearer picture of the behaviour of the system of crown models emerged (Table 2.7). Generally, biases within each basal area and stems ha-1 class were low among the four species tested and tended to be similar to biases calculated in Table 2.6. For all four species, there were no systematic patterns of over- or under-estimation of crown depth or crown radius among the different density classes. Consider, as an example, the seven mean biases of crown radius calculated for basal area classes 25 m2 ha-1 and 35 m2 ha-1 for the four species. Out of these seven mean biases, four were positive while three were negative. Then, looking at the seven mean biases of crown radius calculated for basal area classes 5 m2 ha-1 and 15 m2 ha-1, a similar proportion of mean biases show over- and under-estimation. This balance of over- and underestimation among dense and open stands suggests that systematic patterns of bias were eliminated. Conversely, the fitted standard crown allometry equation (Eq 2.1) tended to overestimate crown radius in the two larger basal area classes. Here, six of the seven mean biases that were calculated for the four species in the 25 m2 ha-1 and 35 m2 ha-1 basal area classes were positive. Systematic biases associated with the standard allometry equations were particularly noticeable for lodgepole pine. For this species, 28 there was also a strong tendency to underestimate crown radius and crown depth among stands with a lower density. Table 2.7 Mean bias (observed - predicted) and RMSE (in parentheses) in metres by density class (stem ha- 1) and by basal area (m2 ha-1) class for N3SLS Stems ha-1 No. trees Mean Bias (RMSE) Basal area No. trees Mean Bias (RMSE) Crown Radius Crown Depth Crown Radius Crown Depth Aspen 250 6 0.062 (0.12) -0.023 (0.31) 5 6 0.062 (0.12) -0.023 (0.31) 750 30 -0.009 (0.31) -0.044 (1.35) 15 26 0.006 (0.31) -0.159 (1.38) 1250 4 -0.020 (0.47) 0.323(1.89) 25 8 -0.060 (0.35) 0.511 (1.37) 1750 0 - - 35 0 - - Douglas-fir 250 0 - - 5 0 - - 750 51 -0.021 (0.54) -0.091 (3.50) 15 18 0.029 (0.35) 0.393 (2.02) 1250 5 0.183 (0.49) 0.772 (2.64) 25 21 0.025 (0.68) -0.186 (4.37) 1750 0 - - 35 17 -0.067 (0.51) -0.233 (3.42) Lodgepole Pine 250 68 -0.106 (0.26) 0.306 (1.99) 5 95 -0.010 (0.26) 0.397 (1.43) 750 237 0.063 (0.29) -0.177 (2.20) 15 207 0.021 (0.30) -0.272 (2.26) 1250 131 -0.049 (0.24) 0.333 (1.83) 25 139 -0.031 (0.25) 0.242 (1.99) 1750 16 -0.102 (0.20) -0.575 (2.29) 35 11 0.074 (0.38) -0.164 (3.87) Hybrid Spruce 250 10 -0.031 (0.19) -0.028 (0.40) 5 8 -0.071 (0.14) -0.148 (0.37) 750 62 0.007 (0.34) 0.084 (2.26) 15 27 0.084 (0.25) -0.067 (1.11) 1250 6 0.081 (0.43) -0.439 (1.57) 25 26 -0.106 (0.34) 0.390 (1.69) 1750 0 - - 35 17 0.053 (0.43) -0.221 (3.69) Figures 2.2 and 2.3 show how the prediction biases for crown radius and crown depth of individual trees varied over four basal area classes (0-10, 10-20, 20-30 and 30-40 m2 ha-1) for the four species tested. 29 Figure 2. 2 Relationship between observed and predicted aspen crown depth (a) and crown radius (b), and Douglas-fir crown depth (c) and crown radius (d). a) c) b) d) 30 Figure 2.3 Relationship between observed and predicted lodgepole pine crown depth (a) and crown radius (b), and hybrid spruce crown depth (c) and crown radius (d). a) b) c) d) 31 To further examine the influence of the stand density measures, partial-F tests were used to test the null hypothesis that stand density variables did not contribute to the predictive abilities of the models used in the system of equations. Results of the partial-F tests suggested that the measures of stand density, in the presence of the tree level variables, tended to improve the predictive abilities of the system of equations. For estimates of crown depth, the F-test values for aspen, Douglas-fir, lodgepole pine and hybrid spruce, where 11.9, 16.6, 94.6 and 53.5, respectively, indicating that the stand density measures significantly improved the predictive abilities (all p values <0.05) of the model. The proportion of variability in crown depth explained by the stand density variables (partial-R2 values) for these were 0.52, 0.50, 0.39 and 0.08 for aspen, Douglas-fir, lodgepole pine and hybrid spruce, respectively. For estimates of crown radius, the density-related variables improved prediction for Douglas-fir (partial-F test = 4.85, p<0.05). Estimates of crown radius for aspen, lodgepole pine and hybrid spruce showed no obvious improvement with the inclusion of stand density measures. For these three species, the effects of competition from neighbouring trees appear to have been accounted for through the use of DBH, height and first stage estimated crown depth as predictor variables. 2.3.2 Model Evaluation and Validation Generally, there was close agreement between the mean of the coefficients obtained by bootstrapping the system of equations and those estimated in the original N3SLS fitting procedure (presented in Tables 2.4 and 2.5). This was particularly true for lodgepole pine and hybrid spruce. For example, the mean of the bootstrapped coefficients for DBH, height and crown depth used to estimate crown radius for lodgepole pine were 0.943, -0.475 and 0.127, respectively. In the original N3SLS fit for lodgepole pine, the estimated coefficients for DBH, height and crown depth were 0.948, -0.499 and 0.109, respectively. In fact, differences between original N3SLS and bootstrapped estimates of coefficients for lodgepole pine and hybrid spruce were rarely greater than 0.1, indicating that original estimators were nearly unbiased. For aspen and Douglas-fir, the means of the coefficients obtained through bootstrap resampling differed from the original N3SLS estimates by a larger amount. For example, the mean of the bootstrapped coefficients for tree height in the equation to predict crown depth in Douglas-fir was 0.339, whereas the original N3SLS estimate was 0.08. The difference suggests that the original estimators for coefficients are likely biased. Increasing the sample size would eliminate this bias as the coefficients for the nonlinear models used in the system of equations are asymptotically unbiased. 32 For aspen, Douglas-fir, lodgepole pine, and hybrid spruce, the crown radius model PRESS statistics were 4.87, 55.62, 37.88, and 11.95, respectively. Comparatively, the sum of squared error (SSE) values for these species obtained in the original fit were 3.58, 15.97, 36.19 and 8.64. The PRESS statistics related to the crown depth model were only marginally larger than the SSE of the original fit for aspen, Douglas-fir and lodgepole pine. For these species, the crown depth PRESS statistics were 94.79, 718.03 and 2082.91, respectively. For the same species, the SSE values were 64.73, 644.25 and 1981.21. For hybrid spruce, the difference was larger, with a value of 753.45 for the PRESS statistic versus an SSE value of 329.31 obtained in the original fit of the system of equations. 2.4 Discussion In the MPB-disturbed forests around Williams Lake BC, trees in dense stands tended to have narrower and shorter crowns, while trees in more open stands tended to have wider and longer crowns. Furthermore, a strong relationship existed between the width and the depth of tree crowns. To account for the effects of stand density, stems ha-1, basal area (m2 ha-1) and BALHT were used as predictor variables in a system of crown allometry equations. Given the presence of other tree size variables, these three measures of density significantly improved predictions of crown depth for the four species examined in this study. Marginal improvements were seen for estimates of crown radius. This suggests that for the study area, stand density plays an important role in determining the size of tree crowns. Thus, in stands with greater crowding, light attenuation from canopy trees may be partially offset by the density- related processes that limit crown size. Conversely, higher levels of light attenuation than expected may occur in open stands, since there are fewer density-related factors limiting the size of crowns. Because regeneration in the understory is often limited by the light environment, establishing the relationship between crown size and density is important from a regeneration modelling perspective. Estimates of regeneration in the forest dynamics model SORTIE-ND are largely a function of light; thus, the crown equations presented in this study are suitable alternatives to the standard crown equations currently in the model. In addition to the effects of density, the crown model described in this study was based on the principle that changes in crown radius and crown depth were integrally related since they are influenced by similar biological processes. As a result, using the opposite axis as a predictor variable on the right-hand-side of a model should significantly improve estimates of the axis on the left-hand-side of the model. Using a system of crown equations and a N3SLS estimator to 33 remove simultaneity bias and take advantage of cross-equation correlation, the models showed a good prediction of crown depth and crown radius for the four species tested. Furthermore, the structural form of the model for crown depth provided built-in assurances that the predictions fell within biologically attainable ranges by using tree height as an asymptote (Eq. 2.19). For the crown radius model, using a power model resulted in biologically realistic predictions. The crown models in SORTIE-ND use only a single tree size variable to predict tree crowns; fitting the standard crown equations to the model dataset resulted in systematic patterns of over- and underestimation among dense and open stands, respectively. To correct for the density- related biases, Eq’s 2.19 and 2.20 explicitly modelled stand density through the use of stems ha-1, basal area (m2 ha-1), and BALHT. From a biological point of view, these simple measures of density each capture a potentially different effect of density on crown depth and radius. Including BALHT as a measure of density provided the models with the ability to capture the effects of shading on the structure of tree crowns, as only taller trees can intercept light before it reaches the target tree (Coomes and Allen 2007). The inclusion of BALHT seemed particularly warranted since the intended application of the crown models was to the light mediated growth model, SORTIE-ND. Both basal area and stems ha-1 in the crown models were included since these variables respond differently as stands age. In the absence of major disturbances, basal area is typically reallocated to larger diameter classes with increasing stand age. Over this same period, stems ha-1 tends to decrease as a result of mortality among older trees and the inability of many understory trees to withstand low levels of light (West 1995; Dean 2004; Keeton 2006). The rate at which basal area increases and stems ha-1 decreases may be highly variable and is often related to factors such as species composition and nutrient availability. Thus, the strength of the relationship between the crown axes and stems ha-1 and basal area will vary with succession stage. From a modelling point of view, it makes sense to include basal area and stems ha-1 as predictor variables for two reasons. First, they corrected the systematic biases in the fitted standard equations observed when the prediction bias was broken down by stems ha-1 and basal area classes. Second, using measures of density that respond differently to forest succession should enable the models to be applied to a wider range of forest conditions. Direct support for this second reason was not obtained in this study. Evidence to support or reject this reason needs to come from external validation of the chosen models. It was somewhat surprising to note that the predictive abilities of the crown radius model were not significantly improved when BALHT, stems ha-1and basal area were included as predictor variables. Prior to fitting the equations, it was hypothesized that the radius of crowns 34 would be smaller in dense stands due to limited growing space and direct physical contact from neighbouring trees. In open stands, crowns would be wider, taking advantage of the available growing space. Furthermore, it was anticipated that stands thinned by repetitive attacks of MPB would allow for trees that were not attacked to expand their crowns as surrounding trees fell. Failure to reject the null hypothesis posed in the partial-F tests suggests that other mechanisms are at work. It is possible that fitting the system of equations with a reduced model consisting of only DBH, tree height and crown depth captured most of the variability in crown radius. For the model dataset, DBH, height and first stage estimated crown depth were moderately to strongly correlated with stems ha-1 and basal area (m2 ha-1) for all four species. Thus, they may each act as a surrogate measure of density. If, as a group, they were strongly related to stand density, then this would explain why adding three explicit measures of stand density failed to significantly improve the predictive abilities of the crown model. It is also possible that the reduction in density that followed MPB attack 25 years ago did not open the stands enough to allow the radius of crowns to significantly increase in size. It is estimated that 25-years after MPB infestation, the average reduction of overstory basal area was 36%. Examining stands where the intensity of attack was higher would be helpful in determining how crowns respond following MPB- disturbance. For aspen, Douglas-fir, lodgepole pine and hybrid spruce in the model dataset, it appeared that crown depth and radius were highly interdependent. This strong relationship has been explained through physiological processes for several tree species (Williams et al. 1999). Most simply, changes in crown depth and radius are a result of the terminal buds exerting control over the depth and orientation of lateral branches, a phenomenon known as epinastic control (Oliver and Larson 1996). The level of control exerted by the terminal bud tends to be stronger among shade-tolerant conifers such as spruce, when access to indirect sunlight is available. Most of the hybrid spruce trees recorded in the model dataset would have had access to direct sunlight following the loss of foliage and fall-down of pine trees infested by MPB. This would have allowed the terminal buds of hybrid spruce trees to assert control over the overall shape of the crown. This pattern of strong excurrent growth form likely contributed to the precise estimates obtained using the fitted system of crown equations. For shade-intolerant species, epinastic control tends to be weaker. It was expected that this weaker relationship among shade intolerant species would translate to poor fits of the model for lodgepole pine and aspen in the model dataset. However, this was not the case, as good fit statistics were achieved for both species. Although the main advantage of fitting the crown models as a system was to take advantage 35 of the inherently strong relationship between crown radius and depth, the approach used to fit the equations also provided some good statistical properties. Using a three-stage least squares estimator for Eq’s 2.19 and 2.20 provided coefficients that were both asymptotically efficient and consistent. Simultaneity biases that could have resulted from using crown depth and radius on the right-hand-side of the crown equations were purged through the use of two first-stage equations, one to estimate crown depth and one to estimate crown radius. The estimated crown radius and depth from the first-stage equations were then passed to the second-stage of the simultaneous fitting process. The instrumental variables in the first-stage equations employed transformations of the tree size and stand density variables used in the second- and third-stage of model fitting, plus the slenderness coefficient (i.e., height/DBH). Adding the slenderness coefficient improved the estimates of crown radius and crown depth in these first-stage equations. Trees with a high slenderness coefficient tend to bend to a greater degree under wind stress, increasing the chances of crown abrasion and influencing crown size (Meng et al. 2007). An examination of the coefficients and their standard errors revealed that some efficiency was gained when switching from a N2SLS to a N3SLS estimator, and the Hausman test indicated that simultaneity biases were removed through the first-stage equations. The resampling methods used to assess the bias in the estimated coefficients provided encouraging results. Since nonlinear models were selected, only asymptotically unbiased and consistent estimates from N2SLS and N3SLS were possible. This was in doubt for aspen, Douglas-fir and hybrid spruce. For hybrid spruce, coefficients obtained through bootstrapping were in close agreement with those estimated in the original model fit. Thus, sample sizes for spruce were large enough to support large sample theory. For aspen and Douglas-fir, differences between the bootstrapped coefficients and those estimated in the original model fits suggested that more crown measurements for these species are needed. For now, the recommendation is that the coefficients reported in Tables 2.4 and 2.5 are useful for the study area, but should be used with some caution. In terms of validation, the PRESS statistics for crown radius suggested that for aspen, lodgepole pine and hybrid spruce, the system of equations had good predictive abilities. The same conclusions were reached regarding estimates of crown depth for aspen, Douglas-fir and lodgepole pine. These results provided added support to the internal validity of the fitted regression model (Kutner et. al, 2005). However, for estimates of crown radius on Douglas-fir and crown depth on hybrid spruce, it appears that the system of equations was somewhat unreliable. Further evaluation of the model using a true validation dataset may shed light on 36 where future improvements to the system of equations could be made to ensure that estimates of both crown axes are reliable for the most common tree species within the study area. 2.5 Conclusions The standard crown allometry models used by SORTIE-ND are based on the structural relationship between the size of the tree bole and crown size. The density-dependent crown models presented in this study go beyond this basic relationship on two levels. First, by using crown depth and crown radius as dependent regressors, the models reflect the dual causality that arises between the two crown axes, which is a result of physiological and mechanical processes. Second, the models explicitly account for the effects of crowding on crown size. The impetus to develop the density-dependent crown allometry models came partially from the need to have a set of models suitable for use in stands that have been disturbed by MPB. The assumption was that following MPB disturbance, crowding among overstory trees decreases, allowing crowns on surviving trees to increase in size. The accurate predictions of crown depth and crown radius that were obtained by fitting the density dependent crown models using data from stands that have been attacked by MPB support this assumption. For these stands in the Williams Lake study area, lodgepole pine is by far the most dominant species by proportion, with hybrid spruce, Douglas-fir and aspen generally forming minor components. Consequently, any of the noted deficiencies related to the estimates of crown depth and crown radius of Douglas-fir and aspen would likely have minimal impact in terms of changes to the total crown area for a stand. At least for lodgepole pine and hybrid spruce, the objective of obtaining accurate estimates of crown depth and radius for a range of stand densities for unmanaged, MPB-disturbed stands in the Williams Lake area of BC appears to have been met. However, the ramifications of using the new equations within SORTIE-ND remain unknown. Thus, the next step is for the equations and the species-specific coefficients to be tested within SORTIE-ND using data from the Williams Lake study area. Once the density dependent crown models have been added to SORTIE-ND, their influence on the light submodel will be assessed through an evaluation of Global Light Index (GLI) values. It will be the subject of further study to determine if the predicted changes in GLI are reflective of the actual light levels for these stands and how the predicted light-mediated establishment and growth of natural regeneration within SORTIE-ND will be altered. Until this next level of testing is complete, it would be difficult to describe a broader set of situations or objectives for which use 37 of the density-dependent equations is appropriate. For now, the density-dependent equations present an alternative to the existing standard crown allometry equations when estimates for MPB-disturbed stands are of interest. 2.6 Literature cited Amemiya, T. 1974. Multivariate regression and simultaneous equation models when dependent variables are truncated normal. Econometrica 42: 999-1012 Amemiya, T. 1977. Maximum likelihood and nonlinear 3-stage least-squares estimator in general nonlinear simultaneous equation model. Econometrica 45: 955-968. Assmann, E. 1970. The Principles of Forest Yield Studies. Pergamon Press, Oxford, 506 pp. Astrup, R. 2006. Modeling growth of understory aspen and spruce in western boreal Canada. Ph.D. Dissertation. University of British Columbia, Vancouver, BC. 156 p. Astrup, R., Coates, K. D, Hall, E., and Trowbridge, A. 2007. Documentation of the SORTIE-ND SBS research parameter file version 1.0. February, 2007. Bulkley Valley Centre for Natural Resources Research and Management, Smithers, B.C.. Aukema, B. H., Carroll, A. L., Zhu, Jun., Raffa, K. F., Sickley, T.A., and Taylor, S.W. 2006. Landscape level analysis of mountain pine beetle in British Columbia, Canada: spatiotemporal development and spatial synchrony within the present outbreak. Ecog. 29: 427-441. Baldwin, V.C., and Peterson, K.D. 1997. Predicting the crown shape of loblolly pine trees. Can. J. For. Res. 27: 102-107. Bechtold, W.A. 2003. Crown-diameter prediction models for 87 species of stand-grown trees in the eastern United States. South. J. Appl. For. 27: 269-278. Biging, G.S., and Gill, S.J. 1997. Stochastic models for conifer tree crown profiles. For. Sci. 43: 25-34. Borders, B. E., and Bailey, R. L. 1986. A compatible system of growth and yield equations for slash pine fitted with restricted three-stage least squares. For. Sci. 32: 185-201. 38 Canham, C.D., Coates, K.D., Bartemucci, P., and Quaglia, S. 1999 Measurement and modeling of spatially explicit variation in light transmission through interior cedar-hemlock forests of British Columbia. Can. J. For. Res. 29: 1775-1783. Canham, D. C., Finzi, A. C., Pacala, S. W., and Burbank, D. H. 1994. Causes and consequences of resource heterogeneity in forest - interspecific variation in light transmission by canopy trees. Can. J. For. Res. 24: 337. Coates, K.D., Canham, C.D., Beaudet, M., Sachs, D.L., and Messier, C. 2003. Use of a spatially explicit individual-tree model (SORTIE/BC) to explore the implications of patchiness in structurally complex forests. For. Ecol. Manage. 186: 297-310. Coomes, D. A. and Allen, R. B. 2007. Effects of size, competition and altitude on tree growth. J. Ecol. 95: 1084-1097. Curtis, R. O. 1982. A simple index of stand density for Douglas-fir. For. Sci. 28: 92-94. Dean, T.J. 2004. Basal area increment and growth efficiency as functions of canopy dynamics and stem mechanics. For. Sci.. 50: 106-116. DiCiccio, T. J. and Efron, B. 1996. Bootstrap confidence intervals. Stat. Sci. 11: 189-228 Efron, B. 1979. Bootstrap methods: Another look at the Jacknife. Ann. Stat. 7: 1-26. Efron, B. 1987. Better bootstrap confidence-intervals. Amer. Stat. Ass. 82: 171-185 Gallant, A. R. 1975. Seemingly unrelated nonlinear regressions. J. Econo. 3: 35-50. Gilmore, D.W. 2001. Equations to describe crown allometry of Larix require local validation. For. Ecol. & Manage. 148: 109-116. Hasenauer, H., and Monserud, R.A. 1996. A crown ratio model for Austrian forests. For. Ecol. & Manage. 84: 49-60. Hasenauer, H., Monserud, R. A. and Gregoire, T. 1997. Using simultaneous regression techniques with individual-tree growth models. For. Sci. 44: 87-95. Hawkes, B., Taylor, S., Stockdale, C., Shore, T.L., Alfaro, R.I., Campbell, R., and Vera, P. 2004. Impacts of mountain pine beetle on stand dynamics in British Columbia. pp. 175-195 in 39 T.L. Shore, J.E. Brooks, and J.E. Stone, editors. Mountain Pine Beetle Symposium: Challenges and Solutions. Kelowna, British Columbia October 30-31, 2003, Nat. Res. Can., Can. For. Serv., Pac. For. Centre, Inform. Rep. BC-X-399, Victoria, BC. 298pp. Heath, R. and R.I. Alfaro. 1990. Growth response in a Douglas-fir/lodgepole pine stand after thinning of lodgepole pine by the mountain pine beetle: a case study. Journal of the Entomological Society of British Columbia 87:16-21. Hillier, G. 2006. Yet more on the exact properties of IV estimators. Econ. Theor. 22: 913-931. Hynynen, J. 1995. Predicting tree crown ratio for unthinned and thinned Scots pine stands. Can. J. For. Res. 25: 57-62. Huang, S. and Titus, S. J. 1999. Estimating a system of nonlinear simultaneous individual tree models for white spruce in boreal mixed-species stands. Can. J. For. Res. 29: 1805-1811. Judge, G. G., Griffiths, R. C., Lutkepohl, H. and Lee, T. C. 1985. The theory and parctic of econometrics. 2nd ed. John Wiley and Sons, New York. Keeton, W. S. 2006. Managing for late-successional/old-growth characteristics in northern hardwood-conifer forests. For. Ecol. & Manage. 235: 129-142. Kellomaki, S., Ikonen, V. and Kolstrom, T. 1999. Modelling the structural growth of Scots pine with implications for wood quality. Ecol. Modell. 122: 117-134. Kutner, M.H., Nachtsheim, C.J., Neter, J., and W. Li. 2005. Applied Linear Statistical Models. International Edition. McGraw-Hill, New York. pp. 214-232. LeMay, V. M. 1990. MSLS: A linear least-squares technique for fitting a simultaneous system of equations with a generalized error structure. Can. J. For. Res. 20: 1830-1839. LeMay, V.M.; Lee, T.; Scott, R.E.; Sattler, D.F; Robinson, D.; Zumrawi, A-A.; Marshall, P. 2006. Modeling Natural Regeneration Following Mountain Pine Beetle Attacks in the Southern and Central Interior of British Columbia: Results for Year 1. Internal report for Nat. Res. Can., MPB Standard Contribution Agreement, PO # 8.35. 70 pp. http://www.forestry.ubc.ca/prognosis/extension.html 40 LeMay, V., T. Lee, D. Sattler, P. Marshall, D. Robinson, and A-A. Zumrawi. 2007. Modeling Natural Regeneration Following Mountain Pine Beetle Attacks in the Southern and Central Interior of British Columbia: Final report. Internal report for Natural Resources Canada, MPB Standard Contribution Agreement, PO # 8.35. 57 pp. http://www.forestry.ubc.ca/prognosis/extension.html accessed November, 2008. Mäkelä, A. 1997 A carbon balance model of growth and self-pruning in trees based on structural relationships. For. Sci. 43: 7-24. Meidinger, D. and Pojar. J. 1991. Ecosystems of British Columbia. B.C. Min. For., Res. Br., Victoria, B.C. Spec. Rep. Ser. No. 6. Meng, S.X., Lieffers, V.J., and Huang, S.M. 2007. Modeling crown volume of lodgepole pine based upon the uniform stress theory. For. Ecol. & Manage. 251: 174-181. Monserud, R.A., and Sterba, H. 1996. A basal area increment model for individual trees growing in even- and uneven-aged forest stands in Austria. For. Ecol. & Manage. 80: 57-80. Oliver, C.D., and Larson, B.C. 1996. Forest Stand Dynamics. John Wiley and Sons, Inc. ISBN 0- 471-13833-9, New York, 520 p. Pacala, S. W., Canham, C. D., Silander Jr., J. A. 1993. Forest models defined by field measurements. 1. The design of a northeastern forest simulator. Can. J. For. Res. 23: 1980-1988. Pacala, S.W., Canham, C.D., Saponara, J., Silander, J.A., Kobe, R.K., and Ribbens, E. 1996. Forest models defined by field measurements: Estimation, error analysis and dynamics. Ecol. Monogr. 66: 1-43. Pinno, B.D., Lieffers, V.J., and Stadt, K.J. 2001. Measuring and modelling the crown and light transmission characteristics of juvenile aspen. Can. J. For. Res. 31: 1930-1939. Rouvinen, S., and Kuuluvainen, T. 1997. Structure and asymmetry of tree crowns in relation to local competition in a natural mature Scots pine forest. Can. J. For. Res. 27: 890-902. Rudnicki, M., Silins, U., and Lieffers, V.J. 2004. Crown cover is correlated with relative density, tree slenderness, and tree height in lodgepole pine. For. Sci. 50: 356-363. 41 Quan, N. T. 1988. The prediction sum of squares as a general measure for regression diagnostics. J. Busin. Econ. Stat. 6: 501-504 SAS Institute Inc. 2003. SAS/ETS user’s guide, version 9.1.3. Cary, NC: SAS Institute Inc. Schabenberger, O. and Pierce, F. J. 2002. Contemporary statistical models for the plant and soil sciences. CRC Press, Boca Raton, Florida. Schwinning, S. and Weiner, J. 1998. Mechanisms determining the degree of size asymmetry in competition among plants. Oecol. 113: 447-455. Soares, P. and Tome, M. 2001. A tree crown ratio prediction equation for eucalypt plantations. Ann. For. Sci. 58: 193-202. Sprugel, D. G. and Hinckley, T. M. 1988. The branch autonomy theory in W. E. Winner and L. B Phelds (eds.), Response of Trees to Air Pollution: the Role of Branch Studies, Proceedings of an Environmental Protection Agency Workshop held November 5-6 1987, in Boulder, Colorado, pp. 1-19. Stockdale, C., Taylor, S. and Hawkes, B. 2004. Incorporating mountain pine beetle impacts on stand dynamics in stand and landscape models: a problem analysis. In Mountain pine beetle symposium: challenges and solutions. October 30-31, 2003, Kelowna, British Columbia. T.L. Shore, J.E. Brooks, and J.E. Stone (editors). Nat. Res. Can., Can. For. Serv., Pac. For. Centre, Inform. Rep. BC-X-399, Victoria, BC. pp. 200-209. Studenmund, A. H. 2006. Using econometrics: a practical guide. Pearson Addison Wesley, New York, USA. Valentine, H. T. and Mäkelä, A. 2005. Bridging process-based and empirical approaches to modeling tree growth. Tree. Physiol. 25: 769-779. Valentine, H. T., Ludlow. A. R. and Furnival, G. M. 1994. Modeling crown rise in even-aged stands of Sitka spruce or loblolly pine. For. Ecol. & Manage. 69: 189-197. Vanclay, J. K., Skovsgaard, J. P. and Hansen, C. P. 1995. Assessing the quality of permanent sample plot databases for growth modeling in forest plantations. For. Ecol. & Manage. 71: 177-186. 42 West, P. W. 1995. Application of regression analysis to inventory data with measurements on successive occasions. For. Ecol. & Manage. 71: 227-234. Williams, H., Messier, C. and Kneeshaw, D. D. 1999. Effects of light availability and sapling size on the growth and crown morphology of understory Douglas-fir and lodgepole pine. Can. J. For. Res. 29: 222-231. Wood, C. S. and Unger, L. S. 1996. Mountain pine beetle: a history of outbreaks in pine forests in British Columbia, 1910 to 1995. Nat. Res. Can., Can. For. Ser., Pac. For. Centre, Victoria, B.C. Wu, D. M. 1973. Alternative tests of independence between stochastic regressors and disturbances. Econo. 41: 733-750. Zellner, A. and Theil, H. 1962. Three-stage least squares: Simultaneous estimation of simultaneous equations. Econo. 30: 54-78. 43 3. A hybrid model to estimate natural recruitment and growth in stands following mountain pine beetle disturbance2 3.1 Introduction The historical presence of Mountain Pine Beetle (Dendroctonus ponderosae Hopkins, MPB) in lodgepole pine (Pinus contorta Dougl. ex Loud. var. latifolia Engelm) forests in the interior of British Columbia (BC) is well documented (Heath and Alfaro 1990; Wood and Unger 1996; Hawkes et al. 2004). Tools to forecast regeneration and growth following MPB-disturbance are lacking, however (Zumrawi et al. 2005). Developing forest growth models suitable for MPB- disturbed forest is currently a priority in BC as an ongoing outbreak of MPB has reached epidemic levels (Griesbauer and Green 2006). Many stands affected by the current outbreak will go unsalvaged and be left to regenerate naturally. It is anticipated that a large proportion of BC’s mid-term timber supply will come from these unsalvaged stands (Patriquin et al. 2008). Thus, there is a particular need for forest growth models capable of modelling natural regeneration in stands affected by MPB. Traditionally, forest managers have relied on empirical models to obtain estimates of timber growth and yield (Mäkelä 2003). In empirical models, growth is estimated through the use of mathematical equations that predict the behaviour of the response variable. A criticism is that there is often no attempt to explain the underlying biological factors influencing growth (Vanclay 1994; Valentine and Mäkelä 2005). Despite this, empirical models can provide a high level of precision if they are well formulated and properly calibrated. Generally, this requires the calibration dataset be collected from permanent sample plots (PSPs) located within a particular population and include records of annual growth over several growth periods (Robinson and Ek 2003; Robinson and Monserud 2003; Valentine and Mäkelä 2005). The robustness of empirical models is questioned when applied to populations outside of the region to which they have been calibrated (Robinson and Monserud 2003). Likewise, MPB outbreaks can change forest conditions to the point that highly empirical models calibrated using pre-disturbance data are no longer suited to post-MPB conditions (Zumrawi et al. 2005). One empirically based growth model that is currently being used to forecast growth in the 2 A version of this chapter will be submitted for publication. Sattler, D.F., LeMay, V.M., Coates, D.K., and Marshall, P.L. Estimation of natural recruitment following mountain pine beetle infestation using SORTIE-ND and PrognosisBC. 44 forests of central and southeastern BC is PrognosisBC (Zumrawi et al. 2003; Hassani et al. 2004). Adapted from the Northern Idaho variant of the Forest Service Forest Vegetation Simulator (FVS) (Stage 1973), PrognosisBC is best suited for projecting existing stands composed of multiple species and age classes. This model can simulate a wide range of silvicultural treatments and is therefore, a valuable tool for forest practitioners interested in testing different harvesting regimes. To simulate regeneration, PrognosisBC uses an empirical imputation method that imputes regeneration from overstory characteristics using the principal of most similar neighbour (MSN) (Hassani et al. 2004). This approach uses a database composed of reference plots containing regeneration information, and overstory tree and site information to match to target plots which have only overstory and site information. Once the best match is found based on overstory and site information, regeneration from the reference plot is imputed to the target plot. For this technique to be successful, a reference database sampled from a wide range of stands is required. Alternatively, regeneration can be added manually to PrognosisBC. This approach is generally used when a specific planting density is prescribed, for example, following harvest. However, neither of these methods of simulating regeneration is ideal for use in MPB-disturbed stands. First, the reference database is limited for MPB-disturbed stands. Second, manually adding regeneration to unsalvaged MPB-disturbed stands is neither feasible nor accurate since the number of stands that will go unsalvaged is extremely large and the expected density of natural regeneration unknown. To extend the capabilities of PrognosisBC into unsalvaged MPB-disturbed stands, an alternative method of estimating natural regeneration is needed. In this study, the regeneration submodel from the forest growth and dynamics model SORTIE-ND (Astrup 2006) is explored as a means of providing PrognosisBC with estimates of natural regeneration. The SORTIE-ND model is capable of simulating many of the processes believed to strongly influence the establishment and survival of natural regeneration in MPB- disturbed stands, including light transmission to the understory, seed dispersal, and the availability of substrates amicable to seedling recruitment (Coates and Hall 2005; Dordel et al. 2008). The light behaviour in SORTIE-ND can predict solar radiation at any point in a stand (Canham et al. 1999). Estimates of solar radiation are used by the juvenile growth behaviour to predict seedling growth. Species-specific preferences for substrate types and canopy cover are used to predict seedling establishment. SORTIE-ND has no default parameters; thus, field data must be collected to parameterize the functions used by each submodel. The model was previously calibrated and tested for use in the Sub-Boreal Spruce (SBS) zone of BC by Astrup et al. (2007). 45 This paper describes the methods used to link SORTIE-ND to PrognosisBC. The linked model is termed a hybrid model since it brings together two different approaches to modelling, an empirical approach used by PrognosisBC, and a more process-oriented approach used by SORTIE- ND. The linkage allowed for estimates of natural regeneration to be passed from SORTIE-ND to PrognosisBC after an elapsed forecast time following MPB-disturbance. Following a description of the linkage, the hybrid model is evaluated in terms of: 1) the feasibility of the linkages between SORTIE-ND and PrognosisBC; and 2) the accuracy of model predictions over a mid-term planning horizon of 25 years. At the outset of testing, it was unclear at what point in the simulation process regeneration estimated in SORTIE-ND should be imputed into PrognosisBC. In theory, the best time to transfer regeneration would coincide with a period during post- disturbance stand development when the majority of new recruits have established and the increased competition for growing resources limits further seedling establishment. Thus, the second objective was to identify a point in the simulation process where the transfer of regeneration from SORTIE-ND to PrognosisBC provided estimates that were accurate over a 25- year planning horizon and which coincided with biological theories related to post-disturbance ingrowth. Diagnostics for evaluating these two objectives were based on comparisons of model outputs with mensurational data recorded in the sampled stands, as well as a comparison of model predictions against more general ecological expectations. The intent is for the hybrid model approach to have practical applications at the stand level. Thus, accuracy was assessed with comparisons to stems ha-1, basal area (m2 ha-1) and mean tree height (m) for different tree size classes. Model outputs for periodic growth (cm year-1) and percent of full sunlight (Global Light Index - GLI) were used to evaluate the hybrid model approach in terms of its ability to meet ecological expectations. 3.2 Methods 3.2.1 Study area The study area included stands located within an approximate radius of 250 km around Williams Lake, BC (52°08’18.59”N 122°08’31.07”W). Many of the sampled stands were located on the Chilcotin Plateau, situated within the former Cariboo Forest Region (Figure 3.1). 46 Figure 3.1 Location of stands sampled for the Williams Lake dataset (LeMay et al. 2007) This area has a long history of disturbance events caused by MPB (Stockdale et al. 2004; Aukema et al. 2006), although evidence suggests that stand replacing and ground fires have played a major role in influencing stand structure across much of the forest region (Hawkes et al. 2004). Historical records gathered from aerial surveys (Wood and Unger 1996), and tree ring analyses (Hawkes et al. 2004) indicate that the most recent significant disturbance caused by MPB, aside from the current outbreak, occurred during a period beginning in the late 1970s stretching until 1984/85. The reemergence of large populations of MPB in the early 1990s and subsequent disturbance brought about by the beetle is documented for some areas of the Chilcotin Plateau (Campbell et al. 2007). Studies within the Chilcotin plateau indicate that the most recent disturbance from large scale wildfires occurred around 1910/20 (Alfaro et al. 2004). Hawkes et al. (2004) noted that many of the stands in the Chilcotin Plateau that were selected for their study of post-MPB stand dynamics showed the effects of multiple disturbances caused by the late 1970s MPB outbreak and low-level ground fires in subsequent years. This combination of events has resulted in the creation of uneven-aged, mixed-species stands throughout much of the Chilcotin plateau (Heath and Alfaro 1990). The topography of the study area is characterized by gentle slopes and an elevation between 47 900 and 1500 m. Lodgepole pine is either present or is the dominant species in most of the upland stands in this area. The sampled stands were located in three major biogeoclimatic ecosystem classification subzones (Meidinger and Pojar 1991): Sub-Boreal Spruce (SBS), Sub-Boreal Pine Spruce (SBPS), and Interior Douglas-fir (IDF). 3.2.2 Data description and preparation Data for the parameterization of SORTIE-ND and testing of the hybrid model approach in this study were collected by LeMay et al. (2007). These data were gathered as part of a BC Forest Science Program project to obtain additional data for imputation of regeneration following MPB attack using PrognosisBC. The sampling design targeted stands where physical evidence and/or historical records indicated that MPB disturbance occurred roughly 25 years ago. A second group of stands were sampled where MPB disturbance had occurred roughly 8 years ago; however, the narrow set of site conditions and short time span between disturbance and sampling did not lend these latter stands to the type of testing carried out in this study and consequently, they were not included in the simulations. However, this second dataset did prove to be useful for the parameterization of SORTIE-ND substrate behaviours, as later described in this paper. Evidence of MPB disturbance was sought through a combination of historical forest inventory records, personal communication with local forest practitioners, and the presence of physical characteristics typical of MPB disturbance at the targeted site. For each selected stand, two to six plot centres, depending on the size of the stand, were established using systematic sampling with a random start. The distance between each plot centre ranged between 50 m and 150 m, with a minimum distance of 50 m from a road or any other significant opening. At each sample point, a nested fixed area plot was used. First, a large fixed-area plot of 11.28 m radius was established and measurements of diameter outside bark at breast height (DBH; 1.3 m above ground), total tree height, and height to the base of the live crown were collected for live trees greater than 7.5 cm DBH. Measures of crown diameter at crown base were recorded from two randomly selected trees of every species tallied in the plot. A 5.64 m radius plot was then established at the same point as the 11.28 m plots, where all trees greater than 2 cm and less than 7.5 cm DBH were measured for DBH and total tree height. On the perimeter of the 11.28 m radius plot, four 2.07 m radius plots were established on each of the four cardinal directions, while a fifth plot was located at the centre of the 11.28 m plot. Within these five small plots, seedlings less than 2 cm DBH and greater than 15 cm height were tallied into four height classes; 15 cm - 49 cm, 50 cm - 100 cm, 100 cm - 150 cm, and greater than 150 cm (but under 2 cm 48 DBH). The proportion of cover by different substrates in each 2.07 m radius plot and the type of substrate used by seedlings was also recorded. 3.2.2.1 Reconstruction of tree-lists Simple tree reconstruction methods, applied to each plot in the 2006 field data, were used to estimate stand structure characteristics (DBH, Height, basal area, stems ha-1) one or two years following MPB-disturbance. These estimated characteristics were used as the starting point for model simulations. Based on sources from the literature (Heath and Alfaro, 1990; Hawkes et al. 2004), the MPB population in the study area last peaked in 1981 (i.e., 25 years ago). The reconstruction of stands was completed in two phases; the first phase involved back-casting live trees, snags and downed trees recorded in the 11.28 m radius plot, while the second phase of back-casting was exclusive to those trees recorded in the 5.64 m radius plot. A terse description of reconstruction methodologies for live trees, snags and downed trees recorded in 11.28 m radius plots is provided here, while full details of the stand reconstruction process are provided in LeMay et al. (2007). Briefly, wildlife tree classification codes recoded for snags and downed trees, based on the level of decay, were used to indicate which trees were alive 25 years ago. From the list of live trees recoded in 2006, tree DBH and height values were obtained by reducing diameters backwards in 10 year cycles with heights re-estimated at each time. Some live trees in the 2006 tree-list were back-casted DBH < 0 and were therefore not present in the 1981 tree list. A modified version of the adult tree growth model in PrognosisBC was used to estimate the 10 year diameter growth increments. Back-casting of trees recorded in the 5.64m radius plot employed the same general methods, except that a modified version of the PrognosisBC small tree 5-year height increment model was used to back-cast tree heights, and diameters were estimated from tree height at each cycle. Parameters for the small tree 5-year height increment model were obtained from Robinson (2004). 3.2.3 The SORTIE-ND Model The SORTIE-ND model contains several behaviours, each corresponding to a biological process. They range from light behaviours to seedling establishment and substrate behaviours. In this paper, the description of the SORTIE-ND model is restricted to the behaviours directly related to growth, mortality, and natural regeneration, as these were the most important processes within the context of this study. Following this is a description the parameterization of SORTIE- ND. 49 3.2.3.1 SORTIE-ND growth and mortality Several behaviours within SORTIE-ND are specific to tree size classes, or life-history stages as they are referred to in the model. For this study, these were: Seedlings: trees ranging in height from 0.1 cm (minimum height allowed by SORTIE) to 1.35 m; Saplings: trees greater than 1.35 m height and between 0.01 and 7.5 cm DBH; Adults: trees greater than or equal to 7.5 cm DBH; and Snags: adult trees that died during the simulation. In the version of SORTIE-ND used for this study, diameter increments for trees ≤7.5 cm DBH were estimated using the following functions: ε+ + = ×− )(1 10 GLIcbe aGrowthD (3.1) ε+ + ×+ = ×− )(1 1010 GLIdce DbaGrowthD (3.2) ε+××+= cGLIDbaGrowthD )100/()10(10 (3.3) where D10Growth is the estimated diameter increment 10 cm from the base of the tree. Estimates for aspen (Populus tremuloides Michx.) and hybrid spruce (Picea engelmannii x glauca (Moench) Voss) were derived from Eq. 3.1 and Eq. 3.2, respectively. Lodgepole pine and Douglas-fir (Pseudotsuga menziesii var. glauca (Beissn.) Franco) estimates were derived from Eq. 3.3. GLI is the Gap Light Index for the tree in consideration, D10 is the diameter at 10 cm and a, b, c and d are species-specific parameters. Incremented diameters at 10cm were used in an allometric function to estimate diameter at breast height (DBH, 1.3 m). Adult tree diameter increments were calculated from: fectCrowdingEfectShadingEffSizeEffectMaxGrowthDBHGrowth ×××= (3.4) where DBHGrowth is estimated diameter increment at breast height, and MaxGrowth, SizeEffect, ShadingEffect and CrowdingEffect are species-specific parameters between 0 and 1 that are used 50 to reduce the maximum growth rate of a tree. The CrowdingEffect parameter is based on a neighbourhood competition index (NCI) that measures the potential reduction in growth based on the proximity of other trees with a minimum of height of 2 metres. In SORTIE-ND, tree vigour, expressed through annual radial growth, is used to predict the probability of mortality for juvenile trees. An additional density dependent mortality function is applied to juvenile aspen. Adult tree mortality due to competition is estimated as a function of the relative increment of a tree and is directly related to adult DBH growth and the MaxGrowth and SizeEffect parameters in Eq. 3.3. Lastly, a stochastic background mortality function is applied to all juvenile and adult trees. 3.2.3.2 SORTIE-ND natural regeneration Natural regeneration in SORTIE-ND is simulated through a regeneration behaviour. First, the simulation of seedling dispersal is accomplished through spatial or non-spatial functions that use species-specific parameters. Spatial and non-spatial functions can operate simultaneously. Next, establishment, which constitutes the progression from seeds to seedlings, is estimated as a function of substrate availability, substrate preference, and proximity of conspecifics. 3.2.4 Parameterization of SORTIE-ND Behaviours Not all parameters for the SORTIE-ND behaviours could be estimated from the field data collected within the Williams Lake study area. The primary source for parameters that could not be estimated from the field data was Astrup et al. (2007). This section outlines those parameters that were obtained from the Williams Lake field data and then describes the refinements to the parameters supplied by Astrup et al. (2007). Using the field data from Williams Lake, sapling and adult tree height was estimated through: ( ) ε+−×−+= ×− DBHbeMaxHeightHeight 1)35.1(35.1 (3.5) where Height is the estimated total height of the tree, MaxHeight is the maximum height attainable by a given species on the sampled sites, and b is a species-specific parameter. The constant 1.35 is based on measurement of DBH at 1.35 m above ground used in SORTIE-ND; this was not altered to the 1.3 m used in the sample data. The MaxHeight parameter in Eq. 3.5 represents the upper asymptote for tree height, and limits the height of individual trees. In theory, this value represents the species’ maximum potential height given the environmental conditions. 51 However, the Williams Lake datasets did not contain a sufficient number of trees that could be considered to be representative of this maximum height. Instead, estimates of MaxHeight for each species were obtained following inspection of the model dataset and peer-reviewed literature, including Nigh (2002) and Temesgen and Gadow (2004). Eq 3.5 was then fit for each species using nonlinear ordinary least squares with a Gaussian iterative search method (PROC NLIN, SAS Institute Inc. 2003). The fitted equations showed no lack of fit and had pseudo-coefficient of determination (R2) values of 0.83, 0.86, 0.80 and 0.84 for aspen, Douglas-fir, lodgepole pine and hybrid spruce, respectively. A nonlinear three-stage least squares (N3SLS) fitting technique was used to estimate crown depth and crown radius (Chapter 2 of this thesis). In SORTIE-ND, crowns are represented as cylinders. A pseudo-measure of crown radius, defined as one quarter of the diameter at the base of the crown, was used in SORTIE-ND to adjust the representation of crowns to the more conical or paraboloid crown shape observed for many conifers. Canham et al. (1999) described the adjustment to one quarter crown diameter as a method of obtaining the “effective” crown radius. The adjustment to the effective crown radius should result in an increased amount of light penetrating through the upper canopy. Substrate parameters were obtained from the Williams Lake dataset from the plots established in stands attacked by MPB 8 years ago. Substrate data collected from these plots were preferred as they were more likely to be reflective of conditions immediately following MPB disturbance. The parameters were estimated based on the proportion of ground occupied by a given substrate type, where values summed to 100%. The substrate types included: 1) a forest floor litter and moss pool; 2) fresh logs, defined as newly fallen trees (fresh logs are then converted into decayed logs by the substrate behaviour); 3) decayed logs; and, finally, 4) exposed mineral soil. Obtaining substrate proportions reflective of conditions following MPB attack were critical to the simulation process. In SORTIE-ND, each substrate has a favourability rating that is specific to each species. These ratings are one of two factors that determine the number of seedlings that successfully establish. The other factor is the amount of canopy cover estimated to be above the various substrates. All other parameters were obtained from Astrup et al. (2007). This included parameters relating to the transmittance of light through crowns, tree growth, mortality, snag fall-down, seed dispersal, and seedling establishment behaviours. To estimate these parameters, Astrup et al. (2007) sampled mainly from stands located on mesic sites in the SBS Moist Cold Subzone, 52 Babine variant (SBSmc2) variant situated around the town of Smithers, located in the northwestern region of the province of BC. Relative to the Williams Lake study site, growing conditions around Smithers are characterized by higher rainfall, and higher mean summer temperatures (Environment Canada 2007). These climatic differences result in slightly higher productivity ratings for spruce, aspen and lodgepole pine (BC Ministry of Forests and Range 2006). Climatic and site conditions are not explicitly modeled in SORTIE-ND; rather, environmental influences are reflected in the parameter estimates. By using parameters from the Smithers area, some minor overestimation of tree growth was expected. Arguably, the most notable difference in growing conditions between the two study sites is the amount forest floor moss cover. The drier conditions of the Chilcotin Plateau tend limit the amount of forest floor moss cover exposing a greater proportion of substrates favourable to the germination of serotinous and non-serotinous pine seeds (BC Forest Practices Board 2007). However, as noted, substrate parameters were estimated using the Williams Lake dataset. Although it was assumed that most of the parameters obtained from Astrup et al. (2007) would provide reasonable estimates, extra attention was paid to the seed dispersal parameters for lodgepole pine. In SORTIE-ND, seeds can be dispersed through a spatial or non-spatial function, or through a combination of the two. Parameters were available for both of these functions, although the non-spatial parameter was intended to only simulate a background level of seed dispersal. This posed a challenge to the hybrid model simulations since the location of trees imported into SORTIE-ND were randomized and it was unknown how this would effect the behaviour of the spatially dependent seed dispersal function. Having estimated the substrate parameters from the Williams Lake dataset, a preliminary round of simplified simulations was used to evaluate the potential effects. Results were compared against simulations that used an alternative non-spatial seed dispersal parameter while the spatial parameters were set to zero. The alternative parameter was derived from records collected within the SBS Dry Warm Subzone, Horsefly variant (SBSdw1) subzone located in the Cariboo Forest Region (BC Ministry of Forests 1999). These records indicated that between 24,000 and 70,000 lodgepole pine seeds ha-1 per year, were counted over a 3-year period. A conservative estimate of 3 seeds m-2 yr-1 was used for the non-spatial seed dispersal parameter. Results seemed to corroborate the use of the alternative non-spatial parameter, as predicted seedling densities fell between the upper and lower boundaries of ecological expectations for the conditions of the study area. Conversely, lodgepole pine seed dispersal parameters from Astrup et al. (2007) seemed ill-suited for the stands sampled in the Williams Lake study site, since predictions of lodgepole pine seedling densities fell far 53 below recognized natural stocking levels for the study area (Alfaro et al. 2004; BC Forest Practices Board 2007; Hawkes et al. 2004). For spruce and aspen, alternative seed dispersal parameters amicable to the environmental conditions of the study site could not be found. Thus, the parameters from Astrup et al. (2007) were used. 3.2.4 The PrognosisBC Model The most important PrognosisBC sub-models related to the hybrid modelling approach are the growth and mortality sub-models. These are described in this section, followed by a description of the parameters used to initialize the model for this study. 3.2.4.1 PrognosisBC Growth and Mortality In PrognosisBC, trees ≤ 2 in (5.08 cm) DBH are incremented first through a height increment model (Robinson, 2004): feGTH =ˆ (3.6) where f is defined as: 1 2 3 4 5 6 7 8 cos( ) sin( ) log( ) 100 f a a ASP SL a ASP SL a SL BAL a H a H a CCF a = + + + + + + + (3.7) and where ASP is the stand aspect (in radians), BAL is the basal area of trees larger in DBH than the subject tree in m2 ha-1, CCF is the crown competition factor, H is tree height in m, SL is the slope ratio for the stand and GTH ˆ is the estimated 5-year height increment in m. A height-based allometric function is then used to estimate the diameters of small trees: ( ) 65.00.2 )75.101232.0( 36 )5.4(ˆ 21 +−× −××+−= RELH RELHCCFAVHHbHBD b (3.8) where AVH is the top height for the stand in feet, and RELH is the relative height calculated as (H - 4.5) / (AVH - 4.5) (0.0 ≤ RELH ≤ 1.0). Increments of DBH are obtained through subtraction. For PrognosisBC, estimates from Eq’s 3.7 and 3.8 are converted from imperial to metric units. When diameters of small trees reach the upper threshold of 4 in (10.2 cm), heights estimated 54 through Eq. 3.6 are combined with diameter-weighted estimates of large tree height growth. This allows for a smooth transition from small to large tree growth models. The diameter growth of large conifer trees is estimated through: feSDD =ˆ (3.9) where f is defined as: 2 1 2 3 4 5 6 2 7 8 9 10 11 12 cos( ) sin( ) 100 100 100 log ( 1) EL EL CCFf a a a a ASP SL a ASP SL a BAL BAL a DBH a DBH a BAL a a a CR DBH DBH = + + + + + + + + + + + + (3.10) and CR is the ratio of crown depth to total tree height, DDS is square of the 10-year DBH growth increment in cm2, EL is the stand elevation in m and 1a to 12a are species-specific parameters. Mortality in juvenile trees is a function of tree size and vigour. The probability of mortality increases for smaller trees and for trees with less vigour. The mortality of larger trees is primarily a function of the maximum occupancy of a site, an upper threshold that is estimated from the maximum basal area for a given stand. As stands approach this threshold, growth rates decline and the probability of mortality increases. 3.4.2 PrognosisBC Parameters Parameters for the growth equations and all other core functions are already built into the model. In the version used for this study, parameters were available for the Nelson, Kamloops and Cariboo forest districts. Within the Cariboo forest district, calibration of the model was limited to the IDF dk3 01 (i.e., dry cool) variant. Although three of the stands tested were from this same variant, the remaining stands were from different variants within the IDF, SBS and SBPS zones. Nevertheless, parameters for the IDF dk3 01 variant were used in all simulations. The rationale behind this choice was that the stands being tested generally developed under mesic and medium soil moisture and nutrient regimes. Also, it was assumed that climatic and site differences would produce only minor differences in growth predictions. The Fire and Fuels Extension (FFE) (Reinhardt and Crookston, 2003) to PrognosisBC was used to control the rate of snag-fall. Since the snag-fall rate among pine trees killed by MPB is known to be lower than the average rate of fall for the species (Hawkes et al. 2005), the SNAGFALL keyword in the FFE was used to simulate conditions where 95% of snags 25cm DBH or larger would fall within a 55 span of 30 to 90 years. 3.2.5 Hybrid Model Structure A schematic of the linkages and overall flow of the hybrid model approach used for this study is shown in Figure 3.2. Figure 3.2 Schematic of the hybrid model approach. Overstory (O/S) and understory (U/S) tree-lists of the same population are imported into SORTIE-ND and PrognosisBC at Time 1, which is shortly after disturbance. At Time 2, PrognosisBC simulations are paused while SORTIE-ND is stopped. Understory trees are exported from SORTIE-ND, formatted and imported into PrognosisBC, replacing trees within the same size class. PrognosisBC is restarted with an updated tree-list and projected to points in time within the mid-term planning horizon (Time 3). The modelling process began by importing the same combined overstory and understory tree tree-list into SORTIE-ND and PrognosisBC shortly after the time of disturbance, noted henceforth as Time 1. The tree-lists were a list of the species, sizes (DBH and height), and stems per ha (i.e., expansion factor) for each tree in the list. For SORTIE-ND, the tree-lists were modified to a file of single trees, with an expansion factor for each tree of one. For PrognosisBC, the input to the model was the tree-lists with the expansion factors. The practical implications of this distinction with regards to linking the two models are discussed in a latter section. Once the models were populated with individual tree records, parameters for the tree population were specified and the simulations began. Time 1 Time 2 Time 3 56 As a simulation proceeded from Time 1 to Time 2, the PrognosisBC model was operated independently from the SORTIE-ND model. During this phase of the simulation, no new seedlings were added to the tree-list projected in PrognosisBC. At Time 2, SORTIE-ND simulations were stopped, and simulations in PrognosisBC were paused. Using predetermined upper and lower size classes, seedlings and advanced regeneration were selected from the SORTIE-ND projections, reformatted and imported into the paused simulation in PrognosisBC. Upon import, the cohort of understory trees selected from SORTIE-ND replaced all trees within the same size class that remained in tree-list forecasted to Time 2 using PrognosisBC. In theory, Time 2 represents the point at which new seedlings and advanced regeneration have occupied the growing space in the understory that became available either through disturbance (e.g., the MPB outbreak), crown rise, or both (Oliver and Larson 1996). Although the ingress of seedlings may continue after this point, competition for resources and growing space limit the number of successful germinants. In practice, the time at which this occurs may be difficult to define as there are several factors that influence the timing, including changing understory light conditions, soil moisture conditions, proximity of seed sources, prevailing wind direction, and species composition (Oliver and Larson 1996; Bainbridge and Strong 2005). Identifying a point in time that provides a reasonable estimate of the onset of stabilized conditions following understory reinitiation is a key determinant in the ability of this hybrid model approach to provide accurate estimates for the mid-term planning horizon. In the current state of this hybrid model, the linkages that have been described are static in the sense that seedling and advanced regeneration populations in PrognosisBC simulations are not dynamically updated through a feedback mechanism with SORTIE-ND. There is only the transfer of tree-list information from SORTIE-ND to PrognosisBC. Conceptually, however, a dynamic feedback system could be achieved by creating an two-way flow of information between the models, where PrognosisBC would simulate the growth of the stand over one cycle (5 or 10 years) following an update of its seedling and advanced regeneration population, then provide SORTIE- ND with a full tree-lists following this growth cycle, which SORTIE-ND then would use to project the next cohort of seedlings and advanced regeneration targeted for transfer to PrognosisBC. 3.2.6 Simulations Thirteen stands were randomly selected from the Williams Lake dataset and projected using the hybrid model. Three of the stands were from the SBPS, two from the SBS and eight from the 57 IDF. The selected stands were considered to be structurally complex, composed of multiple species and multiple age classes. Lodgepole pine was the dominant overstory species in eight of the stands, and the dominant understory species in nine stands. Aspen, Douglas-fir and hybrid spruce occurred mainly as secondary species in the understory, although for one stand Douglas-fir was the only understory species recorded, while for another stand aspen was the dominant species. Douglas-fir was the dominant overstory species in four stands, but was completely absent from the overstory in seven of the stands tested. For each of the 13 stands, the reconstructed plot data for 1981 were used as inputs to the models (Time 1 data). Field measurements from 2006 were used to obtain a stand-level tree-list for each of the 13 stands to assess model projections, as discussed later under model evaluation. Summary statistics for the 13 stands are given in Tables 3.1 and 3.2. Assuming that dead pine trees in the stands backcast to 1981 died as a result of attack from MPB, the proportion of stand basal area (m2 ha-1) that was killed by MPB for these 13 stands ranged from a high of 74% to a low of 19% with an average of 43%. However, these snags were not included in the tree-lists used to initialize SORTIE-ND and PrognosisBC. It was felt that excluding these snags from the starting tree-lists would not significantly alter predictions of understory light conditions since many of the pine trees that were likely killed by MPB in the 3 or 4 years prior to the date which the stands were backcast. Therefore, it is likely that many of these trees would have dropped most of their needles, thus, casting less shade. For each stand, the simulated growing environment in SORTIE-ND was defined by a square 9 ha plot (300 m by 300 m). To create a 9 ha stand from each tree-list, the tree-lists for all plots of a stand were pooled. Then, the trees were replicated to obtain a 1 ha area. For example, if there were five 11.28 m radius (0.04 ha) plots within a stand, then the combination of the five plots represented 0.20 ha. Each tree was then replicated four more times to obtain 1 ha. The 1 ha tree- list was then replicated eight more times to give the equivalent of a tree-list that covered 9 ha. This approach was used to expand the list of trees recorded in the 5.64 m radius with the appropriate adjustments to reach a stand of 9 ha. No trees remained in the 2.07 m radius plots following reconstruction. Since the spatial coordinates of trees were not collected during field sampling, individual trees were assigned a random location coordinate prior to importing tree lists into SORTIE-ND. 58 Table 3.1 Mean stem density (stems ha-1) and basal area (m2 ha-1) by tree size class in 1981 (estimated through reconstruction) for the 13 stands used in the hybrid modelling simulations. < 2 cm DBH 2 ≤ DBH < 7.5 cm 2 ≤ DBH < 7.5 cm Stand# No. plots Species Stems ha-1 Stems ha-1 Basal area Stems ha-1 Basal area 1 4 Aspen 0 0 0 0 0 Douglas-fir 200 55 0.2 545 17.4 Lodgepole pine 0 0 0 114 5.6 2 2 Lodgepole pine 450 170 0.2 114 4.6 Hybrid Spruce 0 0 0 0 0 3 3 Aspen 419 74 0.1 62 1.0 Lodgepole pine 197 80 0.2 456 15.4 Hybrid Spruce 42 9 0.1 8 1.0 4 2 Lodgepole pine 1199 150 0.4 328 10.9 Hybrid Spruce 0 0 0 13 1.0 5 3 Douglas-fir 461 43 0.2 479 31.5 Lodgepole pine 0 0 0 90 4.0 Hybrid Spruce 0 8 0.1 0 0 6 6 Aspen 17 0 0 5 0.1 Douglas-fir 5 0 0 0 0 Lodgepole pine 1629 372 0.8 549 7.5 Hybrid Spruce 0 4 0.1 0 0 7 3 Aspen 0 0 0 8 0.3 Douglas-fir 494 252 0.7 673 15.1 Lodgepole pine 0 18 0.1 149 4.9 8 2 Lodgepole pine 299 90 0.3 324 11.0 9 4 Lodgepole pine 249 40 0.1 298 11.32 10 3 Lodgepole pine 197 115 0.4 706 17.2 Hybrid Spruce 107 151 0.2 127 2.5 11 3 Aspen 329 0 0 17 0.7 Douglas-fir 33 33 0.2 586 17.8 Lodgepole pine 0 0 0 659 33.1 Hybrid Spruce 0 18 0.1 211 6.2 12 3 Aspen 923 144 0.3 71 1.1 Lodgepole pine 1006 75 0.1 152 5.2 Hybrid Spruce 0 0 0 0 0 13 4 Aspen 199 20 0.1 31 0.4 Douglas-fir 6 12 0.1 13 0.1 Lodgepole pine 424 264 0.5 750 13.0 Hybrid Spruce 25 12 0.1 6 0.1 59 Table 3.2 Mean stem density (stems ha-1) and basal area (m2 ha-1) by tree size class in 2006 for the 13 stands used in the hybrid modelling simulations. < 2 cm DBH 2 ≤ DBH < 7.5 cm 2 ≤ DBH < 7.5 cm Stand# No. plots Species Stems ha - 1 Stems ha-1 Basal area Stems ha-1 Basal area 1 4 Aspen 297 0 0 0 0 Douglas-fir 2042 200 0.5 581 23.0 Lodgepole pine 0 0 0 75 3.6 2 2 Aspen 222 0 0 0 0 Lodgepole pine 7577 1550 3.0 287 6.9 Hybrid Spruce 445 0 0 0 0 3 3 Aspen 4902 433 0.7 133 2.2 Lodgepole pine 5150 566 0.9 350 8.4 Hybrid Spruce 247 33 0.1 25 0.3 4 2 Lodgepole pine 5794 2450 3.6 438 13.4 Hybrid Spruce 520 0 0 12 0.4 5 3 Douglas-fir 8270 667 1.2 525 36.7 Lodgepole pine 49 0 0.1 75 3.8 Hybrid Spruce 0 33 0 8 0.1 6 6 Aspen 4060 50 0.1 4 0.1 Douglas-fir 50 0 0 4 0.1 Lodgepole pine 3045 1483 2.5 917 15.9 Hybrid Spruce 0 0 0 4 0.1 7 3 Aspen 198 0 0 8 0.4 Douglas-fir 3268 667 1.4 867 21.9 Lodgepole pine 50 0 0 133 5.8 8 2 Aspen 222 0 0 0 0 Lodgepole pine 6760 1700 2.29 275 8.5 Hybrid Spruce 74 50 0.1 0 0 9 4 Aspen 37 0 0 0 0 Lodgepole pine 6165 975 0.9 287 14.1 Hybrid Spruce 334 0 0 0 0 10 3 Aspen 99 0 0 0 0 Lodgepole pine 990 233 0.6 792 24.2 Hybrid Spruce 544 233 0.5 292 5.7 11 3 Aspen 544 466 0.3 8 0.4 Douglas-fir 247 0 0.1 483 22.4 Lodgepole pine 0 0 0 542 30.7 Hybrid Spruce 990 67 0.1 200 7.0 12 3 Aspen 2129 1133 1.3 200 3.1 Lodgepole pine 7280 2166 2.5 275 7.8 13 4 Aspen 2525 175 0.2 31 0.7 Douglas-fir 148 50 0.1 31 0.6 Lodgepole pine 3194 600 1.3 794 18.1 Hybrid Spruce 0 50 0.1 19 0.3 60 The transfer of tree-lists from SORTIE-ND to PrognosisBC followed the model flow depicted in Figure 3.2. To determine the most appropriate point at which to transfer trees less than 7.5 cm DBH from SORTIE-ND to PrognosisBC (i.e., Time 2), four transfer times were tested. These were 5, 10, 15 and 20 years following the start of the simulations. Since the time between the reconstructed plot data and the measures in 2006 spanned 25 years, Time 3 was defined as 25 years from Time 1. For SORTIE-ND, each growth cycle was defined as a single year, while in PrognosisBC, growth cycles operated on a 5- year basis. For the 5-year hand-off simulation, a list of trees < 7.5 cm DBH was obtained from SORTIE-ND at the end of the fifth growth cycle. Since this list represented the number of trees over a simulated 9 ha plot, an expansion value of 0.111 stems per hectare (i.e., from a 1 ha area) was given to each tree in this list to be compatible with the tree-list inputs of PrognosisBC. Trees selected from SORTIE-ND were inserted to the PrognosisBC tree-list following the initial 5-year growth cycle, and replaced any trees < 7.5 cm DBH. Included in the tree-lists imported from SORTIE-ND, were trees that were less than 1.3 m in height and therefore, had no DBH (in SORTIE-ND, diameters for seedlings are recorded at 10 cm above the ground). For the transfer to PrognosisBC, a DBH of 0.1 was assigned to all imported trees below 1.3 m in height. Making this adjustment did not have implications for the modeling process as PrognosisBC uses a height increment model for trees of this size, and thus the input DBH was not used for these very small trees. The newly updated tree-list was then projected forward for 20 years. The same approach was used for all remaining transfer times (i.e., variations in Time 2) and adjusting the amount of time from Time 2 to Time 3 in order to obtain a 25-year projection. Following the initialization of the simulations in SORTIE-ND for the 13 randomly selected stands, it quickly became apparent that the length of time required by SORTIE-ND to simulate the 25-year period was very long. The average run time for a stand to complete 25 time steps (1 year/time step) was 80 hours. Therefore, the simulated hand-off for the 5 and 10 year periods from Time 1 to Time 2 were completed for all stands. However, only 10 of the 13 stands were simulated for the 15 year hand off, and nine stands completed the 20 year period. For seven of the 13 stands, a projection period of 25 years was also simulated. 3.2.7 Evaluation To evaluate the performance of the hybrid model approach, mensurational data collected in 2006 from the 13 stands selected for testing were compared to predicted values obtained at the end of the 25 year simulation period for the four transfer points (5, 10, 15 and 20 year transfers). 61 Additionally, the 2006 field data were compared to results obtained from running the stands in SORTIE-ND for the entire 25 year period. Similar comparisons were also made to results obtained from using PrognosisBC alone; however, since it was known beforehand that natural regeneration would not be simulated in PrognosisBC, comparisons of predicted results to field data were only useful for the larger tree size classes. The accuracy of the model predictions were assessed through the average differences between observed and estimated values (mean bias) for basal area ha-1, stems ha-1 and tree height. Mean biases for basal area ha-1 were calculated for two tree size classes (2.0 ≤ DBH < 7.5cm; and DBH ≥ 7.5cm), while mean biases for stems ha-1 were calculated for three size classes (< 2.0cm DBH; 2.0 <= DBH < 7.5cm; and DBH > 7.5cm). Mean biases for tree heights were calculated for five size classes (2.0m ≤ height < 5m; 5m ≤ height < 8m; 8m ≤ height < 11m; 11m ≤ height < 14m; and 14m ≤ height < 24m). Since the number of stands simulated for each tree list transfer was unequal, calculations of the mean biases may have been affected. However, evaluation through mean biases seemed reasonable since the range of stand types simulated at each transfer time remained relatively unchanged. As well as mean bias, the root mean square error (RMSE) was calculated for each variable and size class. Further evaluation was performed by comparing model outputs to what Robinson and Ek (2003) refer to as emergent properties. This entailed examining the interactions among the hierarchical processes that are modeled, from individual trees to the stand. Specifically, emergent trends pertaining to periodic radial growth and understory light levels were compared to ecological expectations. This allowed for an evaluation of variables that were not collected in the field yet were critical to the overall performance of the models. Emergent trends also were compared to data from the literature which allowed for the models to be evaluated at points in time other than 25 years after MPB attack. For example, Heath and Alfaro (1990) reported on the growth response of MPB affected stands observed 14 years after disturbance. This allowed for an evaluation of the general trends rather than just the final model estimates. Finally, a comparison of predicted growth rates between SORTIE-ND and PrognosisBC was completed in order to shed some light on the source of the biases. Growth functions derived from PrognosisBC were used to reconstruct the stands and thus, represent what can be considered to be the expected growth rate. Predicted growth rates from SORTIE-ND that exceed or fall below this expected rate indicate where the hybrid model may be underperforming. 62 3.3 Results 3.3.1 Stand Simulation Mean Biases Since lodgepole pine was the dominant species in most of the simulated stands, reporting of results is mainly focused on this species. The smallest mean bias for the basal area of advanced pine regeneration (trees 2 ≤ DBH <7.5cm) was generated through the 15-year tree list transfer point simulations (Table 3.3), hereafter referred to as “the 15-year transfer”. For the simulations run using this transfer point, the predicted mean basal area of advanced pine regeneration was 1.47 m2 ha-1, while the actual mean basal area was 1.59 m2 ha-1. The mean predicted basal area represented a 45% increase in basal area from the beginning of the simulations. The variability in the mean biases for advanced pine regeneration was generally low; however, there were a few stands with large overestimates of basal area. For example, one stand was estimated to have 7.1 m2 ha-1 of advanced pine regeneration, when in fact the stand actually only had 2.98 m2 ha-1. Comparatively, the SORTIE-ND simulations resulted in larger mean biases for basal area of advanced pine regeneration than the 15-year transfer, although the variability in the mean biases were similar to those for the 15-year transfer. For estimates of the basal area of larger trees (trees ≥ 7.5cm DBH), the 15-year transfer also seemed to do a good job and greatly improved upon estimates provided by SORTIE-ND. The 15- year transfer provided an estimated mean basal area of 10.44 m 2 ha -1 for adult pine trees, which was an increase of 26% from the basal area at Time 1. Except for two of the 10 stands used in this transfer, the variability of the prediction biases for adult pine was low, remaining within 4 m 2 ha -1 of actual values. The largest bias for a single stand was an overestimate of 13 m 2 ha -1 , which was alarming since all the stands tested were of similar density and species composition. For comparison, PrognosisBC was run independently for each stand for the entire 25-year projection. Not surprisingly, these runs produced the largest underestimate for the basal area of pine ≥ 7.5 cm DBH, an average difference of 5.07 m2 ha-1 from actual values. It was slightly surprising to see that when SORTIE-ND was run independently, the basal area of pine ≥ 7.5 cm DBH was overestimated by upwards of 9 m2 ha-1. Almost certainly, this is due to the use of SORTIE-ND growth parameters obtained from stands that are generally considered to be more productive than the stands used in the simulations. 63 Table 3.3 Mean bias and root mean square error (in parentheses) of basal area and stem density classes for four transfer simulations and independently run SORTIE-ND simulations Basal Area Stems ha-1 2≤DBH<7.5 DBH≥7.5 DBH<2cm 2≤DBH<7.5 DBH≥7.5 5yr Transfer Aspen 0.1 (0.3) 0.7 (1.1) 1385 (2173) 171 (339) 39 (69) Douglas-fir 0.3 (0.6) 2.1 (3.6) 591 (2043) -494 (1841) 90 (119) Lodgepole pine -10.3 (11.6) 2.3 (4.2) 4187 (4997) -2219 (2859) -16 (188) Hybrid spruce 0.1 (0.1) 0.3 (0.9) 228 (376) 28 (54) 30 (53) 10yr Transfer Aspen 0.2 (0.4) 0.5 (0.8) 1270 (2081) 166 (352) 18 (34) Douglas-fir 0.3 (0.6) 1.1 (2.2) 835 (2431) -1076 (2437) 62 (94) Lodgepole pine -3.3 (3.9) 0.3 (5.0) 4187 (4997) -1475 (1927) -273 (785) Hybrid spruce 0.1 (0.1) 0.1 (0.4) 225 (374) 17 (60) 26 (55) 15yr Transfer Aspen 0.2 (0.4) 0.5 (1.0) 1324 (1983) 167 (403) 27 (59) Douglas-fir 0.2 (0.5) 0.2 (1.0) 1114 (2864) -696 (1926) 40 (53) Lodgepole pine 0.1 (1.7) -0.4 (4.5) 4729 (5365) 31 (1042) -307 (842) Hybrid spruce 0.1 (0.1) -0.1 (0.1) 141 (251) -9 (24) -3 (8) 20yr Transfer Aspen 0.2 (0.5) 0.5 (1.0) 1574 (2393) 192 (404) 30 (59) Douglas-fir 0.3 (0.5) 0.1 (0.7) 1053 (2724) -693 (1926) 30 (45) Lodgepole pine 1.4 (1.8) 0.7 (1.17) 3402 (4438) 1147 (1375) -69 (109) Hybrid spruce -0.1 (0.1) -0.1 (0.14) 110 (208) -1 (23) 1 (6) SORTIE-ND Aspen 0.3 (0.5) -1.5 (2.6) 1657 (2536) 223 (454) -26 (40) Douglas-fir 0.2 (0.3) -0.7 (1.5) 136 (483) 49 (100) 7 (12) Lodgepole pine 1.9 (2.1) -9.0 (10.2) 4411 (4985) 1483 (1616) 6 (94) Hybrid spruce -0.1 (0.1) -0.5 (0.7) 145 (236) -11 (22) -37 (41) Estimates of basal area for aspen, Douglas-fir and hybrid spruce advanced regeneration were also quite good for the 15-year transfer. However, since the actual basal area for these three species tended to be quite low for most stands, it was difficult to assess model performance for these species. Nevertheless, it was encouraging that no exceedingly large overestimates of basal area were predicted for these three secondary understory species. For aspen and Douglas-fir, the basal area of advanced regeneration was predicted to decrease from 0.15 m2 ha-1 and 0.13 m2 ha-1 64 to 0.03 m2 ha-1 and 0.01 m2 ha-1 , respectively. For hybrid spruce, changes from the beginning of the simulations were only marginal. For aspen, Douglas-fir and spruce the mean predicted adult basal area was 0.21 m2 ha-1, 15.00 m2 ha-1, and 0.19 m2 ha-1, respectively. For aspen and Douglas- fir, these estimates amounted to an increase from Time 1 of 0.16 m2 ha-1 and 5.5 m2 ha-1, respectively. Adult spruce basal area was predicted to decrease by 0.47 m2 ha-1 from Time 1. Estimates for stems ha-1 of new regeneration (trees < 2.0cm DBH) were largely underestimated for all species, regardless of tree-list transfer time. The best estimate for the stem density of new pine regeneration was obtained through the 20-year tree-list transfer simulations. In fact, only the 20-year tree-list transfers and the simulations run to completion in SORTIE-ND predicted the presence of trees in the < 2.0cm DBH size class after 25 years. The variability in the biases of new pine regeneration for the 20-year tree-list transfer simulations were quite high, with a maximum overestimate of 2081 stems ha-1 and a maximum underestimate of 6843 stems ha-1. The mean estimate of 756 stems ha-1 is an increase of only 30 stems ha-1 from the mean stem density of pine regeneration at the start of the simulations. The best estimates of stems ha-1 for advanced pine regeneration were also achieved through the 15-year transfer. However, variability in the estimates was high, with a maximum overestimate of 1892 stems ha-1 and a maximum underestimate of 1706 stems ha-1. The mean stems ha-1 of advanced pine regeneration from the 15-year simulations was estimated to be 1191 stems ha-1 and is an increase of 1038 stems ha-1 from the mean stem density at Time 1. Although the 15-year transfer provided good estimates of basal area for larger trees (≥ 7.5 cm DBH), it provided surprisingly poor estimates for stem density for the same size class. Conversely, the SORTIE-ND simulation provided very good estimates of the stem density of larger trees for all four species. Among the hybrid model simulations, the 5- and the 20-year transfer provided the best estimates of stem density for larger trees. With the 5-year transfer, densities for larger pine were projected to increase, starting from an average density of 260 stems ha-1 at Time 1 to 419 stems ha-1 at Time 3. Overall, however, it was somewhat reassuring that the 15-year transfer simulations performed reasonably well in terms of predicting the mean basal area of advanced pine regeneration as well as larger pine trees. The consistency in the biases suggest that this transfer point strikes the best balance between allowing the appropriate amount of understory trees to grow into the overstory and replacing those that did not with a new crop predicted from SORTIE-ND. By far the worst estimate for both basal area and stem density for advanced pine regeneration were obtained using the 5-year transfer. Here, the hybrid model greatly overestimated the basal area of advance pine regeneration. Furthermore, estimates were 65 quite variable. None of the simulations predicted any notable occurrence of new aspen regeneration. The mean biases for new Douglas-fir regeneration were lower; however, the variability was quite high. In some cases, underestimates were close to 8000 stems ha-1 while overestimates were near 3300 stems ha-1. For spruce, the 15 and 20 year tree-list transfer simulations provided the smallest mean bias (158 stems ha-1 and 125 stems ha-1, respectively). Projecting stands in SORTIE-ND for the entire simulation period provided similar results for trees ≤ 2.0 cm DBH. The average projected stem densities of aspen, Douglas-fir and spruce advanced regeneration and adult trees were very near observed densities for all hybrid model simulations, as well as when stands were run to completion in SORTIE-ND. For aspen and spruce, the variability in the biases among most stands was low. However, there were some notable exceptions for advanced aspen regeneration. For example, one stand continued to be largely underestimated over all four transfer simulations, reaching a maximum underestimate of 1132 stems ha-1 with the 20-year tree- list transfer simulations. For a second stand, none of the tree-list transfers provided predictions of advanced aspen regeneration after 25 years, when, in reality, the observed density was 300 stems ha-1. Few stems ha-1of advanced spruce regeneration was estimated throughout the 25 year projection period, regardless of when tree-lists were transferred. Projected changes from 1981 levels were minimal for advanced regeneration and large spruce. Among the stands that were simulated in all four tree-list transfers, advanced Douglas-fir regeneration was observed in only three stands. Thus, improvements resulting from any of the different tree-list transfer points were unclear. For the 5-year tree-list transfer, two of the stands with Douglas-fir had extremely large underestimates in the new regeneration size class and equally large overestimates for advanced regeneration. For these stands, it appeared that the growth rate of small trees was too high, resulting in an overestimate in of advanced regeneration that was juxtaposed against an underestimate of new regeneration. This trend was repeated for Douglas-fir in many of the stands in the 10-, 15- and 20-year transfer simulations. This same line of reasoning may also apply to the contrasting over- and underestimation that was described for new and advanced pine regeneration in the 5- and 10-year transfer simulations. Finally, SORTIE-ND and PrognosisBC both seemed to track the heights of trees quite well (Table 3.4). Within SORTIE-ND, the close tracking of tree heights is likely because the parameters for the height-diameter allometry equation were estimated using the Williams Lake dataset. Similarly, the equations used to relate height and diameter within PrognosisBC used 66 parameters calibrated for the IDF subzone within the Chilcotin Plateau. Looking at the four transfer simulations, there was a tendency to slightly overestimate the heights of smaller trees. No single transfer point provided a clear advantage over the others in terms of predicting heights for smaller and larger trees. Table 3.4 Mean bias and root mean square error (in parentheses) of tree height classes for four transfer simulations and independently run SORTIE-ND simulations Tree Height (Ht) 2m≤Ht<5m 5m≤Ht<8m 8m≤Ht<11m 11m≤Ht<14m 14m≤Ht<24m 5yr Transfer Aspen -0.22 (0.53) -0.23 (0.76) 0.92 (1.02) -0.02 (0.31) -1.25 (1.49) Douglas-fir 0.27 (1.25) -0.28 (0.54) -0.03 (0.97) -0.61 (0.69) 0.20 (0.54) Lodgepole pine -- -0.72 (0.86) -0.25 (0.71) -0.13 (0.66) 0.14 (0.66) Hybrid spruce -0.23 (0.67) -0.17 (1.05) -0.48 (1.10) -0.51 (1.04) -0.96 (1.00) 10yr Transfer Aspen -0.46 (0.57) -0.32 (0.70) -0.31 (0.34) 0.05 (0.43) -1.32 (1.53) Douglas-fir -- 0.44 (0.76) -0.08 (0.49) -0.57 (0.60) 0.16 (0.49) Lodgepole pine -1.38 (1.43) -0.51 (0.82) -0.41 (0.85) 0.23 (0.61) 0.16 (0.69) Hybrid spruce 0.27 (0.61) 0.31 (0.86) -0.51 (1.08) -0.43 (1.03) -1.01 (1.03) 15yr Transfer Aspen -0.43 (1.31) -0.96 (1.12) -0.16 (0.21) 0.11 (0.38) -1.55 (1.67) Douglas-fir -- -0.1 (0.59) 0.16 (0.46) -0.69 (0.76) -0.1 (0.39) Lodgepole pine -0.59 (0.66) -0.31 (0.80) -0.37 (0.81) 0.28 (0.64) 0.1 (0.71) Hybrid spruce -0.11 (0.62) 0.39 (0.84) -- -- -- 20yr Transfer Aspen 0.83 (1.03) -0.84 (1.28) 0.1 (0.14) 0.14 (0.38) -1.56 (1.68) Douglas-fir 1.44 (1.49) -0.46 (0.57) 0.15 (0.43) -0.63 (0.72) -0.1 (0.39) Lodgepole pine -0.66 (0.74) -0.14 (0.57) -0.28 (0.77) 0.41 (0.75) 0.13 (0.72) Hybrid spruce -0.23 (0.24) 1.67 (1.69) -0.55 (0.59) -- -- SORTIE-ND Aspen 0.47 (0.49) -0.41 (0.99) 0.14 (0.16) -0.29 (1.52) -1.87 (1.91) Douglas-fir 1.44 (1.48) -0.1 (0.22) 0.16 (0.18) -0.41 (0.48) -0.44 (0.47) Lodgepole pine -0.66 (0.74) -0.14 (0.58) -0.32 (0.82) 0.45 (0.61) 0.21 (0.79) Hybrid spruce -0.18 (0.22) -- -0.55 (0.87) -- -- 67 3.2 Ecological expectations 3.2.1 Density of seedlings and GLI For all stands simulated using the SORTIE-ND model, the predicted stems ha-1 of pine seedlings peaked between 5- and 10-years after the start of the simulations. The mean predicted density of pine seedlings after five time steps in SORTIE-ND was 3463 and 2711 stems ha-1 after 10 years. The minimum height considered by SORTIE-ND is 10cm. Given that a constant seed rain function was used to add new pine seeds, examining changes in the stems ha-1 of trees below 20cm height should give an indication of how the seedling establishment behaviour is performing. After five years, the density of pine seedlings less than 20 cm in height was just over 1000 stems ha-1 and represented nearly 25% of the total number of trees below 1.35m. This decreased to 350 stems ha-1 and a 50% share of the seedling population from 15 years onward. Not surprisingly, peaks in pine seedling stem density lagged neatly behind peak understory GLI values, which is partially a result of the indirect relationship between seedling mortality and GLI (Kobe and Coates, 1997; Coates and Hall, 2005). Shortly after GLI levels began to decline, periodic growth decreased, thereby increasing the probability of mortality among pine seedlings and saplings. When understory GLI values were below 20%, the predicted periodic growth of pine seedlings was nearly zero for all simulated stands (Figure 3.3). Since seedling mortality is predicted as a function of diameter increment, these seedlings would have been at an increased probability of mortality. Decreasing understory GLI levels had a lower impact on the periodic growth of pine saplings. However, this seemed to only delay slightly the onset of light mediated mortality, as densities for all understory pine dropped rapidly. Conversely, spruce seedlings and saplings were predicted to have higher growth rates, even at GLI levels below 20% . Seedling densities for spruce remained at a low, yet constant level throughout the simulation period within the SORTIE-ND model. After 5 or 10-year time steps, the stem density of spruce saplings began to slowly increase. The simultaneous decrease of pine saplings and increase of spruce saplings was a reflection of the interspecific differences in the competition for light that is built into the SORTIE-ND growth and mortality models. 68 Figure 3.3 Average periodic growth (cm / year) from the 7 simulated stands for pine (a and b, respectively) and spruce (c and d, respectively) seedlings (trees < 2 m tall) versus predicted global light index when simulations were run to completion in SORTIE-ND. Average periodic growth rates from each stand are plotted in 5 year intervals over a 25 year simulation. a) b) c) d) 69 The expectation was that understory GLI values would be closely tied to estimates of crown radius. There was a predicted increase in the mean crown width of adult trees as simulations progressed, from 1.5m in 1981 to 1.75m in 2006. Likely, this was in response to the predicted increase in adult DBH and height. Further increases in crown width were countered by the tendency for stands to experience minor increases in large tree basal area and stems ha-1. Under the new density-dependent large tree crown allometry equations, changes in crown width and crown height are a function of individual tree size dimensions (DBH, height, and crown size) as well as stems ha-1, basal area (m2 ha-1), and basal area of larger trees) (Chapter 2). It seems possible that an average increase of 0.25 m at the crown base of large trees was sufficient to cause predicted understory GLI values to drop by an average of 44% over 25 years. The predicted increase in crown widths would also have reduced the number of pine seeds that were converted to seedlings, since the establishment of seedlings is a partially a function of the amount of forest cover above the seeds that have been dispersed. Using this rationale provides an explanation to the sharp decrease seen in the number of pine seedlings that were below 20cm in height. 3.2.2 Growth rates of regenerated trees Predicted periodic growth rates for pine less than 2cm DBH were higher in SORTIE-ND than PrognosisBC after 5 years (0.49cm/yr vs. 0.14 cm/yr), but were very similar for the remainder of the simulation period (e.g., 0.17cm/yr in SORTIE-ND and 0.11cm/yr in Prognosis at year 15). After five years, trees between 2 and 4 cm DBH were predicted to have a periodic growth rate of 0.26cm/yr in SORTIE-ND and 0.33cm/yr in PrognosisBC. From ten years onward, predicted periodic growth rates were noticeably different (0.12cm/yr for SORTIE-ND and 0.32cm/yr for PrognosisBC). Differences were observed even sooner for pine between 5.5 and 7.5cm DBH (Figure 3.4). Within the hybrid model, these differences resulted in abrupt increases in periodic growth when trees are transferred from SORTIE-ND to PrognosisBC. This abrupt increase is shown in Figure 3.4, which depicts the predicted periodic growth trends observed in the 15 year tree-list transfer simulations for pine trees between 5.5 and 7.5cm DBH. 70 Figure 3.4 Predicted periodic growth (cm/yr) for pine between 5.5 and 7.5cm DBH in SORTIE-ND and PrognosisBC using the 15 year tree-list transfer. The gap in the periodic growth trends marks when understory trees were transferred from SORTIE-ND to PrognosisBC. Each trend line is the predicted growth from a single stand. Symbols on the trend lines identify individual stands, which are run in both SORTIE-ND and PrognosisBC Given the drop in estimated understory GLI, the predictions within SORTIE-ND are as expected. However, the predicted patterns do not seem to match with field measurements collected by Heath and Alfaro (1990). From their records, growth rates of surviving understory pine remained higher than average growth rates prior to attack. This increased growth rate continued to be observed 14 years after MPB disturbance. Conversely, 15 years after the start of the simulations, mean predicted diameter growth rates for understory pine using SORTIE-ND were below expected norms for stands undisturbed by MPB. These discrepancies emphasize the need to parameterize SORTIE-ND’s juvenile growth to the conditions of the study area. One reassuring part about the predicted growth estimates was that peak growth was observed 5 to 10 years after the start of the simulations, regardless of size class. A similar response interval was reported by Heath and Alfaro (1990), who measured a peak increase in pine growth 5 to 9 71 years following MPB disturbance. This lends credence to the assumption that reconstructed stand conditions at the beginning of the hybrid model simulations were representative of conditions that would have been present shortly after MPB attack. 3.4 Discussion The main objective of this study was to examine the functionality of the hybrid model approach and to evaluate performance over a medium-term planning horizon. The importance of evaluating the functionality cannot be understated if the intent is for the model to be used in a forest planning capacity. On this issue, the linkages established between SORTIE-ND and PrognosisBC resulted in a reasonably functional hybrid model. With some additional programming, the process of exporting tree-lists from SORTIE-ND, reformatting them, and importing them into an existing PrognosisBC simulation could be fully automated. One difficulty is the exchange of tree-list information between the models. For SORTIE-ND, the tree-list must be created with an expansion factor of 1 for each tree. Trees transferred to PrognosisBC had to be collapsed into small size classes prior to being imported. Developing this ‘work-around’ has cleared the way to begin the shift from a hybrid model that uses static linkages to a dynamic process where tree-lists are exchanged throughout the simulation process. Testing the hybrid model approach against mensurational data provided some encouraging results. The average biases for the basal area of advanced regeneration and large trees were low for the three secondary species and within reasonable limits for advanced pine regeneration whilst using the 15-year tree-list transfer simulations. The predicted increase in pine basal area from the beginning of the simulations (45%) was slightly higher than the 33% relative increase in basal area reported by Coates and Hall (2005) over a similar time frame following their simulation of understory pine trees in MPB attacked stands in the SBS zone. Average biases were only slightly higher for large lodgepole pine basal area. The relative increase in large tree basal area estimated through the 15-year tree-list transfer simulations were slightly larger than others reported in the literature. For example, the mean predicted basal area from the 15- year tree-list transfer simulations showed an increase of 26% from the basal area of large pine at Time 1 (i.e., a 25 year period). However, only ten years after MPB infestation, Stockdale et al. (2004) already noted a 22% increase in adult pine basal area in stands located on the Chilcotin plateau, with trends suggesting that further increases would likely be observed. Estimates for the stems ha-1 of new pine regeneration showed high mean biases with respect 72 to the mensurational data, but showed both agreement and confliction with reports from literature. For example, the rapid decrease in seedling densities seem to be in agreement with findings from simulations run by Coates and Hall (2005), which suggested that even with understory planting, understory light levels following MPB attack can be limiting and result in extremely high rates of mortality in pine seedlings. However, Coates and Hall (2005) went on to note that there were no studies to corroborate their findings with respect to changing understory light conditions following MPB attack. Results of the simulations run for this study are in stark contrast to the observed densities in the Williams Lake dataset as well as densities reported by Hawkes et al. (2004) who re-measured stands roughly 16 years after MPB attack. In the stands sampled by Hawkes et al. (2004), the density of trees <1.5m height was 4538 stems ha-1 (a decrease from 4687 stems ha-1 14 years prior). From the conflicting evidence presented here, one might suggest that calibrating the SORTIE-ND light transmission parameters for conditions in MPB disturbed stands is needed before further testing of the hybrid model approach is conducted. Estimates of basal area and stems ha-1 of lodgepole pine were generally consistent among the stands used in the simulations; however, the outliers raised some concerns regarding the use of the hybrid model approach at the stand level. In its current state, the hybrid model approach appears to be better suited for application at the landscape level where less resolution is required. For example, the hybrid model approach could be used to test the effects of resource management strategies on the mid-term timber supply. To downscale the hybrid model for application at the stand level would likely require parameterization of SORTIE-ND and PrognosisBC at this smaller scale. For this study, the limited amount of field data meant that most SORTIE-ND parameters had to be obtained from external sources. For the parameters that could be estimated from the Williams Lake dataset, tree records had to be pooled from three different BEC subzones in order to have a sufficiently large sample size. As a result, resolution at the stand level was sacrificed so that a reliable set of parameters useful at the landscape level could be obtained. Parameterizing these equations to specific BEC zones, variants, or even to individual stands would most likely reduce prediction biases. Calibrating the PrognosisBC model to a wider range of BEC units would also be necessary and would have a similar effect on reducing prediction biases. It is important to recall that natural regeneration is highly variable at even the finest scales. Several studies have tried to quantify and characterize patterns of natural regeneration to varying degrees of success (Steijlen et al. 1995; Ehle and Baker 2003; Zagidullina and Tikhodeyeva 2006). The complex interaction of seeds and seedlings with parent trees, substrate availability, soil moisture and light presents significant challenges to the modeler. The regeneration submodel 73 in SORTIE-ND attempts to simulate many of these conditions, and in theory, this should improve the chances of obtaining accurate and reliable estimates of regeneration. Put into practice, the resulting biases and variability suggest that there is room for improvement. However, given that the seed dispersal rates for pine were parameterized from the literature, and most other parameters obtained from an external source, the results seem promising. The evaluation of the hybrid model approach on the basis of ecological expectations proved very insightful. When the models were run independent of each other (i.e., from Time 1 to Time 2), they generally behaved as expected. In SORTIE-ND, peaks in seedling densities matched with peaks in growth rates, which in turn were closely correlated with high understory light levels. Within the first 10 years of the simulations, these trends were in accordance with reports from literature. However, sharp drops in the stems ha-1 of understory trees in the latter half of the simulations were somewhat unexpected, though explainable within a modeling context. That is, within the SORTIE-ND model, predicted decreases in understory light translated to increased rates of mortality among shade intolerant juvenile trees. In reality, pine growing on dry or well- drained sites under an overstory have been shown to display a higher tolerance to low levels of light (Waring and Pitman 1985; Hawkes et al. 2004). It seems likely that parameterizing the juvenile growth function using data collected from drier sites such as those found in the Chilcotin Plateau would result in model predicted growth rates that would be less responsive to moderate decreases in understory GLI, thereby reducing the rate of mortality among juvenile pine trees. Sharp drops in seedling densities may also be explained through the randomization of the location of trees imported into SORTIE-ND. By randomly assigning a coordinate, understory trees that were at one point in a favourable location within the stand may have been relocated to a less favourable location (e.g., one with lower light levels). Using spatially referenced data would almost certainly improve predictions for complex stands experiencing gap-level disturbance events (Rouvinen and Kuuluvainen 1997). However, obtaining spatially reference data is expensive and the full benefits would not be realized within the hybrid model approach since PrognosisBC is spatially independent. There remains some uncertainty with regards to the seed dispersal and seedling establishment parameters. Using a non-spatial seed dispersal parameter obtained from the literature seemed to move model predicted densities of pine seedlings closer to observed values, at least for the first 10 years of the simulations. The current explanation put forth to rationalize the observed decrease in the stems ha-1 of pine seedlings less than 20cm in height is that the increase in average crown width, though moderate, was enough to decrease the rate at which seeds were converted to 74 seedlings. The assumption with this rationale is that seedlings less than 20cm in height likely established in the previous time step and mortality due to slow growth would require at least one or two more growth cycles before its effects could be seen. Therefore, predicted decreases in the stem density of trees less than 20cm could be largely attributed to the establishment parameters. However, until the juvenile growth model is re-parameterized to drier sites, one can only speculate on the degree to which the current set establishment parameters affect seedling densities. When tree-lists were transferred to PrognosisBC, unexpected increases in periodic radial growth were observed, particularly among larger understory trees. These increases are relative to the expected growth rates had the model been allowed to continue. These increases are likely the result of a rearrangement in the social rank of understory trees within the imported tree-lists. In PrognosisBC, the social rank of trees affects the prediction of periodic growth rates. The social rank of a tree is based on the size of a tree, relative to other trees in the stand. Trees with a higher social rank have increased growth and continue to grow well in future growth cycles. Conversely, smaller trees which have a lower social rank are predicted to have slower growth rates. Differential growth rates for seedlings and saplings in SORTIE-ND likely resulted in an understory with a social ranking structure different from that of the understory in PrognosisBC. The result is that minor, but unexpected changes to the growth rate of some trees are likely to occur when understory trees are transferred. For example, in SORTIE-ND the higher growth rates predicted for pine trees between 5.5 and 7.5cm DBH appears to have increased their social rank relative to say, trees less than 2cm DBH. Within PrognosisBC, this gain in social rank would have resulted in a growth rate higher than expected. Improving juvenile growth rates in SORTIE-ND may reduce this problem. The secondary objective of this study was to identify a point in the simulated development of the stands that best reflected conditions characteristic of the end of understory reinitiation. Among the four transfer times tested, the exchange of understory tree-lists after 15 years appeared to provide the best balance in terms of minimizing mensurational biases across all species and size classes. Within the simulations, it appears that the availability of understory growing space and resources following MPB attack has become fully utilized after 15 years. This is an interesting hypothesis generated by this study that would require further field testing. As it is used here, the hybrid model approach was used to provide PrognosisBC with a naturally regenerating understory defined as trees less than 7.5 DBH. This, of course, is a flexible definition of understory trees. Depending on the type of stands being simulations, alternative size 75 classes may be used so long as trees of the same size class are removed from PrognosisBC prior to the import of new trees. For example, stands that experience a high amount of crown rise may have more vertical growing space for advanced regeneration to fill, resulting in the development of larger understory trees at the end of the understory reinitiation phase. In this case, the use of larger size class limits may be desired. 3.5 Conclusions The hybrid modelling approach presented and tested here is a work in progress and further testing is required. One of the priorities for the next round of testing should be to automate the linkages between SORTIE-ND and PrognosisBC and evaluate the benefits of using a dynamic feedback system. Developing a dynamic feedback system would distinguish the hybrid model approach from the existing static imputation techniques currently used by PrognosisBC. Based on the results of this study, further investment into the hybrid model approach seems warranted. Silvicultural decisions are often based on stand basal area. Thus, it is of utmost importance that growth models used to inform management decisions at a minimum provide accurate and reliable estimates of stand basal area. Held in this context, the hybrid model approach appears to be a promising and worthwhile tool that forest managers could use to aid in the development of mid- term harvest plans for stands that have been disturbed by MPB. 3.6 Literature cited Alfaro, R.I, Campbell, R., Vera, P., Hawkes, B., and Shore, T.L. 2004. Dendroecological reconstruction of mountain pine beetle outbreaks in the Chilcotin Plateau of British Columbia. In Mountain pine beetle symposium: challenges and solutions. October 30-31, 2003, Kelowna, British Columbia. T.L. Shore, J.E. Brooks, and J.E. Stone (editors). Nat. Res. Can., Can. For. Serv., Pac. For. Centre, Inform. Rep. BC-X-399, Victoria, BC. pp.245-256. Astrup, R. 2006. Modeling growth of understory aspen and spruce in western boreal Canada. Ph.D. Dissertation. University of British Columbia, Vancouver, BC. 156 p. Astrup, R., Coates, D.K., Hall, E., and Trowbridge, A. 2007. Documentation of the SORTIE-ND SBS research parameter file version 1.0. February, 2007. Bulkley Valley Centre for Natural Resources Research and Management, Smithers, B.C. 76 Aukema, B. H., Carroll, A. L., Zhu, Jun., Raffa, K. F., Sickley, T.A., and Taylor, S.W. 2006. Landscape level analysis of mountain pine beetle in British Columbia, Canada: spatiotemporal development and spatial synchrony within the present outbreak. Ecog. 29: 427-441. Bainbridge, E.L. and Stong, Wl.L. 2005. Pinus contorta understory vegetation dynamics following clearcutting in west-central Alberta, Canada. For. Ecol. and Manage. 213:133- 150 British Columbia (BC) Forest Practices Board. 2007. Lodgepole pine stand structure 25 years after mountain pine beetle attack. Forest Practices Board, Victoria, BC FPB/SR/32 British Columbia (BC) Ministry of Forests Research. 1999. Uniform shelterwood systems for Douglas-fir and lodgepole pine stands: Update, year 8. Extension note 28, April 1999. Research section, Ministry of Forests, Cariboo Forest Region. British Columbia (BC) Ministry of Forests and Range. 2006. Site index estimates by site series: Report by Biogeoclimatic unit (2006 Approximation). Forest Science Program. Campbell, E. M., Alfaro, R. I., and Hawkes, B. 2007. Spatial distribution of mountain pine beetle outbreaks in relation to climate and stand characteristics: A dendroecological analysis. J. Integr Plan. Biol. 49: 168-178. Canham, C.D., Coates, K.D., Bartemucci, P., and Quaglia, S. 1999. Measurement and modeling of spatially explicit variation in light transmission through interior cedar-hemlock forests of British Columbia. Can. J. For. Res. 29: 1775-1783. Coates, K.D. and Hall, E.C. 2005. Implications of alternative silviculture strategies in Mountain Pine Beetle damaged stands. Technical report for Forest Science Program Project Y051161. Bulkley Valley Centre for Natural Resource Research and Management, Smithers, BC. 31p. Dordel, J., Feller, M.C., and Simard, S.W. 2008. Effects of mountain pine beetle (Dendroctonus ponderosae Hopkins) infestations on forest stand structure in the southern Canadian Rocky Mountains. For. Ecol. and Manage. 255: 3563-3570 77 Ehle, D.S. and Baker, W.L. 2003. Disturbance and stand dynamics in ponderosa pine forests in Rocky Mountain National Park, USA. Ecol. Monog. 73: 543-566 Environment Canada 2007. Canadian daily climate data CDs. (1 July - Aug 31 2002 - 2007). http://www.climate.weatheroffice.ec.gc.ca/prods_servs/cdcd_iso_e.html, accessed November, 2008. Griesbauer, H. and Green, S. 2006. Examining the utility of advanced regeneration for forest reforestation and timber production in unsalvaged stands killed by the mountain pine beetle: Controlling factors and management implications. BC Journal of Eco. and Manage. 7(2): 81-92. Hassani, B., LeMay, V.M., Marshall, P.L., Temesgen, H., and Zumrawi, A.A. 2004. Regeneration imputation models for complex stands of Southeastern British Columbia. The Forestry Chronicle 80:271-278. Hawkes, B., Taylor, S., Stockdale, C., Shore, T.L., Alfaro, R.I., Campbell, R., and Vera, P. 2004. Impacts of mountain pine beetle on stand dynamics in British Columbia. Pages 175-195 in T.L. Shore, J.E. Brooks, and J.E. Stone, editors. Mountain Pine Beetle Symposium: Challenges and Solutions. Kelowna, British Columbia October 30-31, 2003, Nat. Res. Can., Can. For. Serv., Pac. For. Centre, Inform. Rep. BC-X-399, Victoria, BC. 298p. Hawkes, B.C., Taylor, S.W., Stockdale, C. , Shore, T.L., Beukema, S.J., and Robinson, D.C.E. 2005. Predicting mountain pine beetle impacts on lodgepole pine stands and woody debris characteristics in a mixed severity fire regime using PrognosisBC and the Fire and Fuels Extension. pp. 123-135. In: Lagene, L., J. Zelnik, S. Cadwallader and B. Hughes (eds.). Mixed Severity Fire Regimes: Ecology and Management. November 17-19, 2004, Spokane, Washington. Washington State University Coop Extension Service/The Association for Fire Ecology, Pullman, Washington, Vol. AFE MISC03. Heath, R. and Alfaro, R.I. 1990. Growth response in a Douglas-fir/lodgepole pine stand after thinning of lodgepole pine by the mountain pine beetle: a case study. Journal of the Entomological Society of British Columbia 87:16-21. Kobe, R.K., and Coates, K.D. 1997. Models of sapling mortality as a function of growth to characterize interspecific variation in shade tolerance of eight tree species of northwestern British Columbia. Can. J. For. Res. 27:227-236. 78 LeMay, V., Lee, T., Sattler, D., Marshall, P., Robinson, D., and Zumrawi, A-A. 2007. Modeling Natural Regeneration Following Mountain Pine Beetle Attacks in the Southern and Central Interior of British Columbia: Final report. Internal report for Natural Resources Canada, MPB Standard Contribution Agreement, PO # 8.35. 57 pp. http://www.forestry.ubc.ca/prognosis/extension.html accessed November, 2008. Oliver, C.D. and Larson, B.C. 1996. Forest stand dynamics. John Wiley and Sons, Inc. ISBN 0- 471-13833-9, New York, 520 p. Mäkelä, A. 2003. Process-based modelling of tree and stand growth: towards a hierarchical treatment of multiscale processes. Can. J. of For. Res. 33: 398-409. Meidinger, D. and Pojar. J. 1991. Ecosystems of British Columbia. B.C. Min. For., Res. Br., Victoria, B.C. Spec. Rep. Ser. No. 6. Nigh, G. 2002. Site index conversion equations for mixed trembling aspen and white spruce stands in northern British Columbia. Silva Fenn. 36: 789-797 Patriquin, M.N., Lantz, V.A., Stedman, R.C., and White, W.A. 2008. Working together: a reciprocal wood flow arrangement to mitigate the economic impacts of natural disturbance. Forestry 81: 227-242. Reinhardt, E. and Crookston, N.L. , eds.. 2003. The fire and fuels extension to the forest vegetation simulator. USDA For. Serv. Gen. Tech. Rep. RMRS-GTR-116. 209 pp. Robinson, D.C.E. 2004. PrognosisBC southern interior geographic variant description: Interior Douglas-fir and interior cedar hemlock zones. Prepared by ESSA Technologies Ltd., Vancouver, BC for British Columbia Ministry of Forests, Victoria, BC. 42 pp. Robinson, A.P. and Ek, A.R. 2003. Description and validation of a hybrid model of forest growth and stand dynamics for the Great Lakes region. Ecol. Modelling 170: 73-104. Robinson, A.P. and Monserud, R.A. 2003. Criteria for comparing the adaptability of forest growth models. For. Ecol. and Manage. 172: 53-67. Rouvinen, S. and Kuuluvainen, T. 1997. Structure and asymmetry of tree crowns in relation to local competition in a natural mature Scots pine forest. Can. J. For. Res. 27: 890-902. 79 SAS Institute Inc. 2003. SAS/ETS user’s guide, version 9.1.3. Cary, NC: SAS Institute Inc. Stage, A.R. 1973. Prognosis model for stand development. USDA, For. Serv. Res. Pap. Intermountain Forest and Range Exp. Station. INT-137. Ogden, UT. Steijlen, I., Nilsson, M.C. and Zackrisson, O. 1995. Seed regeneration of Scots pine in boreal forest stands dominated by lichen and feather moss. Can. J. of For. Res. 25: 713-723 Stockdale, C., Taylor, S. and Hawkes, B. 2004. Incorporating mountain pine beetle impacts on stand dynamics in stand and landscape models: a problem analysis. In Mountain pine beetle symposium: challenges and solutions. October 30-31, 2003, Kelowna, British Columbia. T.L. Shore, J.E. Brooks, and J.E. Stone (editors). Nat. Res. Can., Can. For. Serv., Pac. For. Centre, Inform. Rep. BC-X-399, Victoria, BC. pp.200-209. Temesgen, H. and Gadow, K. 2004. Generalized height-diameter models - an application for major tree species in complex stands of interior British Columbia. Euro. J. of For. Res. 123: 45-51 Vanclay, K. J. 1994. Modelling forest growth and yield: Applications to mixed tropical forests. CABI, New York. pp. 103-133. Valentine, H. T. and Mäkelä, A. 2005. Bridging process-based and empirical approaches to modeling tree growth. Tree. Physiol. 25: 769-779. Waring, R.H. and Pitman, G.B. 1985. Modifying lodgepole pine stands to change the susceptibility to mountain pine-beetle attack. Ecology 66: 889-897. Wood, C. S. and Unger, L. S. 1996. Mountain pine beetle: a history of outbreaks in pine forests in British Columbia, 1910 to 1995. Nat. Res. Can., Can. For. Ser., Pac. For. Centre, Victoria, B.C. Zagidullina, A. and Tikhodeyeva, M. 2006. Spatial patterns of tree regeneration and ground cover in dry Scots pine forest in Russian Karelia. Ecoscience 13:203-218. Zumrawi, A.A., LeMay, V., Temesgen, H., Snowdon, B. 2003. Expanding PrognosisBC to the Chilcotin variant of the interior Douglas-fir dry cool subzone: Large tree radial growth models. Report prepared for: Forestry Innovation Investment, Project No: FIIRFP RES2002/03. 80 Zumrawi, A.A., LeMay, V., Marshall, P.M., Hassani, B.T. and Robinson, D. 2005. Implementing a PrognosisBC regeneration sub-model for complex stands of Southeastern and Central British Columbia. Collaborative Project between Research Branch, Ministry of Forests. Forest Resources Management Department, University of British Columbia, and ESSA Technologies, Inc. Vancouver, British Columbia, Canada. Report submitted to: Forest Science Program, Project No: Y051355, 157 p. http://www.for.gov.bc.ca/hre/pubs/azumrawi.htm accessed, June 2007. 81 4. Concluding Remarks 4.1 General discussion and conclusions This study had two main objectives: to develop a model to predict overstory tree crown size in stands of varying density and to build a growth model capable of estimating natural regeneration in unsalvaged MPB-disturbed forests. The first objective was accomplished by predicting crown depth and crown radius through a nonlinear system of equations. The crown models took into account important biological factors that influence crown size, namely structural dependencies between tree size and crown size, and the relationship between stand density and crown size. For the second objective, two disparate forest growth models were linked together to create a hybrid model. The first, SORTIE-ND, was used to estimate natural regeneration as a function of available light, seed dispersal and seed substrate availability. The second, PrognosisBC, received estimates of regeneration for SORTIE-ND, and projected individual tree growth as a function of tree size, site variables and the target tree’s social position within the stand. The two objectives were connected through stand dynamic processes simulated within the hybrid model. For MPB-disturbed stands that go unsalvaged, regeneration will occur under highly variable understory light conditions. For many of these stands, understory light is a key predictor of natural regeneration and growth (Hawkes et al. 2004; Coates and Hall, 2005). The variability arises from changes in overstory tree crown structure. Within SORTIE-ND, simulated light conditions in the understory are a function of crown size, which is defined through estimates of crown depth and crown radius. Thus, the introduction of density-dependent crown allometry models to SORTIE-ND was conducive to the application of the hybrid model for the purpose of predicting natural regeneration in MPB-disturbed stands of varying overstory density. 4.2 Contributions and criticism 4.2.1 Density-dependent crown models The density dependent crown models were developed with three criteria in mind: 1) select an equation form suitable to the biological processes influencing crown size; 2) ensure logical consistency in the estimates; and 3) include indices of plant competition as predictors of crown length and crown radius. Chapter 2 described how these three criteria were met using a nonlinear 82 three-stage least squares (N3SLS) estimator to fit logistic and power statistical models for crown depth and crown radius, respectively. In the field of forest biometrics, fitting equations using a N3SLS estimator is not uncommon. To my knowledge, there were no previous examples in the scientific literature where this parameter estimation method was used to predict tree crown size. Biologically, using a N3SLS estimator made sense as it allowed two strongly interdependent variables, crown depth and crown radius, to be used to improve the prediction of the other axis. Based on the fit statistics reported in Chapter 2, the density-dependent crown models were appropriate for aspen, Douglas-fir, lodgepole pine and hybrid spruce within the study area. A true validation dataset was not available for this study. As an acceptable substitute, n-way validation using a PRESS statistic was used. This provided support for the internal validity of the model. Validation beyond the sample population is still needed (Rykiel 1996; Yang et al., 2004). Furthermore, for small sample sizes since estimators of parameters and associated standard errors are only asymptotically unbiased efficient for nonlinear models, it is recommended that bootstrap resampling be used instead of fitting the equations through a single N3SLS estimation. The variables used to describe stand density in the new crown equations were supported by biological principles and statistical diagnostics. Also, it was argued that using simple measures of stand density would be more reliable over a range of stand densities. However, using simple measures of stand density has potential disadvantages. For example, stem density does not factor in the effects of tree size. This could lead to unexpected or incorrect conclusions about the effect of crowding on crown size. For the sample population, the contributions of the density-related variables to the overall model fit were as expected; that is, crowns were smaller in stands with more crowding and larger in stands with less crowding. External validation of the crown models should reveal if the potential disadvantages associated with the simple measures of density exist among other populations. 4.2.2 A hybrid model constructed by linking SORTIE-ND and PrognosisBC Overall, the results of Chapter 3 showed that the SORTIE-ND and PrognosisBC models were suited to their roles within the hybrid model. The list of biotic and abiotic factors potentially influencing natural regeneration in MPB-disturbed stands is lengthy (Griesbauer and Green 2006). For trees of advance-regeneration size, the list of processes simulated within the hybrid model provided good estimates of basal area over a 25-year horizon. However, for smaller regenerated trees, estimates of density were generally poor. The evaluation of ecological expectations revealed that the poor estimates for smaller-sized regeneration may be due to the 83 parameters used in the small tree growth model in SORTIE-ND, which resulted in predicted periodic growth rates that were lower than expected. Full parameterization of the SORTIE-ND model to the conditions within the sample population is needed to gain further insight into the factors affecting the growth of new regeneration. Few, if any, tools are available to predict natural regeneration following MPB infestation (Zumrawi et al. 2005). Therefore, the advancements made in the development of the hybrid model described in Chapter 3 provide important contributions to the current body of knowledge on growth and yield modelling in stands disturbed by MPB. The strategy of linking SORTIE-ND and PrognosisBC proved to be a viable alternative to using the models on their own. This was particularly true for PrognosisBC, which prior to the hybrid model, did not have a regeneration submodel well-suited for use in MPB-disturbed stands. For SORTIE-ND, the gains are less clear since many of the model parameters were obtained from outside the study area. As was previously noted, full parameterization of SORTIE-ND to the sample population is needed and should clarify this issue. Demonstrating that two disparate growth models can be efficiently linked to take advantage of each model’s strengths hopefully encourages others to view existing growth and yield models not simply as independent units, but as potentially interchangeable parts. The hybrid models developed by Robinson and Ek (2003) and Milner et al. (2003) are two other successful examples of this approach to modelling. Furthermore, linking two models with different approaches to modelling growth and forest dynamics should help to increase communication among researcher and forest managers using these tools. 4.3 Management implications Models that can predict natural regeneration in MPB-disturbed forests have taken on greater importance as more questions arise regarding contributions from unsalvaged stands to the mid- term timber supply. Currently, there is little long-term information on forest growth and dynamics from unsalvaged stands following MPB-disturbance. The hybrid model could be used to simulate a variety of stand conditions following MPB infestation. The results could fill some the knowledge gaps and help forest managers make more informed decisions on several issues. For example, it is anticipated that many unsalvaged stands will not meet current stocking standards (Griesbauer and Green 2006). For these stands, underplanting or overstory thinning may be helpful in stimulating new regeneration and releasing advance regeneration. However, due to the 84 variability in processes affecting natural regeneration, it is unclear which stand types are likely to need silvicultural intervention. The hybrid model is capable of tracking processes affecting advance regeneration and therefore, could be used to identify unsalvaged stands likely to achieve stocking standards and those that will require some form of management. 4.4 Future directions To realize the potential contributions of the hybrid model described in this thesis, further work is needed. One of the main priorities is the parameterization of SORTIE-ND behaviours to MPB-disturbed stands located on the Chilcotin plateau. This area of BC harbours a substantial number of MPB-disturbed stands. Thus, the cost of collecting field measurements to parameterize the model seems warranted. Following this, external validation of the density-dependent crown models and the hybrid model is needed. This step is critical as it will identify any weaknesses and test the robustness of the hybrid model. A second priority for future work is the development of a dynamic feedback system between SORTIE-ND and PrognosisBC. Under the current system, PrognosisBC is supplied with estimates of regeneration from SORTIE-ND at one point in time. This system was appropriate for the objectives of this study. However, for longer projection periods, SORTIE-ND needs to provide estimates of regeneration to PrognosisBC at multiple points in time. Furthermore, tree-lists in SORTIE-ND would need to be updated to match tree lists in PrognosisBC. This way, overstory composition and structure in the two models would remain similar throughout the projection period. Lastly, there are some technical aspects of the hybrid model that need to be addressed. The simulation times for the 25-year projections in SORTIE-ND were very long. The cause of the long simulation times appears to be related to the new crown models; however, the exact reasons remain unknown. Also contributing to lengthy simulation times were snags in the tree lists used to initialize SORTIE-ND. To overcome this, snags were removed. If simulation times could be reduced when snags are included in the initial tree list, then the hybrid model could be applied to a wider set of starting stand conditions. Working with the SORTIE-ND developers will hopefully identify the related problems and reduce the computing time. 85 4.5 Literature cited Coates, K.D. and Hall, E.C. 2005. Implications of alternative silviculture strategies in Mountain Pine Beetle damaged stands. Technical report for Forest Science Program Project Y051161. Bulkley Valley Centre for Natural Resource Research and Management, Smithers, BC. 31p. Griesbauer, H. and Green, S. 2006. Examining the utility of advance regeneration for reforestation and timber production in unsalvaged stands killed by the mountain pine beetle: Controlling factors and management implications. BC J. of Ecosystems and Manage. 7: 81-92. http://www.forrex.org/publications/jem/ISS35/vol7_no2_art9.pdf Hawkes, B., Taylor, S., Stockdale, C., Shore, T.L., Alfaro, R.I., Campbell, R., and Vera, P. 2004. Impacts of mountain pine beetle on stand dynamics in British Columbia. Pages 175-195 in T.L. Shore, J.E. Brooks, and J.E. Stone, editors. Mountain Pine Beetle Symposium: Challenges and Solutions. Kelowna, British Columbia October 30-31, 2003, Nat. Res. Can., Can. For. Serv., Pac. For. Centre, Inform. Rep. BC-X-399, Victoria, BC. 298p. Milner, K.S., Coble, D.W., McMahan, A.J. and Smith, E.L. 2003. FVSBGC: a hybrid of the physiological model STAND-BGC and the forest vegetation simulator. Can. J. of For. Res. 33: 466-479. Robinson, A.P. and Ek, A.R. 2003. Description and validation of a hybrid model of forest growth and stand dynamics for the Great Lakes region. Ecol. Modelling 170: 73-104. Rykeil, E.J.Jr. 1996. Testing ecological models: the meaning of validation. Ecol. Modelling 90: 229-244. Yang. Y., Monserud, R.A., and Huang S. 2004. An evaluation of diagnostics tests and their roles in validating forest biometric models. Can. J. For. Res. 34: 619-629. Zumrawi, A.A., LeMay, V., Marshall, P.M., Hassani, B.T. and Robinson, D. 2005. Implementing a PrognosisBC regeneration sub-model for complex stands of Southeastern and Central British Columbia. Collaborative Project between Research Branch, Ministry of Forests. Forest Resources Management Department, University of British Columbia, and ESSA Technologies, Inc. Vancouver, British Columbia, Canada. Report submitted to: Forest 86 Science Program, Project No: Y051355, 157 p. http://www.for.gov.bc.ca/hre/pubs/azumrawi.htm accessed, June 2007.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A hybrid model to estimate natural recruitment and...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A hybrid model to estimate natural recruitment and growth in stands following mountain pine beetle disturbance Sattler, Derek Felix 2009
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | A hybrid model to estimate natural recruitment and growth in stands following mountain pine beetle disturbance |
Creator |
Sattler, Derek Felix |
Publisher | University of British Columbia |
Date Issued | 2009 |
Description | A method of linking SORTIE-ND and PrognosisBC was developed for the purpose of predicting natural regeneration and forecasting future stand conditions in mountain pine beetle (Dendroctonus ponderosae Hopkins - MPB) attacked stands in the Interior Douglas-fir (IDF) and Sub-Boreal Spruce (SBS) biogeoclimatic ecosystem zones of central and southeastern British Columbia. PrognosisBC, a spatially-implicit growth model, lacked a submodel suitable for predicting natural regeneration in unsalvaged MPB-disturbed stands. To fill this gap, estimates of regeneration (trees <7.5 cm diameter at breast height - DBH) were supplied to PrognosisBC using the light-mediated forest dynamics model SORTIE-ND and the linked model was used to forecast future stand conditions. In order to improve results, a density-dependent system of crown allometry equations to predict crown depth and crown radius was developed and then added to SORTIE-ND. The equations used stand-level measures of stems ha-¹, basal area (m² ha-¹), and the basal area of trees taller than the target tree to explicitly account of the effects of crowding on the crown axes. Additionally, crown radius and crown depth were used as dependent regressors. The equations were fit using a nonlinear three-stage least squares estimator and generally provided good estimates of crown depth and crown radius for lodgepole pine (Pinus contorta var. latifolia), hybrid spruce (Picea engelmannii x glauca (Moench) Voss), Douglas-fir (Pseudotsuga menziesii var. glauca (Beissn.) Franco) and trembling aspen (Populus tremuloides Michx.). Tests of the hybrid model with the improved system of crown allometry equations were performed using reconstructed plot data collected from natural stands disturbed by MPB 25-years ago. The hybrid model provided good estimates (small mean bias and low root mean square error) for the basal area of advance regeneration (2 < DBH < 7.5 cm) for lodgepole pine (Pinus contorta var. latifolia). The best estimates were achieved when trees <7.5 cm DBH were transferred from SORTIE-ND to PrognosisBC 15-years after MPB-disturbance. For trees <2 m in height, poor estimates of stems ha-¹ where obtained. Despite the shortcomings with respect to trees <2 m tall, the results suggest that linking SORTIE-ND and PrognosisBC is an effective method of building a hybrid model capable of being used in MPB-disturbed forests. However, full parameterization of the SORTIE-ND model is likely needed to obtain accurate estimates for all sizes of natural regeneration. |
Extent | 748943 bytes |
Subject |
Natural regeneration Mountain pine beetle Growth modelling Hybrid model |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-01-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0066900 |
URI | http://hdl.handle.net/2429/3819 |
Degree |
Master of Science - MSc |
Program |
Forestry |
Affiliation |
Forestry, Faculty of |
Degree Grantor | University of British Columbia |
GraduationDate | 2009-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2009_spring_sattler_derek.pdf [ 731.39kB ]
- Metadata
- JSON: 24-1.0066900.json
- JSON-LD: 24-1.0066900-ld.json
- RDF/XML (Pretty): 24-1.0066900-rdf.xml
- RDF/JSON: 24-1.0066900-rdf.json
- Turtle: 24-1.0066900-turtle.txt
- N-Triples: 24-1.0066900-rdf-ntriples.txt
- Original Record: 24-1.0066900-source.json
- Full Text
- 24-1.0066900-fulltext.txt
- Citation
- 24-1.0066900.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0066900/manifest