Wildlife species may respond to industrial development with changes in distribution. However, discerning a response to development from differences in habitat selection is challenging. Since the early 1990s, migratory tundra Bathurst caribou Rangifer tarandus groenlandicus in the Canadian Arctic have been exposed to the construction and operation of two adjacent open-pit mines within the herd's summer range. We developed a statistical approach to directly estimate the zone of influence (area of reduced caribou occupancy) of the mines during mid-July-mid-October. We used caribou presence recorded during aerial surveys and locations of satellite-collared cow caribou as inputs to a model to account for patterns in habitat selection as well as mine activities. We then constrained the zone of influence curve to asymptote, such that the average distance from the mine complex where caribou habitat selection was not affected by the mine could be estimated. During the operation period for the two open-pit mines, we detected a 14-km zone of influence from the aerial survey data, and a weaker 11-km zone from the satellite-collar locations. Caribou were about four times more likely to select habitat at distances greater than the zone of influence compared to the two-mine complex, with a gradation of increasing selection up to the estimated zone of influence. Caribou are responding to industrial developments at greater distances than shown in other areas, possibly related to fine dust deposition from mine activities in open, tundra habitats. The methodology we developed provides a standardized approach to estimate the spatial impact of stressors on caribou or other wildlife species.
The impact of industrial development on wildlife is a frequent and worldwide concern, and this is especially true for long-distance migrants whose traditional routes can be threatened by industrial developments (Berger 2004). Those long-distance migrant species include migratory tundra caribou Rangifer tarandus groenlandicus, and shifts in caribou distribution in response to human activities have been the focus of much research (e.g. Nellemann et al. 2003, 2010, Vistnes & Nellemann 2008, Polfus et al. 2011). However, comparing findings among studies is complicated by different methodologies and scales of disturbance (Stankowich 2008). Differences in results from analyses of the same data sets can trigger controversy (Noel et al. 2004, Joly et al. 2006), detracting from effective conservation and mitigation for species that may be impacted by industrial development.
We became interested in measuring displacement of migratory tundra caribou when investigating the impact of mine development on the Bathurst caribou herd on the central Canadian tundra (i.e. Northwest Territories and Nunavut). Migratory tundra caribou is a gregarious and migratory ungulate with ecological similarities to other open habitat, gregarious ungulates in Africa and Asia which face industrial developments on their ranges (e.g. Mongolian gazelles Procapra gutturosa; Ito et al. 2005).
The Bathurst caribou herd has declined since 1996 from an estimated 349,046 (± 204,975 CI) to 31,897 (± 10,932) caribou in 2009 (Nishi et al. 2010, Boulanger et al. 2011). The decline adds urgency toward understanding the cumulative effects of industry, harvest and other stressors. From the early 1990s onward, Bathurst caribou have been exposed to a boom in mining exploration, which culminated in the construction of two open-pit mines and one underground diamond mine on the tundra range of the herd between 1996 and 2005. During and after environmental assessment hearings for the diamond mines, strong concerns were expressed about how the mines would affect caribou movements and distribution (Boulanger et al. 2004, Johnson et al. 2005). The distance at which caribou change their behaviour, habitat selection and distribution relative to disturbance, which we term the ‘zone of influence’, has implications for measuring the cumulative effects of industrial activities on wildlife, especially where there are multiple projects (Duinker & Greig 2007).
Previous estimates of the distances over which caribou were displaced from industrial disturbance were based on recording the frequencies of caribou occurrence relative to distance from the source of disturbance (e.g. Nellemann et al. 2000, Cameron et al. 2005, Joly et al. 2006), and fitting polynomial curves to distances of caribou from disturbances (Boulanger et al. 2004, Johnson et al. 2005, Golder Associates Ltd. 2008a,b). Estimates of displacement for the Bathurst caribou herd using satellite collar and aerial survey data using polynomial-based methods ranged from 17 to 130 km (Boulanger et al. 2004, Johnson et al. 2005, Golder Associates Ltd. 2008a,b). We suspected that the large differences in the apparent zone of influence were the effect of scale (ranges of distances considered in the analysis), and uncertainty in the exact distance due to the curvilinear nature of polynomial curves. A recent literature review and simulation study by Ficetola & Denoel (2009) demonstrated that the method used to detect threshold effect distances in ecology had substantial effect on estimates, which further supported our suspicion.
The primary objectives of our study was to derive a robust statistical measurement of the zone of influence and determine mechanisms for estimated zones of influence We adopted a piecewise regression method that fit the hypothesized asymptotic zone of influence threshold relationship by estimating the exact distance out to which mines affected caribou distribution while accounting for habitat selection. We furthered the piecewise methodology by fitting models that were robust to sample biases in our data sets and by estimating confidence limits on zones of influence and associated statistics. The scale of the zone of influence raises the question of possible mechanisms. Given the zone of influence, we explored dustfall as a potential mechanism that may contribute to the changed caribou distribution beyond behavioural responses to people, physical structures and vehicles. We suggest that this methodology is applicable to the estimation of spatial response of any wildlife species to stressors.
Material and methods
Study area
Our study area is approximately 300 km northeast of Yellowknife, Northwest Territories, Canada (Fig. 1) and is within the southern Arctic ecozone (Ecological Stratification Working Group 1996). The glaciated landscape has esker complexes, boulder moraines, ancient beach ridges and numerous lakes. Shrub communities of willow Salix spp., shrub birch Betula spp. and Labrador tea Ledum decumbens dominate areas with adequate soil development. Mats of lichens, mosses and low shrubs are found across exposed rocky and gravel sites.
The Bathurst herd of migratory tundra caribou annually moves hundreds of kilometres from wintering ranges below the treeline, to calving and summer range on the open tundra (Gunn et al. 2001). It is during the post-calving through summer seasons that the caribou are most likely to be in the vicinity of the mines, so we restricted our analyses to 15 July-15 October. During this period, Bathurst caribou occupy about 100,000 km2 of tundra with a high-use area (70% kernel) of about 53,000 km2 and a core (50% kernel) of about 33,000 km2, all based on satellite-telemetry of cows from 1996 to 2008 (Environment and Natural Resources, Government of Northwest Territories, unpubl. data).
We analyzed caribou distribution relative to two of the three existing diamond mines within the Northwest Territories (see Fig. 1): Ekati (BHP Billiton Canada Inc.) and Diavik (Diavik Diamond Mines Inc.). The main Ekati mine and Diavik are 30 km apart. Both mines are open-pit mines with accommodation complexes, ore-processing buildings and airstrips. Ekati has a separate camp and an open pit (Misery), which is connected by a 29 km all-weather road to the main Ekati site. The Misery camp and pit are 7 km from the Diavik mine, which is restricted to a large island in Lac de Gras. Because of the juxtaposition of the Ekati and Diavik operations, we modelled these mines as a combined unit. Caribou data were available for the mine sites from preconstruction (Diavik: 1996-1999) and construction (Ekati: 1996-1998, Diavik: 2000-2002) through to full operation (Diavik: 2003-2008, Ekati: 1998-2008), which allowed description of caribou distribution relative to mine areas across a range of mine footprints and activities. The total footprint of the Diavik and Ekati mines in 2008 were 9.7 km2 and 29.9 km2, respectively.
General approach
We used a multi-step approach to analyses, which is summarized in Table 1 and described in detail in the following paragraphs. We first organized the aerial survey and satellite collar data and developed methods to confront sampling biases inherent to each of the data sets. Secondly, we developed a base habitat model using logistic regression with habitat covariates used to predict caribou distribution based on caribou locations from aerial survey and satellite collar data. Thirdly, we used the base habitat model to iteratively estimate a zone of influence (area of reduced caribou occurrence) around the mine sites through the use of a piecewise regression procedure with distance from mine site as an additional predictor variable. Separate analyses were conducted for aerial survey and satellite collar data. Finally, we tested dustfall predictions as a possible mechanism for the observed zone of influence distances.
Caribou data sources
We used two sets of caribou location data. First, we obtained caribou sighting data from weekly helicopter surveys during July-October, which used systematically spaced strip transects (BHP Billiton 2009, Diavik Diamond Mines Inc. 2009; Appendix I). Transect route, spacing and width, study area size and frequency of data collection (collectively called survey design) varied within and between the mines. Most transects, however, were systematically placed at 4- or 8-km spacing to provide a coverage of 15-30 km radius study areas out from mine infrastructure and were flown at an altitude of 150 m and at a speed of 145-160 kph. The transect width was 600 m on both sides of the helicopter, and the transects were divided into 1-km long cells. For analysis we considered surveys in which > 1 cell had caribou present (> 0.2% relative occupancy/survey), resulting in 168 useable aerial surveys flown between 1998 and 2008 with a mean of 10.5 surveys/year (range: 8-18 surveys). For these surveys, the mean relative proportion of cells in which caribou were observed was 5.1% (SD = 6.4%; range: 0.3-41.0%).
The second data set of caribou locations was obtained from satellite transmitter collars fitted to adult cow caribou tracked from April 1996 to October 2008 (Gunn et al. 2001, Environment and Natural Resources, Government of Northwest Territories, unpubl. data). The number of collared caribou annually available for analysis that could encounter the mine sites ranged from four to 19 (see the section Treatment of satellite collar data below; see Appendix I). The satellite collars varied from transmitting every seven days beginning in 1996 to every five days beginning in 1998 with the addition of daily duty cycle for mid-July-mid-August beginning in 2002. We used 3,705 point locations during our period of interest (57.1% daily, 36.9% 5-day and 6.0% 7-day) from an annual average of 11.5 (± 1.25 SD) individual cows.
Habitat variables
We used two sets of data to describe habitat. First, we used vegetation classes and landform features from the Land Cover Map of Northern Canada (NLC; Olthof et al. 2009) and Earth Observation for Sustainable Development of Forests land cover classification (EOSD; available at: http://www.geobase.ca/geobase/en/data/landcover/index.html). Esker coverage was extracted from 1: 250,000 scale National Topographic Data Base maps (Natural Resources Canada; available at: http://geogratis.cgdi.gc.ca/geogratis/en/index.html). We used 12 habitat classes pooled between the NLC, EOSD and eskers coverage.
Second, we included data on plant phenology and productivity, which could influence caribou habitat use and movements (Russell et al. 1993). We used Normalized Difference Vegetation Index (NDVI) imagery to track plant phenology and productivity within the study area (Appendix II). We generated NDVI values for 10-day intervals for each year thus accounting for differences in seasonality for each year in the analysis.
Treatment of aerial survey data
We applied resource selection functions (RSF; Manly et al. 2002) to assess habitat selection using caribou locations from both aerial survey and satellite collar data. We treated the aerial survey observations as presence or absence of caribou rather than absolute abundance to minimize the effect of contagious behaviour and group size (Millspaugh et al. 1998). We compiled the observations of presence or absence into successive 1-km cells that were 1.2 km wide and calculated the proportion of habitat classes within each cell. We determined the distance from the mine site for all transect cells using the distance from the centroid of each transect cell to the nearest mine site centroid or mine road. When we added outlying components of the Ekati development (Misery and Fox pits), we used the distance to the nearest infrastructure component or associated road.
To address spatial autocorrelation, we used a generalized estimating equation model (GEE; Ziegler & Ulrike 1998) to estimate correlations and robust empirical standard errors between successive observations on the same transect line for the most supported base habitat model. We used an exchangeable correlation matrix structure to account for spatial autocorrelation. We used type 3 χ2 tests, which are less sensitive to order of parameters in the models, to test for significance (SAS Institute 2000). We used ROC curves to estimate the goodness-of-fit for how well a model predicts presence or absence through a range of probability cutpoints. A cutpoint was the probability level at which presence or absence was declared in each cell. ROC scores vary between 0.5 and 1. A score of 0.5 would correspond to a model with no predictive ability and a score of 1 would correspond to a model with perfect predictive ability. Models with scores of > 0.7 are considered to be of ‘useful’ predictive ability (Boyce et al. 2002). We used SAS PROC GENMOD or PROC LOGISTIC for all analyses (SAS Institute 2000).
The abundance of caribou varied annually and seasonally, which created variation in habitat selection. We therefore used the relative abundance of caribou in the survey area, as indexed by the number of cells in which caribou were detected relative to the number of cells sampled, as a ‘nuisance’ predictor variable. This minimized the influence of abundance on habitat selection by allowing the probability of habitat being selected to increase linearly with abundance.
We explored the effect of varying survey design between mines and over time by estimating the interaction of different designs (as a categorical variable) and the estimated zone of influence predictor variables. We assumed that any issue with detection of caribou during aerial surveys occurred evenly across all habitat classes. This was a reasonable assumption given that all surveys occurred in open tundra areas with minimal topography or vegetation to obscure caribou. We expressed all estimates as odds ratio given that detection probability could bias absolute estimates of occupancy of caribou around mine areas (MacKenzie 2006).
Treatment of satellite collar data
We defined habitat availability for the satellite collar data set based on each caribou location and estimated movement rates. We determined the proportion of habitat classes in a 1-km buffer radius (the maximum error of the satellite collar locations) around collar locations. Then we compared each buffered point with the buffered area around six random points that were within a circle around the previous location of the collared caribou. The circle was the ‘availability radius’ defined by the 95th percentile of the distance moved for caribou for the interval between successive point locations (Arthur et al. 1996, Johnson et al. 2005). We tested the satellite collar data for interactions between availability radius (duration between fixes; i.e. duty cycle), which determined the size of the buffer where available locations were placed, and habitat variables.
Caribou can select habitat at a finer scale than the scale reflected by the availability radius, as the radius depends on the time between successive telemetry fixes. For this reason, we considered the interaction of each habitat variable with the availability radius. This accounted for possible scale effects and allowed all the data to be simultaneously considered in a single analysis. We included locations from caribou which could encounter the mine sites (as indicated by the availability radius) at least once in a given year in the analysis. This filter excluded only 4.8% (seven of 145) of the caribou-year combinations from the analysis.
We compared caribou location points (used) and random points using a conditional logistic regression (Hosmer & Lemeshow 2000). The analysis defined each used and six accompanying random points as a cluster. This cluster centred each comparison on the habitat available to the caribou at the time when the location was taken. This approach avoided issues with pseudoreplication caused by pooling telemetry data from different individual caribou (Pendergast et al. 1996, Johnson et al. 2005). We used k-fold cross validation to test goodness-of-fit of the used-random satellite collar data (Boyce et al. 2002). For this analysis, we subdivided the data into training and testing data sets based on Huberty's rule of thumb (Huberty 1994). We then tested the goodness-of-fit of a model developed with the training data set with the testing data set. We estimated the Pearson correlation (Zar 1996) of successive RSF score bins with the frequency of used locations in each bin (adjusted for availability area of each bin). If the model fitted the data, then the RSF bin score and area-adjusted frequencies should be positively correlated (Boyce et al. 2002). We expressed all estimates from the logistic model as odds ratios given that absolute probability of presence cannot be estimated using used/availability analyses (Boyce et al. 2002).
Procedure for base model habitat variable selection
To build a base habitat model, we applied individual logistic regression analysis tests (as previously described for aerial survey and satellite collar data) to determine the statistical significance of individual habitat predictor variables (Hosmer & Lemeshow 2000). We pooled yearly data for this analysis. The general form of the model was:
response = habitat variable + habitat variable2 + habitat variable*movement rate + habitat variable*season + habitat variable*mean NDVI score + buffer scale*habitat variable (satellite collar analysis only).
The quadratic term (habitat variable2) tested for situations when stronger associations with habitat values were likely to occur in the midpoint of the habitat variable value as opposed to a linear relationship. The interaction between habitat variables and movement rate was tested for cases when a habitat was used transitionally, as indicated by a significant relationship between movement rate and the given habitat variable. We used the interactions among habitat variables and mean NDVI to test for seasonal selection of habitats and account for yearly differences in seasonality. We also tested for seasonal selection of habitats due to factors such as insect harassment and other factors not accounted for by NDVI by modelling the interaction of habitat type and season. Seasons were defined as early summer (15-18 July), mid-late summer (21 July-20 August), fall migration (25 August-5 October) and rut/late fall (6-14 October). We standardized habitat variables to allow easy interpretation of slope coefficients and to minimize issues with varying measurement scales.
We then added significant variables from univariate tests into a multivariate model in the same order as in the univariate model (i.e. linear habitat variable and habitat variable*movement rate). We evaluated the fit of individual terms by Type 3 χ2 tests and empirical standard error estimates (SAS Institute 2000). From this, we derived a base habitat model, which we then used to test for the zone of influence of mine sites.
Estimation of the zone of influence of mine areas
To test for zone of influence, we used logistic regression for analyses with caribou presence/absence as the response variable and the base habitat selection model with distance from mine as the additional predictor variable. The habitat selection model accounted for caribou distribution due to habitat selection, with a ‘zone of influence’ predictor variable (abbreviated to and symbolized as ZOI) and associated regression coefficient (βZOI). We used a procedure analogous to piecewise or segment regression to determine an optimal cutpoint (Hudson 1966). For example, when a 1.5-km distance was tested, all presence or used locations > 1.5 km were set to 1.5 km, regardless of how far out they were. By doing this, the odds ratio of selection relative to the mine site (as estimated by distance from mine*βZOI) was allowed to change linearly up to the hypothesized ZOIs, at which point it would asymptote and remain constant for distances greater than the ZOI (as estimated by ZOI*βZOI; Fig. 2). We assessed the overall fit of each sequential ZOI distance model by its log-likelihood. If fit was improved by the βZOI term, then the log-likelihood should increase to a maximum at the statistically most probable ZOI before decreasing at larger distances (see Fig. 2). If there was no ZOI, then the log-likelihood would remain constant across the range of distances. The distance at which the log-likelihood was maximized is, then, the estimate for the ZOI (i.e. the maximum distance where an influence of the mine on caribou distribution could be detected). In addition, the relative magnitude of the difference in habitat selection caused by the mine could be estimated by the odds ratio of habitat selection at the estimated ZOI (). The odds ratio in this case was the relative increase in habitat selection at distances further than the ZOI relative to habitat selection at mine sites.
The relative shape of the likelihood curve assessed the strength of the ZOI. For example, an irregular shaped likelihood curve, or a curve without a peak indicated that other spatial factors were influencing caribou selection relative to the mine (and which were not already accounted for in the base habitat model). We estimated the confidence intervals for the likelihood curve from the range of ZOI distances in which the log-likelihood was within 1.92 of the maximum likelihood ZOI (Hudson 1971, Hilborn & Mangel 1997).
We accounted for underestimation of standard errors of the odds ratio (ORZOI) due to repeated statistical tests (to determine the optimal cutpoint in Fig. 2) by correcting ORZOI estimates and accompanying confidence intervals using the shrinkage methods of Holländer et al. (2004). In addition, we adjusted P-values for statistical tests of βZOI using the methods of Lausen & Schumaker (1992).
We also analyzed the effect of temporal changes in mine activity by grouping years into phases of mine development. To retain sample size, we combined data for 1996-1999 (1998-1999 for the aerial survey analysis), 2000-2002 and 2003-2008 (when Ekati and Diavik were both in full operation). We accounted for the expanding footprints of the Ekati-Diavik mines by adding the Misery pit and road to the footprint in 2000 and the Fox Pit in 2003.
We did not base habitat models on pre-mine data which potentially could confound habitat selection with the effects of the mine. However, we suspected that our analysis was robust to this issue given the large area used to formulate base habitat models compared to areas affected by the mines. To check for any effects, we reevaluated the base model habitat coefficients after the ZOI term was added to assess the relative sensitivity of the base habitat model habitat coefficients to the estimated impact of mines on habitat selection.
Estimation of dustfall as a possible mechanism for the ZOI
Rescan (2006) applied a CALPUFF atmospheric transport and dispersion model for Ekati and Diavik to predict deposition of finer dust particles (total suspended particles, TSP; mean mass ∼ 10 μm in size). The model generated isopleths of predicted dust deposition based upon wind strength and direction, dust types and other atmospheric factors, and we interpolated the grid values between successive contours (20-5,000 kg/ha/year). A value of 0 was assumed to occur 5 km outside of the 20 kg/ha/year contour based on the average distance between the contours 20 and 50, and adjusted for the interval increment.
We entered predicted fine dustfall (TSP) as a covariate to the base model for the Ekati-Diavik area to generate predictions of the odds ratio of habitat selection relative to TSP levels. We then contrasted these results with ZOI predictions to test whether TSP might explain the larger ZOI distances that we estimated. We used the dustfall model predictions for the period when both mines were operational (i.e. during 2003-2008).
Results
Aerial survey analysis
For the base habitat model, significant predictors of habitat selection included relative caribou population abundance, eskers, low shrub habitat, tundra habitat, tundra*NDVI and water bodies (Appendix III). The base model fit the data with a ROC score of 0.793.
Using the base model with additional ZOI terms, the asymptote of the likelihood curve corresponded to an estimated ZOI of 14 km (CI = 12.0-15.5 km) for all years of data collection (1998-2008; Table 2 and Fig. 3). Odds ratios of the ZOI effect for the Ekati-Diavik mine sites suggested that caribou were 4.2 times (SE = 1.08, CI = 1.4-8.4) more likely to select habitat at distances > 14.0 km from the mine sites compared with the immediate mine site areas, with a gradient of increasing habitat selection between the mine site areas and the zone of influence (see Fig. 3). The ZOI model terms were significant for the pooled Ekati-Diavik complex (Z = 10.94, P < 0.001), and the overall fit of the model was adequate (ROC = 0.795). The significance of base habitat model terms did not change (at α = 0.1) with the addition of ZOI terms, suggesting that our base model variable selection procedure was robust to the effects of the ZOI on habitat variables.
Differences in survey design also affected ZOI estimates as suggested by a significant interaction of survey design and ZOI term (χ2 = 20.3, df = 2, P < 0.0001). We set all predictions to correspond to the aerial design in which both Ekati and Diavik were simultaneously surveyed under the assumption that this was the best data set to estimate the ZOI for the pooled mine complex.
When each phase of development was considered separately, the ZOI changed as mine footprint increased (see Table 2 and Fig. 4). In the initial time period (1998-1999: Ekati construction), a weak ZOI was evident at 4.0 km. In the middle period (2000-2002: Ekati operation and Diavik construction), no ZOI was evident, as indicated by a lack of peak in the likelihood curve. When both mines were in operation (2003-2008 with seven open pits in total), the ZOI was evident at 14.0 km (CI = 13.0-15.0 km) from the mine sites, which was the same as the pooled estimate but with tighter confidence limits (see Table 2 and Fig. 3).
Satellite collar analysis
For the satellite collar base habitat model, significant predictors of habitat selection included bedrock/boulder*season, forest, low shrub, tall shrub, tundra and water as predictors as well as interactions between some of the predictors with season and scale (see Appendix III). The base habitat model displayed an adequate fit to the data as determined by Pearson correlation of area-adjusted frequencies and ordinal odds ratio bins (ρ = 0.902, P < 0.0001).
The proportion of daily fixes for the satellite collar locations increased after 2001, which resulted in higher densities of used points during 2003-2008 (Fig. 5). Although the caribou satellite collar locations were fewer near the mine areas and then peaked from 25-50 km from the mines before decreasing at further distances, habitat influences such as lakes were affecting the distribution as well as the mine activities.
Analysis of the ZOI by mine development phase suggested that initially, there was no effect or attraction to mine areas, but avoidance in later years (Fig. 6 and Table 3). A ZOI of 23.0 km (CI = 19.0-35.0 km) was evident for the early period (1996-1999) of the Ekati-Diavik complex development; however, the odds ratio of the ZOI was considerably < 1, indicating attraction to the mine areas rather than avoidance. Inspection of the raw data revealed large aggregations of caribou near and at mine areas in August-September 1996 and July-August 1999 that likely caused the apparent attraction. A ZOI of 3.0 km (CI = 1.0-39.0) was evident for the middle period (2000-2002) with an odds ratio of 2.26 (CI = 1.32-225.70) suggesting avoidance; however, the confidence limits on the ZOI estimate were large. A tighter ZOI of 11.0 km (CI = 1.0-17.0 km) was evident for 2003-2008 when both mines were in operation, with an odds ratio of 3.22 (CI = 1.64-10.13) also suggesting avoidance of the mine areas. The significance of base habitat model terms did not change (at α = 0.1) with the addition of ZOI terms, again suggesting that our base model variable selection procedure was robust to the effects of the ZOI on habitat variables.
The precision of ZOI estimates and odds ratio estimates were generally lower for satellite collar data (see Table 3) than for aerial survey data (see Table 2). The years 2003-2008 had the highest sample size of collars (see Appendix I) and may be the best representation of the current ZOI of the Ekati-Diavik mine areas based on the collar data.
Predicted dustfall and the ZOI
The CALPUFF model predicted that TSP declined rapidly > 2 km from mine development (Fig. 7; Rescan 2006). Using aerial survey data, the log of TSP as a covariate for the base Ekati-Diavik habitat model was a significant predictor (χ2 = 117.1, df = 1, P < 0.0001); the resulting model had a ROC score of 0.795, which suggested predictive ability. Plots of predictions suggested a steep decline in the odds ratio of caribou occurrence at relatively low levels of TSP (Fig. 8). A similar analysis for the satellite collar data using only caribou locations that were within 50 km of the Ekati-Diavik mine complex indicated that the log of TSP was also a significant predictor (χ2 = 13.9, df = 1, P = 0.0002).
To explore dust fall as a mechanism for the estimated ZOI, we referenced the mean predicted TSP level at 14 km from the CALPUFF model (23.0 kg/ha/year, SD = 11.1; range: 0-43.0) which corresponded to the estimated ZOI from aerial survey analyses during 2003-2008 (see Table 2). For the aerial survey data, the model predicted odds ratios of 0.55 (CI = 0.49-0.62) at TSP levels of 23.0 kg/ha/year, meaning that a caribou was 0.55 times as likely to occur at areas with this TSP level compared to areas with a TSP level of 0. From the satellite collar analysis, a caribou was 0.63 (CI = 0.51-0.81) times as likely to choose areas with TSP levels of 23.0. At TSP levels of < 23.0, odds ratios quickly approached 1 suggesting that negative habitat selection decreased rapidly at lower TSP levels.
Discussion
Our analyses suggest that migratory tundra caribou respond to open-pit mining operations on summer range by reduced probability of occurrence in a ZOI of about 14 km (using the aerial survey data set). The reduced occurrence around the Ekati-Diavik mine complex was most evident during the operation phase of both mines and less evident during initial operation of Ekati and the construction of Diavik. Using the aerial survey data set, caribou were about four times more likely to select habitat at distances > 14 km from the mine complex, with a gradation of increasing selection between the mine site areas and the estimated ZOI (see Table 2 and Fig. 3).
The strengths of our analyses were that first we used two independent data sets (aerial surveys and satellite collars) that generated similar results. Second, our analyses accounted for patterns in habitat selection, as we tested the goodness-of-fit of the base habitat model without the ZOI variables. Third, we constrained the ZOI curve to an asymptote, such that the average distance from mine complex could be estimated along with confidence limits on the ZOI and the strength of the ZOI. Finally, our approach also allowed exploration of mechanistic factors for ZOI while controlling for other factors influencing caribou distribution.
Effect of scale on estimation of ZOI
We ran models to ensure that we had not confounded the different scales of habitat selection. For example, we ran a model using satellite collar data that extended up to 100 km from the Ekati-Diavik area and found that the log likelihoods initially peaked at the estimated mine ZOI (of ∼ 11 km), but then peaked again at larger distance from mine values (of ∼ 70 km) with odds ratios of < 1. For the larger distances, we estimated that the ZOI models the core of summer range, as also indicated by the highest used point densities (see Fig. 5), rather than the ZOI of the mine area.
We suspect that it is the scale of our analyses that allowed us to detect a larger ZOI than previously published response distances. Most regional studies reveal that Rangifer spp. reduce their use of areas within 1-10 km of development (Murphy & Curatolo 1987, Wolfe et al. 2000, Nellemann et al. 2001, Mahoney & Schaeffer 2002, Cameron et al. 2005, Joly et al. 2006, Weir et al. 2007, Vistnes & Nellemann 2008, Polfus et al. 2011; but see also Nellemann et al. 2010). However, our study addressed the effects of large open pit mines (of ∼ 40 km2 cumulative footprint), which is a different configuration of stimuli to caribou than, for example, a road, transmission line or a tourist lodge. The scale at which caribou are selecting habitat relative to the imposed scale of measurement is also likely a mechanistic factor in determining the extent of influence of mines (Vistnes & Nellemann 2008). The open tundra habitat likely allows caribou to respond at greater distances, however, other studies such as at the Prudhoe Bay oilfield were also on tundra post-calving ranges (Wolfe et al. 2000, Vistnes & Nellemann 2008).
Effect of statistical methods on estimation of ZOI
We assumed that the base habitat model accounted for spatial variation in habitat selection, and that the primary factor influencing habitat selection relative to mine sites was the effects of mines. Inspection of the likelihood plots and the associated odds ratios of βZOI assesses the overall adequacy of the ZOI model and the presence of other gradients or factors that confound ZOI estimates. Also, comparison of base model habitat selection coefficients with and without ZOI terms allows a test of the effect of the ZOI on base habitat selection coefficients.
Earlier analyses of the Bathurst herd using polynomial methods suggested larger ZOI around diamond mines (of ∼ 17-30 km, out to 130 km; Boulanger et al. 2004, Johnson et al. 2005, Golder Associates Ltd. 2008a,b). With the polynomial approach, however, other habitat selection gradients, which occur beyond the ZOI, can influence the overall shape of the curve. For example, satellite collar data indicate a steep gradient of habitat use evident at distances > 50 km from mines as indicated by declining point densities (see Fig. 5). A quadratic curve fit to these data would be influenced by both the gradient from mine ZOI but also the other gradients, which would cause the peak of the curve to be shifted to the middle of the gradient. A ZOI based on the peak of the quadratic curve would therefore be overestimated due to the influence of the other gradient. The main issue with polynomial regression is that it can only approximate a ZOI since the asymptote is not clearly defined, as demonstrated in a recent simulation comparison of the estimation of thresholds by polynomial and piecewise regression, which concluded that piecewise regression provided the most robust threshold estimate (Ficetola & Denoel 2009).
Effect of data type on estimation of ZOI
The aerial survey data provided the strongest analysis of ZOI. However, although less influenced by larger summer range selection gradients, these surveys were constrained by the extent of survey area. Our modelling assumed that the areas surveyed encompassed both the zone influenced by the mine and areas beyond the influence of the mine to allow an estimate of the asymptote of the ZOI curve. Even in the early years of the Ekati-Diavik monitoring, there was reasonable coverage out from development (of ∼ 22 km for Ekati-Diavik).
The aerial survey data were not corrected for sightability bias (Buckland et al. 2004), but we assumed this had little impact on the analyses, as we considered the relative change (OR) in habitat selection to estimate ZOI rather than estimating the probability of presence/occupancy within habitat classes near mines. We assumed that sightability was constant across all habitat classes; a reasonable assumption given that all surveys were conducted in open tundra habitat.
The satellite collars provided less precise estimates of ZOI, largely due to limited sample sizes (resulting in fewer data available for areas near the mines) and less frequent duty cycles for the early years of study. Although Vistnes & Nellemann (2008) recommended the use of satellite collars, we suggest that particular attention has to be paid to sample size. In contrast, stringently designed aerial surveys sample areas adjacent to mine sites are a more consistent indication of presence and absence of caribou relative to mine areas.
Mechanisms for ZOI
Responses of wildlife to human activities can be considered as analogous to responses to predation with associated trade-offs in energetic costs (Frid & Dill 2002). Predation risk is often structured in space and time, and prey will shift their distribution at different scales to accommodate predation risk (Tolon et al. 2009). However, the spatial and temporal scales of predation risk for caribou in the vicinity of open-pit mines are unknown as a contributing mechanism to the ZOI.
One factor that correlates to the scale of the response is dustfall. Although dustfall has been described for its effects on vegetation (Myers-Smith et al. 2006), little is known about the response of herbivores to dust on forage. The mines used an atmospheric transport model (CALPUFF; Rescan 2006) to predict TSP deposition rates in excess of 5,000 kg/ha/year (1,360 mg/m2/day) close to mine activity in summer. Deposition rates decreased rapidly with increasing distance from mine activities, and the dust constituents were identified in lichens; however, our analyses suggest that caribou avoid habitats with even low levels of predicted TSP. While caribou distribution around the immediate mine area may also be affected by non-dustfall sensory disturbance, we show that the larger 14-km ZOI for caribou did coincide with the predicted geographic scale of dustfall (see Figs. 7 and 8), suggesting that TSP may be a mechanism for reduced use by caribou of areas within the estimated ZOI.
Implications of our analysis
Our results demonstrate a quantifiable ZOI from open-pit diamond mines on caribou distribution. Our results suggest that researchers studying impacts of industrial development on caribou and other wildlife species should consider a larger range of scales than those caused by immediate behavioural responses to noise or other smaller-scale disturbances. Alternative larger-scale spatial impacts, such as dust deposition on forage, should be considered in addition to behavioural responses that have been the main focus of past ungulate studies (Stankowich 2008). We also suggest that interaction between spatial thresholds of responses to human activities have to be considered in relation to the spatial and temporal scales of predation risk (cf. Tolon et al. 2009).
The area of reduced caribou occurrence from the Ekati-Diavik mine complex is ∼ 6.7% of the core and ∼ 4.2% of the high use area of summer range of the Bathurst herd. Cumulative impacts from other sources of disturbance on the landscape (Johnson et al. 2005) could have wider implications for the herd, but our study was not designed to measure any effects from the reduced use of the summer range (cf. Nellemann et al. 2000, Cameron et al. 2005, Vistnes & Nellemann 2008, Polfus et al. 2011).
Our methods can be further applied to explore the effects of industrial disturbance on other wildlife species by allowing a robust estimate of displacement while accounting for variation in habitat selection and scale effects. We suggest that this standardized robust approach will allow improvement in monitoring and mitigation measures to manage the impact of mines and other developments on wildlife species.
Acknowledgements
the Department of Indian Affairs and Northern Development and Renewable Resources and Environment Program, provided funding for our analysis, and we thank D. Livingstone for his support over the years. Environment and Natural Resources (Government of Northwest Territories, Canada) provided the satellite collar data and we thank A. D'Hont and B. Fournier for their diligence. Mine survey data were generously provided by BHP Billiton Diamonds Inc. (D. Abernethy & A. Nichols) and Diavik Diamonds Mines Inc. (S. Bohnet & C. English). D. Panayi, C. Wood and J. Virgl (all Golder Associates Ltd.) and M. Wen (Rescan Environmental Services Ltd.) provided additional data. C. Johnson (University of Northern British Columbia) and D. Taylor (Taylor Geomatics) provided advice obtaining base data. We thank O. Gimenez and two anonymous reviewers for their helpful suggestions on the manuscript.
References
Appendices
Appendices
Appendix II: Habitat classification and NDVI indices used in the analysis of ZOI
We condensed habitat categories by blending two sources to provide complete coverage of the study areas, based on similarities in descriptions, low frequency of some types and logical assumptions about caribou biology (Appendix II, Table 1). We pooled habitat classes using the Land Cover Map of Northern Canada (NLC; Olthof et al. 2009) and Earth Observation for Sustainable Development of Forests (EOSD; available at: http://cfs.nrcan.gc.ca/subsite/eosd/mapping) land cover classification. The NLC classification coverage was generally north of the treeline and was given precedence where coverage from both products overlapped. Esker coverage was obtained from 1: 250,000 scale National Topographic Data Base maps (Natural Resources Canada; available at: http://geogratis.cgdi.gc.ca/geogratis/en/product/search.do?id=8147). We converted linear eskers into polygons with a standardized width of 100 m.
NDVI is related to the proportion of photosynthetically absorbed radiation and is calculated from atmospherically corrected reflectance from the visible and near infrared channels from Advanced Very High Resolution Radiometer (AVHRR) flown on NOAA-series satellites (James & Kalluri 1994). We used 1-km resolution NDVI amalgamated by 10-day composite periods for 1996-2006 (Latifovic et al. 2005) and calculated the mean values for each 1 × 1 km cell within the study area.
Appendix III: Results of base habitat models using aerial survey and satellite collar
Aerial surveys
For Ekati and Diavik, the univariate tests revealed linear relationships between caribou distribution and relative occupancy, esker, sedge wetland and the water predictor variables (Appendix III, Table 1). Quadratic relationships were suggested between low shrub, tundra and water predictor variables. In addition, an interaction between tundra and NDVI suggested a positive seasonal influence of the use of tundra.
Satellite collars
The base habitat model displayed adequate fit to the data as determined by Pearson correlation of area-adjusted frequencies and ordinal odds ratio bins (ρ = 0.902, P < 0.0001). The base habitat model analysis revealed linear or quadratic selection of forest, tall shrub, tundra and the water habitat variables (Appendix III, Table 2). Seasonal selection was evident for bedrock-boulder, low shrub, treeline herb, tundra and forest (interaction with NDVI) habitat categories. The selection of forest treeline herb and low shrub was also dependent on scale as determined by the availability buffer width and corresponding fix interval.
Figure 1.
Location of the Ekati, Diavik and Snap Lake diamond mines in the Canadian Arctic. The broken line polygon represents the area of high use (70% kernel) of the distribution of collared caribou during 15 July-15 October 1996-2008. The largest extent of the aerial survey study areas is also shown around each mine. The treeline represents the northern extent of continuous forests.
Figure 2.
Model used to estimate the ZOI and the magnitude of the ZOI. If a ZOI exists (grey area), habitat selection (as reflected by odds ratio of selection compared to the immediate mine area) should increase until the distance where the mine has no influence on selection. At this point, the model should best fit the data as indicated by the highest log-likelihood value. The slope of the increase in odds ratio is estimated by βZOI. At distances beyond the ZOI, the ZOI predictor variable was set constant (i.e. all distances < 10 km were set to 10 km), therefore creating an asymptote in the ZOI curve.
Figure 3.
Predicted change in odds ratio (— with confidence limits) and likelihood curve (- - - - -) as a function of distance from the pooled Ekati-Diavik mine complex as determined from aerial survey data during 1998-2008. The confidence limit for the ZOI estimate which is based upon change in the likelihood score is delineated by vertical dashed grey lines. Estimates are modeled upon the aerial survey design that flew both Ekati and Diavik mine sites in the same survey (2004-2005).
Figure 4.
Likelihood curves based on aerial survey data as a function of time periods for the Ekati-Diavik pooled mine complex analysis during 1998-2008.
Figure 5.
Satellite collar point densities used for the Ekati-Diavik mine complex during 1996-2008. The number of collared caribou was different for each time period, therefore each curve should be interpreted in terms of relative distribution rather than the actual densities of caribou near the mines.
Figure 6.
Likelihood curves based on satellite collar analysis for time periods of the Ekati-Diavik mine complex during 1996-2008.
Figure 7.
Mean total suspended particles (TSP) levels (kg/ha/year) as a function of distance from the Ekati-Diavik mine complex during 2003-2008. The mine complex included Misery Road and Fox Pit. Estimates are based on CALPUFF model predictions (Rescan 2006).
Figure 8.
Predicted odds ratio of caribou occurrence as a function of predicted total suspended particles (TSP) level (kg/ha/year) for the Ekati-Diavik mine complex area from aerial survey data (A) and satellite collar data (B). The odds ratios did not change appreciably at concentrations > 500 kg/ha/year. A reference line is shown at TSP levels of 23.0 kg/ha/year, which corresponds to the estimated ZOI distance of 14 km from the mines as determined from aerial survey data.
Table 1.
Summary of hierarchical steps used to estimate zone of influence (ZOI). Each item is covered in detail in the manuscript text.
Table 2.
ZOI estimates for Ekati-Diavik mine areas as a function of time period from aerial survey data. The ZOI estimate, confidence limit (CI), relative precision (CI divided by estimated ZOI), significance of ZOI model term (βZOI), goodness-of-fit (GOF; ROC score) and the magnitude of ZOI effect as described by the odds ratio (ORZOI) are given. ZOI estimates are based upon units of 0.5 km.
Table 3.
Summary of ZOI estimates for the Ekati-Diavik mine complex based on used/random analyses of satellite collar data. The ZOI estimate, confidence limit (CI), relative precision (CI divided by ZOI), significance of ZOI model term (χ2), goodness-of-fit (ρ) and the magnitude of ZOI effect as described by the odds ratio are given. ZOI estimates are based upon 0.5 km units.
Appendix I. Table 1.
Number of aerial surveys where caribou were observed in > 1 cell, and the number of collared caribou used for analysis. Satellite collar data include only caribou which had a mine area within their availability radius at least once in a given year.
Appendix II, Table 1.
Habitat associations used in base habitat models.
Appendix III, Table 1.
Base habitat model for aerial survey analysis for the Ekati and Diavik mine area aerial surveys. Standardized slope estimates are given for habitat variables (see Appendix I).
Appendix III, Table 2.
Base conditional logistic regression habitat model used to estimate ZOI from satellite collar data for the Snap Lake and combined Ekati-Diavik mine sites.
