Development of regeneration standards for sustainable forest management by CRAIG FARNDEN BSF, University of Brtish Columbia, 1985 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Forestry) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April 2010 © Craig Farnden, 2010 ii Abstract Linking regeneration activities to desired future forest conditions is a critical component of sustainable forest management. In much of North America, regeneration standards are an important administrative tool for ensuring this goal, but the links to predicted future outcomes are often weak. There is a need for developing regeneration standards that are solidly based on defensible predictions of future stand attributes and are indicators of a wide range of forest values. This dissertation describes the application of an analysis framework wherein many variations of simulated juvenile stand stem maps are used as the basis for regeneration standards development. The simulated stem maps were generated at the scale of operational cutting units using newly developed software routines. Spatial variation was emphasized at the scale of 5 to 200 m, reflecting underlying terrain and ecosystem patterns. The simulated stands were projected into the future using existing growth models, and sampled using simulated regeneration surveys. Simplified models to mimic predictions of selected stand attributes from the growth models were developed using survey summary data as predictor variables. The analysis system was applied to test relationships between various measures of stocking and yield for conifer monocultures. For the predictive relationships, curve shape and dispersion were found to be highly variable based on various stocking measurement parameters. Stocking was found to imperfectly account for spatial variation in the regenerating stand, and stocking thresholds for achieving a specific relative yield were found to vary by species and the top height at which relative yield was assessed. A second analysis developed simplified models to evaluate the contribution of regenerated stands to landscape level species composition in a vertically stratified, spruce-aspen (Picea glauca and Populus tremuloides) forest type. Various field assessment methods and model formulations were contrasted. Site index was found to have a major impact on absolute yield and on the curve shape for relationships between stocking and species composition. This dissertation provides a greater understanding of how to develop, interpret and apply regeneration standards. Relating standards and assessment methods to forest-level management goals provides a key link that is often missing in stand-level regeneration assessments. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x Chapter 1. General Introduction and Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Stocking Standards and Stocking Surveys . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Linking Standards to Management Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Gaps in Knowledge and Practice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Research Goals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Methodological Approach. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Document Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 References. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 Chapter 2. Generating Stem Maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 Spatial Pattern Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 Filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Assigning Tree Attributes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 Implementation, Memory Limitations and Computational Speed . . . . . . . . . . . . . . . . . 27 Spatial Pattern Evaluation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 L(t) Case Studies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 Conclusions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 References. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 Chapter 3. Observations on the Relationship Between Regeneration Stocking and Yield. . . . . 39 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 Methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 Results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 iv Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 Conclusions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 References. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 Chapter 4. Simplified Models to Predict Future Mixedwood Outcomes from Regeneration Survey Summaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 Methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 Stand Generation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 Growth Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 Simulated Surveys . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 Fitting Simplified Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 Results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 Calibration of Sub-Models in the Stand Generator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 Simplified Models Based on Various Survey Methods and Fixed Site Index – Original Data Set. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 Simplified Models Based on Various Survey Methods and Fixed Site Index – Expanded Data Set. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Quadrat-Based Models Incorporating Spruce and Aspen Site Index . . . . . . . . . . . . . . . 76 Quadrat-Based Models For Discrete Levels of Spruce Site Index . . . . . . . . . . . . . . . . . 77 Comparison of Errors for Models Using Discrete Versus Continuous Site index . . . . . 79 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 Comparison of Survey Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 Options for Inclusion of Site Index. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 General Comments on Simplified models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 Conclusions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 Chapter 5. General Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 Generating Stem Maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 Relationships Between Stocking and Yield . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 Simplified Yield Models for Spruce-Aspen Mixedwoods . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 Additional Future Research Directions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 Use of Distance-Dependent versus Distance-Independent Growth models . . . . . . . . . . 97 Other Management Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 vClosing Remarks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 References. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Appendix A - Silviculture Survey User Guide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 Appendix B - A Silviculture Survey Methodology for Boreal Mixedwoods in Northeastern BC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 vi List of Tables Table 2-1. Model forms available for predicting multipliers (0≤y≤1) for adjusting tree heights relative to overtopping competition and for predicting diameters.. . . . . . . . . . . . . . . . . . . . 26 Table 3-1. Frequency of gaps, approximate affected area and associated random mortality outside of gaps for stands affected by a combination of clustered and random mortality. . . . . . . . 43 Table 3-2. Stand densities by terrain position for simulated stands intended to mimic natural regeneration patterns, along with area breakdown by terrain position. . . . . . . . . . . . . . . . . 44 Table 3-3. Fit statistics for equations predicting relative yield at each of 5 reference ages (z), along with maximum yield (Yz) at each age used to calculate relative yield.. . . . . . . . . . . 54 Table 4-1. Basic description of simplified models used to predict total yield at age 80 (Y80) and the proportion of that yield as conifer volume (Yc). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 Table 4-2. Coefficients and fit statistics for the models used to calculate Y80 by survey method with site index fixed at 20 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Table 4-3. Coefficients and fit statistics for the models used to calculate yield-based species composition (Yc) by survey method with site index fixed at 20 m. . . . . . . . . . . . . . . . . . . . 72 Table 4-4. Coefficients and fit statistics for the models used to calculate Y80 by survey method with site index fixed at 20 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 Table 4-5. Coefficients and fit statistics for the models used to calculate yield-based species composition (Yc) by survey method with site index fixed at 20 m. . . . . . . . . . . . . . . . . . . . 74 Table 4-6. Coefficients and fit statistics for the models used to calculate total yield at age 80 (Y80) and yield-based species composition (Yc) using summary data from a 10 m2 quadrat survey and both aspen and spruce site index. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 Table 4-7. Coefficients and fit statistics for the models used to calculate Y80 assuming discrete site index values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 Table 4-8. Coefficients and fit statistics for the models used to calculate yield-based species composition (Yc) assuming discrete site index values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 Table 4-9. RMSE values from the models incorporating site index, parsed by 1 m intervals of site index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 Table 4-10. Species composition (proportion Sw) at age 80 as simulated in MGM for varying levels of spruce and aspen site index. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 vii List of Figures Figure 1-1. Conceptual diagram of an analysis framework for developing simplified models to predict future stand attributes based on regeneration survey summary data. . . . . . . . . . . . . . 8 Figure 2-1. Stem maps of a stratified mix of trembling aspen overtopping white spruce. . . . . 18 Figure 2-2. Stem map of 15 m by 19 m plot derived from 70 mm aerial photos. . . . . . . . . . . . 19 Figure 2-3. Stem maps are created using either a random or lattice point process that are filtered to create desired stand structures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Figure 2-4. Voronoi polygons (a) are the basis for meso-scale landscape variation. . . . . . . . . . 21 Figure 2-5. Stem maps illustrating use of Voronoi polygons to define meso-scale features. . . . 22 Figure 2-6. Use of micro-scale patches for spatial diversity. . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Figure 2-7. Conifer height growth is affected by overtopping competition from a broadleaved tree patch. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 Figure 2-8. Cumulative beta distribution where alpha = 4 and beta = 1.5.. . . . . . . . . . . . . . . . . 27 Figure 2-9. The calculation of age at breast height (1.3 m) from total age and height assumes a linear approximation of the height growth curve for heights above 1.3 m. . . . . . . . . . . . . . 28 Figure 2-10. Examples of L(t) functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 Figure 2-11. Cumulative distribution for spacing factor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 Figure 2-12. Spatial distributions as characterized by the L(t) function for the 4 large stem maps illustrated in Figure 2-1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 Figure 2-13. Selected stem maps (scale in metres) and corresponding L(t) functions for three selected plots in 8 to 12 year-old aspen stands, illustrating the range in spatial variation observed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Figure 2-14. Stem map and associated L(t) functions evaluated in four localized sampling windows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 Figure 3-1. Examples of stem maps showing three primary variations in spatial pattern . . . . . 43 Figure 3-2. Yield at 80 years for western redcedar, site index 28 m, as predicted by percent stocking values for six different plot sizes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 Figure 3-3. Yield at 80 years for western redcedar, site index 28 m, as predicted by well-spaced values for six combinations of MITD and m-value . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 viii Figure 3-4. Yield at 80 years for lodgepole pine, site index 20 m, as predicted by percent stocking values for six different plot sizes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 Figure 3-5. Yield at 80 years for lodgepole pine, site index 20 m, as predicted by well-spaced values for six combinations of MITD and m-value . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 Figure 3-6. Conceptual diagram of growing space utilization and its measurement in stocking surveys . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 Figure 3-7. The sampling variance for measures of lodgepole pine stocking (a) appears to be a good indicator of differences in spatial pattern, and may be useful in accounting for some of the variation in relationships predicting yield from regneration survey statistics . . . . . . . . 51 Figure 3-8. The relationship between stocking and relative yield varies by age (top height) of assessment and species . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 Figure 4-1. Two-parameter cumulative Weibull fit of relative spruce height (adjustment factor) as predicted by spacing factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 Figure 4-2. Linear fit of diameter as predicted by height for aspen and spruce . . . . . . . . . . . . . 70 Figure 4-3. Fit plots for simplified models using summary statistics from10 m2 quadrat surveys, 5-5-2-2 surveys, and mixedwood well-spaced surveys as parameters to predict total volume and species composition (% spruce volume). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 Figure 4-4. Fit plot of linear model for species composition using untransformed variables from the mixedwood well-spaced survey . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 Figure 4-5. Fit plots for simplified models using the expanded data set and summary statistics from10 m2 quadrat surveys, 5-5-2-2 surveys, and mixedwood well-spaced surveys as parameters to predict total volume and species composition . . . . . . . . . . . . . . . . . . . . . . . . 75 Figure 4-6. Fit plots for simplified models for (a) Y80 and (b) Yc that include site index of both spruce and aspen as predictive parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 Figure 4-7. Comparison of errors (RMSE) between models fit for discrete levels of site index (from Tables 4-7 and 4-8) and those parsed from the fit of a model that includes site index as a predictive variable (Table 4-9). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 Figure 4-8. Mean change in total volume for each 0.25 m increment in spruce site index, based on MGM predictions used in development of the discrete site index models. . . . . . . . . . . . 80 Figure 4-9. Selection of points (bold red) by site index range in the fit plot for the Yc model incorporating site index as continuous variables (Eq. 16) . . . . . . . . . . . . . . . . . . . . . . . . . . 82 ix Figure 4-10. Comparison of fit plots for species composition models using (a) both SISw and SIAt as continuous variables versus where I2 = 0.928 and RMSE = 0.0710 versus b) only SISw. as a continuous variable with SIAt inferred, where I2 0.845 = and RMSE = 0.1077 82 Figure 4-11. Survey assessments that employ competition thresholds for tallying trees will always ignore a component of the stand . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 Figure 4-12. Linear fit of Yc using summary data from the 5-5-2-2 survey summaries . . . . . . 85 xAcknowledgements My journey through the last 4 1/2 years has been one of the most stimulating and exciting periods of my life. Many people have contributed to this experience, and I am grateful to all of them. My supervisor, Dr. Bruce Larson has been a tremendous teacher, friend and fountain of knowledge. His incredible passion for all things related to forestry is an inspiration to all that know him. Our many hours of conversation will long be fond memories. My committee members, Dr. Peter Marshall and Dr. John Nelson, listened to and supported my ideas, provided invaluable advice and guidance, and always coaxed me onward. Over many years, Dr. Gordon Weetman encouraged me to return to UBC for graduate studies. I finally listened. Gordon’s open door and sage counseling are greatly appreciated. Dr. Val LeMay was a great help for all things statistical and was always a friendly ear. An early (and continuing) mentor in my career was Ian Moss – he taught me how to “think like a tree”, the importance of spatial scale, and what forest management is really all about. The collegiality and friendship of many graduate students and staff members around the Forest Science Centre was a great source of inspiration and support. Of particular note are Mariano Amoroso, David Jack, Anna Tikina, Andrew Innerd, Leah Rathburn, Babita Bains, Candace Parsons and Rasmus Astrup. Funding for my research was provided by the BC Government’s Forest Investment Account through a delivery agreement with Canadian Forest Products Ltd., Fort St. John. Jeff Beale (then) at Canfor immediately recognized the value of this research and was an enthusiastic proponent. Wes Neumeier at Canfor shepherded the process through. Other professional colleagues whose ideas notably influenced this work include Pat Martin, Eleanor McWilliams, Al Powelson, John McClarnon, Richard Kabzems and Wendy Bergerud. My parents, David and Audrey Farnden, always encouraged me to find my own path and to follow it with faithfulness and enthusiasm. My wife, Heather, and kids, Cameron and Meaghan, have been loving, patient and true supporters. 1Chapter 1. General Introduction and Background1 Forest management is the art and science of balancing the flow of extractive goods and services with maintaining the forest as a whole in a desirable state. Every management action that foresters undertake should ultimately be linked to a vision of this balance –the forest management objectives. As part of the management process, we must be aware of how individual activities impact on future stand structures, and how local structures are summed at the landscape unit and broader scales (Baskerville 1992, O’Hara et al. 1994, Martin 2005). Assessing the contribution of individual stands to the larger forest condition is often fraught with uncertainty, particularly at an early age. Over the time span from regeneration to maturity, forest stands undergo massive change. Given this uncertainty and the fact that the largest opportunity to shape future forest conditions occurs at the harvest and reforestation stages, it is critical to know the conditions to create in young stands such that future forest landscape objectives are likely to be achieved. A key element of many forest management systems in North America and elsewhere is the application of regeneration standards (Brand and Weetman 1986, Armson 2005). These administrative rules outline minimum expectations for a regenerating crop. Most commonly specified in the standards are species preferences and minimum establishment densities, with factors such as tree health, size, vigour and freedom from competition by non-crop vegetation often included. The standards have historically had a very strong timber production emphasis, although there appears to be a slow trend developing to also incorporate other resource values. Regeneration standards are a major focus of this dissertation, or more specifically, the development of appropriate standards for a given management situation. Motivations for the study originate in deliberations undertaken over at least two decades in which I was frequently involved as a consulting forester. Appropriate regeneration standards were being sought for boreal mixedwood stands – mixtures of white spruce (Picea glauca (Moench) Voss) and trembling aspen (Populus tremuloides Michx.) – in northeastern British Columbia (BC). The typical approach to such standards in BC, based on a conifer model whereby all species being regenerated could largely co-exist in a single canopy stand, was to specify minimum acceptability thresholds for stocking at the stand level. Variations in standards were specified by climate (biogeoclimatic subzones) and bio-physical site types (site series), in particular those relating to species acceptability and stocking. This system worked well for ensuring the conifer stand types for which it was originally designed became suitably occupied by regenerating trees, but for mixedwood stands it had a critical shortcoming that only slowly became apparent. 1 Some portions of this chapter have been previously published. All subsections of the literature review and the subsection on methodological approach were extracted with minor revisions from Farnden (2009). 2A key element of the conifer standards was the concept of “free-growing”, whereby growing space allocation was assessed at the level of individual trees. For a regenerating stand to be acceptable, a minimum number of conifer crop trees had to be assessed as free of competing vegetation, including that from broadleaved trees. Attempting to translate this to the mixedwood context required a threshold number of uniformly distributed aspen overstorey trees such that growth of the more valuable spruce understorey would be largely unimpeded. This one-size-fits-all approach from the conifer standards was initially proposed but never implemented. Operational foresters quickly realized that attempting to create the spatially uniform mixtures suggested by the standards was grossly impractical. Licensees instead chose an approach commonly referred to as “unmixing the mix”; mixed stands of aspen and spruce were regenerated to pure stands of one or the other at varying scales of patch size. Both government and licensee foresters recognized that this situation was less than desirable. However, the absence of appropriate administrative rules and tools for setting goals and assessing regeneration success for intimate spruce-aspen mixtures left few alternatives other than complete avoidance of these stand types. A potential solution to this problem occurred to me as the result of several modeling exercises I undertook that were aimed at simulating proposed operational strategies for regenerating mixedwood stands. These exercises made it apparent that (i) the concept of a free-growing conifer in the mixedwood context was inappropriate and that (ii) allowing for considerable variation in spatial pattern of the mixture was extremely important in understanding how to quantify regeneration success. From discussions I had with both foresters and biologists working in the region, it was apparent that variation in spatial pattern was both desirable and far more practical to achieve than a uniform admixture. The conifer-centric concept of a free- growing tree had to be abandoned, and a new system was required that embraced quantification of growing space allocation at the stand rather than the tree level. Literature Review Stocking Standards and Stocking Surveys The origin of North American stocking standards is found in the early part of the 20th century based on concerns that extensive exploitive harvesting was not followed by sufficient natural regeneration to sustain forest conditions over the long term (Armson 2005). There are many variations on stocking standards combined with many survey methodologies, often leading to confusion when different systems employ similar terminology or when similar systems use different measurement standards (Hosie 1953, Armson 2005). 3The earliest stocking standards in North America were based on population estimates expressed as a number of trees per unit area or density (Lowdermilk 1927, Haig 1931, Hosie 1953), and such standards persist at least to a small extent today. Allen et al. (1951) and Kaltenberg (1980) described the list quadrat method, in which all of the trees present in a plot are tallied to calculate values of density. In some variations on this theme, the number of trees tallied is capped either at a maximum value (Kaltenberg 1980) or using spacing rules whereby counted trees must be a minimum distance apart (Allen et al. 1951). The latter variation is intended to deal with the distribution problem, whereby a simple tally count will de-emphasize the presence of unstocked voids that cause reductions in timber production (Haig 1931). The system predominantly being used in British Columbia is an extension of this philosophy, whereby a stratified count of total trees is tallied along with a subset of “well-spaced” trees defined based on crop tree desirability criteria and minimum inter-tree distances (Wyeth 1984). In the early 1920’s, an alternate system was developed both to account for the distribution problem and to add sampling efficiency (Lowdermilk 1927). The solution was a measure of site occupancy expressed as the percentage of stocked quadrats in a sampled population, a quadrat being evaluated as either containing (stocked) or not containing (unstocked) an acceptable crop tree. Quadrat-based stocking standards have been widely applied and continue in use today for many jurisdictions. However, over the ensuing eight decades since their inception, there has been considerable tinkering and debate over appropriate application. Two inter-related themes have dominated this evolution, focusing on plot size and thresholds for acceptable versus unacceptable stocking. Plot size and thresholds are mathematically linked (Munger 1945, Grant 1951, Ker 1953), with larger plots resulting in higher percentages of stocked plots. Some authors (e.g. Haig 1931) have speculated that plot size should be dictated by desired tree spacing, whereby a stand with theoretically perfect spacing would have exactly one tree in every plot. Others have noted the approximate scalability of plot size versus percent stocking (Allen et al. 1951, Grant 1951), with percent stocking from small plots more easily scaled to the equivalent for large plots (Grant 1951) rather than the reverse. However, Ghent (1969) warned that all scaling risks introduction of bias, particularly when the scale aadjustment is large. Allen et al. (1951), Grant (1951) and Corona et al. (1998) discussed and favoured plot size decisions based on cost and/or statistical efficiency of field sampling. Thresholds for acceptable stocking have also varied considerably even where plot size is held constant. For example, Allen et al. (1951) noted remarkable differences between standards for British Columbia and those used in similar forest types in the Pacific Northwest (Washington 4and Oregon). Hosie (1953) found similar discrepancies among different management agencies in Ontario. Standards other than those based on stocked quadrats have also been employed. Included here are schemes based on point-to-tree distances. Advantages of such a system include improved statistical efficiency (MacLeod and Chaudrhy 1979) and quantification of spatial aggregation (Dennis 1984). Extensive use of such systems has not been observed, possibly due to more complex and slightly more costly field procedures (MacLeod and Chaudrhy 1979), less transparent summaries and the operational inertia of the stocked quadrat system. Linking Standards to Management Objectives The primary purpose of stocking standards is to ensure that adequately stocked stands of trees have replaced stands removed or destroyed by disturbances. However, the definition of what is adequately stocked has rarely been clear (Allen et al. 1951, Hosie 1953, Armson 2005). Many standards based on quadrat surveys are class based, with percent stocking thresholds defining subjective classes such as “fully stocked”, “well stocked”, “moderately stocked” and “understocked” (Armson 2005). Hosie (1953, pg. 20) decried the arbitrary nature of most thresholds, and expressed doubt that the impact of varying stocking percentages for a wide range of forest conditions could be determined without a “… great effort in research over a very long period of time”. More and more, there are demands to make explicit linkages between reforestation standards and management objectives (Martin 2005), most often expressed in terms of future yields. In a notable early example, Staebler (1949) linked stocking percentages to the degree of normality in yield tables for Douglas-fir. Similarly, Braathe (1982) described relationships between percent unstocked areas and yield for Norway spruce, Stage and Ferguson (1984) presented methods for linking regeneration survey data to the Prognosis growth model, Newton (1998) demonstrated a method of linking regeneration standards to yield for black spruce using point estimates of density and a stand density management diagram, and Martin et al. (2002) used survey summary data to directly predict timber yields. In Alberta, an initiative was begun in 2001 to explicitly link regeneration standards to yield curves (Alberta Reforestation Standards Science Council 2001), a process that is not yet complete in 2010. In BC, target stocking levels based on the well-spaced tree concept are intended to achieve maximum volumes at minimum stand densities (Wyeth 1984). Studies by Bergerud (2002), Martin et al. (2005) and others have confirmed reasonably good relationships between numbers of well-spaced trees and yield at or near culmination of mean annual increment, although inconsistent application of minimum inter-tree distances have diminished the potential value of this relationship. Many other examples of links from 5regeneration standards to yields likely exist as locally applied management reports, but have not been noted either in peer reviewed journals or broadly distributed professional literature. Under the paradigm of sustainable forest management, foresters are becoming increasingly aware that simply achieving fully stocked stands after harvest is no longer sufficient to achieve management success (O’Hara et al. 1994). To this end, a new set of regeneration standards is required that go beyond simply stocking, and that are explicitly linked to long-term forest management objectives (Martin 2005). Examples of stand structures or elements that may be desirable under certain circumstances and that extend beyond just stocking include: (i) the need for gaps in a forest canopy, with specific gap sizes and spatial distributions, (ii) the need for certain species compositions and patterns of inter-mixing, (iii) the need for certain vertical structures such as multiple cohorts and/or multiple layers, and (iv) the need for discrete canopy elements such as snags, or crowns of irregular shape or coarse (wolf tree) structure. A few cases of standards with a perspective beyond simply timber production have been developed in BC by adapting existing stocking standards. An example is a set of guidelines for managing lowland grizzly bear forage habitat (Wood 2001). Martin (2005) outlined a new way of thinking about goal-oriented regeneration standards using boreal mixedwoods as an example. He emphasized that the desired future forest condition should be the context within which regeneration standards are developed, with a focus expanded beyond just stocking to include distribution of patch types by species composition. Based on the principles of performance measurement, it was suggested that such a system would require: (i) an explicit statement of the goals of forest management, (ii) identification of the critical factors required to achieve the goals, and (iii) a monitoring system to measure adherence to thresholds based on the critical factors. Gaps in Knowledge and Practice The forest regeneration management systems envisioned by Baskerville (1992) and Martin (2005) require expansion of the perspective for reforestation objectives to the landscape or broader scales using ledger or GIS systems to track cumulative accomplishments. This, in turn, will allow reforestation standards to become active tools for forest management rather than administrative measures of contractual or regulatory compliance. In order to be useful components of a forest management system, the silviculture survey summary statistics upon which standards are based must have predictive ability for desired future outcomes. A set of simple models is required that uses the survey data and/or summary statistics as inputs, and makes predictions about forest attributes or values that are critical indicators of management success. Such a modeling framework should have the following attributes: 61. Reliability: model outcomes should be reasonably accurate under the range of conditions that might be reasonably expected to occur during the normal course of day-to-day operations. 2. Consistency: models should be relatively simple and easy to implement such that several individuals collecting data for the same sample frame for the purpose of making the same predictions would obtain similar outcomes with relative ease. To this end, it is important that models be free of predictive parameters that are difficult to measure and to which model outcomes are particularly sensitive. 3. Transparency: models should not be “black boxes”; model structures should be such that a typical manager can readily understand how and why variations in the model inputs will result in changes in the predicted outcomes. 4. Practicality: the process of collecting data and running models must be simple and inexpensive. In many cases, stand level models already exist that at least partially satisfy the first three of these requirements. However, application of these models is frequently inhibited for day-to-day operations by data requirements for initializing simulations, time requirements for data collection and setting up simulations, and the lack of specialized knowledge required to properly run the models. Major improvements could be achieved either through stronger linkages between routine monitoring programs (such as regeneration surveys) and the data input requirements for existing predictive models, or the construction of new predictive models that are based on existing monitoring summaries for initial model conditions. There are currently two provincial level initiatives in western Canada to address this issue. In Alberta, the stand level GYPSY model (Huang et al. 2009) is being revised to use stocked quadrat summary data as inputs, thus providing a direct link between estimates of stocking to predictions of yield. In British Columbia a software program called TIPSY (BC Ministry of Forests 2009) is commonly used to make yield forecasts for regenerated stands. In this case, yields are interpolated from a set of variable density yield tables, each of which was developed by simulating stand development on a set of reference stem maps using an individual tree, spatially explicit growth model. To date, total stand density and one of up to three spatial pattern classes have been the primary inputs for TIPSY related to stocking. A planned update for 2010 (pers. comm. W. Bergerud, BC Ministry of Forests, December 2009) will provide the option of instead using actual estimates of stocking from regeneration surveys. However, there are many jurisdictions where such modeling efforts are incomplete or missing. In BC, for example, the TIPSY model simulates species mixtures by simply 7weighting the yields from pure species stands based on initial composition by number of stems (BC Ministry of Forests 2007). This makes the model unsuitable for stands where species composition at harvest is significantly affected by successional patterns, as is the case in most mixtures of conifers and hardwoods in BC. In many other jurisdictions, there are no models at all that explicitly link stocking to yield or to other management objectives. Research Goals The intent of this research is to explore relationships between the concept of stocking as commonly applied in forest regeneration surveys and the achievement of long-term management objectives. The current study limits its focus on management objectives to yield and species composition (with a particular focus on spruce-aspen mixtures in the boreal forests of northeastern BC), recognizing that in many cases there may be a much broader set of objectives that need to be addressed. Methodological Approach An analysis framework was conceived to facilitate development of simple and reliable predictive relationships between silviculture survey outcomes and long-term achievement of management objectives. The framework is composed of the following steps (see also Fig. 1-1): 1. Generate a wide range of simulated juvenile stand conditions (stem maps) that is representative of the outcomes of a reasonable set of likely reforestation plans. Stem maps should be generated at the scale of operational cutblocks (i.e. 5 to 30 ha) and should encompass tree size, species and spatial variability induced by factors such as natural regeneration patterns, terrain position and ecosystems, and treatment history (e.g. logging pattern). 2. Predict future stand attributes (yield, species composition, vertical or horizontal structure) using currently or soon-to-be available growth models (e.g. MGM2, SortieND3, TASS III4). Some spatially explicit models such as SortieND and TASS are able to grow 2 Mixedwood Growth Model: an individual tree, distance-independent model for boreal mixed species stands developed at the University of Alberta (University of Alberta 2007) 3 SortieND: an individual tree, distance-dependent, resource mediated model of tree growth developed in the eastern United States (Pacala et al. 1993) and later adapted for use in British Columbia (Coates et al. 2003, Astrup 2006). 4 TASS III: an individual tree, distance-dependent model for single species, even-aged stands originally developed by Mitchell (1975). The model has since undergone extensive development by the BC Ministry of Forests Research Branch, with a major revision incorporating multiple species, multiple canopy layers and a light sub-model due for release in 2010. 8Figure 1-1. Conceptual diagram of an analysis framework for developing simplified models to predict future stand attributes based on regeneration survey summary data. A 36 ha virtual stand (top left) can be sampled with silviculture survey plots (red dots) to generate stand level survey metrics. Future attributes such as yield and/or species composition for the same stand can be predicted using existing growth models (bottom left). Forecasted outcomes can then be contrasted with stand level objectives (in this case represented by the proportion of merchantable yield as conifer volume). If there is a strong correlation between the survey metric and predicted levels of the objective function (bottom right), the survey metric is useful for determining whether a particular juvenile stand is likely to produce a desired outcome, and desirable thresholds can be explored for use as regeneration standards. If the correlation is weak (top right), alternate metrics should be investigated. 0 100 200 300 400 500 600 0 100 200 300 400 500 600 36 ha Mixedwood Stand Distance (m) D is ta nc e (m ) 0 20 40 60 80 100 0 20 40 60 80 100 Distance (m) D is ta nc e (m ) 0 20 40 60 80 100 0 20 40 60 80 100 Distance (m) D is ta nc e (m ) 70 75 80 85 90 45 50 55 60 65 Distance (m) D is ta n ce ( m ) 0 5 10 15 20 0 5 10 15 20 Distance (m) D is ta n ce ( m ) 0 100 200 300 400 500 0 50 100 150 Stand Age (yrs) Merch. Vol. (m3/ha) Spruce Aspen Total M er ch Y ie ld ( m 3 / ha ) 0 0.2 0.4 0.6 0.8 1.0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Survey Metric C on ife r P ro po rt io n of M er ch . Y ie ld 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 0.2 0.4 0.6 0.8 1.0 Survey Metric C on ife r P ro po rt io n of M er ch . Y ie ld ��������� ������� ������� 9entire stands based on stem maps, but distance independent models such as MGM must work from sampled treelists. In the latter case (as illustrated in Fig 1-1), many random samples from the simulated stand are used to initiate discrete projections of future stand conditions, with the final outcomes compiled from the composite statistics. 3. Contrast projected stand conditions with desired conditions. Assess each stand for success/failure or place into classes reflecting a range of desired conditions. 4. Correlate reforestation metrics to measures of achievement of management objectives to provide a simplified model for predictive purposes. The framework described here has close similarities to processes formerly applied successfully for simply structured conifer stands by Bergerud (2002) and Martin et al. (2002). Ultimately, it is intended to form part of a larger program of adaptive management, whereby management practices are continuously vetted and improved as better tools and techniques are found, and ongoing monitoring of outcomes indicates the degree to which initial predictions are accurate. Objectives Based on the conceptual analysis framework described above, the following objectives have been identified for this dissertation: 1. to construct a stand generation system that will simulate spatially appropriate stem maps of juvenile forest stands and facilitate analyses such as simulated regeneration surveys and growth modeling, and 2. to demonstrate the utility of such a system for developing simple predictive models that are sufficiently reliable and practical for use in developing and implementing regeneration standards. Document Structure The first objective listed above is primarily addressed in Chapter 2, wherein a model is presented that is intended to generate realistic stem maps of juvenile forest stands at the scale of operational cutting units. Realism is introduced by the user through application of its inherent flexibility to include many combinations of stand cohorts and spatial patterns within an underlying biophysical pattern of terrain units. Realism is tested through application of spatial statistics and comparison to known values. The stand generation routines are included within a separately described software package called the Silviculture Survey Simulator, which includes routines for simulating surveys and exporting stem maps as initial conditions for growth simulations. A user’s guide to the software has been included as Appendix A. 10 The second objective is addressed in Chapters 3 and 4 with case studies of model development to predict future outcomes based on regeneration survey statistics. Chapter 3 investigates the simple case of single species conifer stands, and explores the relationships between various measures of stocking and volumetric yield as a simple expression of management objectives. A large emphasis is placed on the application of various stocking metrics and their value as predictive parameters. Chapter 4 moves to the more complex case investigating stratified mixtures of two species. This case addresses the operational needs of forest licensees working in boreal mixed spruce- aspen stands in northeastern British Columbia, where commitments have been made to manage landscape scale species composition within the context of a sustainable forest management plan. The predictive models developed in this study allow prediction of stand level yields by species, which in turn can be tracked in a separately developed ledger system to monitor the cumulative contribution of regenerating stands to the landscape level targets. In Chapter 5, I provide a summary of the major outcomes of the three previous chapters and a discussion of their inter-relationships and implications for forest management. Also included is a section on new research directions that will explore concepts that were problematic within this dissertation or are expected to provide further insights into how regeneration standards can better help foresters to achieve their management objectives. 11 References Alberta Reforestation Standards Science Council. 2001. Linking regeneration standards to growth and yield and forest management objectives. Alberta Ministry of Sustainable Resource Development, Edmonton. 52p. Allen, G.S., Griffith, B.G., and Ker, J.W. 1951. A comparison of several regeneration survey methods in terms of cost and usefulness. B. C. Lumberman 35: 66, 68, 89-90, 92, 94, 96. Armson, K.A. 2005. Regeneration standards: what has the past to show us? For. Chron. 81: 781- 784. Astrup, R. 2006. Modeling growth of understory aspen and spruce in western boreal Canada. PhD Thesis, University of British Columbia, Vancouver. Baskervile, G.L. 1992. Forest analysis: Linking the stand and forest levels. In: The Ecology and Silviculture of Mixed-Species Forests: A Festschrift for David M. Smith. Kluwer Academic Pub., Dordrecht, Netherlands. pp. 257-277. BC Ministry of Forests. 2007. TIPSY 4.1 help – prorating multiple species. Help files embedded in TIPSY version 4.1 software. BC Ministry of Forest and Ramsoft Systems Ltd., Victoria BC. BC Ministry of Forests. 2009. TIPSY features/functions. URL: http://www.for.gov.bc.ca/HRE/ gymodels/TIPSY/features.htm. Accessed Jan. 16, 2010. Bergerud, W.A. 2002. The effect of the silviculture survey parameters on the free-growing decision probabilities and projected volume at rotation. Land Man. Hndbk No. 50. BC Ministry of Forests, Victoria. 32p. Braathe, P. 1982. Stocking control and its effect on yield. Forest Industry Lecture Series, University of Alberta, Edmonton. 21p. Brand, D.G., and Weetman, G.F. 1986. Standards for regeneration establishment in Canada: a case study for Douglas-fir. For. Chron. 62: 84-90. 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. Corona, P., Leone, V., and Saracino, A. 1998. Plot size and shape for the early assessment of post-fire regeneration in Aleppo pine stands. New For. 16: 213-220. 12 Dennis, B. 1984. Distance methods for evaluating forest regeneration. In: New forests for a changing world : proceedings of the 1983 Convention of the Society of American Foresters, October 16-20, 1983, Portland, Oregon. Society of American Foresters, Bethesda, Md. pp. 123-128. Farnden, C. 2009. An analysis framework for linking regeneration standards to desired future forest conditions. For. Chron. 85: 285-292. Ghent, A.W. 1969. Studies of regeneration in forest stands devastated by the spruce budworm. IV. Problems of stocked-quadrat sampling. For. Sci. 15: 417-429. Grant, J.A.C. 1951. The relationship between stocking and size of quadrat. Univ. Toronto For. Bull. No. 1. 35p. Haig, I.T. 1931. The stocked quadrat method of sampling reproduction stands. J. For. 29: 747- 749. Hosie, R.C. 1953. Forest regeneration in Ontario based on a review of surveys conducted in the province during the period 1918-1951. Univ. Toronto For. Bull. No. 2. 134p. Huang, S., Meng, S.X. and Yang, Y. 2009. A growth and yield projection system (GYPSY) for natural and post-harvest stands in Alberta. Alberta Sustainable Resource Development, Edmonton AB. 22p. Kaltenberg, M.C. 1980. Evaluation of regeneration sampling methods: a Monte Carlo analysis using simulated stands. DNR Report No. 39. State of Washington, Dept. of Nat. Res. 50p. Ker, J.W. 1953. The relationship between the number of trees [seedlings] per acre and the percentage stocking of reproduction. J. For. 51: 342-344. Lowdermilk, W.C. 1927. A method for rapid surveys of vegetation. J. For. 25: 181-185. MacLeod, D.A., and Chaudhry, M.A. 1979. A field comparison of distance and plot methods for regeneration surveys. For. Chron. 55: 57-61. Martin, P.J. 2005. Design of regeneration standards to sustain boreal mixedwoods in western Canada. Int. For. Rev. 7(2): 135-146. Martin, P.J., Browne-Clayton, S. and McWilliams, E. 2002. A results-based system for regulating reforestation obligations. For. Chron. 78(4): 492-498. Mitchell, K.J. 1975. Dynamics and simulated yield of Douglas-Fir. For. Sci. Mon. No. 17. 39p. Munger, T.T. 1945. Stocked quadrats vs. number of trees as a basis for classifying reforesting land. For. Res. Notes No. 33. USDA For. Serv., Pacific NW For. Exp. Stn. 7p. 13 Newton, P.F. 1998. An integrated approach to deriving site-specific black spruce regeneration standards by management objective. For. Ecol. Manage. 102: 143-156. O’Hara, K.L., Seymour, R.S., Tesch, S.D., and Guldin, J.M. 1994. Silviculture and our changing profession: leadership for shifting paradigms. J. For. 92: 8-13. Pacala, S.W., Canham, C.D., and Silander, J.A.,Jr. 1993. Forest models defined by field measurements: I. The design of a northeastern forest simulator. Can. J. For. Res. 23: 1980-1988. Staebler, G.R. 1949. Predicting the volume and normality of reproduction stands of Douglas-fir. J. For. 47: 828-833. Stage, A.R., and Ferguson, D.E. 1984. Linking regeneration surveys to future yields. In: New forests for a changing world : proceedings of the 1983 Convention of the Society of American Foresters, October 16-20, 1983, Portland, Oregon. Society of American Foresters. pp. 153-157. University of Alberta. 2010. Mixedwood Growth Model (MGM). 2008. URL: http:// www.ales2.ualberta.ca/rr/mgm/mgm.htm. Accessed Jan. 13 2010. Wood, C. 2001. Grizzly bear habitat in managed forests: silviculture treatments to meet habitat and timber objectives. Ext. Note 54. BC Min. For., Res. Br., Victoria BC. 7p. Wyeth, M.H. 1984. British Columbia Ministry of Forests regeneration survey system. In: New forests for a changing world : proceedings of the 1983 Convention of the Society of American Foresters, October 16-20, 1983, Portland, Oregon. Society of American Foresters, Bethesda, Md. pp. 149-152. 14 Chapter 2. Generating Stem Maps5 Introduction The artificial generation of tree sample or census data has been a topic in the forestry literature since computers became widely available in forestry laboratories. Some of the earlier applications focused on studies of sampling schemes (e.g. MacLeod 1977) and exploration of indices or descriptors of randomness (e.g. Stauffer 1977). In recent years, there has been considerable attention to this topic because mapped tree data are required to initiate many types of simulations using complex forest growth models. For distance-independent (or position-independent) models, sources of real data with individual tree records are relatively plentiful and include research plots, inventory plots, and timber cruises. Where only generalized stand-level descriptions are available, or critical attributes are missing, various imputation methods have been developed to derive tree level data, such as those described by Ek et al. (1997), Temesgen et al. (2003) and Hassani et al. (2004). However, interest is growing in distance-dependant individual tree growth models to explore stand development patterns where the spatial arrangement of trees is an important factor in simulated outcomes (e.g. Pommerening and Stoyan 2008, Puettman et al. 2008). Such models require stand descriptions which include spatial coordinates of each tree, data which are costly to collect and rarely available (Pretzsch 1997). A critical step in the creation of simulated stands is the generation of an appropriate point pattern to represent the spatial position of trees. Spatial patterns can be mathematically produced using a set of methods collectively known as point processes, and their application in forestry to generate tree maps has be summarized by Hanus et al. (1998), and Stoyan and Penttinen (2000). Commonly applied point processes include: • the Poisson process, which creates a completely random pattern • the Poisson cluster process, which creates random locations within random clusters, • the doubly stochastic Poisson or Cox process which randomly varies the densities in Poisson clusters, • the lattice process which creates a regulated pattern based on a grid, and • the hard-core process, which is used in combination with other processes to limit the proximity of tree pairs, usually to mimic the spatial inhibition caused by self-thinning. 5 A version of this chapter will be submitted for publication. Farnden, C. 2010. A large-scale stem map generator for juvenile forest stands. 15 Adaptations and combinations of processes have been added by many authors to better match the variability observed in natural complex systems. For example, a stand structure generation model called STRUGEN was first described by Pretzsch (1997) and subsequently incorporated into a software suite associated with the SILVA growth model (Pretzsch et al. 2002). The intent of this model is to generate “… a given number of (trees) with known stem diameters and height distributions … arranged on a certain stand in such a manner that the greatest possible conformity be achieved between actual and generated structure” (Pretzsch (1997, pg. 243). Points representing trees were determined using a nonhomogeneous Poisson process, where potential points were first created using a homogenous Poisson process then screened to achieve a desired spatial pattern. Potential points were eliminated in a first step based on position- dependent probabilities (to create gaps and clumps), and in a second step to eliminate trees that did not satisfy minimum inter-tree distance requirements (a hard-core process). These steps were repeated for additional species, with controls for the degree of intermingling. A further innovation was the addition of marked points (Penttinen et al. 1992). In this approach, values such as tree species or size are associated with each point and affect the generation of spatial pattern. Marked point processes have been present in statistical literature for some time, but have only recently been introduced to the forestry literature. For generating tree maps, these methods are often linked with a Markov Chain Monte Carlo (MCMC) process, in which an arbitrary starting distribution of points is modified by randomly moving individual points and testing to find if the result better satisfies the objective function than the previous condition. A recently emerging trend in the use of marked points is the application of Gibbs process modeling, which offers considerable potential for mapping both non-stationary conditions and spatial variation based on correlations of point marks (MØller and Waagepetersen 2006, Illian et al. 2008). Examples provided by Kokkila et al. (2002) demonstrate the simulation of linear features such as roads, variations in seedbeds as introduced by planting furrows and the presence of stumps, proximity to seed trees, and non-stationarity related to biophysical (ecosystem) gradients. The Gibbs process defines a repulsive field potential for each point based on the mark. Between any two trees, the repulsive force increases with increasing tree size, but decreases with increasing distance. Successive moves are made in an MCMC process, until an acceptably low potential is achieved for the entire system. The stand level attribute of non-stationarity is introduced by adding a mapped background field. This feature is used to adjust pair-wise field potentials, either increasing or decreasing the effect possibly to the point where it becomes an attractive force. 16 A completely different method of generating realistic stem maps was described by Pommerening and Stoyan (2008). Their approach extrapolated spatial structures observed in a set of small circular sample plots (or subwindows) to the entire intervening area using simulated annealing optimization techniques with nearest neighbour summary statistics (NNSS) as the objective function. Sample plot data representing approximately 25% of a given stand were used to determine both the NNSS values (such as distances to neighbors, directional index, species mingling, dbh differentiation and dbh dominance), the expected point intensity for the stand as a whole, and the distribution of mark values to be associated with the simulated points. Within the reconstructed stand, the sample plots were fixed in place, along with a comparable intensity of points with associated marks distributed over the remaining area using either a Poisson or hard-core process. Similar to MCMC algorithms, a sequence of iterative steps was taken in which a single point (or tree) was moved, and the stand condition re-evaluated against the NNSS objective function. Unlike MCMC, the moves at each iteration were increasingly restricted to improve the efficiency of estimating the global optimum solution. In addition to generating spatial patterns, there is also a need for realistic tree attributes such as species, height, diameter and crown dimensions. These are assigned based on summary statistics from real data or, less desirably, from theoretical relationships. In most cases, attributes are assigned in a specific order, starting with a class variable such as species and/or a size variable such as height or diameter. In unmarked point processes, the generation of the spatial pattern is independent of the tree attributes. In marked point processes, one or both of these attributes will be assigned prior to generating the point pattern based on known probability distributions, and will be used in the generation of the spatial pattern. Commonly, additional tree parameters, where required, are based on the initial size attribute combined with the competitive environment as represented by the proximity of surrounding points. Pretzsch (1997) assigned height, crown width and crown base attributes using allometric relationships to diameter by species. Similarly, many individual tree growth models use allometric relationships to generate missing attributes for trees in the initial condition. The need for systems to generate and model complex forest stand conditions is being driven by an increasing desire to manage for diversity in forest landscapes (Puettman et al. 2008). A suitable analysis framework for testing juvenile stand sampling methods and their outcomes as predictors of criteria for management success was presented in Chapter 1. A critical element of this system is a set of routines that can: 1. generate large numbers of stands at the scale of and with the complexity of operational cutblocks (tens of hectares), 17 2. reflect variation in ecosystem conditions that might occur as a result of variation in geophysical factors (terrain and soils), 3. reflect variation in stand conditions that can occur as a result of stochastic or patterned occurrence of non-crop competing vegetation, 4. appropriately reflect stand conditions with multiple tree cohorts (a cohort being a collection of trees of the same species with similar size and age attributes), and 5. reflect variation that occurs as a result of management actions. The objective of this paper is to describe a treelist generator that meets these criteria. Data The primary data set for this study is a set of stem mapped plots from an 18 year-old post harvest stand located 47 km south of Dawson Creek British Columbia (55º20’28”N, 120º20’30”W). The stand was a stratified mix of sucker-origin trembling aspen (Populus tremuloides Michx.) with occasional balsam poplar (Populus balsamifera L.), willow (Salix spp.) and green alder (Alnus viridis ssp. crispa (Ait.) Turrill) overtopping planted white spruce (Picea glauca (Moench) Voss). Four 17 m radius plots were subjectively located to capture a range of aspen densities, with a low of approximately 7500 trees/ha in a clumped distribution to a high of 22,000 trees/ha with near uniform ground coverage (Fig. 2-1). Aspen top heights (mean height of the tallest 100 trees/ha based on the full area of the sample plots) were 8.5 to 9.0 m, while spruce top heights ranged from 2 to 6 m depending on the density of the overtopping aspen. Each mapped plot consisted of a 9 m radius core with fully measured sample trees, and an additional 8 m radius mapped buffer. In the core area, broadleaved species were tallied if above breast height (1.3 m) and conifers if taller than 0.3 m. In the buffer area, broadleaves were tallied if greater than 1 cm dbh, and conifers if greater than breast height. Species, breast height diameter (dbh) and position (bearing and distance from plot center) were recorded for all tallied trees. Within the core area, height, age, breast height age and 3-year height increment were recorded for conifer trees, while height, crown width (longest branch in each of four quadrants) and live crown length were recorded for a selected sample of broadleaved trees, with random sampling of four trees from each of 10 even-width diameter classes. For willow and alder tall shrubs, the center of each cluster of stems was recorded for position, along with cluster height and the number of stems in the cluster. The four large stem maps described above were primarily intended to explore competitive relationships within and between tree cohorts. Two additional sets of stem maps were established 18 to validate competitive relationships within the aspen cohort and to better characterize spatial patterns that might be found within aspen cohorts. The first additional data set was developed using seven smaller stem mapped plots established in four cutblocks occurring within a 70 km radius of Dawson Creek BC. All stands regenerated to aspen suckers, with top heights similar to those found in the original four plots. The core areas of these plots were variable in size, with radii (ranging from 2.2 to 4.0 m) sufficient to encompass approximately 40 aspen trees. A mapped buffer with an additional 4 m radius was also included. Measured attributes were the same as for the larger plots. A final set of 15 stem mapped plots was generated by the consultant Timberline Natural Resource Group using stereoscopic pairs of large scale (approximately 1:1450) 70 mm aerial photographs. The photos were borrowed from a much larger March 2006 set used by Alberta Pacific Forest Industries Ltd. for regeneration surveys on operational areas within their Forest Management Agreement area in northern Alberta. Photo pairs were selected from the oldest stands available (8 to 12 years) based on initial subjective assessment to maximize the range of Figure 2-1. Stem maps of a stratified mix of trembling aspen overtopping white spruce. The circles representing each tree are scaled proportionate to tree diameter-at-breast- height (dbh). Many of the smallest spruce do not show up well, particularly in plot B2, and trees with heights less than 1.3 m (dbh = 0) are not displayed. 19 aspen densities across the sample. Images were digitally scanned (7 micron) and set up in a DiAP soft copy photogrammetric system for analysis. Within a 15 m by 19 m rectangular plot (Fig. 2- 2), the position of every tree and tall shrub (Salix spp. and Alnus spp.) was mapped and classified by species. For each plot, the heights and crown radius of 10 trees were measured, and ocular estimates were made for the remaining trees using the measured trees for reference. Clusters of tall shrubs were treated as a single individual for mapping purposes, with the number of stems in each cluster being recorded. Spatial Pattern Generation Rectangular stem maps are generated using inhomogeneous Poisson and thinned lattice processes (sensu Diggle 1983). Various combinations of filters (Fig. 2-3) are used to vary the Poisson point intensity (probability of tree occurrence) and mortality resulting from factors such as terrain position and aspect, operational history, microsite, seed and sucker patterns, and other site features. For each cohort in the stand, users define a set of rules for how stocking will be affected by each filter. The system currently allows for 20 cohorts and six filters. Ten of the available cohorts utilize Poisson processes to generate spatial patterns (primarily naturally regenerated trees) while the remaining 10 use lattice processes (primarily planted trees). For each cohort, users are required to specify size distributions based on height, likelihood of occurrence as affected by their presence in any of the six filters (filter effects are additive), statistics for tree age and whether heights should be adjusted for the effects of overtopping competition. Figure 2-2. Stem map of 15 m by 19 m plot derived from 70 mm aerial photos. The circles representing each tree are scaled proportionate to tree height. 20 Filters The six filters can be divided into two categories: those that control variation at the meso- scale (roughly 10 to 200 m) and those that control variation at the micro-scale (roughly 0.1 to 30 m). The three meso-scale filters are based on stratification of the entire mapped area using three distinct sets of Voronoi polygons, also known as Dirichlet tessellations (Fig. 2-4a): 1. Terrain (or slope) position is determined by dividing each of the first set of polygons into classes based on concentric zones - any point within a polygon can be classified based on its relative distance along a radial vector running from the polygon centroid to the boundary (Fig. 2-4b). Six zones are used to represent the terrain position classes crest, upper, mid, lower, toe and depression. Subdivision of polygons is controlled first by specifying the area percentage (a i ) that falls within upland (crest, upper and mid-slope positions) versus lowland classes, and estimating a relative boundary distance (r) on each radial vector using a circle as a proxy: where i = filter number (1) Further subdivision is controlled by subdividing r 1 (upland classes) and 1-r 1 (lowland classes) into three classes each based on user defined percentages of the radial distance. �� ���� ������ ������������ ������������ ���������������� ����� �� ������ ����� Figure 2-3. Stem maps are created using either a random or lattice point process that are filtered to create desired stand structures. A reasonable analog for this process is an imagined rain of seedlings onto a surface, where the falling seedlings are partially intercepted by a variably permeable screen. In the illustrated case, two small clonal patches of aspen are created with scattered individuals in the surrounding area. Each component of a stand (e.g. species x cohort) is built using separate filter rules. ri = a i 21 In situations where it is not important to recognize terrain position, the zones can also be used for purposes such as defining meso-scale clusters of trees with graduated edges, or multi-cohort clusters such as where a few large trees are surrounded by a cohort of smaller trees. Terrain position classes are further subdivided into two aspect classes (Fig. 2-5a) to allow for effects such as warm, dry south aspects or cold wet north aspects, based on user defined compass bearings radiating from the centroids. The combination of six terrain positions and two aspect classes provides 12 basic units within which an initial set of stocking values can be specified for each cohort. 2. The second meso-scale filter uses two classes that create a set of patches within a matrix, where upward or downward adjustments to stocking can be made for either or both. Subdivision into classes is based on the same procedure as the first split for terrain classes using equation 1. One of the main uses of this filter is to create linear patterns within a stand, such as those resulting from skid trails (Fig. 2-5b). In this case, a 2 becomes a large percentage of the radial distance, and the matrix is restricted to a narrow band straddling the polygon boundaries. Stocking can be adjusted upward or downward to reflect mortality, planting or other effects related to tree occurrence on machine running surfaces. Figure 2-4. Voronoi polygons (a) are the basis for meso-scale landscape variation. Polygons are generated using randomly located points (red dots) as centroids at an intensity specified by the user over an area that includes the mapped unit (dashed line) plus a 200 m buffer. Polygon boundary segments are defined as straight lines whose points are equidistant from two adjacent centroids. Polygons can be used to represent terrain units (b) by assuming centroids are terrain summits and polygon boundaries are terrain depressions. 1 2 3 4 5 6 Polygon centroid = meso-terrain summit Terrain Positions Polygon boundary= meso-terrain trough a b 22 3. The third meso-scale filter is similar to the second, except that stocking adjustments can be specified only for the patches and not the matrix. The primary purpose of this filter is to create large, discrete patches (Fig. 2-5c) such as aspen suckers originating from a single parent or cluster of parent trees. 4. Three sets of micro-scale filters (numbered 4 through 6) are created using small circular units with randomly located centers. The patch density is set by the user, as is the mean and standard deviation of the normal size distribution. These patches are used to add trees to create clusters, subtract trees to create gaps, or to move trees from one set of patches to another set to rearrange the initially random pattern (Fig. 2-6). Within each of the 3 micro-scale filters, overlapping patches are treated as a single larger patch. Figure 2-5. Stem maps illustrating use of Voronoi polygons to define meso-scale features. In all three figures, tree positions are colour coded by terrain position. Additional stratification is illustrated by reduced stocking on a southerly aspect (a), linear features that could represent effects such as skid roads (b), and the addition of a cohort in large patches (c). Note that terrain position underlies all other filters, and the effects of other filters can be varied based on terrain position. For the aspect case illustrated here, stocking reductions are graduated by terrain position, with more severe effects on the higher, drier sites. b c a 23 Stem maps are created in several steps: 1. Two sets of potential points (P R and P L ) are generated for Poisson-based and lattice-based cohorts respectively. For each set, the number of points (n) is determined as the maximum expected point density (points/ha) multiplied by the mapped area (ha). The initial point density (λ max ) for P R is determined as the maximum density specified by the user for any of the combinations of terrain position and aspect class, while the point density for the lattice-based set is specified directly by the user (i.e. expected planting density). Point positions for P R are based on randomly selected x and y values within the rectangular mapped area. For P L , points are determined using a square grid with adjustments using a positional vector based on a random direction and a distance drawn from a normal distribution with mean and standard deviation set by the user. 2. Each potential point is marked by cohort and by its occurrence within the influence zone(s) of any of the six filters. Cohort assignments are probabilistic based on proportions set by the user. Proportions are set separately for each combination of terrain position and aspect class. 3. Potential points are filtered (accepted or eliminated) based on their probability (p) of occurrence. For P R , this becomes: (2) Figure 2-6. Use of micro-scale patches for spatial diversity. Randomly spaced trees (left) have been re- arranged (right) to create micro-scale spatial diversity. In this case, trees removed from one set of patches were replaced in another set with different random centers. pi = λ j λmax × 1− bk,ge j ,k( ) k=2 6 24 and for P L : (3) where: i = potential point number p i = probability of acceptance j = code value for filter 1 representing the combination of terrain position (6 classes) and aspect (2 classes) that point i falls within; values = 1 to 12 λ j = initial point intensity for terrain/aspect code j λ max = maximum point intensity for all levels of j k = filter number for filters 2 through 6 g = code value (patch = 1 or matrix = 2) for filters 2 through 6 representing the filter class (patch versus matrix) that point i falls within b g = density reduction rate (mortality for plantations) for patches versus matrix; for filters 3 through 6, the value of b 2 must be zero; the values for both b 1 and b 2 will be zero if filter is inactive a j = density reduction rate (i.e mortality) by terrain position and aspect class for filter 1 e j,k = multiplier (0 to 1.00) to modify influence of filter k based on overlap with terrain position j 4. For P R only, additional (cluster) points can be created based on point intensities by filter for filters 3 through 6. For filter 3, the initial point density (λ max ) is specified by the user. The final occurrence of each point is based primarily on its presence or absence within the meso-scale cluster described by the filter, and by interactions with filter 1 (multiplier e): (4) For filters 4 through 6, the point density of added trees λ j (j = 4 to 6) is based on a percentage of whatever the initial density was before adjustments for a particular cohort relative to the underlying terrain position. Adjustments are also made for interactions with filter 1. A similar process is therefore required as was outlined in equation 2: (5) where: f k = percentage of original point density (based on terrain position and aspect) to be added pi = 1− a j( ) × 1− bk,ge j ,k( ) k=2 6 pi = λ j λmax × 1− 1− fkbk,ge j ,k( ) k=3 6⎛ ⎝ ⎜⎜ ⎞ ⎠ ⎟⎟ 25 Assigning Tree Attributes Tree attributes are determined using various combinations of (i) random draws from a normal distribution of heights, (ii) Wilson’s (1946) Spacing Factor (SF) [where SF = mean inter-tree distance / mean height], and (iii) any previously defined attributes. Tree attributes are assigned in the order of: height, diameter, age, breast height age, and crown dimensions. Tree height was chosen as the “driving” variable, as it is believed to be the least affected by community competition (Lanner 1985). Spacing Factor was chosen as the competition index as it is derived from tree locations and height, these being the first two attributes generated in the stem map. Separate normal distributions of height are defined for each cohort by specifying the mean and standard deviation. It is recognized that other and more flexible model forms such as a Weibull may be appropriate in many cases, but the functionality to specify such forms has yet to be implemented. In the calibration data set used for this project, 8 to 12 year old aspen stands in the aerial photo-based plots tended to have height distributions that were either symmetrical or left-skewed (the modal value was smaller than the mean), while the 18 year old aspen cohort in the larger ground based plots tended to be symmetrical or right-skewed. A possible explanation of the difference is the longer period of self-thinning in the older stand, a process which would preferentially remove smaller stems. The spruce cohort in the larger ground-based plots had distributions that ranged from symmetrical to left-skewed (one of four was acceptably represented by a normal distribution). Where height is anticipated to have been affected by overtopping competition by trees from other cohorts (Fig. 2-7), height can be adjusted downward using a multiplier (range = 0 to 1). The multiplier is determined as a function of one of three variants of SF where the trees used to calculate the index can include a) all trees, b) trees from specific cohorts or c) trees taller than the subject tree. Multiplier functions can be based on any of five equation forms (Table 2-1), using calibration parameters set by the user. Figure 2-7. Conifer height growth is affected by overtopping competition from a broadleaved tree patch. Such effects are simulated using a height multiplier determined using a function based on spacing factor (SF) to quantify the magnitude of the effect. 26 Tree diameters are assigned based on allometric functions using height as the independent variable. As with predicting height adjustment multipliers, several functional forms are available (Table 2-1). Total age for planted trees is a fixed value as specified by the user. Total age for natural trees is calculated primarily as a function of height, assuming that the oldest trees are also most likely to be the tallest. The first step, then, is to interpolate age from height using user specified minimum and maximum ages: (6) Given that height distributions are assumed to be normal, this process will result in an initial distribution of ages that is also normal. A second step adjusts this linear distribution upward to generate a distribution that is skewed to older ages (assuming that more trees germinate earlier rather than later, and that older (bigger) trees are more likely to have survived early stages of self-thinning). minminmax minmax min )( )( )( AgeAgeAge HeightHeight HeightHeight Ageinitial +−×− −= Table 2-1. Model forms available for predicting multipliers (0≤y≤1) for adjusting tree heights relative to overtopping competition and for predicting diameters. In all cases, x is a measure of competition as expressed using spacing factor. Spacing factor can be determined using all trees, or it can be restricted to using either certain cohorts or to trees that are taller than the subject tree. Function Form Ht Adj. DBH Weibull where b 1 is fixed at 1.00 √ √ Chapman Richards where b1 is fixed at 1.00 √ √ Generalized Logistic where b 1 = b 4 = 1.00 √ √ Gompertz where b 1 = 1.00 √ √ Segmented Linear for x > b 1 for x < b 1 √ Polynomial √ Exponential Polynomial √ Linear √ y = b1 1− e −b2x b3( ) y = b1 1− e−b2x( )b3 y = b1 b4 + e b2 −b3 x y = b1e−e b2 −b3x y = b2 = 1 y = b2 − b3 (b1 − x) y = b1 + b2 x1 + b3 x2 + b4 x3 + b5 x4 log(y) = b1 + b2 x1 + b3 x2 + b4 x3 + b5x4 y = b1 + b2 x 27 Adjustments are based on random draws from a cumulative beta distribution (Fig. 2-8), such that most trees get a relatively small adjustment, but a few adjustments are quite large: (7) Parameters for the beta distribution are currently fixed in the model code as a crude approximation and are not based on a fit to sample data. Breast height age is calculated from total age and height, assuming that the relationship between height and age for young trees is linear above breast height. While this assumption will not hold over a very wide range of ages, it should be acceptable for the juvenile stands for which the simulator was developed. Breast height age (Age BH ) is then calculated as: (8) where interceptx is the x intercept of the linear fit for the age versus height relationship above breast height (Fig. 2-9). Implementation, Memory Limitations and Computational Speed The simulator is constructed within the Visual Basic for Applications (VBA) environment of a Microsoft Excel® spreadsheet. The VBA environment was attractive in that it offered ready access to all of Excel’s built-in features such as file structure, function controls and charting. Another advantage was that it also provided ready access to the programming code, so that various users can easily make their own modification based on case-specific needs. The combination of Excel and VBA also makes it very easy to set up batched simulations: with a Figure 2-8. Cumulative beta distribution where alpha = 4 and beta = 1.5. 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Random Value C u m u la ti v e B e ta ( 4 ,1 .5 ) )( max initialdrawinitialfinal AgeAgeBetaAgeAge −+= )( )3.1( xTotalBH InterceptAgeHeight Height Age −−= 28 small amount of up-front preparation, a large number of stands can be generated without periodic user inputs. The major disadvantage is limits on memory. Excel 2003 and 2007 are limited to 1 GB and 2 GB of memory respectively for the Excel process (Williams 2009), although both versions appear to have similar effective limits on the maximum cumulative array sizes within the VBA environment. This limit results in a maximum stand size of approximately 150,000 trees. “Out of memory” errors may increase in frequency when multiple simulations are run in succession because of memory fragmentation (also referred to as memory leakage). For many processes, Excel requires large blocks of contiguous memory (Williams 2009), which may become increasingly unavailable after successive model runs. A 64-bit version of the planned Excel 2010 is expected to vastly increase access to physical memory and largely eliminate these limitations. There is also a speed penalty associated with using the VBA environment. VBA code is not pre-compiled into an executable binary and must be compiled on-the-fly. As a result, a comparable model compiled as an executable program would likely run much faster. Actual times to generate stands will vary with computer configuration and the number of points being evaluated in order to generate the desired number of trees. On a computer using an Intel Core Duo processor, 2 GB of RAM, Windows XP Pro and Excel 2003, evaluating 160,000 points to generate 120,000 trees takes approximately 40 minutes. Stands of 2000 trees can be generated in as little as 1 to 3 minutes. Figure 2-9. The calculation of age at breast height (1.3 m) from total age and height assumes a linear approximation of the height growth curve for heights above 1.3 m. 0 1 2 3 4 5 0 5 10 15 Age (yrs) 1.3 m Total Age Years to BH AgeBH H ei g h t (m ) 29 Spatial Pattern Evaluation Spatial patterns produced using the stem map generation routines are evaluated using two graphical functions. The first of these is the L(t) function, a standardized variant of Ripley’s K(t) function (Ripley 1977) that was first proposed by Besag (1977) in comments appended to Ripley’s paper. The K(t) and L(t) functions are frequently used in ecological literature and other fields where the spatial patterns of events are important. The second function is a new application of Wilson’s (1946) spacing factor (SF). The L(t) function is evaluated and presented using a variation as recommended by Fortin and Dale (2005): (9) which contrasts the observed number of points in a circle of radius t with the expected number given complete spatial randomness (CSR). For a range of values of t, a chart can be prepared which evaluates perceptions of spatial randomness at various scales (Fig. 2-10). Up to 10 L(t) charts can be generated using a set of circular observation windows of radius r (as specified by the user where r ≥ 15 m) that evaluate local neighborhoods within the simulated stand. Values of t are presented ranging from 0.05 to 15 m, although caution should be exercised in interpreting outcomes where t ≥ 0.5r, as the number of pairwise distances that can be evaluated for the larger values of t becomes limited. Observation windows can be located either randomly L(t) = t − K (t) ≠ Figure 2-10. Examples of L(t) functions. In the top chart, trees are overdispersed (more evenly spaced than would be expected with CSR) at scales less than 1 m, possibly due to self-thinning. At scales from 1 to 12 m, the stand is clustered, and becomes overdispersed again for larger scales. In the bottom chart, the stand is overdispersed at scales out to 9 m, a pattern that would be typical of regularly spaced plantations. The clustering at larger scales might indicate large unstockable gaps or patches of mortality in the plantation. 1 0.5 0 -0.5 -1 L(t) t 1 0.5 0 -0.5 -1 L(t) t 0 5 10 15 0 5 10 15 30 or on a grid, recognizing that possible locations will be limited to those where the entire window falls within the simulated stand. It should also be recognized that observation windows may overlap, particularly when their locations are random, resulting in discrete depictions of spatial pattern that may not be independent. While the K(t) and L(t) statistics have been very popular for evaluating spatial pattern, one of their limitations is the assumption that all points are equivalent. Where attributes of the points (or marks) are important for answering the questions of interest, these statistics start to lose their value. This is frequently the case in ecology and forestry where competition between and within cohorts of plants is of importance. In these cases, a description based on competition indices may be more useful. While not truly a measure of spatial pattern, histograms and cumulative distributions of a competition index (Fig. 2-11) can provide some valuable additional information on spatial intensity with incorporation of mark (tree size) values. In the current simulator, measures of SF are evaluated at each of 100 points on a grid within the simulated stand, using an observation window with a radius that is 70% of the top height of the dominant stand cohort. Note also that dominant height is replaced by mean height in the calculation of SF to avoid problems associated with the tallest cohort in a stand having a relatively sparse canopy. The main value of spacing factor (or other competition index) histograms and cumulative distributions is not for spatial analysis but for management interpretation. With prior experience in the interpretation of spacing factor, forest stand managers can evaluate the range of competitive conditions that exists within a stand relative to their expectations for stand development toward desired management objectives. The distributions do a much better job Figure 2-11. Cumulative distribution for spacing factor. For any x-value (spacing factor), the corresponding y-value is the percentage of observations that were less than or equal to the x- value. Spacing Factor = Mean Inter-Tree Distance / Mean Height 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 Spacing Factor Cumulative Frequency (%) C u m u la ti ve F re q u en cy ( % ) 31 of indicating possible future outcomes than do the more traditional whole-stand, single-value summaries of these statistics. L(t) Case Studies The L(t) functions for the aspen/cottonwood components of a number of young sampled stands are illustrated in Figs. 2-12 and 2-13. These plots are relatively typical of the range of conditions that was observed in the larger data set. In general, sampled stands produced statistics consistent with either random or clumped distributions over the range of meaningful scales that could be evaluated using the available plots. By contrast, a 1 ha simulated stand is illustrated in Fig. 2-14, along with four L(t) functions for local portions of the stand generated using 20 m radius circular sampling windows. The observation window size was selected to be midway between those provided by the field and aerial plots. In comparing the spatial statistics for the sampled and simulated stands, the following observations are drawn: 1. It appears that close analogues to the sampled aspen stands can be generated assuming an initial random distribution of trees with the addition of either complete or partial voids of varying intensity and size. Figure 2-12. Spatial distributions as characterized by the L(t) function for the 4 large stem maps illustrated in Figure 2-1. The scale of t is in metres. Plot B1 is moderately clumped across a wide range of scales up to 14 m, while plots B2 and B4 exhibit conditions much closer to complete spatial randomness over the same evaluation distances. Plot B6 shows moderate to strong clumping at scales below 10 m. All four plots show a tendency to over-dispersion at broader scales, indicating a relative absence of broader scale gaps and clusters. B1 -0.8 -0.4 0 0.4 0.8 0 2 4 6 8 10 12 14 16 t L (t ) B2 -0.8 -0.4 0 0.4 0.8 0 2 4 6 8 10 12 14 16 t L (t ) B4 -0.8 -0.4 0 0.4 0.8 0 2 4 6 8 10 12 14 16 t L (t ) B6 -0.8 -0.4 0 0.4 0.8 0 2 4 6 8 10 12 14 16 t L (t ) 32 Figure 2-13. Selected stem maps (scale in metres) and corresponding L(t) functions for three selected plots in 8 to 12 year-old aspen stands, illustrating the range in spatial variation observed. The top plot illustrates a distribution with near CSR conditions across all sampled scales (0 to 7.5 m). The middle plot illustrates a condition wherein moderate clumping occurs across all but the most narrowly focused scales, while the bottom plot illustrates a high degree of clumping at scales of 2 to 7 m, and moderate clumping at narrower scales. The scale of t is in metres. 0 22 0 18 0 22 0 18 0 22 0 18 11141_5 -0.8 -0.4 0 0.4 0.8 0 2 4 6 8 t 02951_5 -0.8 -0.4 0 0.4 0.8 0 2 4 6 8 t 36171_4 -0.8 -0.4 0 0.4 0.8 0 2 4 6 8 t L (t ) L (t ) L (t ) 33 2. Several of the sampled stands exhibit a sharp drop in the L(t) statistic for very small distances (< 1 m) that is not evident in the simulated stands – subsequent simulations suggest this can be easily remedied with an additional set of micro-scale (< 2 m ) voids. 3. Neither the sampled nor the simulated stands show any evidence of overdispersion at small distances (i.e. patterns requiring a hard-core process for simulation). It is important to note that the L(t) function in no way can be used to prove similarity of spatial pattern. Cautions regarding evaluation of spatial equivalence using summary statistics were voiced by Hurlbert (1990) for simple assessments of randomness using mean and variance, and are continued by authors such as Illian et al. (2008) for more complex cases. The importance of this principle is illustrated by imagining the reverse exercise of reconstructing a population from sample statistics alone. Pommerening and Stoyan (2008) illustrate a case where spatial structure can be adequately reconstructed using spatial summary statistics, but only by combining the information captured in six different statistics, and only where the summary statistics are generated using a sample that is a large proportion (~ 25%) of the original population. Figure 2-14. Stem map and associated L(t) functions evaluated in four localized sampling windows (red circles). The stem map was created starting with a CSR point density of 15,000/ha. The point density was reduced to 40% of the original in one set of circles (λ=25/ha, radii: mean = 4.0 m, sd = 3.0 m) and to 20% in a second set (λ=20/ha, radii: mean = 4.0 m, sd = 2.0 m). The scale of t is in metres. -0.8 -0.4 0 0.4 0.8 0 2 4 6 8 10 12 t Top Left -0.8 -0.4 0 0.4 0.8 0 2 4 6 8 10 12 t Top Right -0.8 -0.4 0 0.4 0.8 0 2 4 6 8 10 12 t Bottom Left -0.8 -0.4 0 0.4 0.8 0 2 4 6 8 10 12 t Bottom Right 34 Discussion The stem map generation routines in this simulator were developed using a mathematically naïve approach. The basic philosophy was to create a wide range of spatial patterns that might be observed under field conditions, through user-based experiential mimicry of the causal biophysical factors such as soil moisture, nutrient and microclimate conditions (terrain position and aspect), soil disturbance patterns, seed rain and competing vegetation. The approach taken requires users to imagine the type of distributions they wish to create and the factors that might produce them, and achieve those structures through “brute force” manipulation. Once the pattern is generated, it is then quantified and summarized using spatial statistics. This is in contrast to other approaches in the literature such as Kokkila et al. (2002) and Pommerening and Stoyan (2008) where more mathematically elegant approaches are used primarily to achieve conditions mimicking pre-determined summary statistics and/or spatial samples. Pretzsch (1997) takes an approach with elements of both philosophies. One of the primary design criteria that shaped the architecture of the current simulator was the requirement for non-stationarity, primarily to reflect variation in terrain, soils and hydrologic conditions. This in turn led to the adoption of the nonhomogeneous Poisson process and the use of Voronoi polygons to generate simulated maps of such features. The same logic was then carried through to the simulation of large and small clusters. The nonhomogeneous Poisson process used in the current simulator has similarities to that described by Pretzsch (1997) for the STRUGEN stem map generator, but differs in several important ways. Where STRUGEN assumes a largely stationary landscape, the current algorithms were explicitly designed to simulate at least some effects of non-stationarity at scales of tens to hundreds of metres, primarily as imposed by terrain and ecosystem variability. Considerable overlap can be recognized between the systems in the handling of meso- and micro-scale patches, species intermingling and anthropogenic linear features. Three features that the STRUGEN generator includes that have no direct analogues in the current simulator are a hard-core process to limit the probability of pairs of trees occurring in close proximity, explicit control of species intermingling to match a pre-determined index, and probability distributions that affect filter density relative to the center (or center line for linear features) of a cluster. The latter feature allows for “soft” or dispersed edges to spatial features, although the degree to which such a feature is important has not been assessed. The Gibbs process approach to stem map generation such as that illustrated by Kokkila et al. (2002) shows great promise in the future for generating highly realistic stem maps that can account for all of the design criteria for the current project. However, such methods have yet to be demonstrated in the literature at scales beyond a single hectare, and the calibration of 35 the models for each new simulation requires specialized mathematical knowledge regarding the conditions to be simulated (i.e. pairwise field potentials and interactions with marks and background fields). As such, these methods are not yet suitable for the purposes for which the current simulator was primarily developed: the rapid generation of a wide variety of stand conditions at operational scales, primarily for the purpose of simulating operational sampling schemes. The stand reconstruction approach employed by Pommerening and Stoyan (2008) is also capable of producing stem maps that satisfy most of the design criteria for this project. However, that approach is seriously limited by its needs for large quantities of spatial data for each new stand to be generated. The primary application of such reconstruction methods appears to be one of cost savings in the mapping of very large forest plots: a highly realistic stem map can be generated by field sampling only a portion of a stand. All of the tree generation methods discussed here have at least one common feature that is superior to the approach presented in this chapter: they employ either a hard-core process or an alternate process that can, where appropriate, limit the proximity of any pair of trees. For many stand types, such processes may be very important; the process of self-thinning in even-aged stands, for example, leads to surviving trees that are spaced farther and farther apart as stands grow. However, it is debatable whether or not the lack of such a feature in the stand generator presented here is an important omission. This model is primarily intended to generate stands for use in regeneration sampling simulations, an application where hard-core distances are not particularly important (it seldom matters that two trees are 20 cm apart as opposed to 1 m). It could also be argued that there is little need for simulating hard-core distances in young stands. For example, in the L(t) function plots for 18 year-old and 10 to 12 year-old aspen in Figs. 2-12 and 2-13 respectively, there is no evidence of over-dispersion of trees even at scales of 0.05 m. There are at least two possible explanations for the lack of over-dispersion detected: (i) a pattern of over-dispersion at micro-scales has not yet developed because insufficient self-thinning has taken place and (ii) a pattern of over-dispersion exists at micro-scales but was masked by the measurement protocol. In the latter of these explanations, the stem mapping protocol assumes that the center of the stem at breast height is an acceptable measure for the position of the tree as a whole. Given that trees often develop at least some degree of lean and crown asymmetry as they compete for growing space, it is surmised that a map of crown centroids would elicit a greater degree of over-dispersion at micro-scales than would the map of stems at breast height. Therefore, measures of clustering and over-dispersion at micro-scales (less than the mean inter- tree distance) should be interpreted with caution. 36 Extending this simulator to a wide range of stand conditions, particularly older or multi-layered stands, may also be limited by the lack of spatial control by marks. While the intermingling of trees by species and size can be controlled in a crude fashion by varying the point densities using filters for multiple cohorts, it is not difficult to imagine conditions that would be difficult to mimic. One example is the exclusion of small trees from immediately below the crown of a larger competitor (a common occurrence in dry Douglas-fir [Pseudotsuga menziesii (Mirb.) Franco] forests of the interior of British Columbia). Another is the need to apply a hard-core process to one component of the stand (e.g. larger, older trees) while the remainder are distributed more randomly. Additions to the simulator that may be desirable would be a wide range of further spatial statistics for quantifying and contrasting simulator outcomes. One that would be high on the list would be the pair correlation or g function as recommended by Penttinen et al. (1992) and Illian et al. (2008). These authors point out that the g function is related to the K , the g being the density function and the K being the cumulative function for the same distribution. In this sense, the g function is more appropriate for exploratory analysis while the K and l functions are more appropriate for statistical tests. Further statistical summaries that may prove useful include the mark correlation function as applied by Penttinen et al. (1992) and Pielou’s (1977) index of segregation as employed by Pretzsch (1997). Also of interest would be a selection of different competition indices that could be charted in place of SF in the cumulative distribution function from Fig. 2-11. Conclusions The described forest stand stem map generator was designed to generate a large number of realistic juvenile stands with a wide diversity of conditions resulting from biophysical (terrain and ecosystem) conditions, (primarily anthropogenic) disturbance patterns, vegetative competition and other stochastic but unspecified factors. It largely succeeds in its stated goal when the evaluation scales for spatial distribution are greater than the mean inter-tree distance, but may fail if the user is concerned with the subtleties of nearest neighbor distances. It is here that the stand generator makes a compromise by sacrificing spatial precision at micro-scales in favour of computational speed and achieving spatial diversity at meso- to macro- scales. For the purpose that the stand generator was developed (forest sampling simulations in juvenile stands), it is believed that this compromise is acceptable. 37 References Besag, J.E. 1997. Contribution to the discussion of Dr. Ripley’s paper. J. Roy. Stat. Soc. B. 39: 193-195. Diggle, P.J. 1983. Statistical analysis of spatial point patterns. Academic Press, London. 148 p. 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. Manage. 94: 129- 140. Fortin, M.J., and Dale, M.R.T. 2005. Spatial analysis: a guide for ecologists. Cambridge University Press, Cambridge, UK. 365 p. Hanus, M.L., Hann, D.W., and Marshall, D.D. 1998. Reconstructing the spatial pattern of trees from routine stand examination measurements. For. Sci. 44: 125-133. Hassani, B.T., LeMay, V., 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. Hurlbert, S.H. 1990. Spatial distribution of the montane unicorn. Oikos 58: 257-271. Illian, J., Pentinnen, A., Stoyan, H., and Stoyan, D. 2008. Statistical analysis and modelling of spatial point patterns. John Wiley & Sons Ltd., West Sussex UK. 534 p. Kokkila, T., Mäkelä, A., and Nikinmaa, E. 2002. A method for generating stand structures using Gibbs marked point process. Silva Fenn. 36: 265-277. Lanner, R.M. 1985. On the insensitivity of height growth to spacing. For. Ecol. Manage. 13: 143- 148. MacLeod, D. 1977. Accuracy of regeneration stocking estimates: tests with simulated data. For. Chron. 53: 77-81. MØller, J., and Waagepetersen, R.P. 2007. Modern statistics for spatial point processes. Scand. J. Stat. 34: 643-684. Pentinnen, A., Stoyan, D., and Hentonnen, H.M. 1992. Marked point processes in forest statistics. For. Sci. 38: 806-824. Pielou, E.C. 1977. Mathematical ecology. Wiley, New York. Pommerening, A., and Stoyan, D. 2008. Reconstructing spatial tree point patterns from nearest neighbour summary statistics measured in small subwindows. Can. J. For. Res. 38: 1110- 1122. 38 Pretzsch, H. 1997. Analysis and modeling of spatial stand structures. Methodological considerations based on mixed beech-larch stands in Lower Saxony. For. Ecol. Manage. 97: 237-253. Pretzsch, H., Biber, P., and Dursky, J. 2002. The single tree-based stand simulator SILVA: construction, application and evaluation. For. Ecol. Manage. 162: 3-21. Puettmann, K.J., Coates, K.D., and Messier, C. 2008. A critique of silviculture: managing for complexity. Island Press, Washington D.C. 189 p. Ripley, B.D. 1977. Modelling Spatial Patterns. J. Roy. Stat. Soc. B 39: 172-212. Stauffer, H.B. 1977. Application of the indices of nonrandomness of Pielou, Hopkins and Skellam, and Clark and Evans. BC-X-166. Canadian Forest Service, Pacific Forest Research Centre, Victoria BC. 32 p. Stoyan, D., and Penttinen, A. 2000. Recent applications of point process methods in forestry statistics. Stat. Sci. 15: 61-78. 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. Manage. 177: 277-285. Wilson, F.G. 1946. Numerical expression of stocking in terms of height. J. For. 44: 758-761. Williams, C. Excel memory limits. URL: http://www.decisionmodels.com/memlimitsc.htm. Accessed Aug. 12, 2009. 39 Chapter 3. Observations on the Relationship Between Regeneration Stocking and Yield6 Introduction Forest stocking is a measure of the amount of tree vegetation, expressed in any of many different units, relative to an ideal desired condition (Smith et al. 1997). In the current context, it is used to generally describe a condition relative to one where all of the available resources needed by trees to grow and thrive are or will be effectively accessed. Foresters working in many jurisdictions have put great effort into the quantification of stocking in regenerated stands (see Chapter 1). One of the primary reasons for this work is the implicit recognition that stocking has a large impact on yield (Bella 1976, McWilliams and McWilliams 2009). Stands that are not making full use of the available resources (growing space) will produce less than they otherwise could; conversely, unstocked areas in the stand (or voids) are not contributing to stand growth. Despite this recognition, explicit links between regeneration stocking and future yield have remained tenuous for most of the last 100 years (Lutze et al. 2004). Doucet (1991) noted that, despite concerns about yield, most early applications of measured stocking and the setting of thresholds for acceptable stocking (stocking standards) utilized subjective perceptions of desired stand conditions. Throughout the 20th century, stocking has been measured primarily against an idealized number of well distributed trees per unit area. Following on the work of Pound and Clements (1898) who first proposed quantification of plant occurrence using quadrat plots and Lowdermilk (1927) who first described the use of stocked quadrats in forestry, Haig (1931) visualized a gridded stand such that each square was equivalent to the vertical crown projection of a dominant or codominant mature tree that has largely been unaffected by its neighbours. He then recommended (pg. 748) that such an area be used for sampling, following “… the premise that the size of the unit adopted must be such that one established seedling per square would give full yields at maturity”. Other authors (e.g. Cowlin 1932, Wellner 1940, Bella 1976, Boyer 1977) have followed up on this philosophy, relating stocking levels and ideal quadrat size to both minimum numbers of well distributed trees per unit area and total numbers of trees per unit area. In a few cases over the last decade, various jurisdictions have started to develop stocking standards that are explicitly linked to yield (e.g. Martin et al. 2002, 2005). In the absence of explicit links through modeling exercises, the easiest assumption to make about the relationship between stocking and yield is that unstocked areas do not contribute to volume production, 6 A version of this chapter will be submitted for publication. Farnden, C. 2010. Observations on the relationship between regeneration stocking and yield. 40 leading in turn to a 1:1 relationship between percent unstocked area and percent volume losses from the expected value for a fully stocked stand. Such an assumption is analogous to the ratio adjustment procedures commonly applied to predict actual yields from normal yield tables. Grant (1951) warned of the dangers of this approach, citing the fact that there is no unique measure of unstocked area. Measured values for stocking vary with measurement parameters, both in terms of the smallest gap that is detected and the related issue of what proportion of each gap is assumed to be accessed by adjacent trees. One variation on this theme is the application of percent stocking adjustments to variable density predictions of yield. The BC Ministry of Forests (1998a, b) uses Operational Adjustment Factors (OAF’s) to net down yields from uniformly stocked stands to account for large voids, using stocking measures based on quadrats of either 22.9 m2 or 40.7 m2 depending on species. Likewise Newton (1998) suggested using stocking values from standard regeneration stocking surveys (10.1 m2 quadrats) to adjust yields predicted using a stand density management diagram. However, in both of these examples it is unclear the degree to which the assumption of a 1:1 stocking to yield relationship has been tested. Several studies have addressed explicit linkages between stocking and yield by contrasting current stocking in pole-sized to mature stands with backcast estimates of stocking; or estimates of stocking that would have been measured many years previous at an earlier stage of stand development. Perhaps the first of these was described by Staebler (1949). His study measured stocking and standing volume in 15- to 35-year old Douglas-fir stands, and correlated percent stocking to the percent of normal yield expected for those stands based on published yield tables by site class. Using samples of 100 groups of four quadrats, Staebler was able to predict yields to within 15%, 19 times out of 20. However, his study was limited to forecasts of future percent normality only over a fixed 15-year time horizon starting from 5, 10 or 15 years of age. Braathe (1976, 1982) summarized work in Norway that correlated percent unstocked plots (termed zero-square percentage) with yield at 55% of culmination age for Norway spruce (Picea abies (Linn.) Kars.). Similar to Staebler (1949), measured yields were contrasted to those in existing yield tables to provide relative values. In this case, an R2 value of 0.98 was found between percent stocking and relative yield at the 55% reference age. Doucet (1991) stem-mapped 35 semi-mature jack pine (Pinus banksiana Lamb.) and 28 semi-mature black spruce (Picea mariana (Mill.) BSP) stands (0.2 ha each) for a wide range of stand densities and site qualities in Quebec, Canada. Local volume equations were produced estimating yield either at age 50 (jack pine) or 65 (black spruce) as influenced by stocking as 41 estimated for juvenile stands using both 4 m2 and 8 m2 quadrat plots. R2 values in all cases exceeded 0.9. The fit quality was similar using either plot size for estimating stocking, although the curve shape for the larger plots provided a more nearly linear form, suggesting “…that a 8- m2 quadrat is more closely related to the area needed for a tree to reach maturity than is a 4-m2 quadrat” (Doucet 1991, pg. 186). When backcasting stocking, the assumption must be made that early stand conditions can reasonably be determined from those extant at the time of examination at mid-rotation or later ages. Doucet (1991) discussed the risks inherent in this assumption, and suggested that stocking will most commonly increase slightly from establishment to ages 25 or 30, and will most commonly be reduced by at most a total of 3% thereafter to approximately age 50. Similar findings to the latter point, found using modeling of percent stocking over time, were presented by Yang et al. (2008). Doucet’s (1991) main arguments were that mortality occurring in pole- sized stands would remain detectible for many years and that self-thinning mortality mainly occurs in dense clumps and has little effect on stocking, given that only a single tree remaining is sufficient to class a plot as stocked. The final obvious progression to linking stocking and yield is through growth modeling. A stand level approach was taken by Yang et al. (2008) and Yang and Huang (2008) using panel data from permanent sample plots in Alberta, Canada. They developed relationships, by species and site, relating changes in percent stocking to stand age. In turn, this percent stocking model is intended to be a component of a larger growth and yield model (GYPSY) that uses regeneration stocking both as a user-supplied descriptor of initial condition and as a stand attribute affecting growth through time (Huang et al. 2009). In the most recent version of GYPSY, the ability to evaluate regeneration stocking effects on yield is an implicit design parameter. In contrast to using a stand level model, an alternative approach is to simulate yields using an individual tree, distance dependent growth model such as TASS (Mitchell 1975) or SILVA (Pretzsch et al. 2002). In this type of exercise, point pattern generators such as STRUGEN (Pretzsch 1997) or the stand generator described in Chapter 2 are used to establish either realistic or meaningfully abstract tree spatial distributions. Those distributions serve as a stem map both for the initiation of growth simulations and as the basis for regeneration survey simulations (see Chapter 1). In this manner, each survey outcome can be linked explicitly to a yield prediction, and simulations for a large number of stand variations can be used to elicit general trends. Given that different survey parameters (e.g. quadrat size) result in different estimates of stocking for a given stand, it is logical to assume that the shape of any curve representing the relationship between stocking and yield will change with the survey parameters used. If a relative scale for both stocking and yield is assumed, a stand with zero stocking will result in 42 zero yield, and under any reasonable set of useful survey parameters a stand with 100% stocking will produce 100% of potential yield. If any two variations in survey parameters between these extremes result in different stocking values for a given stand, the shape of the entire curve must be different. Explicit or implied links to yield play a large role in the setting of stocking standards used to evaluate statutory or contractual obligations for reforestation in many jurisdictions (Hosie 1953, Lutze et al. 2004). Stocking effects on yield can also be important inputs into forest estate analyses such as timber supply forecasts, carbon budget modeling and scenario planning for sustainable forest management (Alberta Reforestation Standards Science Council 2001, Ontario Ministry of Natural Resources 2001). Therefore, a greater understanding of the links between stocking measures and anticipated yield is desirable. In addressing this problem, the current study follows the decision framework outlined in Chapter 1; a set of simulated juvenile stands is generated which in turn is the basis for both growth simulations to determine yields and sampling simulations to generate silviculture survey summary statistics. Models are then constructed from the outcomes, whereby survey statistics are used to predict future yields. The primary objective of this study is to use simulated case studies to explore important issues in the quantification and interpretation of regeneration stocking with respect to yield prediction. Important aspects of this relationship that will be explored include the impacts of species, spatial pattern of regenerated trees, and growth rates. Methods Species and site productivity effects on the relationship between stocking and yield were explored by selecting two conifer species for simulation, each having very different silvical characteristics and crown forms. Very different sites were also simulated for the two species to further differentiate the cases. Western redcedar (Thuja plicata Donn) is generally considered to be highly shade tolerant (Krajina et al. 1982) with a broad, conical crown shape, while lodgepole pine (Pinus contorta var. latifolia Dougl. ex. Loud) is shade intolerant (Krajina et al. 1982) with a narrow cylindrical crown. Western redcedar was simulated using a site index typical of its native range in coastal British Columbia (28 m), and lodgepole pine was simulated using a roughly median site index (20 m) for inland areas of BC. For each species, 65 stands were generated using the stand generation routines described in Chapter 2. For each species, 51 of 69 stands were generated to mimic planted stand conditions, with three planting densities for each species (500, 800 and 1200 for western redcedar; 1200, 1800 and 2400 for lodgepole pine). A 20% random deviation from perfectly square spacing was utilized to mimic variation in tree position caused by planter selection of the best micro-sites. One case for each species and density combination was created assuming no mortality, nine were created with random mortality (Fig. 3-1a) ranging from 5% to 85% in 10% increments, and eight were created with varying combinations of random and clustered mortality (Fig. 3-1b, Table 3-1). For both species, tree heights were drawn from a normal distribution (mean = 2.0 m, standard deviation = 0.2 m and heights conditioned to a maximum of 3 standard deviations), with a uniform tree age of 9 years. Figure 3-1. Examples of stem maps showing three primary variations in spatial pattern: a) lattice (square spacing) with random mortality, b) lattice with clustered and random mortality and c) filtered Poisson with clustered re-distribution (representing naturally regenerated stands). The three stands illustrated were assessed as having similar stocking: approximately 66% by 16 m2 quadrats. Table 3-1. Frequency of gaps, approximate affected area and associated random mortality outside of gaps for stands affected by a combination of clustered and random mortality. Mortality rate in gaps was set at 90%. Gap Frequency (number/ha) Approximate Area (%) Random Mortality (% of trees) 5 10 5 11 21 10 17 30 15 24 38 20 34 47 25 45 56 30 60 66 35 44 A further 18 stands were generated for each species to mimic patterns of natural regeneration, with stand densities controlled primarily by terrain position (Fig. 3-1c). In the base case for full stocking, stands for both species were generated using a Poisson distribution with 7000 trees/ha, followed by a partial redistribution where 70% of the trees were removed from a set of randomly located circular areas (mean radius = 8 m, sd = 2 m, density = 40 points/ha) and replaced in a similar set of circular patches with different random locations. Subsequent stands were developed using an 85% redistribution of trees in random clusters, a 75% reduction in trees on linear features representing skid trails, and varied stand density changes from the base case by terrain position (Table 3-2) as a proxy for variation in germination rates related to soil edaphic conditions. All 36 stands mimicking natural regeneration used the same pattern of terrain units and skid trails. Tree heights were drawn from a normal distribution (mean = 1.4 m, sd = 0.4 m and heights conditioned to a maximum of 3 standard deviations), and ages from 1 to 9 years (based on height) reflecting an 8-year pattern of germination. Growth of all stands was simulated using TASS, an individual tree, spatially explicit model originally developed by Mitchell (1969, 1975) and further refined and expanded over many years by the BC Ministry of Forests’ Research Branch. This model simulates the growth Table 3-2. Stand densities by terrain position for simulated stands intended to mimic natural regeneration patterns, along with area breakdown by terrain position. Stand Density By Terrain Position Variant Crest Upper Slope Mid- Slope Lower Slope Toe Slope Depression Base 7000 7000 7000 7000 7000 7000 1 7000 7000 7000 7000 7000 7000 2 7000 6997 6918 6422 4684 200 3 6500 6484 6290 5555 3749 200 4 6000 5942 5576 4652 2933 200 5 5500 5346 4792 3775 2254 200 6 5000 4684 3980 2976 1708 200 7 4500 3969 3193 2287 1283 200 8 4000 3240 2480 1720 960 200 9 3500 2543 1869 1272 720 200 10 3000 1920 1372 930 546 200 11 2500 1394 984 677 422 200 12 2000 975 694 495 335 200 13 1500 657 484 367 276 200 14 1000 427 338 280 236 200 15 700 300 235 200 165 140 16 490 210 165 140 115 95 17 343 160 115 100 80 60 Approx. Area Breadown (%)* 6 23 39 19 6 7 * based on standardized map of terrain units derived from zoned Voronoi polygons (see Chapter 2) 45 and interaction of neighboring tree crowns, with the resulting crown volume being a major determinate of growth. The model has been extensively tested against a wide range of thinning and espacement trials (e.g. Goudie 1998) and has been found be accurate for simulating tree growth relative to competitive conditions. All growth simulations were run to a stand age of 100 years with a yield table generated at annual increments and a full tree list including tree position, height and diameter and mean crown width at 20-year increments. For most comparisons of yield, an arbitrary reference age of 80 years was selected as a reasonable approximation of rotation length. All yields were based on inside bark, merchantable cubic volumes net of a 30 cm stump, tops less than 10 cm in diameter and all trees less than 12.5 cm dbh. Surveys were simulated on each stand using two primary stocking assessments: the stocked quadrat system (Lowdermilk 1927) and the BC Ministry of Forests’ well-spaced system (Wyeth 1984, BC Ministry of Forests 2009). The stocked quadrat system utilizes relatively small sample plots to evaluate the presence of trees, with the percent of plots being occupied by at least one acceptable tree being the percent stocking for the stand. Commonly used plot sizes include 1 and 4 mil-acres in the imperial system (4.05 and 16.19 m2), and 4, 10 and 16 m2 in the metric system. The well-spaced system uses a 50 m2 circular plot, and tallies the maximum number of trees that can be separated by a specified minimum inter-tree distance (MITD), which is most commonly set at 2 m. The maximum number of trees that can be tallied at each plot is often capped at an m-value that roughly correlates to target stand densities (e.g. 6 well-spaced trees in a 50 m2 or 1/200th ha plot is associated with a target of 1200 well-spaced trees/ha). Seven variants of the quadrat system were employed, six using variations in plot size (4, 7, 10, 12.5, 16, 20 and 30 m2) and an eighth using 50 m2 plots divided into 12.5 m2 quadrants. In the latter case, termed Mean Stocked Quadrant or MSQ, each quadrant is treated as equivalent to a quadrat for field sampling, but must be summarized at the plot level for statistical purposes (the number of samples is based on the number of points or plots rather than the number of quadrants). Quadrat surveys were run using 1000 quadrats at randomly generated points (or 250 points where plots were split into 4 quadrants); it was assumed that these numbers were sufficient to get precise estimates of the expected population values. Nine variants of the well-spaced system were tested, using MITD values of 1.4, 2.0 and 2.6 m, and m-values of 4, 6 and unlimited. Each survey used 250 randomly located plots to estimate population parameters. The selection of well-spaced trees in each plot was simulated using multiple iterations of a weighted Markov Chain Monte Carlo (MCMC) process. For each iteration, a starting tree was randomly selected, and all trees within the MITD eliminated. The remaining trees were weighted by distance from the last selected tree, such that trees closer to the last selected tree had a higher probability of being selected to facilitate maximum “packing”. 46 Again, all remaining trees within the MITD were eliminated. This process was repeated until all trees had been exhausted. The number of well-spaced trees was recorded, and the entire evaluation repeated through 100 iterative steps in an attempt to find the maximum well-spaced value. In generating the stands and in completing growth and sampling simulations, all trees were considered to be acceptable crop trees. This means that the generation of simulated stands excluded trees that would not be considered as acceptable (diseased, deformed or otherwise not considered capable of reaching maturity). A necessary assumption was that excluded trees would have trivial effects on growth simulations and do not contribute to yields. Predictive models for exploring stocking effects on yield were generated using JMP version 8.0.1 (SAS Institute 2008). Both linear and nonlinear model forms were considered; ultimately the nonlinear Schnute equation (sensu Sit and Poulin-Costello 1994) was selected for its curve shape flexibility and the control it provided over asymptotes and/or extreme convergence points for a family of similar models. Results Graphs of the relationship between stocking measures and yield at age 80 are provided in Figs. 3-2 through 3-5. These graphs illustrate trends in curve shape and dispersion of points as stocking survey parameters are varied (plot size for quadrat surveys; MITD and m-value for well-spaced surveys). Not included in the illustrations are 12.5 m2 quadrats and the MSQ surveys, which have identical expected relationships to yield. The relationships to yield for these surveys can be inferred from the curve shape and dispersion trends for other quadrat plots sizes in Figs. 3-2 and 3-4. One of the first things that is apparent in many of the graphs in Figs 3-2 through 3-5 is the stratification of data points by categories of spatial pattern (see Fig. 3-1). In Fig. 3-2, for example, naturally regenerated redcedar stands assessed using 4 m2 quadrats have a distinctly different yield/stocking relationship than both categories of planted stands. In other cases (e.g. 4 m2 quadrats for pine in Fig. 3-4) it is the planted areas with random mortality that stand apart. It is evident that spatial distribution of trees plays an important role in stocking assessment outcomes and how they are interpreted. It is also evident that any models intended to predict yields using stocking values must be constructed with an awareness of the impacts of spatial pattern variability. A single model fit to most of the relationships in Figs. 3-2 through 3-5 would produce biased predictions for certain categories of stands. For example, models fit to all of the data points for 4, 7 and 10 m2 quadrats 47 Figure 3-2. Yield at 80 years for western redcedar, site index 28 m, as predicted by percent stocking values for six different plot sizes. Stands are stratified by spatial distribution, where E= planted stands with random mortality, F= planted stands with clustered and random mortality, and S = naturally regenerated stands. Figure 3-3. Yield at 80 years for western redcedar, site index 28 m, as predicted by well-spaced values for six combinations of MITD and m-value. Stands are stratified by spatial distribution, where E= planted stands with random mortality, F= planted stands with clustered and random mortality, and S = naturally regenerated stands. 48 Figure 3-5. Yield at 80 years for lodgepole pine, site index 20 m, as predicted by well-spaced values for six combinations of MITD and m-value. Stands are stratified by spatial distribution, where E= planted stands with random mortality, F= planted stands with clustered and random mortality, and S = naturally regenerated stands. Figure 3-4. Yield at 80 years for lodgepole pine, site index 20 m, as predicted by percent stocking values for six different plot sizes. Stands are stratified by spatial distribution, where E= planted stands with random mortality, F= planted stands with clustered and random mortality, and S = naturally regenerated stands. 49 respectively in Fig 3-2 would over-predict yields for the natural stands and under-predict yields for planted stands with random mortality. A likely explanation for these trends lies in the relative frequency of gap sizes (Fig. 3-6). When considering only quadrat sampling, assessments using smaller plots will detect small gaps or voids that surveys using larger plots will ignore. In looking at a uniform yield level (e.g. 250 m3/ha in the graphs for 4 and 7 m2 plots in Fig. 3-4), stocking levels are lower for the planted stands with random mortality than for the other two spatial categories. At the opposite end of the plot size spectrum, the 30 m2 plots in Fig. 3-4 illustrate a case where the planted stands with random mortality have a higher than average assessed stocking for a given level of yield. In the first case, small plots are overestimating unstocked area while in the second case large plots are underestimating unstocked area. A similar effect likely exists within the well-spaced survey system, but is more difficult to visualize. At least three possible solutions can be pursued to address this problem. The simplest is to find relationships for which the problem is minimized and the bias falls within acceptable limits. For the two species studied here, there appears to be a range of quadrat sizes in the 15 to 20 m2 range where the effects of spatial pattern are minimized; possibly with a smaller optimum plot size for lodgepole pine than for redcedar. While the optimum quadrat size (or combination of sample parameters for well-spaced surveys) may vary somewhat by species, it seems likely that Figure 3-6. Conceptual diagram of growing space utilization and its measurement in stocking surveys. Uniformly distributed stands (fine dashed line) have a high proportion of voids between trees in small gap sizes, making them particularly sensitive to survey plot parameters. Small plots applied to such stands will detect far too many gaps that have little or no effect on yield, whereas very large plots will ignore abundant small to medium gaps that cumulatively add to a significant impact on yield. The lower relative frequency of small gaps in clustered stands (coarse dashed line) will decrease the impact of such effects. 0 100 R el at iv e Fr eq u en cy % o f G ro w in g S p ac e U sa b le b y T re es Small Large Effective Gap Radius �� �� �� �� �� � �� �� �� � �� ��� � �� �� �� �� �� �� �� 50 there is sufficient similarity in curve shapes ( as in Figs 3-2 through 3-5) between species that a multi-species compromise can be found that sacrifices very little in terms of bias-by-spatial pattern. Two other approaches require quantification of spatial pattern in order to account for its effects. Figure 3-7a illustrates the stratification of sample variances by spatial pattern class for simulated stands of lodgepole pine surveyed using a well-spaced survey (MITD = 2.0 m; m- value = 6). Based on this chart, it appears that sample variance may be a useful indicator of spatial pattern. If the relationship between well-spaced stocking and yield for lodgepole pine is fit (Fig. 3-7b) using a variant of Schnute’s equation (sensu Sit and Poulin-Costello 1994): (1) where: Y80 = Stand yield at age 80 in m3/ha WS = mean number of well-spaced trees (MITD = 2.0; m-value = 6) WSmax = m-value b 0 and b 1 are shape parameters b 2 is the predicted maximum yield for a fully stocked stand a pattern of residuals is obtained as shown in Fig 3-7c. It is apparent that the results are biased by spatial pattern category: the mean residual value for the planted stands with random mortality is significantly different from zero on the positive side (p < 0.0001), while the mean residuals for the planted stands with clustered mortality and the natural stands are significantly different from zero on the negative side (p = 0.0013 and p = 0.0187 respectively). Adding the sampling variance to Eq. 1 as a linear adjustment: (2) where: Var = sampling variance results in a reduction in RMSE from 9.43 to 7.08 m3/ha and largely eliminates the problems observed in the earlier residuals (Fig. 3-7d). However, it is important to note that not every case worked this well; attempts to make similar corrections with other combinations of species and surveys produced mixed results (well-spaced and MSQ variants only; the binomial variance of quadrat surveys is unsuitable for this purpose). In all cases there was improvement in the RMSE values, but outcomes related to the correction of bias by spatial pattern class were highly 51 variable. It is possible that a better adjustment methodology exists but has yet to be discovered. Other spatial indices that may be useful and which are compatible with quadrat sampling as recommended by Pretzsch (2009) are those proposed by Clapham (1936) and Morisita (1959), although both would require the addition of a total tree count in each quadrat. A final option for addressing variation in predicted yields by spatial pattern class is a separate set of model fits for each class. Such an approach eliminates the issue of finding an appropriately consistent adjustment mechanism, but introduces a separate problem. While the spatial pattern classes used in this study make it appear that there are distinctly different relationships for identifiable classes, this is purely an artifact of the study design. In reality, there is far more variability in spatial patterns than employed here; real stands do not so readily fall into categories such as “planted stands with random mortality”. Great effort would be required to ensure that the spatial pattern of each sampled stand could be appropriately categorized without unduly adding to the cost of regeneration surveys. Figure 3-7. The sampling variance for measures of lodgepole pine stocking (a) appears to be a good indicator of differences in spatial pattern (E= planted stands with random mortality, F= planted stands with clustered and random mortality, S= naturally regenerated stands), and may be useful in accounting for some of the variation in relationships predicting yield from regeneration survey statistics (b). Residuals for the illustrated fit (c), when stratified by spatial distribution class, are biased by spatial pattern class and show patterns of nonlinearity. The addition of sample variances from (a) for the independent variable, which are also stratified by spatial pattern, improved the fit and corrected both the bias and nonlinearity problems (d). 0 1 2 3 4 5 0 1 2 3 4 5 6 -30 -20 -10 0 10 20 30 0 1 2 3 4 5 6 -30 -20 -10 0 10 20 30 0 1 2 3 4 5 6 0 100 200 300 400 0 1 2 3 4 5 6 Y ie ld a t A g e 8 0 ( m 3 /h a) Sa m p lin g V ar ia n ce Y 8 0 R es id u al s: U n ad ju st ed F it Y 8 0 R es id u al s: A d ju st ed F it Mean Number of Well Spaced Trees Mean Number of Well Spaced Trees a c b d 52 Returning to the general nature of the relationships illustrated in Figs 3-2 through 3-5, it is `apparent that the shape of the relationships varies considerably by survey parameters (plot size for quadrat surveys and MITD for well-spaced surveys), with corresponding impacts on the predictive utility. Some predictive relationships have broad ranges for the independent variable where there are only small changes in the dependent variables (shallow slopes), paired with other regions where the opposite case (steep slopes) is prevalent. Extreme examples include the 4 m2 quadrat relationships in both Figs. 3-2 and 3-4, and the combination of a 1.4 m MITD with no m-value in Figs. 3-3 and 3-5. From the perspective of using these relationships as predictive models, this situation is far less desirable than a relationship that is close-to-linear with little vertical dispersion over any given (small) range of the independent variable. Examples include those illustrated for 16 m2 quadrats in Figs. 3-2 and 3-4, and the combination of a 1.4 or 2.0 m MITD and an m-value of 6 as shown in Fig. 3-5. However, minimal curvature is not always a definitive characteristic for a good predictive relationship. A case in point is the addition of an m- value of 4 for the redcedar well-spaced surveys in Fig 3-3 (top row versus bottom row of charts). In this example, adding the m-value primarily compresses the x-axis for stands with higher yields (note change in x-axis scaling), without a change in vertical dispersion (prediction errors). A related issue exists for some relationships where a large percentage of stands are assessed as having a stocking level at or near the maximum possible value. Examples of this condition can be found for the 20 m2 and particularly 30 m2 quadrats in Figs. 3-2 & 3-4, and for the well- spaced values capped with an m-value of 4 in Fig. 3-3 (the same trend appeared for lodgepole pine with an m-value of 4, and is less prevalent or absent for redcedar with an m-value of 6). In these cases, predictive sensitivity is lost at the upper range of volumes where small changes in the value of survey summary statistics again result in relatively large changes in yield. For the case of quadrat surveys, this change results from the use of large plots, where fewer and potentially more widely dispersed trees can still result in high values for percent stocking. In the case of the well-spaced surveys in Fig. 3-3, it appears that the relatively low m-value of 4 allows plots where the growing space is not fully utilized to be erroneously assigned a maximum value. In both cases there are stands that are not achieving maximum volume growth, but are assessed as having full stocking. To this point, all relationships between stocking measures and yield have been evaluated using the single reference age of 80 years for yield. Comparison across ages is difficult when yield is assessed using absolute values, but is more reasonable using relative values (Fig 3-8). To determine these relationships, the predicted yield at full stocking was first determined using variants of Eq. 1 for each species and level of Yz where z = assessment age for yield. In each case, the determined value of the parameter b 2 was set as the maximum yield, and the simulated 53 yield for each stand was divided into the value of b 2 to calculate the relative yield. Models to predict relative yield at each assessment age and for each species, this time based on quadrat stocking, could then be fit using a variant of the Schnute equation (sensu Sit and Poulin-Costello 1994): (3) where: RYz = Relative stand yield at reference age z in m3/ha %Stocking = percent of plots assessed as stocked using 16 m2 quadrats b 0 and b 1 are shape parameters Parameters and fit statistics for the outcomes can be found in Table 3-3. Better fits were found for older ages than for younger ones, as indicated by the decreasing RMSE values with increasing age. The somewhat weaker correlations between stocking and yield at the younger ages were primarily affected by spatial pattern: stratification of data points by spatial pattern class was much more distinct for the younger ages resulting in greater dispersion of data points. The relationships presented in Fig. 3-8 can also be adapted, with one basic assumption, for use across a range of site indices. Applying Eichhorn’s rule (Eichhorn 1904 as cited in Assman 1970), stand yield can be associated with top height and density rather than age and density, thus accounting for the effects of site quality. Following this principle, the ages from Fig. 3-8 have been converted to top height, based on site index curves by Kurucz for western redcedar Figure 3-8. The relationship between stocking and relative yield varies by age (top height) of assessment and species. Intermediate labels for lodgepole pine have been left out for clarity, but follow the same order as for western redcedar. 0.0 0.2 0.4 0.6 0.8 1.0 0.0% 20.0% 40.0% 60.0% 80.0% 100.0% 16 m2 Quadrat Stocking 0.0 0.2 0.4 0.6 0.8 1.0 R el at iv e Y ie ld 0.0% 20.0% 40.0% 60.0% 80.0% 100.0% 16 m2 Quadrat Stocking �������� �������� �������� �������� ��������� �������� ��������� Western Redcedar Lodgepole Pine ���� ������ � ������ ��� � ��� ���� � ��� ���� � ��� ���� � ��� ���� � ���� ���� ���� ������ � ������ ��� � ��� ���� � ��� ���� � ��� ���� � ��� ���� � ���� ���� 54 (Mitchell and Polsson 1988) and J.S. Thrower and Assoc. (1994) for lodgepole pine as expressed in the TASS output tables. This allows the relative yield curves to be applied independent of site using a reference height in place of a reference age.In evaluating the changes in the yield/ stocking relationships over time in Fig. 3-8, it appears at first glance that there is a much greater degree of variation for western redcedar than for lodgepole pine. However, this perception is skewed by differences in growth rate. The redcedar was intentionally modeled with a higher site index (28 m @ 50 yrs breast height age) than the pine (20 m @ 50 yrs breast height age). If the vertical distance between fit lines for 40 and 50 yrs for redcedar (4.3 m increment in top height) are contrasted with that for 60 to 100 yrs for lodgepole pine (4.6 m increment), the line spacing is similar. If a stocking level of 50% is assumed for both species, then application of Eq. 3 along with parameters from Table 3-3 would result in relative yield differences of 0.052 for redcedar and 0.056 for pine. It appears, then, that similar increases in height between these species result in similar increases in relative yield at a given level of stocking. The difference in curve shape between the two species likely results from a combination of two factors. For a given level of stocking, redcedar as modeled in this study is achieving a much higher relative yield than is the pine. A large part of this can be attributed to the faster growth rates of redcedar, owing at least in part to the higher site indices used. When compared on an equivalent height basis, redcedar still achieves a greater relative yield, although the margin is much smaller than if the comparison was based on age. Much of the remaining difference can likely be attributed to the wider crowns of the redcedar, which are able to increase site occupancy for a given stand top height and stand density. The higher shade tolerance of redcedar may also be a factor, allowing it to develop deeper crowns that capture a higher percentage of incoming sunlight. Table 3-3. Fit statistics for equations predicting relative yield at each of 5 reference ages (z), along with maximum yield (Yz) at each age used to calculate relative yield. Western Redcedar Lodgepole Pine z Max Yz b 0 b 1 RMSE Max Yz b 0 b 1 RMSE 40 287.4 0.3894 0.9646 0.089 118.6 -0.1849 0.8489 0.085 50 471.8 0.5503 1.064 0.054 186.1 -0.1751 0.8804 0.048 60 648.3 0.7725 1.1424 0.040 245.8 -0.0607 0.8864 0.033 80 966.8 1.0511 1.2989 0.029 336.5 0.2116 0.8761 0.021 100 1280.9 1.1471 1.4662 0.027 398.2 0.4859 0.8509 0.017 55 Discussion While common definitions of stocking appear as relatively simple concepts, our collective understanding of the concept is less clear. If stocking is defined as the relative ability of a population of trees to access all of the available resources to facilitate growth, then our understanding of stocking is crude at best. A thorough comprehension of stocking would require an integrated awareness of light interception at various sun angles, the spread of root systems to access water and nutrients, the interaction between root and crown systems of adjacent trees, and the influence of climatic and edaphic variables on all of these factors. In lieu of such omniscience, the concept of stocking is simplified to something comparable to the vertical projection of tree crowns onto the planar ground surface. This concept is most easily visualized using stocked quadrats as a measure of stocking. For measurement purposes, a randomly or systematically located quadrat plot is stocked if it contains an acceptable tree and unstocked otherwise. For visualizing the “stocked” area of a stand, the plots are instead centered on each tree. The net area of the overlapping plots projected onto the ground (the percentage of the total area covered by the resultant “shadow”) is the expected value of the stocking survey (pers. comm. E. McWilliams). The total area covered will vary with plot size. If an area devoid of trees is imagined, characterization of that gap will vary with the encroachment of different sized quadrat areas into the gap edges as the amount of resources in the gap that are accessed by adjacent trees is quantified. With the utilization of larger quadrats the gap, in essence, shrinks. The characterization of unstocked area using the well-spaced concept can be viewed in a somewhat similar manner. Within a given plot, a circular shadow projected onto the ground from each well-spaced tree can be imagined using the MITD as the radius, with the sum of the shadowed area representing the stocking proportion. In contrast to the quadrat method, trees that are not well-spaced are not considered to cast their shadows, nor are trees outside of the plot. This characterization is imperfect in that (i) circular shadows are not permitted to overlap, leaving areas between circles that cannot be shadowed even if occupied by a non well-spaced tree, (ii) some of the area shadowed by well-spaced trees is outside of the plot and trees outside the plot cannot cast shadows inside the plot, and (iii) m-values (if applied) limit the number of trees within a plot that can throw shadows. The simpler, more consistent definition of growing space provided by the quadrat surveys as compared to the well-spaced method likely contributes to the better correlations between stocking and yield that are evident for the quadrat method in Figs. 3-2 through 3-5. However, this is not to say that the well-spaced survey is inferior. The intent of the well-spaced survey goes beyond simply stocking; in addition it attempts to characterize the desired stand conditions 56 within the portion of the area that is stocked (Wyeth 1984). This philosophy is carried even further in a regeneration survey methodology for Ontario (White et al. 2005) where 16 m2 plots are used to assess both percent stocking and well-spaced trees (MITD = 1.2 m, m-value = 4). In this case, the well-spaced assessment is intended to be used only to assess stand conditions within the “stocked” area. Regardless of the survey method employed, the characterization of growing space as the vertical projection of a current or future crown is flawed. In reality, there is no fixed boundary between “utilized” and “unutilized” portions of the resources available to a stand. The roots of a tree will frequently extend far beyond the limits of the crown in search of resources (e.g. Göttlicher et al. 2008), and the ability of a tree to capture extra sunlight on the edge of a gap will vary with its position around the gap relative to the sun direction and vertical angle. Physical space may be wholly accessed if nearby a tree, but only partially accessed if further away. Sunlight within gaps, particularly those on the smaller end of the size spectrum, will be increasingly captured as the stand grows taller and crowns on the gap edge become longer, an effect that is accentuated with lower sun angles at higher latitudes. This effect is likely a major cause of the change in curve shape with increases in height (age) in Fig 3-8. The vertical or cylindrical projection of tree canopies as an expression of growing space is further flawed in cases where light is not a primary factor limiting growth. On arid sites where water becomes limiting, as is often the case for ponderosa pine (Pinus ponderosa P. Laws ex C. Laws.) ecosystems in western North America, stands develop where canopies of adjacent trees are often far apart, and full stocking might be characterized where growing space is better represented using a conical projection (the “shadow” area is considerably larger than the cross section of the crown). In this case, the height of the cone would be based on the expected stand top or mean height at a reference age (analogous to the expected crown width for an open grown tree at a reference age for the cylindrical projection). Alternately, since the cylindrical and conical projections will be mathematically linked (although with some bias if the cone base becomes too large and shadows excessively overlap), it might be acceptable to describe this case such that a relative yield of 1.0 is achieved at a stocking level less than 100%. Despite the imperfect characterization of growing space provided by measures of stocking, it appears that reasonable predictive relationships exist between stocking and yield. Two major applications are evident. The first of these is the traditional stocking assessment, where a minimum acceptability threshold is set for site occupancy by regenerating trees following harvest operations. The explicit links to yield demonstrated in this study and others (e.g. Martin et al. 2002) will, with further development, allow for yield-based definitions of these thresholds. If for example, an expectation exists that the minimum acceptable relative yield threshold is 0.85 57 at 80 years (all things other than stocking being equal), relationships such as those illustrated in Figure 3-8 can be used to determine appropriate stocking thresholds. For the site indices used in this study (20 m for lodgepole pine and 28 m for western redcedar), the appropriate stocking thresholds would be 72% and 86% respectively using 16 m2 stocked quadrat plots. The second major (but closely related) application is for adjustments to yield predictions made using either yield tables or non-spatially explicit growth models. The nature of these yield prediction methods makes them unsuitable on their own for dealing with unstocked areas. An individual tree, distance-independent model of forest growth inherently assumes a uniform spatial distribution of all trees in a stand represented by a single sampled tree. This type of model is poorly suited to dealing with unstocked areas, although the scale of voids is a factor (Bruce 1977). In many cases, such models are seriously criticized by operational foresters for being too optimistic in their predictions, at least partly because of differences between the uniform stocking assumption of the model and the non-uniform stocking conditions found in forested landscapes. Such problems are exacerbated when foresters’ perceptions are based on experience primarily from (older) naturally regenerated stands where medium to large stocking gaps may be common (pers. comm. M. Bokalo, Mixedwood Growth Model Development Team, University of Alberta). A stocking correction, or OAF, is essential in these cases to calibrate yield predictions to actual stocking conditions observed in the field. In correlating stocking to yield for the purpose of applying OAF’s, the issues are slightly different than was the primary focus of the current study. Firstly, the age of assessment will not always be juvenile stands, and stocking will change over time (Yang et al. 2008). However, such changes will be small up to approximately 50 to 60% of culmination age (Doucet, 1992, Yang et al. 2008). Of greater concern is the split between small stocking gaps that might be caused by changes in stand density and large stocking gaps that cannot ever be occupied by the existing trees. Variable density yield tables and individual tree, distance-independent growth models account for at least a major portion of the first category, but not the second. The need to account for only the larger voids is likely the reason for the relatively large quadrats (22.9 and 40.7 m2) used in the British Columbia OAF surveys (BC Ministry of Forests 1998a,b). Both regeneration stocking assessments and yield prediction adjustments using OAF’s require well defined and consistent application of stocking measurements in order to extract full value from monitoring systems. Hosie (1953) decried inconsistencies in measurement and reporting of stocking, a problem that continues in many jurisdictions to the present (Lutze et al. 2004). In Alberta, recent government initiatives are pushing forest licensees to develop explicit linkages between regeneration stocking and yield, with corresponding rigour in measurement and reporting systems. By contrast in neighboring British Columbia, current policy is allowing 58 an increase in the diversity of regeneration monitoring systems, with little consistency attached to the definition of objectives, measurement parameters and/or acceptability thresholds. From the perspective of stocking effects on yield, I argue that the former approach will pay large dividends in the long term based on increased knowledge and certainty regarding forest stand inventories and the modeling of forest estate goods and services. Conclusions For the species and spatial patterns in this study, there is a strong relationship between stocking (or growing space occupancy) and yield. The abstract concept of stocking is simplified into a measurable quantity typically using the physical occupation of 2-dimensional space, most often using assumptions about maximum potential crown expansion. Such practices are formalized in the design of regeneration stocking surveys, which in turn are linked to stocking or regeneration standards. This study largely assumes that the actual yield of a stand relative to the maximum achievable yield for a “fully stocked” stand is a useful criterion for assessing reforestation success. Based on this assumption, it is concluded that several factors must be considered that have not been widely recognized previously. First, for a given relative yield threshold, the required stocking will vary by stand height (or age) with higher stocking levels required for shorter rotations. Second, depending on the definition of stocking as expressed in the measurement parameters, different stocking levels may also be required based on different spatial patterns although future studies may find ways of eliminating this problem using sampling variance or other measures of dispersion as predictive parameters. Finally, this study also confirms that species-specific thresholds are important. The selection of sampling parameters is important to simplify the process of assessing stocking based on yield objectives, although there is no absolutely “correct” choice. For the two species studied here, the widely used 16 m2 quadrat or its 4 mil-acre equivalent appear to be suitable. Several previous authors have commented on the importance of consistency in measurement and reporting of stocking values, and these sentiments are echoed here. 59 References Alberta Reforestation Standards Science Council. 2001. Linking regeneration standards to growth and yield and forest management objectives. Alberta Ministry of Sustainable Resource Development, Edmonton AB. 52 p. Armson, K.A. 2005. Regeneration standards: what has the past to show us? For. Chron. 81: 781- 784. Assman, E. 1970. The principles of forest yield study: studies in the organic production, structure, increment and yield of forest stands. Pergamon Press, Oxford UK. BC Ministry of Forests. 2009. Silviculture survey procedures manual: stocking and free growing survey. BC Ministry of Forests, Forest Practices Branch, Victoria BC. 219 p. BC Ministry of Forests. 1998. Ground-based survey method. OAF1 Project Report, BC Ministry of Forests, Victoria BC. 19 p. BC Ministry of Forests. 1998. An overview of stocking gaps and OAF1 estimates for TIPSY. OAF1 Project Report, BC Ministry of Forests, Victoria BC. 7 p. Bella, I.E. 1976. Assessment of regeneration stocking standards used in Alberta. Canadian Forestry Service, Northern Forest Research Centre, Edmonton AB. NOR-X-167. 38 p. Boyer, W.D. 1977. Stocking percent and seedlings per acre in naturally established longleaf pine. J. For. 75: 504-506. Braathe, P. 1982. Stocking control and its effect on yield. Forest Industry Lecture Series No. 10, University of Alberta, Edmonton AB. 21 p. Braathe, P. 1976. Investigations concerning the development of Norway spruce regeneration which is irregularly spaced and of varying density: 2. the stability of the zero-square percentage. Meddelelser fra Norsk Institutt for Skogforskning (Reports of the Norwegian Forest Research Institute) 32.15: 507-520. Bruce, D. 1977. Yield differences between research plots and managed forests. J. For. 75: 14-17. Clapham A.R. 1936. Over-dispersion in grassland communities and the use of statistical methods in plant ecology. J. Ecol. 24: 232–251 Cowlin, R.W. 1932. Sampling Douglas fir reproduction stands by the stocked-quadrat method. J. For. 30: 437-439. 60 Doucet, R. 1991. The influence of stocking of regeneration on the yield of naturally regenerated jack pine and black spruce stands. In Proceedings of the conference on natural regeneration management, March 27-28, 1990, Fredericton NB. Edited by C. M. Simpson. Forestry Canada, Maritimes Region, Fredericton NB. pp. 181-192. Eichhorn, F. 1904. Beziehungen zwischen Bestandshöhe und Bestandsmasse. Allgemeine Forst- und Jagdzeitung 80: 45-49. Göttlicher, S.G., Taylor, A.F.S., Grip H., Betson, N.R., Valinger, E., Högberg M.N., and Högberg, P. 2008. The lateral spread of tree root systems in boreal forests: estimates based on 15N uptake and distribution of sporocarps of ectomycorrhizal fungi. For. Ecol. Man. 255: 75- 81. Goudie, J.W. 1998. Model validation: a search for the magic grove or the magic model. In Stand Density Management: Planning and Implementation, Nov. 6-7 1997, Edmonton AB. Edited by C. Bamsey. Clear Lake Ltd., Edmonton AB. pp. 45-58. Grant, J.A.C. 1951. The relationship between stocking and size of quadrat. Univ. Toronto For. Bull. No. 1. 35 p. Haig, I.T. 1931. The stocked quadrat method of sampling reproduction stands. J. For. 29: 747- 749. Hosie, R.C. 1953. Forest regeneration in Ontario based on a review of surveys conducted in the province during the period 1918-1951. Univ. Toronto For. Bull. No. 2. 134 p. Huang, S., Meng, S.X., and Yang, Y. 2009. A growth and yield projection system (GYPSY) for natural and post-harvest stands in Alberta. Alberta Sustainable Resource Development, Tech. Rep. T/216, Edmonton AB. 22 p. J.S. Thrower and Assoc. 1994. Revised height-age curves for lodgepole pine and interior spruce in British Columbia. Contract report to BC Ministry of Forests, Victoria BC. 27 p. Lowdermilk, W.C. 1927. A method for rapid surveys of vegetation. J. For. 25: 181-185. Lutze, M., Ades, P., and Campbell, R. 2004. Review of measures of site occupancy by regeneration. Austral.. For. 67: 164-171. Krajina, V.J., Klinka, K. and Worrall, J. 1982. Distribution and ecological characteristics of trees and shrubs of British Columbia. University of British Columbia, Faculty of Forestry, Vancouver BC. 131 p. 61 McWilliams, J. and McWilliams, E. 2009. A review and analysis of the effect of BC’s current stocking standards on forest stewardship. Contract report to the Association of BC Forest Professionals, Vancouver BC. 32 p. Martin, P.J., Browne-Clayton, S., and McWilliams, E. 2002. A results-based system for regulating reforestation obligations. For. Chron. 78: 492-498. Martin, P.J., Browne-Clayton, S., Day, K., and Taylor, G. 2005. Improving regeneration performance standards: comments based on early experience with three new approaches in British Columbia. In The thin green line: a symposium on the state-of-the-art in reforestation. July 26 to 28, 2005, Thunder Bay, ON. Edited by S. J. Colombo. Ontario Ministry of Natural Resources, Sault Ste. Marie, ON pp. 62-68. Mitchell, K.J. 1975. Dynamics and simulated yield of Douglas-Fir. Forest Science Monographs 17. 39 p. Mitchell, K.J. 1969. Simulation of the growth of even-aged stands of white spruce. Bull. Yale Sch. For. No. 75. 48 p. Mitchell, K.J., and Polsson, K. 1988. Site index curves and tables for British Columbia: Coastal species. BC Min. For., Res. Br., FRDA Rep. 37. 29 p. Morisita, M. 1959. Measuring of the dispersion of individuals and analysis of the distributional patterns. Memoirs of the Fac. Sci., Kyushu Univ., Series E (Biology) 2: 215–235 Newton, P.F. 1998. An integrated approach to deriving site-specific black spruce regeneration standards by management objective. For. Ecol. Manage. 102: 143-156. Ontario Ministry of Natural Resources. 2001. Silvicultural effectiveness monitoring manual for Ontario. Ontario Mininistry of Natural Resources, Natural Resources Information Centre, Peterborough ON. 42 p. Ponge, J.-F. 2005. Emergent properties from organisms to ecosystems: towards a realistic approach. Biol. Rev. 80: 403-411. Pound, R., and Clements, F.E. 1898. A method of determining the abundance of secondary species. Minn. Bot. Studies 2: 19-24. Pretzsch, H. 1997. Analysis and modeling of spatial stand structures. Methodological considerations based on mixed beech-larch stands in Lower Saxony. For. Ecol. Manage. 97: 237-253. Pretzsch, H., Biber, P., and Dursky, J. 2002. The single tree-based stand simulator SILVA: construction, application and evaluation. For. Ecol. Manage. 162: 3-21. 62 Pretzsch, H. 2009. Forest dynamics, growth and yield: from measurement to model. Springer- Verlag, Berlin. 664 p. SAS Institute. 2008. JMP user guide, release 8. SAS Institure Inc, Cary NC. Sit, V., and Poulin-Costello, M. 1994. Catalogue of curves for curve fitting. BC Ministry of Forests, Victoria BC. Biometrics Information Handbook 4. 110 p. Staebler, G.R. 1949. Predicting the volume and normality of reproduction stands of Douglas-fir. J. For. 47: 828-833. Wellner, C.A. 1940. Relationships between three measures of stocking in natural reproduction of the western white pine type. J. For. 38: 636-638. White, R.G., Bowling, C.L., Parton, W.J., and Towill, W.D. 2005. Well-spaced free-growing regeneration assessment procedure for Ontario. Ontario Mininistry of Natural Resources, Northwest Science and Information. NWSI Technical Manual TM-007. 39 p. Wyeth, M.H. 1984. British Columbia Ministry of Forests regeneration survey system. In New forests for a changing world : proceedings of the 1983 Convention of the Society of American Foresters, October 16-20, 1983, Portland, Oregon. Society of American Foresters, Bethesda, Md. pp. 149-152. Yang, Y., and Huang, S. 2008. Modeling percent stocking changes for lodgepole pine stands in Alberta. Can. J. For. Res. 38: 1042-1052. Yang, Y., Huang, S., and Dempster, W.R. 2008. Percent stocking models for four major Alberta tree species. Forestry 81: 599-615. 63 Chapter 4. Simplified Models to Predict Future Mixedwood Outcomes from Regeneration Survey Summaries7 Introduction Since early in the 20th century, forest regeneration surveys have been employed in assessing reforestation success (see literature review in Chapter 1). Various strategies have been employed to define success, but linking survey results to long term desired outcomes has always been the ultimate goal. With the advent of electronic computers in the middle of the 20th century, increasingly sophisticated simulation techniques have gradually resulted in improvements to these linkages. A number of studies have been undertaken that involve simulated sampling of forest stands to test the efficiency and effectiveness of various sampling methodologies. Many of these (e.g. Mackisack and Wood 1990, Magnussen 1999, Carter 2007) are focused on forest inventory methods and in particular angle-count sampling. Only in a few isolated cases have simulation studies been undertaken to investigate regeneration surveys. MacLeod (1977) evaluated various implementations of a stocked quadrat survey on four simulated populations with two levels of stand density by two levels of clumping. Quadrats were arranged systematically either in groups or individually, and evaluated based on the sampling intensity required to achieve a specified standard error. It was concluded that clusters of plots resulted in the need for a much greater number of quadrats; however, from an operational efficiency perspective this might be offset by reduced plot-to-plot travel time. Kaltenberg (1980) implemented a series of Monte Carlo simulations of regeneration surveys following eight different sampling rules including clustered and unclustered quadrats of varying size using both stocked and list methods (clustered quadrats were in reality quadrants of a larger plot), location to nearest tree, and variable radius by height. Each method was evaluated for accuracy, precision and operational efficiency with regard to estimating stocking percentage and overall stems/acre (SPA). Nearest tree methods were generally the most efficient. The nearest tree and variable radius methods produced biased estimates of SPA if the underlying distributions were unknown, otherwise all methods were unbiased. Precision varied by sampling intensity and method. Martin et al. (2002) developed a survey system linked to yield predictions for even-aged stands of interior spruce (a common hybrid of Picea glauca [Moench] Voss and P. engelmanii 7 A version of this chapter will be submitted for publication. Farnden, C. 2010. Simplified models to predict future mixedwood outcomes from regeneration survey summaries. 64 Parry) and lodgepole pine (Pinus contorta Dougl. ex Loud var. latifolia Engelm.) in British Columbia. They used the mean number of stocked quadrants (MSQ) in a 50 m2 circular plot to predict yield at age 80. This work has been developed through several unpublished reports (J.S. Thrower & Assoc. 2002, 2003) into an operational regeneration survey system that is applied in several areas of BC, and has provided considerable philosophical influence on the early stages of this thesis. The purpose of this chapter is to explore the ability of various regeneration survey systems to provide summary statistics that are useful for predicting future yields in a boreal mixedwood forest. Specifically, the study will: 1. contrast the summary statistics of three survey systems (quadrat, 5-5-2-2 and mixedwood well-spaced) for their utility as parameters for predicting yields by species: white spruce (Picea glauca [Moench] Voss) and trembling aspen (Populus tremuloides Michx.), 2. explore the impact of plot size for quadrat surveys on the utility of quadrat summary statistics for predicting yields by species group, and 3. contrast simplified model formulations for predicting yields by species group, particularly with regard to the inclusion of variation by site index. Methods This study develops a broad set of predictive models for total yield (m3/ha) at age 80 (Y80) and species composition (the proportion of that yield as spruce volume - Yc) in boreal spruce- aspen mixedwoods. Models were built based on (i) generation of a large number of simulated stand conditions representing a wide range of stocking conditions and species mixtures, (ii) predicting yields for those stands using a widely accepted individual tree growth model, (iii) simulating regeneration surveys on those stands and compiling summary statistics and (iv) constructing simplified models that will mimic the yield predictions of the individual tree model using the regeneration survey summary statistics as predictor variables. Stand Generation Using the stand generation simulator described in Chapter 2, an initial set of 108 stands was generated representing three levels of planted spruce stocking (1200 trees/ha with 10%, 40% and 70% random mortality), nine levels of area coverage by aspen clones (10 to 90% coverage in increments of 10%), and four levels of aspen density for areas not covered by clonal clumps (0, 50, 200 and 800 trees/ha). For the area covered by aspen clones, stand density ranged randomly by stand from 10,000 to 16,000 trees/ha. Each stand covered 9 ha. 65 To augment the initial populations, “extra” stands were developed to cover conditions representative of very low spruce stocking. The five highest levels of aspen clonal coverage (50 to 90%) were used coupled with 800 trees/ha of aspen in the inter-clonal areas and a random spruce density ranging from 0 to 200 trees/ha. Size distributions were developed for spruce and aspen based on 13-year-old stands generated using the Mixedwood Growth Model (MGM) version 2007A (University of Alberta 2008)8 tempered with observed values from a mixedwood field experiment (Kabzems et al. 2007). Aspen heights were assumed to follow a normal distribution with a mean of 3.5 m and a standard deviation of 0.5 m. Spruce heights before adjustments for suppression were assumed to follow a normal distribution with a mean of 2.0 m and a standard deviation of 0.4 m. Two of the sub-models in the stand generator required fitting to data specific to the stand types being generated: the height adjustment model for overtopping competition, and the diameter prediction models for white spruce and trembling aspen. Data for calibrating these functions was obtained from the large stem mapped plots described in Chapter 2. The sub- models were fit using JMP version 8.0.1 software (SAS Institute 2008). Growth Simulations Future outcomes for each simulated stand were predicted using MGM. The distance- independent nature of this model necessitated discrete simulations using sample plots from each stand to capture spatial variability in stocking. Thirty-six samples were selected using 5.64 m radius plots (100 m2) centered on randomly selected points. Each sample treelist was simulated independently in MGM, with yield tables for the 36 samples in a stand blended to generate a composite yield table. Stands had initially been generated assuming a site index (breast height age 50) of 20 m for both spruce and aspen. Where other site indices were required for growth simulations, adjustments to the size of all trees were employed, assuming that the needed size adjustment would be the same as the ratio of site indices: (1) This adjustment factor was applied equally to both heights and diameters, leading to small shifts in allometry given that diameter reaches zero when tree height is 1.3 m. Several MGM simulations were conducted to test the effects of these shifts, and the differences were found to be negligible. FactorAdjustmentSize = SI 50 20 8 The current version of the model, as described in the online documentation, is MGM 2008A. Changes to the model from version 2007A to 2008A affect only the user interface and bug fixes – no differences should exist in model outcomes 66 Simulated Surveys It was necessary to find survey summary statistics that could be used as predictive parameters to forecast future mixedwood yields. In order to achieve this, the summary statistics would need to adequately quantify the effects of growing space allocation to the conifer and broadleaved stand components. Three survey methods were tested as described below. Quadrat Surveys Quadrat surveys typically classify individual plots as stocked or not stocked based on specified acceptability criteria for crop trees. For most forest types and most applications of the system, a single acceptable tree in a plot is sufficient to classify the plot as stocked. In the boreal mixedwood variant of this survey, plots are classified into four categories: unstocked (0), stocked with at least one acceptable broadleaved tree (1), stocked with at least one acceptable conifer tree (2) or stocked with at least one of each (3). Acceptability standards for individual trees are typically based on size and health criteria; for the purpose of simulation all trees were assumed to be acceptable. Six variants of the quadrat survey were simulated, each with a different plot size ranging from 5 to 20 m2 in 2. 5 m2 increments (Q5, Q7.5, Q10, Q12.5, Q15, Q17.5 and Q20). Each simulated survey used 144 quadrats. Summary statistics consisted of the percentage of plots falling into each of the four stocking classes (C 0 , C 1 , C 2 and C 3 ). 5-5-2-2 Surveys The 5-5-2-2 survey (BC Ministry of Forests 2009) is a sampling protocol currently proposed for use in boreal mixedwoods. This protocol classifies cardinally-oriented quadrants (as opposed to quadrats) of a 50 m2 circular plot in much the same manner as quadrats are classified as described above (C 0 though C 3 ), but with strict acceptability standards for conifer trees based on proximity to competing broadleaved trees. As implemented in simulations here, each conifer tree is assessed using a tree-centered sub-plot with a 5 m radius. Sub-plots are again divided into quadrants, with quadrant orientation varied such that the likelihood of the tree being accepted is maximized. Acceptable conifers are those in which the distance to the nearest broadleaved tree is at least 2 m in all directions, and the sum of the distances to the nearest broadleaved tree in all four quadrants is at least 14 m. This is a slightly different variant than is described in Ministry of Forests (2009) whereby the nearest deciduous tree must be at least 5 m away in two adjacent quadrants (5 m + 5 m + 2 m + 2 m = 14 m). For the purpose of simulation, the orientation of sub-plot quadrants was varied through 90º in 1º increments to find the optimum orientation. Assuming that a quadrant is roughly analogous to a quadrat for sampling purposes, 36 plots were simulated for each survey (one quarter the number of plots as was used for the quadrat surveys). Summary statistics were as for quadrat surveys: the percentage of quadrants falling within each of the four stocking classes (C 0 , C 1 , C 2 and C 3 ). Mixedwood Well-Spaced Surveys The well-spaced concept, widely used in regeneration surveys in British Columbia, is one whereby the maximum number of acceptable trees is tallied that can be separated by a minimum inter-tree distance (MITD) within a 50 m2 circular plot (BC Ministry of Forests 2009). The intent of this methodology is to assess occupancy of growing space. For this study, an adaptation of this procedure has been developed for species-stratified boreal mixedwood stands, wherein three independent assessments of well-spaced trees are carried out in each plot using an MITD of 2.5 m. The first well-spaced tally considers only broadleaved trees, the second considers only conifer trees, and the third considers all trees regardless of species. Summary statistics consist of the mean number of well-spaced broadleaved trees (WS 1 ), the mean number of well-spaced conifer trees (WS 2 ), the mean number of total trees (WS 3 ), and the mean number of broadleaved trees by plot for each tallied well-spaced conifer (WS r ) as a measure of overtopping competition. In each plot tally, trees are selected such that the number of well-spaced trees is maximized. For simulation, this was accomplished through multiple iterations of a weighted Markov Chain Monte Carlo (MCMC) process. For each iteration, a starting tree was randomly selected, and all trees within the MITD were eliminated. The remaining trees were weighted by distance from the last selected tree, such that trees closer to the last selected tree had a higher probability of being selected as the next well-spaced tree to facilitate maximum ”packing”. Again, all remaining trees within the MITD were eliminated. This process was repeated until all trees had been exhausted. The number of well-spaced trees was recorded, and the entire evaluation repeated through 1000 iterative steps in an attempt to find the maximum well-spaced value. Fitting Simplified Models Initially, two sets of simplified models were created, using as predictive parameters the various sets of survey summary statistics (Xs j where j = survey method) and, where appropriate, site index (Table 4-1). The first of these was intended to explore the questions of survey method (quadrat vs. 5-5-2-2 vs. mixedwood well-spaced) and quadrat plot size (5 m2 to 20 m2 in 2.5 m2 increments). In these cases, site index was fixed at 20 m for both spruce (SI Sw ) and aspen (SI At ), roughly approximating the median site productivity for mixedwood stands in northeastern BC. The second set of models used a fixed survey method (10 m2 quadrats) and incorporated site index of both spruce and aspen to provide broad applicability across a wide range of site productivity. The default option for fitting all models was a linear combination of the untransformed survey summary statistics and, where appropriate, SI Sw and SI At . Very simple linear models using such untransformed statistics worked well in early exploratory studies (Farnden 2009) where site index was not a factor, but in this study such models were often unable to handle the increasing complexity. In a few cases, models were sufficiently linearized using transformations of the survey statistics, but more often a nonlinear model form was required. In these cases, a linear fit was first attempted, and the shape of the data cloud on the fit plots was used to indicate an appropriate nonlinear adaptation. For most models, the fit plots indicated that the relationships continued to be very nearly linear over a wide range of predicted and actual values. Where the lower tail of the point cloud tended to a horizontal asymptote, a sigmoidal model form such as a Weibull or Chapman-Richards (sensu Sit and Poulin-Costello 1994) was suggested. Where the lower tail of the point cloud tended toward a vertical or steeply sloped asymptote, a logarithmic or mixed-type function (sensu Sit and Poulin-Costello 1994) was indicated. In all of the nonlinear model forms, a complex variant of the model was required whereby the independent variable “X” was replaced by a linear combination of the survey summary statistics and, where appropriate, SI Sw and SI At . In the process of fitting the models that included site index, two concerns arose. The first was that aspen site index is difficult to measure (Jones et al. 1985)and is often unavailable, primarily Table 4-1. Basic description of simplified models used to predict total yield at age 80 (Y80) and the proportion of that yield as conifer volume (Yc). Xs is a linear vector of survey summary statistics, and j is a code indicating various survey methodologies. Model Scope Generalized Model Parameters Description Various Surveys; Site Index Fixed Y80 j = f(Xs j ) Yc j = f(Xs j ) Separate models developed based on each survey method; site index for both spruce and aspen is fixed at 20 m. Single Survey (10 m2 Quadrat); Site Index as Variable Y80 = f(Xs, SI Sw , SI At ) Yc = f(Xs, SI Sw , SI At ) Single pair of models developed using both spruce and aspen site index as independent predictor variables; SI At drawn randomly from range of 13 to 27; for SI Sw there were 118 cases for each of a) conversion equation (SI Sw = 3.805 + 4.786 x SI At ) without adjustment (Nigh 2002), b) conversion equation + Z r and c) conversion equation + Z r *3 where Z r is a random value between -1 and 1. Y80 i = f(Xs) Yc i = f(Xs) Separate models developed for each 1 m increment of spruce site index (i) from 13 to 27. Aspen site index determined as a function of spruce site index (Nigh 2002): SI At = -4.768 + 1.253 x SI Sw 69 due to indistinct growth rings. This factor can be dealt with by assuming that aspen and spruce site indices are highly correlated, and that aspen site index can be adequately predicted from spruce site index. A conversion equation based on this assumption is provided by Nigh (2002). The second concern was that curve shape may vary with site productivity. This factor was at least partially dealt with by making the curve shape parameters to be functions of site index. However, suspicion remained that variation in curve shape was not adequately addressed. To further evaluate the impacts of curve shape, a set of models was created for 15 discrete levels of spruce site index (13 to 27 m) using the survey outcomes for 10 m2 quadrats as the sole predictive parameters. Aspen site index was fixed as a function of spruce site index using the Nigh (2002) conversion function. Most of the simplified models used the 113 stands as discrete data points. For the models where spruce and aspen site index vary independently, each of the original 108 stands were duplicated twice and the extra stands four times, producing three groups of stands with one copy of each of the original stands and two copies each of the extras. Within each group, each stand was randomly assigned a site index value for aspen between 13.0 and 27.0 m, and an initial value of spruce site index was calculated based on the Nigh (2002) conversion equations. Final spruce site index values were determined by adding a random variable ranging from -1 to 1 for one of the groups, and from -3 to 3 for a second group. In the third group, spruce site index was unaltered. It was assumed here that the data used by Nigh (2002) for his conversion equations were a reasonable representation of variation that might be expected between spruce and aspen site index on the same site; the random variates then were intended to mimic the limits of Nigh’s (2002) SI AT by SI Sw data. All models were generated using JMP v. 8.0.1 software (SAS Institute 2008). Linear models were fit using ordinary least squares (OLS). In cases where there was no theory available to linearize relationships, multiple transformations were tested on one or more parameters, with parsimonious combinations explored through stepwise methods. Nonlinear models were developed using a Gauss-Newton search algorithm and least squares solutions where satisfactory linear models were not found. Where several model forms appeared feasible, (as with sigmoidal shapes), the option which produced the best root mean squared error (RMSE) was chosen for cases where the number of model parameters were the same, or AICc (Hurvich and Tsai 1989) where the number of parameters were different. 70 Results Calibration of Sub-Models in the Stand Generator Data for calibration of the sub-model controlling height adjustments related to overtopping competition are illustrated in Fig. 4-1, along with accompanying model parameters and fit statistics. Similar data, parameters and fit statistics for the diameter prediction sub-models are illustrated in Fig. 4-2. Figure 4-1. Two-parameter cumulative Weibull fit of relative spruce height (adjustment factor) as predicted by spacing factor. Spacing Factor was calculated using trees in the overtopping aspen canopy. Relative spruce height was determined by first fitting the curve using absolute spruce height, then dividing each height value by the asymptote. Figure 4-2. Linear fit of diameter as predicted by height for aspen (left) and spruce (right). For aspen, diameter was logarithmically transformed prior to fitting, leading to the curved appearance of the line on standard axes, and the use of the I2 statistic (1-SSE/SSY using y values in original units and back- transformed residuals) in place of the R2 statistic. y = 1− e−377.2x 3.228 RMSE = 0.19 xey 2797.04106.0 +−= I2 = 0.87 RMSE = 0.67 R2 = 0.94 RMSE = 0.42 71 Simplified Models Based on Various Survey Methods and Fixed Site Index – Original Data Set A set of models was constructed based on the original 108 stands to contrast various survey methods for their ability to produce summary statistics useful for predicting Y80 and Yc. Site index was fixed at 20 m for both spruce and aspen. Acceptable linear models were found in all cases except one, where a complex Weibull was required (Eq. 5): All Quadrat Surveys (2) (3) 5-5-2-2- Survey (4) (5) Mixedwood WS Survey (6) (7) In all cases, b 0 through b 6 are separate sets of coefficients for each model. Model parameters and fit statistics are provided in Tables 4-2 and 4-3, with a sample of fit plots in Fig. 4-3. For the nonlinear model, the R2 statistic has been replaced by the pseudo-R2 (or I2) statistic (1-SSE/ SSY) which has similar characteristics but is not necessarily bounded by 0 and 1. All linear models produce residuals that show no obvious pattern of nonlinearity, and residuals were homoschedastic and normally distributed. Based on assessments using RMSE, the best fits for both Y80 and Yc are produced using the mixedwood well-spaced survey summaries, followed by those using the smaller sized quadrats (≤ 10 m2). Fit quality generally declined for increasing quadrat sizes above 10 m2. The 5-5- 2-2 summary data provided a fit for Y80 that was almost as good as the better quadrat fits, but provided an inferior fit for Yc. For the models using the original 108 stands, it is evident from the bottom row in Fig. 4-3 that the scope of usage would be limited to stands wherein the predicted proportion of spruce is greater than approximately 30%. For practical applications this was viewed as a serious flaw, thus necessitating expansion of the data set for further model construction. Y80 = b0 + b1C0 + b2C1 + b3C2 Yc = b0 + b1C0 + b2C1 + b3C2 Y80 = b0 + b1 ln(1+ C0) + b2 C1 + b3C2 + b4 C2 + b5 C3 Yc = 1− e−b5 b0 +b1 C 0 + b2 C 2 +b3 C 3 +b4 ln C 0 +1( )( ) b6 Y80 = b0 + b1WS2 + b2WS3 Yc = b0 + b1WS2 + b2WS3 + b3WSr 72 Simplified Models Based on Various Survey Methods and Fixed Site Index – Expanded Data Set A second set of models was constructed to predict Y80 and Yc using survey summary statistics from the expanded data set. As before, site index was fixed at 20 m for both spruce and aspen. Outcomes were based either on linear fits using ordinary least squares (OLS) regression or a complex variant of either a Weibull, Chapman-Richards, mixed type or logarithmic function where the dependent variable X is replaced by a linear function of the survey summary variables: Table 4-2. Coefficients and fit statistics for the models used to calculate Y80 by survey method with site index fixed at 20 m. Survey Y80 b 0 b 1 b 2 b 3 b 4 b 5 R2 RMSE (m3/ha) Q5 344.5 -276.3 -97.8 37.2 0.895 13.59 Q7.5 313.7 -262.5 -79.0 -7.1 0.901 13.25 Q10 297.6 -256.3 -68.3 -27.9 0.903 13.08 Q12.5 298.0 -278.9 -96.0 -52.9 0.879 14.65 Q15 293.4 -274.9 -105.0 -55.28 0.854 16.08 Q17.5 288.0 -310.4 -106.9 -51.5 0.856 15.96 Q20 285.5 -326.57 -124.45 -58.4 0.816 18.02 5522 512.5 -669.4 -274.0 -505.5 318.4 -73.9 0.893 13.87 WS 57.6 11.7 24.5 0.927 11.28 Table 4-3. Coefficients and fit statistics for the models used to calculate yield-based species composition (Yc) by survey method with site index fixed at 20 m. Survey Yc b 0 b 1 b 2 b 3 b 4 b 5 b 6 R2 or I2 RMSE (%) Q5 1.133 -0.331 -1.040 0.080 0.883 0.061 Q7.5 0.907 -0.112 -0.828 0.187 0.900 0.057 Q10 0.824 -0.035 -0.774 0.222 0.893 0.059 Q12.5 0.777 -0.023 -0.757 0.250 0.868 0.065 Q15 0.758 -0.076 -0.727 0.279 0.850 0.069 Q17.5 0.749 -0.047 -0.795 0.278 0.840 0.072 Q20 0.750 -0.057 -0.945 0.258 0.846 0.070 5522 0.616 -1.664 0.620 0.670 2.167 2.672 2.096 0.762 0.088 WS 0.774 0.121 -0.076 -0.020 0.913 0.053 73 10 m2 Quadrats 5-5-2-2 Mixedwood WS Figure 4-3. Fit plots for simplified models using summary statistics from10 m2 quadrat surveys, 5-5-2-2 surveys, and mixedwood well-spaced surveys as parameters to predict total volume and species composition (% spruce volume). Quadrat Surveys for Q5, Q7.5 (8) all others (9) for Q 17.5, Q 20 (10) all others (11) 5-5-2-2- Survey (12) (13) Mixedwood WS Survey (14) (15) Y80 = b6 1− e −b4 b0 +b1C0 + b2C1 + b3C 2( ) b5 Y80 = b0 + b1C0 + b2C1 + b3C2 Yc = b0 + b1C0 + b2C1 + b3C2 Yc = b4 b0 + b1C0 + b2C1 + b3C2 + b5 b0 + b1C0 + b2C1 + b3C2( ) + b6 Y80 = b0 + b1 ln(C 0 +1) + b2 C 0 + b3 C1 + b4C 2 + b5 C 2 Yc = 1− e −b 6 b0 +b1C1 + b2 C1 +b3C 2 +b4 C2 +b5 C 3( ) b7 Y80 = b0 + b1WS1 + b2WS3 + b3WS4 Yc = b0 + b1 WS1 + b2WS2 + b3 WS2 + b4 ln WS3( ) 74 In all cases, b0 through b6 are separate sets of coefficients for each model. Model parameters and fit statistics are provided in Tables 4-4 and 4-5, with a sample of fit plots in Fig. 4-4. For nonlinear models, the R2 statistic has again been replaced by the I2 statistic. All models produce residuals that show no gross pattern of nonlinearity. While most models were fit using a linear function with the original data set, most models for the expanded data set required nonlinear forms. For Y80, it became evident that the additional data points accentuated a small but previously trivial nonlinear pattern in the original data set. In comparing the nonlinear forms to the linear forms using AICc, the gains over using a linear form were small for surveys Q5, Q7.5, 5522 and WS, and insufficient to justify a nonlinear model form for the remainder. Table 4-4. Coefficients and fit statistics for the models used to calculate Y80 by survey method with site index fixed at 20 m. Survey Y80 b 0 b 1 b 2 b 3 b 4 b 5 b 6 R2 or I2 RMSE (m3/ha) Q5 19.38 -0.554 -0.400 -0.136 9.298e-152 118.15 387.5 0.921 15.87 Q7.5 131.54 -63.00 -36.69 -10.67 6.247e-8 3.605 433.1 0.943 13.4 Q10 402.9 -312.4 -161.4 -55.9 - - - 0.925 15.2 Q12.5 392.6 -339.7 -171.8 -70.0 - - - 0.924 15.4 Q15 382.6 -333.6 -169.3 -67.5 - - - 0.884 19.0 Q17.5 374.7 -364.2 -170.7 -64.4 - - - 0.871 20.0 Q20 367.0 -404.3 -171.9 -62.6 - - - 0.835 22.6 5522 520.5 -202.4 -283.3 -584.8 440.9 -202.9 - 0.879 19.6 WS 11.87 0.8311 0.3810 0.2141 7.568e-5 3.279 498.2 0.951 12.5 Table 4-5. Coefficients and fit statistics for the models used to calculate yield-based species composition (Yc) by survey method with site index fixed at 20 m. Survey Yc b 0 b 1 b 2 b 3 b 4 b 5 b 6 b 7 R2 or I2 RMSE (%) Q5 0.2136 -0.008452 -0.1795 0.09308 -0.02506 2.766 0.3173 - 0.863 0.07109 Q7.5 1.6563 0.2519 -1.2694 1.1698 -0.2615 0.2657 0.3608 - 0.889 0.0643 Q10 0.3089 0.0598 -0.1842 0.1876 -0.0984 1.3530 0.5014 - 0.888 0.0645 Q12.5 1.0082 0.2848 -0.8165 0.7033 -0.0972 0.4666 0.2192 - 0.875 0.0679 Q15 0.3769 0.0834 -0.3107 0.2790 -0.0214 1.2797 0.1623 - 0.862 0.0716 Q17.5 0.5943 0.1790 -0.5571 0.3756 - - - - 0.866 0.0697 Q20 0.5939 0.1195 -0.6070 0.3685 - - - - 0.868 0.0692 5522 2.983 5.497 -8.318 -0.2696 1.389 1.424 1.175 0.8408 0.835 0.0789 WS 0.4554 -0.1985 -0.1066 0.6269 -0.1321 - - 0.925 0.0526 75 For Yc, it became evident that the linear models used for the original data set would lead to heavily biased results when applied to stands with low levels of spruce composition (Fig. 4- 4). For these models, a logarithmic or mixed type of model was employed to retain the largely linear nature of the data in the original set along with the downward curve for stands with species compositions below approximately 30%. These model forms provided better fits at the bottom end of %Sw composition for some surveys than for others, as evidenced by the fit plots in the bottom row of Fig. 4-5. Figure 4-4. Fit plot of linear model for species composition using untransformed variables from the mixedwood well- spaced survey. The 5 points below the diagonal line at the bottom of the chart are the additional points for stands with low spruce presence. All of the data sets showed this trend whereby fits using the original data set would over-predict spruce volume where it is a minor component of the entire stand. 10 m2 Quadrats 5-5-2-2 Mixedwood WS Figure 4-5. Fit plots for simplified models using the expanded data set and summary statistics from10 m2 quadrat surveys, 5-5-2-2 surveys, and mixedwood well-spaced surveys as parameters to predict total volume and species composition (spruce volume proportion). 76 For the most part, trends in model fit by survey type using the expanded data set follow those for the original data set. Mixewood well-spaced surveys produced the best fits, and the 5-5-2-2 survey the worst. Within the range of quadrat plot sizes, surveys using smaller plots performed better than those using larger plots. Trends were somewhat more difficult to discern with clarity due to the mixture of model forms. Quadrat-Based Models Incorporating Spruce and Aspen Site Index Models incorporating site index of spruce and aspen as continuous variables were based on summary statistics from the 10 m2 quadrat surveys. Models for Y80 and Yc were created based on nonlinear least squares fits of equations based on a Chapman-Richards function. For both Y80 and Yc, the X variable in the Chapman-Richards function was replaced by a linear combination of variables (survey summary statistics and site index). In the equation for Yc (Eq. 17), the Chapman-Richards shape parameters become linear functions of the average of site index for spruce and aspen. (16) (17) where: b 0 through b 10 are separate sets of coefficients for each model and SI Sw = White spruce site index, SI At = Aspen site index, and = Mean of spruce and aspen site indices. Model parameters and fit statistics are provided in Table 4-6. Fit plots are provided in Fig. 4-6. Models were also fit using a Weibull function in place of the Chapman-Richards function, but were found to produce less precise fits based on RMSE and I2 statistics. Y80 = b8 * 1− e −b6 (b0 +b1C0 +b2C1 +b3C 2 +b4 SI Sw +b5SI At )( )b7 Yc = 1− e − b6 + b7 *SI x�+ b8 *SI x� 2( ) b0 +b1C 0 +b2C1 +b3C 2 +b4 SI Sw +b5SI At( ) b9 + b10 *SI x�( ) Table 4-6. Coefficients and fit statistics for the models used to calculate total yield at age 80 (Y80) and yield-based species composition (Yc) using summary data from a 10 m2 quadrat survey and both aspen and spruce site index. Model b 0 b 1 b 2 b 3 b 4 b 5 I2 RMSE (m3/ha) Y80 1492 -140.9 -59.71 -31.65 8.317 11.03 0.985 23.88b6 b7 b8 - - - 0.004634 8734.8 1739.7 - - - Model b 0 b 1 b 2 b 3 b 4 b 5 I2 RMSE (%) Yc 1.1586 0.08177 -0.9303 0.4847 0.06522 -0.07809 0.928 0.0710b6 b7 b8 b9 b10 - 10.46 -0.8375 0.01836 0.7298 -0.004347 - 77 Quadrat-Based Models For Discrete Levels of Spruce Site Index In the cases where separate simplified models for Y80 and Yc were developed for discrete 1 m intervals of SI Sw , fits were based either on ordinary least squares (OLS), a complex variant of a Weibull function or a complex variant of a logarithmic function: for SI = 13 to 21 (18) for SI = 22 to 27 (19) for SI = 13 to 18 (20) for SI = 19 to 22 (21) for SI = 23 to 27 (22) In all cases, b 0 through b 6 are separate sets of coefficients for each model. Model parameters and fit statistics are provided in Tables 4-7 and 4-8. All linear models produced residuals that showed no obvious pattern of nonlinearity, were homoschedastic and normally distributed. In general, fits were better based on RMSE values for Y80 at low site indices. However, RMSE values for Y80 cannot be directly compared as each model predicted a different range of volumes (better sites have higher yields). To correct for this problem, RMSE has been standardized by dividing it by the mean predicted volume for all stands. For Yc, model fits were poorer for the mid-range site indices and better toward the extremes. A possible explanation is that stands are tending towards single species prevalence at the extremes, with spruce becoming dominant at low site indices and aspen becoming dominant on the most productive sites. Figure 4-6. Fit plots for simplified models for (a) Y80 and (b) Yc that include site index of both spruce and aspen as predictive parameters. a b Y80 = b0 + b1C0 + b2C1 + b3C2 Y80 = b6 * 1− e−b4 (b0 +b1C0 +b2C1+b3C2 ) b5 Yc = 1− e−b4 (b0 + b1C0 + b2C1 + b3C 2 ) b5 Yc = b4 + b5 log b0 + b1C0 + b2C1 + b3C2( ) Yc = b0 + b1C0 + b2C1 + b3C2 78 Table 4-7. Coefficients and fit statistics for the models used to calculate Y80 assuming discrete site index values. Separate models have been generated for different levels of site productivity based on 1 m increments of white spruce site index. Sw SI 50 Y80 b 0 b 1 b 2 b 3 b 4 b 5 b 6 RMSE (m3/ha) 13 89.38 -74.63 -68.79 9.52 3.80 0.061 14 109.47 -90.22 -72.35 12.52 4.45 0.056 15 135.28 -110.34 -77.09 12.77 5.28 0.052 16 168.83 -136.33 -84.62 8.94* 6.49 0.051 17 212.86 -170.26 -94.54 -0.48* 8.31 0.051 18 270.41 -214.11 -113.56 -17.50 10.51 0.051 19 338.82 -265.93 -136.47 -40.43 13.04 0.051 20 415.78 -323.60 -161.10 -68.23 15.85 0.051 21 496.45 -382.77 -178.63 -97.47 18.84 0.051 22 36.53 -15.413 -9.318 -7.419 1.772e-7 4.734 533.84 21.77 0.051 23 26.61 -10.244 -5.861 -5.0277 1.076e-7 5.372 589.8 23.58 0.048 24 8.544 -0.1417 -0.07340 -0.07096 1.336e-151 162.8 647.5 28.43 0.052 25 5.750 -1.0310 -0.4970 -0.5597 1.357e-10 14.08 702.6 27.24 0.044 26 2.829 -0.3355 -0.1802 -0.1834 2.516e-11 25.74 742.4 29.5 0.044 27 1.518 -0.001953 -0.001184 -0.001104 1.1953e-106 1741 787.56 30.47 0.043 Table 4-8. Coefficients and fit statistics for the models used to calculate yield- based species composition (Yc) assuming discrete site index values. Separate models have been generated for different levels of site productivity based on 1 m increments of white spruce site index. Sw SI 50 Yc b 0 b 1 b 2 b 3 b 4 b 5 RMSE (%) 13 1.1240 -0.6415 -1.1225 -0.2422 8.868 1.814 0.0349 14 0.8877 -0.5186 -0.9103 0.0938 6.024 1.360 0.0421 15 0.9351 -0.4746 -0.9621 0.4181 3.567 1.194 0.0498 16 1.1150 -0.4773 -1.1511 0.8132 2.059 1.087 0.0568 17 0.6057 -0.1881 -0.6284 0.6201 2.829 1.027 0.0628 18 0.5691 -0.0594 -0.5955 0.7979 2.101 1.040 0.0650 19 0.6052 0.0906 -0.4225 0.4346 0.9468 0.6249 0.0647 20 0.5753 0.1180 -0.3072 0.3447 1.020 0.8275 0.0650 21 0.5567 0.0724 -0.1515 0.1588 1.514 1.702 0.0645 22 0.7754 0.0124 -0.0264 0.0241 4.207 14.57 0.0618 23 0.5011 0.2052 -0.5115 0.4461 0.0580 24 0.5147 0.1647 -0.5400 0.4301 0.0550 25 0.5399 0.1074 -0.5778 0.4021 0.0520 26 0.5757 0.0587 -0.6158 0.3618 0.0500 27 0.6273* 0.0268 -0.6551 0.3016 0.0496 79 Comparison of Errors for Models Using Discrete Versus Continuous Site index A comparison of the single pair of predictive models that incorporates spruce and aspen site index as continuous variables (Eqs. 16 and 17) with the 15 pairs of models for discrete levels of site index (Eqs. 18 through 22) can be accomplished by contrasting their prediction errors (RMSE values). In doing so, the errors for the continuous SI models must be parsed by site index to match those for the discrete models. To accomplish this, the spruce site index for each data point in the continuous model was rounded to the nearest whole number and used as a bin value to group the data points. The root mean squared error was then recalculated for the values in each bin. The outcomes are presented in Table 4-9 and Fig. 4-7. In comparing the models for predicting Y80 (Fig. 4-7), it initially appears as if the discrete models have smaller errors, with the mean difference between the continuous and discrete models being 6.9 m3/ha. However, it is important to recognize the artificial error term introduced into the discrete models by assuming site index values to be whole numbers. In comparison to the Y80 model using site index as a continuous variable, the discrete models require a mean input error for site index of ±0.25 m. For most levels of site index, this leads to an implied additional prediction error (Fig 4-8) that in most cases exceeds the differences in RMSE. This factor shifts the advantage in predictive precision to the model using SI as a continuous variable rather than the discrete models. Table 4-9. RMSE values from the models incorporating site index, parsed by 1 m intervals of site index SI Y80 (m3/ha) Yc (proportion Sw) 13 7.8 0.085 14 14.7 0.076 15 12.2 0.077 16 16.8 0.071 17 16.5 0.078 18 19.9 0.077 19 15.6 0.055 20 21.2 0.063 21 21.5 0.059 22 27.9 0.079 23 32.5 0.077 24 27.9 0.062 25 33.0 0.047 26 26.9 0.062 27 56.3 0.054 80 Similar comparisons can be made for the models predicting Yc. As with Y80, the RMSE values are higher for the continuous SI models than for the discrete SI models, although in this case the concern is the interspecies difference in site index rather than the absolute level. For predictions of Yc on spruce site indices above 18, the mean difference in RMSE for the discrete versus continuous SI models is 0.0048 (Fig 4-7). However, dropping below SI 18 this difference increases dramatically to as much as 0.05. To assess the impact and possible causes of these changes, a controlled set of MGM runs was conducted to illustrate the changes in Yc that must be accommodated by the continuous SI model, with the outcomes presented in Table 4-10. At site index 24, there is a range in Yc values from 0.44 to 0.68 that must be explained by the continuous model but not the discrete models, so an increase in RMSE of 0.0048 seems entirely reasonable. At site index 16, there is an increased range in Yc ranging from 0.45 to 1.00 that must be explained, so again the increase in RMSE of Figure 4-7. Comparison of errors (RMSE) between models fit for discrete levels of site index (from Tables 4-7 and 4-8) and those parsed from the fit of a model that includes site index as a predictive variable (Table 4-9). 0 10 20 30 40 50 60 10 15 20 25 30 Site Index (m) R M S E ( m 3/ h a) Parsed RMSE Discrete RMSE Y80 0 0.02 0.04 0.06 0.08 0.1 10 15 20 25 30 Site Index (m) R M S E (p ro po rt io n S w ) Parsed RMSE Discrete RMSE Yc Figure 4-8. Mean change in total volume for each 0.25 m increment in spruce site index, based on MGM predictions used in development of the discrete site index models. 0 5 10 15 20 10 15 20 25 30 Site Index (m) C h an g e in T o ta l V o lu m e (m 3 /h a) p er 0 .2 5 m S I I n cr em en t 81 approximately 0.01 seems entirely reasonable. Overall, the increased RMSE for the continuous models appears to be correlated with an increased range of values that must be explained, with the increased error being a very small proportion of the increased range. It appears that the continuous models do a good job of explaining the increased variation in Yc that results from changing the balance of SI Sw verus SI At . It is likely that the divergence in RMSE pattern illustrated in Fig. 4-6 is indicative of a shift in the competitive balance between aspen and spruce on poorer sites. As overall site index decreases, there is a tendency in the simulation outcomes toward higher proportions of spruce volume, particularly where the ratio of spruce to aspen site index increases. In the discrete models, this trend necessitated a change in model forms (Eqs. 20 to 22), which in turn led to concerns of bias at low site indices for the continuous model which is limited to a single model form. Highlighting specific site index ranges in the fit plot for this model (Fig 4-9) provides a hint that such a bias may exist for stands with low spruce proportions, but the evidence is not strong. There is also a hint of a nonlinear pattern in the residuals for the high site index stands. A final comparison of the continuous versus discrete models is warranted wherein the continuous models are modified and refit to work within the conditions for which the discrete models were designed; that is, where SI At is unknown. In this case, models 16 and 17 become: (22) (23) Table 4-10. Species composition (proportion Sw) at age 80 as simulated in MGM for varying levels of spruce and aspen site index. For each level of spruce site index, aspen site index was determined using the Nigh (2002) conversion equation then adjusted upward or downward in 1 m increments. The common stand used for all simulations had 7% coverage with aspen clones, 800 trees/ha in the inter-clonal areas, and approximately 1000 trees/ha of planted spruce. SI At Proportion Sw Volume by SI Sw 16 20 24 +3 0.45 0.36 0.44 +2 0.56 0.41 0.46 +1 0.71 0.47 0.48 0 0.88 0.56 0.51 -1 0.97 0.68 0.55 -2 0.99 0.83 0.61 -3 1.00 0.94 0.68 Y80 = b7 * 1− e−b5( b0 +b1C0 +b2C1+b3C2 +b4 SISw )( )b6 Yc = 1− e − b5 +b6 *SISw +b7 *SISw 2( ) b0 +b1C0 +b2C1+b3C2 +b4 SISw( ) b8 +b9 *SISw( ) 82 Only very small differences in the predictive value were found for Y80, with equation 16 performing slightly better than Eq. 23; such should not be surprising given the generally accurate fits found for almost all of the models predicting Y80. A comparison of fits for Eqs. 17 and 23 for Yc produces a different picture (Fig 4-10), where there is considerable loss in prediction precision with the loss of SI At as an independent variable. In comparing this single SI fit to the discrete models, the RMSE value of 0.1077 is considerably poorer than all of the RMSE values for the discrete SI models as presented in Table 4-8. SI Sw < 15 SI Sw 19 to 21 SI Sw > 24 Figure 4-9. Selection of points (bold red) by site index range in the fit plot for the Yc model incorporating site index as continuous variables (Eq. 16). There are no obvious signs of bias, except possibly for the lowest spruce proportions over the lowest range of site indices. a b Figure 4-10. Comparison of fit plots for species composition models using (a) both SI Sw and SI At as continuous variables versus where I2 = 0.928 and RMSE = 0.0710 versus b) only SI Sw . as a continuous variable with SI At inferred, where I2 0.845 = and RMSE = 0.1077. 83 Discussion The approach used to predict species composition in this study arises from the need to manage landscape level species composition as described in the Sustainable Forest Management Plan for the Fort St. John Pilot Project (Fort St. John Pilot Project 2004). In that initiative, landscape level species composition is defined based on the percentage of land area in each of four classes (≤20% conifer, 20%80%). Classes are loosely defined based on volume or basal area for stands ≥ 20 years of age. Given that species composition may change dramatically over time given mid-rotation aspen harvests and the strong pattern of succession that is expected, it was decided early in this study to refine the definition to using a reference age of 80 years (Y80), and to predict species composition assuming no post establishment management interventions. The implication of this strategy is that the spruce component (Yc) will be overestimated for stands younger than the reference age, and underestimated for stands older than the reference age. However, it is important to remember that it is not the condition of any one stand that is important but the overall condition of the landscape. Comparison of Survey Systems For both Y80 and Yc, the best model fits resulted using input data from the mixedwood well- spaced survey, followed closely by surveys using small quadrat plots. The poorest model fits resulted from using the larger quadrat plots and the 5-5-2-2 survey. There was less variation between survey methods for predicting Y80 than for Yc, and greater overall precision. Simulation (e.g. Bergerud 2002) and empirical studies (e.g. Braathe 1982) have demonstrated a high correlation between stocking measures and yield, which is echoed in the high correlations found here, despite large variations in species composition in the initial conditions for various stands. It is surmised that overall stocking measures are far more critical for estimating yield than are species specific measures, and that all of the survey methods employed here used satisfactory estimates of overall stocking, or conversely, methods to detect unoccupied growing space. The better fits for Y80 models using small (≤ 10 m2) quadrats as opposed to larger quadrats suggests that small plots are more effective at quantifying variation in stocking than are large plots, at least for the mixedwood stands in this study. However, these results are different than those presented for single-species conifer stands presented in Chapter 3, where larger plots in the neighborhood of 16 m2 provided better predictive parameters. At least two possible explanations exist. First, it is possible that the stand growth modeling strategy employed in this study masked the stocking variation that existed in the simulated stands. Such a case could result 84 if the localized treelists used for MGM growth simulations were capturing too many trees around small gaps, resulting in half void / half tree patch treelists. This in turn might have skewed the relationship between the then erroneous measures of stocking and the simulated yields. The possibility of such an occurrence was anticipated in designing the study, and a plot size (100 m2) was selected that was a compromise between too large a plot that would excessively negate small voids and too small a plot that would not capture a sufficient number of trees to represent a realistic treelist. At the present time there is no easy way to test the validity of such an explanation for these stand types, although a new version of the spatially explicit TASS9 model is expected shortly that may offer opportunities for further testing. A second explanation for the recommended quadrat size differences between this study and that presented in Chapter 3 relates to the differences in stand composition. It seems feasible that for the mixed species stand, a larger quadrat size increases the liklihood that quadrats that include a spruce tree will also include an aspen as well. This shift in the balance between spruce (C 2 ) and mixed (C 3 ) plots may contribute to a loss in predictive precision. However, if such were the case one might also expect similar or greater loss in predictive precision with increasing quadrat size for species composition, which was not the case. In predicting Yc, the relatively poor performance of the 5-5-2-2 survey was initially surprising, given that it was specifically designed for the type of stand condition being evaluated in this study. In retrospect, the complexity of the survey process and the fact that it was intended to address the mixedwood issue from a somewhat different perspective are its primary shortcomings when it comes to predicting species composition. The survey evolved from some preliminary work done by Dr. Phil Comeau at the University of Alberta, where aspen gap size was related to minimum light thresholds required for unimpeded understorey spruce height growth. In the BC survey methodology applying this concept, the objective is to identify and tally the largely unimpeded spruce trees, ignoring the rest (Fig. 4-11). The result is that the survey does a reasonably good job of discerning species composition when a large percentage of spruce are tallied, but a poor job when a large percentage fail the 5-5-2-2 test. As a result, the method does a much better job of predicting species composition where the unimpeded spruce component is high, and a much poorer job where it is low (Fig. 4-12). Conversely, the mixedwood well-spaced survey appears to be particularly effective for use in predicting species composition. Again, many trees may not be tallied, but in this case only if they 9 The Tree And Stand Simulator (TASS) is an individual tre, spatially explicit growth model originally developed by Mitchell (1969, 1975) and further refined and expanded over many years by the BC Ministry of Forests’ Research Branch. Version 3 of this model, scheduled for release in 2010, is expected to be able to simulate stratified mixtures of spruce and aspen. 85 are redundant (they are deemed to be occupying the same growing space as another tree). In this survey, the competition gradient is explicitly recognized: each understorey spruce is tagged with the number of well-spaced aspen occurring in the same plot. The mean tag value (WS r in Eq. 6) becomes a measure of the overall competition experienced by spruce throughout the stand, and appears to be effective at helping to quantify the species-specific yield impacts of the stratified competition. This becomes a much better method of quantifying the apportionment of growing space between the spruce and aspen components at the stand level than was provided by the 5-5- 2-2 survey. The quadrat survey also appears to be effective at quantifying growing space apportionment between species, but does so in a simpler manner. In this case, it is unimportant whether any particular tree is experiencing a high or low level of competition – the system works instead by using probability and a large number of observations. Growing space apportionment is quantified using the relative number of plots in each of the stocked categories (broadleaved vs. conifer vs. both). The probability of any one plot being classed as conifer is related to stand level broadleafed (aspen) numbers; given a relatively consistent distribution of spruce the ratio of conifer to mixed plots will increase with decreasing levels of overstorey aspen competition. Figure 4-11. Survey assessments that employ competition thresholds for tallying trees will always ignore a component of the stand. Such effects can be critical in assessments of species composition, particularly when a large percentage of one species or species group are not tallied. ����������� ���������� ���������� ����������������������� ������������������� �������������������� ������������� ���������� Figure 4-12. Linear fit of Yc using summary data from the 5-5-2-2 survey summaries. In stands expected to have low spruce volume components, very few spruce trees are tallied as acceptable, and models fit to such data tend to interpret very low levels of detected spruce as predictors of considerable spruce volume. 86 In the final analysis, the choice of survey systems to be used under operational conditions is dependent not only on the predictive ability of the survey summary statistics, but also on the ease and practicality of application. The mixedwood well-spaced method produces summary statistics that are excellent predictors of yield and species composition, but the survey method is relatively complex and has three separate well-spaced assessments each with a high degree of subjectivity (which are the well-spaced trees?). While not assessed in this study, actual field application of this method could realistically result in sufficient measurement errors such that all predictive advantages provided by superior model fits would be lost. In contrast, the field application of quadrat plots is simple, with few potential sources of measurement error. This makes the choice of a quadrat-based survey system very attractive. A quadrat survey using 10 m2 plots would not only provide very good predictions of Y80 and Yc, but would be reasonably consistent with many survey system already in use in many jurisdictions. Options for Inclusion of Site Index Two methods have been presented for including site index as a predictive parameter: using both spruce and aspen site index (SI Sw and SI At respectively) as predictive parameters in a single pair of models, and constructing a separate pair of models for each discrete level of spruce site index. In the latter case, there is an assumption of a linear relationship between spruce and aspen site indices following Nigh (2002). In the single pair of continuous site index models, the major advantage is that SI Sw and SI At can be varied independently. Green (2004) and Green and Hawkins (2005) cautioned that variations in environmental gradients can cause changes in the relative productivity of multiple species, which can have important implications for successional dynamics. While the degree of inter-specific site index variation that might exist was not addressed by this study, the potential implications were. It is clear that a 1 to 2 m variation in the site index of one species while the other is held constant can have large effects on predictions of Yc (changes of 30% or more). However, a major problem exists in using models that incorporate aspen site index. Aspen trees are notoriously difficult to accurately age (Jones et al. 1985), and an error of a single year for a young tree can introduce a very large error in the site index estimate. For this reason, aspen site index is often not available or is unreliable. If this is the case, the major advantage of the models employing both site indices may be lost. Where only spruce site index is available there are two choices. First, the aspen site index can be inferred using Nigh’s (2002) conversion equation for use in the single pair of continuous models. Second, the discrete models can be employed. Each offers separate advantages. 87 Discrete models likely offer an advantage in that fits over the range of site indices are not restricted to a single model form, allowing for greater flexibility in curve shape. The curve shape problem was at least partially addressed in the continuous model for Yc by replacing the Chapman-Richards shape parameter with functions based on site index, but the same basic sigmoidal function persisted. In the discrete models for Yc, the functional forms included three general shapes for different site index ranges: Weibull, logarithmic and linear. While complex forms of both the Weibull and logarithmic models might have been adequately fit to encompass the range of site indices that used the linear forms, they are not readily interchangeable themselves. It appears, then, that using a single model form for the entire range of site indices introduces a compromise in fit over that provided by the discrete models. The continuous models offer the advantage of avoiding errors introduced by rounding spruce site index to even 1 m increments, but only if site index as a predictive parameter is realistically measured to a precision better than to the nearest metre. General Comments on Simplified models Beyond contrasting the various options for simplified models, the simplified models should be contrasted with the alternative: the original model from which they were derived. The uncertainty with which predictions can be made using a model such as MGM may or may not be well documented, but managers applying such models generally take it on faith that the model predictions are reasonable representations of the real world. In building and applying simplified models, the assumption is that a complex model such as MGM is impractical to apply on an extensive basis. In doing so, precision and usually a large degree of flexibility are sacrificed (Farnden 2009). For the simplified models, estimates of the additional uncertainty added to the predictive system through the unexplained errors are summarized in the RMSE values. Depending on the chosen strategy, the simplified models add a mean error ranging from 4 to 30 m3/ha for total volume, and from 0.035 to 0.071 for the proportion of spruce volume. From an application stand point, the loss in flexibility may be the more serious limitation. While attempts were made within this study to cover as broad a range of stand conditions as possible that might be experienced in the real world, there will always be conditions found that do not match the assumptions used in the modeling framework. Of immediate concern is the presence of significant components of other species (such as lodgepole pine) or large residual trees. It is important that users are aware that the effects of such variations in stand conditions are not evaluated in the simplified models, and that missing these effects may lead to erroneous conclusions. In my experience, many managers are all too willing to push application of an 88 available model far beyond the conditions and assumptions under which it was created, simply because it is the only tool available. The general process of modeling used in this study to develop regeneration standards also likely has another limitation. The simplified models developed here cover explicitly only two species, various combinations of which are sufficiently common in the boreal forest that the models should have wide applicability. However, extrapolation of the modeling framework used here to more complex stand conditions seems unlikely. The addition of a third species, for example, would greatly increase the complexity of the model forms and the cost of building the models, and would likely also increase considerably the magnitude of the error terms. Coupled with greater cost in complexity would be a considerable decrease in the frequency of real world cases to which the model could be applied. Beyond the relatively simple case of two-species stands, it is strongly suggested that the best strategy is to stick with the original model from which the simplified models would be developed (such as MGM). Conclusions In forest stands where a two species mixture occurs with silvically disparate species such as spruce and aspen, the interactions between the species and the resulting long term outcomes can be very important for the achievement of management objectives. Such interactions have historically been difficult or operationally expensive to predict. This study demonstrates that it is possible to develop simple models using silviculture survey summary statistics as input parameters to predict future yields by species for this type of two-species mixture. The quality of fits for such models will vary with the type of survey being used. Desirable survey characteristics include an ability to address the frequency of unstocked patches and the competitive status of trees in lower layers of a stratified mixture (quantified allocation of growing space). The primary undesirable characteristic is an on/off acceptability criteria for tallying trees based on competitive condition. For the three surveys evaluated in this study, the general performance ranking was: Mixedwood Well-spaced > Small Plot Quadrats >> (5-5-2-2 and Large Plot Quadrats) While the mixedwood well-spaced survey provided summary statistics that were slightly better as predictive parameters, the simplicity of small plot quadrat surveys likely makes them a better choice for operational application. Two options were developed for predicting total yield and the proportion of spruce in that yield based on different modeling strategies. The option where site index of both species is included as a continuous variable is best suited for situations where satisfactory estimates of site index are available for both species. However, it is recognized that aspen ages are difficult to 89 measure making it problematic to acquire reliable site indices for that species. The option where separate models have been developed for discrete 1 metre increments of spruce site index is best suited for conditions where aspen site index is unavailable and the spruce site index estimate is of relatively low precision and/or accuracy. For conditions where only spruce site index is available and the estimate of that value is believed to have good accuracy and precision, neither modeling approach appears to be clearly superior. 90 References BC Ministry of Forests. 2009. Silviculture survey procedures manual: stocking and free growing survey. BC Ministry of Forests, Forest Practices Branch, Victoria BC. 219p. Bergerud, W.A. 2002. The effect of the silviculture survey parameters on the free-growing decision probabilities and projected volume at rotation. Land Management Handbook No. 50. BC Ministry of Forests, Victoria BC. 32p. Braathe, P. 1982. Stocking control and its effect on yield. Forest Industry Lecture Series No. 10. University of Alberta, Edmonton AB. 21p. Carter, D.H.H. 2007. An assessment of variable radius plot sampling techniques for measuring change over time: a simulation study. Master of Science Thesis, University of British Columbia. 111p. Farnden, C. 2009. An analysis framework for linking regeneration standards to desired future forest conditions. For. Chron. 85: 285-292. Fort St. John Pilot Project. 2004. Sustainable forest management plan. Fort St. John Pilot Project, Fort St. John BC. URL: http://fsjpilotproject.com/sfmp.html, accessed Mar. 29, 2010. 274p. Green, D.S. 2004. Describing condition-specific determinants of competition in boreal and sub- boreal mixedwood stands. For. Chron. 80: 736-742. Green, D.S., and Hawkins, C.D.B. 2005. Competitive interactions in sub-boreal birch-spruce forests differ on opposing slope aspects. For. Ecol. Manage. 214: 1-10. Hurvich, C.M., and Tsai, C. 1989. Regression and time series model selection in small samples. Biometrika 76: 297-307. J.S. Thrower and Assoc. 2003. Stand survey and growth modeling for the Fort St. John TSA. Contract report to Canadian Forest Products Ltd., Chetwynd BC. 32p. J.S. Thrower and Assoc. 2002. Stand survey and growth modelling for the TFL 49 results-based pilot project. Contract report to Riverside Forest Products Ltd., Kelowna BC BC. 45p. Jones, J.R., DeByle, N.V., and Winokur, R.P. 1985. Wood resource. In Aspen: Ecology and management in the western United States. Edited by N.V. DeByle and R.P. Winokur. General Technical Report RM-119. USDA Forest Service, Rocky Mountain Forest and Range Experiment Station, pp. 161-167. 91 Kabzems, R., Nemec, A.L., and Farnden, C. 2007. Growing trembling aspen and white spruce intimate mixtures: early results (13-17 years) and future projections. BC Journal of Ecosystems and Management 8: 1-14. Kaltenberg, M.C. 1980. Evaluation of regeneration sampling methods: a Monte Carlo analysis using simulated stands. DNR Report No. 39. State of Washington, Dept. of Nat. Res. 50p. Mackisack, M.S., and Wood, G.B. 1990. Simulating the forest and the point-sampling process as an aid in designing forest inventories. For. Ecol. Manage. 38: 79-103. MacLeod, D. 1977. Accuracy of regeneration stocking estimates: tests with simulated data. For. Chron. 53: 77-81. Magnussen, S. 1999. Effect of plot size on estimates of top height in Douglas-fir. W. J. Appl. For. 14: 17-27. Martin, P.J., Browne-Clayton, S., and McWilliams, E. 2002. A results-based system for regulating reforestation obligations. For. Chron. 78: 492-498. Mitchell, K.J. 1975. Dynamics and simulated yield of Douglas-Fir. Forest Science Monographs 17. 39 p. Mitchell, K.J. 1969. Simulation of the growth of even-aged stands of White Spruce. Bull. Yale Sch. For. No. 75. 48 p. Nigh, G. 2002. Site index conversion equations for mixed trembling aspen and white spruce stands in Northern British Columbia. Silva Fenn. 36: 789-797. SAS Institute. 2008. JMP user guide, release 8. SAS Institure Inc, Cary NC. Sit, V., and Poulin-Costello, M. 1994. Catalogue of curves for curve fitting. Biometrics Information Handbook 4. BC Ministry of Forests, Victoria BC. Rep. 110p. Smith, D.M., Larson, B.C., Kelty, M.J. and Ashton, P.M. 1997. The practice of silviculture: applied ecology. 9th Ed. John Wiley & Sons, New York. 537 p. University of Alberta. 2008. Mixedwood Growth Model (MGM). URL: http://ales2.ualberta.ca/ rr/mgm/mgm.htm. Accessed August 18 2009. . 92 Chapter 5. General Discussion and Conclusions There has long been recognition that spatial pattern in regenerating forest stands has an important impact on yield. The primary aspect of this relationship is the frequency and size class distribution not of the trees themselves but of the gaps between trees. Measures of density alone are poor indicators of yield unless, as may often be the case in intensively managed forests, the frequency and size of larger canopy gaps is consistently low and the distribution of smaller gaps is relatively uniform. A large percentage of forest growth and yield models assume this condition, and use stand density as a major descriptor of initial condition. As a result, crude adjustment multipliers are needed to “correct” predictions for stands with less than perfect stocking. The appropriate characterization of species mixtures is also important. This trait in forest stands may be critical in cases where the species present have very different competitive strategies and economic values (such as in boreal spruce-aspen mixes), but less so where various species grow with similar patterns of growing space usage (Kelty 1992) and have more similar values (such as with many temperate conifer forests). In cases similar to the former example, it is important to quantify the growing space allocation to the various species, and to know how that allocation will affect future outcomes such as canopy pattern and harvest volumes by species with respect to management objectives (e.g. faunal habitat quality and timber revenues). Foresters have been collecting suitable data to address stocking and growing space allocation issues for most of the last century through the commonly applied practice of forest regeneration surveys. However, direct links to the effects on long term outcomes such as yield by species have been tenuous at best. Largely missing is the ability to use stocking data directly in modeling. While jurisdictions such as British Columbia and Alberta in western Canada are taking steps with their TIPSY and GYPSY models respectively to address this problem (pers. comm. W. Bergerud, BC Ministry of Forests, December 2009; Huang et al. 2009), there are many opportunities to address alternate stand structures, species compositions or entire jurisdictions not covered by existing efforts. For example, the current version of the TIPSY model is unsuitable for mixed species stands where strong patterns of ecological succession will influence changes in species composition over the period of interest (BC Ministry of Forests 2007). Generating Stem Maps The recognition that spatial pattern is important when evaluating relationships between forest stocking and yield suggests that such patterns should also be incorporated into modeling efforts. Following on this assumption, the analysis framework I described in Chapter 1 uses a wide range of simulated stem maps as the basis for subsequent model construction. These simulated stem maps are used both as the initial condition for growth models and for simulating regeneration surveys. 93 I described in Chapter 2 a methodology for generating simulated stem maps for juvenile forest stands. I designed the system to mimic the variation in stand conditions that might be expected in actual regenerating stands at the scale of operational cutting units. The result reflects my perception of factors that are important for characterizing forest stands based on two decades of sampling regenerating stands, prescribing treatments, using growth models to project future outcomes, and watching (albeit slowly) stands grow. The software into which the system is embedded (see Appendix A) provides the user with tools to help imagine and specify appropriate stand conditions, a set of spatial analysis statistics to evaluate the outcomes, and a set of routines to simulate surveys and to export treelists for use in growth models. The approach I selected to generate stem maps uses a nonhomogeneous (or filtered) Poisson or lattice process, reflecting spatial patterns expected from natural regeneration and planting respectively. The initial set of potential points is based on maximum stand densities as specified by the user. Before points are accepted as tree positions, they are subjected to various filters that reduce densities to account for factors such as random or patch mortality in plantations, and variation in establishment densities of natural regeneration caused by variations in biophysical factors such as meso-scale terrain position, patch or linear features imposed by treatment, or larger patch features imposed by other biotic factors. Spatial deviations from the Poisson process are then achieved through adding or subtracting trees from numerous small circular areas representing clusters and (semi-) voids. In comparison to other approaches to generating stem maps, my method is closest to that of Pretzsch (1997) who also employed a nonhomogeneous Poisson process. These approaches differ primarily in their emphasis on meso-scale (5 to 200 m) versus micro-scale (0.1 to 20 m) spatial variation, with my methods putting increased importance on the former. I suspected at the outset that meso- to macro-scale variation (or non-stationarity) would be an important factor in evaluating the predictive utility of regeneration survey summaries, and I incorporated it as a primary design criteria. The ability to easily incorporate variation at this scale became an important factor in generating stands with various spatial patterns for both the conifer yield prediction study in Chapter 3 and the boreal mixedwood study in Chapter 4. Conversely, the lack of specific micro-scale control features such as that for species intermingling and minimum inter- tree distances provided by Pretzsch (1997) did not appear to be a hindrance for my studies. I contrasted my method of generating stem maps with several other recent approaches including Gibbs processes (e.g. Kokkila et al. 2002) and a method for extrapolating intensive sample data using simulated annealing and nearest neighbour summary statistics (Pommerening and Stoyan 2008). Both of these approaches offered greater certainty of achieving a simulated stand with known statistical properties, but suffer in terms of flexibility and practicality – neither 94 would be suitable for generating a large number of stands at the scale desired in this dissertation within a reasonable time frame. Broad application of these methods may be restricted in that the former requires a user with sophisticated spatial statistical skills, while the latter requires a large amount of expensive field data for each stand to be generated. Overall, the methods I developed to generate simulated stem maps appear to work well where the scale of spatial pattern that is important is greater than the mean inter-tree distance, and where strict achievement of pre-determined summary statistics is not necessary; fine control at very small scales (< 2 m) is sacrificed for computational speed. The success of the model is demonstrated by the diversity of stand structures that could be generated for the subsequent studies, described in Chapters 3 and 4, and particularly by the distinction in spatial effects that could be detected in the correlations between survey summary statistics and monoculture yields in Chapter 3. Relationships Between Stocking and Yield In Chapter 3, I demonstrated a number of characteristics for the relationship between stocking and yield that I was unable to find previously documented in peer reviewed literature. For a desired level of relative yield, different stocking thresholds are required based on the expected stand height at harvest and the spatial pattern (frequency and size distribution) of voids. I was also able to support the previous but weakly documented belief that stocking thresholds should vary by species. The need for different stocking thresholds by stand height has important implications for jurisdictions where either (i) there is a wide range of site quality or (ii) harvest age will vary considerably to satisfy different management objectives. However, in the latter case, appropriate threshold variation can only be applied where expected harvest timing for a particular stand can be identified at stand establishment. An example for the former case is lodgepole pine in British Columbia, where site indices for harvested stands range at least from 14 to 26 m (breast height age 50), resulting in top heights at culmination of mean annual increment ranging from 19 to 25 m (based on TIPSY yield tables). Extrapolating from Figure 3-8, in order to achieve a minimum relative yield of 0.8, a stocking threshold (based on 16 m2 quadrats) of approximately 70% would be required for the better sites, where a value closer to 85% would be required for the poorer sites. Early in the 20th century, measures of stocking replaced measures of density in assessing regeneration success. This was largely based on the recognition that density measures were unable to account for the large degree of yield variation that would result from variations in 95 spatial pattern. The characterization of under-utilized growing space using stocking measures has been assumed ever since to account for this problem, but that assumption had never been adequately tested. I demonstrated in Chapter 3 that some effects of spatial pattern will frequently remain unaccounted for, but the magnitude of this problem will vary considerably based on selection of sampling parameters. Failure to account for this factor can lead to considerable bias in relative yield predictions for specific spatial pattern strata (i.e. planted versus naturally regenerated stands), and to reduced overall precision in predicting relative yield. Such problems will be particularly evident in jurisdictions where a wide range of regeneration spatial patterns might be expected. In setting up the species and spatial pattern contrasts for Chapter 3, I deliberately chose two conifer species with very different silvical characteristics, and assigned to each a very different site index. As this study was much closer to an inductive rather than a deductive exercise, my intent was to elicit previously unknown or poorly understood effects. From a science perspective, this approach provides many ideas to further test the bounds of the relationships that became apparent. Future analyses can now focus on relationships for additional species and species combinations, and the fluctuations in these relationships for various spatial patterns and reference heights (ages) for evaluation. From a management perspective, the relationships I developed in Chapter 3 provide a considerable increment to our collective understanding of the reliability and predictive value (for yield) of stocking estimates. There had previously been some indications of the shapes of these relationships for specific species and survey parameters (e.g. Braathe 1982), but options for improvement based on alternate survey parameters and uncertainty induced by spatial pattern had been largely unexplored. The results of this study should provide considerable guidance for investment in targeted research to further explore these issues. My current study also provides some hints on limits to the growth models used in the case studies. For example, I suggested that a relative yield of 1.0 may be achieved at a stocking level of less than 100% (using conventional quadrat stocking) for situations where soil water is the primary limiting factor. In these cases where a closed canopy will not normally develop, stocking measures based on closed canopy assumptions may start to develop unforeseen weaknesses. Testing them likely requires a growth model that accounts for below-ground competitive effects; the TASS model used in this study does not explicitly do that. The same is also true for most other forest growth models. Simplified Yield Models for Spruce-Aspen Mixedwoods In Chapter 4, I presented the results of applying the analysis framework from Chapter 1 to boreal mixedwoods, a slightly more complex problem than was done for the single species 96 conifer stands in Chapter 3. Intimate mixtures of spruce and aspen in such mixedwoods undergo potentially dramatic changes in species composition as a result of ecologic succession. For example, a stand with a composition of 95% aspen and 5% spruce (by species or basal area) at age 10 might reasonably be expected to produce 25 to 40% (or more) of its volume as spruce by age 80. Given the very different economic values (spruce >> aspen) and contributions to ecologic function of these species, it is important to account for such patterns in forest planning. I developed and contrasted several sets of models for predicting yield at age 80 (a fixed arbitrary reference age as an approximation of harvest age) by species from regeneration survey summary statistics. An important additional parameter in these models was site index, as this factor has a large influence not only on total yield but also on successional dynamics. In models that incorporated site index of both spruce and aspen as continuous variables, small changes in the relative magnitude of these variables had a large influence on predicted species composition at age 80. At the low end of the site index range, growth of spruce was generally favoured and a small reduction in aspen site index could result in complete eradication of that species from the simulated stands. At the high end of the site index range, the reverse appeared to be true with aspen gaining a complete competitive advantage with a small reduction in spruce site index. In the middle of the range there were equally dramatic shifts, but rarely with complete loss of either species. The large influence of site index on predicted outcomes played a major role in recommending one set of the models over the others for management application. A major problem in application of the models is the general lack of reliable site index for aspen, a species for which age is notoriously difficult to estimate. Given this problem, I elected to recommend against using a model that required aspen site index, preferring instead a model set that assumes a fixed linear relationship between spruce and aspen and separate models for each discrete level of spruce site index (in 1 m increments). In making this choice, I traded off a potentially large and potentially biased measurement error for a large but more quantifiable prediction error as a preference for accuracy over precision. The separate models for discrete levels of site index allowed for more flexible variation in curve shape, which in turn led to better predictions of species composition at the very high and very low ends of the site index range. Development of the simplified models for predicting yield and species composition in spruce- aspen mixedwood stands was my original motivation for undertaking the research described in this thesis. These models address a real need for management tools that has been expressed by government and licensee foresters in northeastern BC since they first began attempting to manage mixedwood forest types on an extensive basis in the early 1990’s. Based on this need (and as a contractual obligation of my research funding), I have adapted my recommended set of 97 models as the core components of a set of mixedwood regeneration standards for northeastern BC, accompanied by a regeneration survey system for collecting appropriate input data (see Appendix B). In developing the regeneration standards and associated survey system, there were several data and modeling parameters that needed to be addressed that were not factored into the simplified models from Chapter 4. These primarily included the impacts of non-crop competing vegetation and individual tree health, where assumptions or adjustments were required outside of the research outcomes for full implemention. Non-crop competing vegetation was thought to be unlikely to have any overtopping effect on the aspen trees, and no adjustments were necessary. However, there was believed to be significant risk associated with spruce, where overtopping non-tree vegetation could slow growth rates relative to those predicted by the models. It was assumed that there would be a 10% growth impact on affected spruce, and the growth effect was approximated by stratifying out the affected stocked quadrat plots, and reducing the spruce stocking in those plots (10% of affected mixed plots would become aspen plots, and 10% of spruce plots would become unstocked). This crude adjustment was based on best guess estimates of the growth effect and the need for a very simple scheme to implement it. The ideal situation would be to have non-crop vegetation explicitly included in forest growth models, thus eliminating the need for adjustment factors. Few growth models are approaching this level of sophistication. The alternative is to develop better estimates of growth impacts imposed by various non-crop vegetation species which would help provide a higher degree of confidence in any applied adjustment factors. In developing the predictive models, it was also assumed that all trees were healthy and free of defects that would impede their growth. In reality, there may be many trees that fail this assumption. In operational surveys, the typical approach to this issue is to stratify the population into acceptable and unacceptable trees based on known or suspected risk factors, and to essentially assume that the unacceptable trees do not exist. This strategy will likely lead to an underestimate of stocking, as some of these trees will recover and contribute to yield. Continued efforts to identify and monitor forest pest effects on growth, and to incorporate these factors into forest growth models will help alleviate this problem in the long run. Additional Future Research Directions Use of Distance-Dependent versus Distance-Independent Growth models One issue that I raised briefly in the discussion for Chapter 4 that deserves further exploration is the general utility of distance-dependent versus distance-independent models for the type of 98 analyses included in this dissertation. With distance-dependent models such as TASS, spatial variation within a stand is explicitly accounted for by mapping the planar Cartesian coordinates of each tree and evaluating the competitive neighbourhood accordingly. With distance- independent models such as MGM, a sample of trees from a given plot is used to represent a stand, and the distances between individual trees is not known. Stand components, as represented by the individual sampled trees, are assumed to be unclustered at stand level scales. Competitive conditions for individual trees are not evaluated based on immediate mapped neighbours, but instead based on overall density measures represented in the entire sample and on social rank by tree size. In Chapter 4, it was assumed that variation in spatial pattern within a stand could be accounted for using a distance-independent model by sampling many local stand conditions, simulating growth for each sample independently, and finally combining the results into a composite yield curve. The risk with this type of approach is that the effects of spatial scales for distances less than that of the sample plot size are lost – all trees within a sample plot are essentially assumed to be uniformly distributed. In designing a sample scheme, the analyst must compromise between inexpensive small plots and larger plots that are more likely to have enough trees for a representative treelist. There are no hard and fast rules for a minimum number of trees, but model developers I have talked with often recommend a minimum of approximately 10 trees, with an ideal target of at least 20 to 30 trees. The issue of influence trees from just outside a sample plot should also be raised, such as a cluster of large trees that falls just outside of a plot. In reality, the competitive influence of those trees might be quite large. One solution to this problem is found in the Parallel Processing Extension to the Forest Vegetation Simulator or FVS (Crookston and Stage 1993). In this case, plots are arranged as hexagonal tiles, with large trees in one tile affecting smaller trees in adjacent tiles. For stands of relatively uniform height, such effects are assumed to be minimal and have no effects on growth projections, but they can play a much larger role in the simulations of stands where tree sizes are distinctly clustered. The overall effect of such issues on the outcomes portrayed in Chapter 4 continues as a source of uncertainty. Further work on quantifying the risks and potential errors associated with the application of distance-dependent models in spatially variable stands is warranted. Other Management Objectives While this dissertation has largely focused on yield and to some extent species composition as management objectives, there are increasingly other objectives that need to be addressed. Restricting the discussion to even-aged stands, there are a wide range of forest values that are 99 impacted by combinations of species composition and spatial pattern. Perhaps the most critical of these is habitat for a wide range of faunal species, where species composition and the presence or absence of varying sized canopy gaps can influence thermal regimes, understorey development and the resulting influences on hiding cover, food sources (seeds, browse, insect variety and abundance etc.), perches, snow interception and other factors. The same factors can have large influences on snow pack management, visual quality, First Nations traditional food sources and other issues generally related to ecosystem health and biodiversity. If these values are to be actively pursued as valid objectives for forest management, they will need to be explicitly recognized in the same stocking standards that are used to ensure regenerating stands are suitable for timber production. While most timber-oriented stocking standards have ignored scale and focused strictly at minimizing unstocked gaps at the stand level, many other values will require a multi-scalar approach. For example, in conifer-dominated landscapes where hardwood stocking is generally low, the contribution of the hardwoods to values such as habitat can be disproportionately large relative to their overall abundance by adding considerable niche diversity (Bunnell et al. 1999). The contribution of such hardwood components to different values can be highly varied, and a one-size-fits-all approach to densities and spatial pattern may be highly inappropriate. Such cases may benefit from a landscape level perspective, where the spatial pattern and density of the hardwood component is intentionally varied (from low density, highly dispersed to higher density and highly aggregated). The same approach could be taken for management of canopy gap distribution objectives. Such strategies are a logical extension of the landscape level approach to species composition management that was partly addressed in Chapter 4, and the proposed application of the Chapter 4 results to the development of a survey system and regeneration standards for boreal mixedwoods in Appendix B. Closing Remarks The major contribution of this dissertation to forest management is a greater understanding of how to develop, interpret and apply regeneration standards. While this dissertation is not the first study to provide explicit links from stocking measures to future forest objectives, it is certainly the most comprehensive to date. Even so, the results presented here barely scratch the surface of a very large topic, particularly given the diversity of species, environmental conditions and management objectives that exist in North American forests in particular but also around the globe. In the current era of sustainable forest management, foresters are challenged with developing and managing for a clear vision of future forests that satisfy a wide range of values. The ability to ensure stands are contributing to those goals at an early stage of development is a critical logistical step in meeting that challenge. 100 References BC Ministry of Forests. 2007. TIPSY 4.1 help – prorating multiple species. Help files embedded in TIPSY version 4.1 software. BC Ministry of Forest and Ramsoft Systems Ltd., Victoria BC. BC Ministry of Forests. 2009. TIPSY features/functions. URL: http://www.for.gov.bc.ca/HRE/ gymodels/TIPSY/features.htm. Accessed Jan. 16, 2010. Bunnell, F.L., Kremsater, L.L. and Wind, E. 1999. Managing to sustain vertebrate richness in forests of the Pacific Northwest: relationships within stands. Environmental Reviews 7: 97-146. Crookston, N.L. and Stage, A.R. 1993. User’s guide to the parallel processing extension of the Prognosis model. Gen. Tech. Rep. GTR-INT-281. USDA Forest Service, Intermountain Research Station, Ogden UT. 13 p. Huang, S., Meng, S.X. and Yang, Y. 2009. A growth and yield projection system (GYPSY) for natural and post-harvest stands in Alberta. Alberta Sustainable Resource Development, Edmonton AB. 22p. Jozsa, L.A. and Middleton, G.R. 1994. A discussion of wood quality attributes and their practical implications. Special Publication SP-34. Forintek Canada Corp., Vancouver. 42 p. Kelty, M.J. 1992. Comparative productivity of monocultures and mixed-species stands. In: The Ecology and Silviculture of Mixed-Species Forests: A Festschrift for David M. Smith. Kluwer Academic Pub, Dordrecht, Netherlands. pp. 125-141. Appendix A Silviculture Survey Simulator User Guide 101 Appendix A - Silviculture Survey User Guide Craig Farnden Note: Minor formatting changes have been made to this document from a previous version released in March 2009 Conditions This document and the associated Excel workbook are provided free of charge on an as-is basis. Usage implies acceptance of the following conditions: 1. No warranty can be assumed. 2. Users are solely responsible for the outcomes. 3. All modifications to the macro code and fit statistics for local calibrations must be made publicly available on a free-of-charge basis. 4. All reports and publications using outcomes from the workbook must include an explicit and forthright acknowledgement of usage. Acknowledgements The Silviculture Survey Simulator was developed as part of a PhD project at the University of British Columbia. Funding for the project was provided by the BC Government’s Forest Investment Account through a delivery agreement with Canadian Forest Products Ltd., Fort St. John. Gratefully acknowledged are the patience, guidance and advice of my supervisory committee (Bruce Larson, Peter Marshall and John Nelson), staff at the BC Ministry of Forests and Range - Forest Practices Branch, and staff of Canadian Forest Products Ltd, BC Timber Sales and the BC Ministry of Forests in the Peace River region of British Columbia Craig Farnden, March 2009 105 Silviculture Survey Simulator - User Guide Introduction The Silviculture Survey Simulator is an MS Excel workbook that is used to generate realistic stem maps of juvenile stands, and to simulate silviculture surveys on those stands. The simulator is designed to have the following attributes: • it will operate at the scale of typical harvest openings, • it will have the capability to include a wide diversity of commonly used and proposed survey plot options, sample designs and methods for evaluating individual crop tree suitability; it will be flexible to produce a wide range of stand level metrics for a diversity of management objectives, • it will have the capability to encompass a wide range of tree species and spatial diversity, and • it will allow for effects of brush species, terrain and ecosystem units on stocking and competitive position of individual trees. In its current format (Version 4.00), the simulator is parameterized to generate young stands of spruce and aspen based on conditions commonly found in northeastern British Columbia, and to simulate surveys that are commonly used or proposed for use in that area. The simulator can be functionally divided into several modules. These include: • a set of functions for generating stem maps, • a set of routines for implementing silviculture surveys on simulated tree lists, • a set of routines for exporting treelists to be used as initial conditions for growth modeling (MGM1 is the only model currently supported), • a set of routines for generating spatial statistics for the generated stem map (Ripley’s K and L, and Spacing Factor at grid-based sample points are currently supported) , and • a routine for displaying the stem map as an x-y scatter chart, including the locations of survey plots, MGM treelist plots and Ripley’s K sample plots. Stand structures in the model are controlled using a combination of multi-scale horizontal complexity layers (micro- and meso-site variation), species composition, and vertical stratification. In various combinations, these features can be used to generate a wide variety of structural types. System Requirements The minimum recommended system requirements to run the simulator are a Pentium 4 processor (or equivalent) and 2 GB of RAM. The simulator will run much faster on a computer with a dual core (or better) processor. The simulator has been run extensively on Excel 2002 and 2003 on Windows XP. Minimal testing has occurred using either Excel 2007 or Windows Vista. One known issue exists when using Excel 2007, related to the drawing of scatter plots with large numbers of data points. At the time of writing, a “Hot Fix” is available from Microsoft to fix this 1. Mixedwood Growth Model, produced at the University of Alberta: http://www.ales2.ualberta.ca/rr/mgm/mgm.htm (as accessed Feb 24, 2009). 106 Silviculture Survey Simulator - User Guide issue. A discussion paper and links to the downloadable patch can be found on the Microsoft support page by searching the reference KB938538. Note that at a future date this patch may also be incorporated into a Microsoft Office Service pack. Creating Stands Stands are created in the survey simulator by generating x-y coordinates and tree attributes such as species and size. Within this process, the user is responsible for defining two sets of parameters: one describing the underlying ecosystem structure and horizontal variability, and another describing various cohorts of trees defined by species and size. These parameters in combination can be used to generate a very wide range of structural conditions. Stem maps are generated using a filtered Poisson process (Figure A-1). At the highest level of structure the process is similar to that used by Pretzsch (1997), although the details of generating tree locations differ considerably. Various combinations of filters, called complexity layers within the software, are used to generate variations in tree patterns resulting from factors such as terrain position and aspect, operational history, microsite, seed and sucker patterns, and other site features. Each cohort in the stand (e.g. an aspen overstorey resulting from post-harvest suckering versus a planted spruce understorey) is generated using a separate “rain” of points, often using its own set of distributional rules or filters. Figure A-1. Stem maps are created using a random distribution of points that are filtered to create desired stand structures. A reasonable analog for this process is a rain of seedlings onto a surface, where the falling seedlings are partially intercepted by a variably permeable screen. In the illustrated case, two small clonal patches of aspen are created using regions in the filter with high permeability, and scattered individuals in the surrounding area are created using a region in the filter with low permeability. . 0% 100% Filter Permeability Aspen Filter Random Seedling “Rain” Simulated Stand 107 Silviculture Survey Simulator - User Guide Horizontal Complexity Layers To generate stem maps that reflect species and spatial diversity at a cutblock scale, 7 layers of information are utilized, each reflecting a different pattern or scale of variation: 1. Terrain Position Terrain position is the primary layer of information underlying all others and is used to control the influence of all others. Terrain positions are mapped using Voronoi polygons (Figure A-2). Points representing terrain summits within the stem map and a 200 m surround are assigned randomly or at user selected positions. Polygon boundaries representing valleys or draws are defined using lines whose points are equidistant to the nearest two summit points. Figure A-2. Voronoi polygons as the basis for terrain units. Voronoi polygons are based on random points within the stem mapped area (dashed line) and a surrounding zone of influence. Red dots (polygon centroids) represent terrain summits, and polygon boundaries (lines for which all points are equidistant to the two nearest centroids) represent valleys (troughs) between. Figure A-3. Derivation of Terrain positions. Within a Voronoi polygon, any radius can be subdivided to generated concentric zones reflecting terrain position (left), based on user defined radial percentages. The result is a zonation of the stem map (right), where points representing tree locations are colour coded by terrain position. 1 2 3 4 5 6 Polygon centroid = meso-terrain summit Polygon boundary= meso-terrain trough Terrain Positions 108 Silviculture Survey Simulator - User Guide Voronoi polygons can be further subdivided to reflect slope position (Figure A-3). Any point within a polygon can be classified based on the relative radial distance from the summit to the polygon boundary. Six terrain positions are recognized. While this zonation is intended to mimic the classification of crest, upper, mid, lower, toe and depressional slope positions, it can also be used for other purposes such as defining large scale clusters of trees with graduated edges. 2. Affected Aspect Within terrain position classes, portions of the Voronoi polygons having an aspect within a specified range can be identified and treated as a separate stratum (Figure A-4). Any of the stocking variations identified for the main portion of the polygon can be adjusted independently. Figure A-4. Stem map with reduced stocking on a south aspect. 3. Large Patches with Matrix - 1º Horizontal Structure Similar to terrain position, large-scale patches and linear features can be defined using Voronoi polygons. Instead of six zones, this layer uses only two. Polygon centers can be used as large-scale patches (to add or eliminate trees), and linear features (matrix) can be identified using the polygon boundaries (Figure A-5). Patch size versus linear feature width is controlled by specifying the percentage of the total area covered by each feature type. 109 Silviculture Survey Simulator - User Guide Figure A-5. Stem map illustrating use of Voronoi polygons to define linear features. In this case, the linear features have a different stand composition than the remainder of the stand. Note that parameters can be specified such that effects of linear features on composition can be varied by terrain position. 4. Large Discrete Patches - 2º Horizontal Structure Similar to terrain position and the large patches with matrix, large discrete patches are generated using Voronoi polygons. In this case, however, once the polygons have been divided into two zones, only the central zone is utilized. The result is a set of discrete patches (Figure A-6) that are non-overlapping. Large patches are defined based on a mean number of centers in a 10 ha area and the percentage of the total area covered. Figure A-6. Use of large patches for spatial diversity. Large patches can be superimposed over a stem map, and trees added or subtracted by species. In this case, the dense clumps might represent clonal aspen patches in a conifer plantation. If desired, the occurence of one species can be reduced while another is increased. 110 Silviculture Survey Simulator - User Guide 5. Small Scale Patches (3 layers) - 3º Horizontal Structure Small-scale patches are circular units with randomly located centers. Patch radii are drawn from a truncated normal distribution with the mean and standard deviation specified by the user, and tails cut off at +/- 3 standard deviations. These patches are used to add and/or subtract trees from a simulated stand to create spatial diversity (Figure A-7). Within each of the 3 small-scale layers, patches can overlap, but the effects of overlap within a discrete layer are not additive. Figure A-7. Use of small-scale patches for spatial diversity. Randomly spaced trees (left) have been re-arranged (right) to create small-scale spatial diversity. In this case, trees removed from one set of patches were replaced in another set with different random centers. The seven horizontal complexity layers are controlled through an input form on an Excel worksheet (Figure A-8). Included on this form are entries to specify cutblock size, relative area occupied by each terrain position, definition of the “affected aspect”, and size/distribution of the various horizontal structure features. A final segment of the form (lower right) allows the effect of any of the horizontal structure layers to be partially or wholly attenuated by terrain position (i.e. if 4000 trees/ha are added to the stand in large discrete patches, this can be reduced to a lower number on lower slopes and in depressions to reflect poorer suckering success on wetter sites). 111 Silviculture Survey Simulator - User Guide Figure A-8. Input form for specifying horizontal variability features. Stand Size: set x and y scales for rectangular stand; area is calculated automatically. Parameters for controlling terrain position mapping: The first field controls the density of Voronoi polygons, with the polygon centroids referred to as nodes. Voronoi polygons are initially split into upland (1, 2 and 3) and lowland (4, 5 and 6) terrain positions using the percentage of the total area in each as the controlling factor. Further subdivisions are based on the width of the band (see coloured zones in Fig. 9) making up each terrain position. Values for terrain positions 3 and 6 are calculated automatically based on entries for the other positions. Aspect Limits: Where there is a certain portion of a landscape defined by aspect that requires alternate stocking conditions, the limits of that aspect are entered here. The affected area will start at the minimum values and progress clockwise to the maximum value. Parameters for controlling large patches with matrix: The first field controls the density of Voronoi polygons, while the second specifies the percent of the area that occurs in the large patches. The matrix area between patches is the remainder. An estimate of the mean patch area is calculated automatically. Parameters for controlling large patches with no matrix: The first field controls the density of Voronoi polygons, while the second specifies the percent of the area that occurs in the large patches. An estimate of the mean patch area is calculated automatically. Parameters for controlling small scale horizontal structure: deviations from a random horizontal structure at small scales (typically less than 20 m) are controlled by adding and removing trees in small circular areas. The density and size distribution of the circles is determined by the user. A particularly useful strategy is to remove a certain percentage of trees from one set of circles and replacing them in a separate set with the same density and size distribution. Cluster occurrence: As influenced by the underlying terrain positions, all other horizontal complexity layers can be turned on or off to gain greater control of their influence on horizontal structure. Partial effects can be obtained using decimal values between 1 and zero. For example, aspen clonal patches can be specified to occur using 2º Horizontal Structure parameters, with full stocking where the patch overlaps on upland terrain positions and decreasing stocking at successively lower terrain positions (a value of 0.4 will result in 40% of the trees that would have occurred given a value of 1). Use Existing Nodes: random terrain and cluster locations can be kept for multiple stem maps – a value of 1 will result in previous locations being re-used, while a value of zero will result in new locations being generated. Stand Creation Parameters 3˚ Horizontal Structure - A (small, no matrix) 200 Maximum "X" Distance (m) 4.00 ha 200 Maximum "Y" Distance (m) 25 Number of nodes per ha 8 Mean Radius (m) Terrain Position (Ecosystems, Terrain) 2.5 Radius Std Dev* (m) 10 Mean number of nodes per 10 ha 3˚ Horizontal Structure - B (small, no matrix) 80% Percent area in Upland (1, 2 & 3) versus Lowland (4, 5 & 6) 20.0% Percent radius of Upland in Core Type (Position 1) 25 Number of nodes per ha 50.0% Percent radius of Upland in Mid Type (Position 2) 8 Mean Radius (m) 30.0% Percent radius of Upland in Edge Type (Position 3) 2.5 Radius Std Dev* (m) 50.0% Percent radius of Lowland in Edge Type (Position 4) 10.0% Percent radius of Lowland in Mid Type (Position 5) 3˚ Horizontal Structure - C (small, no matrix) 40.0% Percent radius of Lowland in Core Type (Position 6) 25 Number of nodes per ha Affected Aspect 10 Mean Radius (m) 3 Radius Std Dev* (m) 135 Min Compass Bearing Limits for Aspect Effect 225 Max Compass Bearing Limits for Aspect Effect 1˚ Horizontal Structure (large; with matrix) 2 Number of Nodes per 10 ha 95% Percent area in Cluster versus Matrix Cluster Occurrence* by Terrain Position 4.75 Mean Cluster Area 1 2 3 4 5 6 1˚ 0 0 0 0 0 0 2˚ Horizontal Structure (large, no matrix) 2˚ 0 0 0 0 0 0 3˚ - A 1 1 1 1 1 1 1 Number of nodes per 10 ha 3˚ - B 1 1 1 1 1 1 3% Percent area in Clusters 3˚ - C 1 1 1 0 0 0 0.3 Mean Cluster Area * use 1 for occurs; 0 for doesn't occur; gradations between for partial effects 0 Use existing Nodes (1 = Yes, 0 = No) * Cluster Radii will be constrained to 3 standard deviations and a minimum of 1 m Å Ç É Ñ Ö Ü á à 112 Silviculture Survey Simulator - User Guide Vertical Structure and Species Composition Species composition and vertical structure are also controlled through a series of input forms. Separate forms are used to specify planted or other layers with a semi-regular distribution (Figure A-9) versus those with a random or clumped distribution (Figure A-10) typical of natural regeneration. In both cases, ten records (rows) are available, which can be used either for different species, different species specifications in different terrain positions and/or different layers or cohorts within a species. For multiple layers, more than one row can be used for each species with different tree size, density and spatial criteria. For planted trees, the entire stand is assumed to be planted to the same density, and is apportioned by species using percentages by terrain position and aspect. Adjustments to density due to mortality can be specified in a random pattern by terrain position or in a clumped pattern using complexity layers. For natural trees, density is specified as trees/ha by species, terrain position and aspect, and is initially specified by terrain position and aspect. Further adjustments are made by adding or subtracting trees by terrain position, aspect and/or complexity layers. 113 Silviculture Survey Simulator - User Guide Figure A-9. Input form for planted trees with regular spacing. Each row in this table represents a different stand cohort, broken down primarily by species but potentially also by size or age. It is assumed that the same total density (trees/ha as specified in top left corner) will be initially established throughout the stem map, but that species composition and post-planting mortality can vary based on various combinations of horizontal complexity layers. Species: Any species codes as a text string can be used; codes will be carried through to the tree list and survey reports. Preference Code: This numerical value is intended for use in silviculture surveys to allow preferential selection of one species over another where two adjacent trees can be selected for a stocking assessment and reporting of overall species composition by stocking. Height: For each layer, specify the mean height and a standard deviation. Height distributions are assumed to be normal. Percent of Total Trees by Terrain Position (and Aspect): Species composition of planted trees must be specified separately for each combination of terrain position and aspect class. In the example above, the species composition is 50% white spruce (SW) and 50% lodgepole pine (PL) for terrain position 1 over most of the area, but changes to 20% spruce and 80% pine on the “affected” aspect. Each column in this section must sum to 100%. Age and Best Years To Breast Height: The age for each layer should include the age of trees at the time of planting (e.g. 2-year old seedling + 8 years since planting = 10 years). Best years to breast height is the total age of the tallest seedlings upon reaching a height of 1.3 m assuming no overtopping suppression. Random Mortality by Terrain Position: Mortality expected to have occurred from the time of planting to the current age of the simulated stand must be specified separately for each combination of terrain position and aspect class. In the example above, 5% of the spruce and 10% of the pine is expected to have died over most of terrain position 1, switching to no pine mortality in the “affected aspect”. Additional Random Mortality by Cluster Type: Mortality in addition to that specified by terrain position can be specified for each of the 1º, 2º, and 3º horizontal complexity layers. For planted trees, this section can be used to generate effects such as losses in stocking due to unstockable areas (e.g. rock and swamp), intense competition from brush patches or root diseases. It can also be used in situations where planting would not have occurred in the presence of patches of advance regeneration. In this latter case, the ‘mortality’ for the planted trees would be 100%, and the same horizontal complexity layer would also be used within the “Stand - Natural Trees” worksheet to create the advance regeneration. Height Affected by Overtopping Cohorts?: Where the height growth of planted trees is likely to be affected by patches of overtopping trees from other layers, such effects can be triggered by entering a value of 1 for the appropriate layer in this section. Correspondingly, calibration parameters for a height reduction model will also be required (see Appendix A-2). Planted Trees 1200 Total Number of Planted Trees per Hectare 20% Spacing Tolerance (percent of mean inter-tree distance) Species Mean StdDev 1 2 3 4 5 6 1 2 3 4 5 6 101 SW 1 3 0.59 50% 50% 50% 80% 80% 80% 20% 20% 20% 80% 80% 80% 9 6 102 PL 1 4.3 0.42 40% 50% 50% 20% 20% 20% 80% 80% 80% 20% 20% 20% 9 4 103 104 105 106 107 108 109 110 Total: 100% 100% 100% 100% 100% 100% 100% 100% 100% 100% 100% 100% 2 Number of Species Species 1 2 3 4 5 6 1 2 3 4 5 6 (1=Yes, 0=No) 101 SW 5% 5% 5% 5% 5% 5% 5% 5% 5% 5% 5% 5% 0% 0% 0% 0% 0% 0% 1 102 PL 10% 10% 10% 10% 10% 10% 0% 0% 0% 0% 0% 0% 0% 0% 0% 0% 0% 0% 1 103 0 104 0 105 0 106 0 107 0 108 0 109 0 110 0 Age Best Years to BH Height affected by overtopping cohorts?3 ̊- C Clusters Random Mortality by Terrain Position (%) Random Mortality by Terrain Position Within Affected Aspect (%) 1 ˚ Clusters 1 ˚ Matrix 3 ̊- A Clusters 3 ̊- B Clusters Additional Random Mortality By: 2 ˚ Clusters Preference Code Height Percent of Total Trees by Terrain Position Within Affected Aspect (each column must sum to 100%) Percent of Total Trees by Terrian Position (each column must sum to 100%) 114 Silviculture Survey Simulator - User Guide Figure A-10. Input form for naturally regenerated trees with variations on a random distribution. Each row in this table represents a different stand component, broke down primarily by species but also by size or age. Unlike the situation for planted trees, the density is not constant across all terrain positions and aspect classes, but instead is the sum of the values specified for each layer. Species: Any species codes as a text string can be used; codes will be carried through to the tree list and survey reports. Preference Code: This numerical value is intended for use in silviculture surveys to allow preferential selection of one species over another where two adjacent trees can be selected for a stocking assessment and reporting of overall species composition by stocking. Height: For each layer, specify the mean height and a standard deviation. Height distributions are assumed to be normal. Total Trees by Terrain Position (and Aspect): Species composition of natural trees must be specified separately for each combination of terrain position and aspect class. In the example above, the overall density (trees/ha) declines from upland terrain positions moving to the wetter lowland terrain positions. At the same time, the prevalence of aspen (AT) declines while the prevalence of cottonwood (AC) increases. Also in this example, there is no change specified by aspect class. Percentage Stand Density Decrease By Cluster Type: The initial stand density specified by terrain position and aspect class can be altered based on overlaps with the horizontal structure features. For example, a narrow matrix for the 1º horizontal structure clusters can be used to simulate reduced natural regeneration on skid trails. Additional Trees By Cluster Type: The initial stand densities specified by terrain position and aspect class can also be augmented based on overlap with horizontal complexity features. This ability is useful for creating discrete patches of natural regeneration. It is also extremely useful in combination with ‘Decreases By Cluster Type’ for creating fine scale (< 10 m) clustering. In the example illustrated above, 70% of the trees are removed from one set of 3º clusters and replaced in a second set. Such re-arrangement of trees from the initial Poisson distribution can be used to create a large degree of variability in spatial pattern, although considerable trial and error may be required to get the desired results. Age and Best Years To Breast Height: Specify the maximum and minimum ages that might be measured for trees in a particular layer, and the age of the fastest growing trees upon reaching a height of 1.3 m assuming no overtopping suppression. Height Affected by Overtopping Cohorts?: Where the height growth of trees is likely to be affected by patches of overtopping trees from other layers, such effects can be triggered by entering a value of 1 for the appropriate layer in this section. Correspondingly, calibration parameters for a height reduction model will also be required (see Appendix A-2). Natural Trees Species Mean StdDev 1 2 3 4 5 6 1 2 3 4 5 6 201 AT 1 3.5 0.5 10000 10000 10000 7000 5000 4000 10000 10000 10000 7000 5000 4000 202 AC 1 3.7 0.5 500 500 500 1000 2000 2000 500 500 500 1000 2000 2000 203 204 205 206 207 208 209 210 Total 10500 10500 10500 8000 7000 6000 10500 10500 10500 8000 7000 6000 2 Number of Species Additional Trees By: (% of original number by terrain position) Species Min Max (1=Yes, 0=No) 201 AT 0% 0% 0% 70% 0% 0% 0% 0% 70% 0% 6 10 4 0 202 AC 0% 0% 0% 70% 0% 0% 0% 0% 70% 0% 6 10 4 0 203 0 204 0 205 0 206 0 207 0 208 0 209 0 210 0 Height affected by overtopping cohorts?3 ˚ - B Clusters Percentage Stand Density Decrease By: 3 ˚ - C Clusters 3 ˚ - C Clusters 1 ˚ Clusters 2 ˚ Clusters 3 ˚ - A Clusters 3 ˚ - B Clusters Age Total trees/ha by Terrain Position Due to Aspect HeightPreference Code 1 ˚ Matrix 2 ˚ Clusters 3 ˚ - A Clusters Total Trees/ha by Terrain Position (distribution but not numbers adjusted by lower level clusters) Best Yrs to BH 115 Silviculture Survey Simulator - User Guide Generating the TreeList Once all of the required tree and horizontal complexity parameters have been entered, the stem map can be generated. This process is initiated using the button located at the top of the “TreeList” worksheet (Figure A-11). Generating a new stand can be a time intensive process – up to 30 minutes or more depending on stand size and computing power. A progress indicator is provided on the status bar at the bottom of the Excel window. Once generated, the treelist will be displayed on the “TreeList” worksheet. Individual trees each occupy a row, with attribute names listed across the top in worksheet row 4. A description of the attributes can be found in Appendix A-1. If more than 60,000 trees are generated, the overflow will reside on a set of hidden worksheets (labeled TreeList 2 through TreeList 10). If needed, these can be revealed using the menu command Format > Sheets > Unhide; they are not password protected. Figure A-11. Selected area from the TreeList worksheet. All of the information on this worksheet is generated automatically – no entries are required by the user. Evaluating Generated Stands Stem Maps The generated treelist can be plotted as a stem map (Fig. A-12) on the “Stem Map” worksheet. Several options exist for controlling how the stem map is displayed. Points on the map representing individual trees can be colour coded by species, cohort or terrain position. In addition to tree locations, plot locations used for calculating Ripley’s Statistics, for simulating silviculture surveys or for subsampling the stand to generate treelistsfor growth modeling can also be displayed. It is important to note, however, that only location is shown, and that plot boundaries cannot be scaled to show their true size. 116 Silviculture Survey Simulator - User Guide Figure A-12. Stem maps of simulated stands are plotted as scatter charts in Excel. Radio buttons in the top left corner of the input form allow specification of what classification scheme is used to colour code individual trees. In cases where treelists are very large, a subset of the simulated treelist can be specified by limiting the minimum and maximum values of x and y to be displayed. Where x and y values have been limited, a button is available to reset these values for displaying the entire stem map. The green box below the input form will inform users of the current status of the displayed stem map. If a new stand has been generated but the stem map has not yet been updated, an appropriate warning will be displayed in this location. Ripley’s Spatial Indices There are a wide range of indices used to evaluate the spatial pattern produced by mapped populations of plants. One of the most commonly used of these indices is Ripley’s K(t) function and its associated L(t) function (Ripley 1976, Fortin and Dale 2005). The K(t) function describes the clumpiness of a spatial pattern by evaluating the cumulative frequency distribution for the distances between all pairs of plants for each distance t. For complete spatial randomness (CSR): K(t) = πt2 so values of K(t) > πt2 indicate that most pairs of trees are closer together than suggested by CSR and the distribution is clustered for the evaluation distance t. If the opposite is true, the spatial pattern is over-dispersed for that particular scale. The K(t) function can be standardized to generate the L(t) function, where the expected value for CSR is always zero: L(t) = t - (K(t)/π)0.5 Based on this statistic, spatial pattern can be evaluated graphically across scales (Fig. A-13). Such charts are generated for up to 9 subsets of the overall stem map. 117 Silviculture Survey Simulator - User Guide Figure A-13. Examples of L(t) functions. In the top chart, trees are overdispersed (more evenly spaced than would be expected with CSR) at scales less than 1 m, possibly due to self-thinning. At scales from 1 to 12 m, the stand is clustered, and becomes overdispersed again for larger scales. In the bottom chart, the stand is overdispersed at scales out to 9, a pattern that would be typically of regularly spaced plantations. The clustering at larger scales might indicate large unstockable gaps or patches of mortality in the plantation. The Ripley spatial statistics and the L(t) charts are found on a worksheet labeled “Ripley Stats”. The top left corner of this sheet contains input forms and control buttons for generating the statistics and charts (Fig. A-14). Caution: Generating spatial statistics is a slow process – expect to wait anywhere from several minutes to half an hour, depending on the number of sub-samples, the sub-sample radius and the processor speed of your computer. Further information on Ripley’s statistics is available from a number of sources including Ripley (1976), Diggle (1983), Dale (1999), and Fortin and Dale (2005). -1 -0.5 0 0.5 1 0 5 10 15 t -1 -0.5 0 0.5 1 0 5 10 15 t L (t ) L (t ) 118 Silviculture Survey Simulator - User Guide Figure A-14. Input form for Ripley’s Statistics. Ripley’s statistics in this usage are based on samples of the population using circular plots. Plots can be randomly located or arranged on a grid, with up to 9 samples per stand. Users need to enter the sample design (grid versus random), the number of sample plots (1 to 9), the number of plots in the x direction (grid systems only), and the sample plot radius. For grid systems, the number of plots in the y direction will be calculated automatically, and the total number of plots specified must be evenly divisible by the number of plots in the x direction. The plot radius should be at least 20 m. The top green box below the input form will inform users of the current status of the displayed charts. If a new stand has been generated but the L(t) charts have not yet been updated, an appropriate warning will be displayed in this location. The bottom green box reminds users that there are two different conventions for calculating L(t) statistics, and indicates which variation has been used in this case. Spacing Factor In some cases, the relative frequency of sample points with different values for spacing factor may provide a good indication of spatial variability or, more particularly, the distribution of competitive environments within a stand. From a silviculture perspective, spacing factor can be used to indicate locations where trees are excessively crowded for achieving a given set of management objectives versus other areas that are within a desired range or are understocked. In boreal mixed stands of juvenile aspen and spruce, for example, spruce growing below an aspen canopy starts to lose height increment below spacing factors of approximately 0.25 (C. Farnden, unpublished data). To provide an estimate of the distribution of competitive conditions within simulated stands, two distribution charts are provided in Figure A-15. Sample Design Grid Number of Plots 4 No. Plots - x Direction 2 No. Plots - y Direction 2 Plot Radius (m) 25 Note: The L(t) charts follow a convention whereby K(t) is subtracted from the expected value given complete spatial randomness. Values of L(t) greater than zero indicate overdispersion, while values less than zero indicate clumpiness at the given scale of evaluation (x-axis in metres). Dispayed L(t) charts represent the current stand Generate Ripley's K Statistics 119 Silviculture Survey Simulator - User Guide Figure A-15. Histogram and cumulative distribution for Spacing Factor. The distribution and variation in competitive conditions is displayed based on 100 points in the simulated stand using a grid sampling scheme. Spacing factor is calculated using all trees in the stand: SF = mean inter-tree distance ÷ mean height. Figure A-16. Silviculture Survey Control Form Run Silviculture Survey: Use this button to initiate a survey. Survey Design: select random or grid (systematic) from the drop-down list. Random designs will select random points from the stand area conditioned to exclude an edge strip equal in width to the maximum radius of the selected survey plot. Grid Angle: not currently implemented – this option is intended to allow for systematic samples on a grid that is oriented other than at right angles to the stand edges. Grid Interva/Number f Plots: for random designs, this field requests the number of plots. For a systematic (grid) sample, this field requests a grid interval between plots. For a systematic design, the simulator will use a random start point and add plots in the x and y directions at the specified interval until no more plots will fit. Plot Type: select the survey type from the drop-down list. Plot Area/ MITD: for quadrat plot surveys, this field requires entry of the quadrat size in m2. For Mixedwood Well Spaced surveys, this field requires entry of the minimum inter-tree distance (MITD). 0 10 20 30 40 50 60 <0.1 0.1 to 0.2 0.2 to 0.3 0.3 to 0.4 0.4 to 0.5 0.5 to 0.6 0.6 to 0.7 0.7 to 0.8 0.8 to 0.9 > 0.9 Spacing Factor 0 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Spacing Factor R el at iv e F re q u en cy ( % ) C u m u la ti ve F re q u en cy ( % ) 120 Silviculture Survey Simulator - User Guide Surveys Silviculture surveys are run from the “Survey” worksheet using an input form (Fig A-16). Three survey methods designed for use in boreal mixedwood stands are currently implemented: 1. 5-5-2-2 (Boreal Mixedwood Survey) The 5-5-2-2 survey starts with a 50 m2 plot that is divided into quadrants (BC Ministry of Forests 2009, pages 112 to 117; see “variant of survey for macro, meso and micro-scale patched mixedwoods”), each of which is roughly analogous to a stocked quadrat plot. Each quadrant is evaluated as being either unstocked, stocked with broadleaved trees (aspen or cottonwood), stocked with spruce or stocked with both broadleaves and spruce. In order for a spruce to be tallied, a further assessment is required using a tree-centered subplot. The subplot is again divided into quadrants, and the spruce must be free of broadleaved competition for a distance of at least 5 m in two quadrants and 2 m in the remaining quadrants in order for the spruce tree to be “well growing” and tallied. In the BC MoF version of this survey, the subplot quadrants with a competition-free radius of 5 m must be adjacent. In its original version and that implemented in the simulator, the 5 m radii need not be in adjacent quadrants. Summary statistics for this survey are comprised of the percentage of plot quadrants that fall within each of the four stocking categories. 2. Boreal Mixedwood Stocked Quadrats Stocked quadrat surveys use small fixed area plots (or quadrats) each of which represents a sample of potentially occupiable growing space. Each quadrat is tallied as being unstocked, stocked with broadleaved trees, stocked with spruce, or stocked with both broadleaves and spruce. There is no fixed plot size for the quadrats – plot area can be set by the user. Summary statistics for this survey are comprised of the percentage of quadrats that fall within each of the four stocking categories, as for quadrants in the 5-5-2-2 survey. 3. Mixedwood Well-Spaced The mixedwood well-spaced survey uses an extension of the standard well spaced survey procedure (BC Ministry of Forests 2009, pages 23 to 27) utilized for a high percentage of regenerating stands in BC. In the mixedwood variation used here, three separate tallies of well-spaced trees are required: • Well-spaced spruce, using a minimum inter-tree distance of 2 m, ignoring the presence of broad-leaved trees • Well-spaced deciduous trees (aspen and cottonwood), ignoring the presence of spruce • Total well spaced trees, ignoring the concept of species In each tally, tree selections are made to maximize the number of well spaced trees; for the final tally of total trees this will usually require selection of different trees than were chosen for the species specific tallies. 121 Silviculture Survey Simulator - User Guide The tally of total well spaced trees is intended to primarily capture the total usage of growing space (or inversely, what percentage is unoccupied). The species specific tallies help to apportion that growing space between the spruce and broadleaves. Summary statistics for this survey include: 1) mean number of total well-spaced trees, 2) mean number of broadleaved well-spaced trees, 3) mean number of spruce well-spaced trees and 4) mean number of aspen well spaced trees associated with each well spaced spruce: where: n = number of plots p = number of well spaced spruce in plot i C = number of well spaced broadleaves in plot i This last statistic is intended as a measure of the mean level of competition experienced by the spruce in the stand as imposed by the overtopping broadleaved trees. Survey summary statistics are displayed on the “Survey” worksheet (Fig. A-17), while individual plot outcomes can be found on a hidden worksheet labeled “Survey Outcomes”. The outcomes worksheet is hidden simply to avoid clutter and is not password protected. Figure A-17. Silviculture Survey Summary Statistics The summary statistics tables can illustrate one survey outcome for each of the 5-5-2-2 and Mixedwood Well Spaced surveys, and three for the Quadrat Surveys. = = == n i n i p j n p C SwperAtMean 1 1 1 No. Plots Unstocked Quadrants Decid Quadrants Conifer Quadrants Mixed Quadrants 36 8.3% 91.7% 0.0% 0.0% Quadrat Survey 1 Quadrat Area (m 2 ): 7.50 No. Plots Unstocked Quadrants Decid Quadrants Conifer Quadrants Mixed Quadrants 100 4.0% 25.0% 13.0% 58.0% Quadrat Survey 2 Quadrat Area (m 2 ): 5.00 No. Plots Unstocked Quadrants Decid Quadrants Conifer Quadrants Mixed Quadrants 100 17.0% 36.0% 13.0% 34.0% Quadrat Survey 3 Quadrat Area (m 2 ): 10.00 No. Plots Unstocked Quadrants Decid Quadrants Conifer Quadrants Mixed Quadrants 100 0.0% 16.0% 9.0% 75.0% WS Mixedwood Survey No. Plots Mean Conif WS Mean Decid WS Mean Total WS Mean Decid WS per Conifer WS 49 5.63 6.49 7.88 5.60 5-5-2-2 Mixedwood Survey 122 Silviculture Survey Simulator - User Guide TreeLists For Growth Modeling Currently, treelists can be generated only for the MGM growth model. Treelist generation is controlled from the data entry form (Fig. A-18) on the “Growth Modeling” worksheet. Figure A-18. Growth Model TreeList Control Form Growth Model: Use this button to initiate a survey. Survey Design: select random or grid (systematic) from the drop-down list. Random designs will select random points from the stand area conditioned to exclude an edge strip equal in width to the maximum radius of the selected survey plot. Grid Angle: not currently implemented – this option is intended to allow for systematic samples on a grid that is oriented other than at right angles to the stand edges. Grid Interva/Number f Plots: for random designs, this field requests the number of plots. For a systematic (grid) sample, this field requests a grid interval between plots. For a systematic design, the simulator will use a random start point and add plots in the x and y directions at the specified interval until no more plots will fit. Plot Area: area of the sample plots (m2) used to generate individual treelists. The plot needs to be large enough to capture 10 to 50 locally competing trees. The plot should be small enough that it represents a relatively uniform set of competitive conditions within the stand. A plot area of 100 m2 is a good place to start. MGM works best with a treelist that is at least 30 to 40 trees. If a smaller treelist is created at a particular sample location, it will be multiplied to generate addition tree records until at least 40 trees are present. The plot expansion factors are automatically adjusted so that the stand level values represented by the plot (e.g. trees/ha) are accurately retained. Treelist files created in this manner are placed in a user defined directory. A form is provided on the “Growth Modeling” worksheet, where the user can select either the same directory as the survey simulator workbook, or browse for any other directory. 123 Silviculture Survey Simulator - User Guide References BC Ministry of Forests, 2009. Silviculture Survey Procedures Manual. http://www.for.gov.bc.ca/hfp/publications/00099/Surveys/Silviculture%20Survey%20Pro cedures%20Manual-April%201%202009.pdf, accessed March 29, 2009. Dale, M.R.T. 1999. Spatial pattern analysis in plant ecology. Cambridge University Press, Cambridge UK. Diggle, P.J. 1983. Statistical analysis of point patterns. Academic Press, London UK. Fortin, M-J and Dale, M. 2005, Spatial analysis, a guide for ecologists. Cambridge University Press, Cambridge UK. Pretzsch, H. 1997. Analysis and modeling of spatial stand structures. Methodological considerations based on mixed beech-larch stands in Lower Saxony. Forest Ecology and Management, 97: 237-253. 124 Silviculture Survey Simulator - User Guide Appendix A-1. Tree Attributes In order of column number, the following attributes are provided on the “TreeList” worksheet: Tree Number: A unique numerical account for each tree. X: Horizontal distance of the tree from the left margin of the stand. Y: Vertical distance of the tree from the bottom edge of the stand. Species: An alphanumeric code, specified by the user, to indicate discrete species. Height: Tree height in metres. An initial height is randomly drawn from a normal distribution where the mean and standard deviation are specified by the user. Where necessary, height is adjusted downward to reflect the effects of overtopping competition (see Appendix A-2) DBH: Bole diameter in centimetres measured at a height of 1.3 m from the ground. Diameter is calculated as a function of height and spacing factor (see Appendix A-2). Crown Width: Mean radial diameter of the tree crown in metres. Currently, crown width is generated as a crude linear function of DBH (CW = 0.31 + 0.148*DBH). Crown width has been included as an attribute to facilitate its use in surveys where it is needed to evaluate of overtopping competition. Total Age: Years since germination. For planted trees, age is specified by the user. For natural trees, age is initially estimated as a function of height, and adjusted based on a random draw from a skewed beta distribution (see Appendix A-2). BHage: Breast Height Age is the number of years since reaching a height of 1.3 m; also described as the ring count of a tree cored with an increment borer at a height of 1.3 m. This value is calculated as a function of height and total age, assuming a linear fit between height and age for trees above breast height (see Appendix A-2). SF All Trees: Spacing Factor (mean inter-tree distance as a percentage of mean height) calculated using all trees within a calculated search radius of the subject tree. The search radius is determined as 70% of the maximum tree height in the stand. SF Taller Trees: Spacing Factor (mean inter-tree distance as a percentage of mean height) calculated using all trees taller than the subject tree within a calculated search radius of the subject tree. The search radius is determined as 70% of the maximum tree height in the stand. SF Selected: Spacing Factor (mean inter-tree distance as a percentage of mean height) calculated using all trees of user selected species set within a calculated search radius of the subject tree. The search radius is determined as 70% of the maximum tree height in the stand. Relative Height: Height of the subject tree as a proportion of the height of the maximum potential height a given cohort. Terrain Position: Record of terrain position associated with subject tree. Possible values are integers 1 through 6. 1˚ Cluster Position: Record of position within primary horizontal complexity layer – value of 1 represents cluster, value of 2 represents matrix. 125 Silviculture Survey Simulator - User Guide 2˚ and 3˚ Cluster Positions: Records of position within secondary and tertiary clusters. A value of 1 indicates tree is within a given cluster set, a value of 0 indicates tree is not within a given cluster set. Affected Aspect: Record of position within affected aspect. A value of 1 indicates tree lies within the aspect range defined as the “affected aspect”, while a value of 0 indicates that it lies outside that range. Nearest Cluster & X & Y: Number and x-y coordinates of nearest terrain position centroid to the subject tree. 2nd Nearest Cluster & X & Y: Number and x-y coordinates of second nearest terrain position centroid to the subject tree. Survey Preference Code: Ranking factor provided by the user in the specification of stand cohorts to indicate a preference of one species over another when choosing trees for stocking assessments in silviculture surveys. Competition Map Pixel: Parameter used to subset the stem map into small grid squares for evaluation of local competition. Pixel dimensions are based on rounding up the spacing factor search radius to the nearest integer that divides evenly into the x and y stand dimensions respectively. A cluster of 9 pixels centered on the one containing the subject tree will then contain all possible trees needed to assess local competition. A similar process is used to subset the stand for survey plots and plots used to generate treelists for initiating growth model simulations (where pixel size is based on rounding up the plot radius). Cohort/Layer: Record of the cohort from which a particular tree was derived, as specified by the user. 126 Silviculture Survey Simulator - User Guide Appendix A-2: Sub-Models Tree Height Adjustments Tree heights are initially drawn from a normal distribution with mean and standard deviation specified by the user (see Figures 9 and 10 in main text). Where height is anticipated to have been affected by overtopping competition by trees from other cohorts2, height can be adjusted downward using a multiplier (range = 0 to 1). The multiplier is determined as a function of one of three variations on spacing factor, using one of eight different functional shapes and using calibration parameters as set by the user. An example of a function for the impact of an aspen canopy on understorey spruce is illustrated in Figure A-2-1. Figure A-2-1. Weibull fit of relative spruce height as predicted by spacing factor calculated using the overtopping deciduous canopy. Relative spruce height was determined by first fitting the curve using absolute spruce height, then dividing each height value by the asymptote. RMSE = 0.19 Built-in functions available to represent this type of relationship include: Weibull: where b 1 is fixed at 1.00 Chapman Richards: where b1 is fixed at 1.00 Generalized Logistic: where b 1 = b 4 = 1.00 Gompertz: where b 1 = 1.00 Segmented Linear: y = b 2 = 1 for x > b 1 y = b 2 -b 3 (b 1 -x) for x < b 1 Specifying the function to use along with the appropriate parameters is done using a form in columns X through AH of the ‘Stand Creation Parameters’ worksheet (Fig. A-2-2). 0.2 0.4 0.6 0.8 1.0 1.2 1.4 R el at iv e H ei g h t 0.10 0.15 0.20 0.25 0.30 0.35 Spacing Factor (Aspen) 228.32.3771 xey −−= y = b1 1− e −b2x b3( ) ( ) 3211 bxbeby −−= xbbeb b y 32 4 1 −+ = xbbeeby 32 1 −−= 1. Cohorts needing adjustment must be indicated on the cohort description forms - see Figures 9 & 10. 127 Silviculture Survey Simulator - User Guide Figure A-2-2. Height adjustment models can be specified separately for each cohort. The model form is specified by choosing the appropriate model Code (top left) and entering it in the Model Code column for the appropriate cohort. The form of spacing factor being utilized (all trees versus trees larger than the subject tree versus trees in selected cohorts) is specified by selecting the appropriate code in the bottom left corner and entering it in the field labeled “Ht Adj Model Uses”. Parameters for previously calibrated models can be listed below the input form, and can be copied into the input form when needed. For cases where spacing factor is to be based on selected cohorts, these can be specified in column AY of the ‘Stand Creation Parameters’ worksheet (Fig. A-2-3). Note that the same set of cohorts must be used for all sub- models using this variation of spacing factor. Figure A-2-3. Where spacing factor is to be calculated using only certain cohorts, those cohorts are specified by entering a value of 1, otherwise enter 0. Height Adjustment Parameters Cohort Model Code b1 b2 b3 b4 b5 Model Form* Model Code 101 1 1 377.2 3.2278 0 0 Weibull (3 param) 1 102 1 1 377.2 3.2278 0 0 Chapman-Richards 2 103 1 1 377.2 3.2278 0 0 Gen. Logistic 3 104 1 1 377.2 3.2278 0 0 Gompertz 4 105 1 1 377.2 3.2278 0 0 Segmented Linear 5 106 1 1 377.2 3.2278 0 0 Polynomial 6 107 1 1 377.2 3.2278 0 0 Exp Polynomial* 7 108 1 1 377.2 3.2278 0 0 Linear 8 109 1 1 377.2 3.2278 0 0 * Hover cursor over model name for format 110 1 1 377.2 3.2278 0 0 Measure of Competition Used in Predictions: 201 1 0 0 0 0 0 Spacing Factor 1 202 1 0 0 0 0 0 Spacing Factor Taller 2 203 1 0 0 0 0 0 Spacing Factor Selected Cohorts 3 204 1 0 0 0 0 0 205 1 0 0 0 0 0 Ht Adj Model Uses: 3 206 1 0 0 0 0 0 DBH Model Uses: 3 207 1 0 0 0 0 0 208 1 0 0 0 0 0 209 1 0 0 0 0 0 210 1 0 0 0 0 0 Calibrated Models: Planted white spruce under aspen (Farnden 2009) 1 1 377.2 3.2278 0 0 Cohort 101 0 102 0 103 0 104 0 105 0 106 0 107 0 108 0 109 0 110 0 201 1 202 1 203 1 204 1 205 1 206 1 207 1 208 1 209 1 210 1 Include in Selected Cohort Calculations? 128 Silviculture Survey Simulator - User Guide Diameter Diameters are calculated as a function of height (Figs. A-2-4, A-2-5). Separate diameter versus height functions can be utilized for each cohort. Figure A-2-4. Linear fit of logarithmically transformed dbh for aspen diameters predicted by height. I2 = 0.87 RMSE = 0.67 A number of functional forms are available including: Weibull: Chapman Richards: Generalized Logistic: Gompertz: Polynomial: Exponential Polynomial: Linear: Specifying the function to use along with the appropriate parameters is done using a form in columns X, Y and AK through AU of the ‘Stand Creation Parameters’ worksheet (Fig. A-2-2). 0 1 2 3 4 5 6 7 8 9 10 D b h ( cm ) 2 3 4 5 6 7 8 9 Height (m) xey 2797.04106.0 +−= y = b1 1− e −b2x b3( ) ( ) 3211 bxbeby −−= xbbeb b y 32 4 1 −+ = xbbeeby 32 1 −−= y = b1 + b2 x + b3 x 2 + b4 x 3 + b5 x 4 4 5 3 4 2 321)log( xbxbxbxbby ++++= xbby 21 += 129 Silviculture Survey Simulator - User Guide Figure A-2-5. Diameter models can be specified separately for each cohort. The model form is specified by choosing the appropriate model Code (columns X & Y at left) and entering it in the Model Code column for the appropriate cohort. The form of spacing factor being utilized (all trees versus trees larger than the subject tree versus trees in selected cohorts) is specified by selecting the appropriate code in the bottom left corner and entering it in the field labeled “DBH Model Uses”. Parameters for previously calibrated models can be listed below the input form, and can be copied into the input form when needed. Columns for parameters b6 through b9 are not currently used, but have been provided to allow for future inclusion of random variation. Diameter Prediction Parameters Cohort Model Code b1 b2 b3 b4 b5 b6 b7 b8 b9 Model Form* Model Code 101 8 -1.503 1.4193 0 0 0 Weibull (3 param) 1 102 8 -1.503 1.4193 0 0 0 Chapman-Richards 2 103 8 -1.503 1.4193 0 0 0 Gen. Logistic 3 104 8 -1.503 1.4193 0 0 0 Gompertz 4 105 8 -1.503 1.4193 0 0 0 Segmented Linear 5 106 8 -1.503 1.4193 0 0 0 Polynomial 6 107 8 -1.503 1.4193 0 0 0 Exp Polynomial* 7 108 8 -1.503 1.4193 0 0 0 Linear 8 109 8 -1.503 1.4193 0 0 0 * Hover cursor over model name for format 110 8 -1.503 1.4193 0 0 0 Measure of Competition Used in Predictions: 201 7 -0.416 0.2797 0 0 0 Spacing Factor 1 202 7 -0.416 0.2797 0 0 0 Spacing Factor Taller 2 203 7 -0.416 0.2797 0 0 0 Spacing Factor Selected Cohorts 3 204 7 -0.416 0.2797 0 0 0 205 7 -0.416 0.2797 0 0 0 Ht Adj Model Uses: 3 206 7 -0.416 0.2797 0 0 0 DBH Model Uses: 3 207 7 -0.416 0.2797 0 0 0 208 7 -0.416 0.2797 0 0 0 209 7 -0.416 0.2797 0 0 0 210 Calibrated Models: Planted white spruce under aspen (Farnden 2009) 8 -1.503 1.4193 0 0 0 Aspen - natural sucker stand (Farnden 2009) 7 -0.4106 0.2797 0 0 0 130 Silviculture Survey Simulator - User Guide Total Age Total age for natural trees is calculated primarily as a function of height, assuming that the oldest trees are also most likely to be the tallest. The first step, then, is to interpolate age from height: where Height max , Height min , Age max and Age min are determined as 3 standard deviations from the mean for a normal distribution as specified by the user (see Figs. A-9 and A-10 in the main text). The next step adjusts this linear distribution upward to generate a distribution that is skewed to older ages (assuming that more trees germinate earlier rather than later, and that older trees are more likely to have survived early stages of self-thinning). Adjustments are based on random draws from a skewed cumulative beta distribution (Fig. A-2-6), such that most trees get a relatively small adjustment, but a few are quite large: Figure A-2-6. Cumulative beta distribution where alpha = 4 and beta = 1.5. Based on a uniform random value, a draw is made from the cumulative distribution to obtain and adjustment factor for total age. minminmax minmax min )( )( )( AgeAgeAge HeightHeight HeightHeight Ageinitial +−×− −= )( max initialdrawinitialfinal AgeAgeBetaAgeAge −+= 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Random Value C u m u la ti v e B e ta ( 1 .5 ,4 ) 131 Silviculture Survey Simulator - User Guide Breast Height Age Breast height age is calculated from total age and height, assuming that the relationship between height and age for young trees is linear above breast height. Breast height age is then calculated as: where intercept x is the x intercept of the linear fit for the age versus height relationship above breast height (Fig. A-2-7). Figure A-2-7. Age BH determination utilizes the assumption that a height versus age relationship (red line) is close to linear (black dashed line) above breast height for young trees. )( )3.1( xTotalBH InterceptAgeHeight Height Age −−= 0 1 2 3 4 5 0 5 10 15 Age (yrs) Height (m) 1.3 m Total Age Years to BH AgeBH H ei g h t (m ) 132 Silviculture Survey Simulator - User Guide Appendix A-3: Batch Runs The simulator contains no built-in routine for batch runs. It has, however, been constructed to facilitate batch runs. Some knowledge of the VBA programming language is required. There are three steps to the process: 1. Build an excel worksheet containing the parameters that need to be changed in each successive run of the simulator, with one row for each simulation and a run ID number or name in the first cell. 2. Construct a VBA macro that will copy the parameters into the appropriate locations on the simulator worksheets, run the simulator, and save the results. In pseusdo-code: Sub Batch_Macro() For i = 1 to Number_Of_Simulations Copy parameters from Row i of batch workbook to simulator, where the first column contains a Run ID Call the appropriate simulator subroutine(s) Save the simulator workbook using the Run ID from cell (i,1) Next i End Sub The simulator workbook must be open before the batch run is started, and should not be closed between successive simulations. Subroutines that may need to be called to initiate batch runs include: Create_Stand2 Equivalent to selecting the “Create New Stand” button on the “TreeList” worksheet. Run_Survey Equivalent to selecting the “Run Silviculture Survey” button on the “Survey” worksheet. Run_Modeling_TreeLists Equivalent to selecting the “Generate TreeLists For Growth Models” button on the “TreeList” worksheet Generate_Map Equivalent to selecting the “Generate Stem Map” button on the “Stem Map” worksheet Generate_Spatial_Stats Equivalent to selecting the “Generate Ripley’s K Statistics” button on the “Ripley Stats” worksheet 133 Silviculture Survey Simulator - User Guide Appendix A-4: Coding New Surveys In order to code in a new survey methodology there are several actions that you need to take with regard to worksheet locations in the simulator workbook: 1. Add the survey name to the Plot Type list in Column BA of the “Survey” worksheet. 2. Ensure that the Data Validation option for cell D10 on the “Survey” worksheet covers the range of cells that includes your entry in Column BA. 3. Determine the number of statistics that need to be recorded for each plot in the survey. 4. Unhide the Sheet labeled “Survey Outcomes” and identify a range of unused columns that is sufficiently wide enough to encompass all of your plot level statistics plus Plot Number, X-coordinate and y-coordinate (see examples already present on this worksheet). Label your columns. Note that this worksheet has been set up with many blank columns between successive survey methods to allow for future adjustments. 5. Create a block of cells on the “Survey” worksheet (starting in column G) to summarize the plot level statistics from the “Survey Outcomes” worksheet into stand level statistics as desired. Use Excel formulas to sum, average or otherwise compile the plot level statistics into stand level statistics. 6. Determine user defined plot variables (such as plot area or minimum inter-tree distance for well spaced trees) and assign then to input fields in the user form on the top left corner of the “Survey” worksheet. In doing so, note that: • The first four fields (rows 4, 6, 8 and 10) should not be changed • The fifth field (row 12) can be used for any purpose, and the field label set to change (via a cell formula) according to the choice made for Plot Type (row 10). • Additional fields can be added below the current ones as needed. The next step is to develop an algorithm for your survey and translate it into VBA code. All code for silviculture surveys is found in Module 1, and follows the flow chart illustrated in Figure A-4-1. 134 Silviculture Survey Simulator - User Guide Figure A-4-1. Flow chart for routines used to simulate silviculture surveys. There are a number of steps in this algorithm that are common to all survey methods. These include: 1. Sub Run_Survey: This subroutine controls the flow of the algorithm. The only change required is an additional case for the “Select Case” set of statements. The case number must match the numerical value from column BB of the “Survey” worksheet beside the new survey name in column BA. The statements that are run within the new case will perform 2 functions: • call the appropriate subroutine to implement the desired survey type, and • write plot summary statistics to the “Survey Outcomes” worksheet, where each plot receives one line within a specific range of columns. 2. Sub Survey_Design: This subroutine reads in stand size and general survey parameters, and assigns plot centers within the stand. The only change required is an additional case for the “Select Case” set of statements. The case number must match the numerical value from column BB of the “Survey” worksheet beside the new survey name in column BA. The statements that are run within the new case will perform 3 functions: • Clear the range of cells in which plot outcomes will be written on the “Survey Outcomes” worksheet • Assign constant values required for the plot (this may include the plot radius if it is fixed) • Read in additional plot variables that are required from column D of the “Survey” worksheet. Note: it is important to ensure that the PlotRadius variable is sufficiently large to encompass all trees that need to be evaluated, including any trees outside of the sample plot that will influence decisions within the plot. If a second or additional radii is/are required, either for the actual plot radius or for subplots, new variables will need to be assigned. Sub Survey_Design • Read in survey parameters • Create plot centers Sub Read_TreeList • Read in stem map from Excel worksheet • Create subsets of treelist by grid squares (pixels) Sub Create_Plot_TreeList • Compile subset of main treelist specific tp plot Various Survey Subroutines • Determine which trees are “in” • Compile plot statistics into array Call appropriate survey subroutine Write plot statistics to “Survey Outcomes” worksheet Count Through plots: More required? Start End Yes No 135 Silviculture Survey Simulator - User Guide 3. Sub Read_TreeList: This subroutine reads the treelist for the entire stand and subdivides it into grid rectangles (or pixels). The dimensions of each pixel are as small as possible such that the x and y dimensions are larger than the value of the PlotRadius variable, and will divide evenly into the respective x and y stand dimensions (MaxX and MaxY variables). Each pixel is assigned a vector within a 2-D array which contains a list of all trees that occur within that pixel. No changes are required in this subroutine. 4. Sub Create_Plot_TreeList: This subroutine finds the pixel that contains a plot center along with the eight surrounding pixels to form a 3 x 3 grid of pixels. The list of trees in each of the nine pixels is then used to compile a new treelist that is a subset of the main treelist and contains all of the trees that need to be evaluated for a given plot. No changes are required in this subroutine. The remaining subroutines in Module 1 are unique to specific surveys, and are called from the Run_Survey subroutine as required. These subroutines are called separately for each plot. The data available to be used by the survey method subroutines include: a) A 2-D array labeled PlotTrees containing all trees to be evaluated in the plot. Each row in the array is an individual tree, and each column is an attribute as described in Appendix A-1. The order of attributes is the same as can be found on the “TreeList” worksheet. b) A set of variables including PlotRadius and any other variables or constants that are specific to the plot type. Also available are a set of built in functions which can be found in Module 3. Of particular interest is a function that calculates the distance between two points: Dist(x 1 , y 1 ,x 2 ,y 2 ) The new survey subroutine should calculate plot level statistics (e.g. stocking class, count of trees by species, count of well spaced trees by species) based on rules imposed by the survey methodology, with the statistics stored as variables that are accessible in the Run_Survey subroutine. 136 Silviculture Survey Simulator - User Guide Appendix A-5: Coding New Modeling TreeLists In order to code in a new sub-sampling methodology for creating tree-lists to initial simulations in accessory growth and yield models, there are several actions that you need to take with regard to worksheet locations in the simulator workbook: 1. Add the survey name to the Plot Type list in Column BI of the “Growth Modeling” worksheet. 2. Ensure that the Data Validation option for cell E5 on the “Growth Modeling” worksheet covers the range of cells that includes your new entry in Column BI. 3. Ensure that the user form in the top left corner of the “Growth Modeling” worksheet contains all of the require plot parameters.. In doing so, note that: • The first four fields (rows 4, 6, 8 and 10) should not be changed • The fifth field (row 12) can be used for any purpose, and the field label set to change (via a cell formula) according to the choice made for Growth Model (row 5). • Additional fields can be added below the current ones as needed. For instance, if a nested plot is required, additional plot radii or basal area factors may be required. The next step is to develop an algorithm for your sampling scheme and translate it into VBA code. All code for growth model sub-sampling is found in Module 5, and follows the flow chart illustrated in Figure A-5-1. Figure A-5-1. Flow chart for routines used to sub- sample simulated stands to generate treelists as the initial conditions for growth modeling. Sub Model_Plots • Read in sample design parameters • Create plot centers Sub Read_TreeList_For_Models • Read in stem map from Excel worksheet • Create subsets of treelist by grid squares (pixels) Create subset treelist from pixels Determine which tree ar “in” the sample plot Count Through plots: More required? Start End Yes No Determine for which models to create treelists Various Subroutines - Start - Write trees to output file with: • Any required additional attributes • Appropriate column order • Appropriate header information 137 Silviculture Survey Simulator - User Guide There are a number of steps in this algorithm that need either no change or minimal change for adding new modeling methods. These include: 1. Sub Run_Modeling_TreeLists: This subroutine controls the flow of the algorithm. The only change required is an additional case for the Select Case set of statements. The case number must match the numerical value from column BJ of the “Survey” worksheet beside the new survey name in column BI. The only statement that is run from the Select Case set of statements is a call to another subroutine that contains algorithms for defining the contents of the sub-sample treelist. 2. Model_Plots: This subroutine reads in stand size and general sub-sampling parameters, and assigns plot centers within the stand. The only change that may be required is the addition of statements to read in any additional plot parameters that may have been added. These statements are near the top of the subroutine following a comment that reads Read in Survey Design Parameters. Note: The variable PlotRadius is calculated from a user defined plot area. 3. Sub Read_TreeList_For_Models: This subroutine reads the treelist for the entire stand and subdivides it into grid rectangles (or pixels). The dimensions of each pixel are as small as possible such that the x and y dimensions are larger than the value of the PlotRadius variable, and will divide evenly into the respective x and y stand dimensions (MaxX and MaxY variables). Each pixel is assigned a vector within a 2-D array which contains a list of all trees that occur within that pixel. No changes are required in this subroutine. There is also a block of code identified in the Create_MGM_TreeList subroutine that will be required for each new sample algorithm. This block of code creates a subset of the main treelist that covers a mapped area slightly larger than each sample plot. The remaining subroutines in Module 5 are unique to specific surveys, and are called from the Run_Modeling_TreeLists subroutine as required. For MGM, there is a main subroutine called Create_MGM_Stands that controls the flow and output, and a second routine called Create_ MGM_TreeList that carries out instruction specific to individual sample plots. The functions of these routines include: 1. Count through the desired number of sample plots 2. Read in a subset of the main treelist specific to each plot 3. Create a 2-D array containing trees from the subset are in the sample plot 4. Write the 2_D array to a new file, with a column order suited to the specific model’s input requirements, alon with any other required formatting and header requirements. The data available to be used by the sampling subroutines include: 1. A 2-D array labeled PlotTrees containing all trees to be evaluated in the plot. Each row in the array is an individual tree, and each column is an attribute as described in Appendix A-1. The order of attributes is the same as can be found on the “TreeList” worksheet. 2. A set of variables including PlotRadius and any other variables or constants defined earlier that are specific to the plot type. 138 Silviculture Survey Simulator - User Guide Also available are a set of built in functions which can be found in Module 3. Of particular interest is a function that calculates the distance between two points: Dist(x 1 , y 1 ,x 2 ,y 2 ) 139 Appendix B A Silviculture Survey Methodology for Boreal Mixedwoods in Northeastern BC Appendix B - A Silviculture Survey Methodology for Boreal Mixedwoods in Northeastern BC 140 A Silviculture Survey Methodology For Boreal Mixedwoods In Northeastern BC Craig Farnden RPF Note: Minor formatting changes have been made to this document from a previous version released in October 2009 A Survey Silviculture Methodology for Boreal Mixedwoods in Northeastern BC 141 ACKNOWLEDGEMENTS This document has been prepared under contract to Canadian Forest Products Ltd. as part of a PhD project by the author in the Faculty of Forestry at the University of British Columbia. Gratefully acknowledged are the patience, guidance and advice of Dr. Bruce Larson, Dr. Peter Marshall and Dr. John Nelson, staff at the BC Ministry of Forests and Range - Forest Practices Branch, and staff of Canadian Forest Products Ltd, BC Timber Sales and the BC Ministry of Forests in the Peace River region of British Columbia. Reviews of this document have been provided by Patrick Smook (Canfor), Maureen Schulting (Pathocon Consulting) and Dave Weaver (MoF, Forest Practices Branch). Their contributions have greatly improved the outcome. However, any remaining errors or inconsistencies are solely the responsibility of the author. The layout of this document follows closely that of an earlier document produced by J.S. Thrower and Assoc. Ltd. that describes a similar survey procedure. Some portions of the surveys are identical, and in such cases text describing field procedures has been copied verbatim. Even where text has not been copied, it should be recognized that the current survey procedure builds upon the work of these and other authors. Their precedents are appreciated. A Survey Silviculture Methodology for Boreal Mixedwoods in Northeastern BC 142 Table of Contents 1.0 Introduction................................................................................................................ 5 2.0 Mixedwood Survey .................................................................................................... 5 2.1 Mixedwood Survey Objectives............................................................................. 5 2.2 Mixedwood Survey Overview............................................................................... 5 2.3 Target Population ................................................................................................. 6 2.4 Office Procedures ................................................................................................ 6 2.4.1 Map and Previous Data ............................................................................... 6 2.4.2 Office Stratification ...................................................................................... 6 2.4.3 Plot Locations .............................................................................................. 7 2.5 Field Sampling ..................................................................................................... 8 2.5.1 Quadrats versus Quadrants ........................................................................ 8 2.5.2 Stratification................................................................................................. 8 2.5.3 Quadrats ...................................................................................................... 9 2.5.4 Enhanced Plots .........................................................................................11 3.0 Survey Outcomes – Yield and Stand type................................................................11 3.1 Overview .............................................................................................................11 3.2 Objectives .......................................................................................................... 12 3.3 Model Development ........................................................................................... 12 3.3.1 MGM Simulations ...................................................................................... 13 3.3.2 Simulated Surveys..................................................................................... 13 3.3.4 Spruce Site Index ...................................................................................... 15 3.3.5 Brush Impacts............................................................................................ 15 3.3.6 Non-Productive Quadrats .......................................................................... 15 3.3.7 TMV’s......................................................................................................... 15 4.0 Survey Outcomes – Inventory Labels ..................................................................... 16 5.0 References .............................................................................................................. 16 A Survey Silviculture Methodology for Boreal Mixedwoods in Northeastern BC 143 1.0 Introduction The purpose of this document is to outline a survey and compilation procedure for the assessment of post-harvest spruce-aspen mixedwood stands where the spruce component consists primarily of planted trees. The procedure is based on extensive modeling of mixedwood stand development, and uses survey outcomes to predict future yields and the partitioning of that yield into conifer and broadleaved components. The procedure is initially targeted at satisfying the management criteria for intimate mixtures of spruce and aspen as identified in the Sustainable Forest Management Plan for the Fort St. John Pilot Project. However, it should be generally applicable to any post- harvest juvenile, even-aged mixture of white spruce and trembling aspen across the boreal forest of BC and nearby regions of other jurisdictions. 2.0 Mixedwood Survey 2.1 Mixedwood Survey Objectives The goals of the survey are to: 1. Describe stand characteristics in sufficient detail to estimate the Predicted Mean Volume (PMV) at age 80 along with the volume breakdown (composition) by species group (conifer versus broadleaved). 2. Contrast PMV to a threshold percentage of the theoretical maximum achievable volume – the target mean volume (TMV) for a stand of the same species composition. 3. Produce inventory labels 4. Identify areas that would potentially benefit from silvicultural treatments. 2.2 Mixedwood Survey Overview The key components of the mixedwood survey are: 1. The survey in its current form can be applied to even-aged mixtures of white spruce and trembling aspen, with aspen also being an analogue for cottonwood. Minor components of other species such as lodgepole pine and paper birch can be tolerated. 2. Stands are surveyed approximately 15 years after harvest to estimate the predicted mean volume (PMV) at age 80, and the percentage of that volume produced by conifer species (mainly spruce). 3. Sample plots for stocking (quadrats) are located at 25 m intervals on strip lines spaced 100 m apart, based on UTM coordinates (NAD 83). Quadrats are 10 m2 in area (1.78 m radius). 4. Quadrat plots are classified based on the presence or absence of trees by species class and (for conifers) free growing status relative to non-tree competing vegetation: a) Null - plot is outside the net area to be reforested (NAR) A Survey Silviculture Methodology for Boreal Mixedwoods in Northeastern BC 144 b) NP c) Unstocked, stockable d) Conifer, well growing e) Conifer, impeded f) Broadleaved tree g) Both conifer and broadleaved; conifer well growing h) Both conifer and broadleaved; conifer impeded 5. Enhanced plots are located on a 200 m grid (every eighth sample point on alternating strip lines). Enhanced plots are 100 m2 (5.64 m radius) in size and are sampled in addition to the quadrat on the same point. 6. Measurements in enhanced plots include: a) Site Series b) % Cover of the tree layer and composition of the understorey c) Height and age of the dominant and co-dominant trees d) Site index 2.3 Target Population The target population to sample in a given year is the NAR created from harvesting 15 years previously, where the intention was to reforest with intimate or micro-scale mixtures of spruce and aspen. While the intent is to always survey stands at age 15, the predictive models are not significantly affected by assessment ages that deviate by at least up to 5 years. 2.4 Office Procedures 2.4.1 Map and Previous Data A Silviculture Prescription (SP) map (or equivalent) should be used to develop the plot locations of the stand survey and should be updated following each survey. This map should show block boundaries, NP areas, non-commercial cover (NCC), wildlife tree patches (WTP’s), riparian management areas (RMA’s), permanent access structures and temporary roads. The surveyor should be familiar with the block history. 2.4.2 Office Stratification Prior to field sampling, the following information should be transferred to the survey map: 1. NAR boundary 2. Standards unit boundaries from the SP (or equivalent) map Standards units can be combined if they have the same site index and the same preferred and acceptable (p+a) species. A Survey Silviculture Methodology for Boreal Mixedwoods in Northeastern BC 145 2.4.3 Plot Locations Strip lines are located at 100 m intervals. Under most conditions, strip lines can be positioned on a UTM grid (NAD 83) at even 100 m intervals (last two digits of either the northing or easting are 00, depending on whether the strips are oriented E-W or N-S respectively). Quadrats are spaced at 25 m intervals on the strip lines, with enhanced plots at 200 m intervals on alternating strips (Fig B-1). In keeping with conventions for other surveys, enhanced plots should be placed where both the northing and easting coordinates are even multiples of 200 m (000, 200, 400, 600 or 800). For both quadrat and full measure plots, allowance should be made at the office stage for the possibility of grid points that appear to be outside the NAR according to mapped boundaries, but are found to be within the NAR in the field. In some cases, it will be inappropriate to run strip lines on a N-S or E-W bearing. The classic case is where linear features affecting stocking have been created through management actions. Such cases can lead to significant bias in survey outcomes, as the probability of any one plot falling within a particular stocking condition are affected by the pattern. There are two options in such cases, depending on the nature of the pattern: 1. If the striped pattern has irregular stripe widths or the stripes do not follow straight lines, strip lines can be run in either the N-S or E-W direction, depending which is closest to being perpendicular to the stripes. Plots are placed at the standard intervals. Figure B-1. Strips and Plots. Quadrat plots are located at 25 m intervals on strips lines spaced 100 m apart. Full measure plots are located at 200 m x 200 m grid points, corresponding to every eighth quadrat on alternating strip lines. A Survey Silviculture Methodology for Boreal Mixedwoods in Northeastern BC 146 2. If the striped pattern has regular stripe widths and stripes follow straight lines, strip lines can be run as in option 1 above, but plots must be spaced at irregular intervals. This can be achieved either by a. plotting each successive quadrat interval as a random distance between 4 and 46 m (a range from 0 to 50 would allow for overlapping plots which is less than desirable), or b. starting with a regular 25 m interval, and adding a random distance along the strip line ranging from –23 to +23 m. In both of these cases, enhanced plots should be placed not on regular grid points but at every eighth quadrat on alternating strips. 2.5 Field Sampling 2.5.1 Quadrats versus Quadrants An unfortunate source of potential confusion in this field procedure arises from the similar spelling and sound of two technical terms used in this versus other similar survey methods (e.g. Mean Stocked Quadrant or MSQ surveys): Quadrat: a fixed area plot, typically quite small, used as a measure of habitat in which individuals in a population are tallied to assess local distribution. The word quadrat originally referred to a rigid square frame used to define the limits of the plot, which in turn is derived from the squarish shape of the quadrate bone in the skulls of birds or reptiles. In practice quadrats may be any shape but are most commonly square or circular (Oxford University Press 2005). Quadrant: one fourth of a circular plot that is subdivided by two diameters passing through the plot center and oriented to each other at right angles. In practice, each quadrant of a circular plot is considered as roughly equivalent to a quadrat (as in MSQ surveys). One important difference is that quadrant data is usually summarized at the plot level first before calculating stand level statistics. In such cases, the number of samples, n, is based on the number of plots rather than the number of quadrants. If quadrant level statistics are treated as discrete samples, problems are encountered with hypothesis tests because the quadrants are spatially correlated and are not independent (independence is a key assumption of many statistical tests). 2.5.2 Stratification Prior to the survey, a walkthrough is required to complete the following: 1. Review block in the field and update the map. This may require identification and GPS location of unmapped natural features or new developments such as well sites and seismic lines. 2. Preliminary identification or confirmation of strata. Strata are areas greater than 2 ha with the same forest cover, site type, stocking class and treatment opportunity. A Survey Silviculture Methodology for Boreal Mixedwoods in Northeastern BC 147 By the completion of the survey, the cutblock must be divided into confirmed strata (polygons), based on common local and Provincial criteria and subject to the following requirements: 1. Inventory labels for the polygon must accurately specify species composition and site index, as per section 4.7.1 of the Surveys Manual1. 2. Areas with different TMV (if any) due to differences in moisture regime or other rationale must be delineated. 3. Areas with different stocking classes (NSR, SR and WG) must be delineated. 4. Areas with different mixedwood stand type classes (D, DC, CD and C) must be delineated 5. Areas with treatment opportunities must be delineated. If, during the survey, the surveyor identifies and stratifies an area 2 ha or larger that would benefit from further treatment (i.e. fill planting, control of non-crop vegetation), then the boundaries of the area must be transferred to the map. A description and rationale for the treatment must be provided. Note: post-stratification based on grouping like plots on a map is not an acceptable practice and will not be tolerated under any circumstances. 2.5.3 Quadrats Quadrats are milhectare (10 m2 or 1.78 m radius) plots that are classified based on the presence of certain combinations of trees, productive capability and competitive environment. The classes are: Null Plot is outside the NAR; plot is not included in any subsequent summaries Unstockable No acceptable trees are present; plot is non-productive ground Unstocked, stockable No acceptable trees are present but plot is potentially stockable. Comments are required regarding potential reasons for no acceptable trees being present and feasibility of adding trees. Conifer, well growing The only acceptable trees present are conifers. At least one conifer is of an acceptable species, an acceptable size and condition, and is free from overtopping non-crop vegetation (see well growing definition below). 1. The Surveys Manual is available at: http://www.for.gov.bc.ca/hfp/publications/00099/Surveys/Silviculture%20Survey%20Proce dures%20Manual-April%201%202009.pdf A Survey Silviculture Methodology for Boreal Mixedwoods in Northeastern BC 148 Conifer, impeded The only acceptable trees present are conifers. At least one conifer is of an acceptable species and is of acceptable size and condition, but all such trees are impeded by competition from non-crop vegetation (see well growing definition below). Broadleaved tree The only acceptable trees present are broadleaved. Mixed; conifer well growing Both acceptable conifers and acceptable broadleaved trees are present, and at least one conifer is free from overtopping tall shrubs (see free growing definition below). Mixed; conifer impeded Both acceptable conifers and acceptable broadleaved trees are present; all acceptable conifers are impeded by tall shrubs (see free growing definition below). In addition to classifying the plot, record: 1. pest codes for any pests present, along with the number of live and dead trees affected. 2. an ocular estimate of the total number of trees by species (Sx, Pl, Bl, Sb, Lt, At, Ac and Ep). Conifer Tree Acceptability – trees must be at least 30 cm tall to be tallied (unless otherwise defined in an SP or equivalent document), and must be free of insect or disease conditions that will prevent them from becoming merchantable-sized trees (minor growth limiting pests and deformities do not limit acceptability) - see FS 660. Broadleaved Tree Acceptability – trees must be at least 150 cm tall and a) the tree pith must not be laterally displaced more than 30 cm from the root crown pith, and b) the tree must not originate from a cut stump (aspen and cottonwood), and c) the tree must have at least one live leader, and d) the tree must not have a wound that is greater than 10% of the stem circumference or greater than 10% of the stem length, and e) the tree must not have any fungal or insect infection affecting tissues below the bark surface (visible without destructive sampling), and f) the tree must not be browsed so as to limit its ability to become a crop tree. Damage to trees where the tree is likely to overcome the damage and continue growing as a dominant or co-dominant tree should not limit the tree from being tallied. A Survey Silviculture Methodology for Boreal Mixedwoods in Northeastern BC 149 Conifer well growing assessments –conifer trees are well growing if: a) they meet the acceptability criteria above, and b) they are at least 100% of the height of all herbaceous and grass competition within a 1 m radius cylinder, and c) they are taller than shrub competition in three of four quadrants of a 1 m radius cylinder It is important to recognize that broadleaved crop trees (aspen, cottonwood and birch) are ignored in well growing assessments for conifers – the competitive effects of these species are accounted for in a different fashion. 2.5.4 Enhanced Plots In addition to the quadrat, enhanced sample points include a 100 m2 (5.64 m radius) plot to collect BEC , site index and additional inventory label data: ß assess and record site series, complete with zone, subzone and variant ß estimate percent cover of the tree layer ß record percent cover and average height of non-crop vegetation by species; include any species that might compete with crop trees ß record mean height and age of dominant and co-dominant crop trees by species ß for each of the dominant conifer and broadleaved crop tree species in the stratum for which growth intercept equations are available, locate the fattest tree in the plot for potential use as a growth intercept sample tree. If a selected tree meets the acceptability criteria listed below, record the height and breast height age. If the fattest tree does not meet the acceptability criteria, use the second fattest tree. If the second fattest tree is also unsuitable, attempt to find a suitable site tree in a 100 m2 plot centered on each successive quadrat until a suitable replacement is found. To be suitable for growth intercept evaluations of site index, a tree must: ß have three or more years height growth above breast height, and ß have an undamaged stem with vigorous, relatively uniform height growth above breast height, and ß have vigorous and uniform ring width at breast height, and ß must not be directly overtopped by other trees or brush, and ß must be in a dominant or co-dominant position 3.0 Survey Outcomes – Yield and Stand type 3.1 Overview This free growing assessment procedure for intimate and micro-scale boreal mixedwoods is based on an analysis framework described by Farnden (2009a). This framework uses a combination of simulated surveys and projections of yield from A Survey Silviculture Methodology for Boreal Mixedwoods in Northeastern BC 150 existing growth models to build new simplified models using input parameters derived from silviculture survey data. Model outcomes include Predicted Mean Volume (PMV), Target Mean Volume (TMV) and Species Composition (Conifer Percentage of PMV). PMV is the predicted total yield at age 80, and TMV is a threshold percentage of the maximum potential yield for a stand of similar species composition. 3.2 Objectives There are two primary objectives for predicting yields for mixedwood stands using silviculture survey data as predictive parameters: 1. Assessing compliance with statutory or contractual reforestation obligations 2. Assessing stand level contributions to landscape level management objectives If the TMV exceeds the PMV for a stand or an area weighted collection of stands, the reforestation obligations are deemed to have been met. Compliance with landscape level species composition objectives are tracked by classifying each stand into a Mixedwood Stand Type class2, and tracking the area reforested into each class in a ledger system (see Martin et al. 2002). 3.3 Model Development Details on the generation and functioning of the models for calculating TMV and projected species composition are provided here for reference. While an under- standing of the system is useful, in most cases field staff will rely on spreadsheets or custom programs to compile survey results and calculate the required outcomes. Separate prediction models were developed for each 1 m increment of white spruce site index, covering a productivity range from SI50 = 13 to SI50 = 27. For each model, a set of 113 mixedwood stands (9 ha each) was generated using the following criteria: ß 9 levels of area coverage by aspen clones from 10 to 90% in 10% increments, where the aspen density varied randomly fro 10,000 to 16,000 trees/ha ß 4 levels of aspen density between clonal patches at 0, 50, 200 and 800 trees/ha ß spruce planted at 1400 trees/ha and three levels of survival: 30, 60 and 90% ß aspen site index was calculated from the spruce site index using a conversion equation provided by Nigh(2002): SIAt = -4.768 + 1.253 x SISw For each stand, yields were predicted using the Mixedwood Growth Model3 (MGM), and surveys were simulated using the Silviculture Survey Simulator developed as part of the PhD work of Farnden (2009b,c) at the University of BC. 2. Stand Type Classes: D = 0 to 20% conifer volume, DC = 20 to 50% conifer volume, CD = 50 to 80% conifer volume and C = 80 to 100% conifer volume. 3. Mixedwood Growth Model (MGM): an individual tree, distance independent growth model developed at the University of Alberta for use in boreal mixedwoods. Information on the model can be accessed at http://www.ales2.ualberta.ca/rr/mgm/mgm.htm (as accessed April 2, 2009). A Survey Silviculture Methodology for Boreal Mixedwoods in Northeastern BC 151 3.3.1 MGM Simulations For each simulated 9 ha stand, a set of 36 sampled treelists was generated using randomly located 100 m2 plots and used to generate a discrete yield curve. The 36 yield curves were then combined to generate a composite yield curve for the stand. Each sample treelist was assumed to represent local growing conditions, such that variation in competitive conditions throughout the stand would be appropriately represented in the predicted yields. Stands were grown from the initial age of 12 years to a projected age of 80 years using BC site curves and merchantability limits set at a 12.5 cm dbh limit, 30 cm stump height and 10 cm minimum top diameter. No yield deductions were made for waste or decay. 3.3.2 Simulated Surveys Surveys were simulated using survey algorithms built into the Silviculture Survey Simulator. For each stand, survey statistics were compiled from a set of 144 randomly located 10 m2 quadrats, each of which was classified as: 1. Unstocked: no trees present 2. Aspen: at least 1 broadleaved tree present; no conifers 3. Spruce: at least one conifer tree present; no broadleaves 4. Mixed: at least one broadleaved tree and at least one conifer present Survey summary statistics are simply the percentages of plots in each class (C0, C1, C2 and C3 respectively). 3.3.3 Model Fitting Models for PMV and species composition (% spruce) are based on either on multiple linear regression using ordinary least squares (OLS), a complex variant of a Weibull function or a complex variant of a logarithmic function (Farnden 2009c): for SI = 13 to 21 for SI = 22 to 27 for SI = 13 to 18 for SI = 19 to 22 for SI = 23 to 27 Where b0 through b6 are separate sets of coefficients for each model. C3 is not utilized as it is redundant information. All models were fit using JMP version 8.0.1. Model parameters and fit statistics are provided in Tables B-1 and B-2. % Sw PMV = b0 + b1C0 + b2C1 + b3C2 = b6 * 1 e b4 (b0 +b1C0 +b2C1+b3C2 ) b5 =1 e b4 (b0 +b1C0 +b2C1+b3C2 ) b5 = b4 + b5 log b0 + b1C0 + b2C1 + b3C2( ) = b0 + b1C0 + b2C1 + b3C2 A Survey Silviculture Methodology for Boreal Mixedwoods in Northeastern BC 152 Table B-1. Coefficients and fit statistics for the models used to calculate PMV assuming discrete site index values. Separate models have been generated for different levels of site productivity based on 1 m increments of white spruce site index. Table B-2. Coefficients and fit statistics for the models used to calculate yield-based species composition (% Spruce) assuming discrete site index values. Separate models have been generated for different levels of site productivity based on 1 m increments of white spruce site index. b0 b1 b2 b3 b4 b5 b6 R 2/I2 RMSE 13 89.38 -74.63 -68.79 9.52 0.961 3.8 14 109.47 -90.22 -72.35 12.52 0.957 4.45 15 135.28 -110.34 -77.09 12.77 0.953 5.28 16 168.83 -136.33 -84.62 8.94* 0.946 6.49 17 212.86 -170.26 -94.54 -0.48* 0.9348 8.31 18 270.41 -214.11 -113.56 -17.5 0.928 10.51 19 338.82 -265.93 -136.47 -40.43 0.925 13.04 20 415.78 -323.6 -161.1 -68.23 0.9224 15.85 21 496.45 -382.77 -178.63 -97.47 0.9187 18.84 22 36.53 -15.413 -9.318 -7.419 1.772E-07 4.734 533.84 0.913 21.77 23 26.61 -10.244 -5.861 -5.0277 1.076E-07 5.372 589.8 0.909 23.58 24 8.544 -0.1417 -0.0734 -0.07096 1.336E-151 162.8 647.5 0.905 28.43 25 5.75 -1.031 -0.497 -0.5597 1.357E-10 14.08 702.6 0.892 27.24 26 2.829 -0.3355 -0.1802 -0.1834 2.516E-11 25.74 742.4 0.877 29.5 27 1.518 -0.001953 -0.001184 -0.0011 1.1953E-106 1741.7 787.56 0.875 30.47 SI Sw 50 Y80 b0 b1 b2 b3 b4 b5 R 2/I2 RMSE 13 1.124 -0.6415 -1.1225 -0.2422 8.868 1.814 0.955 0.0349 14 0.8877 -0.5186 -0.9103 0.0938 6.024 1.36 0.945 0.0421 15 0.9351 -0.4746 -0.9621 0.4181 3.567 1.194 0.931 0.0498 16 1.115 -0.4773 -1.1511 0.8132 2.059 1.087 0.914 0.0568 17 0.6057 -0.1881 -0.6284 0.6201 2.829 1.027 0.898 0.0628 18 0.5691 -0.0594 -0.5955 0.7979 2.101 1.04 0.888 0.065 19 0.6052 0.0906 -0.4225 0.4346 0.9468 0.6249 0.888 0.0647 20 0.5753 0.118 -0.3072 0.3447 1.02 0.8275 0.887 0.065 21 0.5567 0.0724 -0.1515 0.1588 1.514 1.702 0.887 0.0645 22 0.7754 0.0124 -0.0264 0.0241 4.207 14.57 0.898 0.0618 23 0.5011 0.2052 -0.5115 0.4461 0.908 0.058 24 0.5147 0.1647 -0.54 0.4301 0.9181 0.055 25 0.5399 0.1074 -0.5778 0.4021 0.9271 0.052 26 0.5757 0.0587 -0.6158 0.3618 0.9324 0.05 27 0.6273* 0.0268 -0.6551 0.3016 0.932 0.0496 SI Sw 50 Yc A Survey Silviculture Methodology for Boreal Mixedwoods in Northeastern BC 153 3.3.4 Spruce Site Index Spruce site index is required for model selection. Several sources of site index are available, listed in order of preference: 1. growth intercept estimate for spruce using data from surveyed block 2. 2nd or later generation SIBEC estimate4 3. growth intercept estimate for spruce using data from nearby block on same site series and soil type 4. 1st generation SIBEC estimate 5. Pre-harvest inventory estimate 3.3.5 Brush Impacts The growth model used to predict yields (MGM) is unable to explicitly model the effects of non-tree brush competition. The effects of these competitors must therefore be added in after the generation of the fitted models. Adjustments are based on the assumption that tall shrubs overtopping spruce will have a temporary effect on spruce growth, thus reducing spruce yields. This is handled within the survey system by: a) Subdividing the classes of plots where conifers are present (C2 and C3) into those where at least one conifer is free of impeding competition from tall shrubs (mainly alder and willow) versus those where all conifers are impeded b) Assuming that the impeding vegetation will, on average, result in a 10% reduction of yield for the affected trees For the impeded spruce plots, 10% are then transferred to the unstocked (C0) category, and the remainder to the unimpeded spruce (C2) category For the impeded mixed plots, 10% are transferred to the aspen (C1) category, and the remainder to the mixed (C3) category. 3.3.6 Non-Productive Quadrats Non-productive (NP) and null (plot center is outside of NAR) quadrats are treated as exclusive to the sampled population and are not included in the calculation of the percentage of plots in each stocking category. 3.3.7 TMV’s TMV is not fit as a model on its own, but is calculated using the PMV model. For TMV, the number of unstocked plots in excess of 5% are redistributed proportionately to the stocked classes (where 95% stocked plots is a theoretical upper ceiling) to determine a theoretical maximum yield. The threshold volume for acceptable regeneration performance is then set as some proportion of the theoretical maximum yield. For 4. Second generation and later estimates are recognized in the SIBEC tables by the inclusion of a standard error value. A Survey Silviculture Methodology for Boreal Mixedwoods in Northeastern BC 154 example, the value set in the 2003 Sustainable Forest Management Plan for the Fort St. John Pilot Project is 0.95. Example: A stand has survey outcomes with C0=8.1%, C1=21.4%, C2=15.4% and C3=55.1%. In this case, there is 3.1% excess in the unstocked (C0) class that must be redistributed to the other classes: C0(new) = 0.05 C1(new) = C1 + C1*0.031 ÷ (1-0.081) = 0.221 C2(new) = C2 + C2*0.031 ÷ (1-0.081) = 0.159 C3(new) = C3 + C3*0.031 ÷ (1-0.081) = 0.570 To determine TMV using the 0.95 threshold from the Fort St. John SFMP, the calculation becomes 0.95 x PMV(0.05, 0.221, 0.159) 4.0 Survey Outcomes – Inventory Labels Elements of the inventory label for reporting to the BC Ministry of Forests’ Results database are determined as follows: Species Composition: based on the total tally of trees by species in all quadrats: Height and Age: based on the estimates of the mean height and age of the dominant and co-dominant trees of the 1st and 2nd leading species from the enhanced plots. Site Index: based on growth intercept estimates of the dominant species where available, otherwise using the best available alternate method. Crown Closure: based on mean crown closure as estimated for each of the enhanced plots Total Trees: based on the tally of trees within the quadrats 5.0 References Coates, K.D. and Burton, P.J. 1999. Growth of planted tree seedlings in response to ambient light levels in northwestern interior cedar-hemlock forests of British Columbia. Can. J. For. Res. 29:1374-1382. Farnden, C. 2009a. An analysis framework for linking regeneration standards to desired future forest conditions. Forestry Chronicle 85(2): 285-292. Species % i = Speciesi All �Trees A Survey Silviculture Methodology for Boreal Mixedwoods in Northeastern BC 155 Farnden, C. 2009b. Silviculture survey simulator user guide. Contract report to Canadian Forest Products Ltd., Fort St. John, BC. 34p. Farnden, C. 2009c. Development of SFM-based regeneration standards for spruce- aspen mixedwoods in northeastern BC. PhD Thesis, university of British Columbia. In Prep. J.S. Thrower & Assoc., 2003. Stand survey & growth modeling for the Fort St. John TSA. Contract report to Canadian Forest Products Ltd., Chetwynd BC. 32p. Martin, P.J., Browne-Clayton, S., and McWilliams, E. 2002. A results-based system for regulating reforestation obligations. For. Chron. 78: 492-498. Nigh, G. 2002. Site index conversion equations for mixed trembling aspen and white spruce stands in northern British Columbia. Silva Fennica 36(4): 789–797. Oxford University Press, 2005. New Oxford American dictionary. 2nd Ed., embedded in Mac OS X v. 10.5.