Growth of British Columbia Coastal Species in Response to Thinning and Fertilization Treatments by LEAH C. RATHBUN Bachelor of Geological-Engineering, University of Minnesota, 2004 B.Sc., Forest Resources Management, University of Minnesota, 2004 M.Sc., University of Minnesota, 2007 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) November 2010 © LEAH C. RATHBUN, 2010 ii Abstract The successional processes of the mixed-species Pacific coastal temperate rain forests of British Columbia (BC), Canada, are defined by gap dynamics, where small- scale disturbances, mainly due to windthrow, create openings in the canopy necessary for regeneration. Douglas-fir (Pseudotsuga menziesii var. menziesii (Mirb.) Franco) is the dominant, pioneer species in this area and western hemlock (Tsuga heterophylla (Raf.) Sarg.) and western redcedar (Thuja plicata Donn) are the late-successional, shade-tolerant species. Silvicultural systems such as variable retention systems have been applied to many of the secondary growth mixed-conifer forests. Variable retention in this area is designed to differ dramatically from stand to stand. This approach differs from the traditional even-aged management applied to forests of the Pacific Northwest coast. In this study, a model-based approach was used to investigate how multiple treatment interventions as a part of active management across a landscape affect mortality and growth within actively managed stands. There is a need for this information as current growth and yield models used in this area are limited by either the number of species which can co-exist in a stand (e.g., the model TASS of BC) or are limited by the need for data not commonly obtained in inventory databases (e.g., the models FVS and ORGANON of USA). Additionally, no growth and yield models have been developed to include variable retention systems, where a variety of thinning intensities and spatial patterns, timing of thinning and fertilization treatments, and number of treatments are used. Mortality, diameter increment, and height increment models were developed and the effects of thinning, fertilization and the combination of thinning and fertilization were iii examined for Douglas-fir, western redcedar, and western hemlock. For each species, shade-tolerance was found to impact the possible predictor variables included in model development. The use of a generalized logistic survival model resulted in accurate estimates for larger trees, but poor results for smaller trees. To model the effects of fertilization, additional fertilization effect variables were included in the models; conversely, thinning effects were modeled using the immediate change in state variables such as basal area of larger trees which occurred immediately following thinning. iv Table of contents Abstract ............................................................................................................................... ii Table of contents ................................................................................................................ iv List of tables ..................................................................................................................... viii List of figures ...................................................................................................................... x Acknowledgements .......................................................................................................... xiii Dedication ........................................................................................................................ xiv Co-authorship statement ................................................................................................... xv 1. Introduction ................................................................................................................. 1 1.1 Background .......................................................................................................... 1 1.2 Literature review .................................................................................................. 2 1.2.1 Forest practices ............................................................................................. 2 1.2.2 Growth and yield models .............................................................................. 3 1.2.3 Mortality ....................................................................................................... 6 1.2.4 Growth .......................................................................................................... 9 1.3 Gaps in knowledge and practice......................................................................... 10 1.4 Research goals .................................................................................................... 12 1.5 Thesis structure .................................................................................................. 16 1.6 References .......................................................................................................... 19 2. Modeling mortality for mixed-species stands located in coastal British Columbia . 26 2.1 Introduction ........................................................................................................ 26 2.2 Materials and methods ....................................................................................... 29 2.2.1 Study area.................................................................................................... 29 v 2.2.2 Data ............................................................................................................. 30 2.2.3 Model selection ........................................................................................... 31 2.2.4 Model implementation ................................................................................ 35 2.3 Results ................................................................................................................ 38 2.3.1 Model selection ........................................................................................... 38 2.3.2 Model implementation ................................................................................ 39 2.4 Discussion .......................................................................................................... 41 2.5 Conclusion .......................................................................................................... 45 2.6 References .......................................................................................................... 55 3. Modeling the effects of thinning and fertilization on mortality of coastal Douglas-fir (Pseudostuga menziesii var. menziesii. (Mirb.) Franco) ................................................... 59 3.1 Introduction ........................................................................................................ 59 3.2 Materials and methods ....................................................................................... 63 3.2.1 Study area.................................................................................................... 63 3.2.2 Data ............................................................................................................. 63 3.2.3 Survival model ............................................................................................ 65 3.2.4 Treatment effects ........................................................................................ 67 3.2.5 Model accuracy ........................................................................................... 70 3.3 Results ................................................................................................................ 71 3.3.1 Treatment effects ........................................................................................ 71 3.3.2 Model accuracy ........................................................................................... 72 3.4 Discussion .......................................................................................................... 73 3.5 Conclusion .......................................................................................................... 77 vi 3.6 References .......................................................................................................... 89 4. Diameter growth models for mixed-species stands of Coastal British Columbia including thinning and fertilization effects ....................................................................... 94 4.1 Introduction ........................................................................................................ 94 4.2 Materials and methods ....................................................................................... 97 4.2.1 Study area.................................................................................................... 97 4.2.2 Data ............................................................................................................. 97 4.2.3 Box-Lucas model form ............................................................................... 99 4.2.4 Model development for untreated plots .................................................... 101 4.2.5 Empirical model for comparison .............................................................. 105 4.2.6 Treatment effects ...................................................................................... 106 4.2.7 Model accuracy ......................................................................................... 108 4.3 Results .............................................................................................................. 109 4.3.1 Model development for untreated plots .................................................... 109 4.3.2 Empirical model for comparison .............................................................. 110 4.3.3 Treatment effects ...................................................................................... 111 4.3.4 Model accuracy ......................................................................................... 112 4.4 Discussion ........................................................................................................ 114 4.5 Conclusion ........................................................................................................ 119 4.6 References ........................................................................................................ 134 5. Height growth models: Potential, average, and responses to silvicultural treatments ……………………………………………………………………………………..140 5.1 Introduction ...................................................................................................... 140 vii 5.2 Materials and methods ..................................................................................... 143 5.2.1 Study area.................................................................................................. 143 5.2.2 Data ........................................................................................................... 145 5.2.3 Box-Lucas model ...................................................................................... 146 5.2.4 Potential height growth model .................................................................. 147 5.2.5 Height growth model ................................................................................ 151 5.2.6 Treatment effects ...................................................................................... 154 5.2.7 Model accuracy ......................................................................................... 157 5.3 Results .............................................................................................................. 158 5.3.1 Potential height growth model .................................................................. 158 5.3.2 Height growth model ................................................................................ 159 5.3.3 Treatment effects ...................................................................................... 162 5.4 Discussion ........................................................................................................ 165 5.5 Conclusion ........................................................................................................ 168 5.6 References ........................................................................................................ 188 6. Conclusions ............................................................................................................. 194 6.1 Research conclusions ....................................................................................... 194 6.1.1 Mortality ................................................................................................... 194 6.1.2 Growth ...................................................................................................... 195 6.1.3 Experimental design vs. managed stands.................................................. 196 6.2 Management implications ................................................................................ 196 6.3 Future research ................................................................................................. 197 6.4 References ........................................................................................................ 200 viii List of tables Table 2-1. Plot summary statistics at plot establishment and the end of the first measurement period. ......................................................................................................... 50 Table 2-2. Tree summary statistics at plot establishment by species. ............................. 51 Table 2-3. AIC values for selection models averaged across three model subsets. ......... 52 Table 2-4. AIC values averaged over three model datasets and concordance averaged over three test datasets by model and species. .................................................................. 53 Table 2-5. Bias and root mean-squared error for live stems per ha using the three implementation approaches for each species. ................................................................... 54 Table 3-1. Plot summary statistics at plot establishment and at time of treatment .......... 86 Table 3-2. Means and standard deviations for tree-level variables at plot establishment by silvicultural treatment. ................................................................................................. 87 Table 3-3. Parameter estimates for the selected survival model for untreated plots. ...... 88 Table 4-1. Plot summary statistics at plot establishment and at time of treatment. ....... 127 Table 4-2. Means and standard deviations for tree-level variables at plot establishment by silvicultural treatment and species. ............................................................................ 128 Table 4-3. Correlations between possible predictor variables and estimated 1 or 2 for trees with more than three measurements. ...................................................................... 129 Table 4-4. AIC values for augmented Box-Lucas model with varying error structures for trees on untreated plots. .................................................................................................. 130 Table 4-5. Parameter estimates for the augmented Box-Lucas model for untreated plots ......................................................................................................................................... 131 ix Table 4-6. AIC, mean bias, and root mean-squared error values by species for the augmented Box-Lucas model vs. The Empirical model for untreated plots. .................. 132 Table 4-7. AIC values for augmented Box-Lucas model modified to include the fertilizer variable for fertilized plots. ............................................................................................. 133 Table 5-1. Plot summary statistics at plot establishment and at time of treatment. ....... 183 Table 5-2. Means and standard deviations for tree-level variables at plot establishment by silvicultural treatment and species. ............................................................................ 184 Table 5-3. Correlations between possible estimateor variables and predicted 1 or 2 for trees with more than three measurements. ................................................................ 185 Table 5-4. AIC values for augmented Box-Lucas model with varying error structures for trees on untreated plots. .................................................................................................. 186 Table 5-5. AIC values for augmented Box-Lucas model modified to include the fertilizer variable for fertilized plots. ............................................................................................. 187 x List of figures Figure 2-1. Plot of sensitivity and specificity by threshold value for the selected model and the two validation models for each species by small and large diameter sizes. ......... 47 Figure 2-2. Plot of bias and root mean-squared error by DBH class for the final model using all implementation approaches for each species. .................................................... 48 Figure 2-3. Plot of bias and root mean-squared error by incidence of mortality class for the final model using all implementation approaches for each species. ........................... 49 Figure 3-1. Average actual and predicted annual percent mortality for thinned data by time since thinning. ........................................................................................................... 78 Figure 3-2. Average predicted annual percent mortality by thinning intensity for thinned data by time since thinning. .............................................................................................. 79 Figure 3-3. Average predicted annual percent mortality by stand development for thinned data by time since thinning. ................................................................................. 80 Figure 3-4. Average predicted annual percent mortality by number of thinnings for thinned data by time since thinning. ................................................................................. 81 Figure 3-5. Average actual and predicted annual percent mortality for fertilized data by time since fertilization....................................................................................................... 82 Figure 3-6. Average predicted annual percent mortality by site index class for fertilized data by time since fertilization .......................................................................................... 83 Figure 3-7. Average actual and predicted annual percent mortality fertilized and thinned in combination by time since treatment ............................................................................ 84 Figure 3-8. Bias and RMSE values by 5 cm DBH classes by treatment. ....................... 85 Figure 4-1. Box-Lucas model for varying 1 and 2 values for diameter increment .. 120 xi Figure 4-2. Mean bias values by 5 cm DBH classes for diameter increment ............... 121 Figure 4-3. RMSE values by 5 cm DBH classes for diameter increment ..................... 122 Figure 4-4. Average actual and predicted annual diameter increments for the untreated data by 5 cm DBH classes ............................................................................................. 123 Figure 4-5. Average actual and predicted annual diameter increments for the fertilized plots by time since fertilization. ...................................................................................... 124 Figure 4-6. Average actual and predicted diameter increments for the thinned plots by time since thinning. ......................................................................................................... 125 Figure 4-7. Average actual and predicted annual diameter increment for the multiple treatment data .................................................................................................................. 126 Figure 5-1. Box-Lucas model for varying 1 and 2 values for height increment. ..... 170 Figure 5-2. Average annual potential height by 1m height classes using DBH and height as estimator variables. ..................................................................................................... 171 Figure 5-3. Average annual potential height by 1m height classes using DBH, site index, and height as estimator variables. ................................................................................... 172 Figure 5-4. Average annual average height by 1m height classes using the growth potential independent and dependent methods ............................................................... 173 Figure 5-5. Mean Bias values by 5 m height classes by treatment for height increment 174 Figure 5-6. RMSE values by 5 m height classes by treatment for height increment. ..... 175 Figure 5-7. Average annual mean untreated, fertilized, and potential height growth by 2m height classes for Douglas-fir ......................................................................................... 176 Figure 5-8. Average annual mean untreated, fertilized, and potential height growth by 2m height classes for western hemlock. ................................................................................ 177 xii Figure 5-9. Average annual mean untreated, thinned, and potential height growth by 2m height classes for Douglas-fir ......................................................................................... 178 Figure 5-10. Average annual mean untreated, thinned, and potential height growth by 2m height classes for western redcedar. ................................................................................ 179 Figure 5-11. Average annual mean untreated, thinned, and potential height growth by 1m height classes for western hemlock. ................................................................................ 180 Figure 5-12. Average annual mean untreated, fertilized & thinned, and potential height growth by 2m height classes for Douglas-fir .................................................................. 181 Figure 5-13. Average annual mean untreated, fertilized & thinned, and potential height growth by 1m height classes for western hemlock ......................................................... 182 xiii Acknowledgements I would like to begin by thanking: Dr. Valerie LeMay without whom this research would not have been conducted, her insightful observations regarding both research and life have guided me through this journey; Dr. Nick Smith whose critical thinking and vast knowledge of forestry assured my success with this project; my committee members, Drs. Peter Marshall and Lori Daniels who listened and encouraged me throughout this process; Dr. Antel Kozak who was available whenever I had questions, his knowledge of taper equations was invaluable; Dr. David Hann who kindly allowed me to visit and generously shared his research with me; Dr. Cindy Prescott who while supporting all graduate students within the faculty, also provided a roof over my head during the most critical part of this thesis; and Drs. Tom Burk and Alan Ek who helped build the foundation needed for this venture. Funding for this research was provided by the Forest Science Program (FSP) of British Columbia. The data and additional funding critical to this project was provided by Ken Epps and Island Timberland, LP. Thanks go out to my dear friends in the faculty: Andrew, Derik, Mariano, Craig, AnaElia, Babita, Sam, and Suzi; and to my dear friends in life: Wendy, the Kims, Emily, Liz, and all the lovely Germans. For those who watched over my pup when I needed to travel I am eternally grateful: Andrew, Jen, Janelle, Liam, Heike, and Ernie. xiv Dedication To my family who have supported me spiritually, financially, and emotionally throughout all my adventures. xv Co-authorship statement The tables below indicate the percentages of co-author’s contributions at each stage of research for each chapter. Chapter 2. Modeling mortality for mixed-species stands located in coastal British Columbia. Leah Rathbun Valerie LeMay Nick Smith Peter Marshall Lori Daniels Total Problem identification & research design 70 10 10 5 5 100 Performing the research 85 10 5 0 0 100 Data analyses 90 10 0 0 0 100 Manuscript preparation 70 15 5 5 5 100 Chapter 3. Modeling the effects of thinning and fertilization on mortality of coastal Douglas-fir (Pseudostuga menziesii var. menziesii (Mirb.) Franco). Leah Rathbun Valerie LeMay Nick Smith Peter Marshall Lori Daniels Total Problem identification & research design 80 5 5 5 5 100 Performing the research 95 5 0 0 0 100 Data analyses 95 5 0 0 0 100 Manuscript preparation 95 5 0 0 0 100 2 Chapter 4. Diameter growth models for mixed-species stands of coastal British Columbia including thinning and fertilization effects. Leah Rathbun Valerie LeMay Nick Smith Peter Marshall Lori Daniels Total Problem identification & research design 75 5 10 5 5 100 Performing the research 85 10 5 0 0 100 Data analyses 90 10 0 0 0 100 Manuscript preparation 70 15 5 5 5 100 Chapter 5. Height growth models: Potential, average, and responses to silvicultural treatments. Leah Rathbun Valerie LeMay Nick Smith Peter Marshall Lori Daniels Total Problem identification & research design 75 5 10 5 5 100 Performing the research 85 10 5 0 0 100 Data analyses 90 10 0 0 0 100 Manuscript preparation 70 15 5 5 5 100 1 1. Introduction 1.1 Background Forest systems are complex and management of these systems is continually changing. Management and silvicultural decisions are made not only on current stand conditions but also on anticipated future stand conditions. The ability of growth and yield models to update inventories and predict future yield make them a key component for exploring management and silvicultural alternatives and ultimately for making management decisions (Vanclay 1994). In the coastal temperate rain forests of British Columbia (BC), Canada and of the United States, variable retention of live trees in harvest units is a relatively recent forest management practice (Bunnell and Dunsworth 2009, Peterson and Anderson 2009). As a result, few studies exist on the impacts on growth and other forest variables (Aubry et al. 2009). Also, since the level of retention and the spatial distributions of residual trees vary greatly, it is difficult to characterize the impacts of retention on residual trees and on regeneration following harvest. As a result, growth and yield models are needed to simulate alternative variable retention scenarios and evaluate outcomes. The goal of this study was to develop a distance-independent individual tree model for second-growth managed coastal Douglas-fir (Pseudotsuga menziesii var. menziesii (Mirb.) Franco), western hemlock (Tsuga heterophylla (Raf.) Sarg.), and western redcedar (Thuja plicata Donn) of coastal BC. The models developed in this research and described in this dissertation were created for use in managed stands where variable retention systems may occur. In conjunction with a regeneration model (Smith 2010), these models were implemented within FORest Growth Engine (FORGE), a 2 model designed to examine the effects of variable retention systems on the growth of trees on Vancouver Island (Smith 2005). FORGE allows for the modification of height growth through changes in light availability. The final product, an individual tree distance-independent model is now available for industry (Island Timberlands, LP) and will be used to model short-term as well as long term timber supply. FORGE will assist forest managers by providing information needed for silvicultural investment decisions, harvest planning, and mill planning (i.e. by providing tree and log profiles). 1.2 Literature review 1.2.1 Forest practices The mixed-species coastal temperate rain forest of BC is noted for its lack of recent fires (Lertzman et al. 1996). The successional processes of the Pacific Northwest forests are defined by small-scale disturbances, mainly due to windthrow (Gavin et al. 2003). Hence, gap dynamics define the forest structure and natural successional processes on Vancouver Island (Lertzman et al. 1996), with Douglas-fir being the dominant, pioneer species and western hemlock and western redcedar being the late- successional species aggregated within gaps (Getzin et al. 2006). While mortality in younger stands is typically due to self-thinning processes, as a stand matures, gaps formed from single-tree mortality drive the structural and compositional characteristics instrumental in stand development (Lertzman et al. 1996), creating complex mixed- species stands. To mimic the natural disturbance regime in this area, silvicultural systems such as variable retention systems have been applied to many of the secondary growth mixed-conifer forests. Variable retention systems are implemented by applying varying intensities and sizes of disturbances across a gradient resulting in complex stand 3 structures and broad diversity of stands, thereby mimicking the natural successional process of the area (Mitchell and Beese 2002). This approach differs from the traditional even-aged management applied on the Pacific Northwest coast. 1.2.2 Growth and yield models There are typically three types of growth and yield models used to estimate future stand productivity: whole stand models, size class models, and individual tree models. Whole stand models utilize population parameters such as stand basal area and/or stand density to forecast stand growth or yield. These models are generally the least complex and provide information at the stand level only. Whole stand models can be adapted to incorporate stand structure such as size class distribution, termed size class models. These models do not include individual tree information, but rather aggregate trees within categories such as size and/or species classifications. One example of a size class distribution model created for uneven-aged stands is a stand table projection (e.g., model by Ek 1974). Individual tree-level models are the most detailed, where measures for each tree are projected forward in time. This tree-level information is then amassed to obtain stand-level information. Individual tree-level models can be classified as either distance- independent (aspatial) or distance-dependent (spatial information provided for each tree’s location). Gap models are a special case of distance-independent tree-level models, which use multiple regeneration patches to study natural stand dynamics, including long- term changes in species composition, rather than growth and yield (Bartelink 2000). Growth models can also include process components that drive growth and mortality. Process-based models use information regarding physiological processes such 4 as respiration, photosynthesis, and nutrient uptake to alter growth and mortality (e.g., 3PG, Landsberg and Waring 1997; Sands and Landsberg 2002). Recently, hybrid models have been developed by combining process-based models and either stand-level or tree- level empirical growth models. These hybrid models provide the ability to examine possible impacts of changes in physiological processes while using conventional tree allometry measures as model inputs. For example, TREE-BGC, a derivative of FOREST-BGC, calculates the cycling of carbon, water, and nitrogen within trees and has been applied to a small sample of uneven-aged interior Douglas-fir (Pseudotsuga menziesii var. glauca (Mirb.) Franco) stands near Kamloops, BC (Korol et al. 1996). Few hybrid models exist within the Pacific Northwest. Management processes such as thinning and fertilization have been incorporated in growth and yield models using categorical variables or continuous variables such as thinning ratios (e.g., d-ratio, defined as the ratio of the quadratic mean diameter after thinning to the quadratic mean diameter before thinning), residual basal area per hectare, residual trees per hectare, and/or individual tree characteristics. Currently, in the Pacific Northwest, growth and yield models incorporating thinning and/or fertilization practices exist. Douglas-fir Simulator (DFSIM) is an example of a whole stand model which incorporates different thinning regimes and was created for naturally regenerated even- aged Douglas-fir stands located in coastal Oregon and Washington, USA (Curtis et al. 1981). One notable problem with DFSIM is the lack of an ingrowth component for the model (Ritchie 1999). Individual-tree, distant-independent models also exist for the Pacific Northwest of the USA. The Forest Vegetation Simulator (FVS), otherwise named Prognosis, was 5 developed originally for northern Idaho and western Montana and has been adapted to currently include nine variants in the USA (Ritchie 1999). The Pacific Northwest coast variant incorporates different thinning regimes and is used along the coast of Oregon and Washington. While FVS can be applied to mixed-species stands, it requires tree crown measures, a variable not commonly measured in forest inventories. However, these measures can be estimated from other variables. California Conifer Timber Output Simulator (CACTOS) is a mixed-species conifer model developed for northern California (Wensel et al. 1986). CACTOS simulates information for one plot at a time; an additional program can be used to aggregate information across plots. CACTOS was created as a commercial product, but has been made available to the public free of charge. Stand Projection System (SPS), another commercial product, was developed for Douglas-fir stands in Oregon and Washington (Arney 1985). Unlike the other individual tree level models listed, SPS uses age to drive all projections in the simulator; if age is not known it is calculated using the tallest trees and the site index provided (Ritchie 1999). In uneven-aged stands, where age may be poorly defined or undefined, SPS may not provide accurate estimates. ORGANON was created for use in mixed-conifer forests located in Oregon (Hann et al. 2003). ORGANON contains components for thinning, fertilization, and pruning. ORGANON also requires tree crown ratio measurements and the fertilization component can only be simulated on even-aged Douglas-fir stands which are less than 80 years old (Ritchie 1999). The tree and stand simulator (TASS; Mitchell 1975) was created for coastal BC trees and is maintained by the British Columbia Ministry of Forests (Wang and Russell 2006). TASS is an individual tree distance- 6 dependent (ie, spatially explicit), three-dimensional growth simulator for even-aged, single-species coniferous stands. Growth and yield models are generally comprised of three main components: growth, mortality, and recruitment. The probability of tree-level mortality is commonly modelled using using logistic regression. (e.g., Hamilton 1986; Flewelling and Monserud 2002; Zhao et al. 2007). These models can then be applied using either a deterministic or stochastic approach. For individual tree-level models, growth is usually modeled as diameter, height, or basal area growth, or a combination of these. Diameter and height growth are commonly modelled using a log-linear model. Height growth can be modeled using either the growth potential independent method, where height is expressed directly as a function of tree, stand, and site variables, or the growth potential dependent method, where height is expressed as potential height reduced through the addition of a modifier function (Ek and Monserud 1974). Recruitment is often very difficult to predict. In the FORGE model, recruitment is based on light availability. Recruitment models are not specifically addressed in this study, but were addressed in an affiliated research project. 1.2.3 Mortality Tree mortality is a complex process resulting from interactions among biotic and abiotic factors (Franklin et al. 1987). Mortality can be broadly divided into catastrophic mortality, resulting in the death of many trees at one time (e.g., fire, epidemic insect/disease attacks, wind damage, etc.), versus regular mortality, where single or few trees die at one time (Lee 1971). Regular mortality can occur as a result of competition with other trees and plants (i.e., competition-based) or as a result of changes in microsite conditions, endemic insect and/or disease attacks, damage by animals, and windthrow. 7 Competition has been categorized as one-sided or asymmetric, when competition for above-ground resources such as light availability affects smaller trees more aggressively than larger trees, or two-sided or symmetric, when all trees regardless of size compete for below-ground resources such as water and soil nutrients (Cannell et al. 1984; Weiner 1990; Schwinning and Weiner 1998). Some commonly used measures of competition, such as basal area of live trees larger than the subject tree (BAL) measure only competition “from above”, in that only changes in larger trees result in changes in measured competition. Other measures, particularly tree-level measures, reflect competition with trees both larger and smaller than the subject tree, and may better reflect competition for some forest environments. For the mixed-conifer stands of the coastal Pacific Northwest, the horizontal structure is multi-layered and competition for light plays an important role in the mortality process (Spies et al. 1990). In terms of regular, competition-based mortality, models have been developed at either the stand or individual tree level. Often a logistic model is used to predict the probability of tree-level mortality or survival (e.g., Hamilton 1986; Flewelling and Monserud 2002; Zhao et al. 2007). Some studies have included tree-level characteristics to indicate competitive status such as the ratio of live crown length to tree height, termed crown ratio, and crown surface area as useful predictors of mortality (e.g., Monserud 1976; Silkstrom et al. 2005; Voelker et al. 2008). Basal area of larger trees has been used in other studies (e.g., Monserud and Sterba 1999; Bravo-Oviedo et al. 2006). Stand management can affect the rates of regular mortality through the alteration of both competition and site characteristics. Thinning results in an immediate decrease of inter-tree competition from the removal of nearby trees, thereby improving individual 8 tree vigor and reducing the rate of mortality for the residual trees (e.g., Curtis et al. 1981; Hann et al. 2003). Bailey and Tappeiner (1998) concluded that understorey trees had significantly higher live crown ratios compared to unthinned stands in thinned coastal Douglas-fir stands. O’Hara (1988) found that Douglas-fir trees with a medium-sized crown in thinned stands utilized their growing space as efficiently as trees with a large- sized crown in unthinned stands. The impacts of thinning on competition-based mortality have been incorporated into mortality models. For a single thinning event, thinning impacts have been incorporated through the development of two models, one for unthinned stands and one for thinned stands (Avila and Burkhart 1992), or by the inclusion of a thinning indicator variable in the model (e.g., Hann et al. 2003; Karlsson and Norell 2005). These methods do not account for multiple thinning events, high variation among thinning intensities, nor for differences in stage of stand development when thinning was applied, as is so often found with variable retention systems. To do so, mortality models developed with the inclusion of additional stand and/or tree-level variables which indicate changes in competition as a result of thinning are needed (e.g., Hamilton 1986). Fertilization is a common practice in the coastal forests of the Pacific Northwest of Canada and the USA, with the growth of Douglas-fir responding well to nitrogen fertilization (e.g., Brix 1981; Weetman et al. 1997). Shen et al. (2001) found that for Douglas-fir stands located in the inland Northwest, increasing rates of N application increased the probability of mortality for individual trees. In coastal BC, fertilization has been shown to accelerate mortality for Douglas-fir (Brix 1992). This could be due in part to increasing individual tree growth which in turn increases competition and accelerates 9 stand development, or due to increased ingrowth which also increases competition. Fertilization has commonly been incorporated into mortality models through the addition of a categorical variable (e.g., Hann et al. 2003). 1.2.4 Growth Diameter increment has been commonly modeled using a log-linear function (e.g., Cao 2000, Hann et al. 2003). Height growth has commonly been modeled using two methods: the growth potential independent method, where height is expressed directly as a function of tree, stand, and site variables, or the growth potential dependent method, where height is expressed as a function of the potential height which is then reduced through the addition of a modifier (Ek and Monserud 1974). Whereas the growth potential independent method uses all height data to model mean height for all trees, the growth potential dependent method provides estimates of population-averaged potential height values as defined by a subset of dominant and codominant trees. Similar to diameter increment, height increment has previously been modelled using as a log-linear function (e.g., Wykoff et al. 1982). Commonly, the effects of fertilization have been included in growth models through additional categorical variables. This approach was used by Carlson et al. (2008) to model mean relative growth rate of loblolly pine (Pinus taeda L.) and by Huber et al. (2009) for other growth models. Thinning has often been modeled through the inclusion of thinning parameters designed to describe the type of treatment received. For example, Wimberly and Bare (1996) added a parameter to their basal area increment model defined as the amount of basal area removed. Others have modeled thinning through the inclusion of state variables which change instantaneously after thinning as input variables 10 to predict growth following thinning. In particular, state variables which changed as a result of the thinning were basal area per ha, basal area of larger trees, and stems per ha. This approach was used for young Scots pine (Pinus sylvestris L.) stands by Pukkala et al. (2002), who found that models without additional thinning variables or parameters can produce unbiased predictions for thinned stands. The Von Bertalanffy (1957) function was developed to model the relationship between animal growth and metabolic functions. The function describes the net gain in an organism’s weight (W) as the difference in the processes which build up and break down materials used for growth. This relationship was reparameterized by Richards (1959) to describe cumulative growth which follows a sigmoidal curve, where growth begins at zero, increases to a maximum which occurs at the inflection point, begins to decrease and finally levels off to an asymptote. The Box-Lucas (1959) model, developed from the Von Bertalanffy equation, has been used to describe growth in fish (e.g., Carlson and Baremore 2005) and was used by Huang and Titus to model both diameter (Huang and Titus 1995) and height (Huang and Titus 1999) increment of white spruce (Pinus glauca (Moench) Voss). The flexibility of the Box-Lucas (1959) model makes it useful for modeling both diameter and height increment. 1.3 Gaps in knowledge and practice Current growth and yield models used in coastal BC are limited by either the number of species which can co-exist in a stand (e.g., the model TASS) or are limited by the need for data not commonly obtained in inventory databases such as tree crown ratio (e.g., the models FVS and ORGENON). Additionally, no growth and yield models have been developed to include variable retention systems as implemented on coastal BC 11 stands, where a variety of thinning intensities and spatial patterns, timing of thinning and fertilization treatments, and numbers of treatments are used. No modeling context is currently available to provide important information on various issues such as the mortality and growth impacts of variable retention systems. One reason this void exists is due to the lack of data. The dataset used in this study was quite extensive and allowed for analysis of actively managed- stands which vary greatly in age classes, species components, and treatment regimes. Unlike experimentally designed studies, this study was based on managed stands. This was a unique dataset in that the same methods and information were available for trees defined as ingrowth as well as trees defined as old-growth, for plots which were regenerated from plantings to plots which were regenerated naturally, and for stands which received no treatments to stands which received multiple treatments. To account for these variations, other growth and yield models have been created from a patchwork of different datasets collected at different times and locations (Hann et al. 2003), leading to issues of missing information and the creation of separate models and/or variables for unique situations (e.g., old-growth). The dataset used in this study allowed for the investigation of how different treatment regimes, including no treatment, affect the mortality and growth of trees located in highly variable stands throughout coastal BC. Most other studies in this area that have focused on treatment effects have utilized an experimental design approach, limiting the number of treatment regimes and creating specific scenarios which are difficult to repeat using variable retention systems (e.g., He and Barclay 2000). In this study, I used data from actively managed stands rather than experimental data to examine multiple interventions due to active management across a landscape. I 12 investigated how a variety of thinning and fertilization treatments can be modelled and how a model based approach can be used to examine how these treatments affect mortality and growth within a managed setting. This information is not found within the current literature. 1.4 Research goals In this research project, the growth and mortality of regenerating and older trees for coastal temperate rain forests of British Columbia were examined with the overall objective of modeling the impacts of fertilization, thinning, and variable retention on tree growth and mortality. This thesis focused on the following three topics: I. Mortality The first objective was to develop mortality models for second-growth western hemlock, Douglas-fir, and western redcedar trees growing on a variety of sites in coastal BC. The hypotheses were: a. A model developed solely for coastal BC would outperform existing models developed for similar species in other areas; and b. Inter-tree competition variables, such as basal area of larger trees, would be highly related to the probability of mortality for trees located in mixed-species stands with high structural diversity. The second objective was to evaluate the developed mortality models using three alternative implementation methods. The hypotheses associated with this objective were: 13 a. Little differences between approaches would be seen when results are aggregated to the stand level; however, differences would become more evident when analyzed within diameter classes; b. Differences would be more apparent for plots with higher mortality rates; and c. For the threshold probability method, the optimum threshold value would return the best results and should be determined individually for each species within a mixed-species stand. The third objective of this study was to model the impacts of thinning using the developed mortality models. The hypotheses associated with this objective were: a. The time since thinning would affect the reduction in mortality, which would increase to a maximum following thinning and then decline; b. Given a single thinning event, mortality would decline with thinning intensity, as a consequence of reduced competition; c. A greater mortality response would occur when thinning occurs at an earlier time in stand development; and d. Repeated thinning events would not result in repeated reductions in mortality. The fourth objective was to modify the developed mortality model to include the effects of fertilization. The hypotheses associated with this objective were: a. Given a single fertilization application, a greater mortality response would occur compared to untreated stands; 14 b. Multiple fertilization applications would not result in repeated reductions in mortality; and c. Impacts of thinning and fertilization on mortality can be represented additively. II. Diameter increment The first objective was to develop diameter increment models for second- growth western hemlock, Douglas-fir, and western redcedar trees growing on a variety of sites in coastal BC using the Box-Lucas model (1959) based on metabolic processes. The hypothesis was: a. A process-based diameter increment model would outperform an empirical diameter increment model and allow for the direct interpretation of explanatory variables. The second objective was to estimate the impacts of thinning, fertilization, and the combination of thinning and fertilization on diameter increment. The hypotheses were: a. Plots which were fertilized would experience an initial increase in diameter increment growth during the first few years following application only; b. Modifying only the state variables for inter-tree competition would reflect the initial increase in diameter increment growth during the first few years following thinning; and c. The increase in diameter increment growth due to the combination of thinning and fertilization can be represented additively. 15 III. Height increment The first objective of this study was to evaluate potential height increment and average height increment models using the Box-Lucas (1959) model. The hypotheses were: a. For potential height growth, the prediction interval approach would define the upper height growth values better than the binning approach where the tallest trees were selected within diameter classes; and b. Little difference would be seen between the growth potential dependent and independent methods for height growth models The second objective was to estimate the impacts of thinning, fertilization, and the combination of thinning and fertilization on height increment. The hypotheses were: a. Height growth for the fertilized plots would be greater than the untreated plots, and the magnitude of the difference would depend upon species and site characteristics; b. Modifying only the state variables for inter-tree competition would reflect the effect of thinning on height growth; c. Height growth for the thinned plots would be the same or less than the untreated plots, and the magnitude of the difference would depend upon species and site characteristics; d. The largest difference in height growth between treated and untreated stands would be found for stands which received fertilization and thinning in combination; and 16 e. No matter the treatment, height growth would approach but not reach the potential height growth. 1.5 Dissertation structure To address the objectives found within these three topics, four separate research chapters and a concluding chapter are presented in the following order: Chapter 2 – Mortality was addressed in the context of single and mixed- species stands. Because of irregular measurement intervals, survival was estimated using a generalized logistic model and mortality was calculated by subtraction. Models were developed for Douglas-fir, western redcedar, and western hemlock using untreated plots and a random coefficients modeling approach. For each species, an alternative model based on Hann et al. (2003) was also examined as a comparison. Implementation methods for survival models were also explored. A version of Chapter 2 was accepted on April 6, 2010 for publication in the Canadian Journal of Forest Research. Chapter 3 –This chapter investigates the impacts of thinning and fertilization on competition-based mortality for Douglas-fir using the mortality model developed in Chapter 2. The model was modified to include the effects of fertilization and thinning. The hypotheses for this chapter were examined using a modeling approach. A version of this chapter will be submitted for possible publication in Forest Ecology and Management. Chapter 4 – Diameter increment models were developed for Douglas-fir, western redcedar, and western hemlock using the Box-Lucas (1959) model. This Box-Lucas (1959) model was developed to describe the net gain in an 17 organism’s weight as the difference in the processes which build up and break down materials used for growth. A parameter prediction approach was used to develop models for Douglas-fir, western hemlock, and western redcedar for untreated stands. These models were then modified to examine the effect of fertilization and thinning. Different modifications were analyzed and no changes to the models were made for thinned stands. However, for fertilized stands, the models were modified to include categorical variables for fertilization. A version of this chapter has been submitted for possible publication in Ecological Modeling. Chapter 5 – In this chapter, the results for different methods for modeling mean height increment and potential height increment and how these two measures were related are presented. First, potential height increment models were developed for Douglas-fir, western redcedar, and western hemlock using the Box-Lucas (1959) model under alternative definitions of potential height. For each species, model accuracy was calculated and the different definitions and methods used in model development were explored. Secondly, different methods for determining mean height increment were also evaluated and a model was selected for further evaluation. Lastly, the selected mean height increment models were modified to include the effects of fertilization and thinning. A version of this chapter will be submitted for possible publication in Forest Science. Chapter 6 – In this concluding chapter, I summarize the major outcomes of this thesis and provide a discussion of the management and methodological 18 implications of this research. Included in the chapter is a discussion of future research opportunities. 19 1.6 References Arney, J.D. 1985. A modeling strategy for the growth projection of managed stands. Can. J. For. Res. 15(3): 511-518. Aubry, K.B., Halpern, C.B., and Peterson, C.E. 2009. Variable-retention harvests in the Pacific Northwest: a review if short-term findings from the DEMO study. For. Ecol. Manage. 258(4): 398-408. Avila, O.B. and Burkhart, H.E. 1992. Modeling survival of loblolly pine trees in thinned and unthinned plantations. Can. J. For. Res. 22(12): 1878-1882. Bailey, J.D. and Tappeiner, J.C. 1998. Effects of thinning on structural development in 40- to 100-year-old Douglas-fir stands in western Oregon. For. Ecol. Manage. 108(1- 2): 99-113. Bartelink, H. H. 2000. A growth model for mixed forest stands. For. Ecol. Manage. 13(1-3)4: 29-43. Box, G.E.P. and Lucas, H.L. 1959. Design of experiments in non-linear situations. Biometrika 46(1/2): 77-90. Bravo-Oviedo, A., Sterba, H., del Rio, M., and Bravo, F. 2006. Competition-induced mortality for Mediterranean Pinus pinaster Ait. and P. sylvestris L. For. Ecol. Manage. 222(1-3): 88-98. Brix, H. 1981. Effects of thinning and nitrogen fertilization on branch and foliage production in Douglas-fir. Can. J. For. Res. 11(3): 502-511. Brix, H. 1992. Fertilization and thinning effect on a Douglas-fir ecosystem at Shawnigan Lake: A synthesis of project results. For. Can. and B.C. Min. For., Victoria, BC FRDA Rep. No. 196. 20 Bunnell, F.L and Dunsworth, G.B. (eds). 2009. Forest biodiversity: learning how to sustain biodiversity in managed forests. UBC Press, Vancouver, Toronto 349 pp. Cannell, M.G.R., Rothery, P., Ford, E.D. 1984. Competition within stands of Picea sitchensis and Pinus contorta. Ann. Bot. 53: 349-362. Cao, Q.V. 2000. Prediction of annual diameter growth and survival for individual trees from periodic measurements. For. Sci. 46(1): 127-131. Carlson, C.A., Burkhart, H.E., Allen, H.L., and Fox, T.R. 2008. Absolute and relative changes in tree growth rates and changes to the stand diameter distribution of Pinus taeda as a result of midrotation fertilizer applications. Can. J. For. Res. 38(7): 2063- 2071. Carlson, J.K., and Baremore, I.E. 2005. Growth dynamics of the spinner shark (Carcharhinus brevipinna) off the United States southeast and Gulf of Mexico coasts: a comparison of methods. Fish Bull. 103: 280-291. Curtis, R.O., Clendenen, G.W., and DeMars, D.J. 1981. A new stand simulator for Douglas fir-DFSIM user’s guide. USDA For. Serv. Gen. Tech. Rep. PNW-GTR-128. Ek, A.R. 1974. Nonlinear models for stand table projection in Northern hardwood stands. Can. J. For. Res. 4(1): 23-27. Ek, A. R., and Monserud, R. A. 1974. FOR-EST: A Computer Model for Simulating the Growth and Reproduction of Mixed Forest Stands. Research Report A2635, College of Agricultural and Life Sciences, University of Wisconsin, Madison. Pages 1-14 + 3 Appendices. 21 Flewelling, J. W. and Monserud, R. A. 2002. Comparing methods for modeling tree mortality. In: USDA Forest Service Proceedings from the 2 nd Forest Vegetation Simulator Conference, RMRS-P-25, Fort Collins, CO, pp. 168-175. Franklin, J.F., Shugart, H.H., and Harmon, M.E. 1987. Tree death as an ecological process. BioScience 37(8): 550-556. Gavin, D.G., Brubaker, L., and Lertzman, K.P. 2003. Holocene fire history of a coastal temperate rain forest based on soil charcoal radiocarbon dates. Ecology 84(1) 186- 204. Getzin, S., Dean, C., He, F., Trofymow, J.A., Wiegand, K. and Wiegand, T. 2006. Spatial patterns and competition of tree species in a Douglas-fir chronosequence on Vancouver Island. Ecography 29(5): 671-682. Hamilton, D.A., Jr. 1986. A logistic model of mortality in thinned and unthinned mixed conifer stands of northern Idaho. For. Sci. 32(4): 989-1000. Hann, D.W., Marshall, D.D., and Hanus, M.L. 2003. Equations for predicting height-to- crown-base, 5-year diameter growth rate, 5-year height growth rate, 5-year mortality rate, and maximum size-density trajectory for Douglas-fir in the coastal region of the Pacific Northwest. Forest Research Laboratory, Oregon State University, Corvallis, Oreg. Res. Contr. 40. Haung, S. and Titus, S.J. 1995. An individual tree diameter increment model for white spruce in Alberta. Can. J. For. Res. 25(9): 1455-1465. Huang, S. and Titus, S.J. 1999. An individual tree height increment model for mixed white spruce-aspen stands in Alberta, Canada. For. Ecol. Manage. 123(1): 41-53. 22 He, F. and Barclay, H.J. 2000. Long-term response of understory plant species to thinning and fertilization in a Douglas-fir plantation on southern Vancouver Island, British Columbia. Can. J. For. Res. 30(4): 566-572. Huber, M., Halmschlager, E., and Sterba, H. 2009. The impact of forest fertilization on growth of mature Norway spruce affected by Sirococcus shoot blight. For. Ecol. Manage. 257(6): 1489-1495. Karlsson, K. and Norell, L. 2005. Modeling survival probability of individual trees in Norway spruce stands under different thinning regimes. Can. J. For. Res. 35(1): 113- 121. Korol, R.L., Milner, K.S., and Running, S.W. 1996. Testing a mechanistic model for predicting stand and tree growth. For. Sci. 42(2): 139-153. Landsberg, J.J., Waring, R.H., 1997. A generalised model of forest productivity using simplified concepts of radiation-use efficiency, carbon balance and partitioning. For. Ecol. Manage. 95(3): 209–228. Lee, Y. 1971. Predicting mortality for even-aged stands of lodgepole pine. For. Chron. 47(1): 29-32. Lertzman, K.P., Sutherland, G.D., Inselberg, A. and Saunders, S.C. 1996. Canopy gaps and the landscape mosaic in a coastal temperate rain forest. Ecology 77(4): 1254- 1270. Mitchell, K.J. 1975. Dynamics and simulated yield of Douglas fir. Forest Science. Monograph 17, Supplement to Volume 21, Number 4, pp. a0001-z0001(1) Mitchell, S.J. and Beese, W.J. 2002. The retention system: Reconciling variable retention with the principles of silvicultural systems. For. Chron. 78(3):397–403. 23 Monserud, R.A. 1976. Simulation of forest tree mortality. For. Sci. 22(4): 438-444. Monserud, R.A. and Sterba, H. 1999. Modeling individual tree mortality for Austrian forest species. For. Ecol. Manage. 113(2-3): 109-123. O’Hara, K.L. 1988. Stand structure and growing space efficiency following thinning in an even-aged Douglas-fir stand. Can. J. For. Res. 18(7): 859-866. Peterson, C.E. and Anderson P.D. 2009. Large-scale interdisciplinary experiments inform current and future forestry management options in the U.S. Pacific Northwest. For. Ecol. Manage. 258(4): 409-414. Pukkala, T., Miina, J., and Palahi, M. 2002. Thinning response and thinning bias in a young scots pine stand. Silvica Fennica. 36(4): 827-840. Richards, F.J. 1959. A flexible growth function for empirical use. J. Exp. Bot. 10(29): 290-300. Ritchie, M.W. 1999. A compendium of forest growth and yield simulators for the Pacific coast states. USDA For. Serv. Gen. Tech. Rep. PSW-GTR-174. Sands, P.J., and Landsberg, J.J. 2002. Parameterisation of 3-PG for plantation grown Eucalyptus globulus. For. Ecol. Manage. 163(1-3): 273–292. Schwinning, S. and Weiner, J. 1998. Mechanisms determining the degree of size asymmetry in competition among plants. Oecologia 113(4): 447-455. Shen, G., Moore, J.A., and Hatch, C.R. 2001. The effect of nitrogen fertilization, rock type, and habitat type on individual tree mortality. For. Sci. 47(2): 203-213. Silkstrom, U., Jansson, G., and Weslien, J. 2005. Predicting the mortality of Pinus sylvestris attacked by Gremmeniella abietina and occurrence of Tomicus piniperda colonization. Can. J. For. Res. 35(4): 860-867. 24 Smith, N.J. 2005. Variable retention, microclimate, experimental and modeling projects in TFL 39: Modeling mortality in variable retention stands using the FORGE model. Unpublished internal report. Nanaimo, BC, Canada. Smith, N.J. 2010. Young tree growth dynamics in variable retention and clearcut settings in coastal British Columbia in FSP: Y10326 Final Grant Report. Nick Smith Forest Consulting, Nanaimo, BC. Spies, T.A., Franklin, J.F., and Klopsch, M. 1990. Canopy gaps in Douglas-fir forests of the Cascade Mountains. Can. J. For. Res. 20(5): 649-658. Vanclay, J. K. 1994. Modeling forest growth and yield, applications to mixed tropical forests. CAB International Publishing. Wallingford, UK. Voelker, S.L., Muzika, R., and Guyette, R.P. 2008. Individual tree and stand level influences on the growth, vigor, and decline of red oaks in the Ozarks. For. Sci. 54(1): 8-20. Von Bertalanffy, L. 1957. Quantitative laws in metabolism and growth. The Q. Rev. Biol. 32(3): 217-231. Wang, T. and Russell, J.H. 2006. Evaluation of selfing effects on western redcedar growth and yield in operational plantations using the tree and stand simulator (TASS). For. Sci. 52(3): 281-289. Weetman, G.F., Prescott, C.E., Kohlberger, F.L., and Fournier, R.M. 1997. Ten-year growth response of coastal Douglas-fir on Vancouver Island to N and S fertilization in an optimum nutrition trial. Can. J. For. Res. 27(9): 1478-1482. Weiner, J. 1990. Asymmetric competition in plant populations. Trends Ecol. Evol. 5(11): 360-364. 25 Wensel, L.C., Daugherty, P.J., and Meerschaert, W.J. 1986. CACTOS user’s guide: the California Conifer Timber Output Simulator. Calif. Agric. Exp. Sm. Bull. 1920. Wimberly, M.C. and Bare, B.B. 1996. Distance-dependent and distance-independent models of Douglas-fir and western hemlock basal area growth following silvicultural treatment. For. Ecol. Manage. 89(1-3): 1-11. Wykoff, W.R., Crookston, N.L., Stage, A.R. 1982. User’s guide to the stand prognosis model. USDA For. Serv. Gen. Tech. Rep. INT-133. Zhao D., Borders, B., Wang, M., and Kane, M. 2007. Modeling mortality of second- rotation loblolly pine plantations in the Piedmont/Upper Coastal Plain and Lower Coastal Plain of the southern United States. For. Ecol. Manage. 252(1-3): 132-143. 26 2. Modeling mortality for mixed-species stands located in coastal British Columbia 1 2.1 Introduction The mixed-species coastal temperate rain forest of British Columbia (BC) is noted for its lack of recent fires (Lertzman et al. 1996). The successional processes of the Pacific Northwest forests are defined by small-scale disturbances, mainly due to windthrow (Gavin et al. 2003). Hence, gap dynamics define the forest structure and natural successional processes on Vancouver Island (Lertzman et al. 1996), with Douglas-fir (Pseudotsuga menziesii var. menziesii (Mirb.) Franco) being the dominant, pioneer species and western hemlock (Tsuga heterophylla (Raf.) Sarg.) and western redcedar (Thuja plicata Donn) being the late-successional species aggregated within gaps (Getzin et al. 2006). While mortality in younger stands is typically due to self-thinning processes, as a stand matures, gaps formed from single-tree mortality drive the structural and compositional characteristics instrumental in stand development (Lertzman et al. 1996). A survival model, or conversely a mortality model, is an important component within growth and yield models (Avila and Burkhart 1992). Tree mortality is a complex process resulting from interactions among biotic and abiotic factors and many of the ecological processes and systems which contribute to tree mortality are not easily understood (Franklin et al. 1987). Factors affecting tree mortality include: microsite nutrient and climatic conditions, stand-level site and climatic conditions, inter-tree 1 A version of this chapter has been accepted for publication. Rathbun, L.C., LeMay, V., and Smith, N.J. 2010. Modeling mortality for mixed-species stands located in coastal British Columbia. Canadian Journal of Forest Research 40(8): 1517-1528. 27 competition, age, incidence of insect and/or disease attack, and fire. Mortality has classically been divided into catastrophic mortality, involving the death of many trees at one time (e.g., due to fire, epidemic insect/disease attacks, wind damage, etc.), and regular mortality, where single or few trees die at one time (Lee 1971). Regular mortality can occur as a result of competition with other trees and plants (i.e., competition-based), as a result of endemic pests or diseases, or as a result of changes in microsite conditions. Commonly, the mortality component of a growth and yield model considers only regular mortality; however, it should be noted that regular and catastrophic mortality may not always be independent (Dobbertin and Biging 1998). The status of a tree as live or dead is a binary response variable. The probability of individual tree survival has commonly been modeled using a logistic model (e.g., Zhao et al. 2004). The advantage of predicting tree survival probability rather than mortality is that tree survival over a q-year period is the annual probability to the q th power, requiring no interpolation of data with varied measurement intervals (Hamilton and Edwards 1976). Once the probability of survival has been estimated, a deterministic or a stochastic approach can be implemented to predict the status of trees at the end of a growth period. Each measured tree in a sample plot represents a number of stems per ha. For both approaches, the stems per ha represented by each tree at the end of a growth period is reduced based on the survival probability of each tree. There are two commonly used deterministic approaches. The first is a probability multiplier approach, where the predicted probability of survival is multiplied by the stems per ha at the beginning of the growth period to obtain the surviving stems per ha at the end of the period. The second 28 deterministic approach uses threshold probability values, where the stems per ha represented by the tree are predicted as live if the predicted survival probability is greater than a specified threshold, else they are predicted as dead (Monserud 1976). Threshold values can be subjective. However, an objective approach can be used to optimize the balance between the percent correctly predicted for the event (sensitivity) versus the percent correctly predicted for the non-event (specificity). For the stochastic approach, the predicted probability of survival is compared to a generated uniform random number from the open interval 0 to 1. If the estimated survival probability for the tree is greater than the random number, the stems per ha associated with the tree are predicted as live; otherwise they are predicted as dead at the end of the growth period. The stochastic approach is generally applied to an individual tree multiple times, called trials, and the average surviving stems per ha is calculated. As the number of trials is increased, the stochastic approach results in the same surviving stems per ha as the probability multiplier approach, given a true random number generator (Ek 1980, Weber et al. 1986). The first objective of this study was to develop survival models for second-growth western hemlock, Douglas-fir, and western redcedar trees growing on a variety of sites in coastal BC. The hypotheses associated with this first objective were: 1) a model developed solely for coastal BC would outperform existing models developed for similar species in other areas; 2) inter-tree competition variables, such as basal area of larger trees, would be highly related to the probability of mortality for trees located in mixed- species stands with high structural diversity. The second objective was to evaluate the developed mortality models using the three alternative implementation methods. The hypotheses associated with this objective were: 1) little differences between approaches 29 would be seen when results are aggregated to the stand level; however, differences would become more evident when analyzed within diameter classes; 2) differences would be more apparent for plots with higher mortality rates; and 3) of the threshold probabilities, the optimum threshold value would return the best results and should be determined individually for each species within a mixed-species stand. 2.2 Materials and methods 2.2.1 Study area The study area was located on Vancouver Island, Haida Gwaii (the Queen Charlotte Islands), and along the coastal mainland within the Coastal Western Hemlock (CWH) and Coastal Douglas-fir (CDF) Biogeoclimatic Ecological Classification (BEC) zones of BC (Meidinger and Pojar 1991). The CWH BEC zone is divided into multiple subzones: a very dry maritime subzone in the east, a moist maritime subzone in the central area, and a very wet maritime and hypermaritime subzone in the west (Meidinger and Pojar 1991). The temperatures in the CWH BEC zone range from 5.2 to 10.5˚ C, with a mean annual temperature of 8˚ C (Meidinger and Pojar 1991). The CDF BEC zone includes the rainshadow of Vancouver Island and the Olympic mountains (Meidinger and Pojar 1991). It is defined by warm, dry summers and mild, wet winters with mean annual temperatures ranging from 9.2 to 10.5˚ C, and a minimum temperature ranging from -21.1 to -11.7˚ C (Meidinger and Pojar 1991). Study plots were located within latitudes ranging from 48.37 to 53.53˚ N and longitudes ranging from 122.36 to 132.59˚ W. The majority of this area is comprised of second growth multi-species stands, regenerated both naturally and from planting. Common tree species include: Douglas-fir, western redcedar, western hemlock, red alder (Alnus rubra Bong.), Sitka 30 spruce (Picea sitchensis (Bong.) Carr.), and yellow cedar (Chamaecyparis nootkatensis (D.Don) Spach). 2.2.2 Data Permanent sample plot (PSP) data were provided by Island Timberland, LP. The database contained 1,455 untreated plots, ranging in size from 0.008 to 0.253 ha (Table 2-1). Measurements were collected between 1932 to 2002, with measurement intervals varying from 1 to 17 years with an average of 5.0 years. The majority of plots, 1,390, were defined as second growth, less than 150 years since stand establishment; 20 plots were defined as old growth, greater than 150 years since stand establishment; and 45 plots had no available age information. At the time of plot establishment, 528 plots were predominately western hemlock, 515 plots were predominantly Douglas-fir, and 97 plots were predominantly western redcedar, as defined by greater than 50% stems per hectare. Based on more than 50% of basal area per ha, 474 plots contained predominately western hemlock, 580 plots predominantly Douglas-fir, and 96 plots predominantly western redcedar. For the remaining plots, species composition averaged 18, 6, and 5% western hemlock, Douglas-fir, and western redcedar, respectively, based on basal area per ha. Densities at plot establishment ranged from 13 to 11,750 live trees per ha and site index, using coastal Douglas-fir site index at a base age 50 years measured at breast height (1.3 m above ground), ranged from 6.2 to 52.8 m. Average annual plot mortality for the first measurement period following plot establishment ranged from 0 to 11.11% with an average of 1.37%. The average annual mortality for the second measurement period was 1.31%, with a slightly wider range of 0 to 12.07%. 31 At each measurement, the species, tree status (i.e., live or dead), and diameter outside bark at 1.3 m above ground ( DBH, cm) were recorded for all trees with DBH greater than 4 cm and within plot boundaries. Height (ht, m) was measured for a subset of trees, and remaining heights were estimated using height-diameter functions. At the time of plot establishment, the average DBH values were 13.9, 12.5, and 12.6 cm for western hemlock, Douglas-fir and western redcedar, respectively, with corresponding maximum values of 168.3, 149.9, and 217.1 cm (Table 2-2). The average heights for western hemlock, Douglas-fir and western redcedar were 14.2, 11.7, and 11.2 m, respectively. Since the heights were a mixture of measured and estimated heights, the minimum height was 0 rather than the expected value of ≥ 1.3m. For model fitting, validation, and testing, all non-overlapping repeated measures for all plots were used. 2.2.3 Model selection To model the probability of individual tree survival at the end of a growth period, the status of a tree was defined as 1 for live and 0 for dead. A logistic model was used to predict the annual probability of survival. To allow for a variety of measurement periods, the logistic equation was extended for a variable number of years, q, to obtain the following equation for the predicted probability of survival at the end of the projection period (e.g., Zhao et al. 2007): [1] q Xf-p(s) exp1 where f(X) is a linear function of the predictor variables (X) at the beginning of the projection period and q is the number of years of the projection period. When selecting predictor variables, biological processes which affect the probability of mortality were considered. As a result, predictor variables representing 32 tree size and stage of development, site productivity, and inter-tree competition were included. Tree size and stage of development were represented by DBH (cm) and height (m). To represent site productivity, coastal Douglas-fir site index (SI, m) at a base age 50 was included. Also, the variable Growth Effective Age (GEA, years) was used to indicate tree productivity, where GEA is the age of a dominant tree of the same height as the tree of interest, for a given site index (Hann and Ritchie 1988). For these data, GEA was estimated from site index equations for the leading tree species in each plot. A number of stand level competition measures were considered. Curtis’ (1982) Relative Density, an index measure of density, was calculated as: qd G RD Curtis n BHD d n 1i 2 i q where DBHi (cm) is the DBH for tree i, dq (cm) is the quadratic mean diameter for the plot, G (m 2 /ha) is the stand basal area per hectare, and n is the number of live trees across all species within a plot. Stand basal area has been commonly used to measure competition for below ground resources (Fan et al. 2006). In addition, tree-level competition measures such as the basal area of larger trees (BAL, m 2 /ha), relative DBH (RDBH), and crown competition factor of larger trees (CCFL) were evaluated for the model. BAL was calculated as: n i ii BAtree δ 1 x BAL where δi is an indicator variable (1 if tree i has a DBH greater than the tree of interest, 0 otherwise), BAtreei (m 2 /ha) is the basal area per hectare value for tree i, and n is the 33 number of all live trees across all species within a plot. RDBH is a ratio of the diameter of the tree of interest to the average diameter within a plot and was calculated as: DBH BHD iRDBH where DBHi (cm) is the DBH of tree i and dbh (cm) is the average DBH for the plot. CCFL (expressed as a percent) was calculated as: n i ii δ 1 MCA xCCFL 400 / )WĈπ(MCA 2ii ii DBH 12.0 ˆ b2.54a WC where MCAi is the maximum crown area for tree i expressed as a percentage of a hectare that can be occupied by the maximum crown of tree i with DBHi in cm (Avery and Burkhart 2002, pp. 327-328), iWĈ (m) is the crown width of tree i, and the parameters a and b are species-specific constants in imperial units found in Smith (1966). All other constants are used for unit conversion. Using subsets of these possible predictor variables and Eq. [1], twelve candidate models (Models 1 through 12) were formed. The data for each species were divided into three subsets. Then, PROC NLMIXED of SAS, version 9.1.3 (SAS Institute Inc. 2004) was used to fit each candidate model using each subset by species. Initial parameter estimates were varied in order to achieve a global minimum. For each fitted model, Akaike’s Information Criteria (AIC) value was calculated using as: p22logL(AIC )ˆ 34 where )logL( ˆ is the log likelihood function, p is the number of parameters in the model, and smaller AIC values indicate a more accurate and/or more parsimonious model. Model selection for each species was based upon AIC values averaged for the three data subsets. The selected model was then compared to two models from the literature, developed for mixed-conifer stands comprised of similar species. The first model was developed by Temesgen (2002) for cedar-hemlock and Douglas-fir found in the interior of British Columbia. The Temesgen model is: SPHbGbBALbRDBHbDBHbDBHb DBH 1 bbf(X) 7654 2 3210 where SPH is the surviving stems per hectare at the beginning of the period, 70 b tob are parameters to be estimated, and the other variables are as previously defined. The second model was developed by Hamilton (1986) for mixed-conifer stands in Idaho: DBH dinc bRDBHb DBH 1 bdincbGbDBHbbf(X) 6543210 where 60 b to b are parameters to be estimated, dinc (cm) is the incremental diameter growth in the previous measurement period, and all other variables are as previously defined. To compare the three models, the three subsets of data were combined into model (2/3) versus test (1/3) datasets by species. For the selected, Temesgen, and Hamilton models, PROC NLMIXED of SAS, version 9.1.3 (SAS Institute Inc. 2004) was used to fit the models using the model dataset for each species. Each model was then validated using the test data. This process was repeated three times varying the one-third of the data reserved for testing. For each species, AIC values from model fitting and percent concordance using the reserved test data were then compared among the three 35 models. Concordance compares the predicted probabilities of two trees of different status (one live and one dead) to determine whether the order of the predicted probabilities corresponds with that expected based upon tree status. Each pair of one dead tree with one live tree was counted as concordant if the predicted survival probability for the live tree was greater than that for the dead tree. The number of concordant pairs was compared to the total number of pairs to obtain percent concordance. Each model was then fitted using the entire dataset by species, and the optimum probability threshold was calculated for all trees and also for small trees (< 7.5cm DBH) and large trees (≥ 7.5cm DBH) separately. Using the optimum, and two subjectively selected probability threshold values, sensitivity and specificity values were calculated. One model was finally selected for each species. This model was then used for the second objective of evaluating three possible model implementation approaches for accuracy in estimating the surviving stems per ha. 2.2.4 Model implementation Tree and plot measurements at plot establishment (i.e., the beginning of the first measurement period) were used as inputs to the selected model to predict the probability of survival to the end of the first measurement period (i.e., the beginning of the second measurement period). Once the probability of survival was estimated, three implementation approaches were used to predict tree status at the end of the first measurement period as follows: 1. Probability multiplier approach. The probability of survival was used to estimate the surviving stems per ha at the end of the period by: SPHq = p(s)* SPH0 36 where SPHq is the stems per ha surviving to the end of the measurement period represented by the tree, p(s) is the predicted probability of survival for the tree, and SPH0 is the stems per ha represented by the tree at the beginning of the measurement period. 2. Threshold approach. Two subjective threshold values were chosen, and an optimum threshold value was calculated. Many different methods for determining optimum threshold values exist. Cohen’s (1960) optimal threshold, k , was used and calculated as: c c0 p pp 1 k where po is the proportion of units in which the predicted and actual outcomes agree and pc is the proportion of units for which agreement is expected by chance (Fleiss 1975). The optimum threshold value, k , was determined using all plot measurements by species. The probability of survival and the specified threshold value were then used to estimate the surviving stems per ha at the end of the period by: otherwise 0 thresholdp(s) if SPH SPH 0 q where all variables are as previously defined. 3. Stochastic approach, the surviving stems per ha for a single trial were estimated by: otherwise 0 [0,1]number randomp(s) if SPH SPH 0 q 37 where all variables are as previously defined. For more than one trial, the estimated surviving stems per ha was averaged. This approach was run using 1, 5, 20, and 1,000 trials. For each of combination of species and implementation approach, the mean bias and root mean-squared error (RMSE) values were calculated for stems per ha that survived, summed for all trees in a plot. Mean bias was calculated as: m yy m 1j jj ˆ BiasMean where jy is the actual number of trees per ha that survived within the q period for plot j, jŷ is the predicted number of trees per ha that survived within the q period for plot j, and m is the number of plots. Root mean-squared error was calculated as: 1pm yy m 1j 2 jj ˆ RMSE where p is the number of predictor variables in the model and all other variables are as previously defined. Mean bias and RMSE values were also calculated by 5 cm DBH classes and by mortality classes for each species. Four mortality classes were defined based upon the actual average annual percent mortality that occurred at the second measurement period: 1, defined as plots where no mortality occurred; 2, defined as plots with less than 1% average annual mortality; 3, defined as plots with an average annual percent mortality between 1 and 4%; and 4, defined as plots with an average annual percent mortality greater than 4%. 38 2.3 Results 2.3.1 Model selection The best combination of predictor variables for western hemlock and western redcedar was DBH, height, basal area of larger trees, Curtis’ (1982) relative density, site index, growth effective age, and basal area per hectare (Model 1): GbGEAbSIbCurtisRDbBALbheightbdbhbbf(X) 76543210 This model produced the smallest averaged AIC values for western hemlock and western redcedar (Table 2-3). While Models 2 and 3 returned smaller AIC values for Douglas-fir, these AIC values were similar to that of Model 1. Models 8, 10, and 12 were the simplest models, containing DBH, height, and site index. Model 8 also included BAL and returned smaller AIC values for western hemlock and Douglas-fir than Models 10 and 12. Model 10, which included G, outperformed Models 8 and 12 for western redcedar. As well as having the best overall AIC values for western hemlock and western redcedar of the models examined, Model 1 was biologically logical in that the predictor variables included measures of tree size and stage of development, site productivity, and inter-tree competition. Model 1 was then compared to the Temesgen and Hamilton models. The AIC values for Model 1 were the smallest for all three species (Table 2-4). The average AIC value returned for Model 1 (-8,745) was twice as small as those for the Temesgen and Hamilton models (-4,256 and -4,089 respectively). The average AIC value found from Model 1 for western redcedar was negative as compared to the positive values found for the Temesgen and Hamilton models. The percent concordances for Douglas-fir and 39 western hemlock were similar for all three models, but Model 1 had the highest percent percent concordance. For all three species, the threshold values where specificity and sensitivity intersected were quite different for small diameter trees than for large diameter trees (Figure 2-1). A threshold value between 0.79 and 0.82 produced a balance between specificity and sensitivity for small diameter western hemlock; for Douglas-fir, this threshold value was between 0.80 and 0.85; and for western redcedar, this threshold value was between 0.85 and 0.90. All three models indicated that the intersection between specificity and sensitivity was greater than 0.90 for large Douglas-fir and western redcedar. For large western hemlock trees, the Temesgen model returned an intersection value slightly less than 0.90, but Model 1 and the Hamilton model returned a value between 0.90 and 0.95. Based on the model comparisons, Model 1 was selected and used in assessing accuracy of alternative model implementation approaches. 2.3.2 Model implementation Regardless of the implementation approach, the lowest RMSE values were obtained for western redcedar (Table 2-5). Douglas-fir and western hemlock returned similar mean bias values for the stochastic approach, but Douglas-fir had larger RMSE values. For each species, the stochastic approach returned varying results based on the number of trials used. As the number of trials increased, the stochastic approach resulted in mean bias and RMSE values that were similar to the probability multiplier approach, as expected. The two subjective threshold values chosen were 0.50 and 0.75. Using small and large trees combined, Cohen’s optimum threshold value was 0.76 for western 40 hemlock, 0.78 for Douglas-fir, and 0.79 for western redcedar. For western hemlock, the threshold approach resulted in the greatest mean biases and RMSE values, regardless of the threshold value. For all three species, the mean bias was negative for a threshold value of 0.50, meaning the surviving stems per ha were overestimated, and positive for the optimum threshold value for all but western redcedar, indicating underestimation. When comparing implementation approaches, the greatest variation in mean bias and RMSE values by DBH class was obtained for western hemlock and western redcedar, whereas the least variation in mean bias and RMSE was obtained for Douglas- fir (Figure 2-2). Large differences in both mean bias and RMSE were obtained for diameter classes less than 10 cm. The model under-predicted survival (i.e., positive mean biases) for diameter classes less than 7.5cm for all methods and species except when a threshold value of 0.50 was used. For small trees, the threshold approach produced large mean bias and RMSE values which varied greatly depending upon the threshold value, and differed from the other methods (as represented by the three lines furthest from zero). All approaches converged and predicted accurately for trees greater than 15cm DBH. Even with just one trial, the stochastic approach returned results similar to the probability multiplier approach for these trees. When plots were classified by mortality rate, all three approaches predicted well for Douglas-fir located on plots with an average annual mortality less than 4% (Figure 3b); for an average annual mortality greater than 4%, the model over-predicted the surviving stems per ha for all approaches except the optimum threshold. Similar patterns were found for western hemlock (Figure 2-3a) and western redcedar (Figure 2-3c), with the most variation occurring within the threshold approaches. Using the 0.50 and 0.75 41 threshold values resulted in over-predicting the surviving stems per ha for both Douglas- fir and western redcedar on plots with an average annual mortality greater than 4%, while the optimum threshold value under-predicted. Similar to the other two species, differences among the approaches were large for western redcedar on plots where average annual mortality was greater than 4%. 2.4 Discussion Stand development following large-scale disturbances is generally defined by four stages: stand initiation, stem exclusion, understorey re-initiation, and old growth (Oliver 1981). Stands within the first three stages of stand development are often used to develop survival models in areas where fire or other stand-replacing disturbances are prevalent, since there is often a lack of old growth data. Within coastal BC, stand structures vary depending upon the level of human and natural disturbances. For our study sites, regeneration on managed sites has occurred both naturally and from planting, resulting in varied species composition from site to site. Survival is affected by this diversity in structure and composition, and consequently, survival models developed solely for these areas are needed to forecast future conditions. While statistical properties and fit statistics are important in model selection, if a model is not biologically sound it will be difficult to make future predictions or to apply to other datasets (Yang et al. 2003). Both the Temesgen and the Hamilton models developed for similar species were largely driven by DBH. While DBH is a good indicator of growth and vigor, given site quality and age, other variables which represent canopy location and competition for resources become important in areas with complex structural and compositional components. In this study, the selected model includes 42 variables which represent tree size and stage of development, site productivity, and inter- tree competition. Douglas-fir on the BC coast is a shade-intolerant pioneer species, whereas western hemlock and western redcedar are shade-tolerant and can establish at the same time as Douglas-fir, or as late successional species (Carter and Klinka 1992). Coastal Douglas-fir seldom regenerates under a closed canopy, but can regenerate in large gaps in the driest subzones in the CWH zone and the CDF zone (Getzin et al. 2006). Where Douglas-fir does regenerate in smaller gaps, trees tend to have poor growth and small crowns, indicating suppression. The data for this study were gathered from a majority of stands following harvesting. However, the data represent a wide range of mixtures of the three species, including plots in that were primarily Douglas-fir, western hemlock, and western redcedar, as well plots that were mixtures of the three species. Differences in shade tolerance and location within the canopy structure impact the level of competition for both above and below ground resources such as light, water, and nutrients (Carter and Klinka 1992). To indicate this, variables such as BAL and stand basal area were included in the survival model. BAL has been used as a measure of one- sided competition, competition for above ground resources (Bravo-Oviedo et al. 2006). This combination of competition variables performed the best for western hemlock and western redcedar but not for Douglas-fir. For Douglas-fir, a more accurate model, as defined by AIC, was found when BAL, a competition index, was replaced by relative diameter, a measure of size. BAL values for these dominant trees tended to be quite small and varied little from tree to tree. This suggests that benefit of competition variables within a mortality model may be species dependent in mixed-species stands. 43 The selected survival model outperformed the Temesgen and Hamilton models developed originally for interior BC and Northern Idaho, respectively. Unlike more continental regions in North America, coastal BC receives high amounts of rainfall. Carter and Klinka (1992) found that a difference in moisture affects the level of shade tolerance for a species. For example, Douglas-fir is classified as semi-shade tolerant on moderately dry sites (such as in the interior of BC east of the coastal mountain ranges) and as shade intolerant on fresh and moist sites. Shade tolerance levels directly affect canopy position, thereby indirectly affecting rates of mortality. All three models yielded poor results for trees less than 7.5 cm in DBH. These small trees are generally located within the understorey in mixed-species stands. They tend to have higher mortality rates than their dominant counterparts because of increased competition for both above and below ground resources. For Douglas-fir in mature and old growth forests located in the Cascade Mountains of Washington and Oregon, Spies et al. (1990) found that the annual percent mortality for mature trees located in the upper canopy was 0.78% compared to 2.54% for those in the understorey. Differences in shade tolerance impact the level to which mortality increases for these small trees. Shade- tolerant species, such as western hemlock and western redcedar, can persist under suppressed growth for long periods of time, whereas shade-intolerant species, such as Douglas-fir, cannot (Kobe and Coates 1997). These differences in mortality for small trees result in poor accuracies for mortality models if all sizes are included. The development of separate mortality models may improve predictions for these smaller trees. 44 The mortality models did not work well for the 90 plots which experienced larger actual average annual mortality (greater than 4%). The average size of the trees for these higher mortality plots were generally larger than on plots with lower mortality, with average DBH values of 14.7, 15.7, and 13.3 cm and mean heights, 16.4, 16.4, and 12.7 m for western hemlock, Douglas-fir, and western redcedar, respectively. These plots also averaged more stems per hectare and basal area per hectare; however, there were no large differences in site index relative to plots with lower mortality rates. The high mortality plots were comprised mainly of western hemlock (61.4% on average), and the majority (71.1 %) were established in the 1950’s and 1960’s. These plots were remeasured during the early 1970’s. In 1966 through 1968, coastal BC experienced lower than average rainfall for the summer months. Studies have shown that mortality increases for western hemlock when moisture levels are low in the summer months (Kotar 1972). These dramatic changes in site condition due to droughts or other catastrophic events increase the regular rate of mortality, making them difficult to include in mortality models. It is for this reason that the majority of mortality models only model regular mortality and predict poorly when applied to abnormal situations. Differences between the three implementation approaches were apparent at all levels: by plot, by DBH class, and by mortality class. As expected, the probability multiplier and stochastic approaches had nearly identical values for larger numbers of trials, while the threshold approach produced the greatest variability. A threshold approach is similar to a one-run stochastic approach in that an individual tree and all the trees it represents on a per hectare basis are predicted as either all dead or all alive at the 45 end of the growth period. For this reason, different threshold values, whether chosen subjectively or mathematically, produce large differences in outcomes. Cohen’s optimum threshold value was similar, but different, for all three species. Threshold values identify the point below which mortality (or survival) is relatively unaffected by variation in the factors which contribute to it (Deubner et al. 1980). Since variables which contribute to mortality are different for different species, threshold values should be set for each species. Predicting mortality for small trees less than 7.5 cm in DBH produced the largest differences among the approaches overall and the greatest differences among the threshold values for all three species. This indicates that the threshold value should decrease with decreasing diameter classes regardless of the species. Small trees tend to experience greater mortality rates, which should be reflected in the threshold values chosen, another indication that the development of separate mortality models for smaller trees is needed. Mean bias and RSME values increased greatly for plots which experienced an average annual mortality greater than 4% no matter what method was used. The probability multiplier approach returned the best results overall and is recommended. This is equivalent to a stochastic approach averaged over many repetitions, but with much less processing time. 2.5 Conclusion Individual tree survival models are a fundamental component within any tree- based growth and yield model. An understanding of the factors affecting mortality should be considered when developing mortality models; this is especially true for mixed-species stands with complex structural diversity. When implementing survival 46 models, different implementation approaches such as the probability multiplier, stochastic, or threshold approach likely yield different results. Given the importance of growth and yield modeling in designing silvicultural practices, a basic understanding of how a survival model works is critical and will ultimately determine which implementation approach is used. In this study, the use of a generalized logistic survival model resulted in accurate estimates for larger trees, but poor results for smaller trees. Further research on modeling survival and mortality across all size classes, perhaps using a process-modeling approach, is needed. Also, using a probability multiplier approach is recommended for implementing a mortality model. This is equivalent to many repetitions of the stochastic approach but with much less processing time. 47 Figure 2-1. Sensitivity and specificity by threshold value for the selected model and the two validation models for each species by small (<7.5cm) and large (≥7.5cm) diameter sizes: a) small western hemlock, b) large western hemlock, c) small Douglas-fir, d) large Douglas-fir, e) small western redcedar, and f) large western redcedar. 48 Figure 2-2. Bias and root mean-squared error (RMSE) by DBH class for the final model using all implementation approaches for each species: a) and b) western hemlock, c) and d) Douglas-fir, and e) and f) western redcedar. 49 Figure 2-3. Bias and root mean-squared error (RMSE) by mortality class for the final model using all implementation approaches for each species: a) and b) western hemlock, c) and d) Douglas-fir, and e) and f) western redcedar. 50 Table 2-1. Plot summary statistics at plot establishment and the end of the first measurement period (m=1,455 plots). Variable a Establishment End of 1 st Measurement Period Mean Minimum Maximum Mean Minimum Maximum G (m 2 /ha) 38.7 0.1 185.0 43.2 0.5 189.2 SPH 2,044 13 11,750 2,103 247 11,725 CurtisRD 8.7 0.1 24.2 9.6 0.3 24.1 SI (m) 29.3 6.2 52.8 29.4 10.6 51.4 a G is basal area per hectare, SPH is stems per hectare, CurtisRD is Curtis’ Relative Density, and SI is Site Index 51 Table 2-2. Tree summary statistics at plot establishment by species. Variable a Western hemlock Douglas -fir Western redcedar Number of trees 49,092 33,960 11,106 DBH (cm) Mean 13.9 12.5 12.6 Minimum 4 4 4 Maximum 168.3 149.9 217.1 Height (m) Mean 14.2 11.7 11.2 Minimum 0.0 0.0 0.0 Maximum 58.4 55.8 54.1 BAL (m 2 /ha) Mean 34.6 18.4 33.1 Minimum 0.0 0.0 0.0 Maximum 184.9 165.4 134.0 GEA (years) Mean 21.8 17.4 19.5 Minimum 0.0 0.0 0.0 Maximum 148.9 150.0 146.1 a DBH is diameter at breast height, BAL is basal area of larger trees, and GEA is growth effective age 52 Table 2-3. AIC values for selection models averaged across the three model- building data subsets. Model Variables included in model a Western hemlock Douglas- fir Western redcedar 1 DBH, height, CurtisRD, SI, GEA, G, BAL 16,744 -3,900 -892 2 DBH, height, CurtisRD, SI, GEA, G, RDBH 17,065 -4,202 -854 3 DBH, height, CurtisRD, SI, GEA, G, BAL, RDBH 16,825 -4,404 -813 4 DBH, height, CurtisRD, SI, GEA, RDBH 18,332 -3,809 -107 5 DBH, height, CurtisRD, SI, GEA, CCFL 19,322 -2,852 -6 6 DBH, height, CurtisRD, SI, GEA, SPH 19,832 -2,146 -47 7 DBH, height, SI, GEA, BAL 19,795 -2,802 518 8 DBH, height, SI, BAL 19,800 -2,798 551 9 DBH, height, SI, GEA, G 20,227 -1,961 546 10 DBH, height, SI, G 20,228 -1,956 573 11 DBH, height, SI, GEA, RDBH 21,517 -2,104 1,048 12 DBH, height, SI, RDBH 21,559 -2,101 1,088 a DBH is diameter at breast height, CurtisRD is Curtis’ Relative Density, G is basal area per hectare, RDBH is relative diameter at breast height, SPH is stems per hectare, CCFL is crown competition factor of larger trees, SI is Site Index, BAL is basal area of larger trees, and GEA is growth effective age 53 Table 2-4. AIC values averaged over the three model-building data subsets and concordance averaged over three test datasets by model and species. Model Subset AIC Values Percent Concordance Western hemlock Douglas- fir Western redcedar Western hemlock Douglas- fir Western redcedar Model 1 33,121 -8,675 -2,066 0.825 0.833 0.840 1 2 33,112 -8,580 -2,107 0.822 0.836 0.843 3 32,490 -8,980 -2,945 0.824 0.834 0.849 Average 32,908 -8,745 -2,373 0.824 0.835 0.844 Temesgen 1 44,027 -4,447 -81 0.823 0.831 0.822 (2002) 2 44,717 -3,909 613 0.822 0.836 0.826 3 44,007 -4,413 575 0.823 0.833 0.826 Average 44,250 -4,256 369 0.822 0.834 0.825 Hamilton 1 46,923 -4,293 802 0.821 0.831 0.811 (1986) 2 46,094 -3,722 1,417 0.820 0.835 0.820 3 46,237 -4,252 1,380 0.822 0.833 0.817 Average 46,418 -4,089 1,200 0.821 0.833 0.816 54 Table 2-5. Bias and root mean-squared error (RMSE in parentheses) for live stems per ha using the three implementation approaches for each species. Western hemlock Douglas- fir Western redcedar Probability Multiplier -7.1 (105.5) 7.4 (88.8) -3.4 (57.5) Stochastic (Number of trials) 1 -7.1 (106.7) 7.6 (141.4) -2.4 (57.1) 5 -6.8 (106.7) 7.3 (141.4) -3.6 (57.1) 20 -7.0 (105.4) 7.2 (140.1) -3.4 (57.0) 1,000 -7.2 (105.4) 7.2 (141.0) -3.5 (57.4) Threshold (Probabilities) 0.50 -69.1 (165.8) -37.7 (171.6) -23.2 (86.6) 0.75 11.3 (173.1) -0.2 (157.4) -4.8 (88.0) Optimum 17.7 (182.1) 12.3 (162.2) 2.9(101.5) 55 2.6 References Avery, T.E. and Burkhart, H.E. 2002. Forest Measurements, 5th Edition. McGraw-Hill Companies Inc. New York, N.Y. Avila, O.B. and Burkhart, H.E. 1992. Modeling survival of loblolly pine trees in thinned and unthinned plantations. Can. J. For. Res. 22(12): 1878-1882. Bravo-Oviedo, A., Sterba, H., del Rio, M., and Bravo, F. 2006. Competition-induced mortality for Mediterranean Pinus pinaster Ait. and P. sylvestris L. For. Ecol. Manage. 222(1-3): 88-98. Carter, R.E. and Klinka, K. 1992. Variation in shade tolerance of Douglas fir, western hemlock, and western red cedar in coastal British Columbia. For. Ecol. Manage. 55(1-4): 87-105. Cohen, J. 1960. A coefficient of agreement for nominal scales. Educ. Psychol. Meas. 20(1): 37-46. Curtis, R.O. 1982. A simple index of stand density for Douglas-fir. For. Sci. 28(1):92-94. Deubner, D.C., Wilkinson, W.E., Helms, M.J., Tyroler, H.A., and Hames, C.G. 1980. Logistic model estimation of death attributable to risk factors for cardiovascular disease in Evans County, Georgia. Am. J. Epidemiol. 112(1): 135-143. Dobbertin, M. and Biging, G.S. 1998. Using the non-parametric classifier CART to model forest tree mortality. For. Sci. 44(4): 507-516. Ek, A.R. 1980. A preliminary trial of alternative methods for treating mortality in multipurpose forest projection system (MFPS) model. 1980. Univ. Minn., Coll. For., Dep. For. Resour., Staff Pap. Ser. No. 8. 56 Fan, Z., Kabrick, J.M., and Shifley, S.R. 2006. Classification and regression tree based survival analysis in oak-dominated forests of Missouri’s Ozark highlands. Can. J. For. Res. 36(7): 1740-1748. Fleiss, J.L. 1975. Measuring agreement between two judges on the presence or absence of a trait. Biometrics 31(3): 651-659. Franklin, J.F., Shugart, H.H., and Harmon, M.E. 1987. Tree death as an ecological process. BioScience 37(8): 550-556. Gavin, D.G., Brubaker, L., and Lertzman, K.P. 2003. Holocene fire history of a coastal temperate rain forest based on soil charcoal radiocarbon dates. Ecology 84(1) 186- 204. Getzin, S., Dean, C., He, F., Trofymow, J.A., Wiegand, K. and Wiegand, T. 2006. Spatial patterns and competition of tree species in a Douglas-fir chronosequence on Vancouver Island. Ecography 29(5): 671-682. Hamilton, D.A., Jr. 1986. A logistic model of mortality in thinned and unthinned mixed conifer stands of northern Idaho. For. Sci. 32(4): 989-1000. Hamilton, D.A., Jr. and Edwards, B.M. 1976. Modeling the probability of individual tree mortality. USDA For. Serv. Res. Pap. INT-185. 22 pp. Hann, D.W. and Ritchie, M.W. 1988. Height growth rate of Douglas fir: a comparison of model forms. For. Sci. 34(1): 165-175. Kobe, R.K. and Coates, K.D. 1997. Models of sapling mortality as a function of growth to characterize interspecific variation in shade tolerance of eight tree species of northwestern British Columbia. Can. J. For. Res. 27(2): 227-236. 57 Kotar, J. 1972. Ecology of Abies amabilis in relation to its altitudinal distribution in contrast to its common associate Tsuga heterophylla. PhD Thesis. University of Washington, Seattle, WA. Lee, Y. 1971. Predicting mortality for even-aged stands of lodgepole pine. For. Chron. 47(1): 29-32. Lertzman, K.P., Sutherland, G.D., Inselberg, A. and Saunders, S.C. 1996. Canopy gaps and the landscape mosaic in a coastal temperate rain forest. Ecology 77(4): 1254- 1270. Meidinger, D. and Pojar, J., eds. 1991. Ecosystems of British Columbia. The Research Branch, BC Ministry of Forests, Victoria, BC. pp. 95-111. Monserud, R.A. 1976. Simulation of forest tree mortality. For. Sci. 22(4): 438-444. Oliver, C.D. 1981. Forest development in North America following major disturbances. For. Ecol. Manage. 3: 153-168. SAS Institute Inc. 2004. SAS/STAT 9.1 user’s guide. SAS Institute Inc., Cary, N.C. Smith, J.H.G., 1966. Studies of crown development are improving Canadian forest management. In: Proc. VI World Forestry Congress, Madrid. Vol. 2, pp. 2309–2315. Spies, T.A., Franklin, J.F., and Klopsch, M. 1990. Canopy gaps in Douglas-fir forests of the Cascade Mountains. Can. J. For. Res. 20(5): 649-658. Temesgen, H. 2002. Individual-tree mortality models for the interior cedar-hemlock and interior Douglas-fir biogeoclimatic ecosystem classification zones. Internal report for J.S. Thrower and Associates, Vancouver, BC, pp 1-17. 58 Weber, L.A., Ek, A.R., and Droessler, T.D. 1986. Comparison of stochastic and deterministic mortality estimation in an individual tree based stand growth model. Can. J. For. Res. 16(5): 1139-1141. Yang, Y., Titus, S., and Haung, S. 2003. Modeling individual tree mortality for white spruce in Alberta. Ecol. Model. 163(3): 209-222. Zhao, D., Borders, B., and Wilson, M. 2004. Individual-tree diameter growth and mortality models for bottomland mixed-species hardwood stands in the lower Mississippi alluvial valley. For. Ecol. Manage. 199(2-3): 307-322. Zhao D., Borders, B., Wang, M., and Kane, M. 2007. Modeling mortality of second- rotation loblolly pine plantations in the Piedmont/Upper Coastal Plain and Lower Coastal Plain of the southern United States. For. Ecol. Manage. 252(1-3): 132-143. 59 3. Modeling the effects of thinning and fertilization on mortality of coastal Douglas-fir (Pseudostuga menziesii var. menziesii (Mirb.) Franco) 2 3.1 Introduction The mixed-species coastal temperate rain forest of coastal British Columbia (BC) is noted for its lack of recent fires (Lertzman et al. 1996). The successional processes of the Pacific Northwest forests are defined by small-scale disturbances, mainly due to windthrow (Gavin et al. 2003). Hence, gap dynamics define the forest structure and natural successional processes on Vancouver Island (Lertzman et al. 1996), with Douglas-fir (Pseudotsuga menziesii var. menziesii (Mirb.) Franco) being the dominate, pioneer species and western hemlock (Tsuga heterophylla (Raf.) Sarg.) and western redcedar (Thuja plicata Donn) being the late-successional species aggregated within gaps (Getzin et al. 2006). While mortality in younger stands is typically due to self-thinning processes, as a stand matures, gaps formed from single-tree mortality drive the structural and compositional characteristics instrumental in stand development (Lertzman et al. 1996). A survival model, or conversely a mortality model, is an important component within growth and yield models (Avila and Burkhart 1992). Tree mortality is a complex process resulting from interactions among biotic and abiotic factors and many of the ecological processes and systems which contribute to tree mortality are not easily understood (Franklin et al. 1987). Factors affecting tree mortality include: microsite 2 A version of this chapter will be submitted for publication. Rathbun, L.C., and LeMay, V. 2010. Modeling the effects of thinning and fertilization on mortality of coastal Douglas-fir (Pseudostuga menziesii var. menziesii (Mirb.) Franco). 60 nutrient and climatic conditions, stand-level site and climatic conditions, inter-tree competition, age, incidence of insect and/or disease attack, and fire. Mortality has classically been divided into catastrophic mortality, involving the death of many trees at one time (e.g., due to fire, epidemic insect/disease attacks, wind damage, etc.), and regular mortality, where single or few trees die at one time (Lee 1971). Regular mortality can occur as a result of competition with other trees and plants (i.e., competition-based), as a result of endemic insects or pathogens, or as a result of changes in microsite conditions. Commonly, the mortality component of a growth and yield model considers only regular mortality; however, it should be noted that regular and catastrophic mortality may not always be independent (Dobbertin and Biging 1998). Stand management can affect the rates of regular mortality through the alteration of both competition and site characteristics. Thinning results in an immediate decrease of inter-tree competition from the removal of nearby trees, thereby improving individual tree vigor and hence reducing the rate of mortality for the residual trees (e.g., Curtis et al. 1981, Hann et al. 2003). Bailey and Tappeiner (1998) concluded that understorey trees in thinned coastal Douglas-fir stands, had significantly longer live crown ratios compared to unthinned stands and O’Hara (1988) found that Douglas-fir trees with a medium-sized crown in thinned stands utilized their growing space as efficiently as trees with a large- sized crown in unthinned stands, both indications of increased tree vigor. Recently, the impacts of thinning on competition-based mortality have been incorporated into mortality models. For a single thinning event, thinning impacts have been incorporated through the development of two models, one for unthinned stands and one for thinned stands (Avila 61 and Burkhart 1992), or by the inclusion of a thinning indicator variable in the model (e.g., Hann et al. 2003, Karlsson and Norell 2005). In the past decade, management practices in coastal BC have moved from clearcutting toward variable-retention systems which provide varying intensities and sizes of disturbance across a gradient (Beese et al. 2003). The methods commonly used for incorporating thinning into mortality models do not account for multiple thinning events, high variation among thinning intensities, or differences in stage of stand development when thinning is applied, as is so often found with variable retention systems. To do so, mortality models developed with the inclusion of additional stand and/or tree-level variables which indicate changes in competition as a result of thinning are needed (e.g., Hamilton 1986). Fertilization is a common practice in this area as well, with the growth of Douglas-fir responding well to nitrogen fertilization (e.g., Brix 1981, Weetman et al. 1997). Shen et al. (2001) found that increasing rates of N application increased the probability of mortality for individual trees for Douglas-fir stands located in the inland Northwest; in coastal BC, fertilization has been shown to accelerate mortality for Douglas-fir (Brix 1993). In this study, impacts of thinning and fertilization on competition-based mortality (hereafter termed “mortality”) of coastal Douglas-fir within the mixed-conifer forests of coastal BC were studied. The first objective of this study was to modify an existing survival model (Rathbun et al. 2010) developed for untreated second-growth Douglas-fir trees to thinned trees from the same study site in order to model the impacts of thinning. Variations in thinning applications included the number of thinning events (including no 62 thinning), time of thinning, and thinning intensity. The hypotheses associated with this first objective were: 1) the time since thinning will affect the reduction in mortality, which will increase to a maximum following thinning and then decline; 2) given a single thinning event, mortality declines with thinning intensity, as a consequence of reduced competition; 3) a greater mortality response (decrease in the percent mortality) will occur when thinning occurs at an earlier time in stand development; and 4) repeated thinning events do not necessarily result in repeated reductions in mortality. The second objective was to modify the same survival model to include the effects of fertilization. The hypotheses were: 1) given a single fertilization application, a greater mortality response will occur compared to untreated stands; 2) fertilization applied to higher quality sites will increase mortality; and 3) impacts of thinning and fertilization on mortality can be represented using a two-step additive approach, in that mortality was modified for thinning, and then again for fertilization. These hypotheses were examined using a modeling approach. Initially, the existing survival model for unthinned plots was fitted. The model was applied directly to the thinned stands. An additional variable to define the effects or fertilization was added to the model and applied to the fertilized stands and stands which were both thinned and fertilized using measures both prior to and after thinning to forecast the expected mortality. The forecast was compared to actual mortality, and differences were used to assess variations in treatment events on mortality. 63 3.2 Materials and methods 3.2.1 Study area The study area was located on Vancouver Island, Haida Gwaii (the Queen Charlotte Islands), and along the coastal mainland within the Coastal Western Hemlock (CWH) and Coastal Douglas-fir (CDF) Biogeoclimatic Ecological Classification (BEC) zones of BC (Meidinger and Pojar 1991). The CWH BEC zone is divided into multiple subzones: a very dry maritime subzone in the east, a moist maritime subzone in the central area, and a very wet maritime and hypermaritime subzone in the west (Meidinger and Pojar 1991). The temperatures in the CWH BEC zone range from 5.2 to 10.5˚ C, with a mean annual temperature of 8˚ C (Meidinger and Pojar 1991). The CDF BEC zone includes the rainshadow of Vancouver Island and the Olympic mountains (Meidinger and Pojar 1991). It is defined by warm, dry summers and mild, wet winters with mean annual temperatures ranging from 9.2 to 10.5˚ C, and a minimum temperature ranging from -21.1 to -11.7˚ C (Meidinger and Pojar 1991). Study plots were located within latitudes ranging from 48.37 to 53.53˚ N and longitudes ranging from 122.36 to 132.59˚ W. The majority of this area is comprised of second growth multi-species stands, regenerated both naturally and from planting. Common tree species include: Douglas-fir, western redcedar, western hemlock, red alder (Alnus rubra Bong.), Sitka spruce (Picea sitchensis (Bong.) Carr.), and yellow cedar (Chamaecyparis nootkatensis (D.Don) Spach). 3.2.2 Data Permanent sample plot (PSP) data were provided by Island Timberlands, LP. The database for Douglas-fir as the primary species contained 867 untreated plots, ranging in 64 size from 0.008 to 0.253 ha. Measurements were collected between 1932 to 2002, with measurement intervals varying from 1 to 17 years with an average of 4.9 years. The majority of plots, 830, were defined as second growth, less than 150 years since stand establishment, 18 plots were defined as old growth, greater than 150 years since stand establishment, and 19 plots had no available age information. At the time of plot establishment, 515 plots were predominantly Douglas-fir, 40 plots were predominantly western redcedar, and 202 plots were predominately western hemlock, as defined by greater than 50% stems per ha. Based on more than 50% of basal area per ha, 580 plots contained predominantly Douglas-fir, 29 plots predominantly western redcedar, and 168 plots predominately western hemlock. For the remaining plots, species composition averaged 16, 7, and 17% Douglas-fir, western redcedar, and western hemlock, respectively, based on basal area per ha. Densities at plot establishment ranged from 222 to 8,860 live trees per ha and site index, using coastal Douglas-fir site index at a base age of 50 years measured at breast height (1.3 m above ground), ranged from 8.7 to 49.3 m (Table 3-1). Average annual plot mortality for the first measurement period following plot establishment ranged from 0 to 10.29% with an average of 1.20%. The average annual mortality for the second measurement period was 1.16%, with a slightly wider range of 0 to 12.07%. At each measurement, the species, tree status (i.e., live or dead), and diameter outside bark at 1.3 m above ground (DBH, cm) were recorded for all trees with DBH greater than 4 cm and within plot boundaries. Height (ht, m) was measured for a subset of trees, and remaining heights were estimated using height-diameter functions. At the time of plot establishment, the average DBH was 12.5 cm, with corresponding minimum 65 and maximum values of 4.0 and 149.9 cm (Table 3-2). The average height was 11.7 m, with corresponding minimum and maximum values of 0 and 55.8 m. Since the heights were a mixture of measured and estimated heights, the minimum height was 0 rather than the expected value of ≥ 1.3m. For model fitting, validation, and testing, all non- overlapping repeated measures for all plots were used. An additional 395 plots received silvicultural treatments; treatments included fertilization (44 plots), thinnings (294 plots), and a combination of fertilization and thinning (57 plots). The majority of fertilized plots received one application of nitrogen ranging in concentration from 50 to 400 kg per ha and a few plots additionally received ammonium phosphate. Of the 294 plots which were thinned, 218 received one thinning and 76 received multiple thinnings. Basal area removed varied from 0.2 to 45.2 m 2 per ha and stems per ha removed varied from 25 to 9,012 trees per ha, with corresponding averages of 7.6 m 2 per ha and 654 trees per ha. For the plots which received a combination of thinnings and fertilization, fertilization was applied within the measurement period immediately following the thinning treatment. It should be noted that the data used in this study primarily represent silvicultural treatments applied to managed stands rather than treatments applied in a experimental setting. Therefore, the results in this study include the interactions that are often removed from experiments through careful selection of experimental plots. 3.2.3 Survival model The survival model used for the untreated plots was developed and reported in Rathbun et al. (2010) for coastal Douglas-fir. The model contains predictor variables representing tree size and stage of development, site productivity, and inter-tree 66 competition. Tree size and stage of development were represented by DBH (cm) and height (m). To represent site productivity, coastal Douglas-fir site index (SI, m) at a base age or 50 years was included. Also, Growth Effective Age (GEA, years) was used to indicate tree productivity. GEA is the age of a dominant tree of the same height and as the tree of interest, for a given site index (Hann and Ritchie 1988). For these data, GEA was estimated from site index equations for Douglas-fir. Two stand level competition measures were included. Stand basal area has been commonly used to measure competition for below ground resources (Fan et al. 2006). Curtis’ (1982) Relative Density, an index measure of density, calculated as: qd G RD Curtis n DBH d n 1i 2 i q where G is the stand basal area per ha (m 2 /ha) ; DBHi (cm) is the DBH for tree i, dq (cm) is the quadratic mean diameter for the plot, n is the number of live trees across all species within a plot and all other variables are as previously defined. In addition, a tree-level competition measure, basal area of larger trees (BAL, m 2 /ha) was included in the model. BAL was calculated as: n i iBAtreeδ 1 x BAL where δ is an indicator variable (1 if tree i has a DBH greater than the tree of interest, 0 otherwise), BAtreei (m 2 /ha) is the basal area per hectare value for tree i, and n is the number of all live trees across all species within a plot. 67 The selected survival model (Rathbun et al. 2010) was: [1] q Xf-p(s) exp1 GaGEAaSIaCurtisRDaBALaheightaDBHaaf(X) 76543210 where p(s) is the probability of survival at the end of the projection period, f(X) is a linear function of the predictor variables at the beginning of the projection period, q is the number of years within a projection period, a0 through a7 are parameters to be estimated, and all other variables are as previously defined. Parameter estimates for the selected survival model (Eq. [1]) are given in Table 3-3. 3.2.4 Treatment effects 3.2.4.1 Thinning Thinning has often been modeled through the inclusion of thinning parameters designed to describe the type of treatment received. For example, Ritchie et al. (2007) included interactions between categorical variables describing treatment types and the original state variables in the model. In this study, no modifications were made to the survival selected for the untreated plots for thinning. Instead, the selected model given as Eq. [1] was used directly, such that the state variables changed instantaneously after thinning as input variables to predict mortality following thinning. In particular, state variables which changed as a result of the thinning were basal area per ha, basal area of larger trees per ha, and CurtisRD. This approach was used for Douglas-fir stands located in Idaho, USA by Hamilton (1986). To assess the thinning hypotheses, the average predicted annual diameter increments were graphed versus time since treatment using trees on which measurements 68 prior to treatment were taken, and repeated for trees which did not have measures prior to treatment. This was done using four different methods: 1) across all plots, 2) by thinning intensity, 3) by stage of stand development when thinning occurred, 4) and by the number of thinnings. Thinning intensity was defined as low (less than 20% basal area removed), medium (between 20 and 40% basal area removed), and high (greater than 40% basal area removed). Stage of stand development was defined at the time of thinning as stand initiation (stands less than 20 years basal height age), stem exclusion (stands between 20 and 100 years breast height age), and understorey re-initiation (stands between 100 and 250 years breast height age) (Franklin et al. 2002). Old growth was defined as those stands greater than 250 years for breast height age (Trofymow et al. 1997), but so few measurements were available for old growth stands that they were not included in graphing. 3.2.4.2 Fertilization The selected survival model was modified to include the effects of fertilization. An additional fertilization effect variable, defined as: Ftime F effect.fert where Ftime is the time (years) since fertilization, and F is an indicator variable, which was 1 for fertilized plots and 0 otherwise, was added to f(X) to obtain f(Z). This approach was used by Hann et al. (2003) in modeling 5-year mortality rates of Douglas- fir and western hemlock in the coastal region of Oregon and Washington State. The model for the fertilized plots was: [2] q Zf-p(s) exp1 69 F 876543210 time F aGaGEAaSIaCurtisRDaBALaheightaDBHaaf(Z) where a8 is a parameter to be estimated and all other variables are as previously defined. To assess the fertilization hypotheses, the average predicted annual diameter increments were graphed versus time since treatment using trees on which measurements prior to treatment were taken, and repeated for trees which did not have measurements prior to treatment. This was done using four different methods: 1) across all plots, and 2) by site quality. Site quality was defined using the BC Ministry of Forests and Range (1999) site index classifications as low (site index values less than 22 m), medium (site index values between 22 and 30 m), and high (site index values greater than 30 m). 3.2.4.3 Thinning and fertilization For plots that received a combination of thinning and fertilizer treatments, the selected survival model [2] modified for fertilized plots was applied directly. As for the plots which only received a thinning, some state variables for these plots immediately changed following thinning. This approach assumed a greater response would be seen for plots which received both thinning and fertilization than for those which received only one treatment. To assess the fertilization and thinning in combination hypothesis, the average predicted annual percent mortality was graphed versus time since treatment using trees on which measurements prior to treatment were taken, and repeated for trees which did not have measurements prior to treatment. 70 3.2.5 Model accuracy Tree and plot measurements at plot establishment (i.e., the beginning of the first measurement period) were used as inputs to the selected model to predict the probability of survival to the end of the first measurement period (i.e., beginning of the second measurement period). Once the probability of survival was estimated, the model was implemented using the probability multiplier approach (Rathbun et al. 2010), where the probability of survival was used to estimate the surviving stems per ha at the end of the period by: SPHq = p(s)* SPH0 where SPHq is the stems per ha surviving to the end of the measurement period represented by the tree, p(s) is the predicted probability of survival for the tree, and SPH0 is the stems per ha represented by the tree at the beginning of the measurement period. To assess accuracy for all models by treatment type, mean bias and root mean- squared error (RMSE) values by 5 cm diameter classes were calculated for measures after treatment occurred. Mean bias was calculated as: m yy m 1j jj ˆ BiasMean where jy is the actual number of trees per ha that survived within the q period for plot j, jŷ is the predicted number of trees per ha that survived within the q period for plot j, and m is the number of plots. Root mean-squared error was calculated as: 1pm yy m 1j 2 jj ˆ RMSE 71 where p is the number of predictor variables in the model and all other variables are as previously defined. 3.3 Results 3.3.1 Treatment effects 3.3.1.1 Thinning The differences between average actual and predicted annual percent mortality found in Figures 3-1a and b are quite small, especially for Figure 3-1b, plots where measurement were not available prior to thinning. Figure 3-1a shows a slight decrease in average annual mortality following thinning which levels off around 10 years post thinning. Figure 3-1b also shows a slight decrease in average annual mortality for up to 10 years following thinning, after which it increases slightly. For plots with measurements prior to thinning, the high intensity thinned plots had the smallest average predicted annual percent mortality following a thinning (Figure 3-2a), with the low and medium intensity thinning maintaining similar average annual percent mortality. Similar results were obtained for plots without measurements prior to thinning up to 15 years following thinning. The highest average annual percent mortality was found 10 to 15 years following thinning, for plots which were thinned during stand establishment (Figures 3-3a and b). It should be noted that there were only two measurements available for this time period. When plots were thinned during the stem exclusion phase of stand development, the highest average annual percent mortality was found (Figure 3-3a and b). Average annual percent mortality following thinning was the lowest for stands which were in the understorey re-initiation phase of stand development (Figure 3-3a and b); 72 both figures showed a constant rate in actual average annual percent mortality following thinning. The highest average annual plot mortality was found for plots following a second thinning (Figure 3-4). The lowest mortality was found for plots following a third thinning. 3.3.1.2 Fertilization Both average actual and predicted annual percent mortality increased immediately following fertilization on unthinned plots (Figure 3-5a and b). For plots with measurements available prior to fertilization, average annual percent morality increased from just less than 2% to 3.5%. An increase in average annual percent mortality was seen immediately following fertilization for the medium and high productivity sites (Figure 3-6a). For low site index stands, measurements were only available following the application of fertilizer; these stands had the smallest average predicated and actual annual percent mortality immediately following fertilization (Figure 3-6b). 3.3.1.3 Thinning and fertilization Thinning combined with fertilization resulted in a large increase in average actual annual percent mortality which is not modeled well (Figures 3-7a and b). This increase was only seen for the first five-year measurement period immediately following treatment. 3.3.2 Model accuracy All mean bias values were negative except for the smallest diameter trees, less than 5cm (Figure 3-8a). The fertilized plots had the largest mean bias value for the smallest diameter trees, and the mean bias was similar across treatments for trees greater than 30 cm DBH. No matter the treatment type, RMSE values decreased with increasing 73 diameter class (Figure 3-8b). Fertilized trees had the largest RMSE values for trees less than 25cm DBH; while trees which were both thinned and fertilized had the smallest RMSE values for the same size trees. 3.4 Discussion The mortality model incorporated competition-based mortality through the inclusion of additional stand (basal area per hectare and Curtis’ relative density) and tree- level (basal area of larger trees) competition variables (Rathbun et al. 2010). The inclusion of these competition measures in the model capture the change in annual percent mortality resulting from thinning events (Figure 3-1a and b). As was hypothesized, the model predicted a decrease in mortality following thinning for the first 10 years, followed by a leveling off period to year 15, with a subsequent slight increase in mortality following year 15 (Figure 3-1b only). The first 5-year period immediately following thinning produced the greatest post-harvest mortality. During thinning, damage to residual trees can occur. Thinning shock is related to thinning intensity, site quality, and species (Harrington and Reukema 1983). Kneeshaw et al. (2002) found that stem growth did not increase significantly until two years following treatment, and that larger trees experienced a greater growth shock in the first year following treatment for Douglas-fir located on Vancouver Island. Both of these factors may explain the increase in mortality found within the first five years after treatment in this study. Management practices in this area have moved from clearcutting toward variable- retention harvesting (Beese et al. 2003). This shift further complicates stand structure by diversity within treatments through variation in thinning intensity, timing of thinnings, and number of thinning, as represented in this study. Variable retention harvesting in 74 these areas has been shown to increase mortality of overstorey trees, mainly due to increased damage from windthrow and weather on residual trees (Maguire et al. 2006). Following a single thinning event and as a consequence of reduced competition, mortality did not necessarily decline with thinning intensity as was hypothesized. The level of removal was an important determinant of the level of mortality. Stands where greater than 40% of the basal area was removed had much lower mortality than those where less than 40% was removed (Figure 3-2a and b). Zenner (1995) found that low intensity thinning caused short-term increases in understorey light, but did not create sustained increases in belowground resources for Douglas-fir located in Oregon, USA. To do so, controlling the amount of understorey vegetation may be needed (Harrington 2006). Initially after thinning, stands where 20-40% of basal area was removed had higher mortality than those where less than 20% was removed; this changed after 10 years following treatment (Figure 3-2a). This may be due in part to thinning shock. Plots which were thinned within the first 20 years following stand establishment showed no change in mortality following thinning (Figures 3-3a and b). During natural regeneration, when stands are establishing, canopy closure may be gradual and the period of density-dependent mortality may postponed (Franklin et al. 2002). A sustained decrease was seen for plots which were between 20 to 100 years since plot establishment (stem exclusion phase) when thinning occurred. For mixed-conifer stands located in Oregon and Washington States, USA, Aubry et al. (2009) found that mortality was more frequent for suppressed trees following thinning. Thinning during the stem exclusion phase reduces the density dependent mortality which often occurs (Franklin et al. 2002). 75 Stands harvested during the understorey re-initiation phase, between 100 to 250 years since stand establishment, showed the greatest decrease in mortality within the first five to ten years following thinning, after which mortality increased. This was contrary to the hypothesis that thinning at an earlier time in stand development would decrease mortality. This phase of stand development sees a shift from density dependent to density independent mortality for overstorey trees (Franklin et al. 2002). As hypothesized, repeated thinning events did not result in repeated reductions in mortality. Plots which were thinned twice actually had higher mortality following the second thinning than those which were thinned only once (Figure 3-4), while plots which received a third thinning had the lowest mortality. This trend was maintained through time. While the second thinning may have increased the level of resources available, trees may not have been able to utilize them if they were not available or shared. Competitors may have been more efficient at utilizing the newly available resources (Cole and Newton 1986). Damage to residual trees due to a second thinning may also have been a factor as well as the time between thinning events. As was hypothesized, an increase in mortality occurred following fertilization (Figure 3-5a). Mortality almost doubled from pre-fertilization to post-fertilization. Shen et al. (2001) found that heavier nitrogen applications increased the probability of mortality for interior Douglas-fir (Pseudotsuga menziesii var. glauca (Mirb.) Franco). This could be due in part to increasing individual tree growth which in turn increases competition, or due to increased ingrowth which also increases competition. Miller et al. (1986) found that fertilization accelerated mortality in smaller than average trees as a result of increased basal area growth of larger trees in unthinned fertilized stands of 76 Douglas-fir. Mortality was found to decrease slightly after the first five years following treatment. Barclay and Layton (1990) found that much of the mortality in unthinned fertilized plots occurred soon after treatment. As hypothesized, low quality sites had the lowest mortality of fertilized stands and high quality sites had the largest mortality (Figures 3-6a and b). Mortality more than doubled for medium and high quality sites immediately following treatment, with little differences in mortality seen between the two. For stands with the same density of tree regeneration, canopy closure will be more rapid on the higher quality site and therefore increase mortality (Franklin et al. 2002). Miller et al. (1986) found that Douglas-fir stands on poor-quality sites showed increased growth compared to stands on medium and good-quality sites. The increased tree growth, without a reduction in density, increases competition for resources and hence increases mortality. Mortality, while still remaining higher than pretreatment values, decreased after five years following treatment. Fertilization response has been shown to last for 6-8 years (Chappell et al. 1992). It should be noted that Brix (1993) found repeated applications of fertilization accelerated mortality greatly in unthinned plots of Douglas-fir on Vancouver Island. A large difference can be seen between the actual and predicted values of mortality which occur within the first five years following treatment for plots which were both thinned and fertilized (Figures 3-7a and b). Past five years following treatment, the actual and predicted values are similar. This indicates that the impacts of thinning and fertilization on mortality can be represented using a two-step additive approach as the hypothesis indicated, but not for the first five years following treatment where thinning 77 stress may occur. It should be noted that the reported interactions between fertilizer and thinning processes are inconsistent within the literature (Li et al. 2007). As noted previously, the data used in this study primarily represent silvicultural treatments applied to managed stands rather than treatments applied in a experimental setting; the intent was to obtain a mortality model as a component of a growth model for use on these stands. As a result, the impacts of fertilization and thinning are more variable than in studies based on experiments. 3.5 Conclusion The mortality model previously developed for untreated stands of Douglas-fir in coastal BC was modified to reflect the effects of thinning and fertilization. The thinned data included variation due to thinning intensity, timing of thinnings, and number of thinnings, variations often found in variable retention systems. The inclusion of competition variables such as stand basal area per hectare, basal area per hectare of larger trees, and Curtis’ relative density modeled the effects due to thinning well. The additional variable added to the model for fertilized stands modeled the increase in mortality seen after application. If the desire is to reduce mortality through thinning, the model suggests that thinning be done at a high intensity and prior to understorey re-initiation. If stands are fertilized without thinning, low quality sites have less mortality after thinning than high quality sites. Additional research is needed for stands which are both thinned and fertilized and for stands where variable retention systems have occurred. 78 Figure 3-1. Average actual (♦) and predicted (●) annual percent mortality for thinned data by time since thinning for a) plots with measurements prior to thinning and for b) plots without measurements prior to thinning. The time of thinning is given as 0, and the number of plot measures used to determine average values are given across the top of each graph. 79 Figure 3-2. Average predicted annual percent mortality by thinning intensity for thinned data by time since thinning for a) plots with measurements prior to thinning and for b) plots without measurements prior to thinning. The time of thinning is given as 0, and thinning intensity is identified as low (●), medium (♦), and high (■). 80 Figure 3-3. Average predicted annual percent mortality by stand development for thinned data by time since thinning for a) plots with measurements prior to thinning and for b) plots without measurements prior to thinning. The time of thinning is given as 0, and stand development is identified as stand establishment (●), stem exclusion (♦), and understorey re-initiation (■). 81 Figure 3-4. Average predicted annual percent mortality by number of thinnings for thinned data by time since thinning. The time of thinning is given as 0, and the number of thinnings is identified as one (●), two (♦), or three (■). 82 Figure 3-5. Average actual (♦) and predicted (●) annual percent mortality for fertilized data by time since fertilization for a) plots with measurements prior to fertilization and for b) plots without measurements prior to fertilization. The time of fertilization is given as 0, and the number of plot measures used to determine the average values are given across the top of each graph. 83 Figure 3-6. Average predicted annual percent mortality by site index class for fertilized data by time since fertilization for a) plots with measurements prior to fertilization and for b) plots without measurements prior to fertilization. The time of fertilization is given as 0, and site index class is identified as low (●), medium (♦), and high (■). 84 Figure 3-7. Average actual (♦) and predicted (●) annual percent mortality fertilized and thinned stands by time since treatment for a) plots with measurements prior to treatment and for b) plots without measurements prior to treatment. The time of treatment is given as 0, and the number of plot measures used to determine the average values are given across the top of each graph. 85 Figure 3-8. Bias a) and b) RMSE (cm) values by 5 cm DBH classes (defined at the midpoint) by treatment. Treatment was defined as untreated (●), thinned (■), fertilized (♦), and both thinning and fertilization (+). 86 Table 3-1. Plot summary statistics at plot establishment and at time of treatment. (n= number of plots) Treatment Type Variable a Plot Establishment Time of Treatment Mean Min b Max c Mean Min b Max c Untreated n=867 SPH 1993 222 8860 G (m 2 per ha) 21.1 0.1 185.0 SI (m) 29.2 8.7 49.3 CurtisRD 7.6 0.8 24.2 Fertilized n=44 SPH 2,559 765 7,383 2,448 765 6,840 G (m 2 per ha) 33.6 1.4 80.6 37.9 1.4 80.6 SI (m) 30.5 11.2 42.4 30.7 11.2 42.3 CurtisRD 8.6 0.7 14.5 9.3 0.7 15.1 Thinned n=294 SPH 1,853 370.8 9,109 1.847 371 9,235 G (m 2 per ha) 32.0 0.5 110.5 33.7 1.6 110.5 SI (m) 30.7 9.2 50.0 30.9 9.2 50.0 CurtisRD 7.6 0.3 18.4 7.8 0.7 18.4 Multiple Treatments n=57 SPH 2,401 790 7,333 2,361 790 6,815 G (m 2 per ha) 43.5 5.0 76.2 45.1 5.0 76.2 SI (m) 31.6 8,9 41.7 31.6 8.9 41.7 CurtisRD 10.6 2.2 16.1 10.8 2.2 16.1 * SPH is stems per ha, G is basal area per ha, SI is Site Index, and CurtisRD is Curtis’ Relative Density b Min is Minimum c Max is Maximum. 87 Table 3-2. Means and standard deviations (in brackets) for tree-level variables at plot establishment by silvicultural treatment. Variable a Untreated Fertilized Thinned Multiple Treatments Number of trees 33,960 5,040 14,903 4,492 DBH (cm) 12.5 (10.1) 11.3 (7.6) 16.2 (10.4) 15.6 (10.0) Height (m) 11.7 (7.5) 11.4 (6.2) 15.5 (8.3) 15.2 (7.2) BAL (m 2 per ha) 18.4 (15.4) 22.0 (14.6) 20.6 (14.6) 24.9 (13.4) RDBH 1.1 (0.5) 1.0 (0.4) 1.1 (0.4) 1.2 (0.5) GEA (years) 17.4 (16.5) 15.3 (9.5) 21.8 (16.2) 18.9 (12.4) Dinc (cm) 0.42 (0.42) 0.34 (0.30) 0.46 (0.50) 0.36 (0.38) a DBH is diameter at breast height, BAL is basal area per ha of larger trees, RDBH is relative diameter at breast height, GEA is growth effective age, and Dinc is the annual diameter increment. 88 Table 3-3. Parameter estimates for the selected survival model (Eq. 1) for untreated plots (Rathbun et al. 2010). Parameter Estimates a0 7.067 a1 0.270 a2 -0.255 a3 -0.220 a4 -0.202 a5 -0.055 a6 0.030 a7 0.218 89 3.6 References Aubry, K.B., Halpern, C.B., and Peterson, C.E. 2009. Variable-retention harvests in the Pacific Northwest: a review if short-term findings from the DEMO study. For. Ecol. Manage. 258(4): 398-408. Avila, O.B. and Burkhart, H.E. 1992. Modeling survival of loblolly pine trees in thinned and unthinned plantations. Can. J. For. Res. 22(12): 1878-1882. Bailey, J.D. and Tappeiner, J.C. 1998. Effects of thinning on structural development in 40- to 100-year-old Douglas-fir stands in western Oregon. For. Ecol. Manage. 108(1- 2): 99-113. Barclay, H.J., and Layton, C.R. 1990. Growth and mortality in managed Douglas fir: Relation to a competition index. For. Ecol. Manage. 36(2-4): 187-204. B.C. Ministry of Forests and Range. 1999. Timber supply review: Strathcona timber supply area analysis report. Province of British Columbia Ministry of Forests. Victoria, BC. 91 pp. Beese, W.J., Dunsworth, B.G., Zielke, K., and Bancroft, B. 2003. Maintaining attributes of old-growth forests in coastal B.C. through variable retention. For. Chron. 79(3) 570-578. Brix, H. 1981. Effects of thinning and nitrogen fertilization on branch and foliage production in Douglas-fir. Can. J. For. Res. 11(3): 502-511. Brix, H. 1993. Fertilization and thinning effect on a Douglas-fir ecosystem at Shawnigan Lake: A synthesis of project results. For. Can. and B.C. Min. For., Victoria, BC FRDA Rep. No 196. 90 Chappell, H.N., Omule, S.A.Y., and Gessel, S.P. 1992. Fertilization in coastal northwest forests: Using response information in developing stand-level tactics. In: H.N. Chappell, G.F. Weetman and R.E. Miller (Ed.), Forest Fertilization: Sustaining and improving nutrition and growth of western forests. Institute of Forest Resources, Contrib. 73, College of Forest Resources, University of Washington, Seattle, WA, pp. 98-113. Cole, E.C., and Newton, M. 1986. Fifth-year responses of Douglas-fir to crowding and nonconiferous competition. Can. J. For. Res. 17(2): 181-186. Curtis, R.O., Clendenen, G.W., and DeMars, D.J. 1981. A new stand simulator for Douglas fir-DFSIM user’s guide. USDA For. Serv. Gen. Tech. Rep. PNW-GTR-128. Curtis, R.O. 1982. A simple index of stand density for Douglas-fir. For. Sci. 28(1):92-94. Dobbertin, M. and Biging, G.S. 1998. Using the non-parametric classifier CART to model forest tree mortality. For. Sci. 44(4): 507-516. Fan, Z., Kabrick, J.M., and Shifley, S.R. 2006. Classification and regression tree based survival analysis in oak-dominated forests of Missouri’s Ozark highlands. Can. J. For. Res. 36(7): 1740-1748. Franklin, J.F., Shugart, H. H., and Harmon, M.E. 1987. Tree death as an ecological process. BioScience 37(8): 550-556. Franklin, J.F., Spies, T.A., Van Pelt, R., Carey, A.B., Thornburgh, D.A., Berg, D.R., Lindenmayer, D.B., Harmon, M.E., Keeton, W.S., Shaw, D.C., Bible, K., and Chen, J. 2002. Disturbances and structural development of natural forest ecosystems with silvicultural implications, using Douglas-fir as an example. For. Ecol. Manage. 155(1-3): 399-423. 91 Gavin, D.G., Brubaker, L., and Lertzman, K.P. 2003. Holocene fire history of a coastal temperate rain forest based on soil charcoal radiocarbon dates. Ecology 84(1) 186- 204. Getzin, S., Dean, C., He, F., Trofymow, J.A., Wiegand, K. and Wiegand, T. 2006. Spatial patterns and competition of tree species in a Douglas-fir chronosequence on Vancouver Island. Ecography 29(5): 671-682. Hamilton, D.A., Jr. 1986. A logistic model of mortality in thinned and unthinned mixed conifer stands of northern Idaho. For. Sci. 32(4): 989-1000. Hann, D.W. and Ritchie, M.W. 1988. Height growth rate of Douglas fir: a comparison of model forms. For. Sci. 34(1): 165-175. Hann, D.W., Marshall, D.D., and Hanus, M.L. 2003. Equations for predicting height-to- crown-base, 5-year diameter growth rate, 5-year height growth rate, 5-year mortality rate, and maximum size-density trajectory for Douglas-fir in the coastal region of the Pacific Northwest. Forest Research Laboratory, Oregon State University, Corvallis, Oreg. Res. Contr. 40. Harrington, T.B. 2006. Five-year growth responses of Douglas-fir, western hemlock, and western redcedar seedlings to manipulated levels of overstorey and understorey competition. Can. J. For. Res. 36(10): 2439-2453. Harrington, C.A. and Reukema, D.L. 1983. Initial shock and long-term stand development following thinning in a Douglas-fir plantation. For. Sci. 29(1): 33-46. Karlsson, K. and Norell, L. 2005. Modeling survival probability of individual trees in Norway spruce stands under different thinning regimes. Can. J. For. Res. 35(1): 113- 121. 92 Kneeshaw, D.D., Williams, H., Nikinmaa, E., Messier, C. 2002. Patterns of above- and below-ground response of understorey conifer release 6 years after partial cutting. Can. J. For. Res. 32(2): 255-265. Lee, Y. 1971. Predicting mortality for even-aged stands of lodgepole pine. For. Chron. 47(1): 29-32. Lertzman, K.P., Sutherland, G.D., Inselberg, A., and Saunders, S.C. 1996. Canopy gaps and the landscape mosaic in a coastal temperate rain forest. Ecology. 77(4): 1254- 1270. Li, Y., Turnblom, E.C., and Briggs, D.G. 2007. Effects of density control and fertilization on growth and yield of young Douglas-fir plantations in the Pacific Northwest. Can. J. For. Res. 37(2): 449-461. Maguire, D.A., Mainwaring, D.B., and Halpern, C.B. 2006. Stand dynamics after variable-retention harvesting in mature Douglas-fir forests of Western North America. Allg. Forst Jagdztg. 177(6-7): 120-131. Meidinger, D. and Pojar, J., eds. 1991. Ecosystems of British Columbia. The Research Branch, BC Ministry of Forests, Victoria, BC. pp. 95-111. Miller, R.E., and Barker, P.R., Peterson, C., and Webster, S.R. 1986. Using nitrogen fertilizers in management of coast Douglas-fir. In: Douglas-fir: Stand management for the future. Proc. Symp. 18-20 June 1985, Seattle, WA. Univ. Washington, Seattle, Coll. For. Res., Contribution No. 55, pp. 290-303. O’Hara, K.L. 1988. Stand structure and growing space efficiency following thinning in an even-aged Douglas-fir stand. Can. J. For. Res. 18(7): 859-866. 93 Rathbun, L.C., LeMay, V. , and Smith, N.J. 2010. Modeling mortality for mixed-species stands located in coastal British Columbia. Can. J. For. Res. 40(8): 1517-1528. Ritchie, M.W., Skinner, C.N., and Hamilton, T.A. 2007. Probability of tree survival after wildfire in an interior pine forest of northern California: effects of thinning and prescribed fire. For. Ecol. Manage. 247(1-3): 200-208. Shen, G., Moore, J.A., and Hatch, C.R. 2001. The effect of nitrogen fertilization, rock type, and habitat type on individual tree mortality. For. Sci. 47(2): 203-213. Trofymow, J.A., Porter, G.L, Blackwell, B.A., Arksey, R., Marshall, V.G., and Pollard, D.W.F. 1997. Chronosequences for research into the effects of converting coastal British Columbia old-growth forests to managed forests: an establishment report . Can. For. Serv., Pac. For. Cent., Inf. Rep. BC-X-374. pp. 137. Weetman, G.F., Prescott, C.E., Kohlberger, F.L., and Fournier, R.M. 1997. Ten-year growth response of coastal Douglas-fir on Vancouver Island to N and S fertilization in an optimum nutrition trial. Can. J. For. Res. 27(9): 1478-1482. Zenner, E.K., 1995. Effects of residual trees on growth of young to mature Douglas-fir and western hemlock in the Western Central Oregon Cascades. MS Thesis, Oregon State Univ., Corvallis, OR 94 4. Diameter growth models for mixed-species stands of Coastal British Columbia including thinning and fertilization effects 3 4.1 Introduction The mixed-species coastal temperate rain forest of Vancouver Island, British Columbia (BC), Canada, is noted for its lack of recent fires (Lertzman et al. 1996) causing the successional processes of the area to be defined by small-scale disturbances, mainly due to windthrow (Gavin et al. 2003). Hence, gap dynamics define the forest structure and natural successional processes on Vancouver Island (Lertzman et al. 1996), with Douglas-fir (Pseudotsuga menziesii var. menziesii (Mirb.) Franco), as the dominant, pioneer species and western hemlock (Tsuga heterophylla (Raf.) Sarg.) and western redcedar (Thuja plicata Donn) as late-successional species often aggregated within gaps (Getzin et al. 2006). In juvenile stages of individual tree growth, shade-intolerant species, such as Douglas-fir, allocate more photosynthate for height growth rather than diameter growth (Chen and Klinka 2003), outcompeting western hemlock, a shade-tolerant species. It has also been shown that shade-tolerant species grow faster at lower light levels, while shade- intolerant species grow faster at higher light levels (Kobe and Coates 1997). These differences in growth patterns due to levels of shade-tolerance are further complicated by interactions among species in mixed-species stands. Species within mixed-species stands may utilize resources more completely (Piotto 2008), especially if the species have 3 A version of this chapter has been submitted for publication. Rathbun, L.C., LeMay, V., and Smith, N.J. 2010. Diameter growth models for mixed-species stands of coastal British Columbia including thinning and fertilization effects. 95 dissimilar resource requirements and obtain their resources in different ways (Erickson et al. 2009). Due to the structural diversity of crown and root placement within mixed- species stands, positive plant interactions may show a reduction in competition for resources (Chen et al. 2003). In this area, Douglas-fir grows above shade-intolerant western hemlock and it has been suggested that resources between the two species are used in complement (Erickson et al. 2009). Management practices in this area have moved from clearcutting toward variable- retention harvesting (Beese et al. 2003). This shift away from clearcuts further complicates stand structure by diversity within treatments (e.g., variation in thinning intensity, timing of thinnings, and number of thinnings). Thinning has been shown to increase Douglas-fir growth rates, but this effect is dependent upon initial density levels (Wilson and Oliver 2000). Suppressed western redcedar and western hemlock have shown growth increases following thinning, persisting for 10-20 years with 4-year average radial growth rates less than 0.5 mm (Kobe and Coates 1997). Fertilization is a common practice in this area as well. Stands of primarily Douglas-fir have shown a growth response to fertilization (Weetman et al. 1997) as have stands of primarily western redcedar (Devine and Harrington 2009). This response is often temporary and trees tend to resume pretreatment growth rates within 5-10 years of fertilization (Miller 1981). Forest management planning requires diameter growth models as a component of stand development under varying silvicultural treatments (Palahi et al. 2003). To model diameter increment under different silvicultural treatments found on Vancouver Island, a flexible model is required. Also, the model and its components should be logically 96 consistent and biologically realistic (Vanclay and Skovsgaard 1997). The Box and Lucas model (1959) has been used to model diameter increment growth for other species (e.g., Huang 1992). This model has the advantages of a flexible model form and it was developed for processes similar to those governing tree growth. In this study, annual diameter increment models for main coastal BC species growing in mixed-species stands were developed. Because of the flexibility of the Box and Lucas (1959) model. as well as the similarity of processes used to develop this model to growth processes, this model was selected as the base model. First, the model was fitted for each tree diameter growth series of the three main species: Douglas-fir, western redcedar, and western hemlock. A random coefficients modeling (i.e., parameter prediction) approach was then used to replace the estimated parameters by functions of tree size and stage of development, site productivity, and inter-tree competition variables. To estimate the impacts of thinning and fertilization on diameter increment growth, a two-step additive approach was tested, in that diameter increment was modified for fertilization, and then again for thinning. Our hypotheses were: 1) a process-based diameter increment model would outperform an empirical diameter increment model and allow for improved interpretation of the resulting models; 2) plots which were fertilized would experience an initial increase in diameter increment growth during the first few years following application only; 3) modifying only the state variables representing inter- tree competition would result in a model that reflected the initial increase in diameter increment growth during the first few years following thinning; and 4) the increase in diameter increment growth due to the combination of thinning and fertilization could be represented additively. 97 4.2 Materials and methods 4.2.1 Study area The study area was located on Vancouver Island, Haida Gwaii (the Queen Charlotte Islands), and along the coastal mainland within the Coastal Western Hemlock (CWH) and Coastal Douglas-fir (CDF) Biogeoclimatic Ecological Classification (BEC) zones of BC (Meidinger and Pojar 1991). The CWH BEC zone is divided into multiple subzones: a very dry maritime subzone in the east, a moist maritime subzone in the central area, and a very wet maritime and hypermaritime subzone in the west (Meidinger and Pojar 1991). The temperatures in the CWH BEC zone range from 5.2 to 10.5˚ C, with a mean annual temperature of 8˚ C (Meidinger and Pojar 1991). The CDF BEC zone includes the rainshadow of Vancouver Island and the Olympic mountains (Meidinger and Pojar 1991). It is defined by warm, dry summers and mild, wet winters with mean annual temperatures ranging from 9.2 to 10.5˚ C, and a minimum temperature ranging from -21.1 to -11.7˚ C (Meidinger and Pojar 1991). Study plots were located within latitudes ranging from 48.37 to 53.53˚ N and longitudes ranging from 122.36 to 132.59˚ W. The majority of this area is comprised of second growth multi-species stands, regenerated both naturally and from planting. Common tree species include: Douglas-fir, western redcedar, western hemlock, red alder (Alnus rubra Bong.), Sitka spruce (Picea sitchensis (Bong.) Carr.), and yellow cedar (Chamaecyparis nootkatensis (D.Don) Spach). 4.2.2 Data Permanent sample plot (PSP) data were provided by Island Timberlands, LP. The database contained 1,455 untreated plots, ranging in size from 0.008 to 0.253 ha (Table 4- 98 1). Measurements were collected between 1932 to 2002, with measurement intervals varying from 1 to 17 years with an average of 5.0 years. The majority of plots, 1,390, were defined as second growth, less than 150 years since stand establishment; 20 plots were defined as old growth, greater than 150 years since stand establishment; and 45 plots had no available age information. At the time of plot establishment, 515 plots were predominantly Douglas-fir, 97 plots were predominantly western redcedar, and 528 plots were predominately western hemlock, as defined by greater than 50% stems per ha. Based on more than 50% of basal area per ha, 580 plots contained predominantly Douglas-fir, 96 plots predominantly western redcedar, and 474 plots predominately western hemlock. For the remaining plots, species composition averaged 6, 5, and 18% Douglas-fir, western redcedar, and western hemlock, respectively, based on basal area per ha. Densities at plot establishment ranged from 13 to 11,750 live trees per ha and site index using coastal Douglas-fir site index at a base age of 50 years measured at breast height (1.3 m above ground) ranged from 6.2 to 52.8 m. Average annual plot mortality for the first measurement period following plot establishment ranged from 0 to 11.11% with an average of 1.37%. The average annual mortality for the second measurement period was 1.31%, with a slightly wider range of 0 to 12.07%. At each measurement, the species, tree status (i.e., live or dead), and diameter outside bark at 1.3 m above ground (DBH, cm) were recorded for all trees with DBH greater than 4 cm and within plot boundaries. Height (ht, m) was measured for a subset of trees, and remaining heights were estimated using height-diameter functions. At the time of plot establishment, the average DBH values were 12.5, 12.6, and 13.9 cm for Douglas-fir, western redcedar, and western hemlock, respectively, with corresponding 99 maximum values of 149.9, 217.1, and 168.3 cm (Table 4-2). The average heights for Douglas-fir, western redcedar, and western hemlock were 11.7, 11.2, and 14.2 m, respectively. Since the heights were a mixture of measured and estimated heights, the minimum height was 0 rather than the expected value of ≥ 1.3 m. For model fitting, validation, and testing, all non-overlapping repeated measures for all plots were used. An additional 597 plots received silvicultural treatments; treatments included fertilization (85 plots), thinning (419 plots), and a combination of fertilization and thinning (93 plots). The majority of fertilized plots received one application of nitrogen ranging in concentration from 50 to 400 kg per ha and a few plots additionally received ammonium phosphate. Of the 419 plots which were thinned, 344 received one thinning and 75 received multiple thinnings. Basal area cut varied from 0.2 to 82.5 m 2 per ha and stems per ha cut varied from 25 to 10,360 trees per ha, with corresponding averages of 8.2 m 2 per ha and 892 trees per ha. For the plots which received a combination of thinning and fertilization, fertilization was applied within the next measurement period immediately following the thinning treatment. It should be noted that the data used in this study primarily represent silvicultural treatments applied to managed stands rather than treatments applied in a experimental setting. Therefore, the results in this study include the interactions that are often removed from experiments through careful selection of experimental plots. 4.2.3 Box-Lucas model form The Von Bertalanffy (1957) function was developed to model the relationship between animal growth and metabolic functions. The function describes the net gain in 100 an organism’s weight (W) as the difference in the processes which build up and break down materials used for growth and is defined as: [1] ba WW t W where η and κ are constants representing the anabolic and catabolic rates of metabolism, a and b are power coefficients which indicate that these rates are proportional to some power of weight, and t is time. Assuming that the rate of diameter growth with respect to current DBH, where DBH is a substitute for time, can be expressed as the difference between the available constructive metabolism and the destructive metabolism, Eq. [1] can be modified to: inc21 inc AD)DBH( DBH D where Dinc is diameter increment, )DBH(1 is the potential constructive resources available, and incAD2 is the destructive metabolism assumed to be proportional to the current diameter increment modified by a constant, A (Huang 1992). Box and Lucas (1959) showed that if Eq. [1] can be described as a first order and irreversible reaction, W, can be given by: [2] texptexpW Replacing W in Eq. [2] by Dinc and adding in a constant k1, the equation for diameter increment becomes: [3] DBHDBHkD̂ 1inc 12 21 1 expexp 101 where incD̂ is the predicted annual diameter increment, 1k , 1 , and 2 are parameters to be estimated, and DBH was as previously defined. This equation (termed the Box-Lucas equation in the remainder of this paper) is very flexible and is based on metabolic processes similar to processes resulting in diameter growth. Figure 4-1 illustrates typical curves produced from the Box-Lucas function when 1k equals 1, resulting in bounds for the result from 0 to 1. These curves follow the pattern typical of tree diameter increment: diameter increment increases to a maximum, then decreases asymptotically toward zero. The parameters 1 and 2 work in conjunction to describe the net gain in diameter increment where increases in 1 produce larger diameter increments for trees up to a given diameter followed by a decrease in diameter increment, and increases in 2 produce a decrease in diameter increment across all diameters. An increase in 1k maintains the basic shape of the curve while increasing the diameter increment. Using the Box-Lucas model for diameter increment, the DBH at which the maximum diameter increment occurs is given as: 21 12 incmax )/ln(2 DBH where DBHmax inc is the DBH for which the maximum diameter increment occurs and all variables are as previously defined. 4.2.4 Model development for untreated plots To develop diameter increment models for each species and the untreated plots, a random coefficients modeling (i.e., parameter prediction) approach (Clutter et al. 1983 102 pp. 54-56; Schabenberger and Pierce 2002 Chapter 7; Littell et al. 2006, Chapter 8) was used. In the first step, the Box-Lucas model (Eq. [3]) was fit for each tree diameter increment series containing more than three measurement periods using PROC NLIN of SAS version 9.1.3 (SAS Institute Inc. 2004). Estimated parameters 1k , 1 , and 2 , for each tree were then passed to the second step. In the second step, a parameter prediction model was developed for each of the estimated parameters, 1 , and 2 , by species, using functions of tree size and stage of development, site productivity, and inter-tree competition variables. Although an additional parameter prediction model to estimate 1k was initially considered, this was not included for reasons of model parsimony and also because flexible models are obtained with varying 1 and 2 only. However, 1k was replaced by a function when modeling diameter increment in fertilized plots, as later described. Tree size and stage of development were represented by DBH (cm) and height (m). To represent site productivity, coastal Douglas-fir site index (SI, m) at a base age of 50 years breast height age was included. Also, the variable Growth Effective Age (GEA, years) was used to indicate tree productivity, where GEA is the age of a dominant tree of the same height as the tree of interest for a given site index (Hann and Ritchie, 1988). For these data, GEA was estimated from site index equations for the leading tree species in each plot. A number of stand level competition measures were considered. Curtis’ (1982) Relative Density, an index measure of density, was calculated as: qd G RD Curtis 103 n DBH d n 1i 2 i q where DBHi (cm) is the DBH for tree i, dq (cm) is the quadratic mean diameter for the plot, G (m 2 per ha) is the stand basal area per ha, and n is the number of all live trees across all species within a plot. Stand basal area has been commonly used to measure competition for below ground resources (Fan et al., 2006). In addition, tree-level competition measures such as the basal area of larger trees (BAL, m 2 per ha), relative DBH (RDBH), and crown competition factor of larger trees (CCFL) were evaluated for the model. BAL was calculated as: n i ii BAtree δ 1 x BAL where δi is an indicator variable (1 if tree i has a DBH greater than the tree of interest, 0 otherwise), BAtreei (m 2 per ha) is the basal area per ha value for tree i, and n is the number of all live trees across all species within a plot. RDBH is a ratio of the diameter of the tree of interest to the average diameter within a plot and was calculated as: DBH DBH iRDBH where DBHi (cm) is the DBH of tree i and DBH (cm) is the average DBH for the plot. Initially, simple correlations (r) between 1 and 2 estimates by species and possible predictor variables were used to select predictor variables for the parameter prediction models. Transformations were used to linearize the relationships where necessary. Candidate models to predict 1 , and 2 estimates were then fit by species for further selection of predictor variables, using PROC MODEL of SAS version 9.1.3 (SAS 104 Institute Inc. 2004) and the subset of the untreated data used in fitting the Box-Lucas model to each tree. The selection of parameter prediction models for 1 and 2 estimates were based upon Akaike’s Information Criteria (AIC), calculated as: p2ˆ2logL(AIC ) where )ˆlogL( is the log likelihood for the estimated parameter set represented by ˆ , p is the number of parameters in the model, and smaller AIC values indicate a more accurate and/or more parsimonious model. In the third step, the candidate parameter prediction models for 1 , and 2 replaced these parameters in the Box-Lucas model (Eq. [3]) resulting in the following augmented Box-Lucas model: [4] DBHfDBHf ff f kD̂ 1inc 12 21 1 expexp where 1f and 2f are the selected parameter prediction models from step 2 and all other variables are as previously defined. This augmented Box-Lucas model was fit by species using PROC MODEL of SAS version 9.1.3 (SAS Institute Inc. 2004) and the entire dataset. The resultant parameter estimates from the second step were used as an initial set of starting values. In the final step of developing the diameter increment model for untreated plots, the error covariance matrix structure was specified and the model was refit. As noted by Nord-Larsen (2006), repeated measures taken on individual trees can lead to error correlations which may result in inefficient estimates and underestimated standard errors when correlations are strong. To account for temporal autocorrelation, three different error structures were examined for use with the augmented Box-Lucas model: a 105 continuous autocorrelated error structure (CAR(x)); a heteroskedastic error structure using DBH -1 as the weight; and a heteroskedastic, continuous autocorrelated error structure using CAR(x) in combination with DBH -1 as the weight. The CAR(x) structure has been used in other studies (e.g., Crecente-Campo et al. 2009) since this accounts for differing time between measurements. The error structure for first-order CAR(1) is: itqti hh it ee qtiit )( where ite is the ordinary residual for the t th measure of the i th tree, qtie is the ordinary residual of (t-q) th measure for the i th tree, is the first-order continuous autoregressive parameter, )( qtiit hh is the distance separating the t th from the (t-q) th measure of the i th tree, and it is an independent, identically normally distributed error term with mean of 0. The augmented model was then fit again for each species using PROC MODEL of SAS version 9.1.3 (SAS Institute Inc. 2004) and the entire dataset, using the three alternative error covariance matrix structures. The resulting models and the most appropriate error structures were determined using AIC values. For all steps in the fitting process, parameter estimates for 1k were restricted to positive values less than 100 to restrict maximum diameter increment values. Also, initial parameters were varied in order to ensure a global minimum. 4.2.5 Empirical model for comparison The final augmented Box-Lucas model using parameter prediction and untreated plots was compared to a modified model form developed for similar species located in Oregon, U.S. (Hann et al. 2006). The comparison model, termed the Empirical model is: 106 G )DBHlog( BAL )3.1SIlog()DBHlog(expD̂ 43210inc where all variables are as previously defined. The Empirical model was fit using PROC MODEL of SAS version 9.1.3 (SAS Institute Inc. 2004). The same error structures used for the augmented Box-Lucas model were tested and AIC values were used to determine the most appropriate error structure. 4.2.6 Treatment effects 4.2.6.1 Fertilization The augmented Box-Lucas model developed for the untreated plots was modified to include the effects of fertilization. An additional fertilization effect variable, defined as: Ftime F effectfert. where Ftime is the time (years) since fertilization, and F is an indicator variable, defined as 1 for fertilized plots and 0 otherwise, was added to 1f and 2f to obtain 3f and 4f . This approach was used by Carlson et al. (2008) in modeling mean relative growth rate of loblolly pine (Pinus taeda L.) and by Huber et al. (2009) for other growth models. The categorical variable F identifies measurements taken following fertilizer application. As noted earlier, the parameter 1k was also replaced with a parameter prediction model for fertilized plots, as follows: [5] DBHfexpDBHfexp ff f fD̂ 34 43 3 5inc where 5f is a function of the fertilizer effect variable and all other variables are as previously defined. 107 Since modifying all of the parameter prediction models might not have been needed, all possible combinations of modified versus not modified versions of 3f , 4f , and 5f were fitted by species using PROC MODEL of SAS version 9.1.3 (SAS Institute Inc. 2004) for plots that had been fertilized. As with Eq. [4], Eq. [5] was fit using each of the three error structures and AIC values were used to determine the most appropriate model and error structure. 4.2.6.2 Thinning Thinning has often been modeled through the inclusion of thinning parameters designed to describe the type of treatment received. For example, Wimberly and Bare (1996) added a parameter to their basal area increment model defined as the amount of basal area removed. In this study, no modifications were made to the augmented Box- Lucas model developed from the untreated plots for thinning. Instead, the augmented Box-Lucas model given as Eq. [4] was used, with the state variables changed instantaneously after thinning as input variables to predict growth following thinning. In particular, state variables which changed as a result of the thinning were basal area per ha, basal area of larger trees, and stems per ha. This approach was used for young Scots pine (P. sylvestris L.) stands by Pukkala et al. (2002), who found that models without additional thinning variables or parameters can produce unbiased predictions for thinned stands. 4.2.6.3 Thinning and fertilization For plots that received a combination of thinning and fertilizer treatments, the augmented Box-Lucas model [5] modified for fertilized plots was applied directly. As for the plots which only received a thinning, some state variables for these plots 108 immediately change following thinning. This approach assumed a greater response would be seen for plots which received both thinning and fertilization than for those which received only one treatment. 4.2.7 Model accuracy To assess accuracy for all models by species, mean bias and root mean-squared error (RMSE) values by 5 cm diameter classes were calculated. Mean bias was calculated as: m yy m 1i ii ˆ BiasMean where iy is the actual annual diameter increment for tree i, iŷ is the predicted annual diameter increment for tree i, and m is the total number of measures for all trees in the diameter class. Root mean-squared error was calculated as: pm yy m 1i 2 ii ˆ RMSE where p is the number of parameters in the model and all other variables are as previously defined. Mean bias and RMSE was also calculated for all trees by species for the augmented Box-Lucas model [4] and the Empirical model for the untreated plots. In addition, for the untreated plots, graphs of average actual and predicted annual diameter increment by 5 cm diameter class were produced for each species. For fertilized plots, the average actual and predicted annual diameter increments were graphed versus time since treatment using trees on which measurements prior to treatment were taken, and repeated for trees which did not have measures prior to treatment. For both, only 109 trees with at least five measurements were used. Similar graphs were obtained for plots that were thinned. For most plots that were both fertilized and thinned, measurements before treatment were available and only one set of graphs was produced. 4.3 Results 4.3.1 Model development for untreated plots In the first step of the parameter prediction approach, the Box-Lucas model was fitted for each tree in untreated plots having more than three measures. Initial estimates for 1 and 2 were restricted to positive values less than 10; without this restriction the values tended toward infinity. The absolute values of all correlation values between these estimated parameters and basal area per ha was smaller than 0.10 for Douglas-fir (Table 4-3). In addition, for Douglas-fir, the absolute values of correlations between the estimated parameter 1 and basal area of larger trees per ha and Curtis’ relative density were both smaller than 0.10. For western redcedar and western hemlock, the absolute values of correlations between the estimated parameter 2 and basal area of larger trees per ha, site index, and Curtis’ relative density were less than 0.10; absolute correlations with basal area per ha and stems per ha were less than 0.10 for western redcedar. For western redcedar and western hemlock, the absolute values of correlations between the estimated parameter 1 and basal area per ha, site index, Curtis’ relative density, and stems per ha were less than 0.10. Based on the near-normal distribution for the logarithmic transformation of the estimated parameters, an exponential equation form was selected for each parameter prediction equation. Using the simple correlations as a guide, and after fitting a number 110 of possible parameter prediction equations, the selected parameter prediction equations were: For Douglas-fir: )lnexp 8654101 SPHaRDBHaGEAaSIa(DBH)a(af ))ln(exp( 876543102 SPHbCurtisRDbRDBHbGEAbSIbBALbDBHbbf For western redcedar: )RDBHaGEAaBALaGa)DBHln(aaexp(f 6532101 )RDBHbGEAb)DBHln(bbexp(f 65102 For western hemlock: )RDBHaGEAaBALa)DBHln(aaexp(f 653101 ))ln(exp( 8652102 SPHbRDBHbGEAbGbDBHbbf For all three species, an improvement in model fit over using an independent, identically distributed error structure (i.e., No Modification) was obtained using a CAR(1) error structure as expected using these repeated measures data (Table 4-4). The use of a heteroskedastic error structure did not improve the fit, nor did the use of a CAR(1) plus heteroskedastic error structure for any species. Using the CAR(1) error structure, parameter estimates by species for the modified Box-Lucas model (Eq. [4]) are given in Table 4-5. 4.3.2 Empirical model for comparison The best model fit for the Empirical model was found using the same error structure as the augmented Box-Lucas model, CAR(1). Smaller AIC values were found 111 for the augmented Box-Lucas model than for the Empirical model for all species (Table 4-6). Absolute values of mean biases were lower for the augmented Box-Lucas model, although differences were quite small. For all species, the augmented Box-Lucas model had smaller RMSE values than the Empirical model. 4.3.3 Treatment effects 4.3.3.1 Fertilization While only slight changes were seen in AIC values for Douglas-fir and western redcedar when the model was modified to include the fertilizer effect variable, AIC decreased at least two fold from the unmodified model to the modified models for western hemlock, no matter the modification (Table 4-7). For Douglas-fir and western redcedar, the best AIC value was found when 1f and 2f were modified to include the fertilizer effect variable, with no further improvement by introducing 3f as a replacement for 1k . For western hemlock, the addition of 3f to replace 1k resulted in a lower AIC than just modifications of 1f and 2f only. The augmented Box-Lucas model modified to include the fertilizer effect variable (Eq. [5]) and using a CAR(1) error structure was: Ftime F aff 913 exp Ftime F bff 924 exp for all species, and, for Douglas-fir and western redcedar, 15 kf 112 whereas Ftime F ckf 15 for western hemlock. For these equations, the estimated parameters were: -0.62, -3.09, and 0.55 for 9a , and 0.73, 0.11, and 1.21 for 9b , for Douglas-fir, western redcedar, and western hemlock, respectively. The parameter c was 83.59 for western hemlock. 4.3.3.2 Thinning Because the state variables basal area per ha, basal area per ha of larger trees, relative diameter, Curtis RD, and stems per ha change immediately following thinning, the selected models for the thinned plots were the same as that for the untreated plots, with parameter estimates found in Table 4-5. 4.3.3.3 Thinning and fertilization The same reasoning used for the thinned plots applies to the selected models for plots which received a combination of thinning and fertilization. The selected models for plots which received a combination of thinning and fertilization were the same as that found for the fertilized plots (Eq. [5]). 4.3.4 Model accuracy For the untreated plots, the mean bias values by 5 cm diameter class for Douglas- fir were between -0.02 and 0.03 cm, with small positive values found in the two smallest diameter classes and negative values in the larger diameter classes except for trees greater than 55cm DBH (Figure 4-2). Similar results were found for western redcedar with the change from small positive to negative mean bias occurring for trees greater than 25 cm in diameter. For western hemlock, results were similar but without positive mean biases 113 for trees greater than 55cm DBH. The RMSE values were similar for all three species with slightly larger values found in the larger diameter classes, where the values increased with increasing DBH classes (Figure 4-3). As the diameter class increases, the difference between the average actual and predicted values became greater, albeit only slightly (Figure 4-4). This difference was more pronounced for western redcedar and western hemlock than for Douglas-fir. All three graphs in Figure 4-4 showed the greatest difference for trees in diameter classes 40 cm or greater. The change to a positive bias for trees greater than 55 cm in diameter is reflected by the differences seen in Figure 4-4a and b for Douglas-fir and western redcedar, respectively. For the fertilized plots, the mean bias values for Douglas-fir were positive and small for small trees (<10 cm DBH) and became increasingly larger for larger diameter classes. The mean bias values for western redcedar were negative for all diameter classes except for 37.5 cm, and the mean bias values for western hemlock were negative for all diameter classes except for 52.5 cm (Figure 4-2). RMSE values for the fertilized plots were similar to those for the untreated plots for all three species. Figure 4-5 illustrates the change in average actual and predicted diameter increment based on time since fertilization by species. For Douglas-fir where measurements were available prior to treatment, a slight increase in average annual diameter increment immediately following fertilization is shown (Figure 4-5a, c and e). For Douglas-fir where measurements were not available prior to treatment, a decline in diameter increment from time zero to approximately year 10 was shown, after which the diameter increment levels off (Figure 4-5b). For western redcedar and western hemlock, the average predicted and actual values for plots which had measurements prior to 114 treatment showed no increases following treatment (Figure 4-5c and e). However, for plots which did not have measurements prior to treatment, decreases in diameter increments following treatment were shown, followed by a leveling off of diameter increments (Figure 4-5d and f). The mean bias values for the thinned plots were the smallest regardless of diameter class for Douglas-fir, with positive mean bias values for trees smaller than 45 cm and negative mean values for trees greater than 45cm (Figure 4-2). For western redcedar, all mean bias values were positive regardless of diameter class. Western hemlock returned positive mean bias values for trees less than 40 cm DBH, with negative mean bias values for trees greater than 40 cm. RMSE values were similar for all three species (Figure 4-3). No matter the species, average diameter increment increased following thinning (Figure 4-6). This increase in average annual diameter increment was seen for 5 to 10 years following treatment. Western hemlock showed the largest response of the three species, with Douglas-fir and western redcedar showing similar responses. The predicted and actual values showed the largest differences at 5 years after thinning, coinciding with the largest response to thinning. Mean bias and RMSE values for plots which received fertilization following thinning were greatest for western redcedar (Figures 4-2 and 4-3). A slight increase in average annual diameter increment can be seen following treatment (Figure 4-7). The increase was immediately followed by a decline in the average annual diameter increment. 4.4 Discussion The Box-Lucas model used as the basis for modeling diameter increment in this study was developed as a model to describe the process of change, based on incoming 115 and outgoing resources. The model has been used to model animal growth, including fish growth (e.g., Anderson et al. 2005), and has been used to model tree growth (e.g., Huang and Titus 1995), partly because of the underlying development of the model and also because it is very flexible. While the direct interpretation of tree, stand, or landscape- level variables may not translate to the metabolic processes at the cellular level, the model can be indirectly used to interpret the processes involved in the net gain in diameter increment growth, principally, net photosynthesis and respiration. In this study, the Box-Lucas model was fitted by tree and then a random coefficients modeling (i.e., parameter prediction) approach was used to model the changes in growth for tree and stand level variables. For untreated plots, an alternative Empirical model was also fitted for each species and compared to the augmented Box- Lucas model. The Empirical model assumes a basic exponential form, where diameter increment increases with increasing DBH. However, the Box-Lucas model is more flexible and the use of a parameter prediction approach better reflects the nature of the repeated measured tree data. As a result, the augmented Box-Lucas model outperformed the Empirical model for all three species. For untreated Douglas-fir, the pioneer species in this area, variables representing tree size and stage of development and site productivity were included in the parameter prediction equation for 1 , representing the net photosynthetic process, but the competition variables basal area per ha and basal area per ha of larger trees were not. Shade intolerant species such as Douglas-fir grow to the dominant and codominant canopy positions where tree size varies less and competition for light is less severe on a site than for those species in the understorey (Carter and Klinka 1992). Because these 116 dominant trees are similar in height, differences in competition for light resources may be better represented by the relative diameter variable, RDBH. Basal area per ha was also not included for the parameter prediction equation for 2, representing the respiration process. For western hemlock, density measures such as basal area per ha and stems per ha were not well correlated with 1, but basal area per ha of larger trees was. In this coastal area, western hemlock is a shade-tolerant species which can persist in the understorey for many years with extremely small growth rates but can grow quickly when light becomes available (Kobe and Coates 1997). Hence, the basal area per ha of larger trees, a more direct measure of competition for light for an individual tree, may be a more important variable in photosynthetic growth response for western hemlock. Basal area per ha of larger trees was not well correlated with 2 for western hemlock. Unlike the other two species, variables representing tree size and stage of development, site productivity, and competition were all moderately correlated with 1 for western redcedar. Western redcedar has been shown to be less sensitive to site variations and can persist in areas with low nutrient availability and poor soil drainage (Minore 1990). This ability to grow under varying site conditions may indicate the need to include all measurement levels to describe photosynthetic growth response. Only DBH, GEA, and RDBH were correlated with 2 for western redcedar. Site index was not correlated with 1 nor 2 for western redcedar nor western hemlock; however GEA was significant and is derived from site index equations. A slight increase in Douglas-fir diameter increment was found following fertilization in this study. In other studies, Douglas-fir has been shown to respond to fertilization (e.g., Chappell et al. 1990). It should be noted that the impacts on growth 117 due to fertilization are not only a direct result of nutritional improvement but also an indirect effect due to changes in stand structure through time (McWilliams and Therien 1997) and the magnitude of response to fertilization for Douglas-fir is dependent upon site index and soil characteristics (Brix 1992). Unlike our hypothesis, the applications of fertilizer alone did not show an initial increase in average annual diameter increment after fertilization for western redcedar. Evidence in this area indicates that the application of fertilization for western redcedar on poor-quality sites might sustain an increase in tree growth (Blevins et al. 2006). While a possible increase was seen for western hemlock for those plots which had measurements following treatment, this increase was not seen for plots which had measurements prior to treatment. Previous studies in this area have shown that western hemlock response to fertilizer varies (Omule and Britton 1991). Impacts of thinning on diameter increment were addressed using the same parameter estimates as the untreated plots, allowing only the state variables to change. Since the main impact of thinning is the reduction of competition, this response is reflected in the post-thinning state variables such as basal area per ha, stems per ha, and basal area of larger trees. As was hypothesized, the model reflects the post-thinning response well for Douglas-fir. For western redcedar and western hemlock, the model under-predicts for the first 10 years following thinning, but predicts well afterwards. The largest differences between the average actual and predicted annual diameter increments were found for western hemlock. Western hemlock, a shade tolerant species, has been shown to persist within the understorey, producing relatively small growth rates (Kobe and Coates 1997), with residual trees able to grow rapidly following a thinning (Deal, 2001). This instantaneous change in competition for light resources which occurs after 118 thinning is only reflected in the model with the inclusion of basal area per ha of larger trees, which is specific to individual trees. Other stand-level measures of competition for resources such as basal area per ha and stems per ha are not reflected within the model for western hemlock. Competition measures included in the model were based on correlation values derived from the untreated plots, and may not reflect the correlation values for plots after thinning occurs. While the fertilized plots alone did not show an increase in average diameter increment, a more clearly defined increase was seen for fertilizer following thinning, especially for western redcedar and western hemlock. An increase was shown within the first 10 years after application. Thomson and Barclay (1984) have shown that the combination of thinning and fertilization increases the area increment for Douglas-fir in this area. Brix (1992) noted that nitrogen uptake after fertilization was similar in unthinned and heavily thinned plots over a 9-year period, with uptake taking place in the first year in unthinned plots and later in thinned plots for Douglas-fir. The model predicted well on average for all three species which received fertilization and thinning treatments, indicating support for our hypothesis that the increase can be represented additively. As noted, the data used in this study primarily represent silvicultural treatments applied to managed stands rather than treatments applied in a experimental setting; the intent was to obtain diameter growth models as a component of a growth model for use on these stands. As a result, the impacts of fertilization and thinning were more variable than in studies based on experiments. 119 4.5 Conclusion The Box and Lucas model was used to model tree-level diameter increment for three BC coastal species and a parameter prediction approach was used to determine variables affecting growth parameters. The impact of competition variables in each model may have been influenced heavily by the level of shade tolerance and species response within untreated conditions. When applied directly to thinned trees, the competition variables which were not significant for untreated stands, such as basal are per ha for shade tolerant species, may be significant and should be taken into account in the model building process. For the thinned plots, the model predicted well for Douglas- fir and western redcedar, but did a poorer job for western hemlock. This may be due to the exclusion of basal area per hectare in the model and further illustrates the need to maintain competition variables in the model which may not be significant. The model predicted well for fertilized plots across all species. While the model predicted well for plots which were both fertilized and thinned, further research is needed for plots which receive thinning and fertilization in combination. 120 Figure 4-1. Box-Lucas model for varying 1 and 2 values where k1=1; all forms indicate an increase in diameter increment to a maximum diameter increment followed by a decrease to a lower limit. 121 Figure 4-2. Mean bias (cm) values by 5 cm DBH classes for: a) Douglas-fir; b) western redcedar; and c) western hemlock by treatment. 122 Figure 4-3. RMSE (cm) values by 5 cm DBH classes for: a) Douglas-fir; b) western redcedar; and c) western hemlock by treatment. 123 Figure 4-4. Average actual (●) and predicted (♦) annual diameter increments (Dinc) for the untreated data by 5 cm DBH class for: a) Douglas-fir, b) western redcedar, and c) western hemlock. The numbers of trees used to determine average values are given across the top of each graph. 124 Figure 4-5. Average actual (●) and predicted (♦) annual diameter increments (Dinc) for the fertilized plots by time since fertilization for plots with measurements prior to fertilization: a) Douglas-fir, c) western redcedar, and e) western hemlock and for plots without measurements prior to fertilization: b) Douglas-fir, d) western redcedar, and f) western hemlock. The time of fertilization is given as 0, and the numbers of trees used to determine the averages are given across the top of each graph. 125 Figure 4-6. Average actual (●) and predicted (♦) diameter increments (Dinc) for the thinned plots by time since thinning for plots with measurements prior to thinning: a) Douglas-fir, c) western redcedar, and e) western hemlock and for plots without measurements prior to thinning: b) Douglas-fir, d) western redcedar, and f) western hemlock. The time of thinning is given as 0, and the numbers of trees used to determine the averages are given across the top of each graph. 126 Figure 4-7. Average actual (●) and predicted (♦) annual diameter increment (Dinc) for the multiple treatment data for a) Douglas-fir, b) western redcedar, and c) western hemlock over time for thinned and fertilized plots. The time of thinning is given as 0) and the number of trees used to determine the average value is at the top of the plot. 127 Table 4-1. Plot summary statistics at plot establishment and at time of treatment. (n= number of plots). Treatment Type Variable a Plot Establishment Time of Treatment Mean Min b Max c Mean Min b Max c Untreated n=1,455 SPH 2,044 13 11,750 G (m 2 per ha) 38.7 0.1 185.0 SI (m) 29.3 6.2 52.8 CurtisRD 8.7 0.1 24.2 Fertilized n=85 SPH 2,406 370 8,123 2,243 370 8,123 G (m 2 per ha) 35.3 1.4 96.1 39.3 1.4 96.1 SI (m) 31.9 11.2 42.3 32.1 11.2 42.3 CurtisRD 8.7 0.7 18.8 9.3 0.7 17.9 Thinned n=419 SPH 2,126 309 11,370 2,113 309 11,370 G (m 2 per ha) 29.4 0.5 110.5 31.2 0.8 110.5 SI (m) 31.5 9.2 50.0 31.6 9.2 50.0 CurtisRD 7.1 0.3 21.1 7.5 0.5 21.1 Multiple Treatments n=93 SPH 2,332 530 7,802 2,295 530 7,605 G (m 2 per ha) 45.2 5.0 91.9 47.6 5.0 91.9 SI (m) 32.7 8.9 42.9 32.7 8.9 42.9 CurtisRD 10.7 2.2 17.4 11.1 2.2 17.4 a SPH is stems per ha, G is basal area per ha, SI is Site Index, CurtisRD is Curtis’ Relative Density, b Min is Minimum, and c Max is Maximum. 128 Table 4-2. Means and standard deviations (in brackets) for tree-level variables at plot establishment by silvicultural treatment and species. Species Variable a Untreated Fertilized Thinned Multiple Treatments Douglas-fir Number of trees 33,960 5,040 14,903 4,492 DBH (cm) 12.5 (10.1) 11.3 (7.6) 16.2 (10.4) 15.6 (10.0) Height (m) 11.7 (7.5) 11.4 (6.2) 15.5 (8.3) 15.2 (7.2) BAL (m 2 per ha) 18.4 (15.4) 22.0 (14.6) 20.6 (14.6) 24.9 (13.4) RDBH 1.1 (0.5) 1.0 (0.4) 1.1 (0.4) 1.2 (0.5) GEA (years) 17.4 (16.5) 15.3 (9.5) 21.8 (16.2) 18.9 (12.4) Dinc (cm) 0.42 (0.42) 0.34 (0.30) 0.46 (0.50) 0.36 (0.38) Western redcedar Number of trees 11,106 419 2,904 360 DBH (cm) 12.6 (12.1) 11.3 (10.0) 12.5 (12.0) 17.0 (11.6) Height (m) 11.2 (7.6) 10.4 (6.6) 11.6 (8.3) 15.7 (7.5) BAL (m 2 per ha) 33.1 (25.0) 27.6 (18.6) 32.4 (22.1) 28.2 (14.1) RDBH 0.9 (0.5) 1.0 (0.5) 0.9 (0.6) 1.2 (0.7) GEA (years) 19.5 (17.5) 13.4 (9.8) 15.9 (13.6) 22.0 (17.6) Dinc (cm) 0.26 (0.30) 0.35 (0.33) 0.37 (0.38) 0.40 (0.38) Western hemlock Number of trees 49,092 5,845 15,087 3,062 DBH (cm) 13.9 (11.0) 12.1 (8.2) 12.4 (8.7) 19.7 (10.3) Height (m) 14.2 (9.5) 12.1 (7.6) 11.5 (7.3) 19.9 (8.3) BAL (m 2 per ha) 34.6 (23.5) 25.3 (21.0) 22.0 (20.5) 37.5 (22.5) RDBH 1.0 (0.4) 1.0 (0.4) 1.1 (0.5) 1.2 (0.5) GEA (years) 21.8 (19.1) 15.8 (12.4) 15.3 (12.0) 27.0 (14.2) Dinc (cm) 0.28 (0.30) 0.65 (0.67) 0.68 (0.59) 0.30 (0.34) a DBH is diameter at breast height, BAL is basal area per ha of larger trees, RDBH is relative diameter at breast height, GEA is growth effective age, and Dinc is the annual diameter increment. 129 Table 4-3. Correlations between possible predictor variables and estimated 1 or 2 for trees with more than three measurements. Variables a Douglas-fir (n=17,933) Western redcedar (n=4,167) Western hemlock (n=19,136) 1 2 1 2 1 2 DBH -0.273 -0.300 -0.184 -0.252 -0.226 -0.292 ln(DBH) -0.280 -0.382 -0.193 -0.307 -0.240 -0.328 DBH 2 -0.242 -0.198 -0.134 -0.169 -0.183 -0.232 G -0.083 0.066 0.106 -0.053 0.015 -0.148 G 0.5 -0.076 0.084 0.097 -0.050 0.009 -0.138 BAL 0.070 0.313 0.222 0.063 0.208 0.018 SI -0.165 -0.132 0.025 0.011 -0.082 -0.068 GEA -0.128 -0.188 -0.130 -0.233 -0.152 -0.265 RDBH -0.201 -0.418 -0.245 -0.205 -0.283 -0.200 Curtis’ RD -0.002 0.126 0.090 0.002 0.042 -0.026 SPH 0.132 0.148 -0.035 0.093 0.067 0.196 a DBH is diameter at breast height, G is basal area per ha, BAL is basal area per ha of larger trees, SI is Site Index, GEA is growth effective age, RDBH is relative diameter at breast height, Curtis’ RD is Curtis’ Relative Density, and SPH is stems per ha. 130 Table 4-4. AIC values for augmented Box-Lucas model (Eq. 4) with varying error structures for trees on untreated plots (n=numbers of trees). Error Structure Douglas-fir (n=43,059) Western redcedar (n=17,225) Western hemlock (n=65,563) No modification -174,396 -47,232 -251,792 CAR(1) -197,800 -54,974 -283,202 Weighted using DBH -1 -177,562 -50,208 -228,148 CAR(1) and Weighted using DBH -1 -191,954 -55,830 -248,898 131 Table 4-5. Parameter estimates for the augmented Box-Lucas model (Eq. 4) for untreated plots. a N/A indicates the variable is not included in the model Function Parameter Parameter Estimates a Douglas- fir Western redcedar Western hemlock 1f a0 2.9339 1.7655 2.7256 a1 -1.1064 -1.9839 -1.4380 a2 N/A 0.0032 N/A a3 N/A 0.0437 0.0062 a4 -0.0135 N/A N/A a5 -0.0096 0.0575 0.0005 a6 -0.1674 0.2988 -0.3037 a7 N/A N/A N/A a8 -0.0001 N/A N/A 2f b0 1.6337 1.8315 1.4205 b1 -1.2119 -1.1062 -1.0870 b2 N/A N/A 0.0137 b3 0.0091 N/A N/A b4 -0.0017 N/A N/A b5 0.0152 -0.0023 0.0106 b6 0.1072 -0.0692 0.2438 b7 0.0515 N/A N/A b8 1.66E-6 N/A -7.12E-6 1k 35.5176 6.8175 24.0531 132 Table 4-6. AIC, mean bias (cm), and root mean-squared error (RMSE, cm) values by species for the augmented Box-Lucas model vs. the Empirical model for untreated plots (n= number of trees). Fit Statistic Douglas-fir (n=43,059) Western redcedar (n=17,225) Western hemlock (n=65,563) Box- Lucas Empirica l Box- Lucas Empirica l Box- Lucas Empirica l AIC -197,800 -183,714 -54,974 -53,344 -283,202 -274,174 Mean Bias (cm) 0.0028 0.0035 0.0091 0.0097 -0.0021 0.0022 RMSE (cm) 0.1574 0.1632 0.1629 0.1674 0.1474 0.1478 133 Table 4-7. AIC values for the augmented Box-Lucas model modified to include the fertilizer variable for fertilized plots (n= number of trees). Modification AIC Value Douglas-fir (n=4,231) Western redcedar (n=380) Western hemlock (n=5,081) 3f -28,242 -940 -6,350 4f -28,160 -896 -7,138 3f and 4f -28,548 -948 -7,136 5f -28,178 -898 -7,548 5f and 3f -28,412 -945 -7,558 5f and 4f -28,412 -896 -7,586 5f , 3f , and 4f -28,546 -948 -7,694 No modification -28,187 -888 -2,333 134 4.6 References Anderson, M.J., Millar, R.B., Blom, W.M., Diebel, C.E. 2005. Nonlinear multivariate models of successional change in community structure using the von Bertalanffy curve. Oecologia 146(2):279-286. Beese, W.J., Dunsworth, B.G., Zielke, K., and Bancroft, B. 2003. Maintaining attributes of old-growth forests in coastal B.C. through variable retention. For. Chron. 79(3): 570-578. Blevins, L.L., Prescott, C.E., Van Niejenhuis, A. 2006. The roles of nitrogen and phosphorus in increasing productivity of western hemlock and western redcedar plantations on northern Vancouver Island. For. Ecol. Manage. 234(1-3):116-122. Box, G.E.P. and Lucas, H.L. 1959. Design of experiments in non-linear situations. Biometrika 46(1/2): 77-90. Brix, H. 1992. Fertilization and thinning effect on a Douglas-fir ecosystem at Shawnigan Lake: A synthesis of project results. For.Can. and B.C. Min. For., Victoria, BC FRDA Rep. No 196. Carlson, C.A., Burkhart, H.E., Allen, H.L., and Fox, T.R. 2008. Absolute and relative changes in tree growth rates and changes to the stand diameter distribution of Pinus taeda as a result of midrotation fertilizer applications. Can. J. For. Res. 38(7): 2063- 2071. Carter, R.E. and Klinka, K. 1992. Variation in shade tolerance of Douglas fir, western hemlock, and western red cedar in coastal British Columbia. For. Ecol. Manage. 55(1-4): 87-105. 135 Chappell, H.N., Cole, D.W., Gessel, S.P., Walker, R.B. 1990. Forest fertilization research and practice in the Pacific Northwest. Fert. Res. 27(1):12-140. Chen, H.Y.H. and Klinka, K. 2003. Aboveground productivity of western hemlock and western redcedar mixed-species stands in southern coastal British Columbia. For. Ecol. Manage. 184(1-3): 55-64. Chen, H.Y.H., Klinka, K., Mathey, A.H., Wang, X., Varga, P., Chourmouzis, C. 2003. Are mixed-species stands more productive than single-species stands: an empirical test of three forest types in British Columbia and Alberta. Can. J. For. Res. 33(7):1227-1237. Clutter, J.L., Fortson, J.C., Pienaar, L.V., Brister, G.H., Bailey, R.L. 1983. Timber Management – A Quantitative Approach. John Wiley & Sons, New York. 333 pp. Crecente-Campo, F., Marshall, P., LeMay, V., and Dieguez-Aranda, U. 2009. A crown profile model for Pinus radiata D. Don in northwestern Spain. For. Ecol. Manage. 257(12): 2370-2379. Curtis, R.O. 1982. A simple index of stand density for Douglas-fir. For. Sci. 28(1):92-94. Deal, R.L. 2001. The effects of partial cutting on forest plant communities of western hemlock- Sitka spruce stands in southeast Alaska. Can. J. For. Res. 31(12): 2067- 2079. Devine, W.D. and Harrington, C.A. 2009. Western redcedar response to precommercial thinning and fertilization through 25 years posttreatment. Can. J. For. Res. 39(3): 619-628. 136 Erickson, H.E., Harrington, C.A., and Marshall, D.D. 2009. Tree growth at stand and individual scales in two dual-species mixture experiments in southern Washington State, USA. Can. J. For. Res. 39(6): 1119-1132. Fan, Z., Kabrick, J.M., and Shifley, S.R. 2006. Classification and regression tree based survival analysis in oak-dominated forests of Missouri’s Ozark highlands. Can. J. For. Res. 36(7): 1740-1748. Gavin, D.G., Brubaker, L., and Lertzman, K.P. 2003. Holocene fire history of a coastal temperate rain forest based on soil charcoal radiocarbon dates. Ecology 84(1) 186- 204. Getzin, S., Dean, C., He, F., Trofymow, J.A., Wiegand, K. and Wiegand, T. 2006. Spatial patterns and competition of tree species in a Douglas-fir chronosequence on Vancouver Island. Ecography 29(5): 671-682. Hann, D.W. and Ritchie, M.W. 1988. Height growth rate of Douglas fir: a comparison of model forms. For. Sci. 34(1): 165-175. Hann, D.W., Marshall, D.D., Hanus, M.L. 2006. Reanalysis of the SMC-ORGANON equations for diameter-growth rate, height-growth rate, and mortality rate of Douglas- fir. Research contribution 49, Forest Research Laboratory, Oregon State University, Corvallis, OR. Huang, S. 1992. Nonlinear simultaneous diameter and height growth models for major Alberta tree species. PhD Dissertation, University of Alberta, Edmonton, AB. pp.107-129. Huang, S. and Titus, S.J. 1995. An individual tree diameter increment model for white spruce in Alberta. Can. J. For. Res. 25(9): 1455-1465. 137 Huber, M., Halmschlager, E., and Sterba, H. 2009. The impact of forest fertilization on growth of mature Norway spruce affected by Sirococcus shoot blight. For. Ecol. Manage. 257(6): 1489-1495. Kobe, R.K. and Coates, K.D. 1997. Models of sapling mortality as a function of growth to characterize interspecific variation in shade tolerance of eight tree species of northwestern British Columbia. Can. J. For. Res. 27(2): 227-236. Lertzman, K.P., Sutherland, G.D., Inselberg, A. and Saunders, S.C. 1996. Canopy gaps and the landscape mosaic in a coastal temperate rain forest. Ecology. 77(4): 1254- 1270. Littell, R.C., Milliken, G.A., Stroup, W.W., Wolfinger, R.D., and Schabenberger, O. 2006. SAS for Mixed Models, 2nd ed. SAS Institute Inc. Cary, NC. McWilliams, E.R.G., Therien, G. 1997. Fertilization and thinning effects on a Douglas- fir ecosystem at Shawnigan Lake: 24-year growth response. For. Can. and BC Min. For., Victoria, BC FRDA Rep. No. 269. Meidinger, D. and Pojar, J., eds. 1991. Ecosystems of British Columbia. The Research Branch, BC Ministry of Forests, Victoria, BC. pp. 95-111. Miller, H.G. 1981. Forest fertilization: some guiding concepts. Forestry 54(2):157-167. Minore, D. 1990. Thuja plicata Donn ex. D. Don, Western redcedar. In: Silvics of North America. Vol 1. Conifers. Edited by R.M. Burns and B.H. Honkala (technical coordinators). USDA For. Serv. Agric. Handb. 654. pp. 590-600. Nord-Larsen, T. 2006. Modeling individual-tree growth from data with highly irregular measurement intervals. For. Sci. 52(2): 198-208. 138 Omule, S.A.Y., Britton, G.M. 1991. Basal area response nine years after fertilization and thinning western hemlock. For. Can. and B.C. Min. For., Victoria, BC FRDA Rep. No. 137. Palahi, M., Pukkala, T., Miina, J., and Montero, G. 2003. Individual-tree growth and mortality models for Scots pine (Pinus sylvestris L.) in north-east Spain. Ann. For. Sci. 60(1): 1-10. Piotto, D. 2008. A meta-analysis comparing tree growth in monocultures and mixed plantations. For. Ecol. Manage. 255(3-4): 781-786. Pukkala, T., Miina, J., and Palahi, M. 2002. Thinning response and thinning bias in a young scots pine stand. Silvica Fennica. 36(4): 827-840. SAS Institute Inc. 2004. SAS/STAT 9.1 user’s guide. SAS Institute Inc., Cary, N.C. Schabenberger, O., and Pierce, F.J. 2002. Contemporary Statistical Models for the Plant and Soil Sciences. CRC Press. Boca Raton, Florida. Thomson, A.J., Barclay, H.J. 1984. Effects of thinning and urea fertilization on the distribution of area increment along the boles of Douglas-fir at Shawnigan Lake, British Columbia. Can. J. For. Res. 14(6):879-884. Vanclay, J.K., Skovsgaard, J.P. 1997. Evaluating forest growth models. Ecol. Model. 98(1):1-12. Von Bertalanffy, L. 1957. Quantitative laws in metabolism and growth. Q. Rev. Biol. 32(3): 217-231. Weetman, G.F., Prescott, C.E., Kohlberger, F.L., and Fournier, R.M. 1997. Ten-year growth response of coastal Douglas-fir on Vancouver Island to N and S fertilization in an optimum nutrition trial. Can. J. For. Res. 27(9): 1478-1482. 139 Wilson, J.S., Oliver, C.D. 2000. Stability and density management in Douglas-fir plantations. Can. J. For. Res. 30(6):910-920. Wimberly, M.C. and Bare, B.B. 1996. Distance-dependent and distance-independent models of Douglas-fir and western hemlock basal area growth following silvicultural treatment. For. Ecol. Manage. 89(3): 1-11. 140 5. Height growth models: Potential, average, and responses to silvicultural treatments 4 5.1 Introduction Forest management planning requires height growth models as components in models of stand development under varying silvicultural treatments (Palahi et al. 2003). As well as average height growth, potential height growth, defined as the maximum possible height which can be achieved, is of interest. Potential height is an indicator of site productivity often used in site index equations, and also indicates the maximum achievable height under management. For some growth models, potential height growth is used as the limiting height in calibrating growth models. In applying growth models, a measure of potential height may more often be available, since height measurements in forest inventory may be focused on site trees for determining site index (Temesgen et al. 2008). However, to grow all trees, models of average height growth (hereafter termed “height growth”) are needed. In developing models to predict potential height growth, a subset of sample trees that are perceived to be growing at the maxima need to be selected. Commonly, trees have been selected using the tallest trees as defined by site quality (e.g., Hann and Richie 1986) or a specified definition of dominant/co-dominant status (Ishii et al. 2000). In terms of modeling height growth, two methods have been used: the growth potential dependent method, where height growth is predicted as a function of potential height growth; and the growth potential independent method, where height growth is directly 4 A version of this chapter will be submitted for publication. Rathbun, L.C. 2010. Height growth models: Potential, average, and responses to silvicultural treatments. 141 predicted using tree, stand, and site variables (Ek and Monserud 1974). The growth potential dependent method results in logically consistent estimates of height growth and potential height growth, but may not be as accurate as modeling height growth directly. Management practices, such as fertilization and thinning, impact height growth and are often implemented to maximize diameter growth and reach potential height. Fertilization increases below ground resources through the addition of nitrogen and phosphorous to the soil. When resources are made available, species which are adapted to high-resource habitats will grow faster than species which are adapted to low-resource habitats (Grime and Hunt 1975). Conversely, species adapted to low-resource habitats grow faster than those that are not adapted when resources are limited (Fichtner and Schultze 1992, Walters et al. 1993). The increase in height growth expected from fertilization depends upon a species shade-tolerance and on natural growing conditions. Shade-tolerant species grow faster at lower light levels, whereas shade-intolerant species grow faster at high light levels (Kobe and Coates 1997). Fertilized trees are likely to allocate resources directly to height growth, unless thinning is also performed. Thinning primarily increases light levels (above ground resources) within the canopy. Most species will reduce growth in diameter due to competition well before reducing height growth (Piotto 2008). As a result, trees in thinned stands commonly allocate resources directly to diameter growth. Any increase in height growth from thinning depends upon species shade-tolerance, as with fertilizer impacts, and thinning intensity. Particularly in even-aged stands, height growth may not be affected by stand density except at the extremes of very high or very low densities (Oliver and Larson 1996). 142 In this study, potential height growth and height growth models for major species of coastal British Columbia (BC), Canada were developed. The Box-Lucas model (1959), a flexible model, was used as the base model for developing height growth models for Douglas-fir (Pseudotsuga menziesii var. menziesii (Mirb.) Franco), western hemlock (Tsuga heterophylla (Raf.) Sarg.), and western redcedar (Thuja plicata Donn). Potential height growth models were developed using a subset of the data representing potential growth. Two methods were used to select the subset of data: a binning approach, where the data were binned by diameter and site index classes and then the observations with the largest height growth were selected from each bin; and a prediction interval approach, whereby a simple average height growth model was fitted and observations with height growth larger than an upper percentile of prediction intervals were selected. Height growth was modeled using both the potential height growth dependent and independent methods. The main hypotheses associated with model development were: 1) for potential height growth, the prediction interval approach would better define the upper height growth values than the binning approach; and 2) little difference would be seen between the growth potential dependent and independent methods for developing height growth models. Models were also modified to predict the impacts of thinning and fertilization on height growth. Because this study dealt with silvicultural treatments applied to managed stands (i.e. it is not an experimentally-designed study), treated stands were classified by site quality and density levels and compared to untreated stands within the same classes. The main hypotheses associated with management practices were: 1) height growth for the fertilized plots would be greater than the untreated plots, and the magnitude of the difference would depend upon species and site characteristics; 2) modifying only the state 143 variables for inter-tree competition would reflect the effect of thinning on height growth; 3) height growth for the thinned plots would be the same or less than the untreated plots, and the magnitude of the difference would depend upon species and site characteristics; 4) the largest difference in height growth between treated and untreated stands would be found for stands which received fertilization and thinning in combination; and 5) no matter the treatment, height growth would approach but not reach the potential height growth. 5.2 Materials and methods 5.2.1 Study area The study area was located on Vancouver Island, Haida Gwaii (the Queen Charlotte Islands), and along the coastal mainland within the Coastal Western Hemlock (CWH) and Coastal Douglas-fir (CDF) Biogeoclimatic Ecological Classification (BEC) zones of BC (Meidinger and Pojar 1991). The CWH BEC zone is divided into multiple subzones: a very dry maritime subzone in the east, a moist maritime subzone in the central area, and a very wet maritime and hypermaritime subzone in the west (Meidinger and Pojar 1991). The temperatures in the CWH BEC zone range from 5.2 to 10.5˚ C, with a mean annual temperature of 8˚ C (Meidinger and Pojar 1991). The CDF BEC zone includes the rainshadow of Vancouver Island and the Olympic mountains (Meidinger and Pojar 1991). It is defined by warm, dry summers and mild, wet winters with mean annual temperatures ranging from 9.2 to 10.5˚ C, and a minimum temperature ranging from -21.1 to -11.7˚ C (Meidinger and Pojar 1991). Study plots were located within latitudes ranging from 48.37 to 53.53˚ N and longitudes ranging from 122.36 to 132.59˚ W. The majority of this area is comprised of second growth multi-species stands, regenerated both naturally and 144 from planting. Common tree species include: Douglas-fir, western redcedar, western hemlock, red alder (Alnus rubra Bong.), Sitka spruce (Picea sitchensis (Bong.) Carr.), and yellow cedar (Chamaecyparis nootkatensis (D.Don) Spach). The successional processes of the mixed-species coastal temperate rain forest of British Columbia (BC) are defined by small-scale disturbances, mainly due to windthrow (Gavin et al. 2003). Hence, gap dynamics define the forest structure and natural successional processes on Vancouver Island (Lertzman et al. 1996), with Douglas-fir as the dominant, pioneer species and western hemlock, and western redcedar as late- successional species often aggregated within gaps (Getzin et al. 2006). In gaps, the average height of Douglas-fir has been known to increase up to 15 percent (York et al. 2003). In juvenile stages of growth, shade-intolerant species, such as Douglas-fir, utilize photosynthate to increase height growth rather than diameter growth (Chen and Klinka 2003), outcompeting western hemlock, a shade-tolerant species. Shade-tolerant species such as western hemlock are commonly regenerated under a canopy of mature trees, making suppression and release more common for this species (Franklin et al. 2002). Height growth in both suppressed and released western hemlock has been shown to exceed that for Douglas-fir (Renninger et al. 2007). Drever and Lertzman (2001) noted that western redcedar is more shade tolerant than Douglas-fir on Vancouver Island, typically growing faster at lower levels of light. Management practices in this area have moved from clearcutting toward variable- retention harvesting (Beese et al. 2003). This shift away from clearcuts further complicates stand structure by introducing diversity within treatments (e.g., variation in 145 thinning intensity, timing of thinnings, number of thinnings, and spatial patterns). In addition to thinning, fertilization is also a common practice in this area. 5.2.2 Data Permanent sample plot (PSP) data were provided by Island Timberlands, LP. The database contained 1,316 untreated plots, ranging in size from 0.008 to 0.253 ha (Table 5- 1) 5 . Measurements were collected between 1932 to 1996, with measurement intervals varying from 1 to 30 years with an average of 4.8 years. The majority of plots, 1,277, were defined as second growth less than 150 years since establishment, six plots were defined as old growth greater than 150 years since establishment, and 33 plots had no available age information. At the time of plot establishment, 601 plots were predominantly Douglas-fir, 94 plots were predominantly western redcedar, and 621 plots were predominately western hemlock, as defined by greater than 50% stems per ha. Based on more than 50% of basal area per ha, 618 plots contained predominantly Douglas-fir, 92 plots predominantly western redcedar, and 584 plots predominately western hemlock. For the remaining plots, species composition averaged 35, 28, and 37% Douglas-fir, western redcedar, and western hemlock, respectively, based on basal area per ha. Densities at plot establishment ranged from 222 to 11,750 live trees per ha and site index, using coastal Douglas-fir site index at a base age of 50 years breast height age (1.3 m above ground), ranged from 10.7 to 52.8 m. At each measurement, the species, tree status (i.e., live or dead), height, and diameter outside bark at 1.3 m above ground (DBH, cm) were recorded for all trees with DBH greater than 4 cm and within plot boundaries. At the time of plot establishment, the 5 The numbers of PSPs differ from Chapter 4 since heights were not measured on every tree. 146 average DBH values were 17.1, 19.9, and 22.9 cm for Douglas-fir, western redcedar, and western hemlock, respectively, with corresponding maximum values of 146.6, 135.5, and 168.3 cm (Table 5-2). The average heights for Douglas-fir, western redcedar, and western hemlock were 13.7, 13.9, and 19.8 m, respectively. For model fitting, validation, and testing, all non-overlapping repeated height growth measures of all plots were used. 5.2.3 Box-Lucas model The Box-Lucas model (1959) is based on metabolic processes similar to processes resulting in height growth and is defined as: [1] HHkH inc 12 21 1 1 expexp ˆ where incĤ is the predicted annual height growth; 1k , 1 , and 2 are parameters to be predicted; and H is current height. This equation is very flexible; a more extensive description can be found in Rathbun et al. (2010). Figure 5-1 illustrates typical curves produced from the Box-Lucas function when 1k equals 1, resulting in bounds for the result from 0 to 1. These curves follow the pattern typical of tree height growth where height growth increases to a maximum, then decreases asymptotically toward zero. The parameters 1 and 2 work in conjunction to describe the net gain in height growth: increases in 1 produce larger height growth for trees up to a given height followed by a decrease in height growth values, and increases in 2 produce a decrease in height growth across all heights. An increase in 1k maintains the basic shape of the curve while increasing the height growth. Using the Box-Lucas model for height growth, the height at which the maximum height growth occurs is given as: 147 21 21 max )/ln(2 incH where Hmax inc is the height for which the maximum height growth occurs and all variables are as previously defined. 5.2.4 Potential height growth model 5.2.4.1 Data selection To model potential height growth, trees that were considered as reflecting potential height growth were selected from the untreated plots. Two approaches for selecting potential trees were compared: the binning and the prediction interval approaches. Using the binning approach, data were binned by 5cm DBH classes and alternatively by 5cm DBH classes in combination with site index classes for each species. The first reflected height potential regardless of site quality, whereas the second approach gave different potential heights depending upon site quality. The smallest diameter class was defined as less than 5 cm and the largest as greater than 60 cm. Site index was classified using the BC Ministry of Forests and Range (1999) site index classifications as low (less than 22 m), medium (greater than or equal to 22 and less than 30 m), and high (greater than or equal to 30 m), based on a breast height age of 50 years. Within each bin, trees with height growth greater than the 90 th and then the 99 th percentile were selected for modeling potential height. Using these variations in binning approaches, four subsets of data were selected. For the prediction interval approach, a simple height growth model using height, DBH, and site index as predictor variables was fit using a random coefficient modeling (i.e. parameter prediction) approach (Clutter et al. 1983, pp. 54-56; Schabenberger & 148 Pierce 2002, Chapter 7; Littell et al. 2006, Chapter 8). In the first step, the Box-Lucas model (Eq. [1]) was fit for each tree height growth series with more than three measurement periods using PROC NLIN of SAS version 9.1.3 (SAS Institute Inc. 2004). Estimates of 1k , 1 , and 2 for each tree were then passed to the second step. In the second step, estimated parameters from step 1 were predicted from DBH and transformations of DBH and site index for each species using linear models and PROC REG of SAS version 9.1.3 (SAS Institute Inc. 2004). Two alternatives were selected for parameter prediction in the second step. The first alternative used DBH only: [2] HfHf ff f kH inc 12 21 1 1 expexp ˆ 1f = a0 + a1DBH + a2DBHtransformed 2f = b0 + b1DBH + b2DBHtransformed where DBHtransformed in 1f was defined as DBH 2 for Douglas-fir, and ln(DBH) for western redcedar and western hemlock; and DBHtransformed in 2f was defined as DBH 2 for Douglas-fir and western redcedar, and ln(DBH) for western hemlock. The second alternative used DBH and SI: 1f = c0 + c1DBH + c2DBHtransformed + c3SI 2f = d0 + d1DBH + d2DBHtransformed + d3SI where DBHtransformed in 1f was defined as DBH 2 for Douglas-fir, and ln(DBH) for western redcedar and western hemlock; and DBHtransformed in 2f was defined as DBH 2 for Douglas-fir and western redcedar, and ln(DBH) for western hemlock. The two alternative forms of Eq. [2] were then fit using all trees and PROC MODEL, where the resultant parameter estimates from the second step were used as an initial set of starting 149 values. For all steps in the fitting process, parameter estimates for 1k were restricted to positive values less than 100 to restrict maximum height growth values. Also, initial parameters were varied in order to ensure a global minimum. Finally, using this simple height growth model, all observations with height growth greater than the 90 th and also the 99 th percentiles of prediction intervals were selected as those reflecting potential height increments. This resulted in four subsets of data. For both the binning and predicted value selection approaches, if a tree was selected for any measurement period, the entire measurement sequence for that tree was included. 5.2.4.2 Model fitting Using each subset of data (four subsets using the binning approach and four subsets using the prediction interval approach) by species and a random coefficient modeling approach, the Box-Lucas model for potential annual height growth was fit as follows: [3] HfHf ff f kHPinc 34 43 3 1 expexp ˆ where incP Ĥ is potential annual height growth; 3f and 4f are linear functions using DBH, and site index; and all other variables are as previously defined. As with fitting the simple average height growth model for the prediction interval data selection method, only trees with three or more measures were used in the first step to fit Eq. [1] by tree and in the second step to select the parameter prediction equations. However, all trees in the potential growth subset were used in the final fit of Eq. [3]. 150 To account for possible temporal autocorrelation and heteroskedasticity among tree measures, four different error covariance matrix structures were examined: independent identically distributed (iid); continuous autocorrelated error structure (CAR(x)); a heteroskedastic error structure using height -1 as the weight; and a heteroskedastic, continuous autocorrelated error structure using CAR(x) in combination with height -1 as the weight. The CAR(x) structure has been used in other studies (e.g., Nord-Larson 2006; Crecente-Campo et al. 2009) since it accounts for differing time units between measurements. The error structure for first-order CAR(1) is: itqti hh it vee qtiit )( where ite is the ordinary least squares residual for the t th measure of the i th tree, qtie is the ordinary least squares residual of (t-q) th measure for the i th tree, is the first-order continuous autoregressive parameter, )( qtiit hh is the distance separating the t th from the (t-q) th measure of the i th tree, and itv is an independent, identically normally distributed error term with mean of 0. Eq. [3] was fit for each potential height growth data subset models by species using PROC MODEL of SAS version 9.1.3 (SAS Institute Inc. 2004) under the four alternative error covariance matrix structures. The most appropriate error structure for each fitted potential height growth model was determined using AIC values calculated as: p2ˆ2logL(AIC ) where )ˆlogL( is the log likelihood for the estimated parameter set represented by ˆ , p is the number of parameters in the model, and smaller AIC values indicate a more accurate and/or more parsimonious model. 151 Graphs of predicted versus potential annual height growth averaged by 1 m height classes were produced using each data selection method for each species. The graphs were used to select a final model by species based on how well the model reflected potential height growth. 5.2.5 Height growth model 5.2.5.1 Growth potential dependent method For the growth potential dependent method of modeling height growth, the selected potential annual height growth model (Eq. [3]) was adjusted for average height growth through the inclusion of a multiplicative modifier, following the approach used by Hann et al. (2006): [4] incPinc HMH ˆˆ where incĤ is predicted annual height growth, M is a multiplicative modifier model, and incP Ĥ is the predicted potential annual height growth. Linear and nonlinear model forms were considered for M using possible predictor variables describing tree size and stage of development, site productivity, and inter-tree competition. Tree size and stage of development were represented by DBH (cm) and height (m). To represent site productivity, coastal Douglas-fir site index (SI, m) at a base breast height age of 50 was included. Also, the variable Growth Effective Age (GEA, years) was used to indicate tree productivity, where GEA is the age of a dominant tree of the same height as the tree of interest, for a given site index (Hann and Ritchie 1988). For these data, GEA was predicted from site index equations for the leading tree species in each plot. A number of stand level competition measures were considered, including stems per hectare and basal 152 area per ha. Other measures included Curtis’ (1982) Relative Density, an index measure of density, calculated as: dq G RD Curtis' where dq (cm) is the quadratic mean diameter of all live trees within a plot; and G (m 2 per ha) is the stand basal area per hectare. Stand basal area has been commonly used to measure competition for below ground resources (Fan et al. 2006). In addition, tree-level competition measures, such as the basal area of larger trees (BAL, m 2 per ha) and relative DBH (RDBH), were evaluated for inclusion in the models. BAL was calculated as: n i ii BAtree δ 1 x BAL where δi is an indicator variable (1 if tree i has a DBH greater than the tree of interest, 0 otherwise); BAtreei (m 2 per ha) is the basal area per hectare represented by tree i; and n is the number of all live trees across all species within a plot. RDBH is a ratio of the diameter of the tree of interest to the average diameter within a plot and was calculated as: DBH DBHiRDBH where DBHi (cm) is the DBH of tree i; and DBH (cm) is the average DBH for the plot. After initial model fitting, only linear forms of M were further considered. For each species, linear forms for M were fit using PROC REG of SAS version 9.1.3 (SAS Institute Inc. 2004) and selection was based upon R 2 values. The average annual height growth model was then refit using the selected model form for M using PROC MODEL of SAS version 9.1.3 (SAS Institute Inc. 2004) under the four previously mentioned 153 alternative error covariance matrix structures. The most appropriate error structure for each species was determined using AIC values. 5.2.5.2 Growth potential independent method Using all the untreated data, the random coefficient modeling approach was used to fit the Box-Lucas model [1] for average annual height growth by each species. As with the potential height growth models, only trees with three or more measures were used for the first step in fitting the model by tree, and in the second step in choosing the parameter prediction equations. Initially, simple correlations (r) between 1 , and 2 estimates by species and the possible predictor variables given in Section 5.2.5.1 were used to indicate which predictor variables to use for the parameter prediction equations as follows: [5] HfHf ff f kĤ 1inc 56 65 5 expexp Candidate models for 5f and 6f were then fit by species using PROC MODEL of SAS version 9.1.3 (SAS Institute Inc. 2004). Transformations were used to linearize the relationships where necessary. Although an additional parameter prediction equation for 1k was initially considered, this was not included for reasons of model parsimony and also because very flexible models were obtained with varying 1 and 2 only. The selection of parameter prediction models for 5f and 6f and one of the four aforementioned error structures was based on AIC values calculated using PROC MODEL of SAS version 9.1.3 (SAS Institute Inc. 2004). For all steps in the fitting process, parameter estimates for 1k were restricted to positive values less than 100 to 154 restrict maximum height growth values. Also, initial parameters were varied in order to ensure a global minimum. 5.2.5.3 Comparison of the two methods Graphs of average predicted annual height growth by species and 1 m height classes were produced using both the growth potential independent and dependent methods, along with the average actual annual height growth. The graphs were used to select a model for average annual height growth for each species which was then modified for different treatment effects. 5.2.6 Treatment effects An additional 562 plots received silvicultural treatments; treatments included fertilization (82 plots), thinnings (388 plots), and a combination of fertilization and thinning (92 plots). The majority of fertilized plots received one application of nitrogen ranging in concentration from 50 to 400 kg per ha and a few plots additionally received ammonium phosphate. Of the 388 plots which were thinned, 314 received one thinning and 74 received multiple thinnings. Basal area cut varied from 0.4 to 79.9 m 2 per ha with an average of 7.6 m 2 per ha, and stems per ha cut varied from 1.24 to 10,360 trees per ha, with average of 988 trees per ha. It should be noted that the data used in this study primarily represents silvicultural treatments applied to operationally managed stands rather than treatments applied in an experimental setting. Therefore, the results in this study include the interactions that are often removed from experiments through careful selection of experimental plots. 155 5.2.6.1 Fertilization The selected Box-Lucas model developed for average annual height growth was modified to include the effects of fertilization. An additional fertilization effect variable, defined as: Ftime F effectfert. where Ftime is the time (years) since fertilization, and F is an indicator variable, defined as 1 for fertilized and 0 otherwise, was added to 5f and 6f , which were subsequently relabeled as 7f and 8f . The categorical variable F identifies measurements taken following fertilizer application. This approach was used by Carlson et al. (2008) in modeling the mean relative growth rate of loblolly pine (Pinus taeda L.) and by Huber et al. (2009) for other growth models. The parameter 1k was also replaced with a parameter prediction model for fertilized plots, as follows: [6] HfHf ff f fH inc 78 87 7 9 expexp ˆ where 9f is a function of the fertilizer effect variable. Since modifying all of the parameter prediction models might not have been needed, all possible combinations of modified versus not modified versions of 7f , 8f , and 9f were fit by species using PROC MODEL of SAS version 9.1.3 (SAS Institute Inc. 2004) and plots that had been fertilized. For each species, Eq. [6] was fit using each of the four error structures and AIC values used to determine the most appropriate model and error structure. 156 The fertilized data were subset into each combination of two site quality (low and high) and two density (low and high) classes, based on the measurement when fertilization occurred. For comparison to unthinned stands, the untreated data were similarly divided into site and density classes. Site quality was defined as low for stands where site index was less than or equal to 22, else it was defined as high. Density was defined as low for plots with a quadratic mean diameter value less than or equal to 20, else it was defined as high. For each subset and 2 m height class, average predicted values were calculated for: 1) predicted annual height growth for fertilized plots using Eq. [6] and fertilized plot data; 2) predicted annual height growth using equation [5] and untreated plot data; and 3) predicted potential height growth using Eq. [3] and untreated data. For each species, site quality and density combination, a graph illustrating these values was created. 5.2.6.2 Thinning Thinning has often been modeled through the inclusion of thinning parameters designed to describe the type of treatment received. For example, Wimberly and Bare (1996) added a parameter to their basal area growth model defined as the amount of basal area removed. In this study, no modifications were made to the selected Box-Lucas model developed for average height growth to include thinning. Instead, the selected Box-Lucas model developed for average height growth (Eq. [5]) was used, with the state variables changed instantaneously after thinning as input variables to predict growth following thinning. In particular, state variables which changed as a result of the thinning were basal area per ha, basal area of larger trees, and stems per ha. This approach was used for young Scots pine (Pinus sylvestris L.) stands by Pukkala et al. 157 (2002), who found that models without additional thinning variables or parameters can produce accurate unbiased predictions in thinned stands. As with fertilized plots, a model-based approach was used to indicate thinning impacts. The data following thinning were subset into site quality and density classes as with the fertilized data. Average predicted values using Eq. [5] and the thinned data to represent average growth after thinning, using Eq. [5] and untreated data to represent average growth without thinning, and, finally using Eq. [3] with untreated data to represent potential growth were calculated by 2 m height classes and graphed. 5.2.6.3 Fertilization and thinning For plots that received a combination of thinning and fertilizer treatments, the Box- Lucas model modified for fertilized plots was applied directly. As for the plots which only received a thinning, some state variables for these plots immediately change following thinning. This approach assumed a greater response would be seen for plots which received both thinning and fertilization than for those which received only one treatment. As with the fertilized and thinned only treatments, a model-based approach was used to examine effects of thinning and fertilizing. These data were subset into site quality and density classes, and graphs of average predicted values for treated, untreated, and growth potential were produced. 5.2.7 Model accuracy To assess the accuracy of all height growth models, mean bias and root mean-squared error (RMSE) values by 5 m height classes were calculated for annual height growth by species and treatment type. Mean bias was calculated as: 158 m yy m 1i ii ˆ BiasMean where iy is the actual annual height growth for tree i, iŷ is the predicted annual height growth for tree j, and m is the total number of measures for all trees in the height class. Root mean-squared error was calculated as: pm yy m 1i 2 ii ˆ RMSE where p is the number of parameters in the model and all other variables are as previously defined. 5.3 Results 5.3.1 Potential height growth model Once the eight subsets of the untreated data were determined, the potential annual height growth models (Eq [3]) were fitted for each subset and species. Variables selected for 3f were DBH and DBH 0.5 for Douglas-fir, and DBH and ln(DBH) for western redcedar and western hemlock using the binning by DBH or by DBH and site index. Variables selected for 3f were DBH and DBH 2 for Douglas-fir, and DBH and ln(DBH) for western redcedar and western hemlock using the prediction intervals by DBH or by DBH and site index. For 4f , DBH and DBH 2 were selected for Douglas-fir and for western redcedar, and DBH and ln(DBH) were selected for western hemlock, regardless of the method used to select the data. For all models, no advantage was gained in using site index for 3f and 4f nor in using a more complex error structure over an iid error structure, based on AIC. 159 Using the potential height growth models fit using the data selected using the prediction interval approach and the simple model using DBH only, the upper boundaries using the 90% and 99% prediction intervals were necessarily parallel, with the 99% model better outlining the upper edge of observed height growth measures across height classes for all species (Figure 5-2). Adding site index into the simple model decreased the potential height growth curves (Figure 5-3), since the upper limit over all is no longer represented. Using the binning selection method, a more irregular limiting curve resulted, particularly when DBH alone was used in binning (Figures 5-2 versus 5-3) and also when the top 1% was used rather than the top 10%. For all methods, a more irregular limiting curve resulted for larger height classes with less data. Because the prediction intervals method using the upper value of the 99% prediction interval of the simple model with DBH only best defined the upper boundary for potential annual height growth, the model fit using only those data was selected as the best for estimating potential height growth. 5.3.2 Height growth model 5.3.2.1 Growth potential dependent model The multiplicative modifier models were: For Douglas-fir: CurtisRDcRDBHcGEAcSIcBALcGccM 7654320 For western redcedar: CurtisRDcSIcBALcGcDBHccM 743210 For western hemlock: CurtisRDcRDBHcGEAcSIcBALcGcDBHccM 76543210 160 where c0 to c7 represent parameters that were estimated by species. For all models, no advantage was gained changing the error structure from an iid error structure based on AIC. 5.3.2.2 Growth potential independent model In the first step of the parameter prediction approach, the Box-Lucas model was fit for each tree in the untreated plots having more than three measures. Initial estimates for 1 and 2 were restricted to positive values less than 10; without this restriction the values tended toward infinity. The absolute values of all correlations between estimates of 1 by tree and all possible predictor variables were greater than 0.05 for Douglas-fir and for western redcedar, except for GEA in western redcedar (Table 5-3). Since all the absolute values of all correlations for 1 were less than 0.05 for western hemlock, all possible combinations of predictor variables were investigated. The absolute values of all correlations between 2 and all possible predictor variables were greater than 0.05 for all species, except for GEA for Douglas-fir and site index for western hemlock. Based on the near normal distribution for the log transformation of the predicted parameters, an exponential equation form was selected for each parameter prediction equation. The selected parameter prediction equations by species were: For Douglas-fir: )exp( 86543105 SPHaRDBHaGEAaSIaBALaDBHaaf )exp( 643206 RDBHbSIbBALbGbbf For western redcedar: )exp( 876432105 SPHaCurtisRDaRDBHaSIaBALaGaDBHaaf 161 )exp( 86506 SPHbRDBHbGEAbbf For western hemlock: )exp( 87643205 SPHaCurtisRDaRDBHaSIaBALaGaaf )exp( 8765106 SPHbCurtisRDbRDBHbGEAbDBHbbf where SPH is stems per ha, and a0 to a8, and b0 to b8 are parameters to be estimated. For all three species, no gain was obtained in using an alternative to the iid error structure (Table 5-4). 5.3.2.3 Comparison of the two methods Both the growth potential dependent and independent methods used to obtain predicted annual height growth resulted in averages of predicted values similar to average actual values (Figure 5-4), with the exception of large trees (ie, trees taller than 55 m for Douglas-fir and taller than 40 m for western redcedar) where less data were available. The growth potential independent method was selected as the model for untreated stands and was modified for the effects fertilization and thinning. For the selected height growth model using all untreated data, the mean bias values by 5 m height class for Douglas-fir were between -0.60 and 0.01 m, with larger negative values found for the smallest and the largest diameter classes (Figure 5-5). Similar results were found for western redcedar. Western hemlock returned small negative mean bias values for smaller trees only. The fewest observations were found in the small height classes, less than 5 m, as only trees greater than 4 cm DBH were included in the study. The RMSE values were similar all three species with larger values found in the small diameter classes (Figure 5-6). 162 5.3.3 Treatment effects 5.3.3.1 Fertilization Only slight changes were seen in AIC values for Douglas-fir, western redcedar, and western hemlock between Eq [5], the original untreated height growth model, and Eq [6], the model modified to include the fertilizer effect variable (Table 5-5). For Douglas- fir, the best AIC value was found when 5f and k1 were modified to include the fertilizer effect variable as 7f and 9f . No further improvement was found when 6f was modified to include the fertilizer effect variable for Douglas-fir. For western redcedar, the highest AIC was found when 6f was modified to include the fertilizer effect variable as 8f . The best AIC value found for western hemlock included the modification of 6f and k1 to include the fertilizer effect variable as 8f and 9f . For Douglas-fir, the Box-Lucas model modified to include the fertilizer effect variable (Eq. [6]) was: Ftime Faff 957 exp For western redcedar and western hemlock: Ftime Fbff 968 exp For western redcedar: 19 kf and for Douglas-fir and western hemlock: Ftime F ckf 19 163 For the fertilized plots, the mean bias values for all species were nearly zero (Figure 5-5). RMSE values were also very small (Figure 5-6). It should be noted that the mean bias and RMSE values for untreated plots were based on more trees. For Douglas-fir, the fertilized plots had lower predicted height growth values compared to the untreated predicted values, except for the high site quality and high density class, where fertilized plots had larger predicted height growth values for most height classes (Figure 5-7). For western hemlock, height growth appeared to respond positively to fertilizer on high density sites where the fertilized plots had higher average annual height growth values compared to the untreated plots (Figure 5-8 b and d), with perhaps a higher response on lower site quality stands (Figure 5-8 b). Potential height was always considerably higher than fertilized plots achieved. The number of plots for western redcedar in each site and density classes were too small to use the model-based approach to assess fertilization impacts on potential height growth relative to no fertilization and averages are not displayed. Also, there was little information available for fertilized Douglas-fir for low site quality and high density. 5.3.3.2 Thinning Because the state variables basal area per ha, basal area per ha of larger trees, relative diameter, Curtis RD, and stems per ha change immediately following thinning, the untreated model (Eq. [5]) was also used for thinned plots as previously noted. Using Eq. [5] and the plot data following thinning, the mean bias values were nearly zero for all height classes (Figure 5-5), and RMSE values were very low (Figure 5-6) for all three species. 164 Based on the model, thinned plots had a higher average annual predicted height growth compared to the untreated plots for Douglas-fir for shorter heights, particularly on higher productivity sites (Figure 5-9). There was little data available for thinned western redcedar or western hemlock on stands which were defined as low site quality and high density. For other site and density classes, results for western redcedar (Figure 5-10) and western hemlock (Figure 5-11) were similar, with the largest difference between thinned and untreated plots seen for those located on high productivity sites with higher densities. As with fertilization, potential height was always considerably higher than thinned plots achieved. 5.3.3.3 Fertilization and thinning When plots were both fertilized and thinned, Douglas-fir located on low site quality sites had a lower average annual height growth compared to the untreated plots (Figure 5- 12 a). For treated Douglas-fir, on average, the average annual height growth was greater on sites with higher density, regardless of density levels, compared to the untreated plots (Figure 5-12 b and d). For low quality and high density level sites, the average annual height growth for treated Douglas-fir was similar to the untreated trees (Figure 5-12 c). For low density plots, on average, the average annual height growth for treated western hemlock was lower for untreated trees (Figure 5-13 a and b), while on high density plots, it was larger (Figure 5-13 c). Again, potential height growth was considerably higher than the predicted values for treated plots. Mean bias and RMSE values for plots which received fertilization following thinning were greatest for western hemlock (Figures 5-5 and 5-6). The numbers of plots for western redcedar were too little to assess the impact of fertilization in combination 165 with thinning relative to no treatment on potential height growth and averages are not displayed. Also, there was little information available for fertilization in combination with thinning for western hemlock in low site quality and high density stands. 5.4 Discussion Dominant and codominant trees have often been used to define potential height, with dominant trees defined as those with the largest diameter (eg., Wang et al. 2008; Weiskittel et al. 2009). Height growth of dominant trees has been shown to be relatively independent of stand density (Goelz and Burk 1998). However, for young stands of shade-intolerant Douglas-fir, increasing levels of stand density have been found to increase height growth (Scott et al. 1998). If potential height growth is defined as the theoretical maxima across all stand conditions, then the definition of potential height growth as being growth of the dominant trees may not be applicable in mixed-species stands where vertical stand structure layers according to shade tolerance and the status of dominant trees changes with time (Oliver and Larson 1996). The approaches used to select potential height growth trees gave very different results as hypothesized. When potential height growth trees were defined using the binning method and DBH class alone, results varied greatly between the upper 10% versus 1%, with the greatest difference found for the shade-tolerant species western hemlock. Results were more similar when trees were binned by DBH and site index class. Selecting trees from bins using DBH classes alone represent trees of maximum growth across all site conditions. Using upper percentiles of prediction intervals for a simple growth model resulted in a smoother trend over height classes. Regardless of how potential height growth trees are defined, an assumption was made that these sample trees 166 selected as potential height growth trees have attained the maximum height growth rate possible (Hann and Ritchie 1988). Using the prediction interval percentile approach and the 99% percentile better reflected the rarely achieved maximum height growth (i.e., the boundary of height growth) and therefore best defined this theoretical maximum potential. As was hypothesized, little difference between the growth potential dependent and independent methods for height growth models existed. Both methods returned predicted values close to the actual annual height growth for all species. Using the growth potential dependent approach for height growth, often the modifiers use a multiplicative function, reducing the potential height. This may not be sufficiently flexible to reflect height growth, and will depend greatly upon the approach used to model potential height growth. In this study, results were similar for the two approaches likely because the models used for potential height growth and for height growth were similar in functional form. Consequently, modifying the potential height growth model for height growth gave similar predicted values as modeling height growth directly. Because the definition of potential height growth does impact the model of potential height, separate models for potential height growth and height growth are suggested. Increases in Douglas-fir average height increment on fertilized plots were dependent upon stand density, with higher density increasing average height increment and lower density decreasing average height increment. This would suggest that an increase in below ground nutrients without an increase in light does not necessarily increase average height increment for Douglas-fir, a shade intolerant species. In low density stands, an increase in nutrients may increase diameter growth rather than height 167 growth. For Douglas-fir in the Pacific Northwest, Footen et al. (2009) found that five years following fertilization, the increase in DBH was twice that of the increase found for height. For western hemlock, a shade-tolerant species, an increase in average height growth was seen regardless of stand quality or stand density. For western hemlock and western redcedar, fertilization has been show to increase height growth (e.g., Bennett et al. 2003; Devine and Harrington 2009). Impacts of thinning on height growth were addressed using the same parameter estimates as the untreated plots, allowing only the state variables to change. Since the main impact of thinning is the reduction of competition, this response is reflected in the post-thinning state variables such as basal area per ha, stems per ha, and basal area per ha of larger trees. As was hypothesized, the model reflects the post-thinning response well for Douglas-fir. For western redcedar, the model underestimates height growth for all trees, and for western hemlock the model underestimates for trees less than 40 cm in diameter, otherwise it overestimates. The instantaneous change in competition for light resources which occurs after thinning was reflected in the model by the change in basal area per ha, basal area of larger trees, and stems per ha. Unlike the stated hypothesis, an increase in average height growth was seen for all three species in the thinned plots compared to unthinned plots. Renninger et al. (2007) found that the relative height growth of released versus suppressed trees was higher in Douglas-fir than in western hemlock (65%). It has also been shown that relative neighbor size had a negative effect on height growth for Douglas-fir and western hemlock (Erickson et al. 2009). Western hemlock is a shade-tolerant species that frequently regenerates naturally under a canopy of mature trees. Suppression and release 168 is more common in this species (Franklin et al. 2002) and residual trees are able to grow rapidly following a thinning (Deal 2001). Devine and Harrington (2009) found that western redcedar seemed to show an increase in height following thinning. The largest increase in height growth was found when fertilization and thinning was done in combination, as was hypothesized. This increase was found for Douglas-fir and western hemlock on sites which had a combination of high site quality and high density. However, on average, the potential height growth was not achieved through thinning, fertilization, or through the combination of thinning and fertilizing. Other reasons for greater height growth as represented by potential height growth curves include differences in genetics and local competition and site impacts. As noted previously, the data used in this study primarily represent silvicultural treatments applied to managed stands rather than treatments applied in an experimental setting. As a result, the impacts of fertilization and thinning are more variable than in studies based on experiments. The intent was to obtain height growth models as components of a growth model for use on these stands, and to use a model-based approach to assess treatment effects in an operational context in these mixed-species stands. 5.5 Conclusion Since both potential height growth and height growth forecasts are useful in managing stands, both variables were modeled for three species growing in mixed- species stands of Vancouver Island, BC. The Box and Lucas model was used an effective basis for these models. Several approaches were used to define potential height growth. The use of prediction intervals from a simple height growth model is 169 recommended, since observations with height growth measures greater than an upper percentile appear to reflect growth potential and result in a smooth trend line over height. Both the growth potential dependent and independent methods were used to determine average height growth. Modeling average height growth and potential height growth should be done using two separate equations rather than using the growth potential dependent method, as the definition of potential height is subjective and can produce varying results. Models were modified for growth in thinned and fertilized stands, and this is seen as a possibly useful way of assessing treatment impacts in a management setting. Using this model-based approach, an increase in average height growth was observed for thinned Douglas-fir no matter the site or starting density, whereas an increase was only seen on high density stands when fertilized, and on high site quality and high density stands when fertilization and thinning occurred in combination. Little information was available for western redcedar and further research needs to be done to determine how height growth is affected by fertilization for this species. 170 Figure 5-1. Box-Lucas model for varying 1 and 2 values where k1=1; all forms indicate an increase in height growth to a maximum height growth followed by a decrease to a lower limit. 171 Figure 5-2. Average annual potential height by 1 m height classes for a) Douglas-fir, b) western redcedar, and c) western hemlock, using DBH and height as estimator variables. 172 Figure 5-3. Average annual potential height by 1 m height classes for a) Douglas-fir, b) western redcedar, and c) western hemlock, using DBH, site index, and height as estimator variables. 173 Figure 5-4. Average annual average height by 1 m height classes using the growth potential independent and dependent methods by a) Douglas-fir, b) western redcedar, and c) western hemlock, using DBH and height as estimator variables. All methods are quite similar and difficult to distinguish except for the taller trees. 174 Figure 5-5. Mean bias (cm) values by 5 m height classes for: a) Douglas-fir; b) western redcedar; and c) western hemlock, by treatment. 175 Figure 5-6. RMSE (cm) values by 5 m height classes for: a) Douglas-fir; b) western redcedar; and c) western hemlock, by treatment. 176 Figure 5-7. Average annual mean untreated, fertilized, and potential height growth by 2 m height classes for Douglas-fir in stands with a) low site quality and low density, b) high site quality and low density, and c) high site quality and high density. The numbers of trees used to determine the fertilized values are given across the top of each graph. 177 Figure 5-8. Average annual mean untreated, fertilized, and potential height growth by 2 m height classes for western hemlock in stands with a) low site quality and low density, b) low site quality and high density, c) high site quality and low density, and d) high site quality and high density. The numbers of trees used to determine fertilized values are given across the top of each graph. 178 Figure 5-9. Average annual mean untreated, thinned, and potential height growth by 2 m height classes for Douglas-fir in stands with a) low site quality and low density, b) low site quality and high density, c) high site quality and low density, and d) high site quality and high density. The numbers of trees used to determine thinned values are given across the top of each graph. 179 Figure 5-10. Average annual mean untreated, thinned, and potential height growth by 2 m height classes for western redcedar in stands with a) low site quality and low density, b) high site quality and low density, and c) high site quality and high density. The numbers of trees used to determine thinned values are given across the top of each graph. 180 Figure 5-11. Average annual mean untreated, thinned, and potential height growth by 2 m height classes for western hemlock in stands with a) low site quality and low density, b) high site quality and low density, and c) high site quality and high density. The numbers of trees used to determine thinned values are given across the top of each graph. 181 Figure 5-12. Average annual mean untreated, fertilized and thinned, and potential height growth by 2 m height classes for Douglas-fir in stands with a) low site quality and low density, b) low site quality and high density, c) high site quality and low density, and d) high site quality and high density. The numbers of trees used to determine fertilized & thinned values are given across the top of each graph. 182 Figure 5-13. Average annual mean untreated, fertilized and thinned, and potential height growth by 2 m height classes for western hemlock in stands with a) low site quality and low density, b) high site quality and low density, and c) high site quality and high density. The numbers of trees used to determine fertilized & thinned values are given across the top of each graph. 183 Table 5-1. Plot summary statistics at plot establishment and at time of treatment (n=number of plots). Treatment Type Variable a Plot Establishment Time of Treatment Mean Min Max Mean Min Max Untreated n=1,316 SPH 2,073 222 11,750 G (m 2 per ha) 38.1 0.1 185.0 SI (m) 29.6 10.7 52.8 Curtis' RD 8.7 0.1 24.2 Fertilized n=82 SPH 2,430 370 8,123 2,274 370 8,123 G (m 2 per ha) 35.2 1.4 96.1 38.5 1.4 96.1 SI (m) 31.8 11.2 42.3 32.1 11.2 42.3 Curtis' RD 8.7 0.7 18.8 9.2 0.7 17.9 Thinned n=388 SPH 2,180 309 11,370 2,172 309 11,370 G (m 2 per ha) 29.0 0.5 110.5 31.0 1.0 110.5 SI (m) 31.5 9.2 50.0 31.6 9.2 50.0 Curtis' RD 7.2 0.3 21.1 7.5 0.5 21.1 Multiple Treatments n=92 SPH 2,351 530 7,802 2,314 530 7,605 G (m 2 per ha) 45.1 5.0 91.9 47.4 5.0 91.9 SI (m) 32.7 8.9 42.9 32.7 8.9 42.9 Curtis' RD 10.7 2.2 17.4 11.1 2.2 17.4 a SPH is stems per ha, G is basal area per ha, SI is Site Index, and Curtis’ RD is Curtis’ Relative Density, Min is minimum, and Max is maximum. 184 Table 5-2. Means and standard deviations (in brackets) for tree-level variables at plot establishment by silvicultural treatment and species. Species Variable a Untreated Fertilized Thinned Multiple Treatments Douglas- fir DBH (cm) 17.1 (13.2) 15.8 (9.5) 19.5 (14.1) 24.6 (11.3) Height (m) 13.7 (9.5) 13.5 (7.0) 15.9 (10.5) 19.7 (7.6) BAL (m 2 per ha) 8.8 (11.7) 10.7 (10.6) 9.7 (11.4) 11.1 (9.5) RDBH 1.5 (0.5) 1.4 (0.5) 1.4 (0.4) 1.7 (0.5) GEA (years) 20.2 (20.5) 18.3 (11.7) 23.0 (21.3) 25.6 (14.7) Hinc (cm) 0.59 (0.28) 0.66 (0.24) 0.58 (0.31) 0.61 (0.26) Western redcedar DBH (cm) 19.9 (14.2) 17.6 (11.6) 17.8 (14.6) 28.9 (13.2) Height (m) 13.9 (8.6) 13.3 (9.6) 12.7 (8.4) 20.7 (6.6) BAL (m 2 per ha) 17.5 (19.8) 14.2 (15.4) 16.6 (21.1) 12.5 (11.0) RDBH 1.6 (0.6) 1.6 (0.6) 1.5 (0.6) 2.0 (0.7) GEA (years) 26.8 (19.7) 18.7 (14.7) 20.1 (15.0) 33.2 (19.7) Hinc (cm) 0.36 (0.20) 0.45 (0.18) 0.44 (0.23) 0.60 (0.28) Western hemlock DBH (cm) 22.9 (14.9) 18.3 (11.2) 13.2 (9.1) 26.5 (12.5) Height (m) 19.8 (12.0) 15.3 (9.6) 10.5 (7.3) 21.9 (9.5) BAL (m 2 per ha) 20.1 (20.7) 12.2 (16.5) 8.4 (12.2) 16.0 (14.9) RDBH 1.5 (0.5) 1.5 (0.5) 1.5 (0.5) 1.7 (0.5) GEA (years) 32.2 (25.6) 20.9 (16.9) 14.3 (12.7) 31.1 (15.6) Hinc (cm) 0.46 (0.26) 0.59 (0.25) 0.62 (0.26) 0.51 (0.27) a DBH is diameter at breast height, BAL is basal area per ha of larger trees, RDBH is relative diameter at breast height, GEA is growth effective age, and Hinc is the annual height growth. 185 Table 5-3. Correlations between possible estimator variables and predicted 1 or 2 for trees with more than three measurements (n=number of trees). Variables a Douglas-fir (n=4,050) Western redcedar (n=307) Western hemlock (n=1,402) θ1 θ2 θ1 θ2 θ1 θ2 DBH -0.054 -0.209 -0.056 -0.244 0.017 -0.456 G 0.100 -0.197 0.149 -0.196 0.003 -0.439 BAL 0.212 -0.087 0.312 -0.114 -0.021 -0.220 SI -0.143 0.005 -0.085 0.061 0.012 -0.035 GEA 0.079 -0.228 0.041 -0.226 -0.019 -0.428 RDBH -0.135 -0.109 -0.231 -0.102 0.039 -0.090 Curtis’ RD 0.107 -0.158 0.122 -0.094 0.017 -0.342 SPH 0.054 0.130 -0.124 0.275 -0.017 0.279 a DBH is diameter at breast height, G is basal area per ha, BAL is basal area per ha of larger trees, SI is Site Index, GEA is growth effective age, RDBH is relative diameter at breast height, Curtis’ RD is Curtis’ Relative Density, and SPH is stems per ha. 186 Table 5-4. AIC values for augmented Box-Lucas model (Eq. 4) with varying error structures for trees on untreated plots (n=numbers of trees). Error Structure Douglas-fir (n=10,069) Western redcedar (n=2,016) Western hemlock (n=8,059) No modification -32,218 -4,146 -17,044 CAR(1) -24,762 -3,446 -9,866 Weighted using DBH -1 -5,320 -3,594 -10,690 CAR(1) and Weighted using DBH -1 -31,078 -3,980 -16,204 187 Table 5-5. AIC values for augmented Box-Lucas model modified to include the fertilizer variable for fertilized plots (n= number of plots). Modification AIC Value Douglas-fir (n=674) Western redcedar (n=33) Western hemlock (n=629) f7 -1,491 -43 -487 f8 -1,482 -44 -486 f 7 and f 8 -1,497 -42 -485 f 9 -1,487 -43 -487 f 7 and f 9 -1,518 -42 -485 f 8 and f 9 -1,490 -42 -492 f 7, f 8, and f 9 -1,497 -40 -490 No modification -1,472 -42 -486 188 5.6 References B.C. Ministry of Forests and Range. 1999. Timber supply review: Strathcona timber supply area analysis report. Province of British Columbia Ministry of Forests. Victoria, BC. 91 pp. Bennett, J.N., Blevins, L.L., Barker, J.E., Blevins, D.P., and Prescott, C.E. 2003. Increase in tree growth and nutrient supply still apparent 10 to 13 years following fertilization and vegetation control of salal-dominated cedar-hemlock stands on Vancouver Island. Can. J. For. Res. 33(8): 1516-1524. Beese, W.J., Dunsworth, B.G., Zielke, K., and Bancroft, B. 2003. Maintaining attributes of old-growth forests in coastal B.C. through variable retention. For. Chron. 79(3) 570-578. Box, G.E.P. and Lucas, H.L. 1959. Design of experiments in non-linear situations. Biometrika 46(1/2): 77-90. Carlson, C.A., Burkhart, H.E., Allen, H.L., and Fox, T.R. 2008. Absolute and relative changes in tree growth rates and changes to the stand diameter distribution of Pinus taeda as a result of midrotation fertilizer applications. Can. J. For. Res. 38(7): 2063- 2071. Chen, H.Y.H. and Klinka, K. 2003. Aboveground productivity of western hemlock and western redcedar mixed-species stands in southern coastal British Columbia. For. Ecol. Manage. 184(1-3): 55-64. Clutter, J.L., Fortson, J.C., Pienaar, L.V., Brister, G.H., Bailey, R.L. 1983. Timber Management – A Quantitative Approach. John Wiley & Sons, New York. pp.333 189 Crecente-Campo, F., Marshall, P., LeMay, V., and Dieguez-Aranda, U. 2009. A crown profile model for Pinus radiata D. Don in northwestern Spain. For. Ecol. Manage. 257(12): 2370-2379. Curtis, R.O. 1982. A simple index of stand density for Douglas-fir. For. Sci. 28(1):92-94. Deal, R.L. 2001. The effects of partial cutting on forest plant communities of western hemlock- Sitka spruce stands in southeast Alaska. Can. J. For. Res. 31(12): 2067- 2079. Devine, W.D. and Harrington, C.A. 2009. Western redcedar response to precommercial thinning and fertilization through 25 years posttreatment. Can. J. For. Res. 39(3): 619-628. Drever, C.R. and Lertzman, K.P. 2001. Light-growth responses of coastal Douglas-fir and western redcedar saplingsunder different regimes of soil moisture and nutrients. Can. J. For. Res. 31(12): 2124-2133. Ek, A.R., Monserud, R.A., 1974. Trials with program FOREST: growth and reproduction simulation for mixed species even- or uneven-aged forest stands. In: Fries, J. (Ed.), Growth Models for Tree and Stand Simulation. Research Notes. Royal College of Forestry, Skogshögskolan, Sweden, pp. 56–73. Erickson, H.E., Harrington, C.A., and Marshall, D.D. 2009. Tree growth at stand and individual scales in two dual-species mixture experiments in southern Washington State, USA. Can. J. For. Res. 39(6): 1119-1132. Fan, Z., Kabrick, J.M., and Shifley, S.R. 2006. Classification and regression tree based survival analysis in oak-dominated forests of Missouri’s Ozark highlands. Can. J. For. Res. 36(7): 1740-1748. 190 Fitchner, K. and Schultz, E.D. 1992. The effects of nitrogen nutrition on growth and biomass partitioning of annual plants originating from habitats of different nitrogen availability. Oecologia 92(2): 236-241. Footen, P.W., Harrison, R.B., and Strahm, B.D. 2009. Long-term effects of nitrogen fertilization on the productivity of subsequent stands of Douglas-fir in the Pacific Northwest. For. Ecol. Manage. 258(10): 2194-2198. Franklin, J.F., Spies, T.A., Van Pelt, R., Carey, A.B., Thornburgh, D.A., Berg, D.R., Lindenmayer, D. B., Harmon, M.E., Keeten, W.S., Shaw, D.C., Bible, K., and Chen, J. 2002. Disturbances and structural development of natural forest ecosystems with silvicultural implications, using Douglas-fir forests as an example. For. Ecol. Manage. 155(1-3): 399-423. Gavin, D.G., Brubaker, L., and Lertzman, K.P. 2003. Holocene fire history of a coastal temperate rain forest based on soil charcoal radiocarbon dates. Ecology 84(1) 186- 204. Getzin, S., Dean, C., He, F., Trofymow, J.A., Wiegand, K. and Wiegand, T. 2006. Spatial patterns and competition of tree species in a Douglas-fir chronosequence on Vancouver Island. Ecography 29(5): 671-682. Goelz, J.C.G. and Burk, T.E. 1998. Long-term trends in height growth of jack pine in north central Ontario. For. Sci. 44(1): 158-164. Grime, J.P., and Hunt, R. 1975. Relative growth rate: its range and adaptive significance in a local flora. J. Ecol. 63(2): 393-422. Hann, D.W. and Ritchie, M.W. 1986. Development of a tree height growth model for Douglas-fir. For. Ecol. Manage. 15 (2): 135-145. 191 Hann, D.W. and Ritchie, M.W. 1988. Height growth rate of Douglas-fir: a comparison of model forms. For. Sci. 34(1): 165-175. Hann, D.W., Marshall, D.D., and Hanus, M.L. 2006. Reanalysis of the ORGANON equations for diameter-growth rate, height growth rate, and mortality rate of Douglas- fir. Res. Contribution 49. Oregon State Univ. College of Forestry, Corvallis, OR. 24 pp. Huber, M., Halmschlager, E., and Sterba, H. 2009. The impact of forest fertilization on growth of mature Norway spruce affected by Sirococcus shoot blight. For. Ecol. Manage. 257(6): 1489-1495. Ishii, H., Reynolds, J.H., Ford, D.E., and Shaw, D.C. 2000. Height growth and vertical development of an old-growth Pseudotsuga-Tsuga forest in southwestern Washington State, USA. Can. J. For. Res. 30(1): 17-24. Kobe, R.K. and Coates, K.D. 1997. Models of sapling mortality as a function of growth to characterize interspecific variation in shade tolerance of eight tree species of northwestern British Columbia. Can. J. For. Res. 27(2): 227-236. Lertzman, K.P., Sutherland, G.D., Inselberg, A. and Saunders, S.C. 1996. Canopy gaps and the landscape mosaic in a coastal temperate rain forest. Ecology 77(4): 1254- 1270. Littell, R.C., Milliken, G.A., Stroup, W.W., Wolfinger, R.D., and Schabenberger, O. 2006. SAS for Mixed Models, 2nd ed. SAS Institute Inc. Cary, NC. Meidinger, D. and Pojar, J., eds. 1991. Ecosystems of British Columbia. The Research Branch, BC Ministry of Forests, Victoria, BC. Pp. 95-111. 192 Nord-Larson, T. 2006. Modeling individual-tree growth from data with highly irregular measurement intervals. For. Sci. 52(2): 198-208. Oliver, C.D. and Larson, B.C. 1996. Forest stand dynamics. McGraw-Hill, New York. Palahi, M., Pukkala, T., Miina, J., and Montero, G. 2003. Individual-tree growth and mortality models for Scots pine (Pinus sylvestris L.) in north-east Spain. Ann. For. Sci. 60(1): 1-10. Piotto, D. 2008. A meta-analysis comparing tree growth in monocultures and mixed plantations. For. Ecol. Manage. 255(3-4): 781-786. Pukkala, T., Miina, J., and Palahi, M. 2002. Thinning response and thinning bias in a young scots pine stand. Silvica Fennica 36(4): 827-840. Rathbun, L.C., LeMay, V., and Smith, N. 2010. Diameter growth models for mixed- species stands of Coastal British Columbia including thinning and fertilization effects. Submitted for possible publication. Renninger, H.J., Meinzer, F.C., and Gartner, B.L. 2007. Hydraulic architecture and photosynthetic capacity as contstraints on release from suppression in Douglas-fir and western hemlock. Tree Physiol. 27(1): 33-42. SAS Institute Inc. 2004. SAS/STAT 9.1 user’s guide. SAS Institute Inc., Cary, NC. Schabenberger, O., and Pierce, F.J. 2002. Contemporary Statistical Models for the Plant and Soil Sciences. CRC Press. Boca Raton, FL. Scott, W., Meade, R., Leon, R., Hyink, D., and Miller, R. 1998. Planting density and tree-size relations in coast Douglas-fir. Can. J. For. Res. 28(1): 74-78. Temesgen, H., Monleon, V.J., and Hann, D.W. 2008. Analysis and comparison of nonlinear tree height estimation. Can. J. For. Res. 38(3): 553-565. 193 Walters, M.B., Kruger, E.L., and Reich, P.B. 1993. Growth, biomass distribution and CO2 exchange of northern hardwood seedlings in high and low light: relationships with successional status and shade tolerance. Oecologia 94(1): 7-16. Wang, M., Rennolls, K., and Tang, S. 2008. Bivariate distribution modeling of tree diameters and heights: dependency modeling using copulas. For. Sci. 54(3): 284-293. Weiskittel, A.R., Hann, D.W., Hibbs, D.E., Lam, Y.L., and Bluhm, A.A. 2009. Modeling top height growth of red alder plantations. For. Ecol. Manage. 258(3): 323- 331. Wimberly, M.C. and Bare, B.B. 1996. Distance-dependent and distance-independent models of Douglas-fir and western hemlock basal area growth following silvicultural treatment. For. Ecol. Manage. 89(3): 1-11. York, R.A., Battles, J.J., and Heald, R.C. 2003. Edge effects in mixed conifer group selection openings: Tree height response to resource gradients. For. Ecol. Manage 179(1-3): 107-121. 194 6. Conclusions 6.1 Research conclusions In the coastal temperate rain forests of British Columbia (BC), variable retention of live trees in harvest units is a relatively recent forest management practice (Bunnell and Dunsworth 2009, Peterson and Anderson 2009). As a result, few studies exist on the impacts on the growth and mortality of residual trees following harvest (Aubry et al. 2009). Consequently, models are needed to simulate alternative variable retention scenarios and evaluate outcomes. In addition, the successional process of coastal BC forests is driven by gap dynamics, creating complex mixed-species stands. The aim of this study was to develop growth and mortality models across a wide spectrum of silvicultural treatments in order to forecast growth and to use a model-based approach to evaluate treatment effects for Douglas-fir, western redcedar, and western hemlock. 6.1.1 Mortality An understanding of the factors affecting mortality should be considered when developing survival/mortality models; this is especially true for mixed-species stands with complex structural diversity. Parameter inclusion in the mortality models was driven by a species level of shade-tolerance; this was especially important when evaluating silvicultural effects through the inclusion of competition variables. In this study, the use of a generalized logistic survival model using trees with DBH of 4 cm or greater resulted in accurate estimates for larger trees, but poor results for smaller trees. Often growth and yield models create separate models for small trees (regeneration or ingrowth) as small trees experience higher mortality rates. However, when a separate model is used for small trees, the time at which the small tree model passes data to the 195 large tree model must be determined. This may be determined subjectively. Alternatively, the division between the large and small tree models may be objectively based on the size at which trees pass from the stand initiation stage to the stem exclusion stage of stand development. Even if the optimal time is determined, there may be inconsistencies in the mortality estimates when using separate models. Further research is needed on modeling survival and mortality across all size classes to improve accuracy, and perhaps also using a process-modeling approach. 6.1.2 Growth The Box and Lucas (1959) model was used to model tree-level growth and a parameter prediction approach was used to determine variables affecting growth parameters. While the model provides a great deal of flexibility, care should be taken when fitting the model as it is quite sensitive to the initial starting parameter estimates and requires bounds be put in place to restrict the model in order to provide logical outcomes. Similar to the mortality models, the impacts of competition variables in each model were influenced heavily by the level of shade tolerance and by the species. For untreated stands, competition variables which were not significant, but variables such as basal are per ha may be significant particularly for shade tolerant species, and should be taken into account in the model building process. While the diameter and height growth models were developed independently, the two are not exclusive. The models predicted well for plots which were both fertilized and thinned independently, further research is needed for plots which receive thinning and fertilization in combination. 196 6.1.3 Experimental design vs. managed stands It should be noted that this research was conducted using data for trees on actively managed stands. While prior research has been done to analyze the effects of thinning and fertilization on tree mortality and growth, most studies have been conducted using an experimental design approach. An experimental study provides information regarding how selected variables differ from one specific management practice to another specific management practice under controlled stand-level conditions, but does not provide information across a wide spectrum of management practices and/or stand-level conditions. This study was conducted using data collected across a landscape and models will be applied to the same landscape. Also, variables used in the developed models are commonly available in forest inventory databases, whereas models from experimental design studies often include variables such as crown ratio that are not commonly measured.. 6.2 Management implications Individual tree distance independent models grow an existing forest inventory database tree by tree. This is important for modeling short-term productivity, as well as addressing longer term timber supply. The models developed in this thesis, in conjunction with a regeneration model developed in the affiliated research project, were implemented within FORest Growth Engine (FORGE). FORGE is a hybrid model designed to examine the effects of varying treatment processes on the growth and mortality of trees on Vancouver Island (Smith 2005). This model was designed to: 1) model trees (planted or natural) from bare-ground or already established stands; 2) account for the effects of variable retention on growth; and 3) examine the effects of 197 operational thinning and fertilization. This model will assist in addressing management issues concerning levels of silvicultural investment, harvesting planning, and the need for information on short-term tree and log profiles. Specifically, FORGE was designed for use by Island Timberlands, LP for stand projection and as an aid for management decisions. Thus, the research undertaken here was directly concerned with modeling complex variable retention stands using an extensive existing and newly collected database. These models provide a framework for synthesizing several initiatives and addressing a variety of management issues, including forest growth, potential timber supply and silvicultural investment. 6.3 Future research To meet a variety of societal and conservation needs, more retained trees are being left in harvested areas in the Pacific Northwest. The use of variable retention following harvest is a relatively new forest management practice (Beese et al. 2001, Beese et al. 2005). This practice is believed to have the potential to provide a wide array of forest values. However, since it is a new practice, there is little information on the longer term outcomes following harvest. Also, the variety of retention levels and spatial distributions of retained trees is large, making it difficult to transfer information from one retention area to another area. Models can be used to simulate a variety of variable retention levels and patterns for long forecast periods. Since only early results are currently available for limited ranges of retention and dispersion patterns, models can be useful as a mechanism for generalizing results to other conditions. FORGE is a microclimate based individual tree growth model that is responsive to light levels (Smith et al. 2010). Although the tree growth and mortality models within 198 FORGE are aspatial, this model was developed to examine the impact of leaving more trees, such as under variable retention systems (Mitchell and Beese 2002, Bunnell and Dunsworth 2009), on the growth of adjacent regenerating trees. FORGE models regeneration separately from older trees, where older trees are defined as those older than 20 years to greater than 100 years. Although great care was taken in examining and testing each tree submodel (including random validation tests and careful examination of error structure), the FORGE model requires further testing and validation. To reflect changes in growth and mortality because of changed light conditions under variable retention, light availability is modeled in FORGE based on the height and leaf area of the retained trees. Since variable retention systems are relatively new for Coastal BC and no remeasurement data over approximately 10 years old exists, the light modifier developed using regeneration data was also used for older trees in the FORGE model. This will need to be replaced once improved empirical evidence for older trees becomes available. Once this is replaced, the model can be used to examine the effects of light on older tree growth. In addition, this model was developed for the three main species of interest, Douglas-fir, western redcedar, and western hemlock. While these are the three main species of interest for timber in this area, other species such as red alder (Alnus rubra Bong.) are prolific and exist within the same stands. Models need to be developed for these species and included in the FORGE framework. The models were developed using standard permanent sample plot data. Many of the measurement periods for a plot were highly irregular, with measurement intervals ranging from one to 17 years. While research has been conducted to look at different methods for estimating parameters 199 within survival and diameter growth models developed for irregular measurement periods for single-species plantations (e.g., Cao and Strub 2008), these methods have not been applied to mixed-species stands. While existing taper models were examined for this area for a different project, more research needs to be conducted on the effect of variable retention on tree taper. 200 6.4 References Aubry, K.B., Halpern, C.B., and Peterson, C.E. 2009. Variable-retention harvests in the Pacific Northwest: A review of short-term findings from the DEMO study. For. Ecol. Manage. 258(4): 398–408. Beese, W.J., Dunsworth, B.G., and Smith, N.J. 2005. Variable retention adaptive management experiments: testing new approaches for managing British Columbia‟s coastal forests. Paper In: Eds., C. E. Peterson and D.A. Maguire. Balancing Ecosystem Values: innovative experiments for sustainable forestry, USDA Forest Serv., PNW-GTR-635: 55-4. Beese, W.J., Dunsworth, B.G., and Perry, J. 2001. The Forest Project: three-year review and update. J. Ecoforestry 16(4): 10-17. Box, G.E.P. and Lucas, H.L. 1959. Design of experiments in non-linear situations. Biometrika 46(1/2): 77-90. Bunnell, F.L and Dunsworth, G.B. (eds). 2009. Forest biodiversity: learning how to sustain biodiversity in managed forests. UBC Press, Vancouver, Toronto 349p. Cao, Q. and Strub, M. 2008. Evaluation of four methods to estimate parameters of an annual survival and diameter growth model. For Sci. 54(6): 617-624. Mitchell, S.J. and Beese, W.J. 2002. The retention system: Reconciling variable retention with the principles of silvicultural systems. For. Chron. 78(3):397–403. Peterson, C.E. and Anderson, P.D. 2009. Large-scale interdisciplinary experiments inform current and future forestry management options in the U.S. Pacific Northwest. For. Ecol. Manage. 258(4): 409–414. 201 Smith, N.J. 2005. Variable retention, microclimate, experimental and modeling projects in TFL 39: Modeling mortality in variable retention stands using the FORGE model. Unpublished internal report. Nanaimo, BC, Canada. Smith, N,J, Raynor, K., Grant, N., Black, T.A., and Hum A. In Review.. A forest microclimate model based on cellular automata: modelling photosynthetically active radiation in openings surrounding aggregated retention.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Growth of British Columbia coastal species in response...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Growth of British Columbia coastal species in response to thinning and fertilization treatments Rathbun, Leah C. 2010
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Growth of British Columbia coastal species in response to thinning and fertilization treatments |
Creator |
Rathbun, Leah C. |
Publisher | University of British Columbia |
Date Issued | 2010 |
Description | The successional processes of the mixed-species Pacific coastal temperate rain forests of British Columbia (BC), Canada, are defined by gap dynamics, where small-scale disturbances, mainly due to windthrow, create openings in the canopy necessary for regeneration. Douglas-fir (Pseudotsuga menziesii var. menziesii (Mirb.) Franco) is the dominant, pioneer species in this area and western hemlock (Tsuga heterophylla (Raf.) Sarg.) and western redcedar (Thuja plicata Donn) are the late-successional, shade-tolerant species. Silvicultural systems such as variable retention systems have been applied to many of the secondary growth mixed-conifer forests. Variable retention in this area is designed to differ dramatically from stand to stand. This approach differs from the traditional even-aged management applied to forests of the Pacific Northwest coast. In this study, a model-based approach was used to investigate how multiple treatment interventions as a part of active management across a landscape affect mortality and growth within actively managed stands. There is a need for this information as current growth and yield models used in this area are limited by either the number of species which can co-exist in a stand (e.g., the model TASS of BC) or are limited by the need for data not commonly obtained in inventory databases (e.g., the models FVS and ORGANON of USA). Additionally, no growth and yield models have been developed to include variable retention systems, where a variety of thinning intensities and spatial patterns, timing of thinning and fertilization treatments, and number of treatments are used. Mortality, diameter increment, and height increment models were developed and the effects of thinning, fertilization and the combination of thinning and fertilization were examined for Douglas-fir, western redcedar, and western hemlock. For each species, shade-tolerance was found to impact the possible predictor variables included in model development. The use of a generalized logistic survival model resulted in accurate estimates for larger trees, but poor results for smaller trees. To model the effects of fertilization, additional fertilization effect variables were included in the models; conversely, thinning effects were modeled using the immediate change in state variables such as basal area of larger trees which occurred immediately following thinning. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-11-30 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0071478 |
URI | http://hdl.handle.net/2429/30235 |
Degree |
Doctor of Philosophy - PhD |
Program |
Forestry |
Affiliation |
Forestry, Faculty of |
Degree Grantor | University of British Columbia |
GraduationDate | 2011-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2011_spring_rathbun_leah.pdf [ 4.86MB ]
- Metadata
- JSON: 24-1.0071478.json
- JSON-LD: 24-1.0071478-ld.json
- RDF/XML (Pretty): 24-1.0071478-rdf.xml
- RDF/JSON: 24-1.0071478-rdf.json
- Turtle: 24-1.0071478-turtle.txt
- N-Triples: 24-1.0071478-rdf-ntriples.txt
- Original Record: 24-1.0071478-source.json
- Full Text
- 24-1.0071478-fulltext.txt
- Citation
- 24-1.0071478.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0071478/manifest