AN ANALYTICAL PLATFORM FOR CUMULATIVE IMPACT ASSESSMENT IN NORTHEASTERN BRITISH COLUMBIA by Bogdan Mihai Strimbu B.Sc. Transilvania University, 1992 M.Sc., The University of British Columbia, 2003 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) September 2009 ©Bogdan Mihai Strimbu, 2009 ii ABSTRACT The combined influence on the environment of all projects occurring in a single area is evaluated through cumulative impact assessments (CIA), which consider the consequences of multiple projects, each possibly insignificant on its own, yet important when evaluated collectively. Traditionally, the future human activities are included in CIA using an analytical platform, commonly based on complex models that supply precise predictions but with asymptotically null accuracy. To compensate for the lack of accuracy of the current CIA I have proposed a shift in the paradigm governing the CIA. The paradigm shift involves a change in the focus of CIA investigations from the detailed analysis of one unlikely future to the identification of the patterns describing the future changes in the environment. To illustrate the approach, a set of 144 possible and equally likely futures were developed that aimed to identify the potential impacts of forest harvesting and petroleum drilling on the habitat suitability of moose and American marten. The evolution of two measures of habitat suitability (average HSI and surface of the stands with HSI>0.5) was investigated using univariate and multivariate repeated measures. Both analytical techniques (i.e., univariate or multivariate) revealed that the human activities could induce at least one cycle (with a period larger than 100 years), in the moose and American marten habitat dynamics. The planning period was separated into three or four distinct periods (depending on the investigation methodology) following a sinusoidal pattern (i.e., increase – constant – decrease in the habitat suitability measures). The attributes that could induce significant changes in the environment are the choice of harvesting age and the valued ecosystem component. The choice of the valued ecosystem component is critical to the analysis and could change the conclusions of the CIA. iii TABLE OF CONTENTS ABSTRACT ................................................................................................................................. ii TABLE OF CONTENTS ............................................................................................................. iii LIST OF TABLES ....................................................................................................................... vi LIST OF FIGURES ................................................................................................................... vii ACKNOWLEDGEMENTS ........................................................................................................ viii DEDICATION ............................................................................................................................. ix CO-AUTHORSHIP STATEMENT ............................................................................................... x 1. INTRODUCTION ................................................................................................................. 1 1.1 Objectives ..................................................................................................................... 3 1.2 Hypotheses and assumptions ...................................................................................... 4 1.3 Thesis structure ............................................................................................................ 5 1.4 References ................................................................................................................... 8 2. PARADOXICAL EMPIRICAL INVESTIGATIONS OF COMPLEX ENVIRONMENTAL SYSTEMS ................................................................................................................................. 11 2.1 References ................................................................................................................. 19 3. A DETERMINISTIC HARVEST SCHEDULER USING PERFECT BIN PACKING THEOREM ................................................................................................................................ 21 3.1 Introduction ................................................................................................................. 21 3.2 Methods ...................................................................................................................... 22 3.2.1 Data description .................................................................................................. 22 3.2.2 Optimal harvesting age and mean annual increment .......................................... 24 3.2.3 Maximum volume harvested annually determined using MAI ............................. 26 3.2.4 Adjusting MVHA: First Fit Decreasing algorithm ................................................. 28 3.2.5 Simulated annealing ............................................................................................ 32 3.3 Results ........................................................................................................................ 35 3.4 Discussion .................................................................................................................. 39 3.5 Conclusion .................................................................................................................. 41 iv 3.6 References ................................................................................................................. 43 4. AN ANALYTICAL PLATFORM FOR CUMULATIVE IMPACT ASSESSMENT IN NORTHEASTERN BRITISH COLUMBIA ................................................................................. 47 4.1 Introduction ................................................................................................................. 47 4.2 Methods ...................................................................................................................... 49 4.2.1 Paradigm shift ..................................................................................................... 49 4.2.2 Analytical platform ............................................................................................... 49 4.2.3 Generating futures ............................................................................................... 52 4.2.3.1 Petroleum drilling ......................................................................................... 53 4.2.3.2 Forest harvesting ......................................................................................... 58 4.2.3.3 Moose and American marten habitat ........................................................... 60 4.2.4 Cumulative impact assessment ........................................................................... 62 4.3 Results ........................................................................................................................ 64 4.3.1 The development of the set of futures ................................................................. 64 4.3.2 Habitat suitability index for the set of futures ....................................................... 71 4.4 Discussion .................................................................................................................. 76 4.5 Conclusion .................................................................................................................. 80 4.6 References ................................................................................................................. 82 5. MULTIMODELING ASSESSMENT OF THE IMPACT OF PETROLEUM DRILLING AND FOREST HARVESTING ON WILDLIFE HABITAT ................................................................... 92 5.1 Introduction ................................................................................................................. 92 5.2 Methods ...................................................................................................................... 93 5.2.1 Development of independent and equally likely futures ...................................... 97 5.2.1.1 Generation of futures for oil and gas drilling ................................................ 97 5.2.1.2 Generation of futures for forest harvesting ................................................. 101 5.2.2 Assessment of the cumulative impact of oil drilling and forest harvesting ........ 104 5.2.2.1 Quantification of moose and American marten habitat .............................. 104 5.2.2.2 Analysis of futures induced by oil and gas drilling and forest harvesting ... 105 v 5.3 Results ...................................................................................................................... 108 5.4 Discussion ................................................................................................................ 122 5.5 Conclusions .............................................................................................................. 125 5.6 References ............................................................................................................... 127 6. CONCLUSIONS .............................................................................................................. 133 6.1 Status of the hypotheses and significance of the research ...................................... 137 6.2 Potential applications and future research ............................................................... 138 6.3 References ............................................................................................................... 140 APPENDICES ......................................................................................................................... 142 Appendix A .......................................................................................................................... 142 References ...................................................................................................................... 150 Appendix B .......................................................................................................................... 151 Appendix C .......................................................................................................................... 153 vi LIST OF TABLES Table 3.1. Timber harvesting land base (THLB) ....................................................................... 24 Table 3.2. Comparison between the MVHA and the AACs supplied by LP and the adjusted FFD and SA, for the adjacency delay of 1 year and minimum harvest age 120 years ............ 39 Table 4.1 Models describing the evolution of the number of petroleum wells in the northeastern British Columbia ........................................................................................................................ 65 Table 4.2 Variograms of the study areas .................................................................................. 68 Table 4.3. Timber harvesting land base and AAC for the four scenarios ................................. 70 Table 5.1 The harvest scheduling constraints ........................................................................ 102 Table 5.2 Well drilling .............................................................................................................. 110 Table 5.3 The AAC obtained using the two algorithms (i.e., SA and FFD) and the two harvesting ages (i.e., preset minimum age and Age max(MAI)) .................................................. 113 Table 5.4 Principal components analysis ................................................................................ 119 Table 5.5 Canonical discriminant analysis .............................................................................. 121 Table A.1 Time-series model for the monthly average of number of active rigs.. ................... 144 Table A.2 The models for the wells number in northeastern British Columbia. . .................... 149 vii LIST OF FIGURES Figure 2.1 Coefficient of the autoregressive term of the number of active rigs and the dynamics of the moving average for the error model. T. ........................................................................... 15 Figure 2.2 Confidence interval (α=0.02) for coefficient of the autoregressive term (i.e., number of wells drilled per time-step), and the number of economic variables present in the model describing the total number of wells drilled per time-step.. ....................................................... 17 Figure 3.1 Study areas .............................................................................................................. 23 Figure 3.2 First fit decreasing scheduler algorithm that fulfils the PBPT corollary requirements and queues the stands according to the attribute considered for optimization ......................... 31 Figure 3.3 Current volume and area distribution with respect to the harvest timing.. ............... 36 Figure 3.4. AAC dependency on the algorithm and the green-up/ adjacency delay for Doig area (3.4a), Moberly area (3.4b) and Fort Nelson area (3.4c). ................................................. 37 Figure 4.1 Study areas .............................................................................................................. 50 Figure 4.2 Analytical framework for the shift of the CIA paradigm (a) and the development of set of equally likely and independent futures (b). ...................................................................... 52 Figure 4.3 Evolution of the number of wells drilled annually ..................................................... 67 Figure 4.4 Quadrants within the quarter section grid used to determine the well forecast grid 69 Figure 4.5 The evolution of the AAC and of the land disturbance during the planning period .. 70 Figure 4.6 The evolution of the average HSI during the planning period .................................. 72 Figure 4.7 The evolution of the areas with HSI>0.5 using SA as forest harvesting algorithm .. 75 Figure 5.1 Study areas .............................................................................................................. 94 Figure 5.2 Development of the set of futures ............................................................................ 96 Figure 5.3 Cluster analysis revealing the significant models. ................................................. 118 Figure 5.4 PCA separating the significant periods (a) and significant combination of study area harvesting algorithm (b). ......................................................................................................... 120 Figure 5.5 Canonical discriminant analysis indicating the impact of harvesting algorithm on delineating the significant years during the planning period ................................................... 122 viii ACKNOWLEDGEMENTS This thesis is the culmination of the journey that started a quarter of a century ago with the debates about the purpose of life. In my case, this journey is best described by two people: Aristotle ‘’All men by nature desire knowledge” (Metaphysics) and Sir Isaac Newton ‘’If I have seen further it is by standing on the shoulders of giants” (Letter to Robert Hooke, 1675). The completion of this thesis would have not been possible without the support of the people embarked totally or partially in the quest for knowledge. I am particularly grateful to the following people who have contributed to the conclusion of this research endeavor: My supervisor, Dr. John Innes, who guided and trusted me during the last decade. His capacity of seeing beyond figures amazes me now as it did 10 years ago. This thesis contains his delicate but forcefully touches. To my research committee, consisting of Dr. Vlada Limic and Dr. John Nelson, who supported my environmental modeling approach from the beginning. Their input was invaluable and the long conversations from the Department of Mathematics or from the Faculty of Forestry can be encountered almost everywhere in my dissertation. To Dr. Valerie Lemay, Dr. Peter Marshall, Dr. Younes Alila, Dr. Paul Wood, Dr. Richard Anstee and Dr. Joel Friedman who helped, challenged and focused me. To the people who directly participated in this research, particularly Constantin Visan from the Oil and Gas Commission of British Columbia, Noel Gardner from the Treaty 8 Tribal Association, and Warren Jukes and Don Rosen from Canfor. I would like to thank my Sustainable Forest Management lab colleagues, particularly, Denise Allen, Gordon Hickey, and Craig Nitschke for being there when I needed them. I am grateful to all the people that I have encountered, which one way or another shaped my life, the Ţoiu’s Galleria. On a personal note, I want to thank my wife for her never-ending understanding and support.. This research would have not been possible without the financial support of the Sustainable Forest Management Network. ix DEDICATION To my family x CO-AUTHORSHIP STATEMENT The thesis contained four chapters that integrate the collaboration between Prof. Innes, my research supervisor, and me. The identification and design of the research program, performing of the research, analysis of the data and preparation of all four manuscripts were conceived and executed by Bogdan Strimbu, with feedback from John Innes. The results from Chapter 3 were obtained using a computer program developed by Victor Strimbu. 1 1. INTRODUCTION The impact of human activities on the environment ranges from relatively innocuous industries such as tourism, to more invasive extractive industries such as forestry, mining or oil and gas exploration. The harmful consequences of industrial development has become notable as environmental effects associated with these activities challenge the capacity of an ecosystem to incorporate naturally occurring disturbances (Fuhrer 2000; Marcu 1981). Over time, the accumulation of stresses resulting from human activities could transform an ecosystem so profoundly that recovery to a preexisting desirable state is no longer possible using existing conventional processes (Gore et al. 1990; Toffler 1970). To prevent the environment changing beyond socially acceptable limits, the potential effects of human-induced perturbations are generally evaluated by environmental impact assessment studies (Koornneef et al. 2008). The main difficulty in assessing the environmental response to human developments lies in the way that most economic activities have limited impacts when assessed singly, but generally have additive or synergistic effects when considered in the context of all past, present and foreseeable future activities. The combined influence of a suite of projects occurring in a predefined area is therefore evaluated using cumulative impact assessments studies (CIA). These consider the consequences of multiple projects, each possibly insignificant on its own, yet important when considered collectively (Council on Environmental Quality 1969). The CIA literature addresses the qualitative aspects of the environmental changes induced by the evolution of society (Contini and Servida 1992; Glasson et al. 1994; Marr 1997; Masera and Colombo 1992; Reid 2001). However, fewer studies have elaborated specific methods that can be used in evaluating the cumulative impacts associated with human developments (Spaling and Smit 1994). Most current CIAs are performed using either complex modular models (Dube et al. 2006; Voinov et al. 2004) or different heuristic techniques (Stakhiv 1988). The two CIA methods have been possible as an extensive computational effort has been undertaken for the last 30 years to represent the response of a specific set of environmental attributes to different human-induced transformations of the environment [e.g., DHSVM in hydrology, FORECAST in ecology or FOREPLAN in forest harvesting]. The CIA literature abounds in detailed investigations founded on methods initially developed to assess particular activities, such as land-use change associated with urban development (Dickert and Tuttle 1985; Stakhiv 1988), water resource developments (Dee et al. 1973; Gosselink and Lee 1989) and construction (Leopold et al. 1971). These have been enhanced so that they can be used as analytical 2 support for CIA (Rotmans and van Asselt 2001). Regardless of the method adopted, the current CIA modeling platforms are based on the assumption that small incremental changes are the driver of environmental dynamics (Cherp et al. 2007; Dube et al. 2006; Hegmann et al. 1999), and that these changes can be represented mathematically. The complexity of the CIA investigation and the intricate structure associated with modular modeling has led to CIA based on information supplied by one analytical platform (Voinov et al. 2004): one CIA - one complex model. The mono-analytical approach is preferred for CIA investigations for two reasons. First, the unique platform supplies an accurate assessment of the impact of current activities on the environment, and secondly, the modular structure of the investigation platform allows the adjustment of different modules to address the complexity of the CIA. However, the two strengths of the approach diminish when the temporal dimension is incorporated into CIA, as the accuracy of the assessment decreases exponentially with the number of entities (e.g., independent modules, equations or grid cells). The reduction in accuracy may be sufficiently acute for the CIA to convey erroneous conclusions; as the probability of occurrence of a complex model converges to zero. The asymptotic null probability of the results supplied by the complex mono-analytical methodologies employed by CIA is rooted in the frameworks of the models frameworks, which can be fully empirical (Patil et al. 2002), fully process-based (Leimbach and Jaeger 2005) or a combination of the two (Wu and David 2002), in the sense of Korzukin et al. (1996). CIA using models based exclusively on processes and theories have difficulty representing environmental dynamics when socio-economic, biophysical and land- management variables are used (Verburg et al. 2004). Additionally, process–based CIA does not supply a measure of the accuracy of the results, and assume that the environmental dynamics are fully known and predictable. Therefore, the current complex mono-analytical CIA models incorporate empirical components that lead to precise predictions but with reduced accuracy, as the probability of occurrence of the outcome of a model is: Probability of occurrence = 1/ (the number of equally likely results) ≤ ≤ (the smallest number of equally likely results supplied by each empirical equation)-n ≤ (1.1) ≤ n n − ∞→ + )11(lim = 0 where n is the number of non-overlapping entities used by the model and the two 1s correspond to the smallest (respectively largest) value of the confidence interval supplied by the empirical equations. 3 The null probability of occurrence of the results supplied by the multi-modular approach is enforced when possible future activities are considered in the environmental assessment. The accuracy of the CIA decreases even more when the location of the activities, in addition to the timing, is of interest. The description of the spatial dynamics requires the addition of an extra module or set of equations to describe the spatial dynamics which would increase the n value in Equation 1.1, making the convergence to 0 faster. Therefore, considerable effort is required in reducing the uncertainties determining the future evolution of the environment, as CIA of the present state is focused mainly on fulfilling the environmental constraints. The efforts to increase the accuracy of the predicted environmental response to human induced perturbations consider the CIA an instrument to analyze and assess the past, present and future human induced disturbances (Spaling and Smit 1993). Viewed as a scientific tool, CIA is primarily a generator of information to be used in the planning process. The usage of the CIA as a provider of inputs to the decision making process is the natural result of the multi-modular approach. Within this context, the legislative and administrative framework has a significant influence in the implementation of the results supplied by the CIA investigation. To ensure that economic activities are not transforming the environment to states unacceptable by the society, in North America, Canada and the USA promoted specific legislation to implement the CIA; Canada in 1992 (Government of Canada, 1992) and the USA in 1969 (Council of Environmental Quality, 1969). 1.1 Objectives To address the lack of confidence in the techniques currently used for CIA, the present research proposes a change in the focus of the environmental assessment when the prediction of future activities is of interest. In the case when potential developments could impact the environment, the CIA should focus on identifying patterns in the environmental attributes that could lead to undesirable environmental outcomes. To identify the patterns associated with unacceptable states of the environment, the research presented here promotes the development of a set of future environments, each of them possible and equally likely. The set of futures is similar to the states of the molecules in a fluid, with each future resembling the trajectory of a molecule (Boltzmann 1995). As in statistical thermodynamics, the set of futures would be used to represent globally the state (i.e., the possible evolution of the environment during the length of the future) but, because the futures are objects with the probability of existence asymptotically 0 [i.e., (probability of occurrence of a year during the future) -future duration ≈ 0] compared with the 4 particles from a fluid that are physical entities (i.e., probability of existence is 1), the CIA investigation will focus on describing the similarities or differences existing within the state, rather than quantifying a specific attribute of the state. Therefore, the objectives of the present research are to identify the patterns describing significantly different environmental states and the environmental attributes associated with the respective patterns. 1.2 Hypotheses and assumptions A set of hypotheses and assumptions were developed to guide this research. The assumptions used to delineate the possible patterns associated with an environmental state were: 1. The future is unknown; 2. The theoretical framework used to quantify the attributes representing the future does not change; and 3. The influence of human activities on the environment can be separated from the effects of natural events. The three assumptions served as the basis for the testing of the two hypotheses proposed by the CIA paradigm and further served as the foundation for the development of the set of futures. The two hypotheses tested in the present study were: 1. There is an identifiable pattern associated with each state of the environment; and 2. The chosen analytical tools play a significant role in the description and forecast of future environments. To test the hypotheses and to illustrate the proposed CIA approach the research presented in this thesis considered the impact of two human activities (forest harvesting and petroleum drilling) on two valued ecosystem components (VEC) from northeastern British Columbia [moose (Alces alces) and American marten (Martens Americana], as defined by Duinker (1994). Based on the recommendation of the Treaty 8 Tribal Association, three study areas (i.e., Moberly, Doig and Fort Nelson), totaling 1.7 million ha, were used to represent the northeastern part of British Columbia (i.e., the 17 million ha east of the Kechika – Finlay corridor). The research used data supplied by the Ministry of Sustainable Resource Management (2004) for forestry and wildlife, and by the Oil and Gas Commission (2005), Baker Hughes Inc. (2005) and Energy Information Administration (2005) for petroleum wells. 5 1.3 Thesis structure The objectives of the research were achieved by integrating the results supplied by a newly developed set of models describing petroleum drilling and forest harvesting. The set of models and their integration into one analytical platform, which enabled the identification of the patterns associated with different environmental states, are presented individually as chapters of the thesis. The following enumeration presents the title and outlines the contents of each chapter in addressing the research objectives. Chapter 2. Paradoxical empirical investigations of complex environmental systems Chapter 2 serves two purposes: first, the development of the temporal evolution of the petroleum drilling, and secondly, the testing of the second hypothesis of the study (i.e., the analytical tools influence the description of the environment). Chapter 2, by extending the findings of Seppelt and Richter (2005), Descartes (1850) and Bacon (1855), advocate the usage of a stochastic approach (i.e., set of possible and equally likely futures) as an alternative to the current mono-modeling platform approach. To test the second hypothesis of the study, Chapter 2 challenges the present methodologies of representing the dynamics of complex environmental systems, which use mathematical models whose validity relies on the assumption that knowledge regarding the physical world is drawn from observations or experiments subject to the principles of reasoning. In the chapter it is argued that this assumption should be modified, as empirical relationships describing environmental systems have a haphazard character and fail to converge to an assumed model. Using an approach resembling the Liar Paradox it was proved that the scientific inquiry could be influenced by the analysis itself, and does not depend only on the data. Consequently, the chapter findings indicate either the inability to reveal the true relationships within complex environmental systems or that an analytical factor should be included in the description of those systems. Chapter 3. A deterministic harvest scheduler using perfect bin - packing theorem Chapter 3 forecasts the evolution of the landscape as a result of the forest harvesting occurring in the three study areas. Currently, in northeastern British Columbia forest planning uses heuristic techniques, techniques promoted by the need to solve complex problems that cannot be solved using mixed integer programming. In Chapter 3 it is proved that for merchantability standards ensuring the perfect bin-packing theorem (PBPT), the maximum volume that can be harvested annually equals the sum of the maximum MAI of the stands. The method to compute the maximum volume accommodates optimality criteria at the stand level, regarded as 6 maximum MAI, and at the forest level, regarded as maximum annual allowable cut (AAC). The harvestings were scheduled by adjusting the first fit decreasing algorithm (FFD) to the spatial constraints associated with forest planning. When PBPT conditions were not met, a mixed integer programming solution was developed to adjust the merchantability standards of the stands to the distributional requirements of the PBPT, an adjustment that ensured the optimal performance of the FFD. The adjusted FFD was compared with linear programming and simulated annealing (SA) using two harvesting ages (i.e., one based on MAI maximization and one determined as the minimal age) and the same set of spatial-temporal constraints for three areas. Chapter 4. An analytical platform for cumulative impact assessment in northeastern British Columbia Petroleum drilling and forest harvesting are the main activities changing the landscape in northeastern British Columbia. The two activities are used to represent the combined influence on the environment of all projects occurring in the area; influence evaluated by cumulative impact assessments (CIA) studies. Chapter 4, which starts by integrating the results from Chapters 2 and 3, adds a spatial dimension to the petroleum drilling. To compensate for the lack of accuracy of the current mono-analytical CIA, Chapter 4 proposes a shift in the paradigm governing the CIA. The paradigm shift involves a change in the focus of CIA investigations from the detailed analysis of one unlikely future to the identification of the patterns describing the future changes in the environment. To illustrate the approach, a set of 144 possible and equally likely futures were developed that aimed to identify the potential impacts of forest harvesting and petroleum drilling on the habitat suitability of moose and American marten. The habitat of the two species was assessed using a habitat suitability index (HSI). Two statistics were used to quantify the impact of the two human activities on the habitat of moose and American marten (average HSI and Area HSI>0.5). Based on the two HSI statistics, the distinct environmental states, which describe the patterns associated with the future changes of the environment, were identified using univariate repeated-measures analysis. The univariate analysis supplied an answer to the first research hypothesis, namely the pattern associated with each state of the environment can be identified. Chapter 5. Multimodelling assessment of the impact of the petroleum drilling and forest harvesting on the wildlife habitat in northeastern British Columbia Chapter 5 builds on the results of Chapter 4 and expands the findings of Chapter 4 by providing a multivariate perspective of the environmental dynamics. The analysis in this chapter uses the 7 same approach of describing the environment [i.e., a set of probably and equally likely evolutions of two HSI statistics (average HSI and Area HSI>0.5)]. Chapter 5 enhances the univariate investigation of Chapter 4 by clearly delineating the significant environmental states, and identifying the attributes associated with the respective states. The multivariate perspective on the environmental dynamic complemented and confirmed the main results of the univariate approach as multivariate repeated measures analysis, hierarchical cluster analysis, principal component analysis and canonical discriminant analysis, the main tools used to identify the impact of the forest harvesting and petroleum drilling on the habitat of moose and American marten, supplied similar results. The agreement between the multivariate and univariate investigations supported the first hypothesis of the study (i.e., the human activities induce an identifiable pattern within the environment). Chapter 6. Conclusions Chapter 6 concludes the thesis by relating the previous chapters and summarizing the important results and findings of the research. The chapter also discusses the strengths, weaknesses and overall significance of the research, presents the status of the working hypotheses, and evaluates the current knowledge related to the CIA. The thesis discusses the potential applications of the research findings, specifically the extension of the proposed CIA paradigm change to other environmental investigations that routinely use complex models (such as climate dynamics and land-use planning), and the impact of the new bounds of the objectives of the forest planning on the existing heuristic algorithms. The thesis ends by identifying two directions for future research: one practical (inclusion of other human activities within the proposed CIA framework), and one theoretical (assessment of the sensitivity of the CIA to the violation of the analytical assumptions and the impact of different distributions describing the environment on the analytical platform). 8 1.4 References 1. Bacon, F. 1855. The New Organon. Oxford University Press, Oxford. 338 p. 2. Baker Hughes Inc. 2005. Worldwide Rig Counts - Current & Historical Data. Flaharty, G. and Shiels, G. (eds.). Houston TX, Baker Hughes Inc. 3. Boltzmann, L. 1995. Lectures on gas theory. Dover Publications. 490 p. 4. Cherp, A., Watt, A. and Vinichenko, V. 2007. SEA and strategy formation theories: From three Ps to five Ps. Environmental impact assessment review 27(7):624-644. 5. Contini, S. and Servida, A.. 1992. Risk analysis and environmental impact studies. P. 79- 104 in Environmental impact assessment, Colombo, A.G. (ed.), Kluwer Academic Publisher, Dordrecht. 6. Council on Environmental Quality. The National Environmental Policy Act. 42. 1969. 7. Dee, N., Baker, J., Drobny, N.,Duke, K., Whitman, I. and Fahringe, D. 1973. Environmental Evaluation System for Water Resource Planning. Water Resources Research 9(3):523-535. 8. Descartes, R. 1850. Discourse on method of rightly conducting the reason, and seeking truth in the science. Sutherland and Knox, Edinburgh. 118 p. 9. Dickert, T.G. and Tuttle, A.E. 1985. Cumulative impact assessment in environmental planning; a coastal wetland watershed example. Environmental impact assessment review(5):37-64. 10. Dube, M., Johnson, B. Dunn, G., Culp, J., Cash, K., Munkittrick, K., Wong, I., Hedley, K., Booty, W., Lam, D., Resler, O., and Storey, A. 2006. Development of a New Approach to Cumulative Effects Assessment: A Northern River Ecosystem Example. Environmental Monitoring and Assessment 113(1 - 3):87-115. 11. Duinker, P.N. 1994. Cumulative effects assessment: what's the big deal? P. 11-24 in Cumulative effects assessment in Canada: from concept to practice, Kennedy, A.J. (ed.), Alberta Association of Professional Biologists, Calgary. 12. Energy Information Administration. Spot Prices for Crude Oil and Petroleum Products. 2005. U.S. Department of Energy. 13. Fuhrer, E. 2000. Forest functions, ecosystem stability and management. Forest Ecology and Management 132(1):29-38. 14. Glasson, J., Therivel, R., and Chadwick, A. 1994. Introduction to environmental impact assessment. UCL Press, London. 342 p. 9 15. Gore, J.A., Kelly, J.R. and Yount, J.D. 1990. Application of Ecological Theory to Determining Recovery Potential of Disturbed Lotic Ecosystems - Research Needs and Priorities. Environmental management 14(5):755-762. 15. Government of Canada. 1992. Canadian Environmental Impact Assessment Act. 37. 16. Gosselink, J.G. and Lee, L.C. 1989. Special Issue - Cumulative Impact Assessment in Bottomland Hardwood Forests. Wetlands 9: 93-174. 17. Hegmann, G., Cocklin, C., Creasey, R., Dupuid, S., Kennedy, A., Kingsley, L., Ross, W., Spaling, H., and Stalker, D. 1999. Cumulative Effects Assessment Practitioners Guide. Hull, QC, Canadian Environmental Assessment Agency. 143 p. 18. Koornneef, J., Faaij, A. and Turkenburg, W. 2008. The screening and scoping of Environmental Impact Assessment and Strategic Environmental Assessment of Carbon Capture and Storage in the Netherlands. Environmental impact assessment review 28(6):392-414. 19. Korzukhin, M.D., TerMikaelian, M.T., and Wagner, R.G. 1996. Process versus empirical models: Which approach for forest ecosystem management? Canadian Journal of Forest Research-Revue Canadienne de Recherche Forestiere 26(5):879-887. 20. Leimbach, M. and Jaeger, C. 2005. A modular approach to Integrated Assessment modeling. Environmental Modeling and Assessment 9(4):207-220. 21. Leopold, L.B., Clarke, F.E., Hanshaw, B.B., and Balsley, J.R.. 1971. A procedure for evaluating environmental impact. Washington D.C, U.S. Geological Survey. 15 p. 22. Marcu, M. 1981. Meteorologie cu elemente de climatologie forestiera. Ceres, Bucharest. 291 p. 23. Marr, K. 1997. Environmental impact assessment in the United Kingdom and Germany. Ashgate Publishing, Aldershot. 327 p. 24. Masera, M. and Colombo, A.G.. 1992. Contents and phases of an EIA study. P. 53-78 in Environmental impact assessment, Colombo, A.G. (ed.), Kluver Academic Publishers, Dordrecht. 25. Ministry of Sustainable Resource Management. TRIM/TRIM II positional files - basemap. [2004], http://srmwww.gov.bc.ca/bmgs/catalog/bcgs_indexes.htm. 2004. Victoria BC, Ministry of Sustainable Resource Management. 26. Oil and Gas Commission. Well Surface Location. 2005. Oil and Gas Commission. 27. Patil, A.A., Annachhatre, A.P. and Tripathi, N.K.. 2002. Comparison of conventional and geo-spatial EIA: A shrimp farming case study. Environmental impact assessment review 22(4):361-375. 10 28. Reid, L.M. 2001. Cumulative watershed effects: then and now. Watershed management Council Network 10(1):24-33. 29. Rotmans, J. and van Asselt, M.B.A.. 2001. Uncertainty management in integrated assessment modeling: Towards a pluralistic approach. Environmental Monitoring and Assessment 69(2):101-130. 30. Seppelt, R. and Richter, O. 2005. "It was an artefact not the result": A note on systems dynamic model development tools. Environmental Modelling & Software 20(12):1543-1548. 31. Spaling, H. and Smit, B. 1993. Cumulative Environmental-Change - Conceptual Frameworks, Evaluation Approaches, and Institutional Perspectives. Environmental management 17(5):587-600. 32. Spaling, H. and Smit, B. 1994. Classification and evaluation methods for cumulative effects assessment. P. 47-65 in Cumulative effects assessment in Canada: from concept to practice, Kennedy, A.J. (ed.), Alberta Association of Professional Biologists, Edmonton. 33. Stakhiv, E.Z. 1988. An Evaluation Paradigm for Cumulative Impact Analysis. Environmental management 12(5):725-748. 34. Toffler, A. 1970. Future Shock. Random House, New York. 505 p. 35. Verburg, P.H., Schot, P.P., Dijst, M.J. and Veldkamp, A. 2004. Land use change modelling: current practice and research priorities. GeoJournal 61(4):309-324. 36. Voinov, A., Fitz, C., Boumans, R. and Costanza, R. 2004. Modular ecosystem modeling. Environmental modelling and software 19:285-304. 37. Wu, J.G. and David, J.L. 2002. A spatially explicit hierarchical approach to modeling complex ecological systems: theory and applications. Ecological Modelling 153(1- 2):7-26. 11 2. PARADOXICAL EMPIRICAL INVESTIGATIONS OF COMPLEX ENVIRONMENTAL SYSTEMS1 The need for accurate environmental predictions coupled with current computational flexibility has led to the development of complex multi-modular models in subjects ranging from hydrology (Wigmosta et al. 1994) to forestry (Kimmins et al. 1999) and cumulative impacts (Voinov et al. 2004). The predictions are based on the assumption that models can quantify the relationship between different environmental attributes. Bas van Fraassen (1980) has argued that models connecting measured data with theoretical hypotheses are empirically adequate but fail to reveal the truth as constant adjustments are required. Therefore, environmental models have more of a heuristic value (Oreskes et al. 1994) as their verification is precluded by a logical fallacy (i.e., affirming the consequent) rooted in the incomplete access to natural phenomena (Bacon 1855). Models based only on processes, theories and laws solve the fallacy problem but face difficulties in numerically representing the relationships between environmental dynamics and socio-economic, biophysical or land-management variables (Verburg et al. 2004). Additionally, process-based environmental models encounter the challenge that even simple descriptions of reality can lead to chaotic behaviors (May 1974). Stochastic models enhance process modeling by including a random element in the equations describing the environmental relationships (Winston 1994), but do not solve the two difficulties mentioned above. The basis of stochastic modeling is the assumption that the equations describing the environmental dynamics are converging, in the sense of Grimmet and Stirzaker (2002), an assumption met only when the conditions of the ergodic theorem are met (Grimmett and Stirzaker 2002). Therefore, stochastic modeling adds a set of methodological constraints without solving the two process-based challenges (i.e. quantitative representation of the processes and their chaotic behavior). Alternatively, empirical models are representative of the data but do not reveal the structure or the rules governing the relationships between attributes describing the environment (Korzukhin et al. 1996). Therefore, complex environmental relationships are represented by models that are either fully empirical or a combination of empirically-derived equations with process-based models (Korzukhin et al. 1996). Irrespective of the approach, two challenges face the modeling framework. First, the predictions are precise but have reduced accuracy, as the probability of occurrence of the outcome of a model converges to 0: Probability of occurrence = 1/ (the number of equally likely results) ≤ 1 A version of this chapter has been submitted for publication. Strimbu, B.M. and Innes, J.L. Paradoxical empirical investigations of complex environmental systems. 12 ≤ (the smallest number of equally likely results supplied by each empirical equation (2.1) present in a non-overlapping modeling unit)-n ≤ n n − ∞→ + )11(lim = 0 where n is the number of non-overlapping units used by the model and the two 1s correspond to the smallest (respectively largest) value of the confidence interval supplied by the empirical equations. Secondly, the spatial or temporal scale has a significant impact on the outcomes (Refsgaard 1997; Vazquez et al. 2002), as the calibration of a model is scale dependent and usually impacts the results (Chaubey et al. 2005). However, not all environmental models are scale dependent as some allometric models (Enquist and Niklas 2001) exhibit scale- independent relationships, with unchanged relationships for 21 out of the 27 orders of magnitudes in mass that cover the life processes (West and Brown 2004). Besides supplying results with probability of occurrence asymptotically null, the models describing complex environmental systems implicitly rely on an assumption similar to universal biological scaling laws (West et al. 1997; West and Brown 2004): the models are structurally invariant with scale (i.e., the significance and the type of relationships among the components of a model do not change across scales). I argue that the structure of the models describing the relationships between the components of environmental complex systems depends on analytical details that are not necessarily related to the recorded data. Therefore, either the data drawn from observations and experiments are not sufficient to represent the environment or the current scientific tools are unable to reveal the true relationships between the attributes representing the environment. To support this argument I use an approach resembling the Liar Paradox (Mills 1998), namely the hypothesis that the equations describing the relationships among the attributes quantifying the environment are simultaneously linear and nonlinear. I tested the linearity hypothesis by performing a time series investigation, as for a true linear model ∑ ×= i TiiT XaY | the size of the time-step should not influence the linearity of the equations, as long as the models developed for each time-step are correct, where YT represents the dependent variable calculated using the time step T, Xi|T is the ith independent variable for time step T and ai|T is the coefficient of variable Xi|T. The lack of impact of the time step on a linear model is valid if E(YT)=E(YT’)/n where E(YT ) is the expected value of variable YT,, T and T’ are two times steps such that T’ = nT, and n is a natural number. For a true model, the equality of the two expectations holds irrespective the time steps, as 13 )()()()( )()()()( '|| || T i nT i nTinTi i i TiTiTT YEYEXEXE nXEXnEnYEYnE ==== ==== ∑ ∑ ∑ ∑ (2.2) which proves that time-step has no effect on the model structure. I selected time series models to prove the hypothesis as time series represents the simplest forecast relationship (i.e., one-dimensional). Time series models are sensitive to the fulfillment of the statistical assumptions (particularly independence), therefore to ensure that the results are not biased (Gujarati 1995) or meaningless, I considered that the analysis was completed when the residuals were reduced to white noise (WN) and did not exhibit nonlinear dependencies (Brockwell and Davis 1996). To prove the validity of the argument (i.e., the relationship between the environment components are described by equations that are simultaneously linear and nonlinear) only empirical equations were considered, as complex environmental models contain a significant number of empirically-derived equations (Voinov et al. 2004). The proof used data from the oil and gas industry and developed autoregressive equations connecting the count of active petroleum rigs, a variable representing the complex array of environmental conditions required for drilling, with the New York Mercantile Exchange (NYMEX) oil price (Baker Hughes Inc 2005; Oil and Gas Commission 2005). To investigate the time step influence on the relationship between the number of active petroleum rigs and oil prices I adopted a hybrid approach, in the sense of Walls (1992), which covers both econometrical and geological investigations. Following the recommendation of Walls (1992), the relationship was assessed using the most common variables present in petroleum drilling inquiries: actual (xt1) and historical average (xt2)of the number of rigs and NYMEX crude oil prices at time t (i.e., monthly minimum (xt3), maximum (xt4) and average (xt5)). To demonstrate the hypothesis I used monthly averages for five time-steps (i.e., 1, 2, 3, 6 and 12 months). The monthly average was recommended by the length of time for which an active rig is stationed at a particular drilling location (i.e., commonly four to five weeks). To test whether the argument is valid at different spatial scales, jurisdictions and geological conditions I examined the dynamics of the relationship between oil prices and the number of active rigs for the main petroleum production areas around the world: Africa, Canada, Europe, Far East, Latin America, Middle East and USA. The results (Appendix A) confirmed the findings of previous studies (Walls 1992) and identified a significant relationship between the number of active rigs on one hand, and the past number 14 of rigs and the NYMEX oil prices, on the other. However, the time-step changed the structure of the relationship for all regions except the USA, as the NYMEX oil prices were eliminated from the models based on averages computed for time-steps larger than three months (i.e., significance larger than 0.1). For the USA, the equation had non white-noise residuals (p=0.001). Consequently, for the USA, I built an autoregressive model (AR) that met modeling requirements (p>0.06 for WN and p>0.05 for quadratic dependency) but eliminated all the NYMEX oil prices, regardless of the time-step, as being highly insignificant (p>0.5). The results suggest that the number of active rigs is sensitive to oil prices in all regions except the USA, where variables other than NYMEX oil prices seem to play a more important role in the dynamics of active rigs. The lack of significance of the NYMEX oil prices for time steps greater than three months could be related to the loss of information associated with the increase in the time-step size used to compute the monthly average (Kuo et al. 1999). Therefore, generating structurally invariant models for all time-steps (Refsgaard 1997; Schoorl et al. 2000) led to inefficient results. While it was expected that by increasing the temporal scale of the investigation all the oil prices would lose their importance, it was not expected that for the same region different time-steps will contain different oil-prices (e.g., the model for Europe included NYMEX maximum oil price for the monthly time step and the minimum oil price for the tri- monthly data). This change of the variable (i.e., NYMEX oil prices) within the models (i.e., regional equation) only indicates a lack of structural invariance, but evidence of the absence of structural invariance was supplied by the violation of the Banach-Steinhaus theorem (Banach and Steinhaus 1927) by the autoregressive and the moving average terms of the model (Appendix A). Banach and Steinhaus (1927) proved that the linear equations will converge to the true relationship, in respect to the time-step, provided that there is no evidence that the linear equations are incorrect. The theorem ensures that the linear model ∑ ≥ +== 1 |||0|1 ' i p TiTiT p TTT xaaXAx (2.3) converges to the true model XA' only if pTix | constitute a bounded fundamental sequence on a complete set (where AT’ is the transposed vector of the coefficients ai|T, p is the autoregressive order, p≥1, and XTp is the (xpi|T ) vector with time step T). The convergence is trivial for time steps larger than a critical value, Tc, for which the correlogram contains only terms that are not significantly different from 0 as 0''lim aXAXA TTTcT ==> but it is not necessarily valid for all T < Tc as iTiTcT aa ≠< |lim . The latter case characterizes the models for Africa and Europe (Figures 2.1a and 2.1b), where it seems that 0lim | →→ TiTcT a , or 8.0lim | <→ TiTcT a respectively, even that the 15 expected value is 1| =TiEa (i.e., the time series became random noise for T>Tc). Therefore, the linearity of a variable significant at all scales is either nonexistent or the analysis is data dependent, as fulfillment of the Banach-Steinhaus theorem is contingent on data (i.e., not valid for Europe or Africa). Furthermore, all the models had a break or a significant drop in the autoregressive coefficient at the two-month time-step, indicating a haphazard connection between the theoretical findings and the real data. Consequently, besides theoretical considerations strictly related to the model development, the selection of the model depended on the time-step. However, for each time-step there is no indication that the models are incorrect (i.e., all assumptions and requirements were met) (Appendix A), leading to the conclusion that there is no structural invariance across time-steps. Figure 2.1 Coefficient of the autoregressive term of the number of active rigs and the dynamics of the moving average for the error model. The drop in the magnitude of the coefficient at the time step two-months and the variation in the type of equation describing the error model show the haphazard investigation of the complex environmental system, while the decrease of the coefficient of the autoregressive term with time-step contradicts the Banach-Steinhaus theorem. The lack of consistency across different time-steps for the significance of the variables, and the absence of convergence of the coefficient of the autoregressive term to the expected value, indicated the existence of an insufficient relationship between the theoretical framework and data obtained from observations. This insufficiency was deepened by the fluctuation of the model describing the moving average (MA) part of the relationship. The use of the MA or the 16 autoregressive conditional heteroskedasticity (ARCH) process (Brockwell and Davis 1996) was imposed to ensure that the residuals were organized as WN and were linearly and quadratically independent (Appendix A), a condition essential for valid time-series inferences. A loss of information was present to the MA and ARCH terms, as for time-steps less than three months, 90% of the models had a MA or ARCH term, compared to only 47% for time-steps larger than three months. Additionally, the type of equation describing the moving average depended on the time step (Appendix A), as all the models, except for the USA, changed the degree of the MA equation at least once (Fig 2.1c and 2.1d): from linear to no MA term or the opposite (i.e., Latin America, Canada and Europe), from linear to quadratic and no MA (i.e., Africa and Far East) or from quadratic and no MA (the case of Middle East). The variation of the MA part of the model changes the entire model from linear to nonlinear and back. The same haphazard variation as identified for the autoregressive coefficients was present in the model type, confirming that the equations describing the active numbers of rigs is linear and nonlinear, according to the time- step, consequently the absence of an invariant structure across different time-steps. I complemented the previous analysis performed on areas covering millions of square kilometers, different geo-climatic zones and politico-economical regimes, with a similar investigation but focused on one homogeneous area, northeastern British Columbia (BC), Canada. To substantiate that structural invariance is not confined to the variable representing the count of active rigs, the BC analysis considered the number of wells drilled. To enhance the investigation performed on the number of active rigs, besides the linear model a Gamma distribution model, suitable for count type data, was used to describe the dynamics of the number of well existing on the landscape (Renshaw 1994). The results (Appendix A) showed that regardless of the model type (i.e., linear or Gamma) the smaller and geologically and regulatory homogeneous area is characterized by models with greater structural variability than the larger and more heterogeneous regions. However, the main distinction between the BC and the continental or national models was the absence of a linear model across all time-steps (Fig. 2.2a). The presence of a square root model for the number of wells drilled annually, which replaced the linear model, was imposed by the necessity to have the residuals distributed as WN. For the annual time-step, not only the linear model disappeared but the significant variable also changed, from the annual to the cumulative number of wells drilled (Appendix A). The transformation of the model with each time-step was amplified for the Gamma distribution model, which had a linear denominator only for monthly and bimonthly data and a square root denominator for the remaining time-steps (Fig.2.2 b). 17 Nevertheless, for both the linear and Gamma models, the changes of the significant variables from linear equation to square-root equation did not affect the error term, which had a linear MA regardless of the time-step. The loss of information identified at the models describing the number of active rigs was also observed for BC as NYMEX oil prices were dropped from the models with time-steps larger than two months (Fig 2.2.c and 2.2 d). The same conclusion regarding the lack of structural invariance found for the active number of rigs was reached for the models representing the number of wells, as for each time step and model (i.e., linear and Gamma) there were no indications that the models were incorrect but the evolution of their structure violated the Banach – Steinhaus theorem. Therefore, either the models were nonlinear or the analysis should have included information related to the analysis itself, such as time- steps or scale. However, there was no indication that the linear models were incorrect, leading to the conclusion that the models should incorporate in addition to the mathematical formula information regarding their analytical validity. Figure 2.2 Confidence interval (α=0.02) for coefficient of the autoregressive term (i.e., number of wells drilled per time-step), and the number of economic variables present in the model describing the total number of wells drilled per time-step. The same pattern identified for the rig count model (i.e., decrease in the coefficient of the autoregressive term with the time-step) and the absence of a linear model across all time steps confirm the haphazard character of a scientific investigation. The decrease in the number of economic variables confirms the loss of information mentioned by other studies (Kuo et al. 1999) when time step increases; therefore the inefficiency of structural invariant models. 18 The results demonstrate that the development of complex representations of the environment could lead to models whose variable structure denies the convergence to the true model, regardless of the model type (i.e., linear or nonlinear) and distribution used (i.e., normal or gamma). Therefore, it seems that data obtained through observations or experiments are insufficient to quantify complex environmental relationships. The analysis itself could play a significant role in the investigation process as it is structurally dependent on at least one external factor, namely scale. One must therefore either acknowledge the inability of current scientific methods to reveal the true relationships governing complex environmental systems or should introduce an analytical factor into the description of those systems. 19 2.1 References 1. Bacon, F. 1855. The New Organon. Oxford University Press, Oxford. 338 p. 2. Baker Hughes Inc. International Rig Count. 2005. Baker Hughes Inc. 3. Banach, S. and H. Steinhaus. 1927. Sur le principe de la condensation de singularités. Fundamenta mathematicae. 9:50-61. 4. Brockwell, P.J. and R.A. Davis. 1996. An Introduction to Time Series and Forecasting. Springer Verlag, New York.420 p 5. Chaubey, I., A.S. Cotter, T.A. Costello, and T.S. Soerens. 2005. Effect of DEM data resolution on SWAT output uncertainty. Hydrological Processes 19(3):621-628. 6. Enquist, B.J. and K.J. Niklas. 2001. Invariant scaling relations across tree-dominated communities. Nature 410(6829):655-660. 7. Grimmett, G.D. and D.R. Stirzaker. 2002. Probability and Random Processes. Oxford University Press, New York.1-600 8. Gujarati, D.N. 1995. Basic Econometrics. McGraw-Hill Companies, New York, NY.1-350 9. Kimmins, J.P., D. Mailly, and B. Seely. 1999. Modelling forest ecosystem net primary production: the hybrid simulation approach used in FORECAST. Ecological Modelling 122(3):195-224. 10. Korzukhin, M.D., M.T. TerMikaelian, and R.G. Wagner. 1996. Process versus empirical models: Which approach for forest ecosystem management? Canadian Journal of Forest Research-Revue Canadienne de Recherche Forestiere 26(5):879-887. 11. Kuo, W.L., T.S. Steenhuis, C.E. McCulloch, C.L. Mohler, D.A. Weinstein, S.D. DeGloria, and D.P. Swaney. 1999. Effect of grid size on runoff and soil moisture for a variable-source-area hydrology model. Water Resources Research 35(11):3419- 3428. 12. May, R.M. 1974. Biological Populations with Nonoverlapping Generations - Stable Points, Stable Cycles, and Chaos. Science 186(4164):645-647. 13. Mills, E. 1998. A Simple Solution to the Liar. Philosophical Studies 89(2):197-212. 14. Oil and Gas Commission. Drilling Data for all wells in B.C. 2005. Oil and Gas Commission. 15. Oreskes, N., K. Shrader-Frechette, and K. Belitz. 1994. Verification, Validation, and Confirmation of Numerical-Models in the Earth-Sciences. Science 263(5147):641- 646. 20 16. Refsgaard, J.C. 1997. Parameterisation, calibration and validation of distributed hydrological models. Journal of Hydrology 198(1-4):69-97. 17. Renshaw, A.E. 1994. Modeling the Claims Process in the Presence of Covariates. ASTIN Bulletin 24(2):265-285. 18. Schoorl, J.M., M.P.W. Sonneveld, and A. Veldkamp. 2000. Three-dimensional landscape process modelling: The effect of DEM resolution. Earth Surface Processes and Landforms 25(9):1025-1034. 19. van Fraassen, B.C. 1980. The Scientific Image. Oxford University Press, New York NY. 248 p 20. Vazquez, R.F., L. Feyen, J. Feyen, and J.C. Refsgaard. 2002. Effect of grid size on effective parameters and model performance of the MIKE-SHE code. Hydrological Processes 16(2):355-372. 21. Verburg, P.H., P.P. Schot, M.J. Dijst, and A. Veldkamp. 2004. Land use change modelling: current practice and research priorities. GeoJournal 61(4):309-324. 22. Voinov, A., C. Fitz, R. Boumans, and R. Costanza. 2004. Modular ecosystem modeling. Environmental modelling and software 19:285-304. 23. Walls, M.A. 1992. Modeling and forecasting the supply of oil and gas : A survey of existing approaches. Resources and Energy 14(3):287-309. 24. West, G.B. and J.H. Brown. 2004. Life's universal scaling laws. Physics Today 57(9):36- 42. 25. West, G.B., J.H. Brown, and B.J. Enquist. 1997. A general model for the origin of allometric scaling laws in biology. Science 276(5309):122-126. 26. Wigmosta, M.S., L.W. Vail, and D.P. Lettenmaier. 1994. A Distributed Hydrology- Vegetation Model for Complex Terrain. Water Resources Research 30(6):1665- 1679. 27. Winston, W.L. 1994. Operations Research: Applications and Algorithms. Duxbury Press, Belmont. 1372 p 21 3. A DETERMINISTIC HARVEST SCHEDULER USING PERFECT BIN PACKING THEOREM1 3.1 Introduction To incorporate sustainable forest management into annual allowable cut determinations, forest planners are using techniques that have evolved from one-dimensional models, such as the area- or volume-control approach (Davis et al. 2001, p. 528-545; Bettinger et al. 2009, p. 214- 228) to complex multidimensional perspectives, such as mixed integer programming and heuristic techniques (Bettinger et al. 2002). Mixed integer programming can identify the optimal solution that answers the planning requirements, but has limited applicability in real forest planning applications as it depends on the problem size (Crowe et al. 2003). Sometimes even problems with less than 100 variables can be difficult to solve using mixed integer programming (Anderson et al. 1994, p.345), as the computation time can have an exponential or even factorial dependency on the number of integer variables (Nemhauser and Wolsey 1988; Wolsey 1998). Additionally, scheduling is complicated by the interdependency between the temporal constraints (e.g., each stand has a different harvesting age) and the spatial constraints (e.g., the size of an opening has to be smaller than 60 ha or the forested landscape has to contain wildlife corridors). The introduction of multidimensional approaches into forest planning led to the widespread use of heuristic techniques, techniques with reduced complexity comparing to mixed integer programming but supplying suboptimal results (Zomaya and Kazman 1999). The lack of optimality of heuristic algorithms is present both at stand and forest levels. From a stand- level perspective, the departure from optimality is associated with harvesting at ages different to those determined through marginal analysis (Duerr 1993, p. 93; Amacher et al. 2009, p.38, 47). At the forest level, the sub-optimality of the heuristic techniques follows from the fact that the annual allowable cut (AAC) determined using heuristic algorithms is less than the AAC supplied by mixed integer programming [sometimes the difference is less than 5%, as reported by Boston and Bettinger (1999)], when similar planning constraints are imposed (Crowe and Nelson 2003; Pukkala and Kurttila 2005; Bettinger et al. 2002; Murray and Church 1995). The objective of forest planning (i.e., scheduling the harvestings in time and space) is complicated by the type and spatial repartition of the silvicultural systems distributed across the 1 A version of this chapter has been submitted for publication. Strimbu, B.M., Innes, J.L. and Strimbu, V.F. A deterministic harvest scheduler using perfect bin packing theorem. 22 forest (Heinonen and Pukkala 2004; Baskent and Jordan 2002), which can lead to non- polynomial complexity problems (Kangas et al. 2002). Clearcuts are operationally preferred and are easier to implement in computations than selection or shelterwood methods. However, clearcuts introduce a set of constraints that are not present with other silvicultural systems within an integrated forest management framework: the greenup/adjacency delay. This type of constraint further increases the complexity of the planning problem, promoting even more the use of heuristic techniques in solving the large combinatorial scheduling problems, as defined by Lockwood and Moore (1993). Among the heuristic techniques used in forest planning, stochastic algorithms are the most popular and it is common practice to perform a set of runs with the “best” solution being selected (Bettinger et al. 2002). The necessity of several runs complicates the investigation and diminishes the reduced computation time feature of the heuristic algorithms, as the selected solution can be obtained after several hours of total computing time (Murray and Church 1995; Pukkala and Kurttila 2005). In the case of large scheduling problems, the performances of the heuristic technique are assessed by comparing the magnitude of the AAC supplied by the heuristic algorithms with a value considered an upper bound (Bettinger et al 2009), commonly obtained using linear programming (LP). The present research focuses in determining a majorant for the AAC, in the sense of Shanks and Gambill (1973) that accommodates simultaneously the optimality for stand and forest and provides a tighter bound than LP for the AAC. Additionally, I propose a new algorithm for computing the AAC that is deterministic and supplies results close to the value that optimizes the objectives for both stand and forest. 3.2 Methods 3.2.1 Data description The performance of a new harvest scheduler algorithm or the adjustments of an existing one are normally assessed using either computer-generated data (Bettinger et al. 2002) or real data (Murray and Church 1995). Following Heinonen and Pukkala (2004), who used real data to evaluate two harvest scheduling algorithms, the new AAC majorant and the new harvesting algorithm are presented using forest inventory data from three areas from north-eastern British Columbia, Canada (Fig 3.1). 23 Figure 3.1 Study areas Each area has more than 35,000 stands [in the sense of Kangas and Kangas (1999)] larger than 0.5 ha (Table 3.1) and combined, the three areas cover 1.75 million ha. Depending on the area, the timber harvesting land base (THLB) contains between 4700 and 8000 units larger than 5 ha (the units were obtained by merging the original adjacent stands smaller than 5 ha based on species, site index and age) (Table 3.1). The data were supplied by the former British Columbia Ministry of Forests and the former British Columbia Ministry of Sustainable Resources Management. Three biogeoclimatic zones (Anonymous 1997a; Anonymous 1997b; Anonymous 1999) cover the areas considered in the analysis: the Boreal White-Black Spruce (BWBS) zone, which dominates the western extension of the Great Plains region of Canada , and the Spruce– Willow–Birch (SWB) and Engelmann Spruce–Subalpine Fir (ESSF) zones that are found at elevations above 1300 m. The main tree species are white spruce (Picea glauca), black spruce (Picea mariana), trembling aspen (Populus tremuloides), lodgepole pine (Pinus contorta), balsam poplar (Populus balsamifera), Engelmann spruce (Picea engelmannii) and subalpine fir (Abies lasiocarpa). 24 Table 3.1. Timber harvesting land base (THLB) Study Area Forest cover THLB Surface [ha] # stands (> 0.5 ha) Surface [ha] # units >5 ha obtained by aggregation of adjacent stands smaller than 5 ha Doig 684 138 37042 183 057 7926 Fort Nelson 641 024 86893 91 355 5486 Moberly 410 194 36319 140 222 4729 3.2.2 Optimal harvesting age and mean annual increment The moment when a stand is harvested is determined by the dynamics of the attributes considered in the planning process, such as biomass, net present value or habitat suitability. The set of attributes describing a stands can be represented as a multivariate function f : Rm → R, ))(),..,((),..,( 11 iyiyfyyf mm = , where )(),..,(1 iyiy m are the functions quantifying the attributes y1 to ym at age i into the real set R. The harvest timing ensuring stand – level optimality occurs at the age that maximizes the function f (when optimality is represented by minimum then the inverse transformation would revert the analysis to maximum) during an undefined long planning horizon (i.e., an unspecified large number of harvests); that is ∫ horizonPlanning m dttytyf 0 1 ))(),..,(( . In the case when silvicultural treatments are repeated indefinitely ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ×=⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ∫∫ RotationhorizonPlanningdttytyfdttytyf Rotation m horizonPlanning m 0 1 0 1 ))(),..,((max))(),..,((max (3.1) The optimal harvesting age (OHA) is the age for which the Eq. 3.1 holds irrespective the length of the planning horizon. According to Eq. 3.1 the OHA is the rotation that ensures the maximization of the function f during one rotation: ⎟⎟ ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎜ ⎝ ⎛ ∫ Rotation dttytyf Rotation m 0 1 ))(),..,(( max . The maximum of Rotation dttytyf Rotation m∫ 0 1 ))(),..,(( is reached when its total derivative in respect to the rotation R of the stand is null (Shanks and Gambill, 1973): 25 0)( ))(),..,(( )( ))(),..,(( ))(),..,(( 1 0 1 0 1 0 1 =∂ ∂×∂ ∂+∂ ∂= = ∑ ∫∫ ∫ = R Ry R dttytyf RyR dttytyf R R dttytyf dR d jn j R m j R m R m (3.2) The differentiation of the Eq. 3.2, which considers a broad array of attributes, leads to a complicated equation that could occlude the objective of the current research (i.e., determination of a tighter bound for the harvested volume than linear programming and presentation of a new forest harvesting algorithm). Furthermore, a series of functions quantifying some of the stand’s attributes could contain a significant subjective component (such as visual quality objectives or representation of the wildlife habitat loss in monetary terms) that would increase unnecessarily the complexity of the presentation. Without a significant loss of generality (understood as the unaltered expansion of findings between two dimensional spaces) I use only one attribute (i.e., the stand’s yield), which allows a clear description of the algorithm and an objective quantification of the stand. When the yield is assumed to be harvested only once during the rotation (i.e., clearcut), then the function f reduces to one yj(i) and the Eq. 3.2 has the simpler form 0)( = R Ry dR d , where y(R)= yj(R). The R Ry )( is the mean annual increment (MAI) corresponding to the rotation; therefore, for the simple case of clearcutting, the optimal harvesting age occurs at maximum MAI, max(MAI), when only stand’s yield is considered, which confirms the results of Rucareanu and Leahu (1982, p. 127-140), Davis et al (2001, p. 146) and Bettinger et al (2009, p.108). Harvesting at max(MAI) ensures the optimality of the timber accumulation criterion (Davis et al. 2001, p. 145-147) at the stand but not at the forest level. To integrate the stand within the forest, current schedulers optimize the timber harvesting for the entire forest, so that most stands are harvested at ages different than the age at max(MAI), Age max(MAI), therefore at a suboptimal value for the stand. In the case of harvesting at suboptimal ages, it is preferable for harvesting to occur after the MAI has peaked, as the annual stand growth rate has already passed its maximum and the stand is likely to have reached a steady-state phase (Nyland 1996, p. 204). 26 3.2.3 Maximum volume harvested annually determined using MAI A perfect harvest scheduling model is a model that guarantees the optimality of the established goals at the stand and at the forest level. Therefore, when the goal is the maximization of the even-flow annual harvested volume, a perfect scheduler should harvest each stand at OHA [i.e., the age of max(MAI)] and should ensure the constancy and global maximality of the annual harvested volume. The identification of a solution that accommodates both optima, for the stand and for the forest, could start by establishing the OHA for each stand, which is determined by the merchantability standards when only the volume is of interest. Following the establishment of OHA for each stand the maximum volume to be harvested annually (MVHA) can be determined using the perfect bin-packing theorem (Coffman et al. 2000). The perfect bin-packing theorem (PBPT) states that “for positive integers k, j, and r, with k ≥ j, one can perfectly pack a list L consisting of r×j items, r each of sizes 1 through j, into bins of size k if and only if the sum of the r×j item sizes is a multiple of k” (Coffman et al. 2000). Therefore, a necessary and sufficient condition for a constant MVHA (i.e., the bin size from the PBPT) of n stands (i.e., the list from PBPT) is that the volume of the harvested stands can be organized as an arithmetic series (i.e., the 1 through j from the PBPT) whose terms are uniformly distributed across the series (i.e. the r from the PBPT) and having the sum a multiplier of MVHA. In eventuality that the PBPT conditions are met, then is a multiple of the MVHA. Two constraints are imposed by the perfect bin-packing theorem: first, the volumes of the harvested stands range from 1 through j without gaps, and, second, that the forest contain a constant number of stands (i.e., r from the PBPT) with the same volume (i.e., r of volume 1, r of volume 2 and so on). Usually, real forests do not meet these conditions as the volume of the harvested stands are seldom evenly distributed across an arithmetic series with a step of one, even when a normal forest was the long term goal of the forest management (Rucareanu and Leahu 1982). The volume distribution of the harvested stands can be described by continuous and invertible cumulative functions, which can be translated into uniform distributions, as proved by the inversion principle (Devroye 1986, p.28). For the main functions used to describe the volume distribution, Deveroye (2003, p.29) provided the inversion function transforming the non- uniform random variate into a uniform variate, ensuring the fulfillment of the PBPT. The simplest 2 )1( 1 +=×∑ = jjrir j i 27 case of transforming the series of the harvested volume in an arithmetic series with step one and the term uniformly distributed across the series is when the harvested volumes have a uniform distribution. In this case, r from the PBPT would be the number of stands with the same volume (the case of discrete distributions) or within the same volume class (the case of continuous distributions), j= n/r, and each harvested volume would be represented as a value between 1 and j if the following linear transformation is applied: ⎥⎦ ⎤⎢⎣ ⎡ − ×−− − )//()( )101( minmax minstand rnVV VVceiling ε (3.3) where Vmin is the smallest harvested volume, Vmax the largest harvested volume from the series of stands with volume when harvested V stand and ε is a natural number ensuring that the ceiling function will return value 1when Vstand = Vmin and value j when the Vstand = Vmax. Perfect bin-packing theorem identified the necessary and sufficient conditions for perfect scheduling (Coffman et al. 2000), but supplied no indications regarding the size of MVHA. The magnitude of the MVHA was determined on the next corollary of PBPT: Corollary. For a set of merchantability standards ensuring the fulfillment of the PBPT, the maximum volume that can be harvested annually is equal to the sum of the maximum mean annual increment of the stands. (A proof is presented in the Appendix B). To meet the requirements of the PBPT, the volume of the stands should be uniformly distributed or have a left-skewed distribution with respect to the harvest timing, which is the difference between the OHA and the current age of the stand. The volume distribution with respect to harvest timing depends on the conditions of the forest and the selected set of merchantability standards. The conditions of the forest cannot be changed, therefore when the distributions are skewed to the right, the OHA have to be recalculated such that the volume distribution meets the PBPT. The OHA ensuring the PBPT requirements are determined by adjusting the merchantability standards [e.g., from minimum upper stem diameter 10 cm (the stand could be managed for pulpwood) to 20 cm (the stand could be managed for sawtimber)], as according to Eq. 3.2 the OHA depends on the derivative of the function quantifying the attributes of interest (e.g., MAI for pulpwood differs from MAI for sawtimber). The new merchantability standards can be found by solving the following mixed integer programming problem, which would transform a right-skewed distribution to a uniform or slightly left-skewed distribution: 28 agestandpresentagestandpresentimingharvest t moment)tingmax(harvesand1between', standeachfor1 }1,0{ 'thatsuch,'and,moments harvestingtwoanyfor 1thanlargeryearanyfor)max( tosubject )max(max )max( |, '|, , , , −=−= ≤ ∈ <×≥× ×≥× × ∑ ∑ ∑ ∑ ∑ ∑ g MAI g OHA g g i g i ygi ygi g i g i g i g i gi gi g i g i g i g i gi g i g i AgeAge yy iX X yyyyVXVX yMAIXVX MAIX (3.4) where Age gOHA = Age gmax(MAI) is the OHA determined for the goal g (i.e. the age at max(MAI) determined for the merchantability standard g), g iX is the binary variable indicating whether or not stand i would be harvested at Age g max(MAI), g iV is the volume of stand i at OHA for goal g. The objective function ensures that among all possible merchantability standards (e.g., the combination between the minimum diameter at breast height and the minimum upper stem diameter) only the standards that maximize the Corollary’s results are selected. The first constraint selects the combinations of merchantability standards leading to even-flow, while the second constraint restricts the harvestings only to the merchantability standards supplying uniform or left-skew distributions of the harvested volumes. The last two constraints ensure that each stand has one (the third constraint) and only one (the fourth constraint) merchantability standard (the initial or the adjusted standard). 3.2.4 Adjusting Maximum Volume to be Harvested Annually: First Fit Decreasing algorithm The computation time used to solve forest scheduling problems depends on the problem size, importance of the results and the CPU performances, with several authors providing solution in seconds (McDill et al 2002), minutes (Gun and Richards 2005; Bettinger et al 2002) or hours (Crowe et al 2003). The present research focuses on choosing an efficient and deterministic algorithm that provides a solution to the scheduling of the harvest of more than 4000 stands in 29 less than 10 minutes. One possible such algorithm is the First Fit Decreasing (FFD) algorithm (Coffman et al., 1997) that can be implemented to run on average )log(nnR ×× steps, less than 1 min when rotation R is 100 years, the number of stands n is less than 10000 and the CPU is an Intel 3 GHz Pentium 4. The FFD algorithm sorts the stands in decreasing order in relation to the attributes of interest and allocates stands to the nearest year in the planning period that can accommodate its size (Coffman et al. 1997). Johnson (1973) showed that queuing the stands in relation to the attribute quantifying the planning objective (e.g., volume when the objective is even – flow maximization) supplies results that have an expectation larger than 7/9×MVHA [the result without proof can be found in Coffman et al. (1997)]. The use of a different criterion in queuing and selecting the harvest timing of a stand (such as oldest first or randomness) could lead to results with an expectation 70% smaller than the optimal solution (Garey et al. 1976; Ullman 1971). Therefore, I built an efficient algorithm (Fig 3.2) which, by combining the FFD method with the distributional adjustment required by PBPT (i.e., the harvested volumes are uniformly distributed in respect to the harvest timing), ensured that the even – flow of the AAC is larger than 7/9×MVHA. The spatial constraints associated with sustainable forest management, such as greenup – adjacency delays or wildlife corridors, were not considered in MVHA computations and could lead to results below 7/9×MHVA. The spatial restrictions were incorporated into the FFD algorithm using the unit restriction model, in the sense of Murray (1999). The units are either the individual stands larger than 5 ha or an aggregation of adjacent stands smaller than 5 ha, aggregation based on the harvest timing. The selection of unit restricted model to address the greenup-adjacency was preferred over the area restricted model, as the area restriction model is computationally intense (Goycoolea et al 2005), and the purpose of this chapter is to present the algorithm itself not the response of the algorithm to the modeling of the spatial constraints. The green-up adjacency constraint complicates the planning problem but has little effect on the magnitude of MVHA determined based on PBPT corollary when the difference between harvest timing of neighboring stands is larger than the greenup – adjacency delay. When the availability of stands for harvesting is limited by the adjacency delay then either the harvesting age could be determined by adjusting FFD [i.e., the stand with largest volume and age larger or equal to OHA] or by adding extra constraints to the mixed integer programming problem (3.4). Adjusting the harvesting age using mixed integer programming problem (3.4) could lead to infeasible matrices, from the size perspective (matrices larger than 60000 x 60000 values). Therefore, to avoid size-related issues and to maintain the computational time to minutes, the harvest timing of a stand fulfilling the greenup-adjacency delay was established by adjusting the FFD 30 algorithm, such that the stand with the largest volume would be harvested first (condition imposed by the FFD algorithm) but at any age equal or larger than OHA, which relaxes the PBPT’s Corollary [i.e., age at max(MAI)]. The adjustment of the FFD algorithm was preferred over the mixed integer programming (3.4) as bin-packing heuristic algorithms are self organizing to resemble the PBPT requirements (Csirik et al. 1999) and can supply results having the expectation MVHA (Csirik and Johnson 2001). Following the work of Murray and Church (1995), Bettinger et al. (2002), Crowe and Nelson (2003), Pukkala and Kurttila (2005), I compared the results supplied by the adjusted FFD with the results produced by simulated annealing (SA) and a relaxed LP, using data from the three areas. The set of constraints framing the scheduling problem were similar to those used by Crowe et al. (2003) but enhanced with some of the conditions from the Sustainable Forest Management Plan developed for the Fort St. John Pilot Project (Jukes et al. 2004): A maximum 10% fluctuation of the annual harvested volume in respect with the average volume harvested yearly during the planning period; A maximum opening size of 60 ha; A minimum opening size of 5 ha; Stands are harvested either at ages larger than a preset value, as suggested by Bettinger et al. (2002) and Crowe and Nelson (2003), which for the study areas was recommended by Jukes et al (2004) to be approximately 120 years, or at ages larger or equal to OHA (i.e., Agemax(MAI)). The optimal harvesting age (Appendix C) was determined considering the merchantability standards of 17.5 cm for spruce and sub- alpine fir and 12.5 cm for all other species (Pedersen 2003); Greenup/adjacency delay of 1, 5,10 and 20 years, respectively; A planning period of 100 years. The above constraints determined four comparison scenarios [i.e., two algorithms (i.e., SA and adjusted FFD) × two harvesting ages (i.e., harvesting at OHA or older and harvesting above a minimum age)]. The relaxed LP was not included in the comparisons as it was used mainly to compute a majorant for the SA algorithm. 31 Figure 3.2 First fit decreasing scheduler algorithm that fulfils the PBPT corollary requirements and queues the stands according to the attribute considered for optimization 32 3.2.5 Simulated annealing Simulated annealing (Metropolis et al. 1953) has been identified as one of the heuristic algorithms that supplies results close to the linear programming solution attached to the forest planning problem (Crowe and Nelson 2003; Bettinger et al. 2002). The simple implementation of the simulated annealing (SA) has led to its frequent use as a test algorithm in forest planning. Additionally, SA supplies one of the fastest solutions amongst commonly used heuristic algorithms, as for similar problems tabu search requires three times more time and genetic algorithms require five times more time (Bettinger et al. 2002). The parameters of SA [i.e., initial temperature (Ti), freezing temperature (Tf) and the rate of annealing (ra)] are usually selected experimentally for each problem (Boston and Bettinger 1999; Bettinger et al. 2002; and Crowe and Nelson 2003). The tree parameters, together with the number of solutions computed for each temperature, S, determine the run-time of the algorithm if no stopping rule is enforced [e.g., the solution did not change for a certain ratio between the number of accepted and rejected solutions (Zomaya and Kazman 1999; Baskent and Jordan 2002) or the solution is the best solution obtained after a preset number of iterations (Crowe and Nelson 2003)]. The runtime needed by SA to reach an optimal solution with certainty is greater than the runtime required for complete enumeration and evaluation of all possible solutions (Geman and Geman 1984); which makes the guaranteed optimality pointless with SA. An estimation of the time to obtain SA solution using an exponential cooling schedule (Nourani and Andresen 1998) is: (3.5) Based on Eq. 3.5, the expected time required by SA to reach an optimal solution with the parameters (Ti, Tf, ra) is ∫ ∫ Δ Δ ΔΔ≈>×≈ ≈×××+××= )max( )min( )()/exp()]1,0([ )()()( V V T T f i dTVdTVSNtUpESNt ptNSEtNSEtimeE (3.6) where N is the number of temperatures )log( )log()log( a if r TT N −= p is the probability of accepting a new solution, V’’, worse than the current solution, V’ at temperature T: =exp(ΔV/T) )/()(#)/(# solutiontimestempertureetemperatursolutionstime ××= ]/)"'exp[( TVVp −= 33 t is the time to generate and assess and alternate state S is the number of solutions for each temperature U(0,1) is a uniform random variate taking values between 0 and 1. When actual harvest scheduling values are used in the SA computations, the exponential integral ∫ −− dxxx )exp(1 can be assumed negligible in respect to )/exp( TVT Δ ; therefore, the double integral from Eq. 3.6 )/exp()exp()/exp()()/exp( 1 TVTdxxxTVTdTVdTV Δ≈−+Δ=ΔΔ ∫∫∫ − , and consequently ]/)exp[max( ))min()max(exp)min()max(exp[)( ff i i f f TVSNtT T VVT T VVTSNttimeE Δ×≈ ≈Δ−Δ×−Δ−Δ××≈ (3.7) To identify the set of parameters leading to solutions close to optimality I performed a preliminary investigation for each area by running the SA only 15 min and choosing seven rates of annealing, ranging between 0.70 and 0.99 (i.e., from 0.7 to 0.95 in steps of 0.05 plus 0.99) and 13 initial temperatures (103, 104, 105, and from 106 to 108 in steps of 9.9 x 106). Using several trial and error runs I selected as freezing temperature the value 100, except for initial temperatures 106 or less when 0 was used. The number of solutions at any given temperature was determined experimentally starting with the value that supplied the ratio between the number of accepted and rejected solutions (Hoss and Stutzle, 2005, p.77) of 1:10, as proposed by Zomaya and Kazman (1999). For each study area, 10000 solutions per temperature led to an average ratio of maximum one accepted to 10 rejected solutions, resembling the results of Bettinger et al (2002). Consequently, S = 10000 was used in all the computations, irrespective the initial and freezing temperatures, annealing rate or study area. Based on Eq. 3.7, the selected set of parameters and assuming that the time to evaluate a solution is approximately 10-10 seconds, the expected computing time is of the order of 1018 seconds (approximately 1013 hours), which lead to the termination of the algorithm not when the temperature would reach the cooling point but when the ratio between the number of accepted and rejected solution is smaller than the accepted ratio, as recommended by Zomaya and Kazman (1999). Hastings (1970) linked SA to Markov chains; therefore, when the time of a run is relatively short, the generation of several runs is required. Consequently, for green-up/adjacency of 5 years I produced ten solutions for each initial temperature, rate of annealing, study area and method of selecting the harvesting age (i.e., minimum age or OHA). Consequently, I performed a total of 34 10 runs × 13 initial temperatures × 7 annealing rates ×3 study areas × 2 harvesting ages = 5460 SA runs with a total time of 1365 hours. To investigate the impact of the adjacency delay on the performance of SA, an additional set of runs, each one hour long, were executed. For each study area, harvesting age and adjacency delay and using only the set of parameters (i.e., initial temperature, number of iterations per temperature and annealing rate) that supplied the largest results among the above 5460 runs, a series of five runs were executed, a total of 120 runs. The one hour run supplied larger results than the 15 min run and was comparable with the waiting time of Crowe and Nelson (2003), which allowed 40 min for SA to supply a solution. The greenup-adjacency constraint in SA was addressed with the same algorithm used for the adjusted FFD, the unit restriction model. The results supplied by the two algorithms were compared with the solution supplied by a relaxed linear programming (LP) representation of the scheduling problem (Eq. 3.8). The employment of a relaxed LP in conjunction with different heuristic techniques was advocated by the Lockwood and Moore (1993), Csirik et al (1999), Baskent and Jordan (2002), Bettinger et al (2002, 2009) and Gunn and Richards (2005), who used it to identify the bounds of the optimal solution. The LP constraints relaxed were the greenup-adjacency delay and the harvesting of a stand in one intervention (i.e., the harvest of a stand could occur in several clearcuts), which were not incorporated as restrictions. LP algorithm was used to identify an upper bound for the scenarios based on the preset harvesting age; therefore, harvesting could occur only on stands 120 years or older. The relaxed, LP which aim the maximization of the volume harvested during the planning period, is max V୧,୨ ܺ, ୀଵ ே ୀଵ ݏݑܾ݆݁ܿݐ ݐ ∑ V୧,୨ ܺ, ே ୀଵ ሺ1 െ ߮ሻ ൈ ∑ V୧,୩ ܺ, ே ୀଵ ݂ݎ ܽ݊ݕ ݕ݁ܽݎݏ ݆ ܽ݊݀ ݇ ∑ V୧,୨ ܺ, ே ୀଵ ሺ1 ߮ሻ ൈ ∑ V୧,୩ ܺ, ே ୀଵ ݂ݎ ܽ݊ݕ ݕ݁ܽݎݏ ݆ ܽ݊݀ ݇ (3.8) ∑ ܺ, ୀଵ ܣ ݂ݎ ܽ݊ݕ ݑ݊݅ݐ ݅ ܺ, ൌ 0 ݂݅ ܽ݃݁ ݂ ݑ݊݅ݐ ݅ ൏ 120 ݕ݁ܽݎݏ ሺݐ݄݁ ݎ݁ݏ݁ݐ ݄ܽݎݒ݁ݏݐ ܽ݃݁ሻ where Xi,j is the area of unit i to be harvested in year j of the planning period p Vi,j is the volume per hectare of unit i in year j of the planning period p N is the number of units in the timber THLB ϕ is the maximum fluctuation of the annual harvested volume (in this case ϕ=0.1) Ai is the area of unit i. 35 The first two constraints ensure that the annual harvested volume does not vary with more than 10%, irrespective of the years during the planning period. The last two constraints ensure that a unit is not harvested more than once during the planning period and a unit cannot be harvested at ages less than the preset harvest age (i.e., 120 years). The LP (3.8) did not include constraints regarding the growing stock as for the formulation (3.8) does not lead to a significant reduction of the forest inventory. The non-significant difference between the final and the existing or optimal growing stock is the result of the formulation (3.8) and can be proven by comparing the AAC with the annual growth of the forest estate (e.g., AAC ≤ annual growth or AAC ≥ annual growth). From the enumeration of all possible comparisons, the cases fulfilling the LP constraints (i.e., even-flow is ensured between any two years of the planning period, the harvesting of a unit can occur at most once during the planning period and the length of the planning period is larger than half preset harvest age) have the growing stock non-significant different from the existing or optimal growing stock. The decision variables from (3.8) define the LP approach as a Model I formulation (Bettinger et al 2009). The two algorithms, adjusted FFD and SA, were implemented using Delphi 7.0 (Borland Software Corporation 2001) and all runs were executed on a 3 GHz Pentium 4 processor and Windows XP operating system. The same data manipulation procedure and programming platform was used for both algorithms to eliminate possible coding or compiling differences that could influence the comparison of the algorithms (Pukkala and Kurttila 2005). The LP solution was obtained using the C-Whiz 4.1 solver (Optimal Software 2001), while the matrix quantifying the LP problem was generated using Spectrum 3.0 (USDA Forest Service 2008). The results supplied by the three scheduling algorithms were analyzed with SAS 9.1 (SAS Institute Inc. 2004). 3.3 Results The introduction of a new forest planning method is usually tested for a single area (Mathey et al. 2007; Lockwood and Moore 1993; Clements et al. 1990). I enhanced this approach by using three areas to assess the performances of the adjusted FFD. The longitudinal analysis of the distributions describing the relationship between the volume and age and the land surface and age indicated that the differences are non-significant (p=0.84 and p=0.23, respectively), regardless the area. Having similar age – volume and age – surface distributions, the three areas were suitable for use in a comparison of the adjusted FFD algorithm with simulated 36 annealing. This similarity also offered the opportunity to evaluate the performance of the two methods in relation to the green-up / adjacency delay constraint, as consistent results on different areas fulfill the repeatability requirement, mandatory to any valid scientific investigation. For the selected merchantability standards and the corresponding OHA, all three areas revealed a left skewed distribution (Fig 3.3), indicating the fulfillment of the PBPT requirements. The volume available for harvesting, in the absence of spatial constraints, was larger than MVHA for an entire century (the case of Doig and Moberly areas) or for the next 40 years (the case of Fort Nelson area). The existing surplus of volume compared with MVHA indicates that the merchantability standards were not optimal selected, as other limits would have offered larger AAC. However, the MVHA was determined without spatial constraints; therefore, the availability of volume could partially or totally offset the impact of the green-up/adjacency constraint on the magnitude of AAC. Figure 3.3 Current volume and area distribution with respect to harvest timing. The decrease in volume and area with respect to harvest timing indicate the reduced number of units with low site index available for harvesting in the future. This situation is desired as units with longer OHA are difficult to harvest during a planning period shorter than OHA. The left-skewed distribution of the volume in relation to the harvest timing (Fig 3.3) suggested that the AAC determined considering spatial constraints could reach the MVHA. Based on the PBPT corollary and the proposed 7/9×MVHA threshold identified by Johnson (1973), I determined the AAC by performing several runs using the adjusted FFD. The runs successively 37 decreased the AAC threshold from MVHA to 0.75×MVHA in increments of -0.05×MVHA. This approach was possible since the computational time for each area was less than 3 min. In addition to the fast computation, the adjusted FFD algorithm is deterministic, and no additional runs were required. The incremental approach showed that AAC can exceed a threshold of 0.95×MVHA for all three areas when the adjacency/greenup delay is less than or equal to five years. An increase in the adjacency/greenup delay from 5 to 20 years decreased AAC supplied by adjusted FFD with 4% at the Doig area and with 11% at the Moberly and Fort Nelson areas, when harvest occurred at ages close to Agemax(MAI). Irrespective of the adjacency delay, the adjusted FFD supplied AACs significantly larger than the 7/9×MVHA limit (p<0.0001) when the harvesting occurred at ages equal or larger than OHA. The adjusted FFD algorithm supplied the greatest AACs within the four investigated scenarios for all areas (Fig 3.4 a, b and c) when the harvesting occurred at ages close to Agemax(MAI). The adjusted FFD supplied results constantly larger than SA, when the harvestings were scheduled according to the Corollary of the PBPT. However, the use of the adjusted FFD algorithm did not always ensure better results than SA. When the harvesting age was established above a preset age and the green-up/adjacency delay was less than 20 years, the adjusted FFD algorithm supplied the smallest AACs (Fig 3.4b). Similar to the adjusted FFD, SA produced the highest AAC when the stands were harvested at Agemax(MAI) and not at a preset value, except for the Moberly area (Fig 3.4b). Figure 3.4. AAC dependency on the algorithm and the green-up/ adjacency delay for Doig area (3.4a), Moberly area (3.4b) and Fort Nelson area (3.4c). The solid line represents the adjusted FFD algorithm (───) and the dotted line represents the SA algorithm (·····). The squares (■ ■ ■) indicates that the harvest occurred at Age max(MAI), while the disks are associated with harvestings at ages larger than 120 years (● ● ●) 38 The sensitivity of the AAC to the adjacency/greenup delay, identified at the adjusted FFD, was also present at SA, which dropped between 10% (Fort Nelson area) and 70% (Moberly area) when the adjacency/greenup delay increased from 1 to 20 years. While the adjusted FFD supplied a relatively small reduction in AAC (i.e., 20%) regardless of the harvesting age and adjacency delay, SA seems to be sensitive to both harvesting age (i.e., Agemax(MAI) or preset value) and adjacency delay. The AAC calculated using SA could decrease with adjacency delay more than two-thirds when the harvest occurred at max(MAI)) and between 20% and 40% when the harvest occurred at a preset age. The adjusted FFD yielded a counterintuitive result for the Moberly area as the AAC for an adjacency delay of 20 years was greater than that for 10 years delay. However, this anomaly is not especially rare for bin-packing problems, as indicated by Hoffman (1998, p. 172-173) and Rhee (1991). The variation of the AAC during the planning period plays a significant role in the forest decision process (Pedersen 2003), less variation of the AAC being preferred. The adjusted FFD algorithm supplied a relatively uniform AAC for all three areas, when the harvestings were scheduled according to the Corollary of the PBPT, namely Agemax(MAI). The AAC based on the Corollary of the PBPT reduced its variation with the adjacency delay, as the standard deviation of the AAC decreased from 1 to 20 years, in some cases by almost 50% (i.e., Moberly area). In comparison, SA produced increasingly variable results with the increment of the greenup/adjacency delay for all areas, augmenting the variability of AAC (expressed as the standard deviation) sometimes by more than 50% (i.e., Doig area). The LP supplied the maximal even flow AAC in the absence of spatial constraints and for harvests occurring at ages equal or larger than 120 years (Table 3.2). The LP solution was 16% smaller than the MVHA for all three areas. The existence of a ratio of 5/6 between the LP solution and MVHA for all the study areas suggests the existence of a relationship between the two algorithms that could be rooted in the deterministic character of both algorithms. A LP solution smaller than MVHA was not surprising as MVHA could schedule multiple harvests of the same stand during the planning period (sometimes even three times - the case of the lodgepole pine stands with site index larger than 30 and age larger or equal with OHA in year one of the planning period). The LP supplied AAC at most 25% larger than the corresponding AAC determined using SA (for Fort Nelson area) or 40% for AAC determined using adjusted FFD (for Moberly area). 39 Table 3.2. Comparison between the MVHA and the AACs supplied by LP, adjusted FFD and SA, for the adjacency delay of 1 year and minimum harvest age 120 years Study area MVHA [103 m3] AAC [103 m3] LP Adjusted FFD SA Doig 500 422 398 347 Fort Nelson 280 235 231 177 Moberly 315 263 156 212 3.4 Discussion The PBPT corollary established MVHA as the upper bound of the volume that can be harvested annually but did not indicate that the solution is unique. Coffman et al. (2002) showed that MVHA can be reached by scheduling stands for harvesting in more than one combination. This non-uniqueness feature of the adjusted FFD algorithm is desired, as it offers flexibility when implementing the strategic forest management plan. The comparison of adjusted FFD and SA with MVHA represents the highest level of validating a heuristic algorithm (i.e., level 6), while the comparison with LP supplies a level 5 validity, according to the classification of Bettinger et al (2009). Fig 3.4 (a, b and c) revealed that for preset harvesting ages, SA can supply results 50% smaller than MVHA, the aspatial optimal solution, (for the Fort Nelson or Doig areas with a 20 years green-up delay) but within 25% from the solution supplied by the LP (i.e., 80% for Moberly area, 75% for Fort Nelson area and 83% for Doig area). The relatively good performances of SA compared to LP agrees with the results of Boston and Bettinger (1999), Baskent and Jordan (2002), and Crowe and Nelson (2003), who, using SA and a preset harvesting age, achieved results within 95% the LP solution. The excellent performances of SA in forest scheduling mentioned in the literature could depend on the enforcement of a single harvest (at most once) during the planning period for each stand. By ensuring that a stand is not harvested more than once during the planning period, the length of the planning period was related with the scheduling algorithm and for northeastern British Columbia could lead to the harvesting of a stand almost 70 years before the OHA [e.g., a stand of white spruce with SI = 10 has OHA=190 years (Appendix C)]. The harvest prior to OHA could lead to a reduction of the volume planed for harvesting as much as 50% of the volume at culmination MAI (e.g., for white spruce with SI≤9 or subalpine fir with SI≤10). The selection of 40 the harvesting age had an evident impact on the AAC for the Moberly area, which was the only area with an average OHA greater than 120 years, the preset harvesting age (i.e., average OHA was 130 years and standard error of the mean 0.6 years). Considering the surplus of harvesting volume during the century planning period (Fig 3.3), scheduling the harvestings in Moberly area using SA before the optimality of the timber accumulation criterion supplied results larger than SA using OHA or the adjusted FFD but using the preset harvesting age and a greenup delay less than 20 years. Irrespective of the selection of the harvesting age, the largest AAC supplied by SA was at least 20% smaller than the MVHA (Fort Nelson and greenup of one year) confirming the findings of Wah and Wang (1999) and Baskent and Jordan (2002) and proved by Hoss and Stutzle (2005, p. 77-78) when compared with the known maximal solution. The adjusted FFD supplied results that were more uniform than those from the SA in relation to the greenup/adjacency delay, as the variance was smaller in 39 out of 48 comparisons (81%). A constant AAC is not only operationally and socially preferred but also allows integration into the forest planning process of other activities that tend to increase the variability of the results (such as financial fluctuations or land-use change) without exceeding the regulatory or predefined range of the variability of the AAC (in this case 20%). The benefit associated with the use of the simulated annealing algorithm is that it provides an acceptable solution in relatively short time, from seconds (Liu et al 2006; Ohman and Lamas, 2005) to minutes (Crowe and Nelson, 2001; Baskent and Jordan, 2002). However, its heuristic character recommends several runs, 20 in the case of Wah and Wang (1999) and Heinonen and Pukkala (2004) or 10 in the case of Baskent and Jordan (2002). In the analysis, for each area I performed 1820 runs of 15 min and 40 runs of one hour, and only when the stands were harvested at a preset harvesting age SA supplied larger AAC than the corresponding adjusted FFD but in only 4 out of the 12 scenarios (33%). In contrast, the adjusted FFD supplied the solution in 3 minutes and was larger than the SA solution when the harvestings occurred at max(MAI ) (Fig 3.4 a, b and c): from 9% for Fort Nelson area (i.e., adjacency delay 10 years) to 82% for Moberly area (i.e., adjacency delay 20 years). The adjusted FFD supplied AAC on average 25% greater and 100 times faster than those provided by SA. SA can produce larger AAC than the adjusted FFD if either the number of runs or the annealing time is converging to ∞. However, in practice, the number of runs is usually relatively small (less than 20), and in such cases Wah and Wang (1999) or Baskent and Jordan (2002) showed that SA could supply solutions that are 20% below the optimal solution, results confirmed by the findings. 41 Whichever study area, the greenup/adjacency constraint had a reduced influence on the adjusted FFD results, with less than a 15% decrease in AAC between 1 and 20 years delay (Fig 3. 4 a, b and c), except for one case (i.e., 20% for Doig area when the adjacency delay was 20 years and the harvests were scheduled using the preset harvesting age). In contrast, the greenup/adjacency delay had a major impact on the solution supplied by the SA as all areas showed a reduction in AAC of more than 15% (70% to Moberly area) when the delay increased from 1 to 20 years. From an operational perspective, the downward trend of the AAC when the adjacency delay is increasing might be offset with the use of fertilization, irrigation, or elite or genetically enhanced clones for stand regeneration. The intensive silviculture could have a double impact on the harvestings when SA is used for scheduling: increase the yield and reduce the greenup/adjacency delay. For adjusted FFD additional management prescriptions seeking mainly a decrease of the greenup/adjacency delay seems unjustified as the algorithm supplies results close to the optimal solution for delays less than 20 years. Nevertheless, the intensive silviculture could yield harvesting volumes superior to the MVHA, which was determined assuming that the yield of the present and the future stands are the same. 3.5 Conclusion To accommodate the optimal harvesting criterion at the forest and stand level I based the analysis on the PBPT corollary and demonstrated that the maximum volume that can be harvested annually should equal the sum of the maximum mean annual increment of the stands. To achieve MVHA, the harvested volumes have to be distributed uniformly or left- skewed. When the distributional requirements were not met, I developed a mixed integer programming solution that would ensure the fulfillment of the PBPT conditions. PBPT does not consider the constraints associated with the spatial arrangements of the stands, such as adjacency. To accommodate the set of spatial constraints associated with the greenup/adjacency delay the harvestings were scheduled using the unit restriction model (Murray 1999) within the First Fit Decreasing algorithm (Johnson 1973). The inclusion of the spatial constrains in the FFD algorithm using the unit restriction model, rather than developing a spatial specific algorithm inside FFD was preferred, as the FFD algorithm could self organize to structures resembling the PBPT requirements (Csirik et al. 1999), and supply an AAC close to the MVHA (Csirik and Johnson 2001). I compared the results supplied by the adjusted FFD with the results supplied by SA and a relaxed LP using three study areas in northeastern British 42 Columbia and the constraints enforced by a local Sustainable Forest Management Plan (Jukes et al. 2004). The adjusted FFD (i.e., the FFD algorithm constrained to meet the spatial restrictions and the conditions of the PBPT Corollary) supplied the greatest AAC from all areas, within 5% of the MVHA when the greenup/adjacency delay was less or equal than five years. The adjusted FFD and the relaxed LP produced results in less than three minutes; in comparison, SA required six hours. The AAC determined using the adjusted FFD was less variable than the AAC derived by SA. The adjusted FFD seems to be sensitive to the violation of PBPT distributional requirements, as for one of the areas (i.e., Moberly) it supplied both the smallest and the greatest AAC. The greenup/adjacency delay had little impact on the adjusted FFD algorithm but could reduce the AAC by 70% when SA was used. These results suggest that further research into the ability of the adjusted FFD algorithm to incorporate more complex spatial-temporal constraints, such as patch size distribution or financial marginal analysis, is warranted. A natural extension of this work would be the examination of the initial spatial distribution, which limits the accessibility of the stands for harvesting due to adjacency/greenup constraints. Finally, an examination of the reasons why simulated annealing is outperformed by the adjusted FFD algorithm could lead to amendments enhancing the efficiency of the heuristic techniques commonly used in forest planning. 43 3.6 References 1. Amacher, G. S., Ollikainen, M., and Koskela, E. 2009. Economics of forest resources. MIT Press, Cambridge MA USA. 397 p. 2. Anderson, D. R., Sweeney, D. J., and Williams, T. A. 1994 . An introduction to management science. Quantitative approaches to decision making. West Publishing Company, St. Paul, MN. 879 p. 3. Anonymous. 1997a. Fort Nelson Land and Resources Management Plan. Ministry of Sustainable Resources Management. 312 p. 4. Anonymous. 1997b. Fort St. John land and Resources Management Plan. Rep. 31090- 25-04. Ministry of Sustainable Resources Management. 229 p. 5. Anonymous. 1999. Dawson Creek Land and Resource Management Plan. Rep. 31090- 25-01. Ministry of Sustainable Resources Management. 232 p. 6. Baskent, E.Z. and Jordan, G.A. 2002. Forest landscape management modeling using simulated annealing. Forest Ecology and Management 165(1-3): 29-45. 7. Bettinger, P., Graetz, D., Boston, K., Sessions, J., and Chung, W.D. 2002. Eight heuristic planning techniques applied to three increasingly difficult wildlife planning problems. Silva Fennica 36(2): 561-584. 8. Bettinger, P., J. Sessions, and K. Boston. 2009. A review of the status and use of validation procedures for heuristics used in forest planning. Mathematical and Computational Forestry and Natural-Resource Sciences 1(1):26-37. 9. Bettinger, P., Boston, K., Siry, J., and Grebner, D. 2009. Forest Management and Planning. Academic Press, New York NY USA. 331p. 10. Borland Software Corporation. 2001. Delphi 7.0. Cupertino, CA USA. 11. Boston, K. and Bettinger, P. 1999. An analysis of Monte Carlo integer programming, simulated annealing, and tabu search heuristics for solving spatial harvest scheduling problems. Forest Science 45(2): 292-301. 12. Clements, S.E., Dallain, P.L., and Jamnick, M.S. 1990. An Operational, Spatially Constrained Harvest Scheduling Model. Canadian Journal of Forest Research-Revue Canadienne de Recherche Forestiere 20(9): 1438-1447. 13. Coffman, E.G., Courcoubetis, C., Garey, M.R., Johnson, D.S., Shor, P.W., Weber, R.R., and Yannakakis, M. 2000. Bin packing with discrete item sizes, Part I: perfect packing theorems and the average case behavior of optimal packing. SIAM Journal of Discrete Mathematics 13(3): 384-402. 44 14. Coffman, E.G., Courcoubetis, C., Garey, M.R., Johnson, D.S., Shor, P.W., Weber, R.R., and Yannakakis, M. 2002. Perfect packing theorems and the average-case behavior of optimal and online bin packing. SIAM Review 44(1): 95-108. 15. Coffman, E.G., Garey, M.R., and Johnson, D.S. 1997. Approximation Algorithms for Bin Packing: A Survey. P 46-93 in Approximation Algorithms for NP-Hard Problems, Hochbaum DS (ed). PWS Publishing, Boston. 16. Crowe, K. and Nelson, J. 2003. An indirect search algorithm for harvest-scheduling under adjacency constraints. Forest Science 49(1): 1-11. 17. Crowe, K., Nelson, J., and Boyland, M. 2003. Solving the area-restricted harvest- scheduling model using the branch and bound algorithm. Canadian Journal of Forest Research- Revue Canadienne de Recherche Forestiere 33(9): 1804-1814. 18. Csirik, J. and Johnson, D.S. 2001. Bounded space on-line bin packing: Best is better than first. Algorithmica 31(2): 115-138. 19. Csirik, J., Johnson, D.S., Kenyon, C., Shor, P.W., and Weber, R.R. 1999. A self organizing bin packing heuristic. P 246-265 in Lecture Notes in Computer Science 1619, Goodrich M, McGeoch CC (eds). Springer-Verlag, Berlin Germany. 20. Davis, L.S., Johnson, K.N., Howard, T.E., and Bettinger, P. 2001. Forest management. 4 ed. McGraw-Hill, New York USA 804 p. 21. Devroye, L. 1986. Non-Uniform Random Variate Generation. Springer-Verlag, New York NY USA 843 p. 22. Duerr, W.A. 1993. Introduction to forest resource economics. McGraw-Hill, Toronto Canada 485 p. 23. Garey, M.R., Graham, R.L., Johnson, D.S., and Yao, A.C.C. 1976. Resource Constrained Scheduling As Generalized Bin Packing. Journal of Combinatorial Theory Series A 21(3): 257-298. 24. Geman, S. and D. Geman. 1984. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence 6: 721-741. 25. Goycoolea, M., Murray, A.T., Barahona, F., Epstein, R. and Weintraub, A. 2005 Harvest scheduling subject to maximum area restrictions: Exploring exact approaches. Operations Research; 53(3): 490-500. 26. Gunn, E.A. and Richards, E.W. 2005 Solving the adjacency problem with stand-centred constraints. Canadian Journal of Forest Research-Revue Canadienne de Recherche Forestiere; 35(4): 832-842. 45 27. Hastings,W.K. 1970. Monte-Carlo Sampling Methods Using Markov Chains and Their Applications. Biometrika 57(1): 97-109 28. Heinonen, T. and Pukkala, T. 2004. A comparison of one- and two-compartment neighbourhoods in heuristic search with spatial forest management goals. Silva Fennica 38(3): 319-332. 29. Hoffman, P. 1998. The man who loved only numbers. Hyperion, New York USA. 302 p. 30. Hoss, H. and T. Stutzle. 2005. Stochastic local search. Morgan Kaufmann Publishers, New York USA. 658 p. 31. Johnson, D.S. 1973 Near-optimal bin packing algorithms. Ph.D. thesis, Massachusetts Institute of Technology, Dept. of Mathematics. Cambridge, Massachusetts USA 402 p. 32. Jukes,W., Menzies,D., Rosen,D., Taylor,G., Saint-Jean,R., and Beale,J. 2004. Sustainable Forest Management Plan. Fort St. John, British Columbia Canada 561 p. 33. Kangas, A.S. and Kangas, J. 1999. Optimization bias in forest management planning solutions due to errors in forest variables. Silva Fennica 33(4): 303-315. 34. Kangas, J., Kangas, A.S., Leskinen, P., and Pykäläinen, J. 2002. MCDM methods in strategic planning of forestry on state-owned lands in Finland: applications and experiences. Journal of Multi-Criteria Decision Analysis 10(5): 257-271. 35. Lockwood, C. and Moore, T. 1993. Harvest scheduling with spatial constraints: a simulated annealing approach. Canadian Journal of Forest Research 23(3): 468-478. 36. Mathey, A.H., Krcmar, E., Tait, D., Vertinsky, I., and Innes, J. 2007. Forest planning using co-evolutionary cellular automata. Forest Ecology and Management 239(1-3): 45-56. 37. Mcdill, M.E., Rebain, S.A. and Braze, J. 2002. Harvest scheduling with area-based adjacency constraints. Forest Science; 48(4): 631-642. 38. Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H., and Teller, E. 1953. Equations of State Calculations by Fast Computing Machines. Journal of Chemical Physics 21(6): 1087-1092. 39. Murray, A.T. 1999. Spatial restrictions in harvest scheduling. Forest Science 45(1):45- 52. 40. Murray,A.T. and Church,R.L. 1995. Heuristic solution approaches to operational forest planning problems. OR Spectrum 17(2-3): 193-203. 41. Nemhauser, G.L. and Wolsey, L.A. 1988. Integer and Combinatorial Optimization. John Wiley and Sons, New York NY USA 784 p. 42. Nourani, Y. and B. Andresen. 1998. A comparison of simulated annealing cooling strategies. Journal of Physics A-Mathematical and General 31(41): 8373-8385. 46 43. Nyland, R.D. 1996. Silviculture. Concepts and applications. McGraw-Hill, New York NY USA 633 p. 44. Optimal Software. 2001. C-Whiz. Sterling VA USA 45. Pedersen,L. 2003. AAC Rationale for the Fort St. John TSA. Victoria BC Canada 44p 46. Pukkala, T. and Kurttila, M. 2005. Examining the performance of six heuristic optimisation techniques in different forest planning problems. Silva Fennica 39(1): 67-80. 47. Rhee, W.T. 1991. Stochastic-Analysis of A Modified 1St Fit Decreasing Packing. Mathematics of Operations Research 16(1): 162-175. 48. Rucareanu, N. and Leahu, I. 1982. Amenajarea padurilor. Ceres, Bucharest Romania 438 p. 49. Shanks, M. E. and Gambill, R. 1973. Calculus. Holt, Rinehart and Winston, Inc., New York NY USA 748 p. 50. SAS Institute Inc. 2004 SAS Software 9.1.3. Cary NC USA. 51. USDA Forest Service. 2009. Spectrum. 3.0. Fort Collins CO USA. 52. Ullman, J.D. 1971. The performances of a memory allocation algorithm. Technical Report 100. Princeton NJ USA 53. Wah, B.W. and Wang, T. 1999. Simulated Annealing with Asymptotic Convergence for Nonlinear Constrained Global Optimization. P 461-475 in Principles and Practice of Constraint Programming, Jaffar J (ed) Springer Verlag New York NY USA. 54. Wolsey, L.A. 1998. Integer programming. John Wiley and Sons, New York NY USA 264 p 55. Zomaya, A.Y. and Kazman, R. 1999. Simulated Annealing Techniques. P 37-1-37-19 in Handbook on Algorithms and Theory of Computation, Atallah, M.J. (ed) CRC Press, Boca Raton FL USA. 47 4. AN ANALYTICAL PLATFORM FOR CUMULATIVE IMPACT ASSESSMENT IN NORTHEASTERN BRITISH COLUMBIA1 4.1 Introduction Human activities have continuously influenced the natural landscape in a multitude of ways - from relatively innocuous industries such as tourism and hunting, to more intrusive industries such as forestry and oil and gas extraction. The negative impacts of industrial development have become increasingly apparent as the environmental effects associated with these activities have expanded and challenged the capacity of ecosystems to recover from naturally-occurring perturbations (Fuhrer 2000; Marcu 1981). Over time, the accumulation of stresses resulting from human activities will change an ecosystem so profoundly that recovery is no longer possible using conventional processes (Gore et al. 1990; Toffler 1970). To prevent the degradation of the environment beyond socially accepted limits, the effects of economic developments are evaluated through environmental impact assessments (Koornneef et al. 2008). The main difficulty in assessing the environmental responses to human developments relate to the reduced impact of most projects when considered individually that may develop significant impacts when assessed in the context of past, present and foreseeable future activities. As a result, the combined influence on the environment of all projects occurring in a single area is evaluated through cumulative impact assessments (CIA). These consider the consequences of multiple projects, each insignificant on its own, yet important when considered collectively (Council on Environmental Quality 1969). The CIA literature addresses extensively the qualitative aspects of the environmental changes induced by the development of society (Contini and Servida 1992; Glasson et al. 1994; Marr 1997; Masera and Colombo 1992; Reid 2001). However, only a few studies have elaborated specific methods to be used in evaluating the cumulative impacts associated with human development (Spaling and Smit 1994). Currently the most popular CIA methods are based on complex modular models (Dube et al. 2006; Voinov et al. 2004) or operations research (Stakhiv 1988). The CIA literature abounds in detailed investigations based on methods initially developed to assess particular activities, such as land-use change associated with urban development (Dickert and Tuttle 1985; Stakhiv 1988), water resources development (Dee et al. 1973; Gosselink and Lee 1989) or construction (Leopold et al. 1971), all of which were 1 A version of this chapter has been submitted for publication. Strimbu, B.M. and Innes, J.L. An analytical platform for cumulative impact assessment in northeastern British Columbia. 48 enhanced for use as analytical support for CIA (Rotmans and van Asselt 2001). Regardless of the method adopted, current CIA modeling platforms are based on the assumption that small incremental changes are the driver of the environmental dynamics (Cherp et al. 2007; Dube et al. 2006; Hegmann et al. 1999). The complexity of CIA investigation combined with the intricate structure associated with modular modeling led to CIA based on results supplied by one analytical platform (Voinov et al. 2004). The mono-analytical approach for CIA could convey erroneous conclusions, as from a theoretical perspective the probability of occurrence of a complex model converges on zero. The asymptotic null probability of the results supplied by the complex methodologies employed by CIA is rooted on the model’s framework, which can be fully empirical (Patil et al. 2002), fully process-based (Leimbach and Jaeger 2005) or a combination of the two (Wu and David 2002), in the sense of Korzukin et al. (1996). CIA using models based exclusively on processes, theories and laws encountered difficulties in quantifying the relationships between environmental dynamics and socio-economic, biophysical and land-management variables (Verburg et al. 2004). Process-based CIA has no confidence intervals for the predicted values, casting doubts on the results. The complex models incorporating empirical components supply precise predictions but with reduced accuracy, as the probability of occurrence of the outcome of a model is: Probability of occurrence = 1/ (the number of equally likely results) ≤ (the smallest number of equally likely results supplied by each empirical equation)-n ≤ (4.1) ≤ n n − ∞→ + )11(lim = 0, where n the number of non-overlapping units used by the model and the two 1s correspond to the smallest (respectively largest) value of the confidence interval supplied by the empirical equations of each unit. To address the lack of confidence associated with such results; this study proposes a new CIA analytical framework based on a change of the paradigm governing CIA investigations. The focus addressed here is to reveal the existence or not of pattern(s) in a possible future rather than providing a detailed investigation of one improbable future. 49 4.2 Methods 4.2.1 Paradigm shift A paradigm shift is required because the mono-analytical platform of current CIA is unable to forecast accurately complex environmental systems. CIA methodologies constrained to meet all legislative and regulatory requirements epitomize the development of the environment but provide little insight into the integrated behavior of the environment (Rotmans and van Asselt 2001). For CIA that considers the future of the environment as an integrated and significant part of the assessment, the investigation should focus on identifying the commonalities or the differences throughout the planning period. The commonalities and differences would delineate the patterns within the environmental dynamics. To identify these patterns I focused on the analysis of a multitude of futures rather than one future, an approach initially developed by Boltzman (1995) for statistical thermodynamics. To ensure that the CIA is performed within the framework established by Boltzmann, the futures should be independent and with the same probability of occurrence. Therefore, the paradigm shift promotes the development of a set of futures, each independent and equally likely, which would be used to identify the main environmental attributes affected by human activities as well as the moments when significant environmental changes could occur. 4.2.2 Analytical platform The CIA method proposed here is spatially and temporally explicit and assumes that all activities occur within the regulatory framework. Additionally, because the paradigm shift requires the investigation of a set of futures, the method only considers feasible futures The restriction of human activities to those permitted within the regulatory context does not ensure that the evolution of the environment is on the desired path, only that some attributes used to represent the environment numerically do not exceed the appropriate thresholds. This is a substantial difference from current CIA methodologies that determine whether or not the regulatory thresholds affecting the valued ecosystem components are exceeded and then recommend measures to mitigate or eliminate the undesired effects of the human activities on the respective VECs [e.g., Osowski et al (2001), Dube et al. (2006), Perdicoulis and Piper (2008)). The proposed method has previously been suggested by several authors [such as Duinker and Greig (2007)], who argued for a departure from a specific future to a set of possible futures. 50 To illustrate the proposed method, I performed a CIA for three areas from northeastern British Columbia (Fig 4.1). The predominant forest species in the area are white spruce (Picea glauca), black spruce (Picea mariana), trembling aspen (Populus tremuloides), lodgepole pine (Pinus contorta), balsam poplar (Populus balsamifera), engelman spruce (Picea engelmannii) and subalpine fir (Abies lasiocarpa) (Meidinger and Pojar 1991). The biogeoclimatic (BEC) zones dominating the area are the Boreal White-Black Spruce zone, on the western side of the Western Canadian Sedimentary Basin, and the Spruce-Willow-Birch zone and Engelmann Spruce-Subalpine Fir zone that are found at elevations above 1300 m (Ministry of Sustainable Resource Management 2004). Sedimentary (fine and coarse clastic sedimentary rocks, siltstone, and mudstone) and metamorphic rocks (greenstone and green-schist) are prevalent in the area, leading to the development of regosols, gleysols, cryosols or humo-ferric podzols, depending on the local conditions (Meidinger and Pojar 1991). Figure 4.1 Study areas 51 The study areas were selected to reflect the degree of development of oil and gas development activities: low (Moberly area), medium (Fort Nelson area) and high (Doig area). In addition to oil and gas development, I considered forest harvesting, which is traditionally excluded from CIA investigations but could lead to significant and undesirable changes in the environment when combined with the forest clearance required by oil and gas activities. The assessment of the cumulative impact of the two activities was performed on two valued ecosystem components (VEC), in the sense of Beanlands and Duinker (1983), which were identified by the Treaty 8 Tribal Association as representative for the area: moose (Alces alces) and American marten (Martes americana). These two species appear to have no direct relationships (Franzmann and Schwartz 1997) and are therefore suitable subjects for an investigation of the synergistic properties of interest to CIA (Preston and Bedford 1988). I considered just two human activities and two VECs to ensure the generality of the method, as multi-dimensional investigations (i.e., n activities and p VECs) are similar with the two-dimensional ones (Gilbert 1976). Zhang and Montgomery (1994), Chaubet et al. (2005), Therivel and Ross and Kienzle (2004) identified a significant impact of scale on the results when raster data are used in the environmental investigation. To avoid adding scale as a variable in the CIA, I used scale-free vector data (Green et al. 2006). The vector data were supplied by the BC Ministry of Sustainable Resource Management (2004) for forest and wildlife, and by the Oil and Gas Commission (2005) and Energy Information Administration (2005) for petroleum wells. The computations were performed with SAS 9.1 (SAS Institute Inc. 2004) and ArcGIS 9.3 (Environmental Systems Research Institute 2008) was used for the initial data preparation. Based on Boltzman’s approach (1995), the paradigm shift requires the development of a set of futures, each future being considered a sample unit that would be used to delineate the significantly different periods that could exist along the environmental evolutionary tracks (Fig 4.2a). To ensure valid inferences, the set of futures has to cover a large range of possible developments of the environment, with a minimum 10 futures / area being recommended by Tran (1997). Because each future has to be independent and equally likely, I described each future using different spatially and temporally explicit forecasting techniques. The minimum number of futures can be obtained by developing multiple independent models to describe each human activity considered in the CIA, leading to the seemingly parallel structure from Fig 4.2b. 52 Figure 4.2 Analytical framework for the shift of the CIA paradigm (a) and the development of set of equally likely and independent futures (b). The structure apparent in Fig 4.2b differs from traditional scenario analyses [such as Schnorbus and Alila (2004), Nitschke and Innes (2008), Menard et al. (2002) or Caplat et al. (2008)] which vary some attributes within a single quantitative platform, rather than developing a set of different quantitative platforms (i.e., set of different models describing the environmental processes considered in the analysis). In the study, I focused on generating more than 10 futures/area using techniques and results appropriate for the representation of the two human activities being considered, which were heuristic techniques for forest harvesting (Bettinger et al. 2002) and autoregressive (AR) models for petroleum drilling (Walls 1992). As mineral rights supersede surface rights in most jurisdictions, I considered that petroleum drilling will not be impeded by any ownership rights or forest management activities, except when specifically regulated by law (such as archeological sites, sites associated with endangered species or waterfowl nesting sites). 4.2.3 Generating futures To generate the independent and equally likely futures I adopted a hierarchical approach that mimics the real development: the forest harvest accommodates the harvest associated with the petroleum drilling. The hierarchy of human activities is imposed by oil and gas legislation, which 53 attach greater importance to mineral rights than surface rights. For each future, the activities were modeled in a descending sequence according to the position in the hierarchy (i.e., petroleum drilling first, as the least constrained activity, followed by forest harvesting). CIA investigations commonly do not consider long-term forecasts (Lenzen et al. 2004; Nitschke 2008). However, investigations combining oil and gas activities with forest harvesting require a time period determined by the largest duration of the two activities, in this case forest harvesting (Schneider et al. 2003). To illustrate the proposed CIA method, the futures were 100 years long, as the length of the forest plan is usually 100 years, and the centennial CIA was recommended by Schneider et al. (2003). 4.2.3.1 Petroleum drilling Spatially and temporally explicit petroleum drilling futures were generated using a two step approach. The evolution of the number of wells was forecasted and then the predicted number of wells was distributed across the landscape. This approach was motivated by the large number of models developed to predict the numbers of wells (Iledare and Pulsipher 1997; Ringlund et al. 2008; Walls 1994) and the lack of investigations that predict directly the spatial- temporal dynamics of the well-drilling process (AXYS Environmental Consulting Ltd. 2003; Rao 2000). The temporal evolution of petroleum drilling was represented using either the number of wells drilled in each time period (Cuddington and Moss 2001; Iledare 1995; Ringlund et al. 2008) or the total number of wells existing in the landscape (Walls 1994; Walls 1992). Econometric approaches commonly use short time steps in modeling petroleum drilling activity [such as three months (in the case of Canada and Latin America in Ringlund et al. (2008), months (Sadorsky 2001) or weeks (Lanza et al. 2005; Manning 1991)], but when the discovery processes are added to the drilling model then the time step is increased to a year (Cuddington and Moss 2001; Walls 1994) or two years [for the Middle East, Africa or Asia-Pacific in Ringlund et al. (2008) ]. The proposed CIA platform includes forest harvesting (usually modeled using yearly time steps); therefore the futures were produced using a 12 month time-step. The three areas selected to test the proposed CIA are located in the Western Canadian Sedimentary Basin and reflect the range of oil and gas activities occurring in the region. The equations describing the temporal evolution of the petroleum well count were therefore not developed for each area but for the entire region. Walls (1994), Lanza et al (2005), Rao (2000) and Holland (2008) showed that an autoregressive equation can be used to describe the well count evolution. I enhanced the autoregressive approach of Ringlund et al. (2008) with the data 54 investigation techniques proposed by Balke and Gordon (1986) and added a set of covariates to the auto-regressed variable: ttttt sXLCyLB νμ +++= )()( (4.2) where yt and Xt denote the drilling activity and the matrix of the covariate variables in period t, B(L) and C(L) are the two lag-polynomials in the lag operator L, given by: ∑ ∑ = = = −= q i i i ii u i i i LcLC tscoefficiencandbwithLbLB 0 1 )( 1)( μt is the trend, st is the seasonality and νt is the auto-correlated error term: ∑ = − +−= m i titit 1 ενϕν with tNIIDt ∀),0(~ 2σε Following the models proposed by Ringlund et al (2008), the set of covariates used to model the number of wells included the oil prices at the New York Mercantile Exchange. I enhanced the approach of Ringlund et al (2008) and Lanza et al (2005), which used only one oil price as a covariate (i.e., monthly or trimestrial and annual, respectively) by considering the minimum and the maximum oil price in addition to the average. The autocorrelation, inverse autocorrelation and partial autocorrelation functions for the annual and total number of wells were used to identify the order of the autoregressive model (Cleveland 1972). The method proposed by Cleveland (1972) was also used to identify the presence of a non-linear trend in the model describing the dynamics of the number of wells. The augmented Dickey-Fuller (DF) test (Dickey and Fuller 1979) evaluated whether or not the model describing the well drilling dynamics had an intercept and a linear time trend term. The DF test was supplemented by the Bartlet’s Kolmogorov – Smirnov test which assessed whether or not the difference between two consecutive years was white noise, in the sense of Brockwell and Davis (1996) To enhance the AR model mentioned by Walls (1994), Lanza et al. (2005) and Holland (2008), I used the Smallest Canonical (SCAN) correlation and Extended Sample Autocorrelation Function (ESACF) method to tentatively identify whether or not an ARMA process (Brockwell 55 and Davis 1996) could represent petroleum drilling, as proposed by Tsay and Tiao (1985), Choi (1992) and Box et al (1994). In the case that the SCAN and ESCAF indicated the presence of an ARMA process, the Minimum Information Criterion method was used to identify the order of the possibly stationary and invertible ARMA process. The Bayesian information criterion was employed to determine the autoregressive order required to estimate the error series (Hannan and Rissanen 1982). To model counting processes (e.g., the number of petroleum wells) Lindsey (1995) and Kingman (1995) advocated the use of generalized models as an alternative to time series models. Two models can represent the drilling process: Poisson and negative binomial (Eq. 4.3), as presented by Cameron and Trivedi (1998). The same equation (Eq.4.3) described the drilling activities regardless of the model type: ttt XLCy ε+= ))(exp( (4.3) The negative binomial model differs from the Poisson model by allowing the mean of the process to differ from its variance (a property of a Poisson process) by a variable that is gamma distributed. Indications regarding the possible use of Poisson and negative binomial regression models can be supplied by the DF test, autocorrelation and chi-square test (Rees 2002). The models developed for the entire northeastern British Columbia were implemented for each of the three study areas, after the significance of each equation was tested on each area. Regardless of the model type (i.e., autoregressive or generalized model), the representation of the process describing the number of wells was considered complete when the residuals were white noise, assessed using Bartlet’s Kolmogorov-Smirnov test, and stationary, assessed with the Philips-Peron test (1988). To avoid numerical errors and possible overestimation of the parameters of the variance the conditional index and the variance inflation factor (i.e., the measures of numerical errors and parameters overestimation, respectively) should be less than 10, as recommended by Belsley et al. (1980), and Neter et al. (1996). Within a time series framework, the violation of the assumptions can lead to biased or erroneous results (Phillips 1986) and I therefore checked the normality of the residuals using the Jarque and Bera test, homoskedasticity using White’s test, and serial correlation using Vinod’s version of the Durbin- Watson (DW) autocorrelation test (Vinod 1973). I supplemented the DW autocorrelation test with the McLeod and Li portmanteau test. 56 Once the annual amounts of drilling had been determined, indicator kriging (IK) was used to spatially distribute across the landscape the number of wells predicted to be drilled in each year. Indicator kriging was developed by Journel (1983) to predict variables that operate only on the binary set {0,1}, which in the petroleum drilling case corresponds to presence or absence at a given location of a well. However, kriging supplies values into a set organized as a field (Hungerford 1990), therefore the values supplied by IK could be considered as quantifying the probability of a well being drilled at a specific location. The indicator variogram [Eq. 4.4] used by IK is an adjustment of the classical variogram to a variable that has a Bernoulli distribution. })];();({[5.0);( 2zxIzhxIEzhI −+×=γ (4.4) where h is the distance between wells, x are the coordinates and I(x,z) is the indicator function at location x, and z is a real number quantifying the presence or absence of a well at location x. In addition to the indicator variogram, the robust variogram proposed by Hawkins and Cressie (1984a) was adjusted for the indicator variable (Eq. 4.5) and was used in computations: ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ −+= ∑ 5.0 )( 1 |)()(| )( 1 )(/988.0914.0 1)( j hN i xIxIhNhN hγ (4.5) where N(h) is the number of points quantifying the presence or absence of a well for the h distance class. I tested the directional symmetry of the spatial correlation using the isotropy test statistics proposed by Lu and Zimmerman (1994). The IK requires three attributes in the input dataset: the vector x (which gives the position) and 0 and 1 [that identify the presence or absence of a well at location x]. The data supplied by the Oil and Gas Commission (2005) can be used to assign only the value for 1; therefore a method to identify the 0 within the landscape had to be developed. An intuitive procedure that assigns a value of 0 to areas where no drilling can occur (e.g., inside populated areas, archeological sites, streams or lakes) failed to identify the 0 in areas where drilling can occur (e.g., pastures or forest with no harvesting constraints). An intuitive approach is therefore unable to delineate 0 from 1, and could lead to clusters of 0s, that would associate large lag distances to non- negative values, of little use in the variogram development. Consequently, a method to identify the locations where drilling is allowed but it is not performed had to be developed. A method used to determine the places where no drilling would occur should consider two observations: 57 • for any given year the position of wells that would be drilled in subsequent years was unknown; • the true variogram would be known when no further drilling occurs, as the maximum number of wells in the area is reached. When drilling will cease, the variograms built using data before the maximum number of wells is reached converges to the true variogram (Kolmogorov and Fomin 1999). The above two observations indicated that a set of variograms should be computed by splitting the dataset containing the wells information in two subsets, according to the year [i.e., the year would separate the set of wells in wells with known position (past) and wells with unknown position (future)]. The change of the separation year during the petroleum exploration period would create a set of variograms that would converge to the true variogram. Therefore, for the wells drilled before a certain year the I(x,z)=1 while for wells drilled after the respective year the I(x,z)=0. Based on the two observations and using data supplied by the MSRM, I developed variograms for each of the three areas, as each area has different geomorphologic features that could lead to different empirical variograms. The empirical variograms (i.e., indicator and robust) were used to identify the theoretical variograms representing the spatial distribution of the wells. The variogram fitting was assessed using Q1 and Q2 statistics, as suggested by Kitanidis (1991). The IK computes the probability of petroleum drilling at specific locations, commonly the coordinates of the vertexes within a grid. To determine the grid required by IK, each quarter- section (i.e., quarter of a square mile) was separated in four quadrants, each quadrant being associated with one of the four wells that can be drilled on a quarter-section without the consent of the Agricultural Land Commission of British Columbia (Agricultural Land Commission and Oil and Gas Commission 2002). Within each quarter section, the Cartesian coordinates of the wells were computed. Using the Cartesian coordinates, all the wells were plotted at the quarter-section level. The grid size was determined as the minimum distance between any two wells that is a divisor of a mile. The integration into a mile of the IK grid was advocated to conform to the grid used by the Government of British Columbia. The probability of occurrence of a new well within a quadrant was considered as being proportional to the number of wells already within the quadrant. Of the four values, only the largest two were selected, to ensure the compliance with the regulations constraining drilling 58 activities (Government of British Columbia 2004). The set of locations where drilling could occur was determined using data from the Oil and Gas Commission of British Columbia, the quarter- section grid used by the Government of British Columbia, and the constraints imposed by the Petroleum and Natural Gas Act and Land Act (Government of British Columbia 2004). Based on petroleum developments in Alberta (Information Services AEUB 2005), the Canadian province with the greatest amount of petroleum exploration, the density of oil and gas wells on the landscape will reach the maximum capacity allowed by the regulations (Anonymous 2004). Therefore, the grid for which the IK would be computed contains the center of the two quadrants having the largest number of wells within the quarter-section grid. Indicator kriging supplied the probability of drilling at a certain location, and the drilling moment was determined assuming that larger probabilities are associated with earlier drilling; an assumption based on kriging’s exact interpolation property (Cressie 1993) and strict monotonicity of the cumulative probability function (Grimmett and Stirzaker 2002). The probability of drilling at a specific location, conditioned by the number of wells drilled annually, provided the spatial and temporal dynamics of the wells. 4.2.3.2 Forest harvesting To meet the regulations governing the coordination of activities occurring on the landscape, the CIA method assumed that any forest harvesting would include the areas harvested for petroleum drilling (i.e., well pads). To accommodate the hierarchy of different industries changing the landscape, the current forest schedulers used in northeastern British Columbia to determine the annual allowable cut (AAC) are spatially and temporally explicit heuristic algorithms. I used two heuristic algorithms: simulated annealing (SA) (Metropolis et al. 1953) and an adjusted first fit decreasing (FFD) algorithm (Johnson 1973). For forest scheduling problems, SA (Lockwood and Moore 1993) has been identified as one of the heuristic algorithms with the simplest implementation but with results close to the linear programming solution attached to the forest planning problem, which is considered to be the optimal solution (Crowe and Nelson 2003; Bettinger et al. 2002). The SA is sensitive to the initial selection of the parameters used in computations [i.e., initial temperature (Ti), freezing temperature (Tf) and the rate of annealing (ra)]. The identification of the optimal starting parameters is commonly obtained by experimentation, which is time consuming. To attain solutions in times relatively comparable with adjusted FFD the parameters were selected following the recommendations of Boston and Bettinger (1999), Bettinger et al. (2002) and Crowe and Nelson (2003) and not by 59 experimentation. To obtain a solution I allowed the SA to run on average 60 min, in agreement with the approach of Crowe and Nelson (2003), which considered that SA supplied an efficient solution to the forest planning problems in less than 40 min. The selected rates of annealing were 0.85, 0.9, 0.95 and 0.99, the initial temperatures were 1000, 10000 and 100000, and the freezing temperatures were 1, 50 and 100, in agreement with the results of Bettinger et al (2002) and Crowe and Nelson (2003). The stochastic character of the SA required a series of runs to be performed; out of each the best solution is selected. To ensure that the solution supplied by SA can be compared with the adjusted FFD from the computational time perspective, five runs were executed using the same set of constraints. The adjusted FFD algorithm amends the FFD algorithm (Johnson 1973) to the spatial constraints starting from a determined and non-spatial maximal AAC. Using the perfect bin- packing theorem (Coffman et al. 2002), Strimbu et al. (2008) proved that the maximum volume that can be harvested annually in the absence of spatial constraints but fulfilling regulatory requirements equals the sum of the maximum of the MAI of all the stands. The inclusion of the spatial constraints in the scheduling process could lead to values smaller than the upper AAC boundary, so an AAC that met the regulatory requirements was determined by decreasing the upper AAC boundary in steps of 5% until all the forest planning constraints were fulfilled. The set of constraints used to compute the AAC were based on the Sustainable Forest Management Plan developed for the Fort St. John Pilot Project (Jukes et al. 2004a). Among the large number of constraints used to determine AAC (Jukes et al. 2004b), the main restrictions governing the magnitude of the volume that would be harvested annually and included in the computations were: An annual variation in AAC less than 10%, regardless the pair of years during the planning period; The only silvicultural treatment is clear-cutting; The opening size lies between 5 ha and 60 ha; Greenup/adjacency delay of 5 years; A planning horizon of 100 years. In addition to the two scheduling algorithms, two harvesting ages were used to determine the AAC: a preset minimum harvesting age (Crowe and Nelson 2003) and age at maximum mean annual increment (Davis et al. 2001). AAC maximization was ensured when the maximum volume to be harvested for the entire forest estate coincided with the maximum for each stand 60 from the forest estate. When volume is the attribute of interest in the planning process, the maximum volume to be harvested at the stand level occurs at the age when MAI peaks (Davis et al. 2001). The culmination age of the MAI was determined considering the merchantability standards of 17.5 cm for spruce and sub-alpine fir and 12.5 cm for all other species (Pedersen 2003). Stand harvesting at ages larger than a preset value has been suggested by Bettinger et al. (2002) and Crowe and Nelson (2003), and for the three study areas was approximately 120 years (Jukes et al. 2004b). However, a preset harvesting age does not ensure AAC maximization, and could lead to harvesting being undertaken 70 years before the culmination of MAI (e.g., a stand of white spruce with SI ≤10). Nevertheless, I developed future forested landscapes using a preset harvesting age as being the current method used to schedule the forest harvestings in the area. In all, four futures were developed: two algorithms x two harvesting ages (i.e., FFD/age@max(MAI), FFD/Preset age, SA/age@max(MAI) and SA/Preset age). The forest scheduling algorithms were implemented using Delphi 7.0 (Borland Software Corporation 2001) and all runs were executed on a 3 GHz Pentium 4 processor and Windows XP operating system. The stand growth and yield as well as the mean annual increment was computed using TIPSY 3.1 (Ministry of Forests 2005a) for managed stands and WinVDYP (Ministry of Forests 2005b) for unmanaged stands (a total of 18000 stands and more than 100 harvest ages). 4.2.3.3 Moose and American marten habitat The two VECs considered in the proposed CIA (i.e., moose and American marten) were assessed using habitat suitability index (HSI) models. HSI is the most common method used in CIA studies to quantify the quality and quantity of the habitat of the species (Graham et al. 2000; MacDonald 2000; Nitschke 2008; Osowski et al. 2001). While HSI models do not reflect accurately the response of the species to landscape changes, they do provide indications of habitat quality and availability, a necessary but not sufficient condition. However, a large number of studies have identified a significant correlation between the population dynamics and different indexes describing the habitat (Hirzel et al. 2006; Larson et al. 2004; Ottaviani et al. 2004), including HSI. The HSI used in the study area are restricted to assess the habitat suitability during winter, the limiting season (Allen et al. 1988; Dussault et al. 2006; Romito et al. 1999; Takats et al. 1999). The HSI equations quantifying the habitat suitability of each stand are: 61 For moose: ⎩⎨ ⎧ ≥ <= %251 %2525/&% Saplingsfor SaplingsforSaplingsDeciduousShrubs HSI M (4.6) And for marten: 3214 SSSSHSI AM ×××= (4.7) where ⎪⎪⎩ ⎪⎪⎨ ⎧ ≥×− <≤ <≤−× < = %700023.161.1 %70%301 %30%52.004.0 %50 1 ClosureCanopyClosureCanopy ClosureCanopy ClosureCanopyClosureCanopy ClosureCanopy S ⎩⎨ ⎧ ≥ <×+= %50&1 %50&&016.02.0 2 CanopyinFirSpruce CanopyinFirSpruceCanopyinFirSpruce S ⎪⎩ ⎪⎨ ⎧ ≥ ≤≤−× > = mheightTree mheightTreeheightTree mheightTree S 151 1555.01.0 50 3 ⎪⎩ ⎪⎨ ⎧ ≥ ≤≤−× > = %15&,1 %15&,%55.0&,1.0 %5&,0 4 FirSprucePine FirSprucePineFirSprucePine FirSprucePine S The two HSI models are based on the rather questionable assumption that proximity to human settlements or human activities do not affect the species behavior (Romito et al. 1999; Takats et al. 1999). Additionally, the two models hypothesize that water and mineral resources can be obtained in the area that supplies food and cover and are not limiting during the winter season. The forest attributes from Eq. 4.6 and Eq. 4.7 quantify the changes in habitat suitability induced by both petroleum drilling and forest harvesting and confirmed the absence of a direct dependency between moose and American marten. To assess the cumulative impact of the two activities I computed two statistics for each year during the 100 year planning period, one describing the average change of the HSI at the landscape level and one revealing the availability of superior habitats (i.e., HSI>0.5) across the landscape. I selected a threshold of 0.5 for the HSI following the results of Jager et al. (2006), Kroll and Haufler (2006), Liu et al.(2005) and Bailey et al. (2002). The average HSI is an indication of the habitat dynamics at the landscape level and signals whether or not human activities have significantly changed the landscape (Nitschke 2008). The superior habitats index was used to complement the average HSI, as non-realistic steep declines in species occurrence could be modeled when poor habitats are not excluded from the landscape (Jager et al. 2006). 62 4.2.4 Cumulative impact assessment The framework used to develop the futures (Fig 4.2a and 4.2b) ensured that the regulatory conditions governing human activities were not violated, as each future was hierarchically constrained to meet the applicable environmental, social and economic requirements. This hierarchical approach is an enhancement of the CIA modular modeling proposed by Stakhiv (1988), which considered that all activities have the same priority and optimized the economic development subject to environmental constraints. It also improves on the studies of Voinov et al (2004) and Lam et al (2004), neither of which considered human activities within a hierarchical framework, and the study of Costanza et al (2002), which considered unrealistic human activities (e.g., complete conversion from forest to residential or the opposite). However, the hierarchy does not ensure that the activities situated on the lowest hierarchical level will not supply futures that could have the attributes that quantify the VECs violating the regulatory thresholds. To ensure that all the futures fulfilled the regulatory requirements, I compared the HSI of each future and year during the planning period with a threshold HSI. To identify the HSI thresholds, I considered that the minimum amount of suitable habitat should not interfere with the population dynamics of moose or American marten. From this perspective, both species could present a change in the number of individuals larger than 50% (i.e., increase or decrease), as Peterson (1997) found for moose and Fryxell et al (1999) revealed for American marten. Using simulated data Jager et al. (2006) showed that for American marten a 50% reduction in population size could be induced by a 17% loss in the habitat, a value in agreement with Chaplin et al. (1998) who live-trapped the American marten and found that the population did not change significantly when the surface of the habitat diminished by 20%. While a decrease in the size of the habitat by a fifth could lead to significant changes in the population of American marten, a different perspective was identified for moose, as a reduction with 40% of the habitat due to forest and hunting activities was not reflected in a significant reduction of the moose population (Rempel et al, 1997). In fact, Rempel et al (1997) argue that in the absence of hunting, a 44% habitat disturbance induced by the forest harvesting did not lead to a decrease in the number of individuals. To be consistent with the results of the above studies, I considered that the regulatory thresholds could be met when the two statistics quantifying the habitat (i.e., average HSI and Area HSI>0.5) do not decrease by more than 15% for American marten and 40% for moose from the present values. The present values were selected as the reference values as moose and American marten are considered by the BC Conservation Data Center to be 63 ‘species apparently secure and not at risk of extinction’, and have G5 global conservation status (i.e., widespread, abundant and secure) (2008); a status that should be maintained. The scheduling of the forest harvesting was constrained to include forest removal associated with well pads, with the implication that landscape connectivity would be maintained [one of the requirements of the sustainable forest management plan in the region (Jukes et al. 2004b)]. The lack of fragmentation due to forest harvesting and well drilling is the result of the four-color theorem (Appel and Haken 1976), which is enforced by the green-up constraint. The four-color theorem ensures the existence within the forest of multiple corridors with ages at least 25 (the average age of each cluster of adjacent stands is larger than 8 years), as on average each stand has four neighbors stands. Considering Eq. 4.6 and 4.7 in the context of the four–color theorem, the dynamics of the HSI for moose and American marten would not be caused by the spatial arrangement of the stands. To identify whether or not the HSI is impacted by human developments as well as the moments when significant environmental changes could occur I performed a univariate longitudinal analysis of the set of futures using the Helmert transformation (Crowder and Hand 1990). In longitudinal analyses the validity of the F-test used by univariate tests depends on the fulfillment of the Huynh-Feldt condition by the covariance matrix (Huynh and Feldt 1970). The Huynh-Feldt condition was assessed using the sphericity test proposed by Anderson (1984). In cases where the covariance matrix did not meet the Huynh-Feldt condition, the significance level of the F-test was adjusted using the modifications proposed by Greenhouse and Geisser (1959) and Huynh- Feldt (1970). The existence of moments when the significant changes could occur depends on the models used to produce the set of futures. I therefore included an assessment of the differences between futures, and used the futures to identify periods with similar environments. I used a polynomial transformation to separate the futures, as recommended by Crowder and Hand (1990) for quantitative longitudinal analyses. The differences between futures as well as the identification of the periods between the moments when the environment could change significantly were established using the Scheffe test, as advocated by Zar (1996). The two measures used in the CIA investigation operate on two scales (i.e., the average HSI is between 0 and 1 and the area with HSI>0.5 is larger than 150,000). Rencher (2002) and Neter et al. (1996) mention that longitudinal analyses could be impacted by the difference in the 64 orders of magnitude. In order to avoid problems created by including in the same analysis values larger than 105 and values less than 100, the attributes describing the size of the area with HSI>0.5 were transformed to values between 0 and 1 using the linear function: )]5.0min()5.0/[max()]5.0min(5.0[ >−>>−> HSIAreaHSIAreaHSIAreaHSIArea (4.8) where min and max are the minimum and maximum operators applied to the HSI values of each future. 4.3 Results 4.3.1 The development of the set of futures For the number of wells drilled annually the SCAN correlation and ESACF method supplied completely different results, indicating the likelihood that an ARMA process cannot be used to model annual variations in well drilling. However, SCAN and ESCAF revealed that an ARMA (2, 2) process could represent the series describing the total number of wells. Nevertheless, minimum information criterion did not confirm the SCAN and ESCAF results as it associated one of the largest Bayesian information criterion to the ARMA (2, 2) process. Consequently, an ARMA process was not used to model the evolution of the number of wells. Of the three autocorrelation functions used in the analysis (i.e., autocorrelation, inverse autocorrelation and partial autocorrelation) only the autocorrelation for lag one year appeared to be significant (i.e. p<0.05), suggesting a first order autoregressive model for the series describing the wells drilled annually. The autocorrelation functions for the total number of wells indicated the presence of a non-linear trend and the absence of a relationship for lags greater than one year. Therefore, regardless the series used to describe the evolution of the number of wells (i.e., annual or total), there was no evidence to suggest that there was a relationship between drilling activities more than two years apart. The augmented Dickey-Fuller (DF) test revealed that the series describing the annual number of wells may be described by a difference – stationary process as p =0.999. Bartlet’s Kolmogorov Smirnov test confirmed the result of the DF test, and indicated that the difference between the numbers of wells in two consecutive years is white noise (p=0.74). However, the first model 65 from Table 4.1 (i.e., yt = 1.255 × yt-1 + εt) met all the time-series assumptions and indicated that within a linear regression context the annual number of wells can be predicted from the number of wells drilled in the previous year. A similar situation was encountered for the total number of wells, but in addition to the number of wells present on the landscape in the previous years, the models also included the drilling effort of the previous year (second model from Table 4.1). Table 4.1 Models describing the evolution of the number of petroleum wells in the northeastern British Columbia Model Residuals assessment (p value for … ) Equation Fit assessment Nor. Heter sked. Auto-correlation White noise Stationarity Pr>F Adj R2 Max VIF Max Cond Index DW # lags auto. >2σ ML Zero mean Trend 1 yt = 1.255 × yt-1 + εt 0.001 0.8 7 1 1 0.09 0.12 0.38 0 0.29 0.54 0.69 0.99 2 Totalt=1.01 × totalt-1 +1.1×yt-1+εt 0.001 0.9 9 7.27 5.2 0.05 0.18 0.34 0 0.49 0.98 0.78 0.17 3 y t =0.54× totalt-1 – 1.02×(year- 1948)2 – 0.44×year ×f(t)+ εt 0.001 0.9 0 7.7 8.8 0.13 0.73 0.34 0 0.12 0.33 0.12 0.51 4 yt=exp(3.62 + 5.7 x 10-4 x totalt-1 – 1.6 x 10-8 x total2t-1)+εt 0.47* Deviance = 54.26 0.01 0.93 0.20 0 0.32 0.83 0.28 0.80 5 yt=exp(3.69 +5.3 x 10-4 x totalt-1 – 1.39×10-8 ×total2t-1)+εt 0.75* Deviance = 1.14 0.01 0.01 0.15 0 0.05 0.81 0.82 0.98 Total t – total number of wells drilled at year t y t – number of wells drilled in year t f(t) – is the function sin(π(year-1948)/period), w period=12 before 1988 and 8 after 1987. DW – Durbin Watson autocorrelation test (Vinod 1973) ML – McLeod and Li portmanteau test (McLeod and Li 1983) * - for Poisson and negative binomial models the p-value refers to chi-square goodness of fit for Poisson, respectively negative binomial distribution εt ~ WN (0, σε). 66 Irrespective of the dependent variable (i.e., annual or total number of wells), the residual analysis revealed the presence of a parabolic or exponential function. In addition to the nonlinearity, one of the models for the number of wells drilled annually contains a periodic function with increasing amplitude and decreasing period (third model from Table 4.1), confirming the inappropriateness of an ARMA process (i.e., ARMA process is based on a linear function) to the northeastern British Columbia wells data. The first three models from Table 4.1 are regression based models and were developed using Eq. 4.2. The Poisson model (i.e., the fourth model in Table 4.1) has a deviance larger than 1, indicating possible over-dispersion. Additionally the likelihood ratio test proposed by Cameron and Trivedi (1998) revealed that the mean and variances are not equal (p<0.0001), casting doubt on the appropriateness of a Poisson model in this context and recommended the usage of the negative binomial model (the fifth model from Table 4.1). Models with deviances larger than 10 can be used if the deviance of the residuals do not reveal recording errors or significant outliers impacting the model (Neter et al. 1996). The investigation of the results did not exposed any recording errors or the presence of significant outliers, indicating that negative binomial (i.e., model five in Table 4.1) and Poisson models should be used for the interpretation. The analysis confirms the results of Ringlund et al (2008) or Lanza et al (2005), and reveals that regardless of the model type (i.e., regression or generalized model), annual oil prices do not play a significant role in modeling the drilling effort. All the models meet the time-series assumptions (Table 4.1), as the variance inflation factor was smaller than 8 [a threshold of 10 was recommended by (Neter et al. 1996)], and the conditional index was smaller than 9 [the maximum value recommended by Belsey et al (1980) was 10]. The residuals did not exhibit any significant departure from normality (p>0.05), were not heteroskedastic (p>0.12) except in the case of the negative binomial model, and were not autocorrelated (p>0.15 for DW test and p>0.05 for Mcleod and Li portmanteau test). Furthermore, all the models had the residuals organized as white noise (p>0.33) and were stationary with the mean and trend not significantly different from 0 (p>0.12 and p>0.17, respectively). The fulfillment of the assumptions and the agreement with the results of other studies investigating petroleum drilling (Iledare and Pulsipher 1997; Lanza et al. 2005; Rao 2000; Ringlund et al. 2008; Walls 1994) suggests that the models can be considered as correct and as having the same likelihoods. The futures built using the five models can therefore occur and are equally likely. All five models indicated that for northeastern British Columbia the maximum number of permitted wells, as regulated by current legislation (Agricultural Land Commission and Oil and 67 Gas Commission 2002), would be reached. The main difference in the results supplied by the five models lay in their asymptotic behavior, which is finite for the generalized models and infinite for the regression type models. For each study area, the set of five equations were refitted and compared with the equations determined for the entire northeastern British Columbia area. Most pair-wise parameters did not differ significantly (α=0.05). However, for the Moberly and Doig areas, the model including the periodic term (i.e., model 3 in Table 4.1) and the Poisson model was non-significant. In addition, for all three areas the negative binomial model had the coefficient of the squared total number of wells non-significant different than 0. However, the sign and magnitude of the non-significant variables corresponded to the ones developed for the entire region. Therefore, the evolution of the well drilling effort in the three study areas was modeled using the equations developed for the entire northeastern British Columbia (Fig 4.3), as the non-significant models or variables in an area enhance the landscape developments predicted by the significant models. The results from each study area confirmed the results obtained for the entire northeastern British Columbia region: the time series models will reach the maximum number of wells allowed by the regulations (i.e., 20,000 in the Doig area, 19,000 in the Fort Nelson area and 8,400 in the Moberly area) while the generalized models will not reach that maximum during the century long planning period (Fig. 4.3). Regardless of the model used (i.e. Poisson or negative binomial), the three study areas had different peak moments: in approximately 2025 for the Doig area, 2045 for the Fort Nelson area and after 2050 for the Moberly area. Furthermore, for the Moberly area, the generalized model had a non-asymptotic behavior similar to the autoregressive models, as the peak moment would not be reached before 2050. Figure 4.3 Evolution of the number of wells drilled annually 68 The IK method, used to distribute the temporal evolution of the number of wells across the landscape, did not consider geology as influencing the location of a well. For northeastern British Columbia 62% of the wells were drilled above known petroleum reservoirs. This departure from pure randomness (50%) and large sample size (almost 18,000 wells) could lead to the conclusion that there is a relationship between geology and well position. However, since 1992, almost 13,000 wells have been drilled (i.e., 72% of the total number of wells existing in September 2005), out of which less than 10% were above newly discovered reservoirs and 48% were above known reservoirs. Therefore, at a landscape scale, well positioning does not correlate well with oil and gas discoveries, and the location of known petroleum reservoirs was not incorporated into the IK estimates. For each area an indicator and a robust variogram were determined. The type of variogram did not impact the IK estimations, as kriging is a unique and exact interpolator (Kitanidis 1997), but the confidence interval of the predicted values were affected (Cressie 1993). The robust variogram (Eq. 4.5) presented almost 50% less variability than the indicator variogram (Eq. 4.4) so, following the recommendation of Hawkins and Cressie (1984), I used only the robust variogram in the kriging computations (Table 4.2). With a probability larger than 0.75, two wells were separated by more than 300 m, To ensure that the grid used for kriging is integrated within the grid used by the Government of British Columbia, the side of the grid was selected as 400 m, half the side of the quarter-section. Support for the choice of this distance was provided by the oil and gas development occurring in the Doig area, where 665 quarter- sections have already reached the maximum number of wells recommended by current legislation (2 wells/quarter section). Table 4.2 Variograms of the study areas Study area Variogram type Equation Fit assessment Grid size [m] Isotropy P(|Q1||=0) P(|Q2|=1) p(vario. is isotropic) Moberly Spherical 0.05×(1.5×h/4000- 0.5×(h/4000) 3) 0.17 0.23 400 0.061 Doig area Spherical 0.063×(1.5×h/15000- 0.5×(h/15000) 3) 0.22 0.16 400 0.25 Fort Nelson Linear 2×10-6×h 0.19 0.13 400 0.14 The distribution of the wells within the quarter-section for the entire northeastern British Columbia follows a pattern of 35%-21%-23%-21%, clockwise from the NE quadrant to the SE quadrant. The general pattern existing in the British Columbian part of the Western Canadian 69 Sedimentary Basin was confirmed by one of the study areas (i.e., Doig) that indicated that the evolution of the drilling process within the quarter section seems to follow the 35-21-23-21 pattern. For the Moberly and Fort Nelson areas, areas with reduced oil and gas developments, the spatial distribution was almost equal among the four quadrants (Fig. 4.4). Therefore, the probability of drilling within the quadrant was assigned according to the distribution characterizing the drilling occurrence for the entire northeastern British Columbia. Figure 4.4 Quadrants within the quarter section grid used to determine the well forecast grid The average forest harvesting for the next 100 years varied according to the scheduling algorithm and the harvesting age (Fig. 4.5), with the FFD algorithm supplying the largest AAC for all three areas. The selection of the harvesting age played a significant role in the size of AAC which was not necessarily translated in a significant disturbance at the landscape level (Table 4.3), as for the Doig and Moberly areas the difference between the harvested areas was less than 15% , regardless of the scheduling algorithm. However, the algorithm plays a significant role in determining the magnitude of the AAC, as the FFD algorithm supplied results significantly larger (p<0.001) than SA, from 15% for the Fort Nelson area to more than 20% for the Doig and Moberly areas. The results for Fort Nelson are counter-intuitive as despite the significantly larger AAC (i.e., p<0.0001 for AAC=270000m3/year supplied by FFD algorithm compared to 225000m3/year supplied by SA algorithm) the clearing of the forested landscape was insignificant (i.e., p>0.2 when the clear-cut area of 950 ha/year determined by the FFD algorithm was compared to the clear-cut area of 930 ha/year for the SA algorithm). In fact, in 85% of the futures the larger AAC was not reflected in a significantly larger landscape disturbance. 70 Figure 4.5 The evolution of the AAC and of the land disturbance during the planning period Table 4.3. Timber harvesting land base (THLB) and AAC for the four harvesting scenarios Study area Surface [ha] # stands >0.5 ha THLB [ha] Average AAC [1000 m3] / Average Harvested area [ha] FFD/MAI FFD/Preset SA/MAI SA/Preset Doig 684 138 37042 183 057 495 390 355 330 1500 1350 1300 1200 Fort Nelson 641 024 86893 91 355 270 220 225 120 950 950 930 900 Moberly 410 194 36319 140 222 290 155 185 205 1050 1000 700 670 The coefficient of variation for the AAC lay between 5%, for the Fort Nelson area and 30%, for the Moberly area. The reduced AAC variability, the direct result of the 10% AAC annual variation, is not reflected at the scale of individual clear-cuts, which has a coefficient of variation (CV) almost double than the corresponding coefficient of variation of the AAC (e.g., CV AAC – FFD/Age@max(MAI) for Doig = 29% and CV cleared area – FFD/ Age@max(MAI) for Doig = 53% or CV AAC – SA/Preset age for Moberly = 6% and CV cleared area – SA/Preset age for Moberly = 21%). 71 The required seral stage distribution [i.e., stands having high biodiversity emphasis and less than 40 years old should occupy less than 27% of the study area, and the old seral stages (stands older than 120 years) should occur on at least a third of the area (Ministry of Forests 1995)] was met throughout the planning period for all three areas. The fulfillment of the seral stage distribution requirements is obvious for the futures which have the harvests occurring at ages larger than 120 years. The futures for which the harvests occur at the age when MAI peaks also met the seral stage distribution requirements as more than 30% of the stand have the Agemax(MAI)>120 years. Irrespective the scheduling algorithm and the adopted harvesting age, the clearcuts did not change the landscape structure more than 2% year to year, with harvested stands being regenerated with species appropriate to the respective forest ecosystems. 4.3.2 Habitat suitability index for the set of futures The combination of two forest scheduling algorithms, two harvesting ages, five drilling models, two species and two habitat perspectives (i.e., average HSI and stands with HSI>0.5) led to 80 futures for each study area. The five drilling models had an undifferentiated impact on the average HSI (i.e., the difference among the five models was less than 1%), regardless of which scheduling algorithm or harvesting age was used, or of which study area. Therefore, for each study area and VEC, the futures quantified using average HSI were determined only by the forest harvesting (Fig. 4.6). Different results were obtained for the futures represented by stands with HSI>0.5, as drilling activities could play an important role for the first 20 to 40 years, when significantly different environments could be reached, depending on the study area and the model describing petroleum development (Fig 4.7a and 4.7b). For each study area, forest scheduling algorithm and harvesting age the set of five futures converged to one future once the maximum number of wells was reached. After the moment of convergence, the dynamics of the stands with HSI>0.5 were determined entirely by the forest harvesting activities. 72 Figure 4.6 The evolution of the average HSI during the planning period The FFD algorithm supplied better ranges for the average HSI than did the SA algorithm, regardless of the method of selecting the harvesting age (i.e., age at maximum MAI or preset age). The FFD algorithm also led to a relatively constant range of variation for the average HSI (e.g., 4 % for moose HSI), while SA induced larger changes in the average HSI (e.g., 12% for moose). While the harvesting algorithm seems to play a significant role in the dynamics of the average HSI for moose (i.e., 46%–48% for FFD compared with 32%–44% for SA in the Moberly area, 46%–50% compared to 44%–46% in the Fort Nelson area and 64%–68% compared to 49%–52% in the Doig area), the American martin habitat suitability seems to be little influenced by the algorithm, as the HSI had similar pair-wise ranges for each study area (i.e., 34%–38% for the Moberly area, 41%–47% for the Fort Nelson area and 23%–26% for the Doig area). Besides the differences in ranges (which is a measure of the variation in the L1 space), the algorithms also presented different measures of variation when the L2 space was used to assess the variability of the futures (Grimmett and Stirzaker 2002). For each species, harvesting age and study area, the standard errors of the mean of the nonlinear equations describing the temporal dynamics of the average HSI were at least twice larger for the FFD algorithm than for the SA algorithm (i.e., 3.7 times for the Moberly area, 2.1 times for Fort Nelson area and 4.1 times for 73 the Doig area). The equations used to describe the evolution of the average HSI [the equations were based on Bertalanffy’s family of equations )exp()( 1653210 4 −×−××+×+ yearbbyearbbbb b , where b0 –b6 are coefficients] is of little interest for the proposed CIA framework, as it aims to investigate a set of futures, rather than any one in particular. Nevertheless, each equation fulfilled the time-series modeling assumptions and requirements [the equations were significant (p<0.001), had a coefficient of determination larger than 0.9, and the residuals were organized as white noise (p>0.2) with no evidence of nonlinear trends (p>0.08)]. The method of selecting the harvesting age within the forest scheduling algorithm seems to play a significant role in the dynamics of the average habitat suitability as for almost 85% of the comparisons (10 out of 12) the HSI had different evolutions depending on whether the harvesting occurred at MAI or after a preset age. The method of selecting the harvesting age could also play an important role in the long-term evolution of the habitat suitability, as for a preset harvest age the average HSI had an increasing tendency, regardless of the algorithm, species and study area. For the harvesting age as determined by the MAI, half the futures exhibited a relatively constant HSI, and half decreased. The most significant decrease in the average HSI occurred in areas with a reduced timber harvesting land base relative to the size of the entire study area (14% of the entire Fort Nelson area and 27% for the Doig area), and increased variation in the surface cleared annually. However, regardless of the scheduling algorithm and harvesting age the average HSI did not decrease more than 3% during the entire 100-year planning period. While the average HSI decreased in six out of the 24 futures, the area with HSI>0.5 either remained relatively constant (less than 1% decrease) or increased during the planning period. Like the average HSI, the FFD algorithm supplied larger areas with HSI>0.5 than the SA algorithm, regardless of the harvesting age (Fig 4.7a and 4.7b). The standard error of the mean was 30% smaller for the FFD than for the SA, the reverse of what occurred with the average HSI. For American marten, regardless of the scheduling algorithm or harvesting age, the surface of the stands with HSI>0.5 increased. The development of moose habitat was more nuanced than that of American marten, as the dynamics of the area with HSI>0.5 showed a slow but non-significant decreasing trend (in the Fort Nelson area, p>0.06). A characteristic of all the futures describing the surface of the stands with HSI>0.5 for American marten was the reduced standard error of the mean [the models describing the American marten evolution were also developed using Bertalanfy’s family of equations] when compared with moose (less than 20%). Drilling activities seemed to have no impact on the occurrence of highly suitable stands for American marten but played a major role in moose habitat. Depending on the study area, the availability of stands with HSI>0.5 for moose could decrease by more than 10000 ha (e.g., Doig 74 area) as long as petroleum drilling was active. Furthermore, the petroleum drilling seemed to lead to sharp annual changes in moose habitat, with significant reductions from year to year [up to ten times for the Fort Nelson (Fig. 4.7b) or Moberly (Fig. 4.7a) areas]. For the period when the petroleum drilling influenced the environment, the regression type models seem to have a similar impact on the better habitat of the two species, as all three models produced futures with an analogous trend and relatively reduced magnitude in the variation of the stands with HIS>0.5,. For the same period, the futures supplied by the regression types models were significantly different from futures developed using the generalized models (p<0.01), which are responsible for the largest annual drop in the area of the stands with HSI>0.5 (e.g. Fort Nelson area and the forest harvesting scheduled using FFD algorithm or Moberly area and the forest harvesting scheduled using SA algorithm). Figure 4.7 a. The evolution of the areas with HSI>0.5 using FFD as forest harvesting algorithm 75 Figure 4.7 b The evolution of the areas with HSI>0.5 using SA as forest harvesting algorithm In view of the lack of significance of the drilling activities on the average HSI, only 144 futures [3 study areas x 2 species x (5 drilling models HIS>0.5 + 1 drilling model average HSI) x 2 forest scheduling algorithms x 2 harvesting ages] were used for the CIA investigations. All 144 futures had the average HSI and the surface of the stands with HSI>0.5 not decreasing by more than 15 %, indicating that each future could occur, as the selected thresholds were met. Consequently, all 144 futures were used in the CIA. The sphericity test (p<0.001) showed that the covariance matrix did not satisfy the Huynh – Feldt condition, and the results were therefore assessed using the correction proposed by Greenhouse and Geisser (1959) and Huynh and Feldt (1970). The two tests showed that the environment would change significantly during the next century (p<0.0001), regardless of the study area (p<0.001) or harvesting algorithm (p=0.002). A similar conclusion was reached for the differences between the futures, but with a reduced level of significance (p=0.02). The forest scheduling algorithm seems to be responsible 76 for the differences among the futures (p=0.038) while the study area has little impact (p=0.48) in identifying a significant different future within the set of 144 futures. The scheduling algorithm showed different influence on the species, as moose seems to be associated with both SA and FFD (47% of futures with SA and 97% with FFD) while American marten mainly with FFD (46% of the futures). An adjusted Scheffe test separated the 100 years into three periods: year 1 to 50, year 51 to 93 and year 94 to 100. The periods differed not only in duration but also in variability [seven distinct time intervals spanned the first 50 years (year 1, years 2–11, years 12–16, years 17–26, years 27–36, years 37–39 and years 40–50), while the second period had only one time interval and the third had two (years 94–98 and 99–100)]. The time intervals had significantly different covariance matrixes, but the HSI had similar trends within each period increasing in the first period, constant in the second and decreasing in the last period. The models used to develop the futures were significantly different (p<0.0001), with the set of futures developed for moose based on the SA algorithm and HSI>0.5 in Fort Nelson being separated for the remaining 134 futures. The SA algorithm was also responsible for the further separation of the 134 futures into two groups: one having 124 futures and one having the 10 futures derived for American marten in Fort Nelson area. The delineation of the 144 futures in two groups was maintained regardless of the forest harvesting algorithm: one grouping the models describing habitat suitability for American marten in the Ft. Nelson area with the models for moose in the Moberly area (i.e., 30 futures), and one containing the remaining 114 futures. The investigation of the autocorrelation function for all 144 futures revealed the existence of at least one cycle within each future. The cycle periods exceeded 100 years, ranging from 100 to 140 years, with two exceptions, both of which were 30 years (both quantifying the average HSI for moose for the Doig area using the SA algorithm). About 18% of the futures revealed a second, shorter cycle imbedded into the longer cycle, with period varying from 15 to 30 years. 4.4 Discussion The results supplied by the set of 144 futures confirmed the studies performed in the region by Nitschke (2008) and Schneider et al (2003) that found that the combined environmental impacts of forestry and oil and gas development in the western Sedimentary Canadian Basin were not significant. Nitschke (2008) identified that forest harvesting and petroleum drilling activities (i.e., well drilling, roads and resource exploratory clearing) were the main activities changing the landscape in northeastern British Columbia. Among the 82 species investigated by Nitschke 77 (2008), less than a quarter (i.e., 20 species) could have encompassed a reduction in the surface of their habitat, including American marten and moose. The decrease in the habitat availability for moose that Nitschke suggested might have occurred in northeastern British Columbia contradicts the findings of Ball and Dahlgren (2002), who found that the human activities increased the area of suitable habitat. In this study, I found that habitat suitability for moose did not change significantly (Fig 4.6 and 4.7a and b). The models suggested that an increase in habitat suitability would occur for American marten, regardless of the type of measurement used (i.e., average HSI or stands with HSI>0.5), confirming the results of Chapin et al (1998) and Hargis et al (1999). The habitat suitability of both species responded primarily to forest harvesting rather than petroleum drilling, especially when the average HSI was used to evaluate changes in the environment (Fig 4.6). However, when better habitats were considered in the analysis, petroleum drilling started to induce significant changes in the environment (Fig.4.7). These changes were evident before the number of wells peaked. The set of functions describing the evolution of the drilling effort is separated into two classes with or without a periodic term. The increasing amplitude part of the period term seems to emphasize the rise in oil and gas drilling sensitivity, probably a result of the increase in stock- market dynamics combined with the fluctuating geo-economical situation (Williams 2006). However, the decreasing period (i.e., from 12 years to 8 years), starting in 1988, indicates the fast response of oil and gas industry to the socio-economic environment as well as the rapid implementation of petroleum related research (Cuddington and Moss 2001). Forest planning activities aim to maximize, minimize or maintain an objective, which can be explicitly quantitative, such as AAC (Davis et al. 2001) or net present value (Pukkala and Kurttila 2005), or can be represented numerically, such as wildlife habitat (Bettinger et al. 2002) or recreation (Teeguarden and Werner 1968). To avoid subjective quantifications that could influence the results, the aim of the forest planning process in all three areas was AAC maximization and not habitat enhancement or increase of habitat surface. The results revealed that irrespective the algorithm or the harvesting age, both HSI statistics (i.e., average and areas with HSI>0.5) exhibited a pronounced cyclicity. Furthermore, among the set of 144 futures with a significant long term cycle, 24 futures presented a second cycle, enforcing the synergistic impact of the human activities on the environment, as noncyclical developments could induce cycles on the VECs. While the longer cycle could be associated with the variability in the 78 environment, the cycles with shorter periods indicate that when considered together, timber harvesting and well drilling could periodically place the environment at the risk of reaching undesired states, especially when the attributes quantifying VECs dynamics are close to the critical thresholds. The two human activities considered in the study could lead to a decrease in the amount of better habitat, especially when forest harvesting was scheduled using SA (i.e., a third of the futures); however, larger harvesting volumes do not necessarily translate into a decrease in the habitat suitability (i.e., the case of American marten when harvestings were scheduled using adjusted FFD algorithm). From this perspective, scheduling the harvesting using the FFD algorithm produced at least 10 % more volume / ha than SA (for the Moberly area it was almost 20%) without a significant reduction in the availability of better or average habitats. The decrease apparent in the average HSI for 60% of the futures is consistent with the findings of Jager et al (2006), which argued that the average HSI is a weak indicator of the habitat suitability at the landscape level. A reduction in both HSI statistics only occurs after the cessation of drilling activities, and is not constant. The negative trend of the HSI after 40 years is closer to the observed evolution of the environment than to the forecasts commonly made in forest management plans, which assume that the attributes of interest (e.g., AAC) have a desirable long-term trajectory, but only after an initial unfavorable period. Another difference between the results and many current management plans lies in the reduced impact of the two human activities on the availability of better habitat for moose and American marten. While this seems to contradict several earlier studies in the region (Cizek et al. 2002; Dube et al. 2006; Yamasaki et al. 2008), it is not really surprising considering that the main activity changing the landscape is forest harvesting and not petroleum drilling [in agreement with the results of (Nitschke 2008)]. One of the surprising findings of this study is the reduced impact of forest harvesting on the habitat suitability while the AAC is increasing (Table 4.3), a direct result of the arrangement of the harvests. The three periods identified by the proposed CIA analytical platform are characterized by different dynamics, the first having the largest variation and the second the lowest. The first period, which could last 50 years, comprised seven significantly different sub-periods, three of them less than five years. These three sub-periods were characterized by a relatively constant or insignificant decrease of the HSI and separate the remaining four sub-periods, which are characterized by a significant increase of the HSI. The large number of sub-periods in the first half century, which is when petroleum drilling was occurring, could indicate that the clearing for 79 well pads is not responsible for the degradation of the habitat but does increase the variability of habitat quality. Of particular interest is the short sub-period between year 37 and year 39 when petroleum drilling seems to have no impact in comparison with the disturbance induced by the harvesting scheduled with the FFD algorithm and preset age (i.e., the Moberly area in Fig 4.7a). The weak performance of the FFD algorithm when harvesting occurs at a preset age is the result of the violation of the perfect bin - packing theorem (Coffman et al. 2000), and results in the FFD supplying results worse than the SA (Strimbu et al. 2008). The second period, which is the longest period with no sub-periods, would occur most likely after the petroleum drilling stopped, and is reflected in 80% of the futures (i.e., 116 futures), indicating the importance of forest harvesting in determining environmental dynamics. The second period also has the least variation, indicating that besides the constancy of habitat suitability, the variation of the HSI and the area with better habitats would not present large deviations within the period. In the third period the area of stands with HSI>0.5 decrease in all the futures produced using the SA algorithm and in two third of the futures produced with the FFD algorithm. The delineation of the three periods along the planning period follows a sinusoidal pattern that could be associated with the period cycle. The cyclic feature of the HSI emphasizes the long term and difficult to detect cumulative impact of the two activities. The study indicates that cyclical response of the environment to different human induced disturbances seems to be produced by non-cyclical human activities. This type of response resembles the chaotic behavior identified by May (1974) or Williams (1997). The three periods are separated by two significant moments, one most likely to occur after 50 years, and one after 93 years. The HSI describing the two moments do not differ significantly in magnitude (p=0.24) as both could be associated with the second period, a period characterized by reduced variability of the measures used to quantify the habitat. However, the two moments can also be seen as separating an upward trend (i.e., first period) from a downward trend (i.e., third period). The first moment would occur at least 10 years after the cessation of petroleum drilling (after 30 years in the Doig area) confirming the environmental inertia mentioned by Hannan and Freeman (1984) and Ruef (1997). The first moment can be relatively easy identified as the moment when the increase in habitat stops, but the second moment cannot be clearly delineated, as the downward trend is present in 80% of the futures. Furthermore, the decrease starts after 80 years (sometimes after 70 years [i.e., in the Moberly and Fort Nelson areas when the forest harvesting is scheduled using the FFD algorithm)] and not after 93 years, as the longitudinal analysis indicated the presence of a statistically significant change in the set 80 of 144 futures. The 13-year difference is probably associated with the moment when the effects could be detected (i.e., year 93) and not necessarily when the environment started to change (i.e., year 80). The attributes that could induce significant changes in the environment are harvesting age and the VEC. The increasing harvesting age could significantly reduce the availability of better habitats, especially when SA is used to schedule the forest harvesting. The VEC could also change the CIA investigation as the measure for American marten seems to have a constant increase (i.e. 58% of the futures quantifying the availability of better habitats) while the measure for moose decreased (i.e. 54% for the futures quantifying the availability of better habitats). While the choice of harvesting age could be responsible for significant changes in the environment, the algorithm used to schedule the harvesting seemed to have little impact on landscape disturbance, as five out of six scenarios induced similar clearing footprints on the landscape. Regardless of the development of the petroleum drilling, the clearing associated with the well pads did not produce significant changes in the metrics used to assess the VEC. From this perspective, the futures obtained using the five models (i.e., the three regression-type models and the two generalized models) could have not been delineated for the same study area and within the same forest harvesting algorithm and harvesting age. 4.5 Conclusion Traditionally, the inclusion of future activities in cumulative impact assessments has been done using one analytical platform. For complex CIA investigations, the analytical platform is commonly based on modular modeling (Alila and Beckers 2001; Voinov et al. 2004) or on different operations research techniques (Stakhiv 1988). To compensate for the influence of the model, I have suggested a shift in the paradigm governing the CIA. The paradigm shift involves a change in the focus of CIA investigations from the detailed analysis of one unlikely future to the identification of the patterns describing the future changes in the environment. To illustrate the approach, a set of 144 possible and equally likely futures were developed that aimed to identify the potential impacts of forest harvesting and petroleum drilling on the habitat suitability of moose and American marten. Petroleum drilling will cease by the middle of the 21st century, with areas experiencing higher drilling intensity in the past (e.g., Doig) stopping before areas that are currently less well- developed (e.g., Moberly). The economic and environmental impact of the forest harvesting is 81 dependent on the algorithm used to schedule the harvestings as well as on the choice of harvesting age. The combination of the first-fit decreasing algorithm and harvesting at the age when MAI peaks supplied the largest AAC among all the planning strategies investigated. In 85% of cases, the disturbance associated with larger AAC did not differ significantly from lower AACs. The evolution of two measures of habitat suitability (average HSI and surface of the stands with HSI>0.5) revealed that the human activities could induce at least one cycle, having a period between 100 years and 150 years, in the habitat dynamics of moose and American marten, even though regulatory constraints do not plan for cycles. The existence of a cycle in the habitat of the two VEC was reinforced by the presence of a second but shorter cycle in a sixth of the futures. However, the results support the conclusion that habitat suitability would not vary by more than 15% from the present values throughout the 100-year simulation period. The two human activities considered in the CIA had no negative impact on the availability of better habitat for American marten but could lead to an insignificant reduction for moose. The planning period was separated into three distinct periods: the first characterized by an increase in the habitat suitability measures, which would last almost 50 years, the second lasting approximately 40 years and characterized by constant conditions, and the third (10-year) period marking the start of a decrease in the habitat suitability measures. The three periods were separated by two significant moments, one most likely to occur after 50 years, and one during the last decade of the century. The two moments could not be separated based on the magnitude of the attributes describing the periods, but the moments delineated an upward trend (i.e., first period) from a downward trend (i.e., third period). Both moments would occur after the cessation of petroleum drilling. The attributes that could induce significant changes in the environment are the choice of harvesting age and the VEC. The choice of VEC is critical to the analysis and could change the conclusions of the CIA. The study suggests that the harvesting algorithm and petroleum drilling may not have a significant impact on the VEC used in the CIA investigation. However, I recommend that future CIA investigations be based on the analysis of a greater range of human activities and a greater range of VECs. 82 4.6 References 1. Agricultural Land Commission and Oil and Gas Commission. Agricultural Land Commission Act. 36. 2002. 2. Alila, Y. and J. Beckers. 2001. Using numerical modelling to address hydrologic forest management issues in British Columbia. Hydrological Processes 15(18):3371-3387. 3. Allen, A.J., J.W. Terrell, and P.A. Jordan. 1988. An overview of a habitat suitability index model for moose: Lake superior region. Alces 24:118-125. 4. Anderson, T.W. 1984. An Introduction to Multivariate Statistical Analysis. John Wiley and Sons, New York. 704 p 5. Anonymous. Maximum Disturbance Review Criteria. 2004. Fort St. John, Oil and Gas Commission. 8 p 6. Appel, K. and W. Haken. 1976. Every Planar Map Is 4 Colorable. Bulletin of the American Mathematical Society 82(5):711-712. 7. AXYS Environmental Consulting Ltd. A cumulative effects assessment and management framework for northeast British Columbia. 2003. Fort St. John, AXYS Environmental Consulting Ltd. 363 p 8. B.C.Conservation Data Centre. BC Species and Ecosystems Explorer. http://a100.gov.bc.ca/pub/eswp . 2008. B.C. Ministry of Environment Victoria BC. 10-21-2008. 9. Bailey, S.-A., R.H. Haines-Young, and C. Watkins. 2002. Species presence in fragmented landscapes: modelling of species requirements at the national level. Biological Conservation 108(3):307-316. 10. Balke, N.S. and R.J. Gordon. 1986. Historical data. P. 781-850 in The American Business Cycle, Gordon, R.J. (ed.), The University of Chicago Press, Chicago IL. 11. Ball, J.P. and J. Dahlgren. 2002. Browsing damage on pine (Pinus sylvestris and P- contorta) by a migrating moose (Alces alces) population in winter: Relation to habit at composition and road barriers. Scandinavian Journal of Forest Research 17(5):427-435. 12. Beanlands, G.E. and P.N. Duinker. 1983. An ecological framework for environmental impact assessment in Canada. Institute for Resources and Environmental Studies, Dalhousie University, Halifax. 132 p 13. Belsley, D.A., E. Kuh, and R.E. Welsch. 1980. Regression Diagnostics: Identifying Influential Data and Sources of Collinearity. John Wiley and Sons, New York. 292 p 14. Bettinger, P., D. Graetz, K. Boston, J. Sessions, and W.D. Chung. 2002. Eight heuristic planning techniques applied to three increasingly difficult wildlife planning problems. Silva Fennica 36(2):561-584. 83 15. Boltzmann, L. 1995. Lectures on gas theory. Dover Publications. 490 p 16. Borland Software Corporation. Delphi. [7.0]. 2001. Cupertino, CA. 17. Box, G.E.P., G.M. Jenkins, and G.C. Reinsel. 1994. Time Series Analysis: Forecasting and Control. Prentice Hall, Englewood Cliffs, NJ. 18. Brockwell, P.J. and R.A. Davis. 1996. An Introduction to Time Series and Forecasting. Springer Verlag, New York. 420 p 19. Cameron, A.C. and P.K. Trivedi. 1998. Regression analysis of count data. Cambridge University Press, New York NY. 432 p 20. Caplat, P., M. Anand, and C. Bauch. 2008. Symmetric competition causes population oscillations in an individual-based model of forest dynamics. Ecological Modelling 211(3-4):491- 500. 21. Chapin, T.G., D.J. Harrison, and D.D. Katnik. 1998. Influence of landscape pattern on habitat use by American marten in an industrial forest. Conservation Biology 12(6):1327-1337. 22. Chaubey, I., A.S. Cotter, T.A. Costello, and T.S. Soerens. 2005. Effect of DEM data resolution on SWAT output uncertainty. Hydrological Processes 19(3):621-628. 23. Cherp, A., A. Watt, and V. Vinichenko. 2007. SEA and strategy formation theories: From three Ps to five Ps. Environmental impact assessment review 27(7):624-644. 24. Choi, B. 1992. ARMA Model Identification. Springer-Verlag, New- York. 25. Cizek, P., J. McCullum, and A. Booth. 2002. Cumulative impacts in the Fort Liard area. Yellowknife, Cizek Environmental Services. 16 p 26. Cleveland, W.S. 1972. The Inverse Autocorrelations of a Time Series and Their Applications. Technometrics 14(2):277-293. 27. Coffman, E.G., C. Courcoubetis, M.R. Garey, D.S. Johnson, P.W. Shor, R.R. Weber, and M. Yannakakis. 2000. Bin packing with discrete item sizes, Part I: perfect packing theorems and the average case behavior of optimal packing. SIAM Journal of Discrete Mathematics 13(3):384-402. 28. Coffman, E.G., C. Courcoubetis, M.R. Garey, D.S. Johnson, P.W. Shor, R.R. Weber, and M. Yannakakis. 2002. Perfect packing theorems and the average-case behavior of optimal and online bin packing. Siam Review 44(1):95-108. 29. Contini, S. and A. Servida. 1992. Risk analysis and environmental impact studies. P. 79- 104 in Environmental impact assessment, Colombo, A.G. (ed.), Kluwer Academic Publisher, Dordrecht. 30. Costanza, R., A. Voinov, B. Roelof, T. Maxwell, F. Villa, L. Wainger, and H. Voinov. 2002. Integrated Ecological Economic Modeling of the Patuxent River Watershed, Maryland. Ecological Monographs 72(2):203-231. 84 31. Council on Environmental Quality. The National Environmental Policy Act. 42. 1969. 32. Cressie, N. 1993. Statistics for Spatial Data. John Wiley&Sons, New York NY.1-900 33. Crowder, M.J. and D.J. Hand. 1990. Analysis of repeated measures. Chapman and Hall, London. 257 p 34. Crowe, K. and J. Nelson. 2003. An indirect search algorithm for harvest-scheduling under adjacency constraints. Forest Science 49(1):1-11. 35. Cuddington, J.T. and D.L. Moss. 2001. Technological change, depletion, and the US petroleum industry. American Economic Review 91(4):1135-1148. 36. Davis, L.S., K.N. Johnson, T.E. Howard, and P. Bettinger. 2001. Forest management. McGraw-Hill, New York. 804 p 37. Dee, N., J. Baker, N. Drobny, K. Duke, I. Whitman, and D. Fahringe. 1973. Environmental Evaluation System for Water Resource Planning. Water Resources Research 9(3):523-535. 38. Dickert, T.G. and A.E. Tuttle. 1985. Cumulative impact assessment in environmental planning; a coastal wetland watershed example. Environmental impact assessment review(5):37-64. 39. Dickey, D.A. and W.A. Fuller. 1979. Distribution of the Estimators for Autoregressive Time Series With a Unit Root. Journal of the American Statistical Association 74(366):427-431. 40. Dube, M., B. Johnson, G. Dunn, J. Culp, K. Cash, K. Munkittrick, I. Wong, K. Hedley, W. Booty, D. Lam, O. Resler, and A. Storey. 2006. Development of a New Approach to Cumulative Effects Assessment: A Northern River Ecosystem Example. Environmental Monitoring and Assessment 113(1 - 3):87-115. 41. Duinker, P.N. and L.A. Greig. 2007. Scenario analysis in environmental impact assessment: Improving explorations of the future. Environmental impact assessment review 27(3):206-219. 42. Dussault, C., R. Courtois, and J.P. Ouellet. 2006. A habitat suitability index model to assess moose habitat selection at multiple spatial scales. Canadian Journal of Forest Research- Revue Canadienne de Recherche Forestiere 36(5):1097-1107. 43. Energy Information Administration. Spot Prices for Crude Oil and Petroleum Products. 2005. U.S. Department of Energy. 44. Environmental Systems Research Institute, I. ArcGIS. [9.3]. 2008. Redlands CA, ESRI. 45. Franzmann, A.W. and C.C. Schwartz. 1997. Ecology and Management of the North American Moose. Smithsonian Institution Press, Washington DC. 733 p 46. Fryxell, J.M., J.B. Falls, E.A. Falls, R.J. Brooks, L. Dix, and M.A. Strickland. 1999. Density dependence, prey dependence, and population dynamics of martens in Ontario. Ecology 80(4):1311-1321. 85 47. Fuhrer, E. 2000. Forest functions, ecosystem stability and management. Forest Ecology and Management 132(1):29-38. 48. Gilbert, W.J. 1976. Modern algebra with applications. John Wiley and Sons, New York. 348 p 49. Glasson, J., R. Therivel, and A. Chadwick. 1994. Introduction to environmental impact assessment. UCL Press, London.1-342 50. Gore, J.A., J.R. Kelly, and J.D. Yount. 1990. Application of Ecological Theory to Determining Recovery Potential of Disturbed Lotic Ecosystems - Research Needs and Priorities. Environmental management 14(5):755-762. 51. Gosselink, J.G. and L.C. Lee. 1989. Special Issue - Cumulative Impact Assessment in Bottomland Hardwood Forests. Wetlands 9:93-174. 52. Government of British Columbia. Petroleum and Natural Gas Act. 362. 2004. 53. Graham, R.T., T.M. Quigley, and R. Gravenmier. 2000. An integrated ecosystem assessment of the Interior Columbia Basin. Environmental Monitoring and Assessment 64(1):31-40. 54. Green, D.G., N. Klomp, G. Rimmington, and S. Sadedin. 2006. Complexity in Landscape Ecology. Springer, Dordrecht, The Netherlands.1-208 55. Greenhouse, S. and S. Geisser. 1959. On methods in the analysis of profile data. Psychometrika 24(2):95-112. 56. Grimmett, G.D. and D.R. Stirzaker. 2002. Probability and Random Processes. Oxford University Press, New York.1-600 57. Hannan, E.J. and J. Rissanen. 1982. Recursive Estimation of Mixed Autoregressive- Moving Average Order. Biometrika 69(1):81-94. 58. Hannan, M.T. and J. Freeman. 1984. Structural Inertia and Organizational Change. American Sociological Review 49(2):149-164. 59. Hargis, C.D., J.A. Bissonette, and D.L. Turner. 1999. The influence of forest fragmentation and landscape pattern on American martens. Journal of Applied Ecology 36(1):157-172. 60. Hawkins, D.M. and N. Cressie. 1984. Robust kriging -A proposal. Mathematical Geology 16(1):3-18. 61. Hegmann, G., C. Cocklin, R. Creasey, S. Dupuid, A. Kennedy, L. Kingsley, W. Ross, H. Spaling, and D. Stalker. Cumulative Effects Assessment Practitioners Guide. 1-143. 1999. Hull, QC, Canadian Environemntal Assessment Agency. 62. Hirzel, A.H., G. Le Lay, V. Helfer, C. Randin, and A. Guisan. 2006. Evaluating the ability of habitat suitability models to predict species presences. Ecological Modelling 199(2):142-152. 86 63. Holland, S.P. 2008. Modeling Peak Oil. Energy Journal 29(2):61-79. 64. Hungerford, T.W. 1990. Abstract algebra: an introduction. Saunders College Publishing, Orlando.1-575 65. Huynh, H. and L.S. Feldt. 1970. Conditions Under Which Mean Square Ratios in Repeated Measurements Designs Have Exact F-Distributions. Journal of the American Statistical Association 65(332):1582-1589. 66. Iledare, O.O. and A.G. Pulsipher. 1997. The modelling of the petroleum exploration and extraction process for policy analysis: a case study of the Louisiana onshore region. Pacific and Asian Journal of Energy 7(1):21-38. 67. Iledare, O.O. 1995. Simulating the effect of economic and policy incentives on natural gas drilling and gross reserve additions. Resource and Energy Economics 17(3):261-279. 68. Information Services AEUB. Drilling and well Information. 2005. Calgary AB, Alberta Energy and Utility Board. 69. Jager, H.I., E.A. Carr, and R.A. Efroymson. 2006. Simulated effects of habitat loss and fragmentation on a solitary mustelid predator. Ecological Modelling 191(3-4):416-430. 70. Johnson, D.S. 1973. Near-optimal bin packing algorithms. PhD Massachusetts Institute of Technology, Dept. of Mathematics. 402 p. 71. Journel, A.G. 1983. Nonparametric estimation of spatial distributions. Mathematical Geology 15(3):445-466. 72. Jukes, W., D. Menzies, D. Rosen, G. Taylor, R. Saint-Jean, and J. Beale. 2004a. Sustainable Forest Management Plan. Fort St. John, Ministry of Forests. 402 p 73. Jukes, W., D. Menzies, D. Rosen, G. Taylor, R. Saint-Jean, and J. Beale. 2004b. Sustainable Forest Management Plan - Appendices. Fort St. John, Ministry of Forests. 531 p 74. Kienzle, S. 2004. The Effect of DEM Raster Resolution on First Order, Second Order and Compound Terrain Derivatives. Transactions in GIS 81(1):83-111. 75. Kingman, J.F.C. 1995. Poisson Processes. Clarendon Press, Oxford. 104 p 76. Kitanidis, P.K. 1997. Introduction to geostatistics. Cambridge University Press, Cambridge. 249 p 77. Kitanidis, P.K. 1991. Orthonormal residuals in geostatistics: Model criticism and parameter estimation. Mathematical Geology 23(5):741-758. 78. Koornneef, J., A. Faaij, and W. Turkenburg. 2008. The screening and scoping of Environmental Impact Assessment and Strategic Environmental Assessment of Carbon Capture and Storage in the Netherlands. Environmental impact assessment review 28(6):392-414. 87 79. Korzukhin, M.D., M.T. TerMikaelian, and R.G. Wagner. 1996. Process versus empirical models: Which approach for forest ecosystem management? Canadian Journal of Forest Research-Revue Canadienne de Recherche Forestiere 26(5):879-887. 80. Kroll, A.J. and J.B. Haufler. 2006. Development and evaluation of habitat models at multiple spatial scales: A case study with the dusky flycatcher. Forest Ecology and Management 229(1-3):161-169. 81. Lam, D., L. Leon, S. Hamilton, N. Crookshank, D. Bonnin, and D. Swayne. 2004. Multi- model integration in a decision support system: a technical user interface approach for watershed and lake management scenarios. Environmental modelling and software 19:317-324. 82. Lanza, A., M. Manera, M. Grasso, and M. Giovannini. 2005. Long-run models of oil stock prices. Environmental Modelling & Software 20(11):1423-1430. 83. Larson, M.A., F.R. Thompson, J.J. Millspaugh, W.D. Dijak, and S.R. Shifley. 2004. Linking population viability, habitat suitability, and landscape simulation models for conservation planning. Ecological Modelling 180(1):103-118. 84. Leimbach, M. and C. Jaeger. 2005. A modular approach to Integrated Assessment modeling. Environmental Modeling and Assessment 9(4):207-220. 85. Lenzen, M., C.J. Dey, and S.A. Murray. 2004. Historical accountability and cumulative impacts: the treatment of time in corporate sustainability reporting. Ecological Economics 51(3- 4):237-250. 86. Leopold, L.B., F.E. Clarke, B.B. Hanshaw, and J.R. Balsley. 1971. A procedure for evaluating environmental impact. 645. Washington D.C, U.S. Geological Survey. 15 p 87. Lindsey, J.K. 1995. Fitting Parametric Counting Processes by Using Log-Linear Models. Applied Statistics 44(2):201-212. 88. Liu, C.R., P.M. Berry, T.P. Dawson, and R.G. Pearson. 2005. Selecting thresholds of occurrence in the prediction of species distributions. Ecography 28(3):385-393. 89. Lockwood, C. and T. Moore. 1993. Harvest scheduling with spatial constraints: a simulated annealing approach. Canadian Journal of Forest Research 23(3):468-478. 90. MacDonald, L.H. 2000. Evaluating and managing cumulative effects: Process and constraints. Environmental management 26(3):299-315. 91. Manning, N. 1991. The Uk Oil Industry - Some Inferences from the Efficient Market Hypothesis. Scottish Journal of Political Economy 38(4):324-334. 92. Marcu, M. 1981. Meteorologie cu elemente de climatologie forestiera. Ceres, Bucharest. 291 p 93. Marr, K. 1997. Environmental impact assessment in the United Kingdom and Germany. Ashgate Publishing, Aldershot. 327 p 88 94. Masera, M. and A.G. Colombo. 1992. Contents and phases of an EIA study. P. 53-78 in Environmental impact assessment, Colombo, A.G. (ed.), Kluver Academic Publishers, Dordrecht. 95. May, R.M. 1974. Biological Populations with Nonoverlapping Generations - Stable Points, Stable Cycles, and Chaos. Science 186(4164):645-647. 96. McLeod, A.I. and W.K. Li. 1983. Diagnostic Checking ARMA Time Series Models Using Squared-Residual Autocorrelations. Journal of Time Series Analysis 4:269-273. 97. Meidinger, D. and J. Pojar. 1991. Ecosystems of British Columbia. Victoria, BC Ministry of Forests. 179 p 98. Menard, A., P. Dube, A. Bouchard, C.D. Canham, and D.J. Marceau. 2002. Evaluating the potential of the SORTIE forest succession model for spatio-temporal analysis of small-scale disturbances. Ecological Modelling 153(1-2):81-96. 99. Metropolis, N., A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller. 1953. Equations of State Calculations by Fast Computing Machines. Journal of Chemical Physics 21(6):1087-1092. 100. Ministry of Forests. 1995. Biodiversity Guidebook. Victoria BC, BC Minstry of Forests. 101. Ministry of Forests. TIPSY. [3.2b]. 2005a. Victoria BC. 102. Ministry of Forests. WinVDYP. [7.26c]. 2005b. Victoria BC. 103. Ministry of Sustainable Resource Management. TRIM/TRIM II positional files - basemap. [2004], http://srmwww.gov.bc.ca/bmgs/catalog/bcgs_indexes.htm. 2004. Victoria BC, Ministry of Sustaianble Resource Management. 104. Neter, J., M.H. Kutner, C.J. Nachtsheim, and W. Wasserman. 1996. Applied linear statistical models. WCB McGraw-Hill, Boston. 1408 p 105. Nitschke, C.R. 2008. The cumulative effects of resource development on biodiversity and ecological integrity in the Peace-Moberly region of Northeast British Columbia, Canada. Biodiversity and Conservation 17(7):1715-1740. 106. Nitschke, C.R. and J.L. Innes. 2008. A tree and climate assessment tool for modelling ecosystem response to climate change. Ecological Modelling 210(3):263-277. 107. Oil and Gas Commission. Well Surface Location. 2005. Oil and Gas Commission. 108. Osowski, S.L., J.D. Swick, G.R. Carney, H.B. Pena, J.E. Danielson, and D.A. Parrish. 2001. A watershed-based cumulative risk impact analysis: Environmental vulnerability and impact criteria. Environmental Monitoring and Assessment 66(2):159-185. 109. Ottaviani, D., G.J. Lasinio, and L. Boitani. 2004. Two statistical methods to validate habitat suitability models using presence-only data. Ecological Modelling 179(4):417-443. 89 110. Patil, A.A., A.P. Annachhatre, and N.K. Tripathi. 2002. Comparison of conventional and geo-spatial EIA: A shrimp farming case study. Environmental impact assessment review 22(4):361-375. 111. Perdicoulis, A. and J. Piper. 2008. Network and system diagrams revisited: Satisfying CEA requirements for causality analysis. Environmental impact assessment review 28(7):455- 468. 112. Peterson, R.O. 1997. Ecological Studies of Wolves on Isle Royale. Houghton MI, School of Forestry and Wood Products, Michigan Technological University. 17 p 113. Phillips, P.C.B. 1986. Understanding Spurious Regressions in Econometrics. Journal of Econometrics 33(3):311-340. 114. Phillips, P.C.B. and P. Perron. 1988. Testing for a Unit Root in Time Series Regression. Biometrika 75(2):335-346. 115. Preston, E.M. and B.L. Bedford. 1988. Evaluating cumulative effects on wetland functions: A conceptual overview and generic framework. Environmental management 12(5):565-583. 116. Pukkala, T. and M. Kurttila. 2005. Examining the performance of six heuristic optimisation techniques in different forest planning problems. Silva Fennica 39(1):67-80. 117. Rao, R.D. 2000. An integrated modelling framework for exploration and extraction of petroleum resources. Resources Policy 26(3):133-143. 118. Rees, W.P. 2002. The Arrival Rate of Initial Public Offers in the UK. European Financial Management 3(1):45-62. 119. Reid, L.M. 2001. Cumulative watershed effects:then and now. Watershed management Council Network 10(1):24-33. 120. Rencher, A.C. 2002. Methods of Multivariate Analysis. John Wiley and Sons, New York. 708 p 121. Ringlund, G.B.r., K.E. Rosendahl, and T. Skjerpen. 2008. Does oilrig activity react to oil price changes An empirical investigation. Energy Economics 30(2):371-396. 122. Romito, T., K. Smith, B. Beck, J. Beck, M. Todd, R. Bonar., and R. Quinlan. Moose winter habitat. 27, 1-6. 1999. Hinton AB, Foothils Model Forest. Habitat Suitability. 123. Rotmans, J. and M.B.A. van Asselt. 2001. Uncertainty management in integrated assessment modeling: Towards a pluralistic approach. Environmental Monitoring and Assessment 69(2):101-130. 124. Ruef, M. 1997. Assessing organizational fitness on a dynamic landscape: An empirical test of the relative inertia thesis. Strategic Management Journal 18(11):837-853. 125. Sadorsky, P. 2001. Risk factors in stock returns of Canadian oil and gas companies. Energy Economics 23(1):17-28. 90 126. SAS Institute Inc. SAS Software version 9.1.3. 2004. Cary, SAS Institute Inc. 127. Schneider, R.R., J.B. Stelfox, S. Boutin, and S. Wasel. 2003. Managing the cumulative impacts of land uses in the Western Canadian Sedimentary Basin: A modeling approach. Conservation Ecology 7(1). 128. Schnorbus, M. and Y. Alila. 2004. Forest harvesting impacts on the peak flow regime in the Columbia Mountains of southeastern British Columbia: An investigation using long-term numerical modeling. Water Resources Research 40(5). 129. Spaling, H. and B. Smit. 1994. Classification and evaluation methods for cumulative effects assessment. P. 47-65 in Cumulative effects assessment in Canada: from concept to practice, Kennedy, A.J. (ed.), Alberta Association of Professional Biologists, Edmonton. 130. Stakhiv, E.Z. 1988. An Evaluation Paradigm for Cumulative Impact Analysis. Environmental management 12(5):725-748. 131. Strimbu, B.M., J.L. Innes, and V.F. Strimbu. 2008. A harvest scheduler using perfect bin- packing theorem. P. 57 Hahn, A., T. Knoke, and T. Schneider (eds.). Wald Forst Holz Zentrum, Freising DE. 132. Takats, L., R. Stewart, M. Todd, R. Bonar., J. Beck, B. Beck, and R. Quinlan. 1999. American marten winter habitat. 30. Hinton AB, Foothills Model Forest. 8 p 133. Teeguarden, D.E. and K.R. Werner. 1968. Integrating Forest-Oriented Recreation with Timber Growing - A Case Study of Economic Factors. California Agriculture 22(10):10-13. 134. Toffler, A. 1970. Future Shock. Random House, New York. 505 p 135. Tran, Z.V. 1997. Estimating Sample Size in Repeated-Measures Analysis of Variance. Measurement in Physical Education and Exercise Science 1(1):89-102. 136. Tsay, R.S. and G.C. Tiao. 1985. Use of Canonical Analysis in Time Series Model Identification. Biometrika 72(2):299-315. 137. Verburg, P.H., P.P. Schot, M.J. Dijst, and A. Veldkamp. 2004. Land use change modelling: current practice and research priorities. GeoJournal 61(4):309-324. 138. Vinod, H.D. 1973. Generalization of the Durbin-Watson Statistic for Higher Order Autoregressive Process. Communications in Statistics 2:115-144. 139. Voinov, A., C. Fitz, R. Boumans, and R. Costanza. 2004. Modular ecosystem modeling. Environmental modelling and software 19:285-304. 140. Walls, M.A. 1994. Using A Hybrid Approach to Model Oil and Gas-Supply - A Case- Study of the Gulf-Of-Mexico Outer Continental-Shelf. Land Economics 70(1):1-19. 141. Walls, M.A. 1992. Modeling and forecasting the supply of oil and gas : A survey of existing approaches. Resources and Energy 14(3):287-309. 142. Williams, G.P. 1997. Chaos theory tamed. Joseph Henry Press, Washington, D.C. 499 p 91 143. Williams, J.L. 2006. WTRG Economics. http://www.wtrg.com/index.html . 2-17-2006. 144. Wu, J.G. and J.L. David. 2002. A spatially explicit hierarchical approach to modeling complex ecological systems: theory and applications. Ecological Modelling 153(1-2):7-26. 145. Yamasaki, S.H., R. Duchesneau, F. Doyon, J.S. Russell, and T. Gooding. 2008. Making the case for cumulative impacts assessment: Modelling the potential impacts of climate change, harvesting, oil and gas, and fire. Forestry Chronicle 84(3):349-368. 146. Zar, J.H. 1996. Biostatistical analysis. Prentice - Hall, Upper Saddle River NY. 662 p 147. Zhang, W.H. and D.R. Montgomery. 1994. Digital Elevation Model Grid Size, Landscape Representation, and Hydrologic Simulations. Water Resources Research 30(4):1019-1028. 92 5. MULTIMODELING ASSESSMENT OF THE IMPACT OF PETROLEUM DRILLING AND FOREST HARVESTING ON WILDLIFE HABITAT1 5.1 Introduction The increasing demand for petroleum products that has characterized the first decade of the third millennium has led to the profound transformation of landscapes where oil or gas drilling takes place (Yamasaki et al. 2008). When combined with other natural or anthropogenic disturbances, such as fire, insect outbreaks and forest harvesting, petroleum activities can significantly change the environment (Government of Canada 1992). The impact of different developments on the landscape was and is the focus of a large number of scientific investigations. A query of the ISI database using “forest harvesting” and “oil” as keywords reveals that more than 30 articles have been published on the topic over the last 10 years. The large number of articles has arisen partly because the computational challenges associated with the forecast of complex environmental systems have been overcome by technological advances since 1980. While technological advances have largely solved the computation problems, they have not addressed the theoretical foundation of environmental forecasting. One characteristic of the software used in large-scale environmental investigations is the mono-analytical mathematical framework (i.e., most software is restricted to a single mathematical approach to simulating environmental dynamics), such as DHSVM (Wigmosta et al. 1994), FORECAST (Kimmins et al. 1999), ESRI Hydro Data Model (Environmental Systems Research Institute 2001), ALCES (Stelfox 1999) or SORTIE (Menard et al. 2002). As a result, 97% of the papers identified in the ISI database provide results that rely on only one quantitative framework. The use of a single but complex investigative tool to forecast environmental dynamics leads to results with null probability of occurrence as: Probability of occurrence = 1/ (the number of equally likely results) ≤ (the smallest number of equally likely results (5.1.) supplied by each distinct entity used in the quantification)-n≤ 0)11(lim =+ −∞→ n n 1 A version of this chapter has been submitted for publication. Strimbu, B.M. and Innes, J.L. Multimodeling assessment of the impact of petroleum drilling and forest harvesting on wildlife habitat in northeastern British Columbia. 93 where n is the number of distinct entities used in the forecast and the two 1s correspond to the smallest (respectively largest) value of the confidence interval supplied by the equations used in a model entity. In generating this, I have discounted software that does not provide confidence intervals for each entity used in forecast; as such software is based on the unrealistic assumption that the future is perfectly predictable and known. The aim of the study presented here is to analyze the possible and cumulated impact of specific human activities (i.e., petroleum drilling and forest harvesting) on aspects of the environment (e.g., wildlife habitat) by providing confidence in the results. To generate confidence intervals, I advocate a fundamental change in the methods used in the past to forecast the environment, changes derived from the methodology developed by Boltzman (1995) for thermodynamic systems. 5.2 Methods I used three areas from northeastern British Columbia to determine the impact of petroleum developments and forest harvesting on wildlife habitat (Fig 5.1). The dominant tree species in northeastern British Columbia are white spruce (Picea glauca), black spruce (Picea mariana), trembling aspen (Populus tremuloides), lodgepole pine (Pinus contorta), balsam poplar (Populus balsamifera), Engelmann spruce (Picea engelmannii) and subalpine fir (Abies lasiocarpa) (Meidinger and Pojar 1991). The biogeoclimatic zones associated with these species are the Boreal White-Black Spruce zone, located on the western side of the plains of the Western Canadian Sedimentary Basin, and the Spruce-Willow-Birch zone and Engelmann Spruce-Subalpine Fir zone found on the eastern foothills of the Rocky Mountains at elevations above 1300 m (Ministry of Sustainable Resource Management 2004). 94 Figure 5.1 Study areas The boundaries of the three areas were selected by the Treaty 8 Tribal Council and reflect differences in the intensity of disturbance produced by petroleum drilling activities in the landscape: low (Moberly area), medium (Fort Nelson area) and high (Doig area). The impact of forest harvesting and petroleum drilling was assessed for two valued ecosystem components (VEC), in the sense of Beanlands and Duinker (1983), which were also selected by the Treaty 8 Tribal Council: moose (Alces alces) and American marten (Martes americana). The two human activities and the two VECs ensure that the analytical framework used to assess the impact of landscape disturbance on different environmental attributes can be generalized, as there is no difference between multidimensional investigations (i.e., n activities and p VECs) and two- dimensional investigations when Euclidean spaces represent the investigations (Gilbert 1976). The data that were used to evaluate the responses of the two VECs to human developments were supplied by the British Columbia (BC) Ministry of Sustainable Resource Management (Ministry of Sustainable Resource Management 2004) for forest and wildlife, and by the Oil and 95 Gas Commission (Oil and Gas Commission 2005) and Energy Information Administration (Energy Information Administration 2005) for petroleum wells. The analysis was performed using SAS 9.1 (SAS Institute Inc. 2004) and the input-data were prepared with ArcGIS 9.3 (Environmental Systems Research Institute 2008). To determine the confidence surrounding the results, understood here as the probability of rejecting the null hypothesis, I developed a set of independent and equally likely futures. Each future represented a possible evolution of the environment under potential human developments or natural events. These futures symbolized the trajectories of the different particles used by Boltzman (1995) to describe a fluid. Similar to Boltzman’s trajectories, each future represented a sample unit that expressed the dynamics of the environment. However, rather than providing the average values for each attribute of interest, which was the goal of Boltzman’s approach, I focused on delineating significantly different periods that could exist during the evolution of the environment (Fig.5.2). Therefore, to identify the impact of the human activities on the environment, I determined the commonalities or differences that existed throughout the planning period for the entire set of futures. In comparison to Boltzman’s approach I changed the computation objectives (i.e., from averages to periods) as the particles are real entities while the set of futures is non-existent and an average is irrelevant being based on imaginary inputs. The validity of the methodology required that the set of futures be independent and equally likely, a condition of the Boltzman approach. The independence and equal likelihood conditions recommended the description of each future using different (ensuring independence) and possible (ensuring the likelihood) forecasting methodologies, and led to the seemingly parallel structure from Fig 5.2. 96 Figure 5.2 Development of the set of futures To ensure reliable inferences under the structure proposed in Fig 5.2, I generated at least 10 futures/area, the minimum number recommended by Tran (1997). Each future was spatio- temporally explicit and was obtained by combining the development of forestry [using heuristic techniques (Bettinger et al. 2002)] and oil drilling [using a combination of autoregressive (AR) models (Walls 1992) with indicator kriging (Journel 1983)]. The futures built using the above approach had to meet the regulatory requirements imposed on forest harvesting and oil drilling, but did not necessarily ensure that the evolution of the environment was on a desirable path. While forestry and oil drilling were controlled by regulatory constraints, other attributes quantifying the environment and impacted by human development (such as wildlife habitat) could exceed the appropriate thresholds, making the respective future improbable. Identification of significant periods within the planning period was therefore based on the futures fulfilling not only the requirements of each individual human activity but also on constraints associated with other attributes describing the environment (e.g., wildlife habitat suitability). 97 5.2.1 Development of independent and equally likely futures To develop the independent and equally likely futures, I adopted a hierarchical framework that reflects real developments: forest harvesting could only occur on sites not previously harvested by the oil drilling activities. I operated within the multimodel framework from Fig. 5.2, as in British Columbia mineral rights have priority over surface rights (Government of British Columbia 2004). Oil drilling could therefore not be impeded by any ownership rights or forest management activities, except when specifically regulated by the law (such as archeological sites, sites of endangered species or nesting sites of waterfowl birds). Consequently, the possible development of oil drilling was developed first, with forest harvesting being constrained to incorporate land already cleared for well pads. 5.2.1.1 Generation of futures for oil and gas drilling The spatio-temporal development of oil drilling was represented by first forecasting the evolution in the number of wells and then distributing the wells across the landscape. This approach has been advocated by the large number of models developed to predict the evolution of the oil well count (Iledare and Pulsipher 1997; Ringlund et al. 2008; Walls 1994) and the reduced availability of models to describe directly the spatial-temporal dynamics of the well drilling process (Rao 2000). To accommodate the yearly time step of forest planning, the temporal development of the oil and gas drilling was also represented with an annual time-step: using either the number of annual wells drilled (Cuddington and Moss 2001; Walls 1994) or the total number of wells existing on the landscape (Walls 1994; Walls 1992). I developed equations describing the temporal evolution of the oil and gas wells count for the entire region as the three study areas could be considered as a sample of the entire northeastern British Columbia. Furthermore, the larger amount of information available for the whole of northeastern British Columbia led to more accurate models (Neter et al. 1996). Two types of models were elaborated: one based on autoregressive equations (Brockwell and Davis 1996) and one based on generalized models (Neter et al. 1996). Autoregressive equations are the most common tools describing well count development (Holland 2008; Lanza et al. 2005; Rao 2000; Ringlund et al. 2008; Walls 1994), and I therefore developed three different models to predict the future petroleum drilling using the autoregressive approach. The autoregressive approach was based on the equation proposed by Balke and Gordon (1986) to which a set of covariates were added to the auto-regressed variable: 98 ttttt sXLCyLB νμ +++= )()( (5.2) where yt and Xt denote the well activity and the matrix of the covariate variables in period t, B(L) and C(L) are the two lag-polynomials in the lag operator L which are given by: ∑ ∑ = = = −= q i i i u i i i LcLC LbLB 0 1 )( 1)( μt is the trend, st is the seasonality and νt is the auto-correlated error term: ∑ = − +−= m i titit 1 ενϕν with tNIIDt ∀),0(~ 2σε The covariates considered in the three autoregressive models enhanced the approach of Ringlund and others (2008) and Lanza and others (2005), which considered only one covariate (i.e. the average oil price), by assessing the impact of three prices (i.e., average, minimum and maximum) on the oil and gas drilling. The method proposed by Cleveland (1972) for analysis of the autocorrelation, inverse autocorrelation and partial autocorrelation functions was used to identify the order of the autoregressive model as well as the presence of a non-linear trend in the model describing the dynamic of the number of wells. The augmented Dickey-Fuller (DF) test evaluated the slope and the intercept of the model describing the well drilling dynamics and a Bartlett’s Kolmogorov – Smirnov test was used to assess whether the difference between two consecutive years was white noise. The Smallest Canonical (SCAN) correlation and Extended Sample Autocorrelation Function (ESACF) methods (Box et al. 1994; Choi 1992; Tsay and Tiao 1985) were used to identify whether an ARMA process can represent oil and gas well drilling. Contingent on a stationary and invertible ARMA process, the Minimum Information Criterion (MINIC) method was employed to identify the order of the ARMA process while a Bayesian information criterion (BIC) was used to reveal the autoregressive order of the errors. Lindsey (1980) and Kingman (1996) used generalized models as an alternative to time series models to describe the counting processes (e.g., number of oil and gas wells). Following the 99 recommendation of Cameron and Trivedi (1998), the drilling process was represented using two generalized models: Poisson and negative binomial (Eq. 5.3): ttt XLCy ε+= ))(exp( (5.3) The five models (i.e., three autoregressive and two generalized models) were considered as being correct when the residuals were white noise, assessed using a Bartlett’s Kolmogorov- Smirnov test, and stationary, assessed using a Philps-Peron test (Phillips and Perron 1988). Once the temporal modeling for the entire northeastern British Columbia was finished I tested the validity of the five equations on each of the three study areas. The wells predicted to be drilled in each year were spatially distributed using indicator kriging (IK). Indicator kriging (Journel 1983) operates on the binary set {0,1} that can be used to represent the absence, 0, or the presence, 1, at a given location of a well. The variable representing the presence or absence of a well follows a Bernoulli distribution; the classic variogram was therefore replaced with a robust indicator variogram, as proposed by Hawkins and Cressie (1984): ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ −+= ∑ 5.0 )( 1 |),(),(| )( 1 )(/988.0914.0 1)( jj hN ii yxIyxIhNhN hγ (5.4) where h is the distance between wells, xi, yi are the coordinates of a well, I(xi, yi) is the indicator function at location (xi, yi), and N(h) is the number of points quantifying the presence or absence of a well for the h distance class. The results of the indicator kriging were real numbers not integers; I therefore interpreted the indicator kriging as quantifying the probability of a well being drilled at a specific location. In any estimate of the probability of drilling using IK, the input dataset should contain at least three attributes: x and y (which gives the position) and 0 and 1 [that identify the well’s presence or absence at the (x,y) location]. Data supplied by the BC Oil and Gas Commission (Oil and Gas Commission 2005) contain details of the wells present in the study areas. The dataset is suitable for the assignation of the presence value (1), but does not indicate the absence of drilling (0). I therefore sought to develop a method to identify the places where drilling cannot occur. A possible procedure would be to follow the intuitive approach of assigning the value 0 to 100 areas where no drilling can occur (e.g., inside populated places, archeological sites, streams or lakes). However, this approach failed to identify the 0 in areas where drilling can occur (e.g., pastures or forest with no harvesting constraints) and led to clusters of 0s that associated large lag distances to non-negative values, of little use in the indicator variogram elaboration. The inability to delineate 0 from 1 using such an intuitive approach was based on the violation of two conditions that describe the development of drilling: 1. the true variogram would not be known until the drilling will stop, as the maximum number of wells in the area is reached; 2. for any given year the position of wells that would be drilled the following years is unknown. A consequence of the first condition is that the variograms built using data before the maximum number of wells is reached will converge to the true variogram (Kolmogorov and Fomin 1999). The second condition indicated that splitting the dataset containing the well information into two subsets, according to the year of drilling [i.e., wells drilled before a selected year of drilling would have the position known (past) while the wells drilled after the selected year would have the position unknown (future)] would create a set of variograms that would converge to the true variogram (the variogram for the year when drilling will cease). The dataset provided by the Oil and Gas Commission (Oil and Gas Commission 2005) was therefore divided into two, such that for wells drilled before a specific year the I(x,y)=1 while for wells drilled after the respective year the I(x,y)=0. The three study areas had different geomorphologies that could lead to different distributions of wells on the landscape. I therefore developed individual variograms for each of the three areas. The variogram fitting was assessed using Q1 and Q2 statistics, as suggested by Kitanidis (1997). Besides the variogram, the accuracy of the IK estimates also depends on the directional symmetry of the spatial correlation, which was tested using the isotropy test proposed by Lu and Zimmerman (1994). The IK developed using the above procedure computes the probability of oil and gas drilling at specific locations. To identify the locations relevant to the drilling process I constructed a grid that would supply the coordinates where the IK would be estimated. The grid was based on the quarter section system used by the Government of British Columbia in landscape planning (Ministry of Sustainable Resource Management 2004). Each quarter section was divided into four quadrants, as no more than four wells can be drilled on a quarter-section (Agricultural Land 101 Commission and Oil and Gas Commission 2002). For each well, the Cartesian coordinates relative to the quarter-section in which it was located were determined. Using the Cartesian coordinates, all the wells were plotted on one quarter-section and the probability of occurrence of a new well within a quadrant was considered as being proportional to the number of wells within the respective quadrants. Of the four values, only the largest two were selected as regulations constraining drilling activities give unobstructed authority to the Oil and Gas Commission for a maximum of two wells per quarter section (Government of British Columbia 2004). I assumed that oil and gas development in BC would mirror that in Alberta, the Canadian province with the most experience in oil and gas exploration (Information Services AEUB 2005), and that the density of oil and gas wells on the landscape would reach the maximum capacity allowed by the regulations (Anonymous 2004). The grid used for computing the IK has as vertices the center of the two quadrants with the largest number of wells within the quarter- section grid. Finally, the moment of drilling at a certain location was determined by ranking the probability of drilling (i.e., ranking the results of indicator kriging), and, within the order set, assigning to each year the number of wells supplied by one of the five temporal models. This approach was based on two properties: the exact interpolation of IK (Cressie 1993) and the strict monotonicity of the cumulative probability function (Grimmett and Stirzaker 2002), properties justifying the assumption that larger probabilities are associated with earlier drilling. 5.2.1.2 Generation of futures for forest harvesting In northeastern British Columbia, the annual allowable cut (AAC) is determined using heuristic algorithms constrained to meet different spatial and temporal conditions (Jukes et al. 2004). In addition to the requirements of the sustainable forest management plan used for the area (Jukes et al. 2004) I accommodated the hierarchy of the two industries by constraining the scheduling of the forest harvesting to include the surfaces harvested by the petroleum drilling (i.e., well pads). I generated four independent and equally likely future forest landscapes by using two possible heuristic algorithms [i.e., simulated annealing (SA) (Metropolis et al. 1953) and an adjusted first fit decreasing (FFD) algorithm (Johnson 1973; Strimbu et al. 2008) ] and two harvesting ages (i.e., a preset harvesting age (Bettinger et al. 2002) and age at max(MAI) (Davis et al. 2001)]. I selected SA as it is the main algorithm used for forest scheduling in northeastern British Columbia as well as being one of the heuristic algorithms with the simplest implementation (Lockwood and Moore 1993). Furthermore, SA supplies results close to the linear programming solution describing the forest planning problem, which is considered an optimal solution (Crowe and Nelson 2003; Bettinger et al. 2002). To select the parameters 102 required by the SA (i.e., initial temperature (Ti), freezing temperature (Tf) and the rate of annealing (ra)) I followed the recommendations of Boston and Bettinger (1999), Bettinger et al. (2002) and Crowe and Nelson (2003). To obtain a solution I allowed SA to run on average 60 minutes, adhering to the results of Crowe and Nelson (2003), which suggested that SA could supply an efficient solution in less than one hour. I selected the rates of annealing at 0.7, 0.75, 0.8, 0.85, 0.9, 0.95 and 0.99, the initial temperatures at 1000, 10000 and 100000, and the freezing temperatures at 1, 50 and 100, similar to the values used by Bettinger et al. (2002) and Crowe and Nelson (2003). I adjusted the FFD algorithm in the sense that the FFD algorithm was applied to a changing bin size, downwards, starting from a maximal AAC. The maximum AAC was determined using the perfect bin-packing theorem (Coffman et al. 2000) and represented the upper boundary of the AAC, i.e. the maximum volume that could be harvested annually in the absence of spatial constraints. Strimbu et al. (2008) have shown that the aspatial upper AAC boundary is the sum of the maximum mean annual increment of the stands. The computation of the upper AAC boundary did not consider spatial constraints, and the inclusion of green-up constraints in the scheduling process could therefore lead to values smaller than the upper AAC boundary. To meet the regulations, I determined the AAC by decreasing the upper AAC boundary (i.e., sum of the maximum mean annual increment of all stands) in increments of 5%, until all the forest planning constraints were met. The constraints (Table 5.1) governing the magnitude of the volume that could be harvested annually were based on the Sustainable Forest Management Plan developed for the Fort St. John Pilot Project (Jukes et al. 2004). Table 5.1 The harvest scheduling constraints No. Constraint Value 1 The maximum yearly variation of the AAC from the mean AAC of the planning period 10% 2 The stands are harvested using one silvicultural system. Clear-cutting 3 The minimum opening size created by harvesting 5 ha 4 The maximum opening size created by harvesting 60 ha 9 The length of time until a stand is considered fully stocked after harvesting (the greenup/adjacency constraint) 5 years 10 The planning horizon 100 years The maximum volume to be harvested annually at a forest level occurs when the harvesting coincides with the maximum MAI of each stand. The maximum AAC should therefore ensure 103 that each stand is harvested at the age when MAI peaks (Davis et al. 2001). The harvesting guidelines for northeastern British Columbia recommend merchantability standards of 17.5 cm for spruce and sub-alpine fir and 12.5 cm for all other species (Pedersen 2003), and these standards were used to determine the culmination age of the MAI (Appendix C). In addition to the future forested landscapes developed using an algorithm aiming to achieve the maximization of AAC (i.e., adjusted FFD and harvesting at the peak MAI), I developed futures using a preset harvesting age, which is the method currently being used to schedule forest harvesting in the area. For the three study areas the minimum harvesting age was established 120 years, following the recommendations of Jukes et al. (2004). Harvesting above a preset age could result in stands being harvested 70 years before the culmination of MAI (e.g., a stand of white spruce with SI ≤10) leading to local maximization of the AAC. The constraints that could play a significant role in determining the AAC that were not incorporated in the mathematical model were seral stage distribution and growing stock. The constraints regulating the seral stage distribution were not represented mathematically as either the harvest occurred at ages larger than 120 years (for the future considering that the stands are harvested at a preset minimum age) or more than 30% of the THLB had the harvest age larger than 120 years (Appendix C) (for the futures considering that the stands are harvested when MAI peaks). When adjusted FFD is used to schedule the harvests, then the perfect bin- packing theorem (Coffman et al. 2000) conditions were met; consequently, the growing stock constraint was fulfilled, as AAC was smaller or equal to MVHA (in this case MVHA = annual growth of the forest estate where harvest can occur). When the harvests are scheduled using SA, the growing stock could decrease to undesired levels, therefore the implementation of a set of constraints addressing the reduction of the growing stock became necessary. However, if the AAC supplied by SA is smaller than MVHA and the perfect bin-packing theorem condition are fulfilled then the magnitude of the growing stock would not decline to unacceptable levels. The forest scheduling computations were performed using Delphi 7.0 (Borland Software Corporation 2001) and the runs were executed on a PC with a 3 GHz Pentium 4 processor and Windows XP operating system. The forest yield was determined using TIPSY 3.1 (Ministry of Forests 2005a) for managed stands and WinVDYP (Ministry of Forests 2005b) for unmanaged stands. 104 5.2.2 Assessment of the cumulative impact of oil drilling and forest harvesting 5.2.2.1 Quantification of moose and American marten habitat The evolution of the environment was assessed using habitat suitability indexes (HSI) that measure the quality and quantity of the habitat of the two VECs selected (i.e., moose and American marten). The HSI models do not accurately represent the response of the species to the environmental changes but only the habitat quality and availability, a necessary but not sufficient condition. However, several studies identified a significant correlation between population dynamics and different indices describing habitat (Hirzel et al. 2006; Larson et al. 2004; Ottaviani et al. 2004; Roloff and Haufler 1997; Roloff and Kernohan 1999), recommending HSI as a suitable measure for the trends in moose and American marten populations. The limiting season for the two species is winter, and I therefore used only the HSI quantifying habitat suitability during the cold season (Romito et al. 1999; Takats et al. 1999). Both HSI models relied on the assumption that proximity to human settlements or human activities do not affect the species behavior and there is unrestricted access by the animals to water and mineral resources in the area that supplies food and cover. According to Romito et al (1999), the moose winter habitat in the region is influenced mainly by the amount of shrubs and deciduous saplings [a stands is quantified as being 100% suitable if the shrubs and deciduous saplings cover 25% of the stand]. The requirements of American marten for the winter habitat are more nuanced than moose, and depended on the simultaneous fulfillment of several conditions: presence of coniferous (a stand with more than 15% pine, spruce or fir then is 100% suitable), relatively close canopy (between 30% and 70% a stands is assumed to be 100% suitable), and presence of tall trees (stands with height larger than 15 m are 100% suitable). Two statistics were used to assess the cumulative impact of oil and gas drilling and forest harvesting: one describing the average change of the HSI at the landscape level and one measuring the availability of superior habitats (i.e., HSI>0.5) across the landscape. The 0.5 threshold for the HSI was selected based on the results of Jager et al. (2006), Kroll and Haufler (2006), Liu et al. (2005) and Bailey et al. (2002). The average HSI was selected as indicating the habitat dynamics at the landscape level as well as whether or not human activities significantly changed the landscape from the habitat perspective (2002). The superior habitats complemented the average HSI, as HSI models could misrepresent reality by steep declines in species occurrence; declines predicted when poor habitats are not excluded from the modeled landscape (Jager et al. 2006). 105 5.2.2.2 Analysis of the futures induced by oil and gas drilling and forest harvesting The hierarchy outlined in Fig. 5.2 does not ensure that activities situated on the lowest hierarchical level will not supply futures having the two HSI statistics violating regulatory thresholds. To ensure that all futures met the regulatory requirements, the average HSI and the area with HSI>0.5 of each future and year of the planning period were compared with a threshold HSI. I considered that the threshold HSI was reached when the minimum amount of suitable habitat would affect the population dynamics of moose and American marten. Based on the results of Peterson (1997) and Fryxell et al. (1999), I considered that for both moose and American marten a 15% reduction of the two HSI statistics would ensure that the population change would be within the natural range of variation. The reference value was selected the statistic of the first year of each future, as the current conservation status for both moose and American marten (i.e., yellow listed by BC Conservation data center and G5 by NatureServe) is desirable to be maintained or enhanced during the future. The two statistics used to quantify the evolution of the environment operate on two scales (i.e., the average HSI is between 0 and 1 and the area with HSI>0.5 is larger than 100,000). Rencher (2002) and Neter et al. (1996) mention that within a repeated measures framework, the difference in the orders of magnitude could affect the results. To include values operating on a 105 scale with the values operating on a 100 scale in the same analysis the values representing the size of the area with HSI>0.5 (i.e., between 100 000 and 500 000) were transformed to values compatible with the average HSI (i.e., between 0 and 1). The transformation to subunit values was performed using the linear relationship: )]min()/[max()]min([ 5.05.05.05.0 >>>> −− HSIHSIHSIHSI AreaAreaAreaArea (5.5) where min and max are the minimum, respective maximum, operators applied to the 100 HSI values within each future. I performed a univariate repeated measures analysis of the set of futures to identify the impacts of human developments on the habitat of moose and American marten. I used a Helmert transformation (1990) to reveal the existence of significantly different periods, while the differences between the futures were assessed using univariate two-sample profiles, as advocated by Rencher (2002). When the assumptions required for valid repeated measures analysis were violated [i.e., the Huynh-Feldt condition, which was tested using the sphericity test 106 proposed by Anderson (1963)], the significance level of the F-test was adjusted using the modifications proposed by Greenhouse and Geisser (1996) and Huyhn-Feldt (1970). The differences between the futures and the moments when the habitat suitability could change significantly were assessed using the Scheffe test, as suggested by Zar (1996). To enhance the one-dimensional perspective supplied by the univariate approach, I used multivariate repeated measures to reveal the existence of the moments when significant changes could occur in the wildlife habitat during the planning period. This technique also enabled the assessment of the impact of the models in the delineation of the homogeneous periods (i.e., the intervals between the moments when significant changes in the wildlife habitat could occur). I used Wilks’ likelihood ratio, Pillai’s trace and Lawley-Hotelling tests for the interpretation [‘Roy’s root’ (another common multivariate test used in significance testing) is an extension of ‘Pillai’s trace’ tests (Rencher 2002) and therefore was not used]. To ensure the consistency of the comparison between the univariate and multivariate repeated measures, the homogeneous periods were identified using a Helmert transformation in combination with the multivariate two-sample profiles (Rencher 2002); this also enabled the evaluation of the impact of the models on the analysis. The multivariate repeated measures analysis was supplemented by a hierarchical cluster analysis that provided a different multivariate perspective on the possible commonalities in the set of futures. As with the repeated measures analysis, hierarchical cluster analysis was used to separate the futures forecasting period (i.e., 100 years) into distinct sub-periods. The similarities within the set of futures and the moments when significant changes could occur were identified using centroid and Ward’s minimum variance clustering methods, as recommended by Rencher (2002). The distance between two clusters, DAB, was measured using the increase in the sum of squares for Ward’s method and the distance between the mean vectors of two clusters for the centroid method (Rencher 2002): ∑∑∑ ∈∈∈ −−−−−= BAAB Ci Bi Ci Ai Ci ABiAB xxxxxxD 222 |||||||||||| (Ward’s method) (5.6) ∑ −= 2|||| BAAB xxD (centroid method) where xi is the vector representing the future i, when the common periods are identified, or the year i, when the difference between the models is assessed Ax is the mean vector of the cluster CA, CB or CAB, which are in the relationship BAAB CCC ∪= )()(|||| 2 Ai T AiAi xxxxxx −×−=− 107 The representative clusters were delineated based on the formal procedure proposed by Mojena (1977) and improved by Milligan and Cooper (1985), which selects the number of clusters based on the change in the distance between two adjacent levels of grouping at which njsj ...,,2,125.1 =×+> δδδ (5.7) where: δ1, δ2, ... , δn are the distances values for stages with n, n-1, …, 1 clusters, δ and sδ are the mean and the standard deviation of the distances between two adjacent groups, and 1.25 is a constant proposed by Milligan and Cooper (1985). A more detailed investigation of the changes in habitat suitability was supplied by principal component analysis (PCA), which enhanced the results supplied by the hierarchical cluster analysis by allowing a formal analysis of the similarities within the set of futures. Besides revealing the similarities, I used PCA to identify the moments when the wildlife habitat underwent significant changes. The eigenvectors with a level of variance larger than 5% were selected for interpretation based on Noble et al. (2004), while the attributes with loadings larger than 0.32 were considered significant, as suggested by Tabachnick and Fidell (2001). Strimbu and others (2009) have shown that in landscape investigations the distribution of the residuals can affect the significance of the results, and so when the data were found not to be multi- normal, the PCA provided only an indication of the similarities within the set of futures. The PCA identified the similar futures and the moments when wildlife habitat could change significantly by considering that the study area and the forecasting methodology had no influence on the analysis. To check this assumption, I used the methodology proposed by Mackenzie and others (1999), namely canonical discriminant analysis (CDA). The delineation of the periods was performed based only on the harvesting algorithms (i.e., SA and FFD) and not on the drilling models, as the oil and gas disturbance is incorporated in the AAC computation. I identified the impact of the years from the planning period on the harvesting algorithm by selecting for interpenetration the years with correlation between the years and the canonical variates larger than 0.32, as well as the years with standardized canonical coefficients larger than half the maximum standardized coefficient, as advocated by Rencher (2002). The results supplied by CDA are relatively robust to the violation of the assumptions, as CDA relies on relatively weak assumptions [i.e., a multi-normal data distribution is only required within each categorical attribute (Rao 1973) and, for descriptive investigations, the homogeneity of the covariance matrix is relaxed]. 108 Parametric repeated measures, PCA and CDA rely on the assumption of data normality (Diggle et al. 2002; Hardle and Simar 2003). To test this assumption I used two univariate tests (the Shapiro-Wilks and Kolmogorov-Smirnov tests) and three multivariate tests [the Mardia kurtosis and skewness tests (Mardia 1970) and the Henze-Zirkler test (Henze and Zirkler 1990)]. The two univariate normality tests covered the main mathematical approaches associated with distributional testing: Shapiro-Wilks, which is a regression based test, and Kolmogorov-Smirnov, which is based on the empirical distribution function (Thode 2002). The three multivariate tests cover two of the main procedures used to assess the multivariate normality (Thode 2002): the direct – data assessment (Mardia kurtosis and skewness tests) and the empirical characteristic function (Henze and Zirkel test). For this study, the variance homogeneity assumption, required by the Gauss theorem, was relaxed because I did not use the results for prediction, but to detect and describe the patterns (e.g., similarities, differences or periods existent within the set of futures). Besides normality, the outliers were of particular interest because the computations were performed using the least square method, which is sensitive to outliers (Park and Lee 2004; Rocke and Woodruff 1996; Schwager and Margolin 1982) I used the Schwager – Margolin (1982) approach to identify the outliers as it detects several outliers simultaneously (Rencher 2002) and uses the maximum likelihood method, which seems to be more robust than the least square method (Islam and Tiku 2004). 5.3 Results The augmented Dickey-Fuller (DF) test was the first step in the development of the models describing the evolution of the number of oil wells and indicated that the annual well count series may have a difference – stationary process, as p =0.999. The DF test was confirmed by Bartlet’s Kolmogorov Smirnov test, which showed that the difference between the numbers of wells in two consecutive years could be white noise (p=0.74). However, the investigation of the autocorrelation functions revealed that a first order autoregressive model could be used to represent the number of wells drilled annually, as the autocorrelation function for a one year lag was significantly larger than 0 (p<0.01), while the inverse and partial autocorrelations were insignificant, regardless of the lag size (p>0.3). For the total number of wells, the autocorrelation functions indicated the presence of a non-linear trend and the absence of a relationship for any lag larger than one year (p>0.24). Therefore, regardless of the series used to describe the evolution of the number of wells (i.e., annual or total), there was no evidence of a relationship 109 between drilling activities more than two years apart, in agreement with the results of Lanza and others (2005). The SCAN correlation and ESACF methods supplied different results for the number of wells drilled annually, indicating that an ARMA process cannot be used to model the yearly dynamics of well drilling. The total number of wells could be represented by an ARMA (2, 2) process, as revealed by SCAN and ESCAF. However, MINIC did not confirm the SCAN and ESCAF results, as one of the largest BIC was associated with the ARMA (2, 2) process. Consequently, the ARMA process was not used to model the evolution of the number of wells, whether annual or total. The two types of models (i.e., autoregressive equations and generalized models) led to the development of five possible patterns in the evolution of the number of petroleum wells (Table 5.2). All the models fulfilled the main time series assumptions [i.e., significance (p<0.001), and residuals were organized as white noise (p>0.12) and were stationary (p>0.17)]. The residuals were normally distributed for the models based on autoregressive equations (p>0.5), and non- normal for the Poisson and Negative binomial models. The residuals analysis indicated that the functions describing the number of wells are nonlinear (parabolic or exponential), whether using the annual or total number of wells. The nonlinear evolution of the number of wells was reinforced by the presence of a model with a periodic fluctuation, increasing amplitude and decreasing period (model three in Table 5.2) and reconfirmed the inappropriateness of the ARMA process for the northeastern BC well data. Although the Poisson model had a deviance greater than 1, based on the recommendation that models with deviance smaller than 10 are not necessarily incorrect (Neter et al. 1996), I concluded that the results supplied by both generalized models should be selected for the interpretation. Regardless of the model type (i.e., autoregressive or generalized model), the annual oil prices seemed to have no impact on the drilling effort, confirming the results of Ringlund et al. (2008) and Lanza et al. (2005). The fulfillment of the assumptions and the agreement with other studies investigating the evolution of the number of wells (Iledare and Pulsipher 1997; Lanza et al. 2005; Rao 2000; Ringlund et al. 2008; Walls 1994) suggests that the models are correct and do not have different likelihoods. Therefore, it can be concluded that the futures built using the five models can occur and are equally likely. 110 Table 5.2 Well drilling Temporal dynamics of the number of wells Equation Pr>F Normal. White noise Stationarity Zero mean Trend 1 y t = 1.255 x y t-1 + εt 0.001 0.09 0.69 0.99 0.99 2 Total t = 1.01 x totalt-1 + 1.1 x y t-1 + εt 0.001 0.05 0.78 0.17 0.17 3 y t =0.54 x totalt-1 – 1.02 x (year- 1948)2 – 0.44 x year x f(t)+ εt 0.001 0.13 0.12 0.51 0.51 4 y t=exp(3.62+5.7 x 10-4 x totalt-1 – 1.6 x 10-8 x total2t-1)+ εt 0.47* 0.001 0.28 0.80 0.80 5 y t = exp(3.69+5.3 x 10-4 x totalt-1 – 1.39 x 10-8 x total2t-1)+ εt 0.75* 0.001 0.82 0.98 0.98 Spatial dynamics of the wells Study area Variogram type Equation P(|Q1||=0) P(|Q2|=1) P(vario.is isotropic) Moberly Spherical χ(h)=0.05×(1.5×h/4000- 0.5×(h/4000) 3) 0.17 0.23 0.061 Doig Spherical χ(h)=0.063×(1.5×h/1500 0-0.5×(h/15000) 3) 0.22 0.16 0.25 Fort Nelson Linear χ(h)=2×10-6×h 0.19 0.13 0.14 Total t – total number of wells drilled at year t; y t – number of wells drilled in year t; f(t) – is the periodic function sin(π(year-1948)/period), where period=12 before 1988 and 8 after 1987; * - for Poisson and negative binomial models the p-value refers to chi-square goodness of fit for Poisson, respectively negative binomial distribution; εt ~ WN (0, σε); χ(h)- variogram; h is the distance between wells; ∑ =−= n k kn Q 2 1 1 1 ξ and ∑ =−= n k k n Q 2 2 2 1 1 ξ (5.8) while ζk is the kth orthonormalized residual obtained using n points. The set of five equations were refitted for each study area and compared with the equations determined for the entire northeastern British Columbia region. The pair-wise comparison of the 111 parameters was insignificant (p>0.05), except for the Moberly and Doig areas, where the model including a periodic term and the Poisson model were insignificant, as was the variable representing the squared total number of wells in the Negative Binomial model. However, as the sign and magnitude of the insignificant variables were similar to the ones developed for the entire region, the evolution of the drilling effort in the three study areas was considered to be adequately represented by the equations developed for the entire region. To this end, the models based on autoregressive equations will reach the maximum number of wells allowed by the regulations (i.e., 20,000 in the Doig area, 19,000 in the Fort Nelson area and 8400 in the Moberly area) while the generalized models will not reach that maximum for more than 100 years. The generalized models have different peak moments, depending on the existing drilling intensity: approximate 2025 in the Doig area (the area with the largest number of wells), 2045 in the Fort Nelson area and 2050 in the Moberly area (the area with the lowest number of wells). In northeastern British Columbia, it seems that the geology does not play a significant role in determining the location of a well within the regions known to have oil reserves, as of the 13,000 wells drilled since 1992 (i.e., 72% of the total number of wells existing in September 2005), less than 10% were above newly-discovered reservoirs and only 48% were above known reservoirs. As a result, detailed geological maps do not provide particularly useful information in identifying the location of the disturbances caused by oil well drilling, as at the landscape level, the positioning of the wells does not follow oil and gas discoveries. Consequently, the location, size and shape of the known oil reservoirs were not incorporated into the IK estimates. The three areas have different robust variograms, spherical for the southern areas (the Moberly and Doig areas) and linear for the northern area (Fort Nelson), suggesting a stationary model for the areas with low and high petroleum development and a nonstationary model for the area with intermediate development (Table 5.2). An indication of the appropriate grid size to be used for kriging prediction was provided by the oil and gas developments from Doig area, where 665 quarter sections have already reached the maximum number of wells recommended by current legislation (i.e., 2 wells / quarter section). For more than 75% of the wells in the Doig area, any two wells were separated by more than 300 m, suggesting a grid size of 400 m, which would delineate the quarter section into quadrants and would ensure the integration in the grid used for land use planning in British Columbia. Within the quarter-section, the wells were distributed following a pattern of 35%-21%-23%-21%, clockwise from the NE quadrant to the NW quadrant. The 35-21-23-21 pattern identified for the Doig area was also observed for the British Columbian part of the Western Canadian Sedimentary Basin. The probability of drilling at the 112 center of the quadrant was therefore assigned according to the pattern identified in the Doig area. For example, if a well was drilled in the quadrant where the probability of drilling was 21%, then the next well drilled within the quarter-section would be in the quadrant with a probability 35%, while if there is no well drilled within the quarter-section the two wells would be drilled in the quadrants with probability 35% and 23%. The maximum volume that can be harvested annually in the absence of the spatial constraints varies according to the timber harvesting land base, with Doig being the largest and Fort Nelson the smallest (Table 5.3). The scheduling algorithm and the harvesting age had a significant impact on the volume harvested annually (Table 5.3), but the FFD was notable in supplying the largest AAC of any algorithm for all tree areas. The green-up adjacency constraint of five years seemed to have little effect on the AAC supplied by the FFD algorithm, especially when the harvesting occurred at the culmination of MAI (i.e., the AAC was less than 5% smaller than the maximum harvestable volume), but had a dramatic effect on SA, especially when the harvesting was performed at a preset age (for Fort Nelson, the AAC was less than half the maximum harvestable volume). The scheduling algorithm seems to be responsible for the magnitude of the AAC, as FFD algorithm supplied results significantly larger than SA regardless of the area (p<0.001), and ranged from 15% for the Fort Nelson area to more than 20% for the Doig and Moberly areas. The increase in the AAC supplied by the FFD algorithm was not necessary reflected in a significant increase in landscape disturbance, as for the Fort Nelson area, a significantly larger AAC (p<0.0001 for AAC=270,000 m3/year supplied by the FFD algorithm compared to 225,000 m3/year supplied by the SA algorithm) did not lead to a significantly larger area of forest being cleared (p>0.2 when the clear-cut area of 950 ha/year determined by FFD algorithm is compared to the clear-cut area of 930 ha/year for SA algorithm). A similar situation was encountered for the Doig and Moberly areas for which the difference between the harvested areas was less than 15%, regardless of which scheduling algorithm was used. 113 Table 5.3 The AAC obtained using the two algorithms (i.e., SA and FFD) and the two harvesting ages (i.e., preset minimum age and Age max(MAI)) Study area Surface [ha] THLB [ha] MVHA [m3] Average AAC [103 m3] / Average Harvested area [ha] First Fit Decreasing Simulated Annealing Agemax(MAI) Age preset Agemax(MAI) Age preset Doig 684138 183057 500000 495 390 355 330 1500 1350 1300 1200 Fort Nelson 641024 91355 280000 270 220 225 120 950 950 930 900 Moberly 410194 140222 315000 290 155 185 205 1050 1000 700 670 The variability (expressed as the coefficient of variation) of the AAC was between 5% for the Fort Nelson area and 30% for the Moberly area. The clearcut area had a coefficient of variation at least 50% larger than the corresponding coefficient of variation for the AAC (e.g., CV AAC – FFD/Age@max(MAI) for Doig = 29% vs. CV cleared area – FFD/Age@max(MAI) for Doig = 53% or CV AAC – SA/Preset age for Moberly = 6% vs. CV cleared area – SA/Preset age for Moberly = 21%). The selection of the harvesting age ensured that, whether or not the stands were harvested at the culmination of MAI (at which more than 30% of the THLB would be harvested at ages >120 years) or at the preset minimum age of 120 years, the seral stage distribution required by the provincial guidelines [stands with age < 40 years should occupy less than 27% of the study area and old seral stages (stands >120 years) should occur on more than a third on the area] (Ministry of Forests 1995) was fulfilled during the entire planning period for all three areas. The clearcuts did not change the forest structure more than 2% in any year of the planning period, irrespective the scheduling algorithms or harvesting age. The proposed multimodel structure from Fig. 5.2 led to 80 futures for each study area (two species x two HSI statistics x two forest scheduling algorithms x two harvesting ages x five drilling models). Among the 240 futures, oil drilling played a significant role only when the environment was quantified using the HSI>0.5 (i.e., 120 futures). When the futures were represented by the average HSI, the five drilling models produced similar results, with differences less than 1%. Consequently, only 24 futures based on the average HSI were used in the analysis. All 144 futures (i.e., 120 futures for HSI>0.5 plus 24 futures for the average HSI) had the two statistics fulfilling the selected thresholds [i.e., the area with HSI>0.5 and average HSI either increased or decreased within the preset limits (less than 15 %,)] indicating that all the futures are possible, therefore suitable for analysis. 114 The analysis of the superior habitats (i.e. stand with HSI>0.5) revealed the important role that drilling activities could play for the first 20 to 40 years. For each of the three areas, harvesting algorithm and harvesting age, the set of futures determined by the five drilling models converged to one future, when the maximum number of wells allowed by regulations was reached. After convergence, the development of stands with HSI>0.5 would depend only on forest harvesting activities. The average HSI was larger for the futures produced by the FFD algorithm regardless of the method used to select the harvesting age (i.e., age at maximum MAI or preset age). The magnitude of the average HSI for moose was significantly sensitive (p<0.01) to the harvesting algorithm (46%-48% for FFD compared with 32%-44% for SA in the Moberly area, 46%-50% compared to 44%-46% in the Fort Nelson area and 64%-68% compared to 49%-52% in the Doig area). American marten habitat seemed little influenced by the algorithm (p>0.2), as the average HSI has similar pair-wise ranges for each study area (34%-38% for the Moberly area, 41%-47% for the Fort Nelson area and 23%-26% for the Doig area). The development of the average habitat suitability was influenced by the method of selecting the harvesting age within the forest scheduling algorithm, as for almost 85% of the comparisons (i.e., 10 out of 12) the trajectories of the average HSI were different, depending on whether the harvesting occurred at MAI or after a preset age. The choice of harvesting age seemed to influence the long-term development of habitat suitability, as for the preset harvest age the average HSI had an increasing trend, regardless of the algorithm, species and study area, while for the harvesting age at the culmination of MAI, half of the futures had a relatively constant HSI and half declined. The decrease in the average HSI was noticed for areas with a reduced timber harvesting land base relative to the size of the entire study area (14% from the entire area for Fort Nelson area and 27% for the Doig area). Nevertheless, the decrease in the average HSI was less than 3% over the 100 year planning period, regardless of which scheduling algorithm and harvesting age was used. While the average HSI decreased in 25% of the futures (six of 24), the area of the stands with HSI>0.5 either remained relatively constant (less than 3% decrease) or increased during the planning period. The FFD algorithm supplied larger areas with superior habitat than the SA algorithm, regardless of the harvesting age. The superior availability of better habitats on landscapes forecasted using FFD compared with the landscapes obtained using SA is probably because the FFD would harvest more than once during the planning period high productive stands (therefore suitable for moose) and would leave the stands with low SI to grow old (therefore suitable for American marten). The areas with HSI>0.5 ranged from 6000 ha for Moberly to 16000 ha for Fort Nelson when the FFD algorithm was used and from 4000 ha for 115 the Moberly area to 10000 ha for the Fort Nelson area when the SA algorithm was used. The availability of superior habitat for the two species developed differently, as the surface of the stands with HSI>0.5 for American marten increasing regardless the forest scheduling algorithm and the harvesting age, while the area of superior habitat for moose had a slow but non- significant decline (e.g., the Fort Nelson area, p>0.06). Drilling had little impact on the occurrence of highly suitable marten habitat, but negatively affected moose habitat, which decreased by as much as 10,000 ha in the Doig area while drilling was active. Drilling could also induce dramatic annual changes in the availability of moose habitat, with significant drops from year to year [up to ten times for the Fort Nelson or Moberly areas]. The examination of the three assumptions with a possible impact on parametric longitudinal investigations (e.g., normality, presence of outliers and sphericity) revealed the partial fulfillment of the analytical requirements. The univariate normality tests agreed with the multivariate tests and indicated that habitat suitability index is not normally distributed (p<0.0001). The Schwager- Margolin test did not identify any year during the futures as an outlier. The same conclusion was reached for the futures themselves, as no future was found to be positioned as an outlier in the space determined by the 144 futures. Besides statistical assumptions, the data transformation could play a significant role in the analysis (Tukey 1977). The linear transformations associated with the multivariate analyses and the change in the magnitude of the values representing the superior habitats (i.e., from values of the order of 106 to the order of 100) would not alter the results of the statistical tests, as isomorphic transformations of a Cartesian coordinate systems (such as translation or rotation) ensure findings validity, irrespective the Cartesian coordinate system (Hungerford 1990). The sphericity test (p<0.001) indicated that the Huynh – Feldt condition was not met by the covariance matrix, so the correction proposed by Greenhouse and Geisser (1959) and Huyn and Feldt (1970) was adopted to interpret the results. The univariate repeated measure revealed that during the next century, the habitat would change significantly (p<0.0001) in all three study areas (p<0.001), regardless of the harvesting algorithm (p≤0.002 for H0: “The habitat would not change during the planning period”, H0 performed for each study area and algorithm). The significant difference between the futures (p=0.02) indicated that the 144 futures covered a wide array of possible habitat developments. The univariate analysis suggested that of the two human activities investigated, forest harvesting was the more important in determining the trend in habitat suitability (p=0.038), while the study area had no impact (p=0.48). The scheduling algorithms had different impacts on the two species, with moose being 116 associated with both SA and FFD (47% of futures with SA and 97% with FFD) while American marten was only related to the FFD (46% of the futures). The Scheffe test divided the planning horizon into three periods: years 1 to 50, years 51 to 93 and years 94 to 100. The HSI within each period had a similar trend: increasing in the first period, constant in the second and decreasing during the last period. The univariate repeated measures indicated that the models determining the futures were significantly different (p<0.0001), with the ten futures for moose developed using SA algorithm and area of better habitats for Fort Nelson being separated from the remaining 134 futures. The 134 futures were further separated in two groups by the SA algorithm: one having 124 futures and one having the 10 futures derived for American marten habitat in the Fort Nelson area. The multivariate repeated measure analysis confirmed the results of the univariate analysis, with all three multivariate tests (i.e. Wilk’s, Pillai’s and Lawley-Hotelling) agreeing that the habitat would change significantly during the planning period (p<0.0001). The planning period was separated into four sub-periods: years 1 to 11, years 12 to 49, years 50 to 75, and after year 75. The main attribute characterizing the four sub-periods was the variability within each sub-period, as the first and third sub-period were homogeneous (the transition from one year to another during the period was insignificant), while the second and fourth sub-period were heterogeneous (the yearly transitions were separated into 13 and 10 distinct groups, respectively). The four sub-periods delineated three distinct moments during the planning period that seemed to occur at years 11, 49 and 74. While the three moments mark the transition from one sub-period to another, there are two obvious separations during the planning period, one occurring at the end of the fourth decade (between years 36 and 39) and one at the end of the century (between years 93 and 99), when all years differ significantly from one another. The multivariate tests also revealed that the futures were significantly different (p<0.0001), but the partition of the futures into groups was difficult as the two-sample profiles identified more than 20 groupings among the 144 futures. The futures developed using the average HSI was so different (only 2 out of 24 futures were not significantly different) that no groups were detected, agreeing with the results obtained for superior habitats. The hierarchical clustering algorithms seemed to have no impact in grouping either the futures or the time periods, as Mojena’s (1977) formal procedure identified the same significant numbers of clusters using the centroid and Ward’s minimum variance methods. I therefore restricted the interpretation to the Ward’s minimum variance method as it supplied a clearer delineation of the clusters than the centroid method. The hierarchical cluster analysis confirmed 117 the results of the multivariate repeated measures analysis in that the planning horizon was separated into four distinct sub-periods (Fig 5.3a). Although two different techniques (cluster analysis and multivariate repeated measures) were used to identify the moments when the habitat would significantly change, the results were very analogous, as cluster analysis revealed that significant habitat transformations could start in years 11, 37 and 71 (compared to years 11, 49 and 74 as indicated by multivariate repeated measures). In contrast to the multivariate repeated measures analysis, which failed to identify similar models, the cluster analysis delineated four classes of futures based on the models used to build the futures (Fig 5.3b). The set of futures was initially separated into two classes (one containing 64 futures and the other with 80 futures), each being divided into two further classes. Within the group with 80 futures, the futures developed using the superior American marten habitat in the Doig and Moberly areas and moose habitat in the Doig area with the FFD algorithm (50 futures) were differentiated from the futures of moose habitat in the Moberly area with the FFD algorithm and of marten habitat in the Fort Nelson area (30 futures). Of the class with 64 futures, the futures developed using superior moose habitats in the Doig and Moberly areas using the SA algorithm (only the harvesting at culmination MAI for the Moberly area), and at Fort Nelson with the FFD algorithm (25 futures) were isolated from the futures developed using average HSI (24 futures) and the futures developed for the superior moose habitat using SA for the Moberly (only the preset harvesting age) and Fort Nelson areas (15 futures). Fig 5.3 a. The significant periods revealed by cluster analysis using Ward’s minimum variance method (SSE is the sum of square of errors) 118 Figure 5.3 b. The significant models revealed by cluster analysis using Ward’s minimum variance method (SSE is the sum of square of errors). The PCA partially confirmed the results of the cluster analysis and repeated measure analysis and divided the planning period into three periods (Fig. 5.4a) with the significant separation moments occurring after 11 and 59 years. The separation was mainly the result of the first four eigenvectors, which covered more than 90% of variance (Table 5.4). The first eigenvector, responsible for two thirds of the variance, was associated with 98% of the futures (142 out of 144) while the second eigenvector, responsible for 13% of the variance, was associated with the futures developed using superior moose habitats (except for the Moberly area) and for American marten, based on FFD algorithm for the Fort Nelson area (60 futures). The remaining eigenvectors accounted for 13% of the variance and were associated with the futures developed for superior moose habitats (50 futures for the third eigenvector and 40 futures for the fourth eigenvector). PCA offered a different perspective than cluster analysis, as the first two eigenvectors failed to provide a clear separation of the 144 futures, regardless of the study area or harvesting algorithm (Fig 5.4b). Even the combination study area-harvesting algorithm failed to produce an unambiguous grouping of the futures, as each combination presented at least one set of futures outside the overall pattern (e.g., Fort Nelson and SA). The lack of partitioning 119 was supported by the relative distribution of the variance among the selected eigenvectors (63% for the first eigenvector and 26% for the second). The years associated with the selected eigenvectors (the first three eigenvectors) covered the entire planning period, with the first eigenvector being related to the last two thirds of the century while the second eigenvector was associated with the first third (Table 5.4). Table 5.4 Principal components analysis Statistics Principal component 1 2 3 4 Years Proportion of variance 0.658 0.133 0.077 0.055 Significant All futures except the two developed using average HSI and SA for moose on Doig area 1. All futures for moose developed using HSI>0.5 and SA 2. Futures developed for moose using HIS >0.5 and FFD for Doig and Fort Nelson areas 3. Futures developed for marten using HSI>0.5 and FFD for Fort Nelson area 1. All futures for Moberly area developed using average HSI and SA 2. All futures developed for moose using HSI>0.5 except for Doig area 1. Futures for moose developed using HSI>0.5 for Fort Nelson area 2. Futures for moose developed using HSI>0.5 and FFD for Doig area 3. Three unrelated futures developed for moose using average his Models Proportion of variance 0.629 0.259 0.057 0.0168 Selected years 1 33-100 3-39 2-10 - 120 Figure 5.4 PCA separating the significant periods (a) and the significant combination of study area harvesting algorithm (b). The CDA (Fig 5.5) revealed that the harvesting algorithms play a significant role in shaping the future of the wildlife habitat (p<0.0001 for all three multivariate tests). Although the first four eigenvectors had canonical correlations exceeding 0.99, only the first two were selected for interpretation, as each explained more than 5% of the variance (Table 5.5). However, not all years seem to have had a significant impact on the harvesting algorithm (p>0.1), with the first two thirds of the planning period and three other years being important (years 1 – 63, 69, 85 and 86). Only 32 of the 66 significant years were selected for interpretation; these years had a b a 121 correlation with the canonical variate of >0.32, or the standardized canonical coefficient was larger than half the maximum standardized canonical coefficient (Table 5.5). Among the 32 years, no year had a correlation with the second canonical variate of >0.32, and the interpretation was therefore based only on the standardized canonical coefficient. Regardless of the year, only years 85 and 86 had both the standardized canonical coefficients and the correlation with the first canonical variate fulfilling the selected cut-off values. Table 5.5 Canonical discriminant analysis Statistics Canonical variate 1 Canonical variate 2 Canonical correlation 0.999 0.999 Proportion of variance 0.788 0.186 Year Correlation between year and canonical variate Standardized canonical coefficients Correlation between year and canonical variate Standardized canonical coefficients 7-26 0.32 to 0.50 -23924 to 13726 -0.1 to 0.23 -5927 to 10038 31 0.1881 47061 0.2581 -3493 34 0.217 37497 0.1574 -9640 37 0.1166 101867* 0.1916 11750* 43 -0.179 62186 0.0886 -2796 50 -0.218 100920 -0.001 4191.7 52 -0.221 58141 -0.012 -3425 56 -0.226 59579 0.0023 2693 63 -0.306 50558 0.0434 -4804 85 -0.35 -88893 -0.04 7291 86 -0.33 82041 -0.06 -5418 * largest standardized canonical coefficient 122 Figure 5.5 Canonical discriminant analysis indicating the impact of harvesting algorithm on delineating the significant years during the planning period 5.4 Discussion The fulfillment of all major assumptions related to the development of the oil drilling, of the harvesting requirements imposed by BC forest regulations, and of the habitat thresholds indicated that the investigation of the set of possible and equally likely futures could reveal the similarities or differences between the futures or among the years of each future. The repeated measures analysis (both univariate and multivariate) showed that the futures differed significantly. The large array of diverse and possible habitat development produced by the multimodel framework was desirable from two perspectives: firstly to ensure that the models used to build the measure of habitats was not bounded by a narrow homogeneous set of potential futures and secondly to provide confidence in the results. The difference between the futures was likely not an artifact of the violation of one of the least square method assumptions as no outlier was identified among the 144 futures. The variability of the futures increased the likelihood of occurrence of a significant change in habitat suitability during the planning period, as indicated by the univariate and multivariate methods. Similarly, with the investigation of the futures, the delineation of the planning period in distinct sub-periods was not the product of an 123 outlier, as no year was identified as an outlier. However, the lack of outlier years was not surprising, as the entire research was developed within a repeated measure framework, which emphasized that the habitat did not change significantly from year to year. The normality assumption that could have impacted the results was violated but the high significance of the multivariate analyses (p<0.001) meant that the requirement for a normal distribution did not greatly affect the interpretation, which was based on an exploratory rather than explanatory context. The univariate and multivariate repeated measures analyses revealed the existence of at least three distinct periods during the century planning horizon. The disagreement between the univariate and multivariate results lay in the number of periods (three with the univariate and four with the multivariate analyses), but the mid-century breakpoint was present in both types of analysis. The presence of a mid-century breakpoint in the dynamics of habitat suitability was not confirmed by the cluster analysis (the second breakpoint occurred at year 37) but was roughly confirmed by the PCA (the breakpoint was at year 59). The results of the PCA were also in agreement with the three periods identified by the univariate analysis, but the PCA indicated that the planning period separation would occur almost entirely in the first half of the century (years 11 and 59) while the univariate analysis suggested a partition in years 50 and 93. The change in the habitat suitability that was predicted to occur twice in the first 60 years suggests a need for the reassessment of the environmental structure, and was common to all multivariate investigations. The first significant habitat change was identified by all multivariate technique as taking place at the end of the first decade. The disagreement between different multivariate techniques over the timing of the second habitat change likely results from the large variability existent within the set of 144 futures. The variability in habitat suitability, confirmed by the repeated measures analysis, and the likelihood of change during the increase in HSI (as shown by the Scheffe test) suggested that the transition from the second to the third homogeneous periods would occur between years 37 and 59. The results supplied by the multivariate investigations were considered to be more likely than the univariate investigations, as the multivariate analyses provide a more comprehensive picture of the possible futures (Hardle and Simar 2003). The arrangement of the years within the distinct periods in a consecutive pattern, as revealed by the cluster analysis and PCA, confirmed that a significant transformation in habitat will not occur suddenly but will probably require over half of a decade or more. The interlude between the two periods would be characterized by insignificant yearly changes in the habitat suitability but the end year would be different than the beginning year. The lack of habitat delineation between the years suggests the need for an analysis of the environmental attributes 124 within a multivariate platform, which would allow the simultaneous and integrated assessment of the attributes describing the environment. The indistinguishable changes resemble the chaotic behavior of the environment as a whole (May 1977; Prigogine 1997), and draw attention to the possible existence of nonlinear processes underlining the environmental dynamics. Both univariate and multivariate investigations showed that of the two human activities examined, forest harvesting is more responsible for the changes occurring in habitat during the century long planning period. The relative lack of impact of well drilling on habitat suitability was initially indicated by the development of the futures, as the futures built using the average HSI did not delineate among the models describing the drilling, and the futures built using stands with HSI>0.5 revealed no impact of drilling after the first 40 years. The importance of forest harvesting as a determinant of habitat was also reflected at the procedural level, as the harvesting algorithm provided a significant role in deciding the development of moose and American marten habitat (the CDA results). The perspective of change, understood as the method used to represent the development of the habitat suitability, has an impact on the conclusions even when the species and habitat quantification remained constant. From this perspective, the futures assessed using the average HSI were separated from those evaluated using stands with HSI>0.5 (cluster analysis, PCA and CDA). The findings confirm the results of Seppelt and Richter (2005), who found that the choice of methods could alter the results. This undesirable conclusion reinforced the need for multiple methods, as advocated in the multi- model structure (Fig. 5.2). The habitat of the two species presented a similar evolution; even though locally the trends and variability were different. The analogous response of moose and American marten to the two human activities, as clearly indicated by the CDA and PCA, was unexpected as the habitat requirements of the two species are unrelated, and differs from the results obtained by Nitschke (2008) who found that only American marten had a significant change in the average habitat suitability. I confirmed the difference in the behavior of the two species mentioned by the Nitschke, but the results indicate that the habitat change was related to forest harvesting (more specifically the algorithm used to schedule the harvestings) rather than to other human activities [such as seismic lines or roads, which accounted only for 1% of the land-use change, as indicated by the Oil and Gas Commission of British Columbia database (Oil and Gas Commission 2005)]. 125 The investigation indicated that a significant increase in the annual allowable cut does not necessarily translate directly into a significant decrease in habitat suitability, and may even lead to an increase in habitat suitability (as for FFD in the Fort Nelson and Doig areas). Nevertheless, where human activities are changing the habitat of the two species, the shift seems to concentrate on the better habitat, rather than on the overall landscape, as revealed by PCA and cluster analysis. While the stands with HSI>0.5 appear to be sensitive to the combined effect of the two human activities, the study area seemed to have no influence on the evolution of habitat for either species, as determined by all the multivariate analyses. This suggests that the results could be extended from the 1.7 million ha sampling area to the entire northeastern British Columbia region (17 million ha). 5.5 Conclusions To assess the cumulative impact of drilling and forest harvesting on the habitat of moose and American marten in northeastern British Columbia I developed a multi-model analytical framework that enabled the identification of the moments when the habitat would significantly change and the attributes responsible for those changes. Based on the recommendation of the Treaty 8 First Nations, three areas covering 1.7 million ha were selected to represent the 17 million ha of northeastern British Columbia. The use of two species and two human activities (oil drilling and forest harvesting) ensured that the multi-model framework could be generalized to assess an array of human activities and VECs. The multi-model framework required the development of a set of independent and equally likely futures of the attributes that quantify the VEC (in our case, the habitat suitability indices for the two species). The set of futures, with each future representing a possible evolution of the environment, would provide the confidence interval for the moment when significant environmental changes would occur. The set of futures, developed as a combination of the probable and independent models quantifying the evolution of the human activities, led to 240 futures (5 drilling models x 2 forest harvesting algorithms (simulated annealing or adjusted first-fit decreasing algorithm) x 2 harvesting ages (age at culmination MAI or preset age) x 3 study areas x 2 HIS statistics (average HSI or area with HIS>0.5) x 2 VECs). Of the 240 futures, only 144 (i.e. 120 + 120/5) were considered in the analysis, as the five drilling models had an indistinguishable impact on the average HSI. The univariate and multivariate analyses revealed the existence of at least three distinct periods within the next century, with unanimous agreement that management plans should be re- 126 evaluated close to the mid-century to assess whether or not the path followed by wildlife habitat is on the desired trajectory. The multivariate investigation indicated a significant change in habitat after the first decade, requiring a first revision of the planning strategies to evaluate whether the multi-model framework supplied accurate results and to adjust the land-use management plans according to the preferred goals. Regardless of the investigation technique (univariate or multivariate), the main human activity impacting wildlife habitat is forest harvesting, especially the algorithm used to schedule the harvesting (the same result was obtained by Nitschke (2008) using generated data). The human activities seem to affect the better moose habitats (as indicated by cluster analysis and PCA) but the influence is relatively weak (less than 25%). The results reveal that an insignificant increase in landscape disturbance and habitat suitability can lead to a significant increase in the annual allowable cut. Therefore, the cumulated impact of human activities does not necessarily depend on the intensity of the environmental change, but, within certain limits, on the arrangement of the disturbance induced by the activities, confirming the findings of Strimbu et al. (2005). Besides the significant impact of the spatial organization of the disturbances on the attributes quantifying the VECs, this research also substantiates the results of Seppelt and Richter (2005), in the sense that the methodology itself could play a significant role in the analysis, and of Prigogine (1997) and May (1974), in the sense that the environmental dynamics seems to be driven by nonlinear processes. Future research should focus on two directions: one theoretical and one practical. From the practical perspective, future work needs to expand the proposed multi-model framework to incorporate other human activities and VECs (such as agriculture, tourism and mining). From a theoretical perspective, future research should concentrate on the sensitivity of the framework to the violation of the assumptions and on the impact of data distributions on the results. 127 5.6 References 1. Agricultural Land Commission and Oil and Gas Commission. Agricultural Land Commission Act. 36. 2002. 2. Anonymous. Maximum Disturbance Review Criteria. 1-8. 2004. Fort St. John, Oil and Gas Commission. 3. Bailey, S.-A., R.H. Haines-Young, and C. Watkins. 2002. Species presence in fragmented landscapes: modelling of species requirements at the national level. Biological Conservation 108(3):307-316. 4. Balke, N.S. and R.J. Gordon. 1986. Historical data. P. 781-850 in The American Business Cycle, Gordon, R.J. (ed.), The University of Chicago Press, Chicago IL. 5. Beanlands, G.E. and P.N. Duinker. 1983. An ecological framework for environmental impact assessment in Canada. Institute for Resources and Environmental Studies, Dalhousie University, Halifax. 6. Belsley, D.A., E. Kuh, and R.E. Welsch. 1980. Regression Diagnostics: Identifying Influential Data and Sources of Collinearity. John Wiley and Sons, New York. 292 p 7. Bettinger, P., D. Graetz, K. Boston, J. Sessions, and W.D. Chung. 2002. Eight heuristic planning techniques applied to three increasingly difficult wildlife planning problems. Silva Fennica 36(2):561-584. 8. Boltzmann, L. 1995. Lectures on gas theory. Dover Publications. 490 p 9. Borland Software Corporation. Delphi. [7.0]. 2001. Cupertino, CA. 10. Box, G.E.P., G.M. Jenkins, and G.C. Reinsel. 1994. Time Series Analysis: Forecasting and Control. Prentice Hall, Englewood Cliffs, NJ. 11. Brockwell, P.J. and R.A. Davis. 1996. An Introduction to Time Series and Forecasting. Springer Verlag, New York. 420 p 12. Cameron, A.C. and P.K. Trivedi. 1998. Regression analysis of count data. Cambridge University Press, New York NY. 432 p 13. Choi, B. 1992. ARMA Model Identification. Springer-Verlag, New- York. 14. Cleveland, W.S. 1972. Inverse Autocorrelations of A Time Series and Their Applications. Technometrics 14(2):277-293. 15. Coffman, E.G., C. Courcoubetis, M.R. Garey, D.S. Johnson, P.W. Shor, R.R. Weber, and M. Yannakakis. 2000. Bin packing with discrete item sizes, Part I: perfect packing theorems and the average case behavior of optimal packing. SIAM Journal of Discrete Mathematics 13(3):384-402. 128 16. Costanza, R., A. Voinov, B. Roelof, T. Maxwell, F. Villa, L. Wainger, and H. Voinov. 2002. Integrated Ecological Economic Modeling of the Patuxent River Watershed, Maryland. Ecological Monographs 72(2):203-231. 17. Cressie, N. 1993. Statistics for Spatial Data. John Wiley&Sons, New York NY. 900 p 18. Crowder, M.J. and D.J. Hand. 1990. Analysis of repeated measures. Chapman and Hall, London. 257 p 19. Crowe, K. and J. Nelson. 2003. An indirect search algorithm for harvest-scheduling under adjacency constraints. Forest Science 49(1):1-11. 20. Cuddington, J.T. and D.L. Moss. 2001. Technological change, depletion, and the US petroleum industry. American Economic Review 91(4):1135-1148. 21. Davis, L.S., K.N. Johnson, T.E. Howard, and P. Bettinger. 2001. Forest management. McGraw-Hill, New York. 804 p 22. Diggle, P., P. Heagerty, K.Y. Liang, and S.L. Zeger. 2002. Analysis of longitudinal data. Oxford University Press, New York. 396 p 23. Energy Information Administration. Spot Prices for Crude Oil and Petroleum Products. 2005. U.S. Department of Energy. 24. Environmental Systems Research Institute. ESRI Hydro Data Model. 2001. Redlands CA, Environmental Systems Research Institute. 25. Environmental Systems Research Institute, I. ArcGIS. [9.3]. 2008. Redlands CA, ESRI. 26. Fryxell, J.M., J.B. Falls, E.A. Falls, R.J. Brooks, L. Dix, and M.A. Strickland. 1999. Density dependence, prey dependence, and population dynamics of martens in Ontario. Ecology 80(4):1311-1321. 27. Gilbert, W.J. 1976. Modern algebra with applications. John Wiley and Sons, New York. 348 p 28. Government of British Columbia. 2004. Petroleum and Natural Gas Act. 362. 29. Government of Canada. 1992. Canadian Environmental Assessment Act. 37, 74 p 30. Greenhouse, S. and S. Geisser. 1959. On methods in the analysis of profile data. Psychometrika 24(2):95-112. 31. Grimmett, G.D. and D.R. Stirzaker. 2002. Probability and Random Processes. Oxford University Press, New York. 600 p 32. Hardle, W. and L. Simar. 2003. Applied multivariate statistical analysis. Springer-Verlag, New York. 458 p 33. Hawkins, D.M. and N. Cressie. 1984. Robust krigingÇöA proposal. Mathematical Geology 16(1):3-18. 129 34. Henze, N. and B. Zirkler. 1990. A Class of Invariant Consistent tests for Multivariate Normality. Communications in Statistics - Theory and Methods 19(10):3595-3617. 35. Hirzel, A.H., G. Le Lay, V. Helfer, C. Randin, and A. Guisan. 2006. Evaluating the ability of habitat suitability models to predict species presences. Ecological Modelling 199(2):142-152. 36. Holland, S.P. 2008. Modeling Peak Oil. Energy Journal 29(2):61-79. 37. Hungerford, T.W. 1990. Abstract algebra: an introduction. Saunders College Publishing, Orlando. 575 p 38. Huynh, H. and L.S. Feldt. 1970. Conditions Under Which Mean Square Ratios in Repeated Measurements Designs Have Exact F-Distributions. Journal of the American Statistical Association 65(332):1582-1589. 39. Iledare, O.O. and A.G. Pulsipher. 1997. The modelling of the petroleum exploration and extraction process for policy analysis: a case study of the Louisiana onshore region. Pacific and Asian Journal of Energy 7(1):21-38. 40. Information Services AEUB. Drilling and well information. 2005. Calgary AB, Alberta Energy and Utility Board. 41. Islam, M.Q. and M.L. Tiku. 2004. Multiple Linear Regression Model Under Nonnormality. Communications in Statistics - Theory and Methods 33(10):2443-2467. 42. Jager, H.I., E.A. Carr, and R.A. Efroymson. 2006. Simulated effects of habitat loss and fragmentation on a solitary mustelid predator. Ecological Modelling 191(3-4):416-430. 43. Johnson, D.S. 1973. Near-optimal bin packing algorithms. PhD Massachusetts Institute of Technology, Dept. of Mathematics. 402 p 44. Journel, A.G. 1983. Nonparametric estimation of spatial distributions. Mathematical Geology 15(3):445-466. 45. Jukes, W., D. Menzies, D. Rosen, G. Taylor, R. Saint-Jean, and J. Beale. 2004. Sustainable Forest Management Plan. Fort St. John, Ministry of Forests. 402 p 46. Kimmins, J.P., D. Mailly, and B. Seely. 1999. Modelling forest ecosystem net primary production: the hybrid simulation approach used in FORECAST. Ecological Modelling 122(3):195-224. 47. Kitanidis, P.K. 1997. Introduction to geostatistics. Cambridge University Press, Cambridge. 272 p 48. Kroll, A.J. and J.B. Haufler. 2006. Development and evaluation of habitat models at multiple spatial scales: A case study with the dusky flycatcher. Forest Ecology and Management 229(1-3):161-169. 49. Lanza, A., M. Manera, M. Grasso, and M. Giovannini. 2005. Long-run models of oil stock prices. Environmental Modelling & Software 20(11):1423-1430. 130 50. Larson, M.A., F.R. Thompson, J.J. Millspaugh, W.D. Dijak, and S.R. Shifley. 2004. Linking population viability, habitat suitability, and landscape simulation models for conservation planning. Ecological Modelling 180(1):103-118. 51. Liu, C.R., P.M. Berry, T.P. Dawson, and R.G. Pearson. 2005. Selecting thresholds of occurrence in the prediction of species distributions. Ecography 28(3):385-393. 52. Lockwood, C. and T. Moore. 1993. Harvest scheduling with spatial constraints: a simulated annealing approach. Canadian Journal of Forest Research 23(3):468-478. 53. Lu, N. and D.L. Zimmermann. A nonparametric approach to testing for directional symmetry of spatial correlation. 234. 1994. Iowa City IA, Department of Statistics and Actuarial Science, University of Iowa. 54. Mardia, K.V. 1970. Measures of Multivariate Skewness and Kurtosis with Applications. Biometrika 57(3):519-530. 55. May, R.M. 1974. Biological Populations with Nonoverlapping Generations - Stable Points, Stable Cycles, and Chaos. Science 186(4164):645-647. 56. May, R.M. 1977. Thresholds and Breakpoints in Ecosystems with A Multiplicity of Stable States. Nature 269(5628):471-477. 57. Meidinger, D. and J. Pojar. 1991. Ecosystems of British Columbia. Victoria, BC Ministry of Forests. 179 p 58. Menard, A., P. Dube, A. Bouchard, C.D. Canham, and D.J. Marceau. 2002. Evaluating the potential of the SORTIE forest succession model for spatio-temporal analysis of small-scale disturbances. Ecological Modelling 153(1-2):81-96. 59. Metropolis, N., A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller. 1953. Equations of State Calculations by Fast Computing Machines. Journal of Chemical Physics 21(6):1087-1092. 60. Ministry of Forests. 1995. Biodiversity Guidebook. 1. Victoria BC, BC Minstry of Forests. 61. Ministry of Forests. TIPSY. [3.2b]. 2005a. Victoria BC. 62. Ministry of Forests. WinVDYP. [7.26c]. 2005b. Victoria BC. 63. Ministry of Sustainable Resource Management. TRIM/TRIM II positional files - basemap. [2004], http://srmwww.gov.bc.ca/bmgs/catalog/bcgs_indexes.htm. 2004. Victoria BC, Ministry of Sustaianble Resource Management. 64. Mojena, R. 1977. Hierarchical Grouping Methods and Stopping Rules - Evaluation. Computer Journal 20(4):359-363. 65. Neter, J., M.H. Kutner, C.J. Nachtsheim, and W. Wasserman. 1996. Applied linear statistical models. WCB McGraw-Hill, Boston. 1408 p 131 66. Nitschke, C.R. 2008. The cumulative effects of resource development on biodiversity and ecological integrity in the Peace-Moberly region of Northeast British Columbia, Canada. Biodiversity and Conservation 17(7):1715-1740. 67. Oil and Gas Commission. 2005. Geographic Information System Data. Oil and Gas Commission of British Columbia. 68. Ottaviani, D., G.J. Lasinio, and L. Boitani. 2004. Two statistical methods to validate habitat suitability models using presence-only data. Ecological Modelling 179(4):417-443. 69. Park, T. and S.-Y. Lee. 2004. Model Diagnostic Plots for Repeated Measures Data. Biometrical Journal 46(4):441-452. 70. Peterson, R.O. 1997. Ecological Studies of Wolves on Isle Royale. Houghton MI, School of Forestry and Wood Products, Michigan Technological University. 17 p 71. Phillips, P.C.B. and P. Perron. 1988. Testing for a Unit Root in Time Series Regression. Biometrika 75(2):335-346. 72. Prigogine, I. 1997. The end of certainty. The Free Press, New York. 240 p 73. Rao, R.D. 2000. An integrated modelling framework for exploration and extraction of petroleum resources. Resources Policy 26(3):133-143. 74. Rencher, A.C. 2002. Methods of Multivariate Analysis. John Wiley and Sons, New York. 708 p 75. Ringlund, G.B.r., K.E. Rosendahl, and T. Skjerpen. 2008. Does oilrig activity react to oil price changes An empirical investigation. Energy Economics 30(2):371-396. 76. Rocke, D.M. and D.L. Woodruff. 1996. Identification of outliers in multivariate data. Journal of the American Statistical Association 91(435):1047-1061. 77. Roloff, G.J. and J.B. Haufler. 1997. Establishing Population Viability Planning Objectives Based on Habitat Potentials. Wildlife Society Bulletin 25(4):895-904. 78. Roloff, G.J. and B.J. Kernohan. 1999. Evaluating Reliability of Habitat Suitability Index Models. Wildlife Society Bulletin 27(4):973-985. 79. Romito, T., K. Smith, B. Beck, J. Beck, M. Todd, R. Bonar., and R. Quinlan. 1999. Moose winter habitat. 27, Hinton AB, Foothils Model Forest. Habitat Suitability. 6 p 80. SAS Institute Inc. SAS Software version 9.1.3. 2004. Cary, SAS Institute Inc. 81. Schwager, S.J. and B.H. Margolin. 1982. Detection of Multivariate Normal Outliers. The Annals of Statistics 10(3):943-954. 82. Seppelt, R. and O. Richter. 2005. "It was an artefact not the result": A note on systems dynamic model development tools. Environmental Modelling & Software 20(12):1543-1548. 83. Stelfox, J.B. Alces. [Alces II]. 1999. Bragg Creek, AB, Forem Technologies. 2001. 132 84. Strimbu, B.M., G.M. Hickey, and V.G. Strimbu. 2005. Forest conditions and management under rapid legislation change in Romania. Forestry Chronicle 81(3):350-358. 85. Strimbu, B.M., G.M. Hickey, V.G. Strimbu, and J.L. Innes. 2009. On the Use of Statistical Tests With Non-Normally Distributed Data in Landscape Change Detection. Forest Science 55(1):72-83. 86. Strimbu, B.M., J.L. Innes, and V.F. Strimbu. 2008. A harvest scheduler using perfect bin- packing theorem. P. 57 Hahn, A., T. Knoke, and T. Schneider (eds.). Wald Forst Holz Zentrum, Freising DE. 87. Tabachnick, B.G. and L.S. Fidell. 2001. Using multivariate statistics. Allyn and Bacon, Needham Heights. 966 p 88. Takats, L., R. Stewart, M. Todd, R. Bonar., J. Beck, B. Beck, and R. Quinlan. 1999. American marten winter habitat. 30, Hinton AB, Foothills Model Forest. 8 p 89. Thode, H.C. 2002. Testing for normality. Marcel Dekker, New York. 368 p 90. Tran, Z.V. 1997. Estimating Sample Size in Repeated-Measures Analysis of Variance. Measurement in Physical Education and Exercise Science 1(1):89-102. 91. Tsay, R.S. and G.C. Tiao. 1985. Use of Canonical Analysis in Time Series Model Identification. Biometrika 72(2):299-315. 92. Tukey, J. 1977. Exploratory Data Analysis. Addison Wesley. 688 p 93. Walls, M.A. 1994. Using A Hybrid Approach to Model Oil and Gas-Supply - A Case- Study of the Gulf-Of-Mexico Outer Continental-Shelf. Land Economics 70(1):1-19. 94. Walls, M.A. 1992. Modeling and forecasting the supply of oil and gas : A survey of existing approaches. Resources and Energy 14(3):287-309. 95. Wigmosta, M.S., L.W. Vail, and D.P. Lettenmaier. 1994. A Distributed Hydrology- Vegetation Model for Complex Terrain. Water Resources Research 30(6):1665-1679. 96. Wilks, S.S. 1963. Multivariate statistical outliers. Sankhya, Series A 25(4):407-426. 97. Yamasaki, S.H., R. Duchesneau, F. Doyon, J.S. Russell, and T. Gooding. 2008. Making the case for cumulative impacts assessment: Modelling the potential impacts of climate change, harvesting, oil and gas, and fire. Forestry Chronicle 84(3):349-368. 98. Zar, J.H. 1996. Biostatistical analysis. Prentice - Hall, Upper Saddle River NY. 960 p 133 6. CONCLUSIONS Cumulative impact assessments (CIA) focus on the mitigation of the undesirable effects induced by all human activities occurring in an area (Therivel and Morris, 2001). The actual and potential detrimental consequences of different projects are estimated using a variety of methodologies ranging from checklists to mathematical or statistical models (Therivel and Morris, 2001). Among the mathematical and statistical models used for CIA, the most popular are based on modular modeling (Voinov et al., 2004) and heuristic techniques (Stakhiv, 1988). The two methodologies have the advantage that they supply accurate results and can represent the intricate relationships within the environment. The complex interactions between the different components of the environment have typically led to the use of a single analytical platform (i.e., one theoretical approach) for the CIA. When the interest lies in the disturbance induced by a set of current or planned economic activities in an area, the CIA methodology focuses on providing a constellation of measures that provides an indication of whether or not the total instability will exceed acceptable environmental limits. However, the single platform, while providing technical solutions at a specific point in time, fails to maintain the accuracy associated with mathematical and statistical models when a spatial-temporal forecast is required, as the outcomes have a probability of occurrence asymptotically null. The reduced accuracy requires that the assessment of the impact of future activities should not be treated in the same way as for existing activities. The focus of the research presented in this thesis is to illustrate a new approach for CIA that incorporates the spatial and temporal distribution of future activities that could result in a significant disturbance to the environment. The new approach presented in the thesis suggests a need to change the CIA paradigm: rather than performing a detailed assessment of an improbable future the CIA should focus on identifying the patterns associated with undesired states of environment. The identification of patterns related to undesirable environmental conditions can be performed by analyzing a multitude of futures rather than one future, an approach similar to that used in statistical thermodynamics (Boltzmann, 1995). The analysis of a set of futures translates the question of whether or not the attributes describing an environmental state meet the regulatory conditions, which is the focus of current CIA investigations based on single analytical platforms, to whether or not the environment could reach an unacceptable state. To ensure the valid identification of the patterns delineating significantly different states of the environment, the attributes describing each future have to fulfill the regulatory constraints at any moment and at any location. To use an analytical 134 framework compatible with statistical thermodynamics, each future should be independent and with the same probability of occurrence. The fulfillment of the regulatory constraints ensure that the future could occur (i.e., it does not violate legislation), while the requirements of independence and equal probability ensure that standard statistical techniques can be used to identify the constellation of attributes defining an environmental pattern. Based on the set of futures, the paradigm shift promotes identification of the main environmental attributes affected by human activities and, additionally, the moments when significant environmental changes could occur. I selected three areas from northeastern British Columbia (Moberly, Doig and Fort Nelson) to illustrate the proposed CIA method The human activities investigated were forest harvesting and oil and gas drilling, the two industries dominating the region, while the attributes describing the environment were moose (Alces alces) and American marten (Martens americana), the two valued ecosystem components (VEC) identified by Treaty 8 First Nations as significant from social and environmental perspectives. A set of possible independent and equally likely futures was developed using different spatial and temporal analytical tools to model the two human activities. To this end, Chapters 2 and 3 focus on forecasting the oil and gas developments (Chapter 2) and the forest harvesting (Chapter 3). A total of 20 spatially and temporally explicit futures were developed for each area (5 futures oil and gas drilling x 2 possible harvesting algorithms x 2 harvesting ages). The temporal evolution of the oil and gas drilling seems to depend not only the model and data but also on the time step that is used. While the influence of the size of the time step, or more generally the magnitude of the elementary modeling unit, was noted by previous researchers (Chaubey et al., 2005; Claessens et al., 2005; Vazquez et al., 2002; Lanza et al., 2005), the findings of Chapter 2 reveal that the models do not converge, therefore a true model cannot be identified. Consequently, from a theoretical perspective, the analytical tools are unable to identify an accurate model, even when all the modeling assumptions and requirements are met. This result emphasizes the need for a paradigm shift in cumulative impacts analyses, as statistical investigations of sets of futures are likely to be more reliable than a single probable but likely inaccurate future. To harmonize the forest harvesting methodology, which uses annual plans, with the oil and gas drilling, the CIA was performed using only results based on an annual time step (i.e., for both forest harvesting and prediction of the number of wells). The planning period for the foreseeable forest harvesting and oil and gas drilling activities occurring in the region was 100 years. The results from Chapter 3 reveal that the methodology used to schedule the harvestings has an 135 economic and environmental impact, as significantly greater harvest volumes are not necessarily associated with significantly larger disturbances. Furthermore, the greater harvested volumes did not result in a significant increase in the annual variation, suggesting that it is possible to use algorithms that are faster and supply results that meet both economic and ecological constraints. The combined impact of forest harvesting and oil and gas drilling on the two VECs is investigated using univariate repeated measurements (Chapter 4) and multivariate techniques (Chapter 5). The two measures assessing the VEC dynamics during the planning period (average HSI and area with HSI>0.5) revealed that the maintenance of the present trends in forest planning and petroleum drilling will likely lead to futures that fulfill the regulatory requirements. The fulfillment of the regulations does not necessarily ensure that the attributes describing the environment are maintained in a socially and ecologically acceptable state, as shown by Strimbu et al. (2009). Nevertheless, in the eventuality that the set of futures leads to desirable future environments, an increase in the exploitation of the two natural resources (volume of harvested timber and number of petroleum drilling wells) would not be reflected in a decrease of the HSI by more than 15% of the present values. Therefore, whether or not the trend and intensity of the current forest harvesting and petroleum drilling is maintained or increased, as long as the regulations are met, the populations of the two species would likely not experience a decrease outside the natural range of variation (Chapters 4 and 5). The habitats of the two species evolve in a similar way, despite local differences in the trend and variability during the planning period. The univariate repeated measure analysis (Chapter 4) reveals the existence of three distinct periods within the century long planning period. The separation of the planning period into distinct time intervals is confirmed by the multivariate analysis (Chapter 5), with the emphasis on the half century moment. The middle of the planning period separated the environment into two well-defined periods in term of variability (Chapter 4), number of periods (one for univariate and two for multivariate) and trend (the two measures of the VECs having decreased in the first period vs. increasing or remaining relatively constant in the second period). Forest harvesting seems to have a greater impact on the habitat of the two species than oil and gas drilling. The influence of forest harvesting was reflected not only in the magnitude of the harvesting but in the scheduling of the interventions in the forest, as simulated annealing supplied smaller annual allowable cuts than first-fit decreasing algorithm but with a larger negative footprint on the habitat of moose and American marten (Chapter 5). The impact of human activities appears to be concentrated on the better habitats and not on the overall 136 landscape. Finally, the need for a change in the CIA paradigm was suggested not only by the evolution of the environment that could follow significantly different trajectories (Chapters 4 and 5 indicated the existence of at least four classes of futures making the detailed analysis of one future useless), but also by the presence of cycles within the habitat evolution, indicating the possible chaotic behavior of the environment as a whole (Prigogine, 1997; May, 1977). The present research enhances the environmental impact assessment methodology (Morris and Therivel, 2001) by adding a new dimension to the investigation of the attributes used to describe the environment: the unknown future. The results of the research demonstrated that the spatial and temporal dynamics of the potential human activities should be investigated separately from known environmental processes. The approach presented in this thesis offers an alternative to current CIA methods (Therivel, 2001; Morris and Emberton, 2001; Hodson et al., 2001) that incorporate and analyze future human-induced disturbances using the same methodology for known and predicted activities. The approach proposed here has both strengths and weaknesses. The main strength of the proposed paradigm lies in the identification of the moments when significant decisions should be made and the environmental attributes that should be considered during the decision process. The research has provided a relatively accurate means to determine the attributes responsible for potential environmental changes as well as the confidence intervals for the periods when the attributes could change. However, the multi-model analytical platform has two weaknesses. Firstly, the requirement for a set of independent and equally likely futures leads to polynomial complexity (Sipser, 1996) in terms of the number of futures (i.e., the development of m x n futures when the CIA is investigating m human activities using n different analytical frameworks). The large number of futures requires an intensive theoretical and computational effort, making the creation of the set of independent and equally likely futures difficult, and best done by a large interdisciplinary team of experts. Secondly, while the computational and theoretical difficulties are expected to decrease with further developments in information technology, the analytical framework of the proposed paradigm will not provide the technical recommendations normally associated traditional CIA techniques, but will only indicate when significant decisions should be made that will enable an acceptable environmental trajectory. The main approach used to address the objectives of the CIA (mitigation, remediation and compensation of the impact of human activities occurring in an area) consists of ensuring that the constraints associated with different VECs are met (Stakhiv, 1988; Voinov et al., 2004; 137 Costanza et al., 2002). The analytical tools used to provide solutions to the restrictions imposed by regulations range from operational research (such as linear programming or heuristic techniques) to modular modeling (Voinov et al., 2004; Kimmins et al., 1999). Therefore, the aim of the current CIA methodologies is to supply the technical details that ensure that human- induced disturbances conform to environmental legislation as well as to the requirements defined by social inputs. Current CIA methodologies answer technical questions related to the design of the human activities, but fail to provide accurate information about the spatial and temporal response of the environment to further activities. A natural enhancement of current CIA methods would therefore be the development of an analytical platform that could incorporate the forecasting of different human activities while taking into account the accuracy of the CIA. The development of an analytical platform for an accurate CIA could start by developing several equally likely responses of the environment to possible human activities. Each environmental response would then have a CIA undertaken for it using the current methodologies. 6.1 Status of the hypotheses and significance of the research The present research is based on three assumptions: • The future is unknown; • The theoretical framework used to quantify the attributes representing the future does not change; • The influence of human activities on the environment can be separated from naturally occurring events. These three assumptions served as a basis for testing two hypotheses proposed by the CIA paradigm and further served as the foundation for the development of the set of futures. The two hypotheses tested in the present study were: • There is an identifiable pattern associated with each state of the environment. • The chosen analytical tools play a significant role in the description and forecast of the environment. The first hypothesis has been extensively investigated using single analytical platforms (Baskent and Jordan, 1995; Chapin et al., 1998; Enquist et al., 2002; Roberts and Gilliam, 1995; Caplat et al., 2008) and has revealed the existence of a relationship between an association of different attributes and a specific environmental state. The research presented here has not only 138 enhanced the findings of the single analytical approach by confirming the existence of the patterns within the set of attributes used to describe the state of environment but has also associated a measure of confidence with the presence or absence of the respective patterns. The second hypothesis was initially investigated by Seppelt and Richter (2005), who found that different software supplies different solutions to the simple Lotka-Volterra prey-predator model (Volterra, 1926). My results extend the work of Seppelt and Richter by showing that it is not only the software that can influence the results but also the methodology itself (Chapter 2). The findings of Chapter 2 raise the question of whether or not the current scientific investigations, which assume that knowledge should be drawn only by using valid analytical methods applied to data (Bacon, 1855), should be reevaluated. As different results are reached by using same methodology but different methodological details, a scientific investigation should consider the inclusion of a component describing the analysis itself, not only the conclusions founded on the respective analytical tools. The research presented here has three novel elements. Firstly, a new paradigm for CIA is proposed (Chapter 4 and Chapter 5). A second element is the development of a new forest harvesting scheduler that demonstrates that the harvesting of greater volumes is not necessarily associated with a greater land clearing, and could lead to better overall environments (e.g., average HSI for moose and American marten) or landscapes with a greater amount of better habitat (e.g., areas with HSI>0.5). Furthermore, the research has determined the maximum amount of volume that can be harvested in the absence of spatial constraints, a limit difficult to compute for large size problems using standard forest harvest scheduling techniques, such as mixed integer programming or heuristic methods. The third element of novelty is the new perspective on the impact of the analysis on the results, which demonstrated the necessity of adjustments to the theoretical framework of scientific investigation as established by Descartes (1850), Spinoza (1894), Locke (1823) and Hume (1888). 6.2 Potential applications and future research The present research offers a new perspective on the incorporation of explicit future spatial and temporal changes into complex models. Two potential applications of the proposed framework are a natural extension of the advocated analytical platform and paradigm change. The first application is the extension of the proposed paradigm change (i.e., identification of patterns within the analyzed system rather than detailed behavioral investigations of an unlikely future) to other areas involving complex modeling such as hydrology, climate dynamics and land-use 139 planning. The proposed development of an alternative paradigm for CIA involving forecasts was based on the large uncertainties associated with predictions. The paradigm supporting the identification of structures within a system rather than bijective relationships could be applied to domains involving large uncertainties, such as genetics, population dynamics and microbiology. The second application of the present research is associated with the determination of the upper harvesting limit that could lead to the development of better and faster heuristic algorithms for forest planning, especially the selection of the initial values of the parameters required by different methods. An additional possible application of the present research is in the areas of psychology, sociology, history or anthropology, by providing an alternative methodology to identify significant patterns for the respective disciplines, such as archetypes for psychology (Jung, 1959) or civilizations for history (Toynbee and Somervell, 1947). Future research based on the proposed change of the CIA paradigm could focus on two directions: one theoretical and one practical. From the practical perspective, future work should expand the proposed multimodel framework to incorporate other human activities (such as agriculture, tourism or mining) and VECs (such as water, air quality or visually quality objectives) into the CIA. From a theoretical perspective, future research should concentrate on the sensitivity of the proposed CIA framework to the violation of the analytical assumptions as well as the identification of the assumptions with the greatest impact on the results. A second avenue for further research based on the proposed change of the CIA paradigm is the incorporation of observed distributions of natural processes into the analytical platform. 140 6.3 References 1. Bacon, F. 1855. The New Organon. Oxford University Press, Oxford. 338 p. 2. Baskent,E.Z., Jordan,G.A., 1995. Characterizing Spatial Structure of Forest Landscapes. Canadian Journal of Forest Research-Revue Canadienne de Recherche Forestiere 25, 1830- 1849. 3. Boltzmann, L.1995. Lectures on gas theory. Dover Publications. 490 p 4. Caplat,P., Anand,M., Bauch,C., 2008. Symmetric competition causes population oscillations in an individual-based model of forest dynamics. Ecological Modelling 211, 491-500. 5. Chapin,T.G., Harrison,D.J., Katnik,D.D., 1998. Influence of landscape pattern on habitat use by American marten in an industrial forest. Conservation Biology 12, 1327-1337. 6. Chaubey,I., Cotter,A.S., Costello,T.A., Soerens,T.S., 2005. Effect of DEM data resolution on SWAT output uncertainty. Hydrological Processes 19, 621-628. 7. Claessens,L., Heuvelink,G.B.M., Schoorl,J.M., Veldkamp,A., 2005. DEM resolution effects on shallow landslide hazard and soil redistribution modelling. Earth Surface Processes and Landforms 30, 461-477. 8. Costanza,R., Voinov,A., Roelof,B., Maxwell,T., Villa,F., Wainger,L., Voinov,H., 2002. Integrated Ecological Economic Modeling of the Patuxent River Watershed, Maryland. Ecological Monographs 72, 203-231. 9. Descartes, R.1850. Discourse on Method. Simpkin, Marshall and Co, Edinburgh. 118 p 10. Enquist,B.J., Sanderson,J., Weiser,M.D., 2002. Modeling macroscopic patterns in ecology. Science 295, 1835-1836. 11. Hodson,M., Stapleton,C., Emberton,R., 2001. Soils, Geology and geomorphology. In: Morris,P., Therivel,R. (Eds.), Methods of environmental impact assessment. Spon Press, London, pp. 170-196. 12. Hume, D.1888. A treatise of human nature. Clarendon Press, Oxford. 709 p 13. Jung, C. G.1959. The archetypes and the collective unconscious. Routledge & Kegan Paul, London UK. 14. Kimmins,J.P., Mailly,D., Seely,B., 1999. Modelling forest ecosystem net primary production: the hybrid simulation approach used in FORECAST. Ecological Modelling 122,195- 224. 15. Lanza,A., Manera,M., Grasso,M., Giovannini,M., 2005. Long-run models of oil stock prices. Environmental Modelling & Software 20, 1423-1430. 141 16. Locke, J.1823. An essay concerning human undersrding. Oxford University Press, Oxford UK. 535 p 17. May,R.M., 1977. Thresholds and Breakpoints in Ecosystems with A Multiplicity of Stable States. Nature 269, 471-477. 18. Morris,P., Emberton,R., 2001. Ecology-overview and terrestrial systems. In: Morris,P., Therivel,R. (Eds.), Methods of environmental impact assessment. Spon Press, New York USA, pp. 243-285. 19. Morris, P. and Therivel, R.2001. Methods of environmental impact assessment. Spon Press, New York USA. 20. Prigogine, I.1997. The end of certainty. The Free Press, New York. 240 p 21. Roberts,M.R., Gilliam,F.S., 1995. Patterns and Mechanisms of Plant Diversity in Forested Ecosystems -Implications for Forest Management. Ecological Applications 5, 969-977. 22. Seppelt,R., Richter,O., 2005. "It was an artefact not the result": A note on systems dynamic model development tools. Environmental Modelling & Software 20, 1543-1548. 23. Sipser, M.1996. Introduction to the theory of computation. PWS Publishing company, Boston USA. 256 p 24. Spinoza, B.1894. Ethic. Macmillan and Co. New York 297 p 25. Stakhiv,E.Z., 1988. An Evaluation Paradigm for Cumulative Impact Analysis. Environmental management 12, 725-748. 26. Strimbu,B.M., Hickey,G.M., Strimbu,V.G., Innes,J.L., 2009. On the Use of Statistical Tests With Non-Normally Distributed Data in Landscape Change Detection. Forest Science 55 (1): 72-83 27. Therivel,R., 2001. Landscape. In: Morris,P., Therivel,R. (Eds.), Methods of Environmental Impact assessment. Spon Press, London, pp. 105-121. 28. Therivel,R., Morris,P., 2001. Introduction. In: Morris,P., Therivel,R. (Eds.), Methods of environmental impact assessment. Taylor and Francis Group, New York USA, pp. 3-19. 29. Toynbee, A. J. and Somervell, D. C.1947. A study of history. Oxford University Press, Oxford UK. 30. Vazquez,R.F., Feyen,L., Feyen,J., Refsgaard,J.C., 2002. Effect of grid size on effective parameters and model performance of the MIKE-SHE. Hydrological Processes 16, 355-372. 31. Voinov,A., Fitz,C., Boumans,R., Costanza,R., 2004. Modular ecosystem modeling. Environmental modelling and software 19, 285-304. 32. Volterra,V., 1926. Fluctuations in the abundance of a species considered mathematically. Nature 118, 558-560. 142 APPENDICES Appendix A Time-series modeling is more sensitive to the violation of assumptions than modeling of independent data (Gujarati 1995). Therefore, a model is considered finalized when the significant variables produce residuals distributed as white noise and do not exhibit nonlinear deterministic dependencies (Brockwell and Davis 1996). The two conditions, white noise, tested with Bartlett's Kolmogorov-Smirnov statistic (Bartlett 1978), and quadratic independence, tested using Q statistics (McLeod and Li 1983), ensure that time series models are not spurious (Phillips 1986) and are unlikely to be represented as a quadratic chaotic series (Sprott 2003). All the computations were performed using a maximum likelihood algorithm; the generalized Durbin-Watson test (Vinod 1973) assessed the sequence independence, Philips-Perron test evaluated the stationarity of sequences, normality was investigated using the Bera approach (Bera 1982) and studentized deleted residuals were employed to identify the outliers (Neter et al 1996). The significant variables were selected using stepwise selection criterion with a level of significance of 0.02. I demonstrated that the structure of a model can be dependent on the time step by developing time-series models for the count of active rigs (Table 7.1) and the number of wells drilled during a time-step (Table 7.2). The values from Table 7.1 and 7.2 represent the probability to reject the null hypothesis associated with each of the previously mentioned tests. The symbols from Table 1 and 2 represent the number of the active rigs, or the number of drilled wells, respectively (i.e., actual (xt1) and total (xt2)) and the crude oil prices (i.e., monthly minimum (xt3), maximum (xt4) and average (xt5)), based on New York Mercantile Exchange, at time t. 143 The time-series models built for British Columbia (Table 7.2) has all the residuals distributed as white noise but the Portmanteau test shows a possible quadratic relationship among residuals for the linear model with time-step 1, 3 and 12 months. However, the presence of a quadratic relationship was probably the result of the lack of normality as the GARCH process did not identify any significant nonlinear terms. The normality assumption was not met but because all the distributions were unimodal and did not exhibit significant asymmetry; the inference was considered valid from a distributional perspective, as suggested by Glass et al (1972) and Games and Lucas (1966). The Gamma distribution model, supplementing the linear model, used as a link function the reciprocal function. A Bartletts‘s Kolmogorov – Smirnov test revealed that all the residuals are WN. The Portmanteau test indicated the absence of a quadratic relationship among residuals. Table A.1 Time-series model for the monthly average of number of active rigs. I included an insignificant economic variable for Africa region to ensure the fulfillment of the time series assumptions as the estimator efficiency (reduced by this decision) was not relevant for the analysis. xt1 and xt2 are the monthly average for actual and total (cumulated) number of active rigs, while xt3, xt4 and xt5 are the monthly minimum, maximum and average of the NYMEX crude oil prices. Region Equation R2 Pr