Quantifying Long-Term Bird Population Responses to Simulated Harvest Plans and Cumulative Effects of Disturbance

There is interest in linking outputs from land use simulators to bird species distribution models to project how boreal birds will respond to cumulative effects of caribou (Rangifer tarandus) conservation, harvest, fire, and energy-sector development in Alberta. Our hypotheses were: (1) species associated with older mixed-wood stands would decline more if harvest was shifted away from areas used by caribou to areas with more mixed-wood; and (2) species associated with older forests would be more negatively affected by the combined effects of harvest, fire, and non-forestry footprint than by harvest alone. We used vegetation data from two harvest scenarios produced in Patchworks as inputs for density models of 20 boreal forest songbird and woodpecker species in Alberta. We projected abundance of these species over 50 years under: 1) two scenarios created in Patchworks, without fire but with and without deferral of timber harvest within a caribou conservation zone on lands tenured to Alberta-Pacific Forest Industries Inc.; (2) a scenario with fire but no human footprint; and (3) five scenarios in ALCES Online, in which habitat was affected by Patchworks harvest locations, fire (1–2 × current rate), and energy sector development (present or absent; with or without seismic line reclamation to improve caribou habitat). In the Patchworks scenarios, we found similar projected numbers of each bird species over time, whether harvest deferral occurred or not. Both harvest plans increased habitat and numbers for most species associated with older forests over 50 years, while most species associated with younger forests declined in both harvest plans, because average projected forest age increased over 50 years. Fire and other footprint generally reduced relative amount of habitat for species associated with older forests, which still increased over time, while other species responded positively or less negatively to fire. Seismic restoration created habitat for three-quarters of species that responded negatively to energy sector development over 50 years. As projections depended on whether just harvest, fire or all footprints were analyzed, multiple human impacts over time beyond harvest should be considered in conservation and land use planning based on long-term predictions about wildlife in anthropogenic landscapes.

There is interest in linking outputs from land use simulators to bird species distribution models to project how boreal birds will respond to cumulative effects of caribou (Rangifer tarandus) conservation, harvest, fire, and energy-sector development in Alberta. Our hypotheses were: (1) species associated with older mixed-wood stands would decline more if harvest was shifted away from areas used by caribou to areas with more mixedwood; and (2) species associated with older forests would be more negatively affected by the combined effects of harvest, fire, and non-forestry footprint than by harvest alone. We used vegetation data from two harvest scenarios produced in Patchworks as inputs for density models of 20 boreal forest songbird and woodpecker species in Alberta. We projected abundance of these species over 50 years under: 1) two scenarios created in Patchworks, without fire but with and without deferral of timber harvest within a caribou conservation zone on lands tenured to Alberta-Pacific Forest Industries Inc.; (2) a scenario with fire but no human footprint; and (3) five scenarios in ALCES Online, in which habitat was affected by Patchworks harvest locations, fire (1-2 × current rate), and energy sector development (present or absent; with or without seismic line reclamation to improve caribou habitat). In the Patchworks scenarios, we found similar projected numbers of each bird species over time, whether harvest deferral occurred or not. Both harvest plans increased habitat and numbers for most species associated with older forests over 50 years, while most species associated with younger forests declined in both harvest plans, because average projected forest age increased over 50 years. Fire and other footprint generally reduced relative amount of habitat for species associated with older forests, which still increased over time, while other species responded positively or less negatively to fire. Seismic restoration created habitat for three-quarters of species that responded negatively to energy sector development over

INTRODUCTION
Canada's boreal forest is continually being altered by human activities like forestry, energy sector development, agriculture, and climate change Gauthier et al., 2015). The cumulative effects of these activities are affecting the amount and suitability of habitat for wildlife (Schneider, 2019). Forestry operations are a dominant source of land-use change in Canada's boreal forest and generally shift the age-distribution of forests toward younger successional states (Kuuluvainen and Gauthier, 2018;Lavoie et al., 2019). Fragmentation by roads and other linear features (e.g., pipelines, seismic lines, power transmission lines) are also a concern if these features reduce forest patch size or increase edge effects (DeLong and Tanner, 1996;Dyer et al., 2001Dyer et al., , 2002Schneider, 2019). To reduce such effects on wildlife, many forestry companies have begun adjusting the spatial pattern, size distribution, and timing of harvests to better approximate natural disturbances like forest fire (Hobson and Schieck, 1999;Huggard et al., 2014). The goal of approximating natural disturbance is that it provides a coarse-filter approach that should be better at maintaining habitat for more wildlife species than traditional harvesting (DeLong and Tanner, 1996;Dzus et al., 2009;Kuuluvainen and Grenfell, 2012). However, such strategies come with economic costs, so it is important to assess the environmental benefits of different harvesting strategies. Future simulation tools that project outcomes for forestry yields, harvesting costs, and habitat quantity or quality for different species (e.g., caribou, boreal birds) are a crucial component of such evaluations (Sturtevant et al., 2007;Carlson et al., , 2019. Coarse-filter habitat management by forestry companies can be complicated by the needs of declining species. In western Canada, woodland caribou (Rangifer tarandus) have declined in many areas, leading to calls to defer logging in caribou zones (Dyer et al., 2001(Dyer et al., , 2002Wittmer et al., 2007). In the boreal plains ecoregion, harvesting of caribou habitat (black spruce dominated bogs and fens) is very uncommon as the trees are too small to be commercially valuable. However, within caribou zones, there are patches of upland, mesic, deciduous or white spruce stands, which are of low value to caribou per se, but are valuable commercial timber. Deferring or shifting logging away from upland patches adjacent to or within caribou zones has been proposed as a tool to minimize the risk that wolves (Canis lupus) and their primary prey (e.g., white-tailed deer [Odocoileus virginianus], moose [Alces alces]) are attracted to the early seral habitats created by forestry companies while also minimizing fragmentation caused by road construction (Latham et al., 2011). However, whether this fine-filter habitat management approach creates a conflict with other species at risk is not particularly clear (Villard et al., 1999;Drapeau et al., 2000;Imbeau et al., 2001).
While many qualitative assessments of caribou management on other species have been done, relatively few studies have used quantitative models to test the value of caribou as an umbrella species (however see Bichet et al., 2016;Drever et al., 2019). Previous studies on the value of caribou as an indicator species focus on the co-occurrence of other species with caribou using relatively coarse maps of species distribution, rather than projecting abundance of species based on detailed abundancehabitat relationships. Similarly, past work has tended to rely on relatively simple simulated landscapes or conservation networks rather than harvesting plans that will actually occur on the landscape. Declines in boreal birds may be exacerbated by a caribou-centric harvest strategy if forest harvest is concentrated in more contiguous older mesic upland forests, rather than in the small patches of upland interspersed in the lowland complexes preferred by caribou. To understand how long-term changes in land use to benefit one species affect other species, we can use spatial simulation modeling to create future landscapes. By linking spatial simulation models to species distribution models (SDMs) of abundance, it is possible to make predictions about how the size of bird populations may change under different harvest plans that vary in their objectives.
Quantitative predictions about bird response to forest harvest are only as realistic as the assumptions that go into each scenario. For example, in ecosystem-based management where harvests are designed to emulate natural disturbances, fire is usually assumed to be the most important disturbance influencing forest age and boreal forests can be stratified into areas naturally subject to and manageable with different fire or harvest frequencies ("ASIO model" in Angelstam, 1998;Kuuluvainen and Grenfell, 2012). Depending on the forest region, other disturbances or sources of tree mortality such as drought, wind storms or insect outbreaks may be more important disturbances in other regions (Seidl et al., 2017;Navarro et al., 2018). Forest harvest planning is often also based on the assumption that forest harvest is the sole disturbance setting forest age. Forestry companies design harvest strategies with the goal of approximating the distribution of forest fires, but an assumption underpinning these strategies is that fire control reduces the area burned sufficiently to allow for harvest to occur. It is becoming increasingly evident in western Canada that harvest and forest fires jointly affect forest age and the habitat available to different birds, despite best efforts at suppressing fire (Arienti et al., 2006). Furthermore, there is good evidence that bird abundance often differs between recently burned and recently harvested forests (Hobson and Schieck, 1999;Schieck and Song, 2006;Robertson and Hutto, 2007), particularly in the first 10-20 years post-disturbance (Schieck and Song, 2006;Huggard et al., 2014). Fires leave snags and patches of unburned vegetation and create temporary habitats for some species like woodpeckers and flycatchers (Schieck and Song, 2006) and some of these habitat attributes may be absent from some types of harvest blocks (Huggard et al., 2014). Whether or not there are significantly large differences in responses by species to fire and harvest in meaningful space and time remains a key question (Andison, 2003;Messier et al., 2003). SDMs that account for difference in bird response to harvest versus fire provide a much better way to assess whether forestry plans are maintaining birds within the natural range of variation (NRV) that would be expected when uncontrolled fire is the dominant disturbance agent. Simulating landscapes based on NRV in forest age and structure provides a way of estimating species abundance in the absence of human footprint, although the effects of a lack of harvest are confounded with the effects of increased fire due to a lack of fire suppression by humans.
While forest management plans often treat forestry as the only anthropogenic disturbance, there are an increasing number of other land-users in many areas of the western boreal forest. Assessing the impact of different sectors becomes a key priority in assessing overall risk to species and what the most effective management actions to conserve species might be. In western Canada's sedimentary basin, oil and gas development cumulatively deforest and/or alter vegetation structure of large areas each year (Brownsey and Rainer, 2009;Pickell et al., 2015). Climate change is also increasingly influential, with recent studies suggesting that boreal forests are seeing and will be subject to even higher rate of burning (Stralberg et al., 2015(Stralberg et al., , 2018, drought, wind storms, and insect outbreaks (Seidl et al., 2017;Navarro et al., 2018). Given that landscapes are being transformed by multiple agents of change, wildlife management decisions that focus solely on timber harvest may be ill-equipped to actually achieve management objectives if other disturbances are not accounted for. While integrated landscape management is recommended as a best practice, in reality it is rare in most jurisdictions (Kennett, 2006) as data from planned development by all relevant industries is usually unavailable. Decision support tools that simulate and incorporate multiple land uses (e.g., energy, forestry, agriculture, settlements, transportation, mining) and natural disturbances are needed . By linking SDMs to simulation tools, the consequences of management strategies in the presence of the full suite of drivers can be assessed (Carlson et al., , 2019. We linked detailed bird SDMs with several landscape simulation models in order to: (1) compare how populations of 20 boreal bird species of interest respond to forestry over the next 50 years with and without deferrals of harvest in caribou habitat; (2) predict the NRV in those species' populations in the absence of human footprint but with historical fire rates; and (3) predict cumulative effects to these species from forestry, other land uses (bitumen development, settlements), and fire (present burning rate and a doubled burning rate). We focused on species of conservation interest to Canada or Alberta or of management interest to foresters in Alberta. Thirteen of these species use older coniferous, mixed-wood, or deciduous forests as habitat. We predicted that species associated with older mixed-wood, or deciduous forests would decrease with deferral of harvest in caribou habitat, while species associated with older coniferous forests would increase with harvest deferral. We also expected bird species associated with older forests to respond negatively to forest disturbance by fire and non-forestry land use, and positively to restoration of land use footprint. Seven of our 20 species use younger or open habitats, and we predicted these species would respond positively to harvest, fire, and other land use.

Study Area
The study area was the Alberta-Pacific Forest Industries Inc. (Al-Pac) Forest Management Area (FMA) (∼6,300,000 ha) in northeastern Alberta, extending north from the towns of Athabasca and Lac La Biche to the Birch Mountains (∼340 km) and west from the Alberta/Saskatchewan border to Lesser Slave Lake. The predominant ecosystem is boreal forest, dominated by black spruce (Picea mariana) and tamarack (Larix laricina) in lowlands and by trembling aspen (Populus tremuloides), balsam poplar (P. balsamifera), and white spruce (Picea glauca) in mesic uplands. Another climax species, balsam fir (Abies balsamea), is uncommon due to the frequency of forest fires. Drier uplands are dominated by jack pine (Pinus banksiana). The merchantable land base comprises 31% of the FMA (Alberta-Pacific Forest Industries Inc, 2015). The FMA is divided into 12 Forest Management Units (FMUs) varying in size and the amount of forest stands that are suitable for harvest. A Forest Management Plan is produced every 10 years for the entire FMA but there are separate management targets for each FMU (Alberta-Pacific Forest Industries Inc, 2015, Figure 1).

Bird Species Distribution Models
The SDMs that we used for boreal birds were produced by the Boreal Avian Modelling Project (BAM), Alberta Biodiversity Monitoring Institute (ABMI), and Environment and Climate Change Canada (ECCC). The modeling process has been described in detail elsewhere and has been applied in other simulation and modeling studies (Ball et al., 2016;Sólymos et al., 2020b). Briefly, point count data were collated and standardized (n = 141,557 survey visits from 33,002 unique stations) from multiple boreal bird studies in Alberta's boreal forests . Boreal birds select habitat at multiple spatial scales ; thus, predictor variables in the SDMs were assessed at two spatial scales for each survey station. Local-scale variables were assessed in a 150-m radius of each station. Stand-scale variables were assessed in a 564m radius (1 km 2 ) of each survey station. This stand scale was chosen for pragmatic reasons to match the mapping unit in our predictions and because it roughly corresponds to the scale deemed most appropriate for landscape variables based on smoothing kernel estimates for landscape variables (Chandler and Hepinstall-Cymerman, 2016). At the local scale, land cover was assessed for each survey station using provincial land cover information (Alberta Biodiversity Monitoring Institute, 2017, 2018; Allen et al., 2019). Vegetation type included deciduous, mixed-wood, white spruce, pine, black spruce, tamarack fen, shrub, grass/herb, graminoid fen, marsh, and swamp cover types. Human footprint was assessed at each survey point based on the year of sampling. Footprint type included cultivation, forestry, urban-industrial (mines, well sites, urban areas, industrial, rural residential), hard linear (road and rails), and vegetated soft linear (seismic lines, pipelines, power lines, road verges) features. Proportional area of the land cover types was calculated at the local scale, and the dominant vegetation type was assigned to each survey station based on a simple majority rule. Various data sources (Alberta Biodiversity Monitoring Institute, 2017) were used to estimate the years since last disturbance (i.e., forest age) relative to year of sampling for birds. Age was calculated as the area-weighted average of the forested polygons within 150 m of survey stations. When the dominant land cover was a harvest block, the pre-disturbance vegetation type but not age was assumed based on available forest inventory data in the local 150-m buffer. Doing so treated harvested areas as young forest rather than a separate land cover type. We also created a contrast variable that ranged between 1 (harvest) and 0 (converged to natural stands) to describe the convergence trajectory of forestry cut blocks. We assumed that convergence is complete at 60 years after harvest. This allowed us to differentiate young forests of natural (i.e., fire) versus anthropogenic (i.e., timber harvest) origin. Stand-scale variables included: the amount of open water in a 1-km 2 buffer around each survey location; the proportions of total human footprint, vegetated footprint types, non-vegetated footprint types, linear footprint, non-linear footprint, cultivation, and non-cultivation footprint types; and the proportion of suitable habitat for each species ( Table 1). The suitable land cover classes were determined based on the binary classification of land cover types into suitable and unsuitable classes by maximizing the Youden index (Youden, 1950).
Geographic variation was captured by including latitude, longitude, and climate (mean annual precipitation, mean annual temperature, potential evapotranspiration, annual heat moisture index, frost-free period, mean warmest and coldest month temperature at 0.5 • resolution; Wang et al., 2012) ( Table 1).
The final SDMs were Poisson generalized linear models with a log link function. The response variable was the number of male birds of a species counted per survey. The QPAD approach was used to account for differences in sampling protocol and nuisance parameters affecting detectability (time of day, time of year, tree cover, habitat composition; Sólymos et al., 2013). This approach converts sampling distances and durations to a common standard through statistical offsets and adjusts for differences in detection error and sampling area related to broad vegetation types and timing of surveys. As a result, the results allow us to estimate density of birds (individuals per hectare), in a spatially explicit manner, allowing us to use forest stand type, fire or harvest origin, stand age, and human impact coefficients in scenario modeling.
Bird densities for the various scenarios were estimated based on the local stand, its spatial location, and the characteristics of the surrounding polygons. We modeled the effect of forest age on bird density by using weighted age and a quadratic or square root transformed terms as covariates to fit nonlinear responses to age. We incorporated interactions between forest type and age, climate variables, latitude, and longitude ( Table 1). The stand level effects of suitable habitats and human footprint allowed us to differentiate between locally suitable habitats that are surrounded by suitable vs. unsuitable land cover types. The stand-scale predictors effectively measure patch size in a species-specific manner, i.e., suitable habitat was assessed for each species individually, to best describe their optimal habitat characteristics.
We extracted habitat variables, and depending on the scenario, human footprint and climate variables, which served as inputs to the SDMs for predicting abundance of bird species at specific locations. For the Patchworks harvest scenarios and NRV scenario, these locations were individual quarter-sections throughout the Al-Pac FMA and we extracted proportions of different forest-origin type-ageclasses per quarter-section. The quarter-section IDs were linked to latitude and longitude locations stored with climate variables and interaction terms in the cure4insect R package (R Core Team, 2020; Sólymos et al., 2020a). For the ALCES Online cumulative effects scenarios, we extracted variables from 200-m square raster cells including the proportion of each cell occupied by a specific vegetation-origin type or human footprint, the proportion of land within 1 km 2 of each cell in different vegetation-origin types or human footprints, the weighted-average forest age in that cell, latitude, longitude, climate variables, and any interaction terms among these variables. Instead of importing these variables as inputs to SDMs within the cure4insect package, we constructed individual species indicators and ran each species indicator through the cumulative effects scenarios in ALCES Online (Table 1) using a set of parameter estimates from models that included the 1 km 2 scale predictor variables as described in Ball et al. (2016). Model coefficients and indicator formulae for species are stored online at https: //github.com/borealbirds/ABMI-bird-models-ALCES-Online.
We focused on 20 bird species in three key groups: (1) federally threatened species in Canada 1 ; (2) provincially threatened in Alberta (Alberta Environment and Sustainable Resource Development, 2014), or (3) have specific habitat requirements that make them useful indicator species for forest managers in Alberta. Most of these species are associated with mature or older deciduous or mixedwood forests (

Harvest Planning for Caribou Management
Many forest planners in Alberta create their operating harvesting plans using the harvest-scheduling simulator software Patchworks 2 . Patchworks is used because it optimizes economic and ecological tradeoffs when selecting which stands to harvest, in what areas, and when (Sturtevant et al., 2007). We used Patchworks to create two spatial harvest plans for the Alberta-Pacific FMA in northeastern Alberta that are being submitted as options for harvest to provincial regulators.
The first Patchworks scenario solves for different tradeoffs to maximize harvesting pulpwood and timber with relatively even harvest levels over time, while (1) minimizing costs of road construction and maintenance; (2) adjusting harvest area size, shape, and distribution to approximate the size, shape and distribution of natural forest disturbances like fires; and (3) applying other constraints (Hebert et al., 2003;Dzus et al., 2009). We describe this scenario as the Ecosystem-Based Management/Natural Disturbance Model ("EBM") scenario. 1 | Predictors used in the species distribution models to predict density and abundance of each analyzed bird species within individual Forest Management Units and over the Al-Pac Forest Management Area (FMA) over 50 years within the two Patchworks spatial harvest scenarios, the NRV scenario, and the ALCES Online cumulative effects scenarios described in this paper. For the Patchworks harvest scenarios and NRV scenario, we extracted the proportions of different forest-origin type-age classes from each quarter-section in the Al-Pac FMA, since quarter-section IDs were linked to latitude and longitude locations stored with climate variables and interaction terms in the cure4insect package. For the ALCES Online cumulative effects scenarios, we extracted the proportion of each 200-m square cell occupied by a specific vegetation-origin type or human footprint, the proportion of land within 1 km of each cell in different vegetation-origin types or human footprints, the weighted-averaged forest age in that cell, latitude, longitude, climate variables, and any interaction terms among these variables. We constructed individual species indicators using the extracted data with model coefficients from SDMs to run indicator under different scenarios in ALCES Online.
The second scenario is known as the Preferred Forest Management ("PFM") scenario and corresponded to the spatial harvest strategy used by Al-Pac in their 2015 Forest Management Plan. The PFM scenario also solves for the tradeoffs described in the EBM scenario but it also defers large tracts of habitat within Woodland Caribou range (primarily black spruce/larch bogs and fens) from harvest for the first 20 years of the simulation. After the 20-year deferral model allocation within the Caribou Zone is permitted. Both scenarios were run separately on each of Al-Pac's 12 Forest Management Units (FMUs) as required by provincial planning standards (Government of Alberta, 2006).
From each of the two scenarios, we extracted year 0, year 10, year 20, and year 50 outputs as shapefiles for each of Al-Pac's 12 FMUs in each scenario, describing the polygons harvested, when they were harvested, age at harvest, cutblock size, type of cover class, and harvest volume (Figure 2). We spatially unioned these layers with an Alberta quartersection shapefile layer clipped to each FMU's boundaries. This allowed us to calculate the cumulative area of each forest ageforest type per quarter-section. We classified forest age-classes based on the management age of each forest stand since its origin, based on the stand and age-class categories used in the bird SDMs. The abundance of each species per quartersection based on predicted habitat conditions were computed from cure4insect for the 20 bird species for each time period and scenario.
The Al-Pac FMP model underlying the Patchworks scenarios in this paper assumes that the only action that can change the age of a forest stand is a harvest event. When a stand is harvested, stand age is reset to zero and the forest type remains the same (regenerating stands are the same as the original stand). Stand succession (the endemic process of senescence and renewal) is captured in the long-lived yield curves with the assumption that long-term standing volumes decline and stabilize at 50% of peak culmination.

Natural Range of Variation
Patchworks is an optimization tool for harvest scheduling, but does not explicitly allow for dynamic changes in fire or other human land-uses. Thus, we used a spatially explicit land use simulator called LANDMINE (Andison, 1996(Andison, , 2005Andison and Forest, 2000) to compare the effects of fire to the harvest results from Patchworks. Specifically, we compared the two harvest scenarios to a scenario devoid of human activities whereby forest age-structure was set back by fire instead of harvest. LANDMINE uses a dispersal algorithm to spread fires from one pixel to another probabilistically based on fuel-type and topography. Fire movements were calibrated to create different fire shapes and unburnt island remnants based on empirical data. Fire size was controlled by an equation representing the actual fire size distribution for each landscape. Ignition location probabilities were loosely based on historical lightning probabilities. Finally, the total amount of forest burnt in any single time step was based on historical areas burnt. A burn frequency of 63.5 years was used for the Al-Pac FMA. The "NRV" scenario was a Monte Carlo simulation that replaced human footprint with natural disturbance footprints (fire) through time, in which no logging nor any other human activity occurred. Once human footprint was removed, the NRV scenario was run for approximately 1000 simulated years (Andison and Forest, 2000). The NRV simulation of fire and the resulting forest-age structure was developed by DA at Bandaloop Landscape Ecosystem Services 3 .
From the NRV simulation 100 snapshots -each separated by 10 years -of the Al-Pac FMA were extracted. Each snapshot consisted of a grid of 200-m square cells across the entirety of the Al-Pac FMA. Estimated forest age and dominant tree species or non-forest vegetation in each cell were calculated. The amount of each vegetation age-class varied by snapshot. We converted each snapshot to shapefiles and imported the shapefiles into cure4insect to calculate bird abundance in each cell (Figure 3). Once we had 100 abundance predictions of each species, we estimated the mean, standard deviation, and confidence interval estimates for those predictions.

Cumulative Effects
While Patchworks and LANDMINE are particularly good at spatially representing harvesting and fire strategies and optimizing solutions, they were not designed to track cumulative effects. To explore how boreal birds may respond to the cumulative effect of multiple land uses and fire, we created five scenarios using a spatially explicit land use cumulative effects simulator called ALCES Online 4 . To make these simulations as realistic as possible, we used the actual harvest plans from the PFM scenario in Patchworks (the most likely harvest plan) along with theoretical future trajectories for bitumen development and fire to determine how these drivers affect bird response over the next 50 years. The consequences of simulated anthropogenic and natural processes were assessed by tracking changes in landscape composition, forest age, and forest origin (i.e., burn or harvest). The simulator is cell based, with each cell's composition tracked as proportional coverage by various natural and anthropogenic cover types. ALCES was initialized by calculating the current proportional composition of each 200-m cell by forest types, other cover types, and footprints. Forest area, age, and origin was based on the Alberta Vegetation Inventory as in the Patchworks and LANDMINE tools. Coverage by different non-forest, terrestrial and aquatic vegetation types -which may also influence bird abundance -was calculated using the Earth Observation for Sustainable Development (EOSD) land cover and the AltaLIS hydrology datasets. Current location and extent of anthropogenic footprints other than harvest (roads, well sites, pipelines, seismic liens, mines, rural residences, settlements, transmission lines, industrial features) were obtained from the Alberta Biodiversity Monitoring Institute human footprint inventory as well as other sources such as Alberta Energy Regulator, CanVec, and OpenStreetMap. Proportional composition, age, and origin of each cell was then modified during simulations to track the effect of new developments, reclamation, and fire on landscape composition and age. A baseline "Al-Pac BAU" scenario incorporated forestry, bitumen, settlement, road, and gravel pit development as well as fire. The baseline scenario was modified to create the following additional scenarios: "Seismic Restoration, " which restored seismic lines to improve woodland caribou habitat; "No Energy, " which excluded the effect of future energy development; "No Fire, " which excluded the effect of future fire; and "Increased Fire, " which doubled the fire rate to incorporate the projected effect of climate change (Boulanger et al., 2014).
We used the following assumptions to simulate fire and nonforestry land use. Annual area burned (284.8 km2/year) and the fire size class distribution was based on fires occurring in the study area over the past 75 years  rather than on the NRV scenario. Fire location was influenced by forest type and age (Bernier et al., 2016). Simulation of in situ (i.e., well-based) and mineable bitumen development assumed production trajectories consistent with projections by the Alberta Energy Regulator (2016) and National Energy Board (2016) for the first 25 years, after which production plateaued based on the expectation that bitumen production will stabilize in the long-term (Millington and Murillo, 2015;Straatman and Layzell, 2015). The sequencing of mining projects during the simulation was based on anticipated project start-up dates. For the first 25 years of the simulation, in situ bitumen development occurred at operational, approved, and applied projects. Thereafter, location was influenced by bitumen pay thickness. The number of new production wells and mine area required through time to meet the production trajectory was based on productivity assumptions from previous studies (Wilson et al., 2008;Alberta Energy Regulator, 2015). Production wells and exploration footprint (exploration wells, seismic lines) were aggregated around central industrial plants, and pipelines linked the projects to the existing pipeline network. Energy sector footprints remained for the duration of the 50-year simulation because reclamation trajectories for footprints such as seismic lines (Lee and Boutin, 2006) and mines (Rooney et al., 2012) are slow and uncertain. The exception was the "Seismic Restoration" scenario for which seismic lines were reclaimed to natural land cover after 20 years. Rural residential and urban settlements were simulated to expand as per the Government of Alberta's population growth rate for the region. New rural residences (i.e., acreages) occurred within 1 km of existing rural settlement footprint, and urban expansion occurred at the periphery of Fort McMurray, the largest city in the study area. Simulated roads connected new well sites, acreages and timber harvest areas to the nearest existing road. New gravel pits were simulated based on the current ratio of road to gravel pit area; gravel pits were located within gravel deposits in proximity to new roads.
Rather than compute bird abundance by predicting total abundance using the various shapefiles created by Patchworks FIGURE 3 | Examples of "snapshots" of forest age-structure the Al-Pac FMA from simulations of natural disturbance (a 63.5-year fire cycle) in the total absence of human footprint (NRV scenario). In this program, younger forest stands resulted from fires which were also used to replace human footprint over 1000 years of simulated time. We ran this simulation for 1000 years, then used 100 "snapshots" from the simulation to estimate natural range of variation (NRV) of different forest stands and ages. Pure stands of white spruce were not included in the simulation due to their rarity in the Al-Pac FMA. In contrast to the NRV scenario, harvest was the only source of disturbance setting back forest stand age to 0 in the Patchworks scenario and given the small amount of each forest management unit harvested over 50 years in the Patchworks scenarios, forest age on average increased within individual FMUs over 50 years. and LANDMINE as inputs in cure4insect, we took the underlying model coefficients from the SDMs in cure4insect to directly track birds as indicators in the Alces Online tool with reporting of mean density and population size at each 10-year mark in a 50-year simulation.

Harvest Planning for Caribou Management
Most bird species associated with mature and old deciduous and mixed-wood forests were predicted to increase over 50 years (Figure 4, Table 2, and Appendix I). Bird species associated with mature and old deciduous and mixed-wood forests tended to exhibit a larger population increase in the presence of caribou habitat deferral (PFM scenario) than without (EBM scenario). This included Black-throated Green Warbler (9% under PFM vs. 7% under EBM: 4701 more birds after 50 years), Brown Creeper (58% vs. 56%: 20849 more birds), Canada Warbler (36% vs. 30%: 18736 more birds), and Pileated Woodpecker (15.4% vs. 14.7%: 415 more birds). An exception was Ovenbirds, with 24375 fewer birds under the PFM scenario (−1.0% under PFMS vs. −0.6% under EBM), although the differences were small relative to the initial population of this species (5283326 Ovenbirds in Year 0). The density of species associated with older mixed-wood and deciduous forests tended to be lower in FMU A-14, which experienced a large forest fire in 2016 at the start of the simulation (Figures 4-6 and Appendix I).
Bird species associated with older coniferous forests tended to respond more positively to the PFM scenario than the EBM scenario, although the difference was relatively minor.  (Figure 4, Table 2, and Appendix I).
When we compared forest age-class amounts in Year 50 from both harvest scenarios, we did not find large percent differences in the amount of each forest age-class. Four percent of the forestage classes were ≥50% more abundant by Year 50 under the PFM scenario, particularly black spruce < 20 years old, but also some mixed-wood and white spruce 20-60 years old. Just over one percent of the forest-age classes were >50% less abundant under the PFM scenario than the EBM scenario, primarily pine <20 years old originating from harvest. Nearly 95% of forest ageclasses across all forest management units showed smaller relative differences among the harvest scenarios (Appendix II).
Average forest age increased over time in both harvest scenarios. On average in Year 0, most quarter-sections in the FMA were dominated by deciduous and mixed-wood forests <60 years old and coniferous forests <80 years old. By Year 50 in both harvest scenarios, much of the FMA had become dominated by older forests (Figure 7), with relatively large gains in the percent cover of black spruce and pine >80 years old in each quarter-section. Mixed-wood and white spruce occupied small proportions of the FMA across the whole period and large tracts of older deciduous forest were strongly reduced in both scenarios over 50 years (Figure 8). However, by Year 50, declines in deciduous and mixed-wood forests >60 years old and white spruce >80 years old were offset by new older stands that developed from harvested areas existing prior to Year 0 of the harvest scenarios (Appendices I, II).

Natural Range of Variation
Most species associated with older coniferous, mixed-wood, or deciduous forests were projected to be more abundant under the harvest scenarios than in the NRV scenario (Figures 4-6 and Table 1), or after 50 years were at the higher end of the number predicted in the NRV scenario. These increases are presumably due to greater amounts of older deciduous and coniferous (Figure 8) rather than mixed-wood forests after 50 years in both harvest scenarios, since older mixed-wood forests were still uncommon after 50 years in both harvest scenarios. We generally projected more deciduous and fewer young mixed-wood stands under the two harvest scenarios than in the NRV scenario (Figure 9 and Appendix II). The amounts of deciduous and mixed-wood age-classes in the harvest scenarios were based on the initial amounts in the Alberta Vegetation Inventory, in which many natural mixed-wood stands had already been harvested and converted to either pure coniferous or deciduous stands during replanting (Hobson and Bayne, 2000). Blackpoll Warbler was less abundant under the harvest scenarios than the NRV scenario, possibly because fires in the NRV scenario created enough alternative habitat that this species could use. Most of the analyzed species that were associated with younger or open habitats in boreal forests were more abundant under the NRV scenario than either harvest scenario (Figure 4 and Table 1). We reasoned that as the overall forest age structure became older, forests became less suitable for most of these species. Young pine and black spruce stands may have been less abundant while older pine and black spruce were more abundant under the harvest scenarios due to fire suppression: under the NRV scenario, fires were more likely to burn older pine and black spruce, converting burned stands to young stands (Figure 9). In contrast, American Three-toed Woodpecker and Olive-sided Flycatcher, which are associated with burns within boreal forests, were more abundant under the harvest scenarios than the NRV scenario, despite the fact that there would be fewer forests or open areas after 50 years in the harvest scenarios. Since these two are tree-nesters that forage on or from trees (in contrast to the other species that nest in trees or shrubs but forage in clearings), fire rates associated with the NRV scenario over 1000 simulation-years may be reducing some important habitat features for these species, even if fire is associated with habitat for these species.
While we focused on analyses of 20 bird species for this paper, we have included scripts for running the same Patchworks and NRV analyses on other species at https://github.com/borealbirds/ Patchworks-NRV-cure4insect.

Cumulative Effects
When we considered cumulative effects of harvest, fire, and habitat conversion by non-forestry footprint, 13 of 20 species increased under current rates of fire, harvest, and energy sector (2) the Ecosystem-Based Management/Natural Disturbance Model Scenario (EBM) without harvest deferral (blue), relative to the predicted Natural Range of Variation (NRV) in Black-throated Green Warbler abundance in the absence of human footprint. Red line = mean predicted Black-throated Green Warbler abundance per Forest Management Unit in the NRV Scenario. Light blue dashed lines = 50% confidence intervals for mean predicted Black-throated Green Warbler abundance in the NRV Scenario. Black dashed lines = 95% confidence interval estimates in the NRV Scenario. SDMs predicted that Black-throated Green Warbler densities were highest in very old deciduous and mixed-wood stands (>80 years old) and average deciduous and mixed-wood forest age increased in most FMUs over 50 years in the Patchworks simulations, while average deciduous and mixed-wood forest age was reduced by simulated fires in the NRV scenario. development ("Al-Pac BAU Scenario"), with the largest increases being observed for American Three-toed Woodpecker (126%) and Black-throated Green Warbler (120%). The declining species were Bay-breasted Warbler, Blackpoll Warbler, Ovenbird, Palm Warbler, and Rusty Blackbird (as in the harvest scenarios), and Boreal Chickadee and White-winged Crossbill (unlike in the harvest scenarios). The largest decrease was observed for Rusty Blackbird (42%) (Figure 10 and Table 3).
In general, species that were more abundant under the harvest scenarios than the NRV scenario were relatively less abundant under a higher burn rate. This negative response was measured as a larger decrease or smaller increase in the "Increased Fire" scenario relative to the "Al-Pac BAU" scenario and/or as a larger decrease or smaller increase in the "Al-Pac BAU" scenario relative to the "No Fire" scenario ( Figure 11). Most species associated with older forests, along with American Three-toed Woodpecker and Olive-sided Flycatcher, were more abundant in the harvest scenarios than the NRV scenario and also responded negatively to fire in the cumulative effects scenarios. The negative effect of fire on Olive-sided Flycatcher was small: Olive-sided Flycatcher increased in all scenarios but increased less over 50 years in the "Increased Fire" scenario (5%) than the "Al-Pac BAU" scenario (6%) and "No Fire" scenario (9%). The largest negative responses (≥25% difference between "Al-Pac BAU" and "No Fire" population projections at year 50) were for American Three-toed Woodpecker, Blackthroated Green Warbler, Brown Creeper, Canada Warbler, and Cape May Warbler. An exception to the pattern was Ovenbird, which responded positively to fire in the cumulative effects scenarios. It is worth noting, however, that Ovenbird densities were initially reduced by a higher fire rate until near the end of 50 years in the cumulative effects scenarios (Figure 10). Other species that responded positively to a higher fire rate (Blackpoll Warbler, Black-backed Woodpecker, Northern Flicker, Palm Warbler, and Rusty Blackbird) were more abundant in the NRV scenario than the harvest scenarios (Figure 10 and Table 3).
Most (16 of 20) species responded negatively to energy sector development in the cumulative effects scenarios. This negative response was measured as a larger increase or smaller decrease in the "No Energy" scenario relative to the "Al-Pac BAU" scenario and/or as a larger decrease or smaller increase in the "Al-Pac BAU" scenario relative to the "No Energy" scenario. The largest negative responses (≥25% difference between "Al-Pac BAU" and "No Energy" population projections at year 50) were observed for Black-throated Green Warbler, Canada Warbler, and Western Tanager. Species responding positively to energy sector development were limited to Blackpoll Warbler, Black-backed Woodpecker, Rusty Blackbird, and White-winged Crossbill (Figure 10 and Table 3).
Twelve of sixteen species that responded negatively to simulated levels of energy sector development responded positively to restoration of seismic lines. The largest positive responses to seismic line reclamation (≥5% difference between "Al-Pac BAU" and "Seismic Restoration" population projections at year 50) were observed for Black-throated Green Warbler, Brown Creeper, Canada Warbler, and Western Tanager (Figure 10 and Table 3).

Harvest Planning for Caribou Management
In our study, outputs from Patchworks were used to predict bird species abundance under an ecosystem-based management versus caribou-conservation strategy. Differences in population projections for most birds were small with an absolute difference in percent population change between harvest scenarios < 6% on average for all species). While harvest locations differ considerably between the PFM and EBM scenarios in the first 20 years of the simulation, most forest stands over the entire FMA remained unharvested over 50 years in both scenarios as they have not yet become old enough to be harvested. We expected that the caribou-conservation strategy would have negative effects on birds that rely on large patches of older deciduous and mixed-wood forests. The SDMs on which bird population projections are based on emphasize habitat amount rather than habitat configuration per se. Our standlevel modifier that adjusts local bird density based on the amount of suitable habitat surrounding the survey location does indirectly account for patch configuration, because many metrics landscape fragmentation metrics are correlated with habitat amount (Wang et al., 2014). We also considered not only surrounding suitable habitat but amount of water and different types of human footprint at the landscape scale. We believe that a combination of these variables is predictive and interpretable, which were our main concerns from an application perspective. However, we also recognize that the concept of a patch for boreal birds is a fundamental challenge in these types of models as species with small territories may treat a clump of conifers in an otherwise deciduous-dominated forest as a patch, while a bird with a larger home range may view that same area as a mixed-wood. Taking a species-centric view of patch size is needed to address this issue in future simulations and will involve modeling local and landscape level stand characteristics by accounting for territory size differences among species (Westwood et al., 2019). Importantly, total area harvested was similar in both the ecosystem-based management and caribou-deferral plans. After deferral ended in the PFM scenario, the same locations were eligible for harvest within both scenarios. Thus, even when large changes occurred in amounts of some forest age-classes underwent large changes, there ended up being similar amounts of most forest age-class types in both scenarios at the end of 50 years. As a result, habitat available for and predicted abundance of bird species in the FMA was similar in both scenarios after 50 years. Some studies suggest that habitat configuration and fragmentation effects are insignificant for boreal birds in landscapes where forest harvest is the primary agent of habitat conversion, except in extremely fragmented landscapes (Andren, 1994;Schmiegelow and Mönkkönen, 2002).
Our harvest scenarios explored the influence of harvest deferral on one major harvest strategy underlying EBM the location of harvest areas. Apart from varying size, shape, distribution and location, our harvest areas were all assumed to be clear-cuts based on the harvest practices modeled by our SDMs. Other EBM-based practices in Canada like partial cuts, shelterwood cuts, structural retention, and understory protection have been studied throughout Canada for their effects on tree mortality (Thorpe and Thomas, 2007), subsequent tree growth (Montoro Girona et al., 2017, 2019, understory protection (Burke et al., 2008), and biodiversity (Fenton et al., 2013;Huggard et al., 2014;Charchuk and Bayne, 2018). These harvest strategies may be more appropriate than traditional clear-cuts for emulating natural disturbance in regions where forest fires are less frequent than insect outbreaks and other disturbances. As regional SDMs are developed to account for the effects of these other harvest strategies on birds and other wildlife, it will become possible and desirable to project long-term bird abundance under these other strategies, using programs like Patchworks.
Since forest stand age was not set back by fire, other natural disturbances, or non-forestry human footprint in the Patchworks scenarios, increasing average forest age explains why we projected smaller numbers of bird species associated with younger forests in Patchworks relative to the absence of human footprint (including fire suppression) in the NRV scenario. Some bird species associated with younger forests, like Olive-sided Flycatcher, also use habitats that were not modeled in the harvest scenarios (e.g., larch fens, shrublands) (Robertson and Hutto, 2007). For this reason, simulators that also model non-forested vegetation, unlike Patchworks, may provide more realistic projections of habitat available for species such as these.
Previous studies (Bichet et al., 2016;Drever et al., 2019) have quantitatively assessed if conserving or managing habitat for woodland caribou also protects significant habitat for other species. While we did not explicitly test for the role of caribou as an umbrella species for boreal birds in our study, our harvest scenario results suggest that harvest deferral for 20 years within caribou conservation zones does not have large effects on the populations of bird species across the Al-Pac FMA over 50 years. Incidentally, deferral of harvest to benefit caribou in the Al-Pac FMA resulted in more available habitat created or remaining after 50 years for several species associated with older boreal forests. These species included the federally listed Canada Warbler 5 (see footnote 1) and some species of conservation interest in Alberta like Bay-breasted Warbler, Black-throated Green Warbler, and Cape May Warbler (Alberta Environment and Sustainable Resource Development, 2014). Harvest deferral was associated with reductions in two other federally listed species, Olive-sided Flycatcher (678 fewer under PFM after 50 years) and Rusty TABLE 2 | Initial population size and projected percent change in population of 20 species in the Al-Pac FMA over 50 years, along with projected response of each species, under two scenarios: (1) the Preferred Forest Management Scenario (PFM) that included harvest deferral within caribou habitats; and (2) the Ecosystem-Based Management/Natural Disturbance Model Scenario (EBM) without harvest deferral, relative to the predicted Natural Range of Variation (NRV) in abundance of each species in the absence of human footprint.

Habitat
Year 0  "Increased * " or "Decreased * " = less than a 1 per-cent difference in population after 50 years between scenarios with and without harvest deferral for caribou.
Frontiers in Ecology and Evolution | www.frontiersin.org Light blue dashed lines = 50% confidence intervals for mean predicted Palm Warbler abundance in the NRV Scenario. Black dashed lines = 95% confidence interval estimates in the NRV Scenario. SDMs predicted that Palm Warbler densities were in black spruce stands 40-80 years old, which are preferred as habitat by woodland caribou but are not harvested by Al-Pac; therefore habitat availability for and projected numbers of Palm Warblers were virtually the same under both Patchworks scenarios. Average black-spruce forest age increased in most FMUs over 50 years in the Patchworks simulations, reducing habitat for Palm Warblers, while simulated fires reduced forest age and increased habitat for these species in the NRV scenario.
Blackbird (4420 fewer under PFM after 50 years), but the percent population change over 50 years was small for both species (<5%). Based on the species we examined, shifting harvest pressure away from landscapes containing preferred caribou habitat does not appear to have large negative consequences for other species at risk and may incidentally benefit some declining birds as well. Additional management actions for those species at risk that do decline under the harvest deferral scenario could be considered within individual F.M.U.s where there is less harvest deferral occurring.

Natural Range of Variation
It may be intuitively surprising to expect average forest age and available habitat for birds associated with older forests to increase under harvest scenarios relative to the absence of human footprint. In addition to harvest, human footprint includes active suppression of fires. Recent fires like the Horse River Fire created many newly initiated forest stands in the Al-Pac FMA in the years just prior to this study. In fact, some of these fires may have been more severe in areas with a long history of fire suppression, due to accumulation of flammable material (Arienti et al., 2006). As a result, the average forest age in the Al-Pac FMA was low relative to other boreal forest regions in Alberta, in Year 0 of the Patchworks and ALCES Online scenarios. Since harvest resets forest age for only a small proportion of the total available forest, forest age will on average increase in the absence of other forest disturbances.
An assumption underlying some harvest scenario results (the projected increases bird species associated with older forests) is that fire suppression by humans is completely successful. If that assumption is unmet (Arienti et al., 2006), then projected increases of many bird species with increasing forest age will be smaller or even turn to decreases. It should also be noted that the NRV scenario was based on current burn rates but burn rates in boreal forests are predicted to increase with climate change in Canada (Bhatti et al., 2003;Krawchuk et al., 2009). Finally, it should be noted that the NRV scenario modeled only one natural disturbance, fire, but other disturbances like droughts, wind-throw, beavers, and insect outbreaks could also affect tree mortality and hence forest age structure and habitat availability for birds (Seidl et al., 2017;Navarro et al., 2018;Cadieux et al., 2020). These disturbances could further reduce habitat for species associated with older forests while creating habitat for other species. Furthermore, these disturbances are also expected to increase with climate change (Seidl et al., 2017;Navarro et al., 2018;Cadieux et al., 2020).

Cumulative Effects
Bird species associated with older forests (Black-throated Green Warbler, Brown Creeper, Canada Warbler, Cape May TABLE 3 | Initial population size and projected percent change in population of 20 species in the Al-Pac FMA over 50 years under 5 land use scenarios, along with projected response of each species to simulated current fire rates, doubled fire rates, energy sector development without seismic line restoration, and energy sector development with seismic line restoration.

Habitat
Year 0  "Increased * " or "Decreased * " = less than a 1 per-cent difference in population after 50 years between scenarios with and without a particular disturbance. All scenarios share the same amount and locations of harvest disturbances. "Al-Pac BAU" assumes that some forest habitat is either set back by fire or converted to other land uses by energy sector development but also by agriculture and urbanization. "Seismic restoration" assumes the same amount and location of fire, harvest, and non-forestry footprint as "Al-Pac BAU," but also assumes that seismic lines are successfully reclaimed and start regenerating to forests. "No Fire" assumes that no future forests are burned by fire, but levels of conversion to non-forest habitats are the same as in "Al-Pac BAU". "No Energy" assumes that no new energy sector development occurs, but levels of fire and conversion to other non-forest habitats are the same as in "Al-Pac BAU." "Increased Fire" assumes that the same amount of forest habitat is converted to other land uses like energy sector development, agriculture, and urbanization as "Al-Pac BAU," but twice the amount of forest is burned each decade, reducing older forests relative to "Al-Pac BAU." Frontiers in Ecology and Evolution | www.frontiersin.org Warbler) generally increased over 50 years in both the Patchworks and ALCES Online scenarios emphasizing that the forest age still increased over time throughout the Al-Pac FMA in the cumulative effect scenarios. However, fire and energy sector development generally reduced habitat for these species, resulting in smaller projected increases relative to the "No Fire" and "No Energy" scenarios. Forest disturbance by energy sector development in northeastern Alberta is substantial: in years prior to the global economic downturn of 2008 it was even comparable to the amount of harvest by the forestry sector (Schneider and Dyer, 2006;Brownsey and Rainer, 2009). Surprisingly, fire and energy sector development also reduced habitat for Olive-sided Flycatchers in the NRV and ALCES Online scenarios, despite the species' preference for younger forest stands, burns, shrublands, and open lands as habitat (Robertson and Hutto, 2007). However, the negative response of Olive-sided Flycatchers to fire was small relative to species associated with older forests. Doubling the burn rate reduced the population growth rate of this species relative to the "Al-Pac BAU scenario) from 6% to 5% while a scenario lacking fire had 9% population growth. Since this species nests in coniferous trees, which had a higher burn probability in the scenarios, the simulated fires might have reduced some nesting habitat for  (2) the Ecosystem-Based Management/Natural Disturbance Model (EBM) scenario without harvest deferral (blue), relative to the median amount of those forest age-classes predicted from the Natural Range of Variation (NRV) Scenario in the absence of human footprint. (1) Forest types 0-9 and 10-19 years old were treated separately when predicting bird abundance but were combined to simplify display in these plots; (2) Mixed-wood and white spruce were treated separately when predicting bird abundance in the Patchworks scenarios but were combined in these plots for comparison against the NRV Scenario. The reason for doing so is that so little pure white spruce occurred in simulations of the NRV Scenario that white spruce was treated as older mixed-wood forests when predicting bird density. this species even while theoretically creating more open habitat for this species.
Another surprising result was that Ovenbird responded positively to a higher fire rate, which was contrary to the species' lower abundance in the NRV scenario which also incorporated fire. Ovenbirds are associated with mature rather than old deciduous and mixed-wood forests and it is possible that simulated rates of fire in the ALCES Online scenarios created enough new suitable habitat over 50 years for Ovenbirds to exhibit a positive fire response. In contrast, the NRV scenario simulated rates of fire over 1000 years, which may have been long enough for fires to reduce the amount of suitable habitat FIGURE 10 | Projected mean density (# males/ha) of 20 species in the Al-Pac Forest Management Area (FMA) over 50 years under 5 land use scenarios. Projected population size can be calculated by multiplying mean density by the number of hectares in the Al-Pac FMA (6,563,755). All scenarios share the same amount and locations of harvest disturbances. "Al-Pac BAU" (violet line) assumes that some forest habitat is either set back by fire or converted to other land uses by energy sector development but also by agriculture and urbanization. "Seismic Restoration" (gold line) assumes the same amount and location of fire, harvest, and non-forestry footprint as "Al-Pac BAU," but also assumes that seismic lines are successfully reclaimed and start regenerating to forests. "No Fire" (dark blue line) assumes that no future forests are burned by fire, but levels of conversion to non-forest habitats are the same as in "Al-Pac BAU". "No Energy" (light blue line) assumes that no new energy sector development occurs, but levels of fire and conversion to other non-forest habitats are the same as in "Al-Pac BAU." "Increased Fire" (red line) assumes that the same amount of forest habitat is converted to other land uses like energy sector development, agriculture, and urbanization as "Al-Pac BAU," but twice the amount of forest is burned each decade, reducing older forests relative to "Al-Pac BAU". Generally as the amount of simulated disturbance increases from least ("No Energy," "No Fire") to most ("Increased Fire), species associated with older forests are more likely to decline or less likely to increase, while species associated with younger or mature forests are more likely to increase or less likely to decrease.
available on average to Ovenbirds. This incidental result suggests the importance of considering temporal scale when simulating cumulative effects on boreal birds.
When projecting future populations of species, the assumptions underlying forest disturbance, regrowth, and age matter greatly. Due to the assumptions underlying the FIGURE 11 | Mean density (#birds/ha) and distribution of Cape May Warbler -a songbird that responded negatively to fire and energy sector development -in the Al-Pac Forest Management Area (FMA) under four cumulative effects scenarios in ALCES Online: (1) harvest from the Preferred Forest Management (PFM) scenario + current burn rate + moderate energy sector development over 50 years ("Al-Pac BAU" scenario) (top left); (2) harvest from the PFM scenario + 2*current burn rate + moderate energy sector development over 50 years ("Increased Fire" scenario) (top right); (3) harvest from the PFM scenario + no fire + moderate energy sector development over 50 years ("No Fire" scenario) (bottom left); (4) harvest from the PFM scenario + current burn rate + no energy sector development over 50 years ("No Energy" scenario) (bottom right).
harvest and cumulative effects scenarios, forests aged over time in the simulations and so habitat increased for bird species of older forests in our harvest scenarios. Although most species increased during the 50-year simulation period, the negative effect of a higher fire rate suggests that population declines could occur if anthropogenic climate change in Alberta's boreal forests results in more area burned than what we simulated (Bhatti et al., 2003;Krawchuk et al., 2009). Given that forest fire suppression is imperfect (Arienti et al., 2006), forest fires are likely to increase in the future in Alberta. Further, tree species successfully regenerate at different rates after fires (Lieffers et al., 2003;Johnstone et al., 2010) and different kinds of harvest (Montoro Girona et al., 2017, 2019, so future simulation modeling efforts will need to account for differing levels and additional impacts of forest fires and other climate factors (e.g., increased drought) on stand replacement and habitat available to boreal birds (Cadieux et al., 2020). Exception for the higher fire rate scenario, we did not explicitly model climate change in our ALCES Online scenarios, which occurred over a shorter time frame than in Cadieux et al. (2020). The negative effect of doubling the amount of fire on the species we analyzed was consistent with negative population projections of boreal songbird species in other studies (Mahon et al., 2014;Stralberg et al., 2015;Cadieux et al., 2020).
Increases in non-forestry footprint, due mainly to the energy sector, were associated with larger reductions of bird species associated with older boreal forests, consistent with a previous study . Some energy sector footprints like seismic lines are intended to be temporary and left to regenerate to forest, although regeneration has been variable along seismic lines in boreal forests (Lee and Boutin, 2006). The "Seismic Restoration" scenario created habitat for some federally or provincially listed species (66294 Baybreasted Warblers, 16245 more Black-throated Green Warblers, 81555 more Brown Creepers, 820 more Canada Warblers, and 328 more Olive-sided Flycatchers) in the Al-Pac FMA over 50 years. Seismic restoration is one strategy being explored for improving woodland caribou habitat (Bentham and Coupal, 2015;Kansas et al., 2015). Our results suggest that just as deferral of harvest to benefit caribou did not have strong negative effects on boreal birds, restoring habitat for woodland caribou incidentally may have positive effects for some declining boreal bird species as well.
In theory, we could have tried to use one simulator to model all disturbances. ALCES Online can already do this for harvest, fire, non-forestry footprint. Another program, LANDIS-II, has been used to simulate multiple forest disturbance types like harvest, fire, windthrow, and insect outbreaks in many wildlife studies (e.g., Cadieux et al., 2020). However, these simulators do not yet account for many of the factors (e.g., socioeconomic) that must be considered in harvest plans, whereas Patchworks does. Similarly, Patchworks can now account for other disturbance types besides harvest, but these other disturbance types are not usually of interest to forestry companies, and simulators like ALCES Online are more versatile in modeling multiple disturbance types. ALCES Online can be used to remove human footprint from landscapes to simulate a lack of human footprint as in the NRV scenario; however, the NRV scenario simulator that we ran based on LANDMINE gave us estimates of uncertainty in habitat amounts, which allowed us to estimate uncertainty in bird abundance in the absence of human footprint.

CONCLUSION
This study is one of the first attempts to predict abundance of bird species under two harvest plans in response to habitat management for another species, and to compare those predictions with the likely NRV in those species' numbers in the absence of harvest, as well as with the presence of other human footprint in a northeastern Alberta landscape. By regulation, forestry companies develop forestry harvest plans in the absence of input and knowledge of activities from other industrial sectors in Alberta, because government oversight of industries is sector by sector and integrated land management is still uncommon (Kennett, 2006). In the case of forestry, forest management plans are updated every 10 years, thus allowing for the accounting of needs for both the forestry sector and wildlife such as caribou and boreal birds. Harvest scheduling software like Patchworks, when linked with wildlife SDMs, can be used to project the impact of strategic harvest plans on wildlife. In our study, we found that deferring harvest for 20 years in merchantable forest stands embedded within preferred caribou harvest was unlikely to adversely affect overall populations of boreal bird species associated with those stand types. These projections can make sense in the short-term because the location and extent of other industrial footprints are uncertain. However, given that substantial amounts of forest habitat can be removed by increases in fire (Bhatti et al., 2003;Krawchuk et al., 2009) and nonforestry sectors (Schneider and Dyer, 2006;Brownsey and Rainer, 2009), accounting for these additional non-forestry disturbances and sources of potential habitat loss is needed to evaluate fully how species could respond over the long term to "all" forms of disturbance. We found that inclusion of fire and energy sector development in addition to forestry led to lower habitat projections for most species that we analyzed, and that restoring energy sector footprint (seismic lines) to benefit caribou also benefited bird species. By projecting species abundance under a range of scenarios involving multiple industrial sectors and natural disturbance, cumulative effects simulators could facilitate future integrated landscape management for wildlife in Alberta.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

AUTHOR CONTRIBUTIONS
LL developed the manuscript, figures, and appendix, created three ALCES Online scenarios modified from the baseline scenario created by MC, and generated predictions of bird species abundance for the Patchworks, NRV, and ALCES Online scenarios in this manuscript. EB and ED helped conceive the overall study design and edited the manuscript. PS developed the bird species distribution models used to make predictions in this manuscript, developed the cure4insect package in R for generating predictions from GIS outputs, and edited the manuscript. TM developed the Patchworks scenarios and simulations and edited this manuscript. DA developed the NRV scenario and simulation and edited this manuscript. DC edited this manuscript. MC developed the baseline ALCES Online scenario and Seismic Restoration scenario and simulations in this manuscript and edited this manuscript. All the authors contributed to the article and approved the submitted version.

FUNDING
This research was funded by Alberta-Pacific Forest Industries Inc., and the Government of Canada through Mitacs Accelerate grants (Numbers IT07731 and IT11327). PS and EB developed the SDMs whose predictions were used in this manuscript. TM at Spatial Planning Systems developed the Patchworks harvest scenarios used in this manuscript. DA at Bandaloop Landscape-Ecosystem Services developed the NRV scenario used in this manuscript. MC at the Integral Ecology Group developed the ALCES Online cumulative effects scenarios used in this manuscript.