Simulation modeling accounts for uncertainty while quantifying ecological effects of development alternatives

Wildlife management often involves trade-offs between protecting species and allowing human activities and development. Ideally, these decisions are guided by scientific studies that quantify the impacts of proposed actions on the environment. However, critical information to assess impacts of proposed activities may be lacking, such as certainty in where actions will take place, which may hinder a robust impact assessment. To address this issue, we present the Development Impacts Analysis (DIA), which employs Monte Carlo simulation modeling to quantify the environmental consequences of proposed development scenarios, while accounting for uncertainty in the exact location of future development. We applied the DIA to five proposed oil leasing management scenarios under a revised management plan for the National Petroleum Reserve—Alaska. For each management scenario with differing levels of proposed development (“alternatives”), oil production pads and roads were randomly simulated in proportion to estimated undiscovered oil and following alternative-specific restrictions. We assessed habitat displacement for two caribou (Rangifer tarandus) herds, eight shorebird species, and black brant (Branta bernicla) based on reported responses to development, repeating the process 100 times for each alternative. Some habitat loss was reported for each proposed alternative, but the amount of impact varied by alternative and species. One caribou herd and most bird species indicated greatest effects in the alternative with the least restrictions on development and lesser impacts under more protective alternatives. Our results emphasized the importance of considering spatial variation in development effects and species-specific differences when evaluating management proposals. The DIA quantified potential impacts on a suite of species under proposed management alternatives, while accounting for uncertainty in where development will occur and providing confidence intervals on estimated impacts. This illustrates that uncertainty need not preclude management decisions about establishment of broad land use restrictions prior to submission of project-level proposals but can instead be explicitly incorporated into decision making. While no single management approach will likely benefit all species, use of tools such as the DIA allows managers to quantify trade-offs among species and pursue approaches that balance the needs of various taxa and other management objectives.


INTRODUCTION
Public lands in the United States provide valuable ecosystem services (Ruhl 2010) such as species conservation (Dietz et al. 2020), hunting and fishing opportunities (Titus et al. 2009, Gillespie et al. 2018, carbon sequestration (Olander et al. 2012), and economic revenues (Rasker 2006). While not a panacea for specific conservation goals such as endangered species protection (Cassidy andGrue 2000, Clancy et al. 2020), public lands can nonetheless be expected to contribute to meeting global conservation targets (Balmford et al. 2005, Dinerstein et al. 2019) and addressing impacts of climate change on biodiversity (Heller and Zavaleta 2009, Hole et al. 2009, Mawdsley et al. 2009). Federal public lands cover approximately 259 million ha or 28% of the U.S. land area (Vincent and Hanson 2020). Nearly 70% of this area (about 177 million ha) is managed by two agencies, the Bureau of Land Management (BLM; 99 million ha) and the Forest Service (78 million ha; Vincent and Hanson 2020). These two agencies have multiple-use mandates that inevitably lead to management trade-offs (Nelson et al. 2009, Ruhl 2010. Managers of multiple-use public lands must often balance biodiversity conservation with human activities that can fragment natural habitat, disrupt ecosystem processes, and reduce wildlife populations (Mascia andPailler 2011, Watson et al. 2016). To reduce such risks, environmental assessments-such as under the National Environmental Policy Act in the United States (Mandelker 2010)-may be undertaken to identify and mitigate negative environmental consequences at the planning and approval stages of projects. Environmental impact analyses are often conducted early in planning when uncertainty exists about the precise location of future actions within a broad study area. For example, while the impacts of a proposal to allow for clear-cutting in various drainages of Forest Service lands may be reviewed under the National Environmental Policy Act, the precise location and size of specific cuts themselves are often not specified or analyzed in the review.
Uncertainty surrounding future development makes accurate estimation of impacts difficult, especially for highly mobile and widely dispersed wildlife species. Such taxa are common in Alaska, which contains about 30% of the BLMmanaged land in the United States (nearly 29 million ha; Vincent and Hanson 2020). Many migratory species in Alaska rely on seasonally accessed habitat within public lands subject to ongoing development planning. Alaska is home to over 640,000 caribou (Rangifer tarandus), with the vast majority (approximately 80%) occurring in four arctic herds: the Western Arctic Herd (WAH), Teshekpuk Caribou Herd (TCH), Central Arctic Herd, and Porcupine Caribou Herd (Harper and McCarthy 2014). Globally, caribou and related reindeer (R. t. tarandus) are the most abundant large terrestrial herbivore in the circumpolar arctic (Bråthen et al. 2007), but more than 50% of migratory caribou and reindeer have been lost over the past two decades, likely due to varying climate and anthropogenic landscape change (Vors and Boyce 2009, Joly et al. 2011, Mallory and Boyce 2018, Russell et al. 2018. In addition to their importance for many ecological processes (Ballard et al. 1997, Stark et al. 2015, Heggenes et al. 2018, caribou are also central to food security and the cultural well-being of many indigenous groups across the circumpolar arctic (Wolfe and Walker 1987, Bjørklund 1990, Berkes et al. 1994. Caribou may respond to human development and activity with altered distribution and movement throughout their annual cycle (Johnson et al. 2005, Plante et al. 2018, Johnson et al. 2020), but appear most sensitive to human disturbance during the calving period when females with newborn calves may displace 4 km or more from infrastructure related to oil and gas development (Cameron et al. 2005, Johnson et al. 2020. In addition to caribou, the Alaskan Arctic also supports large aggregations of migratory and resident birds (Bart et al. 2012, Smith et al. 2017, including many species with declining global populations (Smith et al. 2020). Within Alaska, oil and gas development and associated infrastructure directly and indirectly compromise shorebird habitat (Meehan 1986, 1988. The construction of gravel pads fills in tundra, wetlands, and other landscape types that are important for nesting birds (Meehan 1986), while widely dispersed dust fallout associated with gravel pads and roads alters soil chemistry, increases permafrost thaw depth, and significantly changes vegetation communities (Walker and Everett 1987, Auerbach et al. 1997, Myers-Smith et al. 2006, Ackerman and Finlay 2019. Beyond geophysical impacts on habitat quality and extent, specific activities associated with oil and gas development have the potential to disturb individual birds (Miller et al. 1994). Geese particularly are sensitive to disturbance while molting, as molting birds have shed their flight feathers and cannot fly away from predators or other threats (Derksen et al. 1982, Miller et al. 1994, Taylor 1995, Lewis et al. 2011. Aircraft overflights may be especially disturbing for black brant (Branta bernicla), which are facing population declines (Leach et al. 2017, Sedinger et al. 2019. Displacement of brant from aircraft activity (Ward et al. 1994(Ward et al. , 1999 has resulted in spatial restrictions on aircraft activity near brant molting areas (BLM 2019a).
In this study, we estimate and compare the expected impact of various oil development scenarios on two caribou herds, eight species of shorebirds, and molting black brant. We present the Development Impacts Analysis (DIA), which employs a Monte Carlo simulation approach to analyze potential impacts of future development under hypothetical management scenarios. Our approach builds upon that of Wilson et al. (2013), with updated development assumptions reflecting technological advancements and an expanded suite of species for which impacts are assessed. As a test case, we apply the DIA to evaluate impacts under the alternatives proposed for development restrictions and stipulations in the National Petroleum Reserve-Alaska (NPR-A) under a revised Integrated Activity Plan (IAP; BLM 2019a). Specifically, we assess the differences in habitat displacement among proposed alternatives and what trade-offs exist among alternatives for different species. Such information quantifies the environmental consequences of proposed development restrictions under the various alternatives, while accounting for uncertainty in the location of future development, making it useful for informing management decisions about establishment of broad land use restrictions prior to submission of project-level proposals.

Study area and management context
In 1976, in accordance with the Naval Petroleum Reserves Production Act, the Naval Petroleum Reserve No. 4 was renamed as the NPR-A and management was transferred from the United States Navy to the BLM (USGS 2001). Covering approximately 9.1 million terrestrial ha, plus over 173,000 additional ha of bays, lagoons, inlets, and other tidal waters (BLM 2019a), the NPR-A is the largest single unit of federal public lands in the United States (DOI 2017). It encompasses a variety of ecosystems ranging from the mountainous peaks of the Brooks Range to wide coastal plain wetlands to brackish lagoons of the Beaufort and Chukchi seas (Fig. 1). These varied habitats are home to a wide array of wildlife species including over 300,000 caribou in the WAH and TCH. The TCH primarily calves around Teshekpuk Lake ( Fig. 1) and then remains in the NPR-A throughout much of the year (Person et al. 2007, Wilson et al. 2012. The WAH calves in the foothills of the Brooks Range mountains (Fig. 1), foraging in and around the NPR-A during the summer and then largely migrating south for the winter (Dau 2015, Joly and Cameron 2019, Cameron et al. 2020, Fullman et al. 2021. The NPR-A is also home to a diverse community of birds that migrate from as far as Antarctica , Saalfeld et al. 2013, Brown et al. 2017) to breed and molt in lowland habitats throughout the NPR-A (Liebezeit et al. 2011, Amundson et al. 2019). Over 5.3 million aquatic birds rely on habitat in the NPR-A, more than any other site in the circumpolar Arctic (Bart et al. 2012). In particular, the wetland complex surrounding Teshekpuk Lake is recognized by the National Audubon Society as an Important Bird Area of global significance ( Fig. 1) and as the Qupałuk Flyway Network Site by the East Asian-Australasian Flyway Partnership.
There has been little permanent development within the NPR-A (Fig. 1), with the first oil and gas production facilities only starting construction in 2013. Current petroleum facilities consist of Alpine CD 5 and Greater Mooses Tooth 1 (BLM 2019a). Greater Mooses Tooth 2 is under construction, and permitting is approved for a Willow Master Development Plan that would add a new central processing facility and at least four additional satellite well pads (BLM 2019a(BLM , 2020a. Four communities lie within the NPR-A boundary: Atqasuk (2019 population size: 248; v www.esajournals.org U.S. Census Bureau 2020), Nuiqsut (425), Utqia _ gvik (4467, formerly known as Barrow), and Wainwright (587; Fig. 1). These communities are not connected by permanent roads to other Alaskan infrastructure and can only be reached by air or boat for most of the year.
Oil and gas development in the NPR-A is governed by an IAP, finalized by BLM in 2013 (BLM 2013). The IAP seeks to balance the NPR-A's dual mandate of allowing for oil and gas leasing and protecting important surface resources, such as wildlife and habitats, by making some areas available for leasing and development with a series of best management practices and prohibiting leasing in other areas of conservation importance. A draft Environmental Impact Statement for a revised IAP (IAP DEIS), released in 2019, proposed four alternatives (BLM 2019a; Appendix S1: Fig. S1): Alternative A-the no-action alternative-maintained the existing land use and management restrictions from the 2013 IAP (about 52% of the subsurface available for oil and gas development), Alternative B slightly increased the area made unavailable for leasing and altered its configuration (50% available), Alternative C expanded the area available for development, especially around Teshekpuk Lake and in the Brooks Range foothills (75% available), and Alternative D greatly expanded the area available for development (81% available), including nearly all land around Teshekpuk Lake. In addition to the four alternatives  (Person et al. 2007, Cameron et al. 2020). Existing roads are sparse within the NPR-A and include both those already completed (CD5, Greater Mooses Tooth 1 [GMT1]) and those under construction (GMT2). Proposed roads, central processing facility (CPF) and satellite pads (BT1-BT5) are depicted as displayed in the Willow Master Development Plan that was recently permitted (BLM 2020a). Special areas where leasing is not allowed (no leasing) and where both leasing and new non-subsistence infrastructure are prohibited (no leasing, no inf.) reflect the 2013 Integrated Activity Plan (IAP, BLM 2013) and correspond to Alternative A in the IAP DEIS (BLM 2019a). Oil and gas leases shown as of January 2020. T. Lake, Teshekpuk Lake; IBA, Important Bird Area; inf., non-subsistence infrastructure.
v www.esajournals.org proposed in the IAP DEIS, we analyzed an additional alternative proposed by the Western Arctic Caribou Herd Working Group, a group of diverse stakeholders that cooperates with and advises management agencies on the conservation of the WAH and its use by people (Cleveland 2019). The group recommended that the leasing closure in the northern portion of the Brooks Range foothills from Alternative B be added to the closed areas under Alternative A to enhance protection of the WAH (Cleveland 2019). We referred to this updated alternative as Alternative A+ (Appendix S1: Fig. S1).

Simulation of future development
We simulated future development in the NPR-A under the five alternatives described above (Fig. 2), assuming sufficiently long timescales to allow development across the NPR-A. Infrastructure simulation started with central processing facilities (CPFs), which are "the operational center for long-term production" of an oilfield (BLM 2019a:B-6). Central processing facility locations approximated the locations of newly discovered oil accumulations suitable for development in the NPR-A (Wilson et al. 2013). We randomly generated CPF locations proportional to estimated undiscovered oil volume but constrained them to avoid areas off-limits to development under the given alternative and placement on top of waterbodies (Fig. 2a). Relative undiscovered oil followed Houseknecht et al. (2017), which greatly increased the estimated amount of undiscovered oil in the NPR-A relative to estimates used by prior assessments and management plans (BLM 2013, Wilson et al. 2013, due to incorporation of information from recent large oil discoveries. We used the mean estimate of undiscovered oil as this represents a balance of a higher probability of finding many smaller oil accumulations and a decreasing probability of finding larger accumulations (D. Houseknecht, personal communication). Data from Table 2 in Houseknecht et al. (2017) were related to Assessment Unit shapefiles (Garrity et al. 2011) and rasterized to a 120-m spatial resolution.
We multiplied relative undiscovered oil by infrastructure restriction level and a waterbody mask to account for additional restrictions on development (Fig. 2a). We assumed that areas indicated in the IAP DEIS as being closed to development would remain closed for the duration of the plan, without granting of waivers or exceptions (infrastructure restriction = 0). Such locations were identified from descriptions, required operating procedures, and maps provided in the IAP DEIS (BLM 2019a) matched with shapefiles provided by BLM with the IAP DEIS (available from https://eplanning.blm.gov/ eplanning-ui/project/117408/590). Existing leases are not subject to altered stipulations or closures imposed under the revised IAP, but rather maintain the stipulations in place when leases were issued (BLM 2019a). This is not reflected in the BLM maps or shapefiles of closed areas; thus, areas overlapped by existing leases (Fig. 1) were set to their current restriction level under the 2013 IAP (Appendix S1: Fig. S1). Areas where development is restricted but may occur, such as within buffers of certain rivers or with agreement of Native allotment owners, were assigned a reduced probability of development (infrastructure restriction = 0.1; Wilson et al. 2013). Only Alternative B included an area with deferred leasing, and we followed Wilson et al. (2013) in assuming that these tracts would become available for lease after the deferral period ended and so did not restrict such areas from development. All other areas received an infrastructure restriction level of 1 (i.e., no restriction; Appendix S1: Fig. S1). We further constrained CPFs to not occur on top of waterbodies, including both lakes and main rivers. Much of the coastal plain of the NPR-A is covered by wetlands and marshes. Such areas have a relatively high expected probability of undiscovered oil and have previously been developed for oil extraction (BLM 2019a); thus, we did not further restrict development from other land cover types. We matched land cover data from the North Slope Science Initiative (2013) to the 120-m resolution of the other input data and assigned values of 0 to all open water pixels and 1 elsewhere.
Each iteration simulated between three and seven new CPFs in addition to the proposed Willow CPF, which we treated as existing following the IAP DEIS' Reasonably Foreseeable Development Scenario (BLM 2019a; see Appendix S1: Table S1 for summary of model input parameters and sources). This includes equal or greater numbers of CPFs to the three described in BLM's high development scenario (BLM 2019a); however, we v www.esajournals.org felt that revising the IAP to permit only three or fewer CPFs and associated anchor fields was an underestimate, especially over the time frame considered in the DEIS of 70 yr after signing of the record of decision (BLM 2019a). Three to seven new CPFs are lower than those simulated Fig. 2. Development Impacts Analysis model workflow. (a) Infrastructure was simulated 100 times for each of the five alternatives (Appendix S1: Fig. S1) to provide a range of possible development (one example iteration from Alternative D is displayed). Simulated infrastructure under each iteration was used to quantify impacts to (b) caribou of the Teshekpuk Caribou Herd (shown here) and Western Arctic Herd, (c) eight species of shorebirds (semipalmated sandpiper shown here), and (d) black brant molting habitat and disturbed geese. See Appendix S1: Table S1 for model input parameters. R code to conduct all analyses is available in Data S1 or at https://github.com/tfullman/dia. Rel., relative; dev., development.
by Wilson et al. (2013), who included 8-15 oil CPFs (they described these as accumulations) per iteration, plus 24-49 gas CPFs. Development of natural gas deposits was not considered in our model. While there are multiple efforts ongoing to develop pipelines for transporting liquid natural gas from Alaska's North Slope to export terminals in southern Alaska, none have yet been shown to be feasible and there is currently little financial incentive for a stand-alone development of natural gas resources (BLM 2019a). Central processing facility locations were constrained to be a minimum distance of 35 km apart (Appendix S1: Table S1). Once CPF locations were established, the DIA simulated four to eight satellite pads for each new CPF (Fig. 2a, Appendix S1: Table S1). As with the CPFs, satellite pads were randomly located proportional to undiscovered oil and development constraints, and to avoid location on waterbodies. Satellite pad locations were constrained to occur within 56.3 km of their affiliated CPF and at least 6.4 km from other satellites or CPFs (Appendix S1: Table S1).
Roads connecting satellite pads to their CPF were then simulated using least-cost paths (Wilson et al. 2013) with the gdistance package in R (van Etten 2017; Fig. 2a). Cost maps were based on development restrictions and waterbodies and ranged from 0 (no development) to 1 (cost increases only with distance). As with CPF and satellite pad simulation, we assigned infrastructure restriction values of 0 in areas where roads are prohibited and values of 0.1 where development is restricted but may occur (Appendix S1: Fig. S1). Unlike with CPFs and satellite pads, water was not considered an absolute barrier to roads, but instead received a cost value of 0.05 to reflect the possibility, but high cost, of building bridges over water (Wilson et al. 2013). In a departure from Wilson et al. (2013), we connected all anchor fields (a CPF and associated satellites) to each other and existing infrastructure with roads. This followed guidance in the IAP DEIS that, "[o]perators on the North Slope have found that roadless developments present operational and logistical difficulties, so future developments are expected to be connected by gravel roads in most cases" (BLM 2019a:B-7). The one exception to this was simulated development north of Teshekpuk Lake under Alternative D. Stipulations for this alternative require development north of the lake to be isolated from southern development as industrial roads are not permitted through the narrow corridors on either side of Teshekpuk Lake in Alternative D (BLM 2019a). We followed BLM (2019a) in treating Alpine, Greater Mooses Tooth, and proposed Willow developments as existing for the purposes of the DIA simulations. v www.esajournals.org

Caribou impact analysis
We analyzed TCH and WAH displacement from simulated development during the calving period (Fig. 2b). We used a 120-m resolution resource selection function raster to represent relative calving habitat suitability for the TCH (Wilson et al. 2012;Fig. 3a) and identified the top quartile of habitat values as high-quality calving habitat (sensu Johnson et al. 2005; Fig. 3b). Recent work noted some TCH calving further west (Prichard et al. 2019), but the vast majority remained in areas around Teshekpuk Lake (BLM 2019a) indicated as suitable by Wilson et al. (2012). This habitat suitability layer remains the best available representation of TCH calving. The IAP DEIS also used this layer as part of its Fig. 3. Caribou input data for the Development Impacts Analysis model. Relative habitat suitability data from resource selection functions for the (a) Teshekpuk Caribou Herd (TCH) and (c) Western Arctic Herd (WAH) were used to identify high-quality calving habitat areas for the (b) TCH and (d) WAH, which we considered the upper quartile of RSF values. For the WAH, we used additional information on (e) calving area overlap between 2010 and 2017 to produce (f) a weighted set of high-quality pixels where high-quality pixels that had been used more frequently received increased weight. impact assessment (BLM 2019a). For the WAH, we used a relative calving habitat suitability raster, based on a similar resource selection approach as with the TCH (Cameron et al. 2020). We resampled the WAH data from its original 30-m resolution to match the 120-m resolution of the TCH data (Fig. 3c), which greatly reduced model run time while yielding similar results for proportional habitat loss (Appendix S1: Fig. S2). As with the TCH, we considered the upper quartile of relative calving suitability pixels to represent high-quality calving habitat (Fig. 3d). We used 95% contour kernel density estimates of WAH calving areas from 2010 to 2017 ( Fig. 3e; Cameron et al. 2020), derived from GPS-detected calving locations , to further refine the WAH calving suitability raster. Western Arctic Herd calving overlap values, ranging from 0 to 8 yr, were intersected with the WAH high-quality calving pixels to produce a weighted set of high-quality pixels, in which pixels that had been used more frequently for WAH calving received a higher weight (Fig. 3f). Comparable data were not available for the TCH.
For each DIA iteration, we overlaid the simulated CPFs, satellite pads and roads onto the calving rasters for each herd (Fig. 2b). Central processing facility and satellite pad locations were simulated as points and were converted to areas of coverage by generating square buffers around points with areas of 40.5 ha (BLM 2019b) and 6.1 ha (BLM 2019a), respectively. Roads were simulated as lines and were converted to areas of coverage by buffering by 9.45 m (i.e., half of the 18.9 m road footprint width; Appendix S1: Table S1). Habitat values of pixels overlapping existing and simulated infrastructure were assigned a value of 0, indicating complete loss of calving habitat. Other pixels within 4 km of infrastructure had their value degraded according to the following equation for caribou displacement from infrastructure (Wilson et al. 2013, based on data from Cameron et al. 2005): w x ð Þ ¼ expðÀ2:974 þ 1:079 Â xÞ 1 þ expðÀ2:974 þ 1:079 Â xÞ where w(x) is the proportional reduction in use by calving caribou, and x is distance from infrastructure in kilometers. We then calculated the proportion of pixels of high-quality calving habitat lost after accounting for the effects of infrastructure (Wilson et al. 2013).

Shorebird impact analysis
We identified suitable habitat within the NPR-A for eight species of shorebirds: semipalmated sandpiper (Calidris pusilla), American golden-plover (Pluvialis dominica), black-bellied plover (Pluvialis squatarola), dunlin (Calidris alpina), longbilled dowitcher (Limnodromus scolopaceus), pectoral sandpiper (Calidris melanotos), red phalarope (Phalaropus fulicarius), and red-necked phalarope (Phalaropus lobatus). Based on species distribution models and species-specific threshold values from Saalfeld et al. (2013), we classified all habitat suitability data into predicted suitable and unsuitable habitat for each of the eight species (Fig. 4). Within the suitable habitat, we identified high-quality suitable habitat using the upper quartile of values (Fig. 4). This allowed us to take advantage of the species-specific thresholds identified by Saalfeld et al. (2013) while maintaining a consistent accounting of high-quality habitat across taxonomic groups.
We analyzed shorebird habitat loss by intersecting the physical footprint of simulated infrastructure, buffered by 100 m to account for documented proximity effects such as dust fallout (Walker and Everett 1987, Myers-Smith et al. 2006, BLM 2019b, with predicted suitable habitat and set habitat values to zero where overlapped by the buffered area (Fig. 2c). Unlike with the caribou analysis, there was not sufficient information available to discount habitat values with distance from infrastructure, so we only considered the buffered footprint of simulated and existing infrastructure. We calculated the proportion of high-quality suitable shorebird habitat lost after accounting for the effects of infrastructure.

Brant impact analysis
To estimate potential impact to molting black brant, we calculated the area of molting habitat and estimated number of individual geese that would be disturbed by helicopter overflights and surface development (Fig. 2d). Helicopter flights are expected to occur between CPFs and satellites (BLM 2019a). We modeled the area in which helicopters would be at a low enough altitude v www.esajournals.org and a close enough lateral distance to disturb molting birds based on previous studies, which indicate that molting brant in the Teshekpuk Lake area are consistently disturbed by helicopter overflights lower than 500 m and closer than 3570 m laterally (Jensen 1990, Miller et al. 1994. We assumed that helicopters would follow IAP DEIS requirements to maintain a minimum altitude of 457 m above the Teshekpuk Lake Caribou Habitat Area (which spatially contains the Teshekpuk Lake Goose Molting Area; BLM 2019a). We rounded this up to conservatively estimate helicopters would be above 500 m except during take-off or landing from CPFs or satellite facilities. This ignores landings at areas away from facilities, such as for clean up after the ice road season, since locations of such landings could not be estimated reliably. Safe landing Fig. 4. Shorebird input data for the Development Impacts Analysis model. We analyzed eight species of shorebirds: (a, b) SESA, semipalmated sandpiper; (c, d) AMGP, American golden-plover; (e, f) BBPL, black-bellied plover; (g, h) DUNL, dunlin; (i, j) LBDO, long-billed dowitcher; (k, l) PESA, pectoral sandpiper; (m, n) REPH, red phalarope; and (o, p) RNPH, red-necked phalarope. For each species, the left map shows relative habitat suitability as determined by Saalfeld et al. (2013), while the right map uses species-specific threshold values from Saalfeld et al. (2013) to distinguish suitable and unsuitable habitat. We considered high-quality suitable habitat to be the upper quartile of suitable habitat values. trajectories for helicopters in low-relief areas typically involve descent angles between 7°and 12°( FAA 2019). We assumed a 10°descent angle as an intermediate value. Therefore, a helicopter landing at a CPF or satellite would be lower than 500 m within 2836 m from the landing strip (500 m × tan[10°] = 2836 m). Given that helicopters disturb birds 3570 m away (Jensen 1990, Miller et al. 1994, we buffered each CPF and satellite pad by 6406 m (2836 m where helicopter is <500 m altitude + 3570 m of disturbance distance around helicopters) to represent the zone of disturbance. We added to the zone of disturbance the effect of disturbance from the physical footprint of infrastructure, using a 100 m buffer to represent proximity effects, as with shorebirds.
We intersected the combined zones of disturbance from physical infrastructure and helicopter overflights with brant habitat north of Teshekpuk Lake (Fig. 2d). We used annual aerial molting goose surveys from the U.S. Fish & Wildlife Service to define this habitat and to calculate the most recent 5-yr (2014-2018) average abundance and density of geese within 209 lakes, coastal segments, and creeks ( Fig. 5; Shults and Zeller 2019). We excluded observations of waterbodies with partial or no survey effort. Within the combined zone of disturbance, we calculated the area of lake habitat affected and estimated the affected number of brant. Some brant molt in other areas of the NPR-A (Lewis et al. 2010, Flint et al. 2014), but we constrained our analyses to the region north of Teshekpuk Lake (Fig. 5) as it supports the largest waterbird concentrations (Shults and Zeller 2019) and was the only location with sufficiently robust abundance and density data.

Statistical analysis
The DIA modeling procedure, encompassing both infrastructure simulation and species impact analyses described above, was replicated 100 times for each alternative to reflect uncertainty in where oil accumulations may be discovered and where development may occur. We confirmed that 100 iterations of the DIA were sufficient to capture the variability for each alternative by calculating the coefficient of variation for each herd/species' impact results and ensuring it reached an asymptote (Wilson et al. 2013).
If coefficients of variation failed to reach an asymptote within 100 iterations, we ran an additional 400 iterations. The amount of simulated infrastructure was recorded for each iteration and compared across alternatives to estimates given in the IAP DEIS for area of surface disturbance, road length, and area within 4 km of development. This was done using both the base DIA runs and a subset of runs simulating exactly three new CPFs per iteration, to match the IAP DEIS high development scenario (BLM 2019a).
For each species/caribou herd, we calculated the mean habitat loss (caribou and shorebirds) or goose numbers and lake hectares affected (molting brant) under each alternative. Because the distribution of impact results strongly diverged from normality for some alternatives, we evaluated difference of means among alternatives using a Kruskal-Wallis rank sum test (Hollander and Wolfe 1973), with a multiple comparison test to identify which alternatives differed (Siegel and Castellan 1988). All analyses were conducted in R (version 3.6.1, R Core Team 2019). See Data S1 and https://github.com/tfullman/dia for example R code. A summary of model inputs and assumptions, with detailed references, can be found in Appendix S1: Table S1.

RESULTS
The amount of simulated infrastructure was generally similar across alternatives, but greatly exceeded that estimated in the IAP DEIS, whether run under the full DIA assumptions (Table 1), or constraining the number of CPFs to match the IAP DEIS high development scenario (Appendix S1: Table S2). The coefficient of variation for impact metrics reached an asymptote within 100 iterations for all alternatives of the caribou and shorebird analyses. The brant analyses showed no variation for alternatives A+, A, and B, and did not reach an asymptote within 100 iterations for alternatives C and D. Running these latter alternatives for a total of 500 iterations resulted in the coefficient of variation for brant impact metrics reaching an asymptote.
Impact analyses for caribou showed differential loss of high-quality calving habitat across alternatives for both herds (Fig. 6). For the TCH, habitat loss generally followed the area available for development, with greatest loss in alternatives C and D and less loss in alternatives A+, A, and B (Fig. 6a). Losses for the WAH did not follow a consistent direction across alternatives. Unlike for the TCH, development restrictions across WAH calving grounds did not simply decline across alternatives but showed more variable patterns (Appendix S1: Fig. S1). High-quality calving habitat loss for the WAH was lowest in alternatives A+ and B, intermediate in alternatives C and D, and greatest in Alternative A (Fig. 6b).
Shorebird impact analyses revealed variability across species, but with some consistent patterns v www.esajournals.org (Fig. 7). All species featured overlap in statistical significance across many of the alternatives, though some statistical differences were apparent. All but one species indicated statistically significant differences between alternatives B and D with lower habitat loss in B compared with D. For molting brant, patterns of significance across alternatives were the same for both lake area affected and the estimated number of molting geese disturbed (Fig. 8). No impact was observed in any iteration for alternatives A+, A, or B. Impacts were rare for Alternative C, with only 15/500 iterations (3%) showing any effect. The resulting level of impact was statistically indistinguishable from alternatives A+, A, and B (Fig. 8). Impacts were more common for Alternative D, with 186 of 500 (37%) displaying effects. Alternative D significantly differed from the other alternatives for both metrics analyzed (Fig. 8).

DISCUSSION
The importance of accounting for uncertainty when conducting environmental assessments has been previously noted (Wilson et al. 2013), but continues to be under-applied. For example, although the 2013 IAP benefitted from a simulation model approach to quantifying development impacts (BLM 2012), the revised IAP DEIS did not attempt to use such a method (BLM 2019a). We have demonstrated the utility of such an approach and hope that by making the DIA model code available, we will facilitate such analyses as part of future environmental assessments. The DIA model quantifies potential impacts on a suite of species under proposed management alternatives, while accounting for uncertainty in where development will occur and providing confidence intervals on estimated impacts. While our analysis was specific to the proposed developments and conservation delineations within the NPR-A, we believe the DIA approach could be applied to any geographic region, both on public and on private lands, and, as we have shown, a wide array of taxa. There are opportunities to continue to build upon this approach, such as integrating the development simulation and impact results with spatially explicit population or movement models, to better quantify the full range of potential impacts of future development.

Comparison of IAP alternatives
In our application of the DIA to the revised NPR-A IAP, some habitat loss was reported for Fig. 6. Proportion of high-quality caribou calving habitat lost across Integrated Activity Plan alternatives for the (a) Teshekpuk Caribou Herd and (b) Western Arctic Herd in the National Petroleum Reserve-Alaska. Violin plots indicate kernel probability density of impact data. Circles indicate medians, thick vertical lines indicate interquartile range. Alternatives with the same lowercase letter had habitat loss that was not statistically significantly different. Comparisons were only done within each herd; letters are not comparable across panels. each proposed alternative, but the amount of impact varied by alternative and species. This emphasizes the importance of considering spatial variation in development effects and speciesspecific differences when evaluating management proposals. For the TCH, molting brant, and five out of the eight shorebird species considered, mean impacts were greatest under Alternative D compared with the other three alternatives proposed in the IAP DEIS. The WAH was a notable exception to this pattern, featuring the lowest impact levels under the alternatives that precluded development from the largest areas (A+ and B), intermediate impacts under alternatives that made the most area available for development (C and D), and, surprisingly, greatest impacts under Alternative A.
The counterintuitive results for the WAH may stem from offset alignment between areas of concentrated calving and development restrictions. The development-restricted areas in Alternative A did not encompass the northern portion of WAH-concentrated calving areas as fully as in the other alternatives ( Fig. 1; Appendix S1: Fig. S1), increasing the possibility of development in an area of intermediate estimated undiscovered oil. While alternatives C and D make greater area in the southern portion of the WAH calving range available for development, the lower amount of expected undiscovered oil Fig. 7. Proportion of high-quality suitable habitat lost across Integrated Activity Plan alternatives for eight shorebird species in the National Petroleum Reserve-Alaska. Violin plots indicate kernel probability density of impact data. Circles indicate medians, thick vertical lines indicate interquartile range. Alternatives with the same letter above them had habitat loss that was not statistically significantly different. Comparisons were only done within species; letters are not comparable across species. SESA, semipalmated sandpiper; AMGP, American golden-plover; BBPL, black-bellied plover; DUNL, dunlin; LBDO, long-billed dowitcher; PESA, pectoral sandpiper; REPH, red phalarope; RNPH, red-necked phalarope. (Fig. 2a) reduces the likelihood of development impacts in this area. Natural gas potential, however, is much higher in this southern region (Houseknecht et al. 2017), as are hard rock mineral resources (BLM 2019a). If these or other adjacent resources were to be developed, the consequences for the WAH under alternatives C or D could be more pronounced.
In July 2020, after our analyses were completed, BLM released a final environmental impact statement for the IAP (BLM 2020b). This version of the IAP made slight alterations to alternatives B-D and added a new Alternative E, which BLM selected as their preferred alternative. This preferred alternative was finalized in a Record of Decision issued 31 December 2020 (BLM 2020c). The new Alternative E is similar to Alternative D but alters infrastructure restrictions near some rivers and in the southwestern portion of the NPR-A and makes additional area available for development (BLM 2020b). Impact results under this new alternative for the species and seasons we analyzed likely would be similar to, though greater than, those for Alternative D. A quantitative tool like the DIA should be applied to analyze impacts under the updated alternatives.

Ecological consequences and underestimation of impacts
Across all alternatives, our results should be considered minimum estimates of the expected impact of future development. In multiple instances, we intentionally chose conservative assumptions to not overestimate predicted impacts of development. In other cases, reliable information simply was not available to incorporate certain impacts, like disturbance by fixedwing aircraft or traffic volume on gravel roads. Furthermore, when identifying species' responses to development, we typically used the lower end of potential responses. For a detailed discussion of aspects that were conservative representations of future development or impact, see Appendix S2.
One key example of our DIA results reflecting minimum impacts is with caribou. We only quantified impacts to caribou during the calving period in this analysis. However, studies have shown that displacement of caribou mothers and calves away from infrastructure continues through at least the post-calving and mosquito harassment periods, despite the presence of mitigation measures and without clear evidence of habituation (Johnson et al. 2020, Prichard et al. 2020. These periods are crucial for nursing caribou to obtain high-quality forage, allowing them to meet the significant nutrient demands of lactation and replenishment of depleted body reserves used to support migration, winter survival, and subsequent year calving (Cameron et al. 1993, Gerhart et al. 1996, Barboza and Parker 2008, Taillon et al. 2013, Veiberg et al. 2017. Moreover, caribou migration can be stopped, disrupted, or delayed by infrastructure and human activity (Panzacchi et al. 2013, Frenette et al. 2020. Development or disturbances that may hinder or preclude caribou from reaching the calving grounds or other essential habitat could have further impacts beyond those reported here (see also Appendix S2). Finally, development in the NPR-A may have impacts on overwintering caribou, especially the TCH as the majority of animals remain on the coastal plain all winter (Person et al. 2007).
The estimated impacts presented here may suggest consequences at an individual or even population level. For example, molting brant in the NPR-A are limited by specific habitat requirements and are unable to alter foraging behavior to make up for lost opportunities, making them particularly sensitive to disturbance (Lewis et al. 2011). Due to the birds' energetic demands while re-growing feathers (Fox et al. 2014), even a brief startle response to a helicopter overflight can be energetically significant in an individual's life cycle (Jensen 1990, Miller et al. 1994. Similarly, caribou that fail to replenish their body stores, that exert additional energy in winter, or that are delayed or prevented during migration from arriving at their preferred calving areas could see additive impacts to those we modeled, including reduced body condition, survival, and reproductive output (Cameron et al. 1993, Taillon et al. 2013, Veiberg et al. 2017, with consequences for population dynamics.
Cumulative effects and broader-scale phenomena would likely further magnify our estimated impacts. Shorebirds face other stressors such as climate change (Galbraith et al. 2014, Weiser et al. 2018, Saalfeld et al. 2019) and loss of stopover sites throughout migration routes (Melville et al. 2016, Szabo et al. 2016). Development-related habitat loss on arctic breeding grounds, especially under Alternative D, could compound other effects and has the potential to exacerbate population declines in sensitive species , Galbraith et al. 2014. Climate-related stressors such as rain-on-snow events and landscape-scale impacts have been well documented in caribou and reindeer (Vors and Boyce 2009, Joly et al. 2011, Forbes et al. 2016, Mallory and Boyce 2018, making undisturbed access to seasonal habitats key to maintaining individual body condition and overall herd productivity (Griffith et al. 2002, Cameron et al. 2005.
Although mitigation efforts may reduce impacts for some species or in some locations, complete avoidance of negative effects is extremely unlikely and unintended consequences could increase impacts on other areas or species, as has been widely observed elsewhere (Chauvenet et al. 2011, Skogen 2015, Dodd and Sharpley 2016, Schwabe et al. 2020. For example, establishment of extensive no surface occupancy areas could result in a concentration of development along their borders to allow maximum access from directional drilling platforms. This may have the unintended effect of increasing barriers to species movement across these borders. Furthermore, because displacement from infrastructure for species such as caribou extends for several kilometers beyond the physical footprint of development (Cameron et al. 2005, Boulanger et al. 2012, Plante et al. 2018, Johnson et al. 2020, habitat would be lost functionally even within areas set aside to prohibit surface occupancy. Such unintended consequences should be taken into account when creating management prescriptions and additional buffers or other mitigation options should be considered.

Model extensions
As we describe above, the DIA considered impacts to caribou at just one stage in their annual cycle-calving-yet development may impact caribou across seasons. Incorporation of such seasonal effects will require additional submodels be incorporated within the DIA framework, such as movement models parameterized using empirical data (Patterson et al. 2008, Avgar et al. 2013, Hooten et al. 2016 to represent migratory movements and responses to infrastructure simulated by the DIA, or energetic balance models that represent fluctuations in energy gain and body reserves due to altered foraging (Russell et al. 2004). Some steps have been taken toward developing such models (Russell and Gunn 2019), but more work is needed. Nonetheless, the modeling approach we have presented is flexible enough that it could be updated with improved data or to test a wider range of assumptions, allowing the effects on predicted impact to be compared.
We focused our analysis on species of conservation concern for which robust data were readily available. However, there are other important species that should be considered in future assessments. For example, polar bears (Ursus maritimus) are iconic arctic species that are listed as threatened under the U.S. Endangered Species Act. They are affected by climate change, especially through loss of sea ice, which provides important foraging opportunities (Wiig et al. 2008, Bromaghin et al. 2015. Loss of sea ice has led to increased use of land areas by polar bears both in summer and in maternal denning (Atwood et al. 2016, Olson et al. 2017, potentially increasing risk of negative encounters with people in and around the NPR-A. Another example is the yellow-billed loon (Gavia adamsii), the largest species in the loon family and a Red List species numbering only in the thousands of individuals in Alaska (Schmidt et al. 2014). This species of conservation concern has been reviewed for inclusion on the U.S. Endangered Species list, and additional loss of habitat could increase its vulnerability (Schmidt et al. 2014). Other species of interest that could likewise be included in future iterations of the DIA if sufficient data were available include wolverine (Gulo gulo), moose (Alces alces), and spectacled eider (Somateria fisherii).

Management implications
Wildlife management on public lands with a multiple-use mandate, such as the NPR-A, often involves balancing the needs of humans and multiple wildlife species. Many studies have sought to quantify the effects of management actions or infrastructure on species after the action has occurred (Boulanger et al. 2012, Plante et al. 2018, Johnson et al. 2020. We believe it is more efficient and beneficial to devise strategies that incorporate conservation concerns in the planning process than to attempt to remedy situations that have caused conservation issues after infrastructure has already been developed. Quantification of expected impacts of different management alternatives prior to their implementation, while accounting for uncertainty, allows for this (Wilson and Durner 2020), but remains rare (Copeland et al. 2009, Wilson et al. 2013. We built upon prior quantitative modeling work (Copeland et al. 2009, 2013, Wilson et al. 2013, Wilson and Durner 2020, expanding it to consider development effects on ten species across multiple taxonomic groups, as well as two subpopulations reflected in the WAH and TCH. Our approach built directly upon the Monte Carlo simulation framework of Wilson et al. (2013) and accounted for updates in technology and business decision constraints (e.g., altering distance to development, connecting development through roads, and the array of impacts considered (e.g., through inclusion of proximity effects to birds from roads and helicopters). The opensource DIA framework we provide allows for future modifications and enhancements.
Our analyses revealed trade-offs between the amount of allowable development and wildlife conservation under different scenarios, as well as conservation trade-offs among various species within a scenario. For example, Alternative A+ was designed to mitigate impacts to the WAH by making more of the concentrated calving grounds unavailable for development. It did result in the lowest habitat loss for the WAH (Fig. 6), but for many shorebird species (five of eight of those investigated here), habitat loss was comparable to the least-restrictive alternative (Fig. 7). While we constrained ourselves to the alternatives proposed by BLM in the IAP DEIS or by stakeholders (i.e., Alternative A+), management alternatives could be generated using conservation planning tools such as MARXAN (Ball et al. 2009) or ZONATION (Moilanen et al. 2009) that seek to maximize conservation targets while minimizing costs, thus meeting multiple management objectives. This would leverage the strengths of the DIA in quantifying impacts among compared alternatives in combination with opportunities to generate candidate alternatives that minimize species impact while maximizing energy production. It is possible that no single management approach will equally benefit all species. Through the use of tools such as the DIA, however, researchers can explicitly test for trade-offs among species and pursue approaches that balance the needs of various taxa and other management objectives.

ACKNOWLEDGMENTS
Ryan Wilson provided R code allowing replication of the Wilson et al. (2013) approach. Support for TJF was provided by grants from the Wilburforce Foundation. Lois Epstein, Stephen Gerlek, and David Krause offered critical insights into development constraints and assumptions that fed into model parameters. Comments from Jeff Rasic, David Krause, and an anonymous reviewer improved this manuscript. Use of any trade names in this manuscript does not imply endorsement by the U.S. Government. At the time of acceptance, The Wilderness Society was involved in a lawsuit with BLM over the validity of the revised IAP. This had no bearing on the analyses or results reported in this manuscript, which were completed prior to selection of the preferred alternative by BLM and filing of the lawsuit. The authors were not directly involved in the lawsuit.