Boreal Caribou Can Coexist with Natural but Not Industrial Disturbances

For species at risk, it is important that demographic models be consistent with our most recent knowledge because alternate model versions can have di ﬀ ering predictions for wildlife and natural resource management. To establish and maintain this consistency, we can compare predicted model values to current or past observations and demographic knowledge. When novel predictor information becomes available, testing for consistency between modeled and observed values ensures the best models are used for robust, evidence ‐ based, wildlife management. We combine novel in formation on the extent of historical disturbance regimes (industrial and ﬁ re) to an existing demo graphic model and predict historical and projected demographics of woodland caribou ( Rangifer tarandus caribou ). Exploring 6 simulation experiments across 5 populations in Alberta, Canada, we identify the relative importance of industrial disturbance, ﬁ re, and population density to observed population size and growth rate. We con ﬁ rm the onset of signi ﬁ cant declines across all 5 populations began approximately 30 years ago, demonstrate these declines have been consistent, and conclude they are more likely due to industrial disturbance from the oil and gas sector within contemporary pop ulation ranges than historical ﬁ re regimes. These ﬁ ndings reinforce recent research on the cause of woodland caribou declines. Testing for consistency between observations and models prescribed for species recovery is paramount for assessing the cause of declines, projecting population trends, and re ﬁ ning recovery strategies for e ﬀ ective wildlife management. We provide a novel simulation method for conducting these tests. © 2020 The Authors. The Journal of Wildlife Management published by Wiley Periodicals LLC on behalf of The Wildlife Society.

Demographic models are used to understand factors affecting population growth and have a long history in biology (McIntosh 1986, Bacaër 2011. These models reflect the current knowledge of species-habitat relationships and can incorporate environmental predictors of demographic change. The application of demographic models for endangered or at-risk species has been cautioned because of the poor quality of available demographic data, associated agents of decline, and difficulties in estimating variance in demographic rates (Chaudhary and Oli 2020). As a result, testing for consistency between models and observations for species at risk when new data become available is important.
In wildlife management, demographic models commonly have a dual purpose: they aid in inferring the processes behind changing species demographics, and in recommending objectives for species recovery plans. For species with large distributions, demographic models within recovery plans have the potential to affect policy and natural resource economies within multiple jurisdictions, where alternative model structures could result in very different predicted effects of management regimes (Beissinger andWestphal 1998, Brook et al. 2000). Many species at risk are in decline because of human-induced habitat use changes (Laliberte andRipple 2004, Ripple et al. 2014). Boreal populations of woodland caribou (Rangifer tarandus caribou), for example, are at the forefront of recovery priorities under Canada's Species at Risk Act (Government of Canada 2019). These populations have garnered national and international attention, with only 15 of 51 populations (29%) considered self-sustaining as of 2017 (Environment and Climate Change Canada 2019), and a distribution that spans 9 provincial and territorial jurisdictions overlapping Canada's most valuable natural resources: oil, natural gas, and forestry (Hebblewhite 2017, Fortin et al. 2020. Populations in Alberta, Canada, have been intensively studied and most are in steady decline (Bergerud 1974, Hervieux et al. 2013, Government of Alberta 2017.
Historically, the major factor influencing landscape change in Alberta was wildfire. Industrial disturbances in the form of oil and gas extraction and forestry now compound these natural disturbances (Fortin et al. 2020, Serrouya et al. 2020, resulting in an accumulation of disturbance effects (Tattersall et al. 2020). The recent, extensive proliferation of industrial disturbances is due to international demand for Canada's natural resources (Festa-Bianchet et al. 2011, Hebblewhite 2017, Fortin et al. 2020. Oil, gas, and forestry practices cause habitat fragmentation and vegetation change, which induce apparent competition between caribou, moose (Alces alces americanus), and white-tailed deer (Odocoileus virginianus) as prey for wolves (Canis lupus;Seip 1992;Wittmer et al. 2005;Serrouya et al. 2017Serrouya et al. , 2019Tattersall et al. 2020) and black bears (Ursus americanus; Latham et al. 2011, Leblond et al. 2016. Assessment of strategies for preventing caribou population declines in highly disturbed landscapes have concluded there is limited potential for long-term population persistence without intensive, and expensive, interventions (Johnson et al. 2019, Serrouya et al. 2019, Winder et al. 2020. Teasing apart the relative influence of predictors on typically small populations is hampered by challenges associated with environmental stochasticity (Beissinger and Westphal 1998). Resultant limitations can lead to reservations in employing models for species recovery at landscape scales (Sorensen et al. 2008;Schneider et al. 2010Schneider et al. , 2012Sleep and Loehle 2010). Uncertainty tactics similar to those employed by climate change denial campaigns have been used to question the disturbance-based demographic models for woodland caribou (Boan et al. 2018). Validating models that link demography to habitat condition (and extending their temporal window) should bolster confidence in recommendations to support species recovery.
Model validation involves testing for consistency between model predictions and existing knowledge (Mayer andButler 1993, Beissinger andWestphal 1998). In other words, does a model conform to what we have seen in a system of interest in the past and what we observe at present? Demographic models for boreal populations of woodland caribou have been developed by academic and government researchers (Sorensen et al. 2008, Johnson et al. 2020, and could be combined with knowledge of historical landscape disturbance regimes (i.e., predictors of demographic change), inferred knowledge of historical caribou population dynamics, and knowledge of today's caribou demographics. This allows biologists to test the reliability of caribou demographic models, determine if existing demographic models need to be updated to improve model accuracy, and draw conclusions about the sensitivity of population declines to industrial versus fire disturbance. Validation is best achieved by using simple, spatially explicit models, evaluating relative rather than absolute rates of extinction, examining all feasible scenarios, and using relevant time periods for projections (Beissinger and Westphal 1998).
Our objective was to test the consistency between predictions and observations for woodland caribou by combining an influential regional demographic model, the Sorensen model (Sorensen et al. 2008), with extensive data on Alberta's historical disturbance regime (fire and industrial disturbances). Sorensen et al. (2008) identified 2 predictors of lambda (λ), the rate of population change, for Alberta's boreal populations of woodland caribou: the percentage of the population's range within buffered (250 m) industrial disturbances and the percentage of the range burned within the last 50 years. We used 6 simulations of historical industrial and fire disturbance regimes on 5 woodland caribou populations in Alberta, and compared predicted and observed population size, growth rate, and persistence. If the Sorensen model accurately describes the effects of industrial disturbance and fire on population-specific growth rates, then we predicted population growth rates were approximately stable (λ~1.0) prior to the onset of industrial disturbances, and recent and intensive accumulation of industrial disturbances has a larger effect on predicted population growth rates than fire or historical population size, when compared to observed values (Vors et al. 2007, Environment Canada 2011, Johnson et al. 2020).

STUDY AREA
The study area is located in Alberta and includes the ranges of 5 of 12 extant woodland caribou populations: West Side Athabasca River, Little Smoky, Cold Lake Air Weapons Range, Caribou Mountains, and Red Earth (Fig. 1). These populations occupy both lowland and upland boreal forest comprised of black spruce (Picea mariana), white spruce (P. glauca), tamarack (Larix larcina), lodgepole pine (Pinus contorata), and jack pine (Pinus banksiana). The Albertan boreal forest occupies 381,046 km 2 , with a mean annual temperature of 0.5°C, and mean annual precipitation of 480 mm (Alberta Parks 2006). For more details about these populations, including the geography, seasonality, movements, and other land use please see the draft Provincial woodland caribou range plan (Government of Alberta 2017).

Simulating λ from the Sorensen Model
The Sorensen model showed that industrial disturbance ( IND % ) and the percentage of the range burned within the last 50 years ( BURN % ) were negatively associated with population-specific (i) λ for each time period (t), and their coefficients did not significantly differ: More recently, measures of industrial disturbances based on features buffered at only 250 m have been considered conservative; buffers of 250 m to 1 km have since been used in modeling caribou demographic parameters (Environment Canada 2011, Johnson et al. 2020). The Sorensen model has many of the same covariates, definitions, and data sources as more recent, and nationally derived, models (Environment Canada 2011, Johnson et al. 2020 providing strong support for a generally negative relationship between woodland caribou population growth and landscape disturbance, despite the relatively small scale and sample size of their investigation. To backcast and forecast λ from the Sorensen model, we required information on the 2 covariates of this equation through time: fire and industrial disturbance within each population's range. We obtained digital range maps for 5 of the 6 caribou populations modeled in Sorensen et al. (2008) from the Alberta Caribou Committee (2006). We defined the state of caribou habitat within each range using a simple forest age-class model subject to annual disturbance by fire, coupled with a dynamic model of the industrial footprint.
The age-class model tracks the proportional area in 1-year age classes (y; 0 to the maximum age Y ). At each time step t, the range i is subject to a proportional area burned (b i,t ). As explained below, these values are empirically determined for 1940-2006, and otherwise sampled from a range-specific distribution fit to the empirical data. We assumed that each forest age class was equally susceptible to burning ( Johnson et al. 2001) and that caribou population range boundaries did not change over the simulation period. The amount of forest that is zero years old in time t + 1, F i,t+1 (0), is simply the amount burned in time t in that population's range, b i,t . The amount of 1-year-old forest at time t + 1, F i,t+1 (1), will be the amount that was zero in the previous year minus the fraction that was burned in the current year, (1−b i,t ) and so on. The forest age-class model may thus be written: where Y is the maximum, accumulating, age class. The percentage area of a population's range i that is 50 years old or younger at time t is then calculated from the age structure to yield: Fire regime.-We characterized the natural disturbance regime within each caribou range using time-series of digital fire maps published by Alberta's Ministry of Agriculture and Forestry, Wildfire Management Branch. The maps show the boundaries of all known human and lightning-caused fires >12 ha, for 1940-2006. Provincial authorities constructed these maps from contemporary fire management records, archival paper maps, and post-fire aerial photography; they are considered to be nearly complete for areas under forest protection (C. Tymstra, Alberta Wildfire Management Branch, personal communication), which included our 5 populations. Since the early 1960s, lightning-caused fires have predominated in this region (Cumming 2001). We overlaid all fire maps with caribou population range boundaries. From this intersection, we estimated b i,t for each range.
For simulation, we sampled b i,t for population i at time t from a beta distribution b i,t~B eta(α i , β i ). We estimated these parameters using the method of moments: where x i ̅ is the mean and i 2 σ is the variance of the annual proportional area burned for population i (Table 1).  :table 1) originally calculated the percent industrial footprint for each range by buffering all industrial disturbances by 250 m, including linear features (seismic lines, roads, pipelines, and transmission corridors, oil and gas well sites), permanent clearings, and forestry cut-blocks. They counted overlapping buffers only once, and identified these industrial features from remotely sensed imagery taken from 2000. To estimate temporal trajectories over the period of record is challenging because time-series data are not available for most of the feature classes contributing to the index, except for oil and gas wells. For these, we used a comprehensive, proprietary database of energy sector infrastructure in western Canada provided by geoLOGIC Systems, Calgary, Alberta. The database records the location to the nearest quarter section and the date of completion for all exploratory and production oil and natural gas wells drilled in Alberta to the present day. At a spatial resolution of 10,000 ha, the density of drilled oil and gas wells is correlated with the densities of the linear feature classes incorporated in the Sorensen model (Cumming and Cartledge 2004), and we assumed this spatial correlation held throughout the duration of our simulations. We calculated the annual well density (w i,t ) within range i for years t = 1940, …, 2006 as the number of wells drilled within the population range up to year t, divided by the range area in km 2 . Using these values as surrogates for total industrial disturbance, we calculated: where %IND i are the herd level values reported by Sorensen et al. (2008). Our surrogates account for differences in development history (1940-2000) between ranges, and also allow us to extend the analysis to 2006 while accounting for development activity after 2000. We further assumed

Simulated Caribou Population Trajectories
We applied the Sorensen model to simulate caribou population trajectories over 1656-2076. Spin-up consisted of replicate runs over 200 years (1656-1855) prior to simulation. Our simulations are split into 3 stages: pre-industrial (1856-1939), industrial (1940-2006), and future . In the pre-industrial stage, we simulated annual area burned using variations of the fire regime model (eq. 5) to evaluate the probability of population persistence under a natural disturbance regime in the absence of industrial development. The industrial stage used the empirical annual areas burned and the estimated developmental trajectories (eq. 6) to determine if the Sorensen model is consistent with the inferred population declines and persistence up to 2006, given our 67-year reconstruction of the historical disturbance regime. Our future stage evaluated the probability of each population's persistence under an industrial disturbance and fire regime from 2006, with annual areas burned sampled from fitted beta distributions (eq. 5). For each simulation, we initialized the forest age structure within each population's range to uniform values of . From these covariate data, calculations, and forecasts, we predicted i t , λ from equation 1. We know of no estimates for woodland caribou population sizes prior to 1940. We assumed an initial mean caribou population density of ρ = 0.03 adult females/km 2 (0.06 caribou/km 2 ). This value is 1.5 times higher than a contemporary woodland caribou population with limited (2%) landscape disturbance (Mealy Mountain, 0.02 adult females/km 2 ; Environment and Climate Change Canada 2019), and 7.9 times higher than an Alberta population with high (87%) landscape disturbance (Cold Lake, up to 0.0038 adult females/km 2 ; Burgar et al. [2019], Environment and Climate Change Canada 2019). This value also closely aligns with the calculated stabilizing density for forest-dwelling caribou, compiled by Bergerud (1992). We calculated average range carrying capacities as (Table 1). We calculated population trajectories under a model of truncated exponential population growth as: where N i,t is the number of adult female caribou in population i at time step t. Equation 7 imposes a strict form of density dependence applying only to populations at or above carrying capacity. Scenario evaluation.-We used simulation scenarios to test the consistency of the Sorensen model to observed Table 1. Landscape condition, historical disturbance regime parameters, and assumed initial size for 5 caribou population ranges in Alberta, Canada. We included the herd range area (km 2 ), estimated initial population size, mean annual proportion of caribou range burned (x i ̅ ), standard deviation of proportion of range burned annually ( 2 σ ), estimates of beta distribution parameters of herd fire regimes (α and β; 1940 through 2006), and estimates of the proportion of a herd's range within 250 m of an industrial feature (%IND) in 1980 and 2006.  (Table 2). Our purpose was first to determine if the Sorensen model results in plausible outcomes given the assumptions specified as part of our model parameterization, and second to identify which parameters have the largest effect on caribou population demographics (size and growth rate).
In the first scenario (Regular), we address the compounded effect of fire and industrial activity using equations 1, 4, 5, and 6. In our second scenario (Decreased Density), we decreased the population carrying capacity to an estimate of contemporary woodland caribou population under limited landscape disturbance (Mealy Mountain, 0.02 females/km 2 ; Environment Canada 2012). The third scenario (+Fire) increased the pre-industrial (pre-1940) burn rate to 0.01 to reflect the effect of 2 possible factors creating a bias towards a low estimate of burn rate: fire suppression over the period of the empirical dataset (i.e., 1940-2006) and climate change into the future (i.e., post-2006; Boulanger et al. 2017). The observed mean annual burn rates x i ( ̅ ) were all <0.009 (Table 1), which may be very low; rates ≥1% (0.01) may be typical of this region in the long term (Larsen 1996). In the fourth scenario (++Fire), we increased the annual burn rate to 1.6% (0.016), a level equivalent to annual forest removal by the forestry industry within this region. Recommended forest harvest rates under a natural-disturbance paradigm are based on fire return intervals under the assumption that forest harvest provides a landscape disturbance similar to fire; current guidelines allow for forest harvest rotations every 50-60 years (or 0.016 of the forest each year; Hunter 1993). This scenario assumes all polygonal disturbances on the landscape are due to forest harvesting. Our fifth scenario (0Industry) separates the effects of industrial disturbance and fire by reducing industrial disturbance to zero and keeping all other parameters the same as in the Regular scenario. Our sixth scenario of no industry and very high fire (0Industry++Fire) eliminates linear industrial disturbance and assumes fire burn rates solely reflect contemporary forest harvest rotations, not a combination of existing fires and contemporary forest harvest rotations (Table 2).
Simulation scenarios ran for 220 years (pre-industrial through future stages: 1856-2076). Each simulation scenario consists of 300 replicate runs (resulting in 66,000 individual year estimates), repeated for each of the 5 caribou populations (Table 1). For each scenario, and each population, we calculated the average yearly λ and population size. We report the mean and 95% confidence intervals where this average yearly λ first fell below 1.0, the projected mean λ and population size in 2017 because this was the last reported observed value provided by the Alberta Government (Government of Alberta 2017), and the range of years where the mean population size fell below a quasi-extinction threshold of 10 adult females. We conducted all analyses in R version 3.3.5 (R Core Team 2019), and present values as mean ± standard errors unless otherwise specified. We obtained animal data from Government reports and websites; no animal care approval was required for this simulation study. Analysis and manuscript code are available at https://github.com/ StewartResearch/Boo2020.git.

RESULTS
We obtained natural and industrial disturbance data from 5 woodland caribou ranges across 67 years (1940 to 2006). On average 0.004 ± 0.0001 of each population's area burned annually ( Table 1). As of 2001, Sorensen et al. (2008) reported between 8.0% (West Side Athabasca) and 29.2% (Red Earth) of population ranges had burned, between Table 2. Parameter values for each stage of 6 simulation scenarios used to assess the reasonableness of the woodland caribou Sorensen demographic model. We changed parameter values between scenarios to test the effect of the proportion area disturbed by industrial disturbance (%IND), proportion area burned (%BURN), and population assumed density (max. density; adult females/km 2 ). We ran simulation scenarios for 5 woodland caribou populations (i) in Alberta, Canada, and across 3 stages, each representing a different historical period: pre-industrial (1856-1939), industrial (1940-2006), and future   Table 1). The historical disturbance data used here indicated 0% of all population ranges fell within a 250-m buffer of industrial disturbances as of 1940 and by 2006 this value increased on average to 77.1 ± 9.8% (Table 1). Our Regular simulation scenario combined historical disturbance data with the Sorensen model. During the preindustrial stage of this scenario, no population had a λ < 1.0 prior to the onset of industrial disturbance, the exception being Cold Lake whose λ may have approached 1.0 during some years because of its small population size ( Fig. 2; Lande 1993). Quasi-extinction never occurred across the other 4 populations, and population size stayed within 95% of K. These same results persisted across industrial and future stages when we set industrial disturbances to zero (0Industry and 0Industry++Fire); quasi-extinction never occurred, and population sizes stayed within 11% of K (Little Smoky; Fig. 2), even with much higher fire return intervals than occurred in our period with fire data (++Fire). Importantly, prior to incorporating historical disturbance regimes, the Sorensen model appeared to be consistent with the persistence of all 5 populations.
The industrial stage represents the start of data-estimated disturbance events. Four of our 6 scenarios had declining population growth (Fig. 2); only the 2 simulation scenarios where industrial disturbances were set to zero (0Industry and 0Industry++Fire) showed no declines. During this stage, all populations had scenarios where λ values fell below 1.0 (Table 3).
The future stage represents the projection of the Sorensen model from 2007 through to 2076. Across all simulation scenarios where we set industrial disturbance to zero (0Industry and 0Industry++Fire), populations remained above the quasi-extinction threshold during this stage ( Table 3). The only exception was the Cold Lake population, where quasi-extinction could have occurred as early as 1856, given stochasticity, in every scenario (Table 3); we attribute this exception to the high rates of pre-industrial fire disturbance (Tables 1 and 2) and small population size. All other scenarios indicated 2 of 5 populations were potentially quasi-extinct prior to 2017 (Little Smoky, Cold Lake; Table 3). This may be surprising, except that the projected λs for all populations in the Regular scenario were within the confidence limits reported by the latest Government of Alberta estimates for the same year (Government of Alberta 2017; Table 3); there was consistency between the Sorensen model tested with historical disturbance regimes and observed caribou populations in Alberta today. Increasing fire frequency to reflect climate change or fire suppression (+Fire scenario), current forest harvest rates (++Fire scenario), or decreasing density (Decreased Density scenario) generally decreased the Reg refers to our regular data scenario using historical disturbance regime data; −Dens refers to our Decreased Density scenario; +Fire and ++Fire refer to our scenarios where burn rate was increased to reflect either fire suppression, climate change, or annual forest removal; 0Ind refers to our zero industrial disturbance scenario; and 0Ind++Fire refers to our scenario with no industrial disturbance and increased burn rates that reflect contemporary forest harvesting regimes. average date (i.e., the date was earlier) where λ fell below 1.0, λ and population size as of 2017, and potential quasiextinction date for all populations, with the exception of Cold Lake. Compared to the Regular scenario, +Fire and ++Fire scenarios generally had a larger effect on these values than the Decreased Density scenario (Table 3). Given stochasticity all populations may have the potential to persist past 2076, but our simulations demonstrate these probabilities were significantly reduced by industrial disturbance ( Fig. 2; Table 3).

DISCUSSION
Our simulations indicate that all 5 Alberta caribou populations examined are able to persist under a historical natural disturbance regime of fire but not under recent cumulative effects of industrial disturbance. Woodland caribou populations were stable under simulations of a reduced carrying capacity or increased area burned by forest fires prior to the onset of industrial disturbances. This work provides a novel simulation method to test the predictability of demographic models relied upon for species recovery planning. We suggest, with similar historical data, our approach could be applied elsewhere to assess models relating demography to changing landscapes ( Johnson et al. 2020). We examined whether the Sorensen model predicted observed population demography when tested against historical disturbance regimes. We demonstrate significant consistency between observed and predicted caribou demographics in 2017 (Table 3), supporting inferences regarding causes of observed woodland caribou declines. The predicted λ in the Regular scenario was contained within the 95% confidence intervals of observed 10-year λ estimates (Table 3). These mean λ estimates correspond to the observed, widespread declines of Alberta's boreal woodland caribou populations over the last 30 years ( Fig. 2; McLoughlin et al. 2003, Hervieux et al. 2013, concomitant with rapid intensification of industrial development within population ranges (Hebblewhite 2017, Fortin 2020; the estimates do not match with our only fire scenarios (+Fire, ++Fire, or 0Industry++Fire), consistent with findings for Table 3. Predicted and observed demographics of 5 woodland caribou populations across Alberta, Canada. Across each simulation scenario (1856-2076) we report the first year that the average rate of population change (λ) fell below zero, the projected population size (N 2017 ) and λ in the year 2017 (λ 2017 ), and the first date of projected quasi-extinction. Quasi-extinction occurs when populations contain <10 adult females; blank cells represent extant populations as of 2076. Reg refers to our regular data scenario using historical disturbance regime data; −Dens refers to our Decreased Density scenario; +Fire and ++Fire refer to our scenarios where burn rate was increased to reflect either fire suppression, climate change, or annual forest removal; 0Ind refers to our zero industrial disturbance scenario; and 0Ind++Fire refers to our scenario with no industrial disturbance and increased burn rates that reflect contemporary forest harvesting regimes. woodland caribou populations in other parts of Canada's boreal forests (Schaefer 2003, Johnson et al. 2020.

Predictions from the Sorensen Model
We developed our predictions to evaluate the relative contributions of natural and industrial disturbances to the declines of Alberta's boreal populations of woodland caribou. Our results support our predictions; across populations, demographic parameters (population size and λ) were higher when we removed the effects of industrial disturbance. Testing the effect of our assumed carrying capacity by reducing it 1.5 times, and increasing the annual rate of disturbance by fire to 1.6% (up to a 33-fold increase over historical rates, depending on the population) had little effect on projected demographics; the time to quasiextinction decreased by only 0-18 years. The outcomes of population extinction or persistence as of 2017 were affected most by variations to the industrial disturbance regime (Table 3; Fig. 2). Emulating historical fire with contemporary harvesting is foundational to ecosystem-based forest management (Perera et al. 2007, but see Cyr et al. 2009 for an empirical critique). The natural range of variation (Landres et al. 1999) of forest fires in this region has been reported as up to 1.6% of the landscape (i.e., forest fires replace the forest every 60-100 yr; Hunter 1993, Andison 2013. We tested the effect of fireemulated harvesting on woodland caribou by increasing the annual proportion of area burned through simulation scenarios while keeping industrial disturbance consistent. If these values accurately represent the natural range of variation, the Sorensen model predicts the 5 study populations may become extinct between 1995 (Little Smoky) and 2035 (Caribou Mountains). Some caribou populations might coexist with a harvest return interval of 60 years (0Industry++Fire scenario), if there are no decreases in habitat carrying capacity, additional fires, or other industrial disturbances. This situation does not reflect the current state of Canada's boreal ecozone. Our simulations suggest that to sustain the boreal populations of woodland caribou, forest harvest rates should not exceed the historical fire regime of these Alberta populations, 0.9% annually (Table 1), or a harvest rotation of roughly 110 years.

Model Limitations and Edge Cases
Although the Sorensen model fit its data well (R 2 = 0.96), it used only 6 data points. Previous validation against 5 years of independent demographic data suggested it generally over-predicted λ (Sleep and Loehle 2010). Our initial values (Fig. 2, λ~1.2) may therefore be too high, leading to overestimated times to quasi-extinction (Table 3). Sleep and Loehle (2010) also question the equal weights assigned to industrial and natural disturbances in the model, corroborated by the more recent findings of Johnson et al. (2020) that herd demography is more sensitive to industrial disturbances than to fire. This would not alter our main results; indeed, it would further emphasize our conclusion pointing to the greater importance of industrial disturbances than fire ( Fig. 2; Table 3).
Our study neglected several potential sources of temporal variation in the data or in underlying processes. For example, our ability to estimate population densities may change with time (e.g., Government of Alberta [2017] vs. Burgar et al. [2019] density estimates for Cold Lake; Table 3); the effects of industrial disturbances may decrease over time as they return to a natural state, at rates that may vary by disturbance types (Dabros et al. 2018, Dhar et al. 2018; and true adult female mortality and juvenile recruitment rates may have varied with increased wolf densities because of landscape change and reduced trapping (Bergerud 1974). Neither did we account for wildlife management actions recently undertaken to conserve Alberta's caribou (see below); however, our model's simplicity supports rapid assessment of the relative effects of the factors that are included (Beissinger and Westphal 1998), and may be easily translated to other systems having spatially explicit disturbance histories.
Our model predicts very high pre-industrial population growth rates (λ~1.2) for some populations, under the scenario 0Industry (Fig. 2). This occurs when fire disturbance rates are very low. In that case, the predicted λ's approach the intercept of equation 1, which is 1.192. In population ranges with higher historical disturbance rates, or scenarios with elevated disturbance rates, these high growth rates are not predicted.
Our simulations project possible quasi-extinction for the Little Smoky population from the year 2013 onwards and for the Cold Lake population as early as 1856 (Table 3). Both populations may be vulnerable because of low numbers as derived from their small range areas (Table 1; Lande 1993). Although both populations have persisted to the present day, they did exhibit substantial and consequent population declines prior to 2001 (Government of Alberta 2017). The predicted and observed declines were ameliorated by wolf cull programs undertaken in winters 2005-2006 and 2016-2017 for Little Smoky and Cold Lake, respectively (Hervieux et al. 2014, Government of Alberta 2017). Our predictions for Little Smoky corroborate a recently published simulation that reported >68% of undisturbed habitat is needed to maintain a self-sustaining population over a 20-year period (Johnson et al. 2020); this condition has not been satisfied since 1981 according to our estimated development trajectory. With respect to the Cold Lake population, possible early quasi-extinction is compounded by the relatively high rate of fire disturbance on their range and limited cross-border movement due to province-specific altered landscape change with the contiguous SK2 herd in Saskatchewan, Canada (Government of Saskatchewan 2019). We think this points to an oversensitivity of the quasi-extinction index rather than any fundamental weakness in our models.
Predictions of demographic models for species at risk will be most accurate if reliable predictor histories are obtained (Chaudhary and Oli 2020), models can be decomposed into life-stage survivorship and recruitment terms (Johnson et al. 2020), assumptions regarding historical population density estimates are tested for consistency against observed data, and environmental stochasticity of survival and recruitment is incorporated into model predictions-all factors that are generally hard to come by (Williams et al. 2002, Lafontaine et al. 2017). In the present context, more detailed industry-specific disturbance histories would increase model sensitivity relative to the surrogate index of drilled well density that we used.

MANAGEMENT IMPLICATIONS
Our results suggest that many of Alberta's woodland caribou populations in boreal ecosystems will continue to decline as a result of cumulative industrial disturbances if preventative or mitigative management actions are not taken. Our data support the need to limit industrial disturbance within caribou population ranges because caribou may persist even under fire regimes with annual burn rates higher than recently observed if industrial disturbance is absent. In other words, caribou appear resilient to fire even at the most extreme end of our fire history uncertainty. Our results are also consistent with a positive effect of existing wolf culls (Little Smoky and Cold Lake) on population persistence; our simulations imply observed persistence would be unlikely in the absence of such treatments. This further suggests there remain opportunities to conserve other caribou study populations (e.g., Red Earth, West Side Athabasca, and Caribou Mountains) through management action that limits, reduces, or mitigates the presence of industrial features within their ranges.