A Multi-Scale Assessment of Spatial-Temporal Change in
the Movement Ecology and Habitat of a Threatened Grizzly
Bear (Ursus arctos) Population in Alberta, Canada
by
Mathieu Louis Bourbonnais
MSc, University of Victoria, 2013
BSc, University of Victoria, 2011
BA, University of Alberta, 2005
A Dissertation Submitted in Partial Fulfillment
of the Requirements for the Degree of
DOCTOR OF PHILOSOPHY
in the Department of Geography
© Mathieu Louis Bourbonnais, 2018
University of Victoria
All rights reserved. This dissertation may not be reproduced in whole or in part, by
photocopy or other means, without the permission of the author.
Supervisory Committee
A Multi-Scale Assessment of Spatial-Temporal Change in the Movement Ecology
and Habitat of a Threatened Grizzly Bear (Ursus arctos) Population in Alberta,
Canada
by
Mathieu Louis Bourbonnais
MSc, University of Victoria, 2013
BSc, University of Victoria, 2011
BA, University of Alberta, 2005
Supervisory Committee
Dr. Trisalyn Nelson, Supervisor
(School of Geographical Sciences & Urban Planning, Arizona State University)
Dr. Chris Darimont, Co-Supervisor
(Department of Geography, University of Victoria)
Dr. Farouk Nathoo, Outside Member
(Department of Mathematics and Statistics, University of Victoria)
Mr. Gordon Stenhouse, Committee Member
(Grizzly Bear Program, fRI Research)
Abstract
Supervisory Committee
Dr. Trisalyn Nelson, Supervisor
(School of Geographical Sciences & Urban Planning, Arizona State University)
Dr. Chris Darimont, Co-Supervisor
(Department of Geography, University of Victoria)
Dr. Farouk Nathoo, Outside Member
(Department of Mathematics and Statistics, University of Victoria)
Mr. Gordon Stenhouse, Committee Member
(Grizzly Bear Program, fRI Research)
Given current rates of anthropogenic environmental change, combined with the increasing lethal and non-lethal mortality threat that human activities pose, there is a vital need to understand wildlife movement and behaviour in human-dominated landscapes to help inform conservation efforts and wildlife management. As long-term monitoring of wildlife populations using Global Positioning System (GPS) telemetry increases, there are new opportunities to quantify change in wildlife movement and behaviour. The objective of this PhD research is to develop novel methodological approaches for quantifying change in spatial-temporal patterns of wildlife movement and habitat by leveraging long time series of GPS telemetry and remotely sensed data. Analyses were focused on the habitat and movement of individuals in the threatened grizzly bear (Ursus arctos) population of Alberta, Canada, which occupies a human-dominated and heterogeneous landscape. Using methods in functional data analysis, a multivariate regionalization approach was developed that effectively summarizes complex spatial-temporal patterns associated
with landscape disturbance, as well as recovery, which is often left unaccounted in studies quantifying patterns associated with disturbance. Next, the quasi-experimental framework afforded by a hunting moratorium was used to compare the influence of lethal (i.e., hunting) and non-lethal (i.e., anthropogenic disturbance) human-induced risk on antipredator behaviour of an apex predator, the grizzly bear. In support of the predation risk allocation hypothesis, male bears significantly decrease risky daytime behaviours by 122% during periods of high lethal human-induced risk. Rapid behavioural restoration occurred following the end of the hunt, characterized by diel bimodal movement patterns which may promote coexistence of large predators in human-dominated landscapes. A multi-scale approach using hierarchical Bayesian models, combined with
post hoc trend tests and change point detection, was developed to test the influence of landscape
disturbance and conditions on grizzly bear home range and movement selection over time. The results, representing the first longitudinal empirical analysis of grizzly bear habitat selection, revealed selection for habitat security at broad scales and for resource availability and habitat permeability at finer spatial scales, which has influenced potential landscape connectivity over time. Finally, combining approaches in movement ecology and conservation physiology, a body condition index was used to characterize how the physiological condition (i.e., internal state) of grizzly bears influences behavioral patterns due to costs and benefits associated with risk avoidance and resource acquisition. The results demonstrated individuals in poorer condition were more likely to engage in risky behaviour associated with anthropogenic disturbance, which highlights complex challenges for carnivore conservation and management of human-carnivore conflict. In summary, this dissertation contributes 1) a multivariate regionalization approach for quantifying spatial-temporal patterns of landscape disturbance and recovery applicable across diverse natural systems, 2) support for the growing theory that apex predators modify behavioural patterns to account for temporal overlap with lethal and non-lethal human-induced risk associated with humans, 3) an integrated approach for considering multi-scale spatial-temporal change in patterns of wildlife habitat selection and landscape connectivity associated with landscape change, 4) a cross-disciplinary framework for considering the impacts of the internal state on behavioural patterns and risk tolerance.
Table of Contents
Supervisory Committee ... ii Abstract ... iii Table of Contents ...v List of Tables ... ix List of Figures ... xi Acknowledgements ... xvi Dedication ... xviiCo-authorship Statement ... xviii
1. Chapter 1 – Introduction ...1
1.1 Introduction ...1
1.1.2 Approaches in movement ecology ...1
1.1.3 Habitat conditions and landscape disturbance ...5
1.1.4 Objectives ...6
1.1.5 Scope of analysis ...7
2. Chapter 2 – Characterizing spatial-temporal patterns of landscape disturbance and recovery in western Alberta, Canada using a functional data analysis approach and remotely sensed data ....10
2.1 Abstract ...10
2.2 Introduction ...11
2.3 Methods ...13
2.3.1 Study area ...13
2.3.2 Data...14
2.3.3 Disturbance pattern metrics ...16
2.3.4 Smoothing the disturbance pattern metrics ...16
2.3.5 Clustering disturbance regions ...17
2.4 Results ...19
2.5 Discussion ...25
2.6 Conclusion ...28
3. Chapter 3 – Behavioural restoration: relaxation of antipredator movement patterns in a threatened apex predator following a hunting moratorium ...29
3.1 Abstract ...29
3.2 Introduction ...30
3.3 Methods ...33
3.3.1 Study area ...33
3.3.2 Data ...35
3.3.2.1 Grizzly bear telemetry and biological data ...35
3.3.2.2 Landscape disturbance and resource availability ...37
3.3.2.3 Covariates ...38
3.3.3 Statistical analysis ...38
3.4 Results ...40
3.5 Discussion ...45
3.6 Conclusion ...48
4. Chapter 4 – A multi-scale assessment of spatial-temporal change in patterns of habitat selection, movement, and landscape connectivity in a threatened grizzly bear (Ursus arctos) population in Alberta, Canada ...49
4.1 Abstract ...49 4.2 Introduction ...50 4.3 Methods ...53 4.3.1 Study area ...53 4.3.2 Data...54 4.3.2.1 Movement data...54 4.3.2.2 Covariates ...55 4.3.3 Statistical analysis ...56
4.3.3.1 Seasonal home range models ...56
4.3.3.2 Movement models ...58
4.3.3.3 Quantifying spatial-temporal change in home range and movement selection ....60
4.3.3.4 Identifying and mapping potential movement corridors ...61
4.4 Results ...61
4.4.1 Home range and movement selection ...61
4.4.2 Change in spatial-temporal patterns of home range and movement selection ...65
4.5 Discussion ...71
4.6 Conclusion ...76
5. Chapter 5 – Integrating conservation physiology and movement ecology: body condition influences the movement behaviour of a threatened large omnivore ...77
5.1 Abstract ...77
5.2 Introduction ...78
5.3 Methods ...81
5.3.1 Study area ...81
5.3.2 Data...82
5.3.2.1 Movement data and biological metrics ...82
5.3.2.2 Habitat and landscape disturbance ...83
5.3.3 Statistical analysis ...84
5.3.3.1 Classifying grizzly bear behaviour ...84
5.3.3.2 Covariates ...85
5.3.3.3 Hierarchical Bayesian multinomial regression ...86
5.3.4 Results ...87 5.3.5 Discussion...94 5.3.6 Conclusion ...99 6. Chapter 6 – Conclusion ...100 6.1 Conclusion ...100 6.2 Key findings ...105 6.3 Future research ...106 Bibliography ...110 Appendix A ...162
A1.0 Disturbance pattern metrics ...162
A2.0 Functional data analysis (FDA) approaches ...163
A2.1 FDA curve fitting ...163
A2.2 Functional principal component analysis (FPCA) ...164
A2.3 Functional analysis of variance (FANOVA) ...165
A2.4 R code for FDA regionalization ...165
B1.0 Tables ...169
B2.0 Figures ...171
Appendix C ...178
C1.0 Tables ...178
List of Tables
Table 2.1 Summary statistics and variance explained by the functional principal components analysis (FPCA) for the disturbance pattern metric time series and curves. ...21 Table 2.2 Functional analysis of variance (FANOVA†) comparing the mean functional disturbance pattern metric†† curves by cluster and disturbance type. ...24 Table 3.1 Description of covariates used in hierarchical Bayesian regression models examining diurnal movement patterns of grizzly bears in Alberta, Canada. ...36 Table 3.2 Model selection for the hierarchical Bayesian regression models. WAIC is the Watanabe-Akaike information criterion and criterion and ΔWAIC the difference in WAIC between the candidate model and the top model. All models included a group level effect for the individual (bear). ...40 Table 3.3 Tukey’s honestly significant difference test of predicted movement estimates in habitat with high non-lethal human-induced risk associated with anthropogenic disturbance and resource availability by gender-age class and for parks and protected areas. Mean estimates are compared for the hunt (2001 – 2005) and post-hunt (2006 – 2014) periods for each season. Significant (p < 0.05) increasing movement estimates (+) and decreasing movement estimates (-) from the hunt to post-hunt period are indicated in bold†. ...42 Table 3.4 Posterior summaries for model coefficients in the final selected model. Parameters whose 95% credible intervals do not contain zero are indicated in bold. ...43 Table 4.1 Posterior summaries for coefficients, B, of the home range selection hierarchical logistic regression model. Parameters† with 95% credible intervals that did not contain zero were considered influential and are indicated in bold. ...62 Table 4.2 Posterior summaries for coefficients, B, of the movement selection hierarchical conditional logistic regression model. Parameters† with 95% credible intervals that did not contain zero were considered influential and are indicated in bold. ...63 Table 4.3 Summary statistics of the least-cost path analysis for the period 2001 – 2006 and 2007 – 2014, by grizzly bear habitat region and season. ...69 Table 5.1 Proportion of grizzly bear movements classified in the three behavioural states (foraging,
active foraging, and searching) identified using a hidden Markov model, by gender, gender-age
class, and season. ...88 Table 5.2 Posterior summaries for coefficients, B, from the hierarchical multinomial logistic regression model with the greatest support. Parameters with 95% credible intervals that did not contain zero were considered influential and are highlighted in bold. All parameters had 𝑅̂ values < 1.10. ...90 Table B1 Principal component analysis factor loadings (PC1 and PC2) for climate variables and NDVI at the home range and movement scales. ...169
Table B2 Results of the Kruskal-Wallis test comparing median values across regions (parks, core, secondary, and tertiary grizzly bear habitat) and season (hypophagia, early-hyperphagia, and late-hyperphagia) for the home range selection and movement probabilities. Significant values (p < 0.01) of the test statistic are indicated in bold. Refer to Appendix B Figure B8 for the median and range. ...170 Table C1 Model selection for the hierarchical Bayesian multinomial logistic regression models. WAIC is the Watanabe-Akaike information criterion and ΔWAIC the difference in WAIC between the candidate model and the top model. All models included a group level effect for the individual. ...178
List of Figures
Figure 2.1 The study area in western Alberta, Canada. Six provincial grizzly bear management units (BMU) are shown as well as the 223 watersheds which represented the landscape unit of analysis. A network of provincial and federal parks and protected areas (P&PA) exists mostly along the western portion of the study area. ...14 Figure 2.2 Disturbance pattern metric curves by watershed (n = 223, the curves are ordered by increasing latitude and longitude of the watersheds). The multivariate time series of proportion disturbance (Pd), probability disturbance adjacency (Pdd), mean disturbance patch area (Mdpa), and number disturbance patches (Ndp), are modelled as a unique curve for each watershed using B-spline basis functions with the amount of smoothing determined using generalized cross-validation. The disturbance pattern metric curves can be interpreted based on their amplitude (y-axis), phase variability (x-(y-axis), and the shape of the curve. Change in curve amplitude indicates periods of increasing cumulative disturbance and recovery (i.e., decreasing amplitude) in the watershed. ...17 Figure 2.3 Results of the functional principal components analysis (FPCA) for disturbance pattern metric curves of proportion disturbance (Pd) (A), probability of disturbance adjacency (Pdd) (B), mean disturbance patch size (Mdpa) (C), and number of disturbance patches (Ndp) (D). Either two or three functional principal components (FPC) were required to explain 90% of the variance in the watershed disturbance pattern metric curves. The upper plot of each panel maps the FPC with the highest absolute score for each watershed. The lower plot of each panel shows the mean of the fitted disturbance pattern metric curve (solid black line) and how the shape of the mean curve varies if the FPC curve (not shown) is added (+) or subtracted (-) from the mean curve. ...20 Figure 2.4 Bayesian Information Criterion (BIC) values from the fourteen multivariate Gaussian finite mixture models. Models are fitted using the eleven functional principal component (FPC) scores calculated from the proportion disturbance (Pd), probability of disturbance adjacency (Pdd), mean disturbance patch size (Mdpa), and number of disturbance patches (Ndp) curves, as well as the centered and scaled watershed area as inputs. The model with the greatest support (BIC = -2095.35) had an ellipsoidal distribution, with variable volume, and equal shape and orientation (VEE), resulting in eight disturbance pattern clusters. The fourteen models vary according to the volume, shape, and orientation of the covariance ellipsoid which are equal (E) or varying (V), with a spherical or axis-aligned orientation (I) (further details available in Scrucca et al. 2016). ...22 Figure 2.5 Mean curves by cluster for the proportion disturbance (Pd) (A), the probability of disturbance adjacency (Pdd) (B), the mean disturbance patch area (Mdpa) (C), and the number of disturbance patch (Ndp) (D) disturbance pattern metrics. The mean curves represent the average temporal trajectory of each metric for the watersheds belonging to the cluster and can be interpreted based on their amplitude (i.e., range along the y-axis) and phase (i.e., range along the x-axis) variability. Each mean curve is associated with the watersheds mapped by cluster membership (E). Parks and protected area boundaries are shown in green...23 Figure 2.6 Results of the functional analysis of variance (FANOVA). Mean curves of the proportion area disturbed (Pd), probability of adjacency (Pdd), mean disturbance patch area (Mdpa), and number of disturbance patches (Ndp) for the attributed disturbances (forest fires, forest harvest, non-stand replacing, roads, and oil and gas well-sites) by cluster are shown. ...25
Figure 3.1 The study area located in the Yellowhead and Grande Cache bear management areas (BMA) in west-central Alberta, Canada. A large network of federal and provincial parks and protected areas (P&PA) is found in the area. A total of 216 grizzly bear captures (121 individuals) occurred from 2001 – 2014. ...34 Figure 3.2 Predicted diurnal movement rates (m/hr-1) by gender-age class, season, and period, controlling for disturbance density and resource availability using the upper quantile of these parameters. The hunt and post-hunt season estimates correspond with the spring hunt (April 1 – May 31) for the management periods 2001 – 2005 and 2006 – 2014. The hunt/mate and post-hunt/mate seasons correspond with latter portion of the mating season (June) for the two management periods. ...41 Figure 3.3 Predicted movement estimates (m/hr-1) of male and female grizzly bears in association with road density (km/km2) for the pre-sunrise, daylight, and post-sunset diurnal periods. Movement estimates are shown for the hunt and post-hunt season (April 1 – May 31) corresponding to the hunt (2001 – 2005) and post-hunt (2006 – 2014) management periods. The hunt/mate and post-hunt/mate season corresponds with the latter portion of the mating season (June) for the same management periods. All other model parameters are held at their mean. The vertical dashed line represents a road density of 0.75 km/km2, which is positively associated with human-induced grizzly bear mortality (Boulanger and Stenhouse 2014). ...44 Figure 3.4 Predicted movement estimates (m/hr-1) of male and female grizzly bears in association with well-site density (km/km2) for the pre-sunrise, daylight, and post-sunset diurnal periods. Movement estimates are shown for the hunt and post-hunt season (April 1 – May 31) corresponding to the hunt (2001 – 2005) and post-hunt (2006 – 2014) management periods. The hunt/mate and post-hunt/mate season corresponds with the latter portion of the mating season (June) for the same management periods. All other model parameters are held at their mean. ...45 Figure 4.1 Scale-dependent differences in the direction and magnitude of home range (A – density) and movement (B – distance decay; C – edge density) selection of roads. ...64 Figure 4.2 Scale-dependent differences in the direction and magnitude of home range (A – density) and movement (B – edge density) selection of forest harvest blocks. ...65 Figure 4.3 Map of the classified modified Mann-Kendall trend test results for home range and movement probability by season (hypophagia, early-hyperphagia, and late-hyperphagia). Pixels are classified based on the direction of the z-score (positive, negative, or no-trend) and significance of the test statistic (p < 0.01). ...66 Figure 4.4 The proportion of parks, core, secondary, and tertiary regions classified based on the significance of temporal trends in predicted home range selection and movement probability values from 2001 – 2014 by season (hypophagia, early-hyperphagia, and late-hyperphagia). Pixels were classified using a modified Mann-Kendall trend test based on the direction of the trend (positive, negative, or no trend) and the significance of the test statistic (p < 0.01). ...67
Figure 4.5 Mean posterior predictions for the period 2001 – 2006 of home range selection and movement probability for the hypophagia, early-hyperphagia, and late-hyperphagia seasons. ...68 Figure 4.6 Mean posterior predictions for the period 2007 – 2014 of home range selection and movement probability for the hypophagia, early-hyperphagia, and late-hyperphagia seasons. ...68 Figure 4.7 Cumulative least-cost path network for the hypophagia, early-hyperphagia, and late-hyperphagia seasons for the period 2001 – 2006 and 2007 -2014. Least-cost paths were smoothed using a Gaussian kernel and represent the cumulative path count (square-root transformed for visualization purposes) for all years during the period. ...70 Figure 4.8 Density of least-cost paths crossings of major highways and primary industrial roads for the hypophagia, early-hyperphagia, and late-hyperphagia seasons, for the periods 2001 – 2006 and 2007 - 2014. Point locations were smoothed using a Gaussian kernel to identify focal densities of road crossings. ...71 Figure 5.1. Predicted probability of grizzly bear behavioural state (foraging, active foraging, and
searching) by time of day for each gender-age class and season. Gender-age classes are defined as
adults (> 5 years old), subadults (3 – 5 years), and females with dependent cubs. ...89 Figure 5.2 Predicted probability of grizzly bear behavioural state (foraging, active foraging, and
searching) by time of day, gender-age class, and body condition index (BCI). BCI values are
defined here using the lower (Low: BCI ≤ 0.33), interquartile (Moderate: BCI > 0.33 & < 1.33), and upper (High: BCI ≥ 1.33) range. ...92 Figure 5.3 Probability of grizzly bear behavioural state (foraging, active foraging, and searching) by time of day, body condition index (BCI), and edge density of anthropogenic disturbance features in profitable high-risk habitat. We defined profitable high-risk habitat as patches that occur in the upper quartile of anthropogenic disturbance edge density and resource availability based on NDVI and food. BCI classes (low, moderate, and high) represent the interquartile range of BCI values by gender-age class. All other model parameters are held at their mean. ...93 Figure 5.4 Probability of grizzly bear behavioural state (foraging, active foraging, and searching) by time of day, body condition index (BCI), and distance (decay) to anthropogenic disturbance features in profitable high-risk habitat. We defined profitable high-risk habitat as patches that occur in the upper quartile of distance (decay) to anthropogenic disturbance and resource availability based on NDVI and food. BCI classes (low, moderate, and high) represent the interquartile range of BCI values by gender-age class. All other model parameters are held at their mean. ...94 Figure B1 The study area in west-central Alberta, Canada includes the Yellowhead and Grande Cache grizzly bear management areas (BMA). A large network of parks and protected areas (P&PA), and grizzly bear habitat classified as core, secondary, and tertiary. Grizzly bear capture locations (n = 218) are shown in blue. ...171
Figure B2 Probability of home range selection associated with covariates in the hierarchical Bayesian logistic regression. Only parameters deemed influential (95% credible interval did not contain zero) are shown. See Table 4.1 for parameter estimates. All covariates are centered and scaled. ...172 Figure B3 Probability of movement selection associated with covariates in the hierarchical Bayesian conditional logistic regression. Only parameters deemed influential (95% credible intervals did not contain zero) are shown. See Table 4.2 for parameter estimates. All covariates are centered and scaled. ...173 Figure B4 Annual predicted probability of grizzly bear home range selection for the hypophagia (A), early-hyperphagia (B), and late-hyperphagia (C) seasons predicted from the hierarchical Bayesian logistic regression. ...174 Figure B5 Annual probability of grizzly bear movement selection for the hypophagia (A), early-hyperphagia (B), and late-early-hyperphagia (C) seasons predicted from the hierarchical Bayesian conditional logistic regression. ...175 Figure B6 Distribution of the maximum predicted probability of home range and movement selection for the years 2001 – 2014 by season (hypophagia, early-hyperphagia, and late-hyperphagia). The maps show home range selection has generally increased over time, while the maximum observed movement selection generally occurred in the mid-2000’s. ...176 Figure B7 Density of the posterior predictions for home range and movement probability by year (2001 – 2014) and season (hypophagia, early-hyperphagia, and late-hyperphagia). The plots show increasing high probability densities of home range selection over time, while movement selection densities have decreased. ...176 Figure B8 Box and whisker plots of probability values associated with home range (A) and movement (B) selection by region and season. The boxes represent the median, 25th and 75th percentiles, the lines represent 1.5 times the interquartile range, with outliers shown as circles exceeding this range. Mean probability values are represented by the filled point in the box. ..177 Figure B9 Annual distribution of change points in the probability values associated with home range and movement selection from 2001 – 2014 for the hypophagia, early-hyperphagia, and late-hyperphagia seasons. The greatest number of change points in the mean and variance of the probability values occurred circa-2006. ...177 Figure C1 Box-and-whisker plots of body condition index (BCI) values for 60 grizzly bears classified by gender and age. Adult bears represent those greater than five years of age. The boxes represent the median, 25th and 75th percentiles, the lines represent 1.5 times the interquartile range, with outliers shown as points exceeding this range. ...180 Figure C2 Probability of grizzly bear behavioural state (foraging, active foraging, and searching) by time of day, body condition index (BCI), and all anthropogenic disturbance in profitable low-risk habitat. We defined profitable low-low-risk habitat as patches with no anthropogenic disturbance
and high resource availability based on the upper quartile of NDVI and food. BCI classes (low, moderate, and high) represent the interquartile range of BCI values by gender-age class. All other model parameters are held at their mean. ...180
Acknowledgements
I would like to take this opportunity to thank my supervisors Dr. Trisalyn Nelson and Dr. Chris Darimont for their guidance and unwavering support over the course of my PhD. You were always willing to lend a hand, answer a question, or just have a chat, and for that I thank you. To my committee member Gordon Stenhouse, and his team at the Foothills Research Institute Grizzly Bear Program, thank you for the numerous insights into grizzly bear behaviour and for entrusting be with the bear data. It has been a real pleasure working with you. Thanks as well to my committee member Dr. Farouk Nathoo for all your hard work and statistical advice over the years.
To my colleagues in the SPAR lab, a sincere thank you for all the support, advice, and laughs. Jessica Fitterer, Robin Kite, Gillian Harvey, Shanley Thompson, Karen Laberee, Michael Branion-Calles, Ben Jestico, and Jed Long, it was a pleasure getting to know you all. Extra thanks to Dr. Jessica Fitterer for all the whiteboard sessions and advice over the last year of my PhD research.
Special thanks to Dr. Mike Wulder and his colleagues at the Pacific Forestry Center, Natural Resources Canada, for the advice, support, and data throughout my program. I’d like to also thank Dr. Nicholas Coops, and the members of the Integrated Remote Sensing Studio at UBC, for their remote sensing expertise and willing collaboration. Thanks as well to the ACS lab, Megan Adams, and Dr. Kyle Artelle, for discussing all things bear related and making me feel so welcome while I was finishing my research.
To my parents and brothers, thank you for providing the foundation, being so positive, and for all the love and support. Finally, I cannot begin to describe how much I owe to my amazing wife Marcy. Thank you for the long walks on the beach, handling the ups and downs, and for all the sacrifices you have made. None of this would have been possible without you, I love you.
Dedication
I dedicate my dissertation to my daughter, Hadley June Bourbonnais. Who knew someone so young could be so wise. I love you kid.
Co-authorship statement
Chapters 2 through 5 of this dissertation are manuscripts that were co-authored. The following outlines contributions of the doctoral candidate, and each of the authors. A reference representing the publication status of each chapter is provided.
Chapter 2
Bourbonnais, M.L., Nelson, T.A., Stenhouse, G.B., Wulder, M.A., White, J.C., Hobart, G.W., Hermosilla, T., Coops, N.C., Nathoo, F., Darimont, C. (2017). Characterizing spatial-temporal patterns of landscape disturbance and recovery in western Alberta, Canada using a functional data analysis approach and remotely sensed data. Ecological Informatics, 39, 140–150.
MB developed the concept of work, conducted the analysis, and prepared the manuscript for publication. TN, GS, MW, JW, GW, TH, NC, FN, & NC aided in the preparation of the manuscript with comments, edits, and advice on structure. MW, JW, GW, TH, & NC provided remotely sensed data.
Chapter 3
Bourbonnais, M.L., Nelson, T.A., Stenhouse, G.B., Nathoo, F., Darimont, C.T. (In Preparation). Behavioural restoration: relaxation of antipredator movement patterns in a threatened apex predator following a hunting moratorium.
MB developed the concept of work, conducted the analysis, and prepared the manuscript for publication. TN, GB, FN, & CD provided comments, edits, and advice on structure. GB provided grizzly bear telemetry data.
Chapter 4
Bourbonnais, M.L., Nelson, T.A., Stenhouse, G.B., Nathoo, F., Darimont, C.T. (In Preparation). A multi-scale assessment of spatial-temporal change in patterns of movement, habitat selection, and landscape connectivity in a threatened grizzly bear (Ursus arctos) population in Alberta, Canada.
MB developed the concept of work, conducted the analysis, and prepared the manuscript for publication. TN, GB, FN, & CD provided comments, edits, and advice on structure. GB provided grizzly bear telemetry data.
Chapter 5
Bourbonnais, M.L., Nelson, T.A., Stenhouse, G.B., Nathoo, F., Darimont, C.T. (In Preparation). Body condition influences the movement ecology of a threatened grizzly (Ursus arctos) population in Alberta, Canada.
MB developed the concept of work, conducted the analysis, and prepared the manuscript for publication. TN, GB, FN, & CD provided comments, edits, and advice on structure. GB provided grizzly bear telemetry data.
Chapter 1
Introduction
1.1 IntroductionOver the past decade, considerable progress has been made to unify the theoretical and methodological approaches used to study animal movement in the field of movement ecology (Holyoak et al. 2008, Nathan et al. 2008). Broadly, movement ecology focuses on understanding the mechanisms, decisions, timing, and spatial-temporal context driving animal movement and behaviour (Nathan et al. 2008). While the field of movement ecology recently formalized the study of animal movement and behaviour, it also has a long and rich history of theoretical development, as well as qualitative and quantitative inference in ecology (Lima and Zollner 1996, Turchin 1998). Animal movement and behaviour are central in many foundational theoretical models in ecology, including optimal foraging theory (Schoener 1971, Charnov 1976, Pyke 1984), the ideal free distribution (Fretwell and Lucas 1970), predator-prey systems (Huffaker 1958, Holling 1959), island biogeography (MacArthur and Wilson 1967), and metapopulation dynamics (Levins 1969, Hanski and Gilpin 1991). For example, the conceptual foundations of optimal foraging theory and the ideal free distribution are based in part on animal memory and spatial awareness of profitable habitat patches and their ability to locate and move among them (Charnov 1976, Pyke et al. 1977, Kacelnik et al. 1992, Farnsworth and Beecham 1999). While many early theoretical models in ecology incorporating animal movement were focused on population distributions (i.e., Eulerian approaches, Turchin 1998), a major paradigm shift in movement research recognized the importance of individual variability (i.e., Lagrangian approaches) for understanding movement and behaviour, as well as population distributions and dynamics (Nathan et al. 2008, Morales et al. 2010).
1.1.2 Approaches in movement ecology
Advances in technologies for monitoring individual animal movement using global positioning system (GPS) telemetry have been a catalyst for recent theoretical and methodological developments in movement ecology (Cagnacci et al. 2010, Hooten et al. 2017). The improved spatial-temporal resolution and accuracy of GPS telemetry has revealed previously unobserved movement patterns and behaviour of individuals and populations at multiple scales (Nathan et al.
2008, Kays et al. 2015). Quantifying spatial-temporal patterns of animal habitat selection, which reflects the behavioural decisions and choices of individuals, represents one of the earliest and most widely used applications of GPS telemetry data (Morris 2003, Hebblewhite and Haydon 2010). Building on early studies using temporally coarse very high frequency (VHF) radio-telemetry (Craighead and Craighead 1965), contemporary approaches for modelling habitat selection typically employ resource selection functions (RSF; Manly et al. 2002), to quantify habitat selection in a used versus available framework based on GPS locations and associated environmental conditions (Johnson et al. 2006). RSF models have broad applications in wildlife management and conservation as coefficients reflect behavioural choices of individuals and populations (Lele et al. 2013). Further, RSF model predictions are commonly employed to characterize the intensity or probability of selection of a habitat unit (Boyce and McDonald 1999, Boyce et al. 2002, Nielsen et al. 2002, Hebblewhite et al. 2005, McLoughlin et al. 2010), and to estimate landscape resistance to movement and connectivity (Chetkiewicz and Boyce 2009, Shafer et al. 2012, Squires et al. 2013, Elliot et al. 2014, Zeller et al. 2016). Ultimately, habitat selection, characterized using GPS telemetry, represents a fundamental mechanistic component of animal movement that influences individual fitness (Fretwell and Lucas 1970) and the distribution and abundance of populations (Boyce et al. 2016).
Scale-dependence in habitat selection is well-known (Johnson 1980, Mayor et al. 2009, DeCesare et al. 2012a, McGarigal et al. 2016), and statistical inference can help reveal variability in the behavioural response and distribution of animals in relation to scale-based differences in resource availability (Rettie and Messier 2000, Dussault et al. 2005, Boyce 2006), risk and predation (Johnson et al. 2002, Hebblewhite and Merrill 2007, 2009), as well as landscape disturbance and heterogeneity (DeCesare et al. 2012a, Northrup et al. 2016a, Zeller et al. 2017). However, defining scale-based availability in RSFs is contentious (Johnson et al. 2006, Aarts et al. 2008) and improper specification of availability and autocorrelation in the movement process can result in significant bias in regression parameters (Beyer et al. 2010, Northrup et al. 2013). To address potential bias, RSFs have been recast in the framework of spatial and spatial-temporal point process models (Hooten et al. 2013, 2017, Johnson et al. 2013), using hierarchical modelling frameworks (Duchesne et al. 2010, Hooten et al. 2016, Northrup et al. 2016b), and by defining availability based on spatial-temporal characteristics of the movement process through step
selection functions (Arthur et al. 1996, Fortin et al. 2005a, Forester et al. 2009, Thurfjell et al. 2014, Avgar et al. 2016).
Compared to habitat selection, quantifying animal movement requires explicit consideration of temporal dependence, as well as the spatial context and scale over which movement occurs (Hooten et al. 2017). Increasing spatial-temporal resolution of GPS telemetry data (Tomkiewicz et al. 2010), which characterizes the movement process as a near-continuous path (Laube 2014), represents a significant challenge as sequential spatial locations violate the assumptions of independence of many traditional estimators of animal movement and space use (Kie et al. 2010). As such, quantifying animal movement as a dynamic space-time process required the development of increasingly sophisticated and computationally intensive approaches (Schick et al. 2008, Long and Nelson 2012a, Patterson et al. 2017). Further, methodological approaches in movement ecology must also account for hierarchical scale-based complexity in animal movement behaviour over space and time. For example, at fine spatial (e.g., metres) and temporal (e.g., minutes to hours) scales, foraging (Fauchald and Tveraa 2003, Fortin et al. 2005b, Frair et al. 2005, Barraquand and Benhamou 2008, Fryxell et al. 2008, Avgar et al. 2013, Auger-Méthé et al. 2016), predation (Whittington et al. 2004, Hebblewhite and Merrill 2009, Merrill et al. 2010, Gurarie et al. 2011), resting/bedding (Franke et al. 2004, Ordiz et al. 2011, Hertel et al. 2017), and avoidance of disturbance features, such as roads (Forman and Alexander 1998, Beyer et al. 2014, Kite et al. 2016), represent movement responses to resources, risk, and habitat features. At broader spatial-temporal scales (e.g., tens to thousands of kilometres, and days to months), home range dynamics (Burt 1943, Börger et al. 2008, Long and Nelson 2012b, Sorensen et al. 2015), migration (Webster et al. 2002, Bauer and Hoye 2014, Le Corre et al. 2014, Scott et al. 2014), and dispersal (Bowler and Benton 2005, Doerr and Doerr 2005, Clobert et al. 2009, Baguette et al. 2014), connect populations and represent behavioural responses to resource distributions and landscape heterogeneity.
Generally, approaches for quantifying animal movement identify the timing of different behaviours and/or the spatial extent over which they occur. Path segmentation (Benhamou 2014, Gurarie et al. 2015, Edelhoff et al. 2016) and state-switching models (Morales et al. 2004, Jonsen et al. 2005, Patterson et al. 2008, Schick et al. 2008, Langrock et al. 2012) are among the most commonly used approaches for identifying the timing of fine-scale behaviours from movement
paths. Here, the movement path is divided into distinct modes based either on the changes in the mean and variance of the properties (e.g., velocity, direction) of the movement process (Fauchald and Tveraa 2003, Nams 2005, Barraquand and Benhamou 2008, Gurarie et al. 2009), or modelled as characteristic behavioural states based on the properties of defined statistical distributions or variants of a random walk (Morales et al. 2004, Bartumeus et al. 2005, Jonsen et al. 2005, Langrock et al. 2012, Beyer et al. 2013). At the home range scale, broadly defined as the area used by an individual to acquire resources and mate (Burt 1943), contemporary approaches in movement ecology quantify the area used (e.g., the utilization distribution) by an individual using time-geographic approaches based on velocity (Long and Nelson 2012b, 2015, Wall et al. 2014) or movement-based variants of kernels (Horne et al. 2007, Benhamou 2011, Downs et al. 2011, Kranstauber et al. 2012, Lyons et al. 2013).
While the breadth and scope of analytical approaches for quantifying animal movement has increased, there remains a need to characterize animal behaviours within the spatial context of the habitat in which they occur. As there is often a temporal disconnect between the resolution of the movement process and available data for characterizing environmental dynamics and landscape disturbance, opportunities still exist to develop a modelling framework for understanding animal movement responses to rapidly changing landscapes and resource availability (Baguette and Van Dyck 2007, Doherty and Driscoll 2018). As multi-annual GPS movement datasets are increasingly collected, new opportunities exist to quantify long-term change in movement and behaviour of individuals and populations. Contextualizing movement and behavioural responses to external habitat conditions also requires consideration of the internal state of the animal, which is difficult to quantify (Nathan et al. 2008, Holyoak et al. 2008, Martin et al. 2013). Further development of cross-disciplinary approaches holds considerable promise to help address this methodological gap, through integration of movement data and ancillary physiological and heath metrics. For example, wildlife health metrics from the emerging field of conservation physiology, which focuses on identifying mechanisms influencing cause-and-effect relationships between the physiological state of individuals and environmental conditions (Cooke et al. 2013, Coristine et al. 2014), can help reveal unobserved internal factors influencing behavioural responses to external conditions and environmental change (e.g., Creel et al. 2013, Bourbonnais et al. 2014a, Creel 2018).
1.1.3 Habitat conditions and landscape disturbance
Given the complexity of animal movement, placing it in the proper context requires ancillary environmental data that closely matches the spatial-temporal resolution of the behaviour of interest. For far ranging terrestrial mammals, matching the scale and extent of their movements often requires the use of remote sensing data products which vary in terms of their spatial and temporal resolution (Pettorelli et al. 2014, Neumann et al. 2015). A number of remotely sensed indices are available for characterizing spatially continuous vegetation biomass and primary productivity. For example, the normalized difference vegetation index (NDVI), which provides a measure of landscape productivity and vegetation health, is often used to quantify seasonal changes in foraging biomass and resource availability (Pettorelli et al. 2011). NDVI has been linked to foraging opportunities (Pettorelli et al. 2007, Hebblewhite et al. 2008, Boettiger et al. 2011), as well as seasonal range dynamics (Mueller et al. 2011, Pettorelli et al. 2011, Viana et al. 2018) and migrations (Boone et al. 2006, Bischof et al. 2012, Van Moorter et al. 2013, Wall et al. 2013). Similarly, the fraction of photosynthetic radiation (fPAR), a measure of energy absorbed by vegetation driving primary productivity and habitat dynamics (Coops et al. 2009a), has been linked to broad-scale seasonal range size and habitat use (Nilsen et al. 2005, Bourbonnais et al. 2014a). Due to correlations with food and resource availability, remotely sensed vegetation indices are found to positively influence the physiological condition of animals (Couturier et al. 2009, Bourbonnais et al. 2013, 2014a), as well as reproduction (Pettorelli et al. 2006), and survival (Pettorelli et al. 2007, Hurley et al. 2014).
Due to their ease of calculation, remotely sensed vegetation indices are commonly employed in animal movement analyses to reflect the spatial and temporal scales of animal movement processes. While derived remotely sensed data, including land cover classifications and landscape disturbance feature detection, represent advantages for animal movement modelling, they still lack the temporal resolution to match multi-annual GPS animal movement data. For example, nation-wide land cover classifications developed at fine spatial resolutions (Homer et al. 2007, Wulder et al. 2008a), as well as maps of anthropogenic disturbance and human footprint (Sanderson et al. 2002, Ellis et al. 2010) are often temporally static and represent multiple years of landscape data and patterns. While informative, landscape patterns are the result of dynamic processes of disturbance and recovery that may not be accurately represented with static maps
(Turner 2010, Pickell et al. 2015, 2016a). Recently, the availability of free Landsat imagery (Wulder et al. 2008b, 2016, Wulder and Coops 2014), combined with increased computing capacity (Azzari and Lobell 2017), has resulted in the development of land cover change detection and attribution approaches for classifying annual anthropogenic and natural disturbance at broad-scales (Hermosilla et al. 2015a, 2015b), which can be used to inform annual change in land cover over time (Hermosilla et al. 2018). Given global concerns regarding the impacts of landscape disturbance and human activities on animal movement patterns (Mueller et al. 2011, Hardesty-Moore et al. 2018, Tucker et al. 2018), the use of time series of landscape disturbance imagery provide opportunities to quantify how movement patterns of individuals and populations have responded to changing landscape patterns over time.
1.1.4 Objectives
Longitudinal studies in movement ecology that quantify change in movement behaviour of individuals and populations in the context of changing landscape conditions are limited (although see for example Le Corre et al. 2014), due in part to the cost of collecting GPS telemetry data (Hebblewhite and Haydon 2010) and limitations in the spatial-temporal resolution of ancillary environmental data. As such, the objective of my PhD research was to develop approaches for quantifying change in the spatial-temporal patterns of wildlife movement and their habitat. To do so, I use a long-term GPS telemetry database, combined with a novel high spatial resolution landscape disturbance time series, to detect and characterize change in the movement behaviour and habitat of grizzly bears (Ursus arctos) in Alberta, Canada.
Grizzly bears represent an ideal species for testing models and approaches in movement ecology. As a large terrestrial omnivore, grizzly bears have complex seasonal dietary requirements (Munro et al. 2006, Coogan et al. 2014, Erlenbach et al. 2014, Adams et al. 2017), and demonstrate considerable individual behavioural plasticity which drives movement and habitat selection patterns at multiple scales (Nielsen et al. 2002, 2010, Edwards et al. 2011, Sorensen et al. 2015). However, similar to many large carnivores in North America, the range and distribution of grizzly bears has been substantially reduced due to human persecution and habitat loss (Laliberte and Ripple 2004, Estes et al. 2011, Ripple et al. 2014). In Alberta, Canada, a moratorium on the legal hunt of grizzly bears was instituted in 2006 and the population was listed as threatened in 2010, with provincial population estimates of less than 700 individuals (Festa-Bianchet 2010). Spatial
patterns of mortality (Nielsen et al. 2004a, McLellan et al. 2012), habitat selection (Nielsen et al. 2002, 2006, 2010, Stewart et al. 2013, Sorensen et al. 2015), and population density (Boulanger et al. 2018), are strongly dependent on human activities and anthropogenic disturbance, as well as associated resources and food.
Less is known about the movement ecology of the Alberta grizzly bear population, and how it has changed over time. As such, approaches and results from my PhD will provide timely information regarding how the movement and behaviours of grizzly bears have responded to landscape change, as well as to conservation and management efforts. Further, understanding the movement ecology of threatened populations, and those that exist at the periphery of the geographic range of the species, can contribute important context for other populations experiencing demographic declines (Carroll et al. 2001, Cooke 2008).
1.1.5 Scope of analysis
Multivariate regionalization approaches can summarize complex spatial-temporal landscape patterns providing insight into the underlying mechanisms and drivers influencing observed patterns (Hargrove and Hoffman 2004, Long et al. 2010, Powers et al. 2012). In Chapter 2, I developed a novel regionalization approach using functional data analysis for quantifying spatial patterns of landscape disturbance and recovery over time. Using a time series of landscape disturbance data derived from Landsat for Alberta, Canada from 1985 - 2011 (J. C. White et al. 2014, Hermosilla et al. 2015a, 2015b), I applied functional principal components (Ramsay et al. 2009) to characterize temporal trends in spatial patterns of disturbance in watersheds throughout grizzly bear habitat. The scores from the functional principal components analysis are clustered using a Gaussian finite mixture model (Fraley and Raftery 2002), and I summarize the relative influence of dominant landscape disturbance types, including forest fires, non-stand replacing (e.g., insect infestation, drought), forest harvest, roads, and oil and gas well-sites, within the clusters. The resulting regionalization provides a useful characterization of how spatial patterns of landscape disturbance have changed through time in support of grizzly bear management and movement research, and has broad applications to characterize spatial-temporal landscape patterns in other systems.
Hunting moratoriums combined with long-term animal movement datasets present a unique quasi-experimental setting to test theory related to the changing impacts of risk on behaviour of animals (Lima and Dill 1990). In Chapter 3, I took a longitudinal approach to quantify the relative influence of factors related to external risk on diurnal patterns of animal movement. Using 216 grizzly bear movement trajectories collected from 121 individuals from 2001 – 2014, I tested the predation risk allocation hypothesis (Lima and Bednekoff 1999), which suggests that individuals will adjust antipredator behaviour based on their own needs and temporal variation in risk. For grizzly bears, lethal (i.e., hunting) and non-lethal (i.e., anthropogenic disturbance) human activities represent the dominant form of risk (Benn and Herrero 2002, Frid and Dill 2002, Nielsen et al. 2004a, McLellan et al. 2012), and many carnivores adjust behavioral patterns in response to human predation and fear (Ordiz et al. 2011, Darimont et al. 2015, Clinchy et al. 2016, Hertel et al. 2016a, Smith et al. 2017). Accounting for resource availability and non-lethal risk associated with anthropogenic disturbance, I showed that individuals adjust diurnal activity patterns when risk associated with hunting is highest during the spring from 2001 - 2005, with behavioural restoration following the end of hunting season and a moratorium on the legal hunt in Alberta beginning in 2006.
Long-term monitoring of individuals using GPS telemetry affords opportunities to characterize multi-annual patterns of behavioural change in wildlife populations. With Chapter 4, I developed a multi-scale approach for quantifying change in grizzly bear home range and movement behaviour, as well as landscape connectivity, in relation to changing landscape conditions from 2001 – 2014 (Mayor et al. 2009, DeCesare et al. 2012a, Mcgarigal et al. 2016). I fit home range and movement models in a hierarchical Bayesian framework to predict the probability of home range and movement selection (Rhodes et al. 2005, Squires et al. 2013, Thurfjell et al. 2014), equivalent to second and third-order selection proposed by Johnson (1980), over time. I used post-hoc tests to assess temporal trends and change points in grizzly bear home range and movement probability. I identified clear increasing trends in home range selection and decreasing trends in movement selection related to anthropogenic disturbance, such as roads and forest harvest blocks, and topography. Landscape connectivity, which I quantified using a least-cost path analysis informed by the home range and movement probability surfaces (Squires et al. 2013, Zeller et al. 2016), has increased over time, with an increasing number of potential movement pathways connecting grizzly bear populations in west-central Alberta (Proctor et al.
2010, Boulanger et al. 2018). The framework and results demonstrate the utility of considering scale-based differences in response to external drivers, and how temporal trends in behavioural probability can be assessed over time in support of conservation management.
Behavioural responses to external stimuli will depend not only on the nature of external conditions, but also on the internal state of the individual (Ran et al. 2008, Martin et al. 2013). In Chapter 5, I examined how animals respond to external risk, represented by anthropogenic disturbance, as a function of their physiological condition (Holyoak et al. 2008, Martin et al. 2013, Blecha et al. 2018). Characterizing the internal state of an individual in movement models is difficult (Holyoak et al. 2008) and often requires the use of additional technologies (e.g., biologging, Cooke et al. 2004) or estimation of the internal state (e.g., hunger, Blecha et al. 2018). Combining approaches from movement ecology and conservation physiology (Cooke et al. 2013, Tomlinson et al. 2018), I used the body condition index (Cattet et al. 2002) of individuals to test the influence of internal state on grizzly bear behaviours classified using a hidden Markov model (Langrock et al. 2012), accounting for resource availability, anthropogenic disturbance, land cover, and topography. Results demonstrated that the probability of an individual engaging in foraging behaviour in high-quality, but risky, habitat was inversely related to their body condition. The overall approach highlights the cross-disciplinary potential of movement ecology and conservation physiology, and the results have important implications for the conservation of the threatened grizzly bear population in Alberta, Canada.
Finally, with Chapter 6, I present concluding remarks and a summary of methodological contributions to the field of movement ecology, as well as how results contribute to understanding of grizzly bear behaviour and conservation. I discuss the importance, as well as the challenges, of quantifying change in spatial patterns of both habitat and movement over time. Finally, I highlight future opportunities for leveraging long-term GPS telemetry data for quantifying change in spatial-temporal patterns of animal movement.
Chapter 2
Characterizing spatial-temporal patterns of landscape disturbance
and recovery in western Alberta, Canada using a functional data
analysis approach and remotely sensed data
2.1 Abstract
Landscape regionalization approaches are frequently used to summarize and visualize complex spatial patterns, environmental factors, and disturbance regimes. However, landscapes are dynamic and contemporary regionalization approaches based on spatial patterns often do not account for the temporal component that may provide important insight on disturbance, recovery, and how ecological processes change through time. The objective of this research was to quantify spatial patterns of disturbance and recovery over time for use as inputs in a regionalization that characterizes unique spatial-temporal trajectories of disturbance in western Alberta, Canada. Cumulative spatial patterns of disturbance, representing the proportion, arrangement, size, and number of disturbances, and adjusted annually for spectral recovery, were quantified in 223 watersheds using a Landsat time series dataset where disturbance events are detected and classified annually from 1985 to 2011. Using a functional data analysis approach, disturbance patterns metrics were modelled as curves and scores from a functional principal components analysis were clustered using a Gaussian finite mixture model. The resulting eight watershed clusters were mapped with mean curves representing the temporal trajectory of disturbance. The cumulative mean disturbance pattern metric curves for each cluster showed considerable variability in disturbance amplitude which generally increased markedly in the mid-1990’s, while disturbance amplitude remained low in parks and protected areas. A comparison of mean curves by disturbance type (e.g., fires, harvest, non-stand replacing, roads, and well-sites) using a functional analysis of variance showed that anthropogenic disturbance contributed substantially to curve amplitude in all clusters, while curve amplitude of natural disturbances was generally low. These differences enable insights regarding how cumulative spatial disturbance patterns evolve through time on the landscape as a function of the type of disturbance and rates of recovery.
2.2 Introduction
Terrestrial ecosystems are subject to a range of natural and anthropogenic disturbances that influence landscape dynamics and heterogeneity. In North America, the frequency, extent, and severity of natural disturbances, including forest fires (Stocks et al. 2002, Flannigan et al. 2006, Bourbonnais et al. 2014b), insect infestation and disease (Volney and Hirsch 2005, Kurz et al. 2008), and environmental impacts such as wind and drought (Dale et al. 2001), has been increasing due to anthropogenic influences and climate change (Turner 2010). Similarly, anthropogenic activities and pressures on many terrestrial ecosystems are growing (Venter et al. 2016), and disturbances such as forest harvest (Masek et al. 2011, Pickell et al. 2014), road network development (Forman and Alexander 1998, Trombulak and Frissell 2001), and energy development and mining (White et al. 2011, Pickell et al. 2014, 2015), contribute substantially to land use change (Defries et al. 2004, Pickell et al. 2015) and landscape fragmentation (Wulder et al. 2008c, Pickell et al. 2016a). Cumulatively, the extent and severity of disturbance influences ecological processes, and is temporally dynamic given post-disturbance recovery, regeneration, and succession (Frazier et al. 2015, Bartels et al. 2016). As such, monitoring and quantifying how spatial patterns of landscape disturbance vary over time can highlight how ecological processes are influenced by disturbance and recovery, as well as forest and land management (Fraser et al. 2009, Wulder et al. 2009).
Methods for detecting and attributing disturbance from remotely sensed time series datasets provide opportunities to quantify spatial patterns and temporal dynamics of landscape disturbance, change, and recovery. While disturbance detection and land cover change have been quantified using a variety of remote sensing technologies, Landsat data have been used extensively due to the longevity and spatial resolution of the available data (Wulder et al. 2008b, Hansen and Loveland 2012). The 30m spatial resolution of Landsat data allows for the capture of disturbances at human scales; that is, landscape alterations that are the result of a given management or land use decision can be discerned over large areas in a systematic and repeatable fashion (Wulder et al. 2012).
Detection of land cover changes has long been a primary focus of methodological development (Jönsson and Eklundh 2004, Lu et al. 2004). Contemporary change detection approaches have extended bi-temporal scene comparisons (e.g., Coppin et al. 2004) to include dense time series of remotely sensed imagery allowing more detailed land cover change detection
and characterization (Kennedy et al. 2007, Huang et al. 2010, Hilker et al. 2011, Banskota et al. 2014, Hermosilla et al. 2015b). Chiefly, bi-temporal scene comparison allows for capture of binary change, while the capture of change using more than two dates allows for insights on disturbance aspects such as rates and persistence combined with directionality, as well as overall pre- and post-disturbance trends (Gillanders et al. 2008). In part, change detection approaches have been facilitated by the rapid development of image compositing techniques, resulting in gap-free time series of spectral reflectance values (Hilker et al. 2009, Roy et al. 2010, Griffiths et al. 2013, J. C. White et al. 2014, Hermosilla et al. 2015b). The temporal dimension of products generated from remotely sensed data can be leveraged to develop new hypotheses on disturbance recovery and land cover change (Roy et al. 2010, Hansen and Loveland 2012, Franklin et al. 2015, Hermosilla et al. 2015a, 2016, Pickell et al. 2016a).
Landscape regionalization approaches, where geographic entities are grouped based on common factors in order to simplify and explain complex landscape and environmental dynamics (Hargrove and Hoffman 2004), are common in ecology (e.g., Coops et al. 2009b, Fitterer et al. 2012, Thompson et al. 2016). While regionalization approaches have been developed to characterize spatial patterns of landscape disturbance (Long et al. 2010) and land use (Zurlini et al. 2007) using landscape pattern metrics (e.g., Turner 1990, 2005, Boots 2006), the temporal dynamics of disturbance and recovery are often left unaccounted which can influence the interpretation of resulting patterns (Gómez et al. 2011, 2015, Pickell et al. 2016b). Further, characterizing spatial-temporal patterns of landscape disturbance and recovery is difficult as multivariate time series data are frequently high-dimensional (Warren Liao 2005).
Methods in functional data analysis (FDA) are specifically designed to characterize data that are multivariate, temporally structured, and high-dimensional (Ramsay et al. 2009, Morris 2015). In the FDA framework, the fundamental object underlying observed data are functions, typically assumed smooth in some sense. Within this context, discrete time series observations are considered to arise through the regular sampling of a single smooth function (i.e., curve), rather than thought of as a realization from a multivariate distribution, representing the underlying process (Ramsay and Silverman 2005, Ramsay et al. 2009). Statistical inference then involves reconstruction of the underlying function based on the noisy time series data, and subsequently methods including clustering (Jacques and Preda 2014, Hitchcock and Greenwood 2015), principal components analysis (Shang 2014), analysis of variance (Cuevas et al. 2004, Ramsay et al. 2009),
and regression (Morris 2015) have been extended to the functional data setting. By modelling multivariate disturbance time series data using the FDA paradigm, landscape pattern metrics quantified annually can be characterized as functions (i.e., curves) representing disturbance as a temporally continuous process, rather than a discrete state (Frolking et al. 2009, Kennedy et al. 2010, Gómez et al. 2011, Frazier et al. 2015).
The goal of our research is to characterize disturbance as temporally dynamic, allowing us to quantify and map cumulative patterns of disturbance while simultaneously accounting for recovery. To do so, we develop an FDA regionalization of landscape disturbance in western Alberta, Canada from 1985 to 2011 using Landsat disturbance time series data (Hermosilla et al. 2015b, 2015a, 2016). The region represents a complex and shifting mosaic of natural ecosystems and land use strategies influenced by extensive cumulative natural and anthropogenic disturbance (White et al. 2011, Pickell et al. 2015). Characterizing landscape disturbance patterns as curves, our regionalization identifies unique temporal trajectories of cumulative disturbance patterns that represent underlying distributions and temporal trajectories of specific natural and anthropogenic disturbance types, including forest fires, forest harvest and roads. The regionalization provides an effective means for quantifying and visualizing complex spatial-temporal patterns of disturbance and recovery and highlights the utility of FDA approaches when working with multivariate remote sensing time series data.
2.3 Methods 2.3.1 Study area
The study area is approximately 158 000 km2 in western Alberta, Canada (Figure 2.1). Monitoring of landscape disturbance in the region is critical due to the overlap between on-going resource extraction activities and habitat of threatened and endangered species in the region including grizzly bears (Festa-Bianchet 2010) and mountain caribou (COSEWIC 2014). Watersheds (n = 223), were selected as the landscape unit of analysis for the regionalization and defined using heights of land along major watercourses (White et al. 2011), and as such correspond to fifth-level catchments. Watersheds are commonly selected as an environmentally relevant landscape unit for monitoring forest and land cover changes (Wulder et al. 2009), as well as monitoring environmental variability and habitat security (Noss et al. 2002). The study area is topographically complex, with elevations ranging from 450 m to 3500 m above sea level. Major
ecosystems include prairies to the east transitioning to coniferous and mixed-wood forests, and high montane sub-alpine and alpine landscapes in the west. Resource extraction activities in the region include forestry, oil and gas exploration, mining, and agriculture, which are all serviced by an extensive road network. A network of parks and protected areas, where resource extraction is limited, exists primarily in the western mountainous area. Important natural forest disturbances include forest fires (Gralewicz et al. 2012a), insect infestation (Safranyik et al. 2010), and non-stand replacing disturbances of variable magnitude (e.g., wind, drought, stress).
Figure 2.1. The study area in western Alberta, Canada. Six provincial grizzly bear management units (BMU) are shown as well as the 223 watersheds which represented the landscape unit of analysis. A network of provincial and federal parks and protected areas (P&PA) exists mostly along the western portion of the study area.
2.3.2 Data
Landscape disturbance from 1985 – 2011 was characterized by Hermosilla et al. (2015b, 2015a, 2016) in a recently developed disturbance time series data product developed for Canada
from Landsat image composites. As further details regarding the development of these data are available in (J. C. White et al. 2014) and Hermosilla et al. (2015b, 2015a, 2016), we will only briefly describe the methodology. Annual Landsat image composites were developed using a best available pixel (BAP) approach with August 1 as the target date (J. C. White et al. 2014, Hermosilla et al. 2015a). Landsat pixels incorporated in the annual image composites were scored based on their temporal proximity to August 1, the sensor type (i.e., Landsat TM or ETM+), distance to clouds and cloud shadows (Zhu and Woodcock 2012), and atmospheric opacity (the pixel scoring criteria are available in White et al. (2014a)). Pixels with the highest cumulative score are incorporated into the annual Landsat image composite. Image composites where further processed to remove noisy observations (due to unscreened clouds, smoke, and haze). Any pixels in the target annual image composite with missing spectral values were infilled using a break-point informed temporal segmentation approach using available reflectance values for a given pixel series (Hermosilla et al. 2015b).
Spectral change detection in the annual composite imagery was accomplished through a breakpoint and trend analysis using Normalized Burn Ratio (NBR) pixel time series values (Hermosilla et al. 2015b). The NBR of a pixel is calculated as NBR = (B4-B7)/(B4+B7) where B4 and B7 correspond to Landsat bands 4 (near-infrared) and 7 (short-wave infrared) respectively (Key and Benson 2006). Based on the spectral trend analysis, pixel-based disturbances are characterized by the change year and a suite of change metrics detailing the magnitude and persistence of the disturbance (for the full list of change metrics see Hermosilla et al. 2015b). The change metrics are used to classify the forest disturbances as either stand replacing (e.g., harvest, wildfire, roads) or non-stand replacing (e.g., insects, stress) (Hermosilla et al. 2015a), while additional information on persistence and nature of a given change can further inform the attribution. Overall for Canada, change detection accuracy was 89%, 92% of the disturbances were correctly classified, 89% of the changes were labelled to the correct year, and 98% within ±1 year (Hermosilla et al. 2016). For the purposes of our research, we used the year of change to quantify annual spatial patterns of disturbance in the regionalization. We then provide context for the regionalization using the disturbance attribution where detected disturbances are classified as either forest fires, forest harvest, roads, well-sites, or non-stand replacing (e.g., insect infestation, wind, drought, disease).
