Spatio-temporal population change of Arctic-breeding waterbirds on the Arctic Coastal Plain of Alaska

Rapid physical changes that are occurring in the Arctic are primary drivers of landscape change and thus may drive population dynamics of Arctic-breeding birds. Despite the importance of this region to breeding and molting waterbirds, lack of a comprehensive analysis of historic data has hindered quantifying avian population change. We estimated distribution, abundance, and spatially explicit population trend of 20 breeding waterbird species using 25 years (1992–2016) of aerial survey data collected on the Arctic Coastal Plain (ACP), Alaska. The ACP is an extensive wetland complex on Alaska’s North Slope that supports millions of breeding waterbirds and includes portions of the National Petroleum Reserve—Alaska and the Arctic National Wildlife Refuge. We summarized annual counts into approximately 6-km by 6-km grid cells and analyzed data with generalized linear mixed models that accounted for survey timing and spatio-temporal autocorrelation. Geese and swans were most abundant along the coast between Admiralty Bay and Prudhoe Bay. Sea ducks, generalist predators (i.e., jaeger, gulls, terns), and loons were most abundant between Utqiaġvik and Point Lay, Alaska. Important areas for most species included the coastal fringe near Teshekpuk Lake, the Colville River Delta, and Admiralty Bay. The National Petroleum Reserve—Alaska was an important area for all species examined. Conversely, density on the coastal plain of the Arctic National Wildlife Refuge was greater than average for 20% of species. Annual population growth rates over the 25-year survey period were variable: 13 increased (range: 1.4%–13.8%), one decreased (-3.4%), and six were stable. However, even species with no overall population trend had areas of changing population size, suggesting localized conditions affected waterbird distributions on the ACP. Our results can be used to better inform land use decisions, improve monitoring of waterbird populations, and increase understanding of avian response to ecological change in the Arctic. Changements spatio-temporels chez des populations d'oiseaux aquatiques nichant dans la Plaine côtière arctique de l'Alaska RÉSUMÉ. Les changements physiques rapides qui se produisent dans l'Arctique figurent parmi les premiers facteurs responsables des changements qui surviennent dans le paysage, et pourraient donc déterminer la dynamique des populations d'oiseaux nichant dans l'Arctique. Malgré l'importance de cette région pour la nidification et la mue des oiseaux aquatiques, l'absence d'une analyse complète des données historiques a freiné la quantification des changements advenus dans les populations d'oiseaux. Nous avons estimé la répartition, l'abondance et la tendance des populations de façon explicite spatialement chez 20 espèces d'oiseaux aquatiques au moyen de 25 années (1992-2016) de données issues de relevés aériens menés sur la Plaine côtière arctique (PCA), en Alaska. La PCA est un vaste complexe de milieux humides dans la région de North Slope en Alaska et supporte des millions d'oiseaux aquatiques nicheurs; elle comprend des portions de la National Petroleum Reserve Alaska et de l'Arctic National Wildlife Refuge. Nous avons compilé les dénombrements annuels à l'échelle de cellules de grille de 6 km par 6 km environ et analysé les données au moyen de modèles linéaires généralisés à effets mixtes qui tenaient compte de la durée d'inventaire et de l'autocorrélation spatio-temporelle. Les oies et les cygnes étaient le plus abondants le long de la côte entre la baie Admiralty et la baie Prudhoe. Les canards de mer, les prédateurs généralistes (c.-à-d. les labbes, goélands et sternes) et les plongeons étaient le plus abondants entre Utqiaġvik et Point Lay. Les secteurs importants pour la plupart des espèces comprenaient la bordure côtière près du lac Teshekpuk, le delta de la rivière Colville et la baie Admiralty. La National Petroleum Reserve Alaska était un endroit important pour toutes les espèces examinées. À l'opposé, la densité calculée sur la plaine côtière de l'Arctic National Wildlife Refuge était supérieure à la moyenne pour 20 % des espèces. Les taux annuels de croissance des populations durant les 25 années d'inventaire étaient variables : 13 ont augmenté (étendue : 1,4 %-13,8 %), une a diminué (-3,4 %) et six étaient stables. Cependant, même chez les espèces globalement stables, la taille de population avait changé dans certains secteurs, indiquant que les conditions locales influaient sur la répartition des oiseaux aquatiques sur la PCA. Nos résultats peuvent servir à éclairer la prise de décision quant à l'utilisation des terres, à améliorer le suivi des populations d'oiseaux aquatiques et à accroître la compréhension de la réaction des oiseaux face aux changements écologiques dans l'Arctique.


INTRODUCTION
Several studies indicate that the distribution and abundance of some avian populations are changing rapidly in the Arctic (Hines et al. 2003, Alisauskas et al. 2011, Van Hemert et al. 2015. Over a similar period, northern Alaska has experienced warming temperatures, changes in the timing, duration, and warmth of growing seasons, loss of sea ice, permafrost degradation, changes in biogeochemical fluxes, and a net loss in terrestrial surface water (Hinzman et al. 2005, Kittel et al. 2011, Serreze and Stroeve 2015, Van Hemert et al. 2015, Gustine et al. 2017). These changes have the potential to alter terrestrial, marine, and lacustrine food abundance and availability for birds, which in turn may affect avian productivity, survival, and associated population growth. Further, climate-driven environmental impacts to birds will likely vary by trophic level (Van Hemert et al. 2015). For example, terrestrial herbivores may benefit from increased vegetation quality (Tape et al. 2013) or biomass, whereas birds that eat invertebrates (hereafter invertivores), e.g., sea ducks and shorebirds, may experience either increased or reduced food availability because of a phenological mismatch or change in invertebrate populations. Hines et al. (2003) found avian herbivore populations in the central Canadian Arctic increased substantially from the mid-1970s to the mid-1990s, but populations of King Eider (Somateria spectabilis) and Long-tailed Duck (Clangula hyemalis) declined during that same period. These inconsistencies suggest climate change may be acting on marine-and terrestrial-dependent bird populations differently, but few studies have examined broad scale avian population trends in the Arctic (though see Bart et al. 2013).
In addition to changing climate, several of the world's largest oil and gas deposits lie north of the Arctic Circle (Brownfield et al. 2012) and active oil and gas exploration and production are underway (Harsem et al. 2015, Hayes 2015. Natural resource extraction may directly affect avian populations by reducing productivity and avian density through increased predator abundance around human developments and disturbance or avoidance of areas surrounding development (Liebezeit and Zack 2008, Northrup and Wittemyer 2013, but see Meixell and Flint 2017. Additionally, oil and gas development activities and site use (especially roads) may have indirect environmental effects by exacerbating vegetation change, e.g., compaction, subsidence, and dust, altering hydrological cycles, e.g., snowmelt and ice road construction, and increasing the frequency and extent of thermokarst events (Natural Resource Council 2003, Myers-Smith et al. 2006. Inferring causative factors affecting Arcticbreeding bird populations is complicated because nearly all are migratory and therefore subject to stressors in a variety of distant migration, staging, and wintering areas. In addition, although some species show high fidelity to nesting and molting areas, some fraction of young birds are capable of dispersing in response to changing environmental conditions (Flint et al. 2014). As such, population change may be the result of birth/death processes either in the Arctic or elsewhere, and/or local immigration/ emigration driven by landscape change.
Like elsewhere in the Arctic, the Arctic Coastal Plain (ACP), an extensive wetland complex along the northern coast of Alaska, is experiencing both increased interest in oil resources and rapid environmental change (Kittel et al. 2011, BLM 2018. The ACP includes North America's largest oil field (Prudhoe Bay and associated fields), much of the National Petroleum Reserve-Alaska (NPR-A), and the 1002 area of the Arctic National Wildlife Refuge (ANWR). The ACP is also an important breeding area for millions of migratory wetland-dependent birds (Brown et al. 2007, Bart et al. 2013. Understanding abundance and population trends across the ACP is necessary to inform oil and gas leasing decisions, identify important areas for species, and identify hotspots of population change to focus subsequent research into causal mechanisms. Two studies previously summarized waterbird aerial survey data on the ACP and reported overall abundance and trends of species ) and east-west patterns in abundance (Bart et al. 2013), but spatially explicit estimates of abundance, trends, and relative importance remain largely unknown. Here, we estimated spatially explicit relative density and population trends for 20 species of breeding waterfowl and waterbirds across the ACP using 1992-2016 aerial observation index data collected by the U.S. Fish and Wildlife Service (USFWS) on the ACP Aerial Breeding Waterfowl Survey (Stehn et al. 2013). Our study expands previous efforts by accounting for survey timing and spatial covariance in estimates, and mapping density, population trends, and important areas across the ACP. Our objectives include (1) provide estimates of abundance and trends for each species across the ACP, the NPR-A, and the ANWR; (2) provide spatially explicit density and trend maps for each species; and (3) provide maps of the relative importance of an area to each species and guild. Our results can inform site-specific development activities and leasing decisions made by public and private land managers and provide insight into how birds may be responding to changes in the environment.

Study area and species
The ACP covers approximately 90,000 km² bordering the Chukchi and Southern Beaufort Seas from Point Lay, Alaska to the Canadian border just east of the Kongakut River Delta (Fig. 1,upper). It is underlain by continuous permafrost and is vegetated by sedge, graminoid, dwarf shrub, and/or moss tundra communities. The ACP landscape is treeless, topographically flat (< 100 m in elevation) and dotted by millions of shallow thermokarst lakes. This wetland mosaic hosts greater populations of breeding waterbirds (> 5 million) than any other Arctic region in the world (Bart et al. 2013). It is also home to over 10,000 human residents in small, scattered villages (Vadapalli et al. 2015). It has several known oil rich deposits including the greater Prudhoe Bay and Colville River Delta area and the potential for more oil and gas exploration in the NPR-A and coastal plain of the ANWR, e.g., 1002 area.

Aerial survey
The ACP Aerial Breeding Waterfowl Survey began in 1986 and is flown annually by the USFWS in a fixed-wing aircraft with a pilot-observer and passenger-observer. The survey consists of systematic east-west transects flown 38 m above the ground at about 175 km/hr. Procedures generally follow the standard protocol developed for the North American Breeding Waterfowl Survey (USFWS and Canadian Wildlife Service 1987). Observers count birds within 200 m on each side of the plane for a combined 400-m wide strip-transect. The aircraft's GPS location was downloaded each time an observer records an observation. The aircraft flightpath is also recorded by GPS location every 5 seconds. Observation location accuracy, stratum boundaries, transect placement, and survey timing have changed considerably through time. We analyzed data 1992-2016 because prior to 1992 observation locations were only summarized into 25.7-km (16mile) segments, therefore, we could not accurately assign these observations to grid cells.
From 1992 to 2006, the USFWS conducted two surveys annually: (1) an early-season survey in the northern portion of the ACP initiated around 10 June that covered approximately 1558 km² Avian Conservation and Ecology 14(1): 18 http://www.ace-eco.org/vol14/iss1/art18/ and averaged 8.6 days duration, and (2) an extensive survey across the ACP started around 24 June that covered approximately 1183 km², and averaged 7.3 days duration. In 2007, the late-flown survey was eliminated and the USFWS devised a new survey design conducted for the entire ACP in early June (average start date: 9 June; area surveyed: 1815 km²) with four strata flown at varying sampling intensities based on waterfowl density and management concern ( Fig. 1 in . Transects were flown in a 4-year rotating panel design and spacing between surveyed transects in a given year for the low, medium, and high density strata and the high interest stratum around Teshekpuk Lake was approximately 28 km, 18 km,9 km,and 5 km,respectively (Fig. 1,lower). Stratum number and boundaries changed through time and were based on previous estimates of bird density and physiographic characteristics, e.g., wetland density (Stehn et al. 2013). The high interest area represents the core of the BLM designated Teshekpuk Lake Special Area, which is an extremely important area for breeding waterfowl and area with high potential for containing hydrocarbon deposits (Derksen et al. 1981, BLM 2018. The southern boundary of the ACP survey area is defined by both elevation and the associated abrupt change in wetland density. These surveys generate count data, comprising observations of a species along each transect during a specific survey (early or late) in each year. Although pilots attempted to replicate the same flight line for a given transect, actual transect locations varied between the early and late survey and among years because of weather, e.g., patchy fog requiring a deviation from the established transect, and imprecision in navigating lines. Furthermore, most observations are of flying birds that are likely responding to aircraft noise flushing them from their habitat. Because of the uncertainty surrounding on-the-ground locations of birds and inconsistency in the placement and survey frequency of transects, we created 6 km by 6 km grid cells across the survey area and assigned specific observations to the appropriate cell so we could include spatially explicit covariates and temporal trends. This grid size was defined to remain consistent with previous summaries of the data (Bart et al. 2013) and maintain sufficient replication of flight-line sections within most cells. Grid cells overlapping the ACP boundary or Teshekpuk Lake were clipped to onshore area resulting in several cells < 36 km². We determined the area surveyed per cell-year-survey from the area of a 400-m diameter buffer superimposed on the transect flightpath that fell within each grid cell during each survey.
Observations represent an index to true abundance rather than an estimate of the actual population size. Waterbird distribution and detectability may change as the season progresses because of behaviors associated with nesting and incubation. Thus, the timing of a survey in a given year likely influences indices of abundance. Moreover, some waterbirds adjust their timing of breeding relative to local seasonal phenology at specific breeding locations (McKinnon et al. 2012, Liebezeit et al. 2014, making calendar date a poor measure of timing of breeding. We calculated a covariate to assess the influence of annual phenology by subtracting the date of predicted onset of the growing season in each cell-year from the survey date in each cell-year-survey to estimate the number of days a cell was surveyed before or after the onset of spring green-up for each cell-year-survey. We obtained spatially explicit (approximately 5 km x 2.5 km resolution) start-of-season-time (SOST) data for our time series from vegetation phenology indices derived from MODIS satellite imagery , Didan et al. 2015.
For cell-year-surveys with > 1 unique SOST or survey day value, we calculated weighted mean SOST and survey date for each cell. We weighted averaged SOST by the area covered by each value in each year, and weighted survey date by the proportion of the area surveyed each day in each cell-year-survey. Beyond seasonal variation in detection probability, we assume that detection rate of a species is consistent for the combined counting ability of the pilot and observer across years and space.

Model assumptions and data considerations
We made several assumptions when including data. First, counts per area surveyed within a cell are representative of the density of birds in that cell. We only included counts during a survey if ≥ 0.4 km² of a cell was surveyed during a single year-survey to ensure counts were representative of the entire cell. This value represents flying a single transect approximately 2 km within a cell. We chose this cutoff based on preliminary analyses that suggested counts were correlated within cells through time when at least this amount of area was covered. Second, consistent with USFWS survey protocol , we doubled single observations for ducks to account for unseen paired females on nests (hereafter, indicated pairs). We also doubled single observations for dark geese (Greater White-fronted Goose, Cackling Goose, and Black Brant). Other species were presumed visible enough to see nesting females, e.g., loons, Lesser Snow Geese, or are known to have male-skewed sex ratios, e.g., Scaup. We excluded observations > 2 birds, including only pairs or indicated pairs, in analyses to minimize inclusion of presettled flocks of birds or flocks of failed breeders that have moved from original breeding areas, e.g., Black Brant (Taylor 1995, Lewis et al. 2010. Therefore, although counts are restricted to single birds and pairs, we define counts as the number of breeding birds and as an index to their total breeding populations across the ACP.

Model structure
We accounted for the hierarchical nature of the counts through space and time by modeling abundance using negative binomial or zero-inflated negative binomial (ZINB) mixed-effects generalized linear models combined with a nonlinear spatiotemporal smoothing parameter, i.e., general additive mixed model, implemented in a Bayesian framework. Every cell was not surveyed in every year; thus, we used a spatio-temporal smoother in lieu of traditional random effects that allowed for autocorrelated variation in counts through space and time. Each species was not observed in all cells either because the habitat was unsuitable or because of insufficient sampling effort or imperfect detection of birds. Consequently, we included only cells where a species had ever been observed in analyses, i.e., we removed structural zeros (Martin et al. 2005). Despite removing structural zeros, data contained a high proportion of zero counts (range: 0.47-0.95). Therefore, for each species, we first attempted to fit the negative binomial distribution, and if models did not accurately reflect the true proportion of zeros in the data as determined by comparing Monte Carlo realizations of count data based on model parameter estimates to the observed proportion of zero counts in the data (Appendix 1), we refit models using the ZINB distribution.

Abundance and trend
For each species, we modeled the likelihood as the sum of birds observed in spatial grid cells c during a specific year y and survey s either as a generalized linear model using (1) the negative binomial realization of λ where Count c,t,s ~ NegBin(p c,t,s , r c,t ), r is an overdispersion term per year-cell and p c,t,s = r c,t /(r c,t + λ c,t,s ) or (2) as a two-part model where the binomial realization z that a count was zero in a given cell-year-survey was a Bernoulli draw with probability Ψ. Nonzero counts were then a negative binomial realization of the expected count λ. For the ZINB distribution, λ c,t,s = Density c,t,s × area surveyed c,t,s × (1 -z c,t,s ), whereas the negative binomial distribution did not include the z parameter. Density (km²) per cell-year-survey is then a log-linear function of random spatial and temporal variation (CELL and Y 0 , respectively), a spatio-temporal smooth term (ST; described below), late survey replicate offset (SURV), and log-linear and log-quadratic forms of survey date adjusted for the onset of spring in year t (A.DOY). Similarly, Ψ was a logit-linear function of random spatial (Z.CELL) and temporal variation (Z.Y 0 ), late survey replicate offset (Z.SURV), and logit-linear and logitquadratic forms of survey date adjusted for the onset of spring in year t (Z.DOY). Across species, the proportion of a cell surveyed varied during each year-survey. Therefore, we weighted the variance of the dispersion parameter r where log(r c,t,s ) = log (r′ -S.PROP c,t,s ); S.PROP is the normalized relative proportion of a cell surveyed and r′ is the overdispersion value at the mean cell proportion surveyed. We then derived annual indices of abundance based on model parameters for the average A.DOY and the early-flown survey. Total abundance in year t was the sum of cell-specific expected abundances in year t. We calculated overall and cell-specific population change as the geometric mean rate of change of predicted abundance (i.e., averaging over random temporal variation, Y 0 ; Sauer et al. 2014), but also report Ricker intrinsic rates of population growth g t = log(N t+1 /N t ); Dennis and Taper (1994) in the supplemental material.

Spatio-temporal smooth
Waterbird observations could be correlated in space and time because of environmental conditions (e.g., habitat), life-history traits (e.g., coloniality), and population trends that may vary spatially. Permafrost characteristics, soil type, and topography vary east to west on the ACP leading to shifts in habitat at similar latitudes (Fresco et al. 2015). Further, elevation, climate, and habitat vary with proximity to the coast (Fresco et al. 2015). Last, changes in climate varied spatially over the period considered (Appendix 2). To capture dependence in observations in space and time, we modeled spatio-temporal autocorrelation in survey data (ST) by fitting low rank, i.e., reduced dimension, thin plate spline smoothers for space and time parameters (Wood 2003, Conn et al. 2015, Hefley et al. 2017. We transformed longitude (Long) and distance to the coast (Dcoast; highly correlated with latitude, r = -0.7, and more biologically meaningful) using penalized thin plate splines into a 2-dimensional plane (SP) with k inflection points (knots) using a full tensor product smooth (te). Inflection points were spaced approximately every 30 km distance from the coast and 60 km longitude across the range of each species. This corresponds to gradients of habitat composition and climate and is consistent with previous distribution maps that showed some species changed more rapidly as a function of distance to the coast than longitude (Larned et al. 2011). We set the minimum number of inflection points in the tensor product to 3 in either dimension to ensure a sufficient number to compute a smooth for species with narrow distributions (k range: 6-32).
We applied a 1-dimensional nonlinear smooth (s) to year (YR; k = 5) (temporal model), and then created a tensor interaction (ti) term between SP and YR using the same number of k as in lower order smooths. We obtained smoothing coefficients and code for spatio-temporal effects using the jagam function in the mgcv package (Wood et al. 2016). We specified vague normal prior distributions with mean 0 and precision (i.e., 1/σ²) τ = 0.001 for fixed effects while random effects had a common mean (μ) and precision τ = 0.001. The hyperprior on μ was then normally distributed with mean 0 and τ = 0.001. We assigned uniform priors to variances ranging from 0 to 1000 (Kéry 2010). We analyzed each species independently using 20,000 iterations conducted across 5 chains, discarding the first half and retaining every 3rd iteration to reduce output file size. We assessed model convergence by examining effective sample size, ensuring the Gelman-Rubin potential scale reduction parameter was < 1.1 (Gelman and Rubin 1992), and visually inspecting density plots of posterior distributions of each parameter. We evaluated model fit using posterior predictive checks, i.e., Bayesian P-values (Royle and Dorazio 2008). Bayesian P-values from 0.2-0.8, with 0.5 as the optimal value, suggest an acceptable fit to the data. We summarized predictive ability of models by calculating mean square error for each iteration and summarizing results. We report relative breeding bird density per km² by dividing predicted abundance by the area surveyed per cell-year-survey, and mean parameter estimates with 95% credible limits (CL) throughout. We used a statistical significance criterion of whether 95% credible intervals overlapped zero for fixed effects and comparisons of interest.

Population change in space and time
We evaluated spatially explicit population change in three ways. First, we mapped the geometric mean rate of change in predicted abundance per cell for cells with ≥ 3 non-zero counts over the 25year period. We restricted maps to these cells so that trends within cells were representative of the cell and not solely an artifact of smoothing. Second, we summarized density and trends for the entire ACP, as well as the NPR-A and the ANWR. Third, we evaluated whether population growth rate varied over the survey period by conducting a rolling-window analysis where we estimated the geometric mean rate of change over successive 10-^

Important areas
We derived indices of relative importance (RI) across the study period for individual species for each cell by calculating the standardized deviance of annual cell density to annual mean density and then averaging these deviances through time for each cell RI c = mean[(density c,t -density t )/density t ].
Calculating ratios by year ensured overall importance of a cell was independent of any temporal trend in population size. Index values equal to zero represented cells that had mean density equal to overall density. We calculated RI for each iteration of the posterior distribution of predicted density to obtain error estimates. We then classified cells as important = 1, average = 0, or unimportant = -1 by comparing whether 95% credible intervals around standardized deviance overlapped zero. We further summarized relative importance of each cell by guild and all analyzed species by summing importance scores, i.e., -1, 0, 1, across groups. We assigned cells never occupied by a species to RI = -1. This approach to estimating relative importance has two advantages: (1) using model-predicted density estimates and their associated error accounts for unequal survey effort and area surveyed among cell-years, and (2) each species had one value, ensuring more abundant species were not over-represented in maps of relative importance at broader taxonomic scales.

RESULTS
Our data set included 98,610 observations and 329,546 indicated breeding birds of the 20 species we evaluated. At least one species was observed in 1774 approximately 36 km² cells (mean = 32.08 km²; range 1.02-36 km²) that encompassed 56,904 km 2 of the 57,336 km² ACP survey area (99.2%; . The earliest year surveyed per cell ranged from 1992 to 2011, but this distribution was highly skewed with most sampling first occurring 1992-1995. Cells were surveyed 1-24 times each (mean = 13, mode = 19; Fig. 2b). The NPR-A was surveyed more intensely (mean = 15 times; range = 1-24) than the ANWR (mean = 6 times; range = 1-15). The mean onset of phenological spring was [11][12]SD = 6.14;, but the onset of spring advanced approximately 9 days between 1992 and 2016 (Appendix 2). The early-flown survey was conducted on average 0.6 (SD = 6.12) days after the predicted onset of spring (range: -26-25) while the average late-flown survey was conducted 14 (SD = 6.59) days after the predicted onset of spring (range: -15-37). In general, survey dates were better aligned with the SOST spring phenology index at higher latitudes closer to the coast (Fig.  3).

Model results
Posterior predictive plots and values suggested models converged for all species after issues were addressed for several species. We modeled 17 species using the negative binomial distribution and the remaining three species using the ZINB distribution (Appendix 1). We averaged predicted density over surveys for Scaup and Yellow-billed Loon because they nest frequently in the southern extent of the ACP, a region that was not surveyed during the early-flown survey 1992-2006. This led to a strongly positive later survey effect on abundance, which when applied to predictions for the early survey across the entire ACP, led to very low predicted abundance in early years. In general, breeding bird abundance declined with phenologically adjusted survey date and the probability a count was zero increased with relative survey date for the three species modeled using the ZINB distribution suggesting breeding pair abundance (or possibly their detection rate) decreased throughout the season (Table A4.1). Nevertheless, Tundra Swan, Scaup, Jaeger, and Pacific Loon showed an increase in abundance with relative survey date, i.e., positive A.DOY coefficients (Table A4.1). In general, predictive ability of models was consistent with the modeled distribution used for each species but was lower for Lesser Snow Goose than other species, possibly because of their highly clustered distribution, i.e., occurring in discrete colonies, and relatively recent colonization of the ACP (Tables A4.4-A4.24).
Total abundance varied by species (Fig. 4, Appendix 4) and Greater White-fronted Goose, Northern Pintail, Long-tailed Duck, and to a lesser extent, King Eider, Arctic Tern, and Redthroated Loon showed negative associations between population size and growth (Fig. 5). Predicted density varied among administrative units (Fig. 6, Table A4.2). Jaeger, Tundra Swan, Red-throated Loon, and Cackling Goose averaged higher density in the ANWR. Greater White-fronted Goose, Pacific Loon, and Sabine's Gull averaged lower densities in the ANWR. Further, White-winged Scoter, Yellow-billed Loon, and Steller's Eider were never observed in the ANWR.
Densities of 10 species were greatest near the coast including all goose species, Tundra Swan, Northern Pintail, Steller's and Spectacled Eider, and Red-throated Loon (Appendix 5). In contrast, White-winged Scoter, Scaup, and Red-breasted Merganser showed an affinity for the southern extent of the ACP (Appendix 5).
Population trends varied widely by species; 13 increased, 1 declined, and 6 were stable over the study period (Fig. 4). Lesser Snow Goose showed the largest annual growth rate (13.8%; 95% CL: 9.0, 18.4%) while Red-throated Loon declined at 3.4% per year (95% CL: -5.6, -1.3%; Table A4.2). Several species showed some evidence of variable rates of growth over the study period; over the entire study period Arctic Tern population trend was positive, but most recent 10-year windows, i.e., > 2002, suggest populations may be declining (Fig. 7). Similarly, several species showed a peak in population growth from the late 1990s to early 2000s, with decreases in recent years, e.g., Greater White-fronted Goose, Sabine's Gull, Glaucous Gull.
Trends were generally similar among administrative units (Fig. 8 inherently drive population change. Our study represents one of the first steps in evaluating patterns of avian population change across the ACP landscape by describing spatially explicit population trends of waterbirds. Our results support earlier studies that found the ACP supports large breeding waterbird populations (Bart et al. 2013). Furthermore, most species increased or were stable over the study period, except for Redthroated Loon and possibly, Long-tailed Duck, which appeared to experience long-term population declines. Within species, patterns of density and population change were heterogeneous across the ACP, and some stable species experienced areas of local decline.

Abundance and important areas
Our indices of abundance are generally consistent with other analyses of data from the ACP breeding waterfowl survey ). Steller's Eider abundance estimates, however, are higher than previous reports because we included data from the later-timed survey that may have included non-or failed-nesting birds that were not incorporated in other studies (Bart et al. 2013, Dunham and Grand 2017. The NPR-A contained areas of high density for most species (Appendix 5), which is consistent with Bart et al. (2013) who showed highest density in similar areas. Bart et al. (2013) concluded the ANWR had the lowest waterbird densities in the ACP based on aerial and ground survey data from 1992-2010. Conversely, we found greater than average densities of four species within the 1002 area of the ANWR. Densities in the ANWR were high for Cackling Goose, Tundra Swan, Jaeger, and Red-throated Loon (Fig. 6), although the wetland area and total population sizes were small compared to the NPR-A. Because surveys occurred relatively late in the ANWR (Fig. 3), compensating for late survey timing caused estimates of density to be higher compared to Bart et al. (2013). The area also had relatively low survey effort (Fig. 2b), increasing uncertainty in estimates.
Important waterbird areas were concentrated between Utqiaġvik and Prudhoe Bay ( Fig. 9), especially for waterfowl (Appendix 7,8), specifically, within the Teshekpuk Lake Special Area in the NPR-A, which has been identified as the most important breeding area for waterbirds on the ACP (Liebezeit et al. 2011, Andres et al. 2012 and is also an important molting area for Arctic-breeding geese (Derksen et al. 1982, Flint et al. 2008, Shults and Dau 2016. Previous work suggests the area just south of Utqiaġvik, Alaska supports high densities (> 0.75 birds per km²) of breeding Steller's, spectacled, and King Eiders at local scales (Obritschkewitsch and Ritchie 2017). Our results show relatively high densities of eiders in this area, but our importance index did not identify these areas as hotspots. Thus, important areas for these and possibly other species may occur at spatial scales smaller than we evaluated or were not identified because of relatively uniform densities throughout their range or high variance within cell predictions. Further, because our index relies on model-based estimates, estimates with greater precision, e.g., surveyed in more years, could be more likely to demonstrate importance (or lack thereof), as such, estimates of importance from the southern ACP and the ANWR, which was surveyed less intensely over our study period (Fig. 2b) should be interpreted with caution.

Population trends
Most waterbird populations on the ACP increased, which is generally consistent with other published studies of breeding population trends in the ACP and elsewhere (Appendix 9). All goose and swan populations increased on the ACP, especially in the NPR-A near Teshekpuk Lake (Fig. 4). Further, only Greater White-fronted Geese showed slowing of population growth in recent years (Fig. 5), suggesting patterns of increase for geese and swans may continue.
We observed several notable contrasts to other published population trends for some species. The population trends for Scaup (3.3%; 95% CL: 0.9%, 6.2%) and Northern Pintail (0.1%; 95% CL: -1.4%, 2%) differ from the broader continental pattern of gradual decline (USFWS 2018). Part of this divergence is that the ACP Scaup population consists mostly of Greater Scaup, and areas where Scaup comprise mostly Greater Scaup are experiencing less drastic declines than areas occupied by Lesser Scaup (Kessel et al. 2002, Ross et al. 2015. Additionally, Northern Pintail were previously shown to overfly the Midcontinent prairies in dry years when numbers in the Arctic would increase (Henny 1973, Derksen and Eldridge 1980, Nicolai et al. 2005. As such, the numbers seen in the Arctic are not necessarily independent of habitat conditions elsewhere. Multiple sea duck species have declined globally, but most species in the Pacific Flyway were stable or increasing through 2011 (Flint 2013). We found similarly mixed results among sea ducks in our study. White-winged Scoters increased > 10% per year, but the ACP represents the northern edge of their breeding range and the little information available on population trends suggests declines in some areas along their southern extent (Brown and Fredrickson King Eiders initially increased during our study but have stabilized in recent years (Fig. 4). Conversely, King Eiders decreased substantially on the North Slope from the 1970s to 1990s (Suydam et al. 2000) and a population model suggests a stable or slowly declining population in Alaska driven mostly by adult and duckling survival (Bentzen and Powell 2012).
The three loon species present on the ACP show the full range of population processes. The numerically dominant Pacific Loon appears stable, whereas the less numerous Red-throated and Yellow-billed Loons demonstrate declining and increasing trends, respectively. Red-throated Loons were the only species to show a substantial rate of decline (-3.4%yr -1 ) in our study, which is contrary to stable populations reported from recent data summaries for this survey (-0.5% yr -1 , 95% CL: -1.8%, 0.08% yr -1 ; ) and the Yukon-Kuskokwim Delta Aerial Breeding Waterfowl Survey (0.0%, 90% CL: -5%, 10%; Swaim 2017). Nevertheless, our results show a slowing of population decline through the study period ( Fig. 7), which is consistent with trends on the ACP and other parts of the state (Mallek and Groves 2011. These loon species have substantially different foraging behaviors and diets as well as migration and wintering habits. Therefore, variation in population trajectory for sympatric species could be driven by local forage availability (Rizzolo 2017) or factors outside the ACP.
Glaucous Gulls are increasing on the ACP whereas some populations have declined in the Canadian Arctic (Gaston et al. 2012, Petersen et al. 2015. These important predators of eggs and chicks have shown similar population growth patterns to cooccurring waterbirds (Barry and Barry 1990, Gilchrist and Robertson 1999). Previous studies have found positive associations between nesting gulls and Common Eiders (Somateria mollissima) at high latitudes (Götmark andAhlund 1988, Robertson andChaulk 2017). Similarly, mean cell-specific density of Glaucous Gulls was positively associated with densities Sabine's Gull and Arctic Tern populations were positively associated with each other (r = 0.57) and both species increased slowly during the study period. Concordantly, populations of Sabine's Gulls appear to be increasing elsewhere in Alaska, i.e., Yukon-Kuskokwim Delta  while Arctic Tern population trends may be declining in Arctic Canada Robertson 1999, Maftei et al. 2015). Our results further suggest both populations have showed slowed growth or declines over the most recent 10 years (Fig. 7), with Arctic Terns potentially experiencing density-dependent growth (Fig. 5).
Almost all species we examined showed spatial variation in population trend. For most broadly distributed species, population trends varied latitudinally, with five species increasing faster toward the coast, and four species showing the opposite trend (Appendix 6). Concurrently, trends in the onset of phenological spring over our study period were stronger in nearshore areas compared to interior areas, especially near the Chukchi Sea coast in northwestern ACP ( Figure A2.2). This area, one of the latest areas to warm, has shown the largest change in timing of green-up. These patterns suggest waterbirds may be responding differently to earlier springs and the magnitude and distribution of change in the arrival of spring on the ACP. Further research is needed to assess mechanisms underlying observed patterns of change in our study.

Study considerations
Our study employed several methods to address potential issues in data collection including spatio-temporal smoothing to account for intermittent sampling in space and time, weighting variances by area surveyed per cell-year-survey, adjusting predictions for phenologically adjusted survey date, and using an offset for the later-survey to account for some variability in detection rate among observers.
Survey timing relative to the onset of spring was important for 85% of species examined and varied in direction; the number of observations decreased throughout the season for all but four species (Table A4.1). Our results were consistent with trends in survey timing from previous summaries of these data (Larned et al. 2008(Larned et al. , 2009) but showed somewhat stronger associations for some species. Therefore, our results confirm the need to control and/or adjust for survey timing relative to breeding status when monitoring abundance of nesting birds with aerial surveys.  We were unable to formally account for imperfect detection of birds in our study, and as such, assumed that pilot-observer pairs had similar detection among years. To refine this assumption, we included an effect of the latter survey to account for potential observer-pilot pair effects within a year. In general, detection probability may vary by species, group size, time of day, habitat type, flight conditions (e.g., weather, sun glare), season date, observer, or aircraft (Cook and Jacobson 1979, Tracey et al. 2005, Pearse et al. 2008, Koneff et al. 2008, Sauer et al. 2014 Despite accounting for several sources of potential error, our inference was somewhat limited by sampling constraints. Inconsistent sampling effort in space and time may have masked some spatial patterns of density, occupancy, and core use at the spatial scale we used (6 km by 6 km). Increasing the spatial scale examined may illuminate additional patterns at the expense of sampling intensity within each sampling unit, i.e., less area surveyed per cell-year. Further, sampling effort was relatively low in the ANWR and the southern ACP and often these areas were surveyed phenologically late, leading to less power to detect population change and possibly lower or higher detection of some species. For example, Yellow-billed Loon maps suggest relatively high density in the south-central portion of the ACP (Appendix 5), despite low density of large lakes necessary for breeding loon pairs (North and Ryan 1989). Upon further examination, these cells were surveyed infrequently and contained < 3 nonzero counts Avian Conservation and Ecology 14(1): 18 http://www.ace-eco.org/vol14/iss1/art18/ across the study period. Thus, uncertainty in density for these cells is high, but not captured in density maps (cell-specific density SEs are available from the corresponding author). We removed cells with < 3 nonzero counts across the study period from trend maps to avoid making conclusions regarding population trend in areas with low sampling effort but did not apply that criterion to density maps to illustrate the full distribution of each species on the ACP. Regardless, increased sampling effort in these areas is needed to improve estimates and resulting maps. Third, we attempted to reduce counts of nonbreeding birds or birds temporarily in an area before dispersing to breeding locations by only including counts of pairs or indicated pairs. This approach may underestimate population size to some extent for many species, particularly colonial species because recording of pairs versus flocks may have been inconsistent among observers and through time. For example, an observation of birds may have been recorded as a flock of eight or four pairs, only the latter of which would have been included in our analysis. Because we only included individuals or pairs, our indices of abundance may be lower than similar studies and trends conservative for some colonial species like Lesser Snow Geese (see Burgess et al. 2017. We recommend observers create standardized protocol to ensure breeding status of groups > 2 are consistently recorded.
Because Arctic waterbirds are migratory, population regulation could occur at other times in their annual cycle. Changes in survival or body condition, e.g., for capital breeders, during the nonbreeding season may influence abundance or breeding success the following year. If so, the resulting change on the ACP is more likely to be uniform and widespread than spatially restricted, a pattern that emerged for several species including Pacific Loon, Northern Pintail, and Long-tailed Duck (Appendix 6). Lifehistory traits make some species more or less susceptible to environmental conditions in the Arctic (Saalfeld and Lanctot 2017). For example, all the species we evaluated are site faithful to breeding areas, but 17 species, i.e., all but Scaup, Red-breasted Merganser, and Northern Pintail, also exhibit delayed reproduction. That is, for most of these species, surviving young may not return to breeding grounds for two to five years and breeding propensity of adults is relatively low and variable (Sibley 2009, Elphick 2016. Therefore, effects of local breeding conditions on reproduction, survival, and dispersal may result in a lagged response in terms of population change. During that time, survival rates of nonbreeding adults will vary based on conditions away from breeding areas, thus masking the effects of local breeding-area conditions on population growth. Population models that incorporate data estimating vital rates such as survival or productivity could help detangle the source, e.g., breeding or nonbreeding areas, and mechanisms underlying observed change for these species.

Future directions
Although changes to physical environmental drivers are documented and modeled well into the future, responses of avian populations are mediated through complex interactions among life-history traits, environmental conditions, and trophic interactions that are not homogeneous across the landscape or the annual cycle. Because of this, accurate prediction of future population dynamics is difficult (Van Hemert et al. 2015). Our study is the first step to identifying important areas for breeding waterbirds on the ACP and areas exhibiting rapid population change based on long-term, comprehensive survey data of multiple species.
Because the ACP landscape is changing as a result of climatedriven environmental change and anthropogenic disturbance, future research is necessary to identify drivers of observed patterns of change, quantify their effects, and forecast population change into the future. Logical next steps include (1) examining patterns in abundance and trends relative to vegetation and wetland composition to identify environmental characteristics associated with avian density and population change in the Arctic, (2) examining population trends relative to life history traits and nonbreeding ecology that may identify traits or areas outside the Arctic influencing population growth, (3) assessing abundance and trends relative to existing oil and gas infrastructure to inform future leasing decisions, and (4) identifying areas of landscape change, e.g., thermokarst events or vegetation change, in the Arctic to determine wildlife responses and predict how populations will continue to change into the future.
Responses to this article can be read online at: http://www.ace-eco.org/issues/responses.php/1383

Spatio-temporal population change of arctic-breeding waterbirds
We evaluated the ability of our model structure to accurately represent the number of zeros in the data by deriving 3,000 realizations of count data using Monte Carlo simulations of converged model parameter estimates and associated error and comparing the resulting proportion of zeros in the realized counts with the proportion of zeros observed in the real data ( Figure 1). Results suggest negative binomial models were sufficient for most species, but underestimated the proportion of zeros for Greater White-fronted Goose, Northern Pintail, and Long-tailed Duck by to 0.02 to 0.05. These species were the most abundant and broadly distributed in our study and thus, higher mean abundance may have prevented the level of overdispersion necessary in the negative binomial distribution to accurately reflect the number of zeros in the data. Further examination of the spatio-temporal trends in bias showed no pattern in bias through time; estimates were generally consistently lower for realized data ( Figure 2). Patterns in space showed that locations of disagreement were generally peripheral to high density areas ( Figure 3). The negative binomial Northern Pintail model could adequately predict the number of zeros in the early survey replicate, but not the later survey ( Figure 4). As a result of this analysis we reran the three biased species using a zero-inflated negative binomial distribution (ZINB) (see main text).   Purple denotes areas where realized data and actual data both predicted non-zero counts, tealblue areas are areas of disagreement between realized and actual data (i.e., only one group predicted a zero count), and yellow denotes areas where both realized and actual data resulted in zero counts. We report the distribution of realized values and the actual proportion of zeros as a dashed line for each survey replicate (early, late).

Spatio-temporal population change of arctic-breeding waterbirds
We obtained data on the predicted start of the growing season (start of season time, SOST; day of the year) from vegetation indices derived from MODIS satellite imagery from 1992 to 2016 , Didan et al. 2016). Data are available at approximately 2.5 km by 5 km spatial resolution across the Arctic Coastal Plain, Alaska (ACP). We summarized SOST by cellyear by calculating the mean, weighted by area covering each cell for cells with > 1 value. Some SOST values are not well predicted along the coast and particularly in river deltas subject to spring flooding. Therefore, we removed SOST values in cells that were < 122 (2 or 3 May) or > 228 (16 or 17 Aug) and imputed values for those cell-years by interpolating values from the 40 nearest neighbors in each year weighted by the inverse-weighted distance to the target cell. In most cases, this represented cells < 3 cell-widths or 17 km from the target cell boundary (i.e., on the diagonal). The days chosen as cutoffs represented natural break-points in the data. We used the knnimputation function in the R package DMwR (Torgo 2015, R Core Team 2018) to impute omitted values.
We evaluated trends in spring phenology through space and time in two ways. First, we evaluated a generalized additive model (GAM; Wood 2017) that incorporated spatial autocorrelation and an additive temporal trend to estimate changes in SOST through time. Second, we evaluated a space by year interaction to examine spatial patterns of temporal change in SOST. We implemented GAMs using the gam function in the R package mgcv (Wood 2018). We obtained confidence intervals around differences in SOST between 1992 and 2016 using Monte Carlo simulations. In this appendix, we provide R code and results to facilitate replication of our analyses.

RESULTS
Model selection results supported an interaction between space and time, and a smoothed spatial component (Table A2.  Figure A2.1. Predicted start of the growing season (day) from 1992 to 2016 for 1,914 6 km by 6 km cells across the Arctic Coastal Plain, Alaska. We modeled the effect of space and time on start of growing season day using generalized additive models that apply a non-linear smooth to spatio-temporal effects. The gray shaded region denotes 95% confidence interval around predictions. The lower plot is identical, but includes semi-transparent data values across cells for each year (in blue).

Figure A2.2.
The change (in days) in the start of spring growing season (day of year) between 1992 and 2016 within approximately 36 km 2 cells (6 km by 6 km) across the Arctic Coastal Plain, Alaska. We modeled the effect of space and time on start of growing season day using generalized additive models that apply a non-linear smooth to spatio-temporal effects, then subtracted cell-specific values in 1992 from 2016. Negative values denote a decreasing trend in the start of spring (i.e., start of the growing season occurred earlier in recent years).