Harbor Seals as Sentinels of Ice Dynamics in Tidewater Glacier Fjords

Tidewater glaciers calve icebergs into the marine environment which serve as pupping, molting, and resting habitat for some of the largest seasonal aggregations of harbor seals (Phoca vitulina richardii) in the world. Although they are naturally dynamic, advancing and retreating in response to local climatic and fjord conditions, most tidewater glaciers around the world are thinning and retreating. Climate change models predict continued loss of land-based ice with unknown impacts to organisms such as harbor seals that rely on glacier ice as habitat for critical life history events. To understand the impacts of changing ice availability on harbor seals, we quantified seasonal and annual changes in ice habitat in Johns Hopkins Inlet, a tidewater glacier fjord in Glacier Bay National Park in southeastern Alaska. We conducted systematic aerial photographic surveys (n = 55) of seals and ice during the pupping (June; n = 30) and molting (August; n = 25) periods from 2007 to 2014. Object-based image analysis was used to quantify the availability and spatial distribution of floating ice in the fjord. Multivariate spatial models were developed for jointly modeling stage-structured seal location data and ice habitat. Across all years, there was consistently more ice in the fjord during the pupping season in June than during the molting season in August, which was likely driven by seasonal variation in physical processes that influence the calving dynamics of tidewater glaciers. Non-pup harbor seals and ice were correlated during the pupping season, but this correlation was reduced during the molting season suggesting that harbor seals may respond to changes in habitat differently depending upon trade-offs associated with life history events, such as pupping and molting, and energetic costs and constraints associated with the events.


INTRODUCTION
The distribution of abiotic and biotic resources required for animal survival and reproduction is seasonally dynamic and can change in space and time (Fretwell, 1972). As such, organisms living in highly seasonal environments may adapt by timing life history events, such as reproduction and migration, to coincide with favorable environmental conditions (Stearns, 1989). Furthermore, there can be complex interactions between life history events, the physiological state of the animals, and the abiotic and biotic environment which can result in trade-offs during the annual cycle (Stephens et al., 2014). Thus, elucidating how environmental factors, such as habitat and prey availability, interact with intrinsic factors, such as physiological state, and influence the annual cycle of individuals (McNamara and Houston, 2008;Varpe, 2017) is essential for understanding and predicting how populations may respond to changes in rapidly changing subpolar and polar ecosystems Pörtner et al., 2019).
Tidewater glaciers are mountain glaciers that flow to the ocean and calve icebergs into the marine environment (Figure 1). These marine-terminating glaciers are prominent physical features in high-latitude environments that play an important role in linking landscape and marine ecosystem processes in Greenland, Iceland, Svalbard, the Canadian and Russian Arctic, Alaska, Chile, and Antarctica Straneo et al., 2019;Bianchi et al., 2020). Tidewater glaciers undergo substantial seasonal variation which can influence freshwater flux, vertical mixing, nutrient input (Bartholomew et al., 2010;Howat et al., 2010;Moon et al., 2014;Arimitsu et al., 2016;Amundson and Carroll, 2017;Fried et al., 2018), biological productivity, marine food webs, and create important habitats for invertebrates, fish, seabirds, and marine mammals during critical phases of the annual cycle (Etherington et al., 2007;Lydersen et al., 2014;Juul-Pedersen et al., 2015;Meire et al., 2017). For example, subglacial discharge plumes transport nutrients and entrain zooplankton near the glacier terminus which creates productive foraging areas for black-legged kittiwakes (Rissa tridactyla), northern fulmars (Fulmarus glacialis), and other seabirds Urbanski et al., 2017). Belugas (Delphinapterus leucas) and narwhals (Monodon monoceros) use habitats near glacial fronts and tidewater glaciers for foraging (Lydersen et al., 2014;Lucey et al., 2015;Laidre et al., 2016). Similarly, ice that is discharged or calved from the terminus of tidewater glaciers produces icebergs that serve as resting, pupping, and molting habitat for several species of pinnipeds including ringed seals (Pusa hispida), bearded seals (Erignathus barbatus), and harbor seals (Phoca vitulina) (Calambokidis et al., 1987;Lydersen et al., 2014;Hamilton et al., 2016Hamilton et al., , 2019. Harbor seals (Phoca vitulina) are the most widely distributed pinniped in the northern hemisphere and some of the largest seasonal aggregations of harbor seals in the world are associated with floating ice habitat in tidewater glacier fjords in Alaska . Tidewater glacier fjords in Alaska range from LeConte Glacier in southeastern Alaska to Kenai Fjords in south-central Alaska, spanning approximately 1,300 km from 56 • to 61 • latitude. Harbor seals exhibit a high degree of seasonal fidelity to tidewater glacier fjords during energetically demanding life-history phases of the annual cycle (Womble and Gende, 2013). Use of ice habitat by harbor seals likely confers benefits related to providing a stable platform for nursing pups during the brief lactation period and reducing the risk of predation (Calambokidis et al., 1987;Blundell et al., 2011). Harbor seals are also an important cultural and subsistence resource for Alaska Natives (Crowell, 2016) and a highly sought-after viewing experience for visitors to tidewater glacier fjords in Alaska (Jansen et al., 2010;Young et al., 2014).
Although tidewater glaciers are naturally dynamic, advancing and retreating in response to physical drivers and local climatic and fjord conditions (Post et al., 2011), most glaciers around the world are thinning and as a result, many tidewater glaciers are retreating (Arendt et al., 2002;Larsen et al., 2007;Wouters et al., 2019;Zemp et al., 2019). Although models predict continued mass loss for glaciers globally, including tidewater glaciers (Hock et al., 2019;Slater et al., 2019), there is a limited understanding of how organisms, such as harbor seals, that utilize ice habitat in tidewater glacier fjords may respond. However, studies elsewhere in polar regions have demonstrated that changes in sea ice may influence pup survival, disease risk, foraging behavior, distribution, and habitat use of ice-associated pinnipeds (e.g., Kelly, 2001;Ferguson et al., 2005;Johnston et al., 2005;Laidre et al., 2008;Bajzak et al., 2011;Kovacs et al., 2011Kovacs et al., , 2020Hamilton et al., 2015).
One of the largest seasonal aggregations of harbor seals in southeastern Alaska occurs in Glacier Bay (Calambokidis et al., 1987), an expansive tidewater glacier fjord that has retreated over 100 km, and lost more than 3,000 km 3 of ice since 1770 (Cooper, 1937;Field, 1947;Larsen et al., 2005). In addition, the number of actively calving tidewaters glaciers in Glacier Bay has decreased from 12 in 1982 (Molnia, 2007) to five in recent years. Over the last several decades, the distribution of harbor seals in Glacier Bay has also changed significantly due in large part to the rapid retreat and grounding of the Muir Glacier in the early 1990's (Hall et al., 1995), where historically over 1,300 seals used icebergs as habitat in the 1970's and 1980's (Streveler, 1979;Calambokidis et al., 1987). Since the grounding of the Muir Glacier, the only remaining ice habitat in Muir Inlet is at the McBride Glacier which has retreated over 3 km from 1948 to 2012 (McNabb andHock, 2014) and provides habitat for up to 200 seals (Womble et al., 2010).
Currently, over 75% of harbor seals in Glacier Bay occur in Johns Hopkins Inlet, a tidewater glacier fjord that is fed by the Johns Hopkins and Gilman Glaciers in the West Arm of Glacier Bay. After undergoing a retreat that began at the end of the nineteenth century, the Johns Hopkins Glacier has advanced nearly 2 km since the mid-twentieth century (McNabb and Hock, 2014). Though the terminus of the Johns Hopkins Glacier has remained relatively stable since 2012, over the last few decades (1992-2017) the abundance of seals in Johns Hopkins Inlet has declined precipitously (Mathews and Pendleton, 2006;Womble et al., 2010Womble et al., , 2020. Collectively, these relatively largescale changes in tidewater glaciers and the distribution of seals in Glacier Bay suggest that ice is a key environmental variable that may influence the ecology and habitat use of seals. Ultimately, harbor seals may serve as sentinels (e.g., Moore, 2008;Hazen et al., 2019) of ice dynamics in these rapidly changing tidewater glacier fjord ecosystems.
To date, few studies have quantitatively addressed the relationships between harbor seals and ice habitat (e.g., Jansen et al., 2015), due in large part to the expansiveness and remote nature of tidewater glacier fjords as well as the dynamic nature of the ice habitat (Boveng et al., 2003;Bengtson et al., 2007). However, recent advances in survey and analytical methods allow for the systematic quantification of seals and ice using aerial photographic surveys which provide a permanent record McNabb et al., 2016). If tidewater glaciers continue to thin and retreat, the amount of ice habitat available for harbor seals may decrease and seals may spend more time in the water, use terrestrial sites, or move to other areas (Calambokidis et al., 1987;Womble et al., 2010;Hoover-Miller and Armato, 2018). Given the reliance of seals on ice habitat for critical life history events, such as pupping and molting, each of these outcomes could have population-level consequences.
Our primary objective was to quantify how changes in ice habitat influences harbor seal distribution and abundance during two energetically demanding life history phases of the annual cycle: the pupping period in June and the molting period in August. Specifically, we (1) conducted aerial photographic surveys to quantify the spatial distribution and abundance of harbor seals and ice from 2007 to 2014; (2) quantified the availability and characteristics of ice from aerial photographs using object-based image analysis; and (3) developed multivariate conditional autoregressive models to jointly model stagestructured animal location data and ice habitat to describe the spatio-temporal distribution of ice and seals in Johns Hopkins Inlet. Given that ice is a dynamically evolving resource that we were interested in understanding jointly with seal distribution, we used statistical methods that permit joint modeling of animals and habitat to better understand harbor seal and glacier ice dynamics.

Study Area
Johns Hopkins Inlet (58 • 50.89 N, 137 • 06.12 W) (Figure 2) is an expansive (12 km long × 3 km wide) tidewater glacier fjord in the upper West Arm of Glacier Bay National Park, a Biosphere Reserve and World Heritage site encompassing over 242,811 hectares of marine waters in southeastern Alaska. Harbor seals aggregate seasonally in Johns Hopkins Inlet to rest, pup, and molt on icebergs that calve from two tidewater glaciers (Mathews and Pendleton, 2006;Womble et al., 2020); the Johns Hopkins (250 km 2 ) and the Gilman (25 km 2 ) which extend from the Fairweather Mountain Range.

Aerial Photographic Surveys of Seals and Ice (2007-2014)
We conducted aerial photographic surveys (n = 55) of harbor seals and ice in Johns Hopkins Inlet (Figure 2) during the pupping period in June (date range: 12-30 June; n = 30 surveys) and during the molting period in August (date range: 2-25 August; n = 25 surveys; Table 1) from 2007 to 2014. Surveys were conducted between 1200 and 1800 Alaska Daylight time, as higher counts of seals typically occur 1-4 h after solar noon (Mathews and Pendleton, 2006).
Aerial photographic surveys were conducted from a de Havilland Canada DHC-2 Beaver single-engine high-winged aircraft (Ward Air Inc., Juneau, Alaska) following methods developed by Jansen et al. (2006) and Ver Hoef and Jansen (2015). Aerial surveys were flown at approximately 1,000 feet (305 m) at 90-95 kts along 12 established transects. Transects were spaced 200 m apart and spanned the length of the fjord from the terminus of the glacier to the opposite end of the fjord (Figure 2). The transects encompassed an area of approximately 10.8 km 2 of the 22.5 km 2 surface area of the water in Johns Hopkins Inlet. An aerial survey of Johns Hopkins Inlet was completed in approximately 1 h (Womble et al., 2020).
During the aerial surveys, non-overlapping digital photographic images (Figure 3) were taken directly under the plane using a vertically aimed digital single-lens reflex (DSLR) camera (Nikon D2X, 12.4 megapixel; Shinagawa, Tokyo, Japan) with a 60 mm focal length lens (Nikon AF Micro-NIKKOR, 2.8D). The camera was attached to a tripod head and mounted to a plywood platform that was secured in the belly porthole of the aircraft. The camera captured an image every 2 s, using a digital timer (Nikon MC36) operated by an observer. The firing rate and the spacing of the transects allowed for a gap between images of approximately 15 m end-to-end and 70 m side-to-side to ensure that images were separated from one another and so seals were only sampled once. Each digital photo (3,216 × 2,136 pixel JPEGS) (Figure 3) covered approximately 80 × 120 m at the surface of the water with ∼3.7 cm pixel resolution. An onboard global positioning system (Garmin 76 CSX) was used to record the track line and position of the plane along the transects at 2 s intervals (Womble et al., 2020).

Post-processing of Aerial Digital Imagery for Seals
The latitude, longitude, and altitude from the track line were written to the EXIF headers of each digital image using RoboGEO V6.3 (Pretek, Incorporated, Christian, TN, United States). Images from each survey were embedded as a raster layer in an ArcGIS project using ArcGIS (ESRI, version 9.3 and 10.1) and R (R Core Team, 2018). Each photo was examined by a trained observer using digital photographic software (ACDSEE Pro 4) and each seal was marked as a point feature in an ArcGIS shapefile. After all images were reviewed and all seals were marked, the point locations for seals were summed within each image and assigned to the centroid of each photo and exported as shape files for statistical analysis (Womble et al., 2020).

Post-processing of Aerial Digital Imagery for Ice
We used object-based image analysis (Blaschke, 2010;Blaschke et al., 2014) to quantify the amount and characteristics of ice habitat from the high-resolution digital imagery collected during aerial photographic surveys. Specifically, we used Trimble eCognition Developer (version 8.9.0) (Trimble Geospatial Imaging) to classify and extract multiple variables from each image including: (1) iceberg (%) defined as the percent of each photo that is icebergs greater than 1.6 m 2 ; (2) brash ice (%) defined as the percent of each scene that is ice smaller than 1.6 m 2 ; (3) water (%) of each scene that is water; (4) iceberg size (m 2 ) defined as average size of icebergs in each scene; and (5) distance to the terminus of the glacier (km) defined as the distance from the glacier calving face to the center point of each scene or image (McNabb et al., 2016). Minimum iceberg size was based on the premise that an iceberg of approximately 1.6 m 2 would be large enough to support a non-pup seal whereas icebergs smaller than 1.6 m 2 would be too small to support a seal and were thus classified as brash ice. The minimum iceberg size was based upon the average curvilinear length of non-pup harbor seals (1.36 ± 0.15; range 1.00-1.76 m; n = 81 non-pup seals) that were live-captured, measured, and released in Johns Hopkins Inlet from 2004 to 2008 (Blundell et al., 2011;Womble and Gende, 2013).

Statistical Methods
Complete documentation of the statistical analysis procedures, including R Statistical Software (R Core Team, 2018) code, is provided in Supplementary Appendix 1. Preliminary analyses and model checking results suggested spatial autocorrelation in the residuals, violating the conditionally independent error assumption in standard linear regression models. Therefore, we used a multivariate conditional autoregressive (MCAR) model to introduce multiple, dependent spatial random effects associated with areal units (Leroux et al., 1999;Gelfand and Vounatsou, 2003;Banerjee et al., 2014). The MCAR model allowed us to summarize and display graphically the spatial distribution of our response variables, while controlling for interrelations among response variables more accurately than extrapolating response variables independently. To conduct our analyses, we discretized Johns Hopkins Inlet into 62 × 7 equally sized grid cells (i = 1,. . .,434 total grid cells), each 200 × 300 m in resolution. We used an MCAR model because we were interested in l = 1,2,3 response variables that included counts of adults, counts of pups, and ice proportion, and correlations among them. We centered and scaled each response variable in each k = 1,. . ., K ij photograph, taken in grid cell i, during time j = 1,. . ., J using z-standardization with respect to all photographs used in our study. We assumed where y ijk ≡ (adults ijk , pups ijk , ice ijk ) are the centered and scaled counts of adults, pups, and ice proportion from the kth photograph in the ith grid cell during the jth time. The vector µ ij represents the expected values of the transformed response Frontiers in Marine Science | www.frontiersin.org variables in each grid cell i during time j, and the covariance matrix V describes variance and covariance of response variables within a grid cell. We assumed where x ij represents a vector of covariates for grid cell i and β ij are the associated coefficients to be estimated. We used the same design matrix X for each l response variable, but parameter estimates β ij were allowed to vary among responses. We used covariates associated with distance to the terminus of the glacier. The distance covariate we used consisted of the linear distance from each cell to the terminus of the Johns Hopkins Glacier. The covariance matrix ≡ diag(W1)−ρW −1 ⊗ characterizes the multivariate and spatial covariance, where (diag(W1)-ρW) −1 is the conditional autoregressive (CAR) prior for the spatial correlation, with adjacency matrix W, diagonal matrix diag(W1) consisting of the number of neighbors on the diagonal, and spatial correlation parameter ρ. The matrix describes the covariance of the multivariate response. We used inverse Wishart priors for V and , with hyperparameters equal to 4 and I 3×3 , a 3 × 3 identity matrix. We used multivariate normal priors for β, with mean 0 and covariance matrix σ 2 β I pxp , where σ 2 β = 10 2 . We fit our model to the data using the MVS.VARleroux function in the CARBayes R package (Lee, 2017). We drew 1,050,000 iterations from the MCMC algorithm, discarding the first 50,000 of the original MCMC samples as burn-in samples, and thinned the remaining samples by retaining every tenth sample. We used the posterior median of µ, which is optimal with respect to an absolute error loss function (Williams and Hooten, 2016), to summarize and plot expected counts of non-pups, pups, and expected ice proportion in space and time.
We assessed trace plots for each full conditional distribution to assess convergence to the approximated marginal posterior distribution. All parameters appeared to converge; this observation was supported by Gelmin-Rubin diagnostics < 1.05 for all parameters. We assessed the goodness of fit of our model by simulating data from the posterior predictive distribution and comparing our simulated data to the observed data. We compared predicted data to observed data for each response variable by comparing the deviance discrepancy function for observed and predicted data and calculating Bayesian p-values (Conn et al., 2018). Bayesian p-values ranged from 0.11 to 0.67, suggesting no lack of model fit. Many covariates not accounted for in our process model likely contribute to ice and harbor seal spatial dynamics. However, because our model was able to reliably reproduce our observed data, and the associated uncertainty in our observed data, we considered our process model a parsimonious approximation of reality.
Abundance estimates for harbor seals during the pupping (non-pups and pups) and molting (non-pups) periods were generated using counts of seals from aerial photographs and predictions from unsampled areas using a model-based estimator following methods developed by Ver Hoef and Jansen (2015).

RESULTS
Across all years from 2007 to 2014, there was substantial seasonal and interannual variation in the amount of ice and the spatial distribution of ice in Johns Hopkins Inlet. On average there was 7.8 times more ice in the fjord during June than during August. However, there was variability across years, which ranged from 1.3 times more ice in June than August in 2014 to 18.4 times more ice in June than August in 2012. Two years, 2011 and 2012, stood out as particularly icy years (Figure 4). Across all years and during both June and August, there was consistently more ice concentrated near the terminus of the glacier than farther from the glacier (the slope parameter for ice β > 0) and the distribution of ice was also more consistent near the terminus of the glacier (Figure 5). The spatial distribution of ice in the same month was consistent among years, particularly during August (Figure 5). However, the distribution of ice in June was much more variable than in August (Figure 6).
There was substantial seasonal variation in seal distribution and abundance between the pupping and molting periods ( Figure 5). During June, seals were typically dispersed more extensively throughout the fjord and generally tracked ice distribution (Figure 7). However, there was substantial variation among years, with seal distribution extending from the near the terminus of the glacier to Jaw Point during June in most years, except for 2008 through 2010. In contrast, during the molting period in August, the distribution of seals was typically clustered much closer to the terminus of the Johns Hopkins Glacier. However, in 2011 and 2012 when there was a more extensive distribution of ice, the distribution of seals extended farther down the inlet away from the terminus (Figures 7, 8). The estimated abundance of seals (non-pups and pups) was greater during the pupping season than during the molting season in 2007, 2009, 2010) and corresponded to increased availability of ice and more extensive spatial distribution of ice in the fjord during June. However, in other years, estimated abundance of seals (non-pups and pups) was more similar between the pupping and molting seasons (Figure 9). During June when pups were still dependent upon adult females for nutrition, the response variables for non-pups and pups (i.e., the standardized counts of pups and non-pups at each site) were highly correlated (λ 1 , 2 = 0.85; range: 0.75-0.94) (Figure 10). Overall, there was a positive correlation between non-pups and ice during both the pupping (λ 1 , 3 = 0.27; range: 0.23-0.31) and molting (λ 1 , 3 = 0.22; range: 0.17-0.27) periods (Figure 11). However, there was variability in the correlation across seasons and years. As expected, given the strong correlation between non-pups and pups, the correlation between pups and ice was λ 2 , 3 = 0.30 during June (Figure 10). There was no correlation between non-pups and iceberg size, during the pupping or the molting season. Across all years, the proportion of pups during June averaged 0.31 and ranged from 0.24 to 0.36. Summaries for each marginal posterior distribution approximated using the MCMC algorithm are provided in Supplementary Appendix 1.

DISCUSSION
Understanding how ice habitat changes is important for understanding patterns in the distribution and abundance of harbor seals, particularly given the expected future changes to tidewater glaciers and associated fjord ecosystems. Across all years there was substantial seasonal and interannual variation in ice dynamics and seal distribution and abundance. There was consistently more ice and the spatial distribution of ice was more extensive during the pupping period in June than during the molting period in August. During the pupping period, when there was more ice and the distribution of ice was more extensive, non-pup harbor seals and glacier ice were more strongly correlated; however, this correlation was reduced during the molting season.
Animals living in dynamic environments may respond to changes in habitat differently depending upon trade-offs associated with life history events and the energetic costs and constraints associated with these events (Stearns, 1989;Stephens et al., 2014). In most vertebrates, reproduction and molting do not overleap due to the energetic costs associated with these life history events (Beltran et al., 2018). The reproductive period of mammals is energetically costly, with lactation being the costliest aspect of mammalian reproduction (Gittleman and Thompson, 1988). In Alaska, harbor seal pups are born from mid-May through June (Jemison and Kelly, 2001;Mathews and Pendleton, 2006) and are dependent upon adult females for energy, via lactation, for 3-5 weeks until they are weaned. Although some adult female harbor seals may forage during the lactation period   (Boness et al., 1994), for the most part, harbor seals are capital breeders (e.g., Jönsson, 1997) and rely on stored energy to nourish their young during the brief lactation period. Given that energetic demands are high during the reproductive period, having access to ice habitat in tidewater glacier fjords may confer several benefits including providing a stable substrate for nursing pups and resting, an isolated floating platform for naïve pups that reduces the risk of predation from terrestrial and aquatic predators, and a substrate that reduces the risk of disease and parasite transmission (Fay, 1974;Calambokidis et al., 1987;Blundell et al., 2011). In addition, the stable ice substrate likely provides thermoenergetic benefits for young pups that have a limited blubber layer . In contrast, terrestrial haulout sites that are used by seals are subject to tidal fluctuations multiple times per day, may expose seals to increased risk of predation, and in some cases, can become space limited, particularly during high tides when the haulout substrate may become submerged. In most years the correlation between seals and ice was reduced during the molting period compared to the pupping period. Although molting is energetically demanding (Thometz et al., 2021), costs associated with reproduction, particularly lactation, are greater (Gittleman and Thompson, 1988). Adult female harbor seals typically molt after the reproductive period and after their pups have been weaned, hence they are no longer constrained by the presence of nutritionally dependent young. The timing of molt is such that pups molt first, in utero (Boulva, 1975), followed by juvenile harbor seals that are not reproductively active, with adult female and adult males molting last (Daniel et al., 2003). Similarly, studies of Weddell seals (Leptonychotes weddellii), an ice-obligate phocid that is associated with the fast ice in Antarctica, have demonstrated that parturient females may delay molting activities to recover from the energetically demanding aspects of reproduction. In addition, the molting phenology of Weddell seals also covaried with the timing of sea ice break-out (Beltran et al., 2019). Collectively, this suggests that some phocids may exhibit plasticity in timing of lifehistory events in dynamic and highly seasonal habitats which may result in shifts in phenology across the annual cycle.
Several other lines of evidence also support the premise that harbor seals in tidewater glacier fjords modify their behavior and habitat use throughout the annual cycle, due in part to trade-offs associated with life history events and environmental factors, such as prey availability and habitat. Satellite telemetry studies have demonstrated that during the more energetically demanding pupping period, seals exhibit high fidelity to ice habitat in Johns Hopkins Inlet, with up to 78% of seals returning to the ice habitat from the previous year (Womble and Gende, 2013). In addition, seals tend to have reduced travel distances to foraging areas and spend more time out of the water during the reproductive and molting periods. While there may be benefits associated with using ice habitat during energetically demanding periods of the annual cycle, there may also be tradeoffs. When using ice habitat, harbor seals may expend greater effort to forage by diving deeper and/or traveling farther to forage in more shallow areas where prey may be more easily accessible (Womble et al., 2014). In contrast, from September to April, after reproduction and molting have been completed and seals are not as constrained, fidelity to ice habitat is reduced, seals travel much more extensively outside of tidewater glacier fjords, and the proportion of time that seals spend out of water is reduced (Womble and Gende, 2013). Similarly, isotopic signatures from seals in tidewater glacier fjords suggest that during the reproductive period, the foraging niche of harbor seals is narrower and more focused on pelagic prey species (Blundell et al., 2011;Smith et al., 2019). In contrast, during the nonreproductive period, when harbor seals travel more widely and have lower fidelity to ice habitat, there is a shift to a more diverse diet with an increased emphasis on benthic fish species (Smith et al., 2019).
There was consistently more ice in the fjord during the pupping season in June than during the molting season in August, which was likely influenced by seasonal variation in physical processes that influence the calving dynamics of tidewater glaciers (Fried et al., 2018). In Disenchantment Bay, a tidewater glacier fjord in the eastern Gulf of Alaska near Yakutat, Jansen et al. (2015) also documented maximum ice coverage in June during the peak of the harbor seal pupping period and decreased ice coverage in August. In Alaska, rates of frontal ablation of tidewater glaciers tend to be higher in spring and early summer (March-May) and coincide with peak surface velocities of glaciers and increased ice supply to the terminus of the glacier (McNabb et al., 2015). When frontal ablation is greatest during the spring and early summer, iceberg production tends to be increased, thereby resulting in more ice habitat in the fjord for seals. In contrast, during the molting period in August, there is typically less frontal ablation resulting in less ice habitat for seals and more spatially clustered aggregations of seals near the terminus of the glacier, where most of the ice is concentrated. Frontal ablation of tidewater glaciers, which includes mass loss via calving and submarine melt, can be driven by multiple interacting processes including submarine melting via elevated ocean temperatures, iceberg calving, convection through the mixing of warmer sea water with meltwater from the terminus, and interactions with subsurface geometry and depth of the fjord (Motyka et al., 2003;Bartholomaus et al., 2013;Sutherland et al., 2019;Jackson et al., 2020). Factors such as tides, currents, katabatic winds, air and ocean temperatures, and the presence and depth of sills in the fjord may also influence the distribution and persistence of ice within the fjord (Bartholomaus et al., 2015;Spall et al., 2017;Amundson et al., 2020).
There was an increased amount of ice in the fjord during 2011 and 2012 with a corresponding increase in seal distribution throughout the fjord during June and August of both years. It is unknown what factors may have contributed to the increased amount of ice during these years; however, inspection of the frontal position of the Johns Hopkins glacier in 2011 from satellite imagery suggests that there was approximately 40 m of length lost along the terminus between late April and early June, which likely resulted in increased ice calving and iceberg production. Furthermore, when coupled with cooler air temperatures, ice could have persisted for longer in the fjord and resulted in the increased availability of ice for harbor seals during 2011 and 2012. Hence, a more mechanistic understanding of ice availability in the fjord will require understanding not only frontal ablation and ice flux, but also the influence of environmental factors such as water and air temperatures on iceberg melting and persistence in the fjord.
While our results focus on one tidewater glacier fjord, previous studies in Glacier Bay (Calambokidis et al., 1987;Mathews and Kelly, 1996;Young et al., 2014) and in other fjords in Alaska also suggest that ice is a key environmental variable that influences seal distribution and abundance (Bishop, 2011;Jansen et al., 2015;Mathews et al., 2016;Hoover-Miller and Armato, 2018). Our results have implications for the timing of surveys of seals in tidewater glacier fjords, most of which occur during the molting period from late July through early September in Alaska, when ice habitat tends to be reduced. Although there may be substantial variability in tidewater glaciers and factors that influence glacier calving dynamics and ice habitat in fjords, further investigation is warranted in other regions where glaciers are in different states of advance and retreat to better understand the variability of ice habitat in fjords. Additionally, incorporating data related to variation in ice habitat into population models will be essential for understanding the influence of ice on seal distribution and abundance and for predicting future scenarios for seals that use tidewater glacier fjords.
Substantial seasonal and annual variability occurs in ice habitat and is driven by physical processes that influence calving dynamics of tidewater glaciers which can ultimately influence seal distribution, abundance, and behavior. However, our study demonstrates that harbor seals may respond to seasonal changes in ice habitat differently depending upon tradeoffs associated with life history events, such as pupping and molting, and energetic costs and constraints associated with the events. Thus, emphasizing the importance of considering the interaction between abiotic and biotic resources and life history events across the annual cycle in highly seasonal environments that are undergoing rapid environmental change (Varpe, 2017). Tidewater glaciers are naturally dynamic, and change can occur on daily, seasonal, annual, and multi-decadal time scales (Post et al., 2011). Furthermore, climate change models predict continued loss of land-based ice and changes to tidewater glaciers (Hock et al., 2019;Slater et al., 2019). Given the reliance of harbor seals on ice habitat for critical life-history events, harbor seals may serve as sentinels (e.g., Moore, 2008;Hazen et al., 2019) of ice in tidewater glacier fjords and provide insight into seasonal and annual changes associated with physical processes occurring along the ice-ocean boundary. Ultimately, elucidating biophysical linkages along the ice-ocean boundary can be facilitated through interdisciplinary collaborations using a systems-based approach (e.g., Catania et al., 2019;Straneo et al., 2019) to better understand how changes to tidewater glaciers will influence organisms that use these highly dynamic tidewater glacier fjord ecosystems.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.