Impacts of tidewater glacier advance on iceberg habitat

Abstract

Icebergs in proglacial fjords serve as pupping, resting and molting habitat for some of the largest seasonal aggregations of harbor seals (Phoca vitulina richardii) in Alaska. One of the largest aggregations in Southeast Alaska occurs in Johns Hopkins Inlet, Glacier Bay National Park, where up to 2000 seals use icebergs produced by Johns Hopkins Glacier. Like other advancing tidewater glaciers, the advance of Johns Hopkins Glacier over the past century has been facilitated by the growth and continual redistribution of a submarine end moraine, which has limited mass losses from iceberg calving and submarine melting and enabled glacier thickening by providing flow resistance. A 15-year record of aerial surveys reveals (i) a decline in iceberg concentrations concurrent with moraine growth and (ii) that the iceberg size distributions can be approximated as power law distributions, with relatively little variability and no clear trends in the power law exponent despite large changes in ice fluxes over seasonal and interannual timescales. Together, these observations suggest that sustained tidewater glacier advance should typically be associated with reductions in the number of large, habitable icebergs, which may have implications for harbor seals relying on iceberg habitat for critical life-history events.

1. Introduction

Tidewater glaciers are well known to respond nonlinearly to climate due to complex relationships between climate, ice flow and processes occurring at the glacier–ocean interface (i.e. iceberg calving, submarine melting and sediment deposition and erosion) (Post and others, 2011; Brinkerhoff and others, 2017; Robel and others, 2018). Consequently, tidewater glaciers can be out-of-phase with climate and neighboring glaciers (McNabb and Hock, 2014) and can go through cycles of slow advance and rapid retreat independent of climate change. Tidewater glacier advance is facilitated by the erosion and deposition of sediment at the glacier terminus, which provides flow resistance that allows sustained glacier thickening while also limiting mass losses from iceberg calving and submarine melting. Once climate forces a tidewater glacier to retreat from its end moraine, increases in ice velocities and calving rates produce an instability that is irreversible until the glacier retreats to a new pinning point (Pfeffer, 2007). The instability associated with calving retreat causes tidewater glaciers to retreat too far for the given climate; thus, once they reach a pinning point, tidewater glaciers begin to re-advance due to climate forcing. This process of slow advance and rapid retreat is commonly referred to as the tidewater glacier cycle.

The tidewater glacier cycle is associated with significant changes in fjord ecosystems. Freshwater fluxes to the ocean, in the form of both icebergs and subglacial discharge, vary during the tidewater glacier cycle due to feedbacks between climate, glacier geometry and glacier dynamics (e.g. Amundson and Carroll, 2018). Collectively, subglacial discharge and iceberg melting (an important source of freshwater in fjords; Enderlin and Hamilton, 2014; Moon and others, 2017) affect fjord circulation and stratification (e.g. Motyka and others, 2003; Davison and others, 2020; De Andrés and others, 2020), with implications for biogeochemical cycles (Hopwood and others, 2018; Kanna and others, 2022), nutrient distributions (Arimitsu and others, 2016), carbon sequestration (Smith and others, 2015; Hopwood and others, 2020), plankton populations (Arendt and others, 2016; Cuevas and others, 2019) and upper trophic levels (Lydersen and others, 2014; Urbanski and others, 2017; Womble and others, 2021).

Despite the clear links between glacier and ecosystem dynamics, few studies have attempted to link large-scale changes in glacier dynamics to ecosystem changes. Here, we leverage an extensive record of high-resolution aerial photographs of seals and icebergs from Johns Hopkins Inlet, Alaska (Womble and others, 2020) to quantify changes in iceberg habitat and harbor seal concentrations and to relate those changes to glacier dynamic processes.

1.1. Study area

Johns Hopkins Inlet (Tsalxaan Niyaadé Wool’éex'i Yé; 58$^\circ$49’ 31” N, 137$^\circ 8'$ 9” W; Fig. 1) is a tidewater glacier fjord located in the northwest corner of Glacier Bay (Sít’ Eetí Geeyi) National Park and Preserve. The fjord contains icebergs produced by Johns Hopkins Glacier (Tsalxaan Niyaadé Sít’) and adjacent Gilman Glacier, which are two of the few tidewater glaciers in Alaska that are currently advancing. As the glaciers have advanced, their termini have begun to coalesce into a single terminus. Since Gilman Glacier is just 10% of the area of Johns Hopkins Glacier (McNabb and Hock, 2014) and less by volume, its contribution to the total iceberg abundance is relatively small. We therefore focus our analysis on Johns Hopkins Glacier, which stretches ~21 km from its origins within the Fairweather Range to sea level and covers ~255 km² (McNabb and Hock, 2014). Johns Hopkins Glacier reached a minimum extent in the 1920s (Hall and others, 1995) following the disintegration of the Glacier Bay Icefield. Since 1948 the glacier has advanced over 1.6 km (McNabb and Hock, 2014) and thickened by over 100 m in its lower reaches (Larsen and others, 2007). From 1972 to 2009, the fjord filled with sediment at rates of ~0.5 m a⁻¹ near the mouth of the fjord, 1–2 m a⁻¹ in the near-terminus region, and >2 m a⁻¹ within 0.5 km of the glacier terminus (Hodson and others, 2013).

Figure 1. Map of the study site showing (a) Johns Hopkins Inlet and Glacier Bay and their location within (b) Alaska and (c) Glacier Bay National Park and Preserve. In (a), small white boxes indicate image footprints from an aerial survey flown on 9 June 2019 and are representative of imagery obtained during all surveys, the purple line indicates the centerline profile, and blue, orange and green points indicate points 1.5, 3.5 and 5.5 km from the 2021 terminus position (used when plotting velocities and elevations in Fig. 4). (d) Close up of terminus region outlined in the dashed red box in (a). Colored profiles indicate the terminus positions from 1935 to 2021. The background image in (a) and (d) is a Sentinel-2 image from 2018.

Johns Hopkins Inlet hosts the largest aggregation of harbor seals (Phoca vitulina richardii) in Glacier Bay with up to 2 000 seals aggregating seasonally (Calambokidis and others, 1987; Mathews and Pendleton, 2005; Womble and others, 2010, 2020). After extensive movements during the post-breeding season, seals begin to arrive in Johns Hopkins Inlet in late April to mid-May for pupping (Womble and Gende, 2013). The number of seals in Johns Hopkins Inlet peaks during the pupping period in June and the molting period in August (Womble and others, 2021). Icebergs are not subject to tidal inundation and likely provide refuge from predation for young pups and stable platforms for nursing young (Womble and others, 2014). Glacial sites, such as John Hopkins Inlet, may serve as source populations for surrounding regions (Womble and others, 2010), and thus changes in glacier dynamics and associated iceberg habitat may have significant impacts on harbor seals regionally in tidewater glacier fjords.

2. Data and Methods

We build on previous work to quantify the impact of glacier dynamics on iceberg habitat (McNabb and others, 2016). We focus on the time period from 2007 to present, coinciding with the onset of aerial photographic surveys in the fjord that have been used to estimate the abundance and distribution of seals (Womble and others, 2020, 2021). Figure 2 provides a timeline of the datasets that we use in this study.

Figure 2. Timeline of data sources used in this study.

2.1. Terminus position and ice velocity

Glacier terminus positions from 2007 to 2012 were obtained from McNabb and Hock (2014) and positions from 2012 to 2021 were manually delineated from Landsat and Sentinel optical imagery. Glacier length was calculated using the ‘box method’ described in Moon and Joughin (2008).

Velocity data were generated using auto-RIFT (Gardner and others, 2018) and provided by the NASA MEaSUREs ITS_LIVE project (Gardner and others, 2022). Velocity measurements were extracted at three points along the glacier centerline (taken from Kienholz amd others, 2014) at distances of 1.5, 3.5 and 5.5 km from the glacier's 2021 terminus position (Fig. 1a). The ITS_LIVE compilation of mean surface velocities for 2007–2018 were subject to limited image availability with scene pairs ranging from 6 to 546 days. Data quality improved drastically in the later part of our survey when Landsat 8, Sentinel-1 and Sentinel-2 data became incorporated in 2014, 2015 and 2017, respectively, allowing for the analysis of seasonal velocity variations.

2.2. Digital elevation models

We obtained digital elevation models (DEMs) from the Shuttle Radar Topography Mission (SRTM) carried out from 11 to 22 February 2000 (Farr and others, 2007), Interferometric Synthetic Aperture Radar (IfSAR) data collected in summer 2010 (https://doi.org/10.5066/P9C064CO) and a combination of Worldview 1, 2 and 3 satellite stereo pairs available through the Polar Geospatial Center's ArcticDEM portal (Porter and others, 2018). All DEMs were co-registered to the IfSAR DEM with the demcoreg Python module (Shean and others, 2016), using the method of Nuth and Kääb (2011). Areas covered by glacier ice or having a slope <0.1$^\circ$ or >40$^\circ$ were excluded during the co-registration. We take the normalized median absolute deviation (NMAD) to represent the uncertainty in the non-glacierized area within the DEMs. NMAD errors ranged between 2.6 and 5.4 m, with some of the error associated with snow cover. We then used the centerline points 3.5 and 5.5 km upstream from the terminus location in 2021 as representative locations to discuss terminus elevation changes below; we excluded the 1.5 km point from our analysis due to noise associated with its proximity to the glacier front. Additionally, we found an ICESAT2 trackline for ATL06 data from 24 May 2021 that overlapped with the point 3.5 km from the terminus (Smith and others, 2020, accessed through OpenAltimetry (Khalsa and others, 2022)).

2.3. Fjord bathymetry

The bathymetry of Johns Hopkins Inlet was surveyed in 1972, 2009 and 2020 by the National Oceanic and Atmospheric Administration (NOAA; see https://coastalscience.noaa.gov/products/noaa-bathymetric-data-viewer/). The 1972 data consist of point data spanning the full length of the fjord. While the 2009 data also covered the full length of the fjord, the 2020 data were limited to an ~2.6 km² region spanning ~350 m to 2.5 km in front of the glacier terminus because dense ice coverage prevented a more comprehensive survey. The 2009 and 2020 datasets were originally made available as Bathymetric Attributed Grids (BAG) files with variable horizontal resolution ranging from 1 to 16 m; we resampled the data to 10 m and converted to geotiff format. Hodson and others (2013) compared the 1972 and 2009 data and computed sedimentation rates. We expand on their analysis to include the 2020 data and provide a focused assessment on the morphology of Johns Hopkins Inlet.

2.4. Aerial photographic surveys

Aerial photographic surveys (n=91) were conducted in Johns Hopkins Inlet from 2007 to 2019 and targeted the harbor seal pupping (June) and molting (August) periods (Womble and others, 2020, 2021). Due to georectification issues we only use 89 of the surveys here. Surveys were conducted from a de Havilland Canada DHC-2 Beaver single-engine high-winged aircraft (Ward Air Inc., Juneau, Alaska). The surveys were flown along 12 established transects spaced 200 m apart and provided systematic sampling of the entire fjord (Fig. 1a).

Photos were taken directly under the plane at 2 s intervals using a 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 photo sample rate and transect spacing prevented the overlapping of images and ensured seals and icebergs were not double-counted. An onboard global positioning system (Garmin 76 CSX) recorded the position of the plane along the transects at 2 s intervals for the later pairing of latitude, longitude and altitude with each photo. Each image covered about 80 m × 120 m on the ground. A trained observer reviewed each digital image and counted seals located on icebergs and the total number of seals on icebergs in each digital image was summed for each survey day. Here we define the seal concentration as the total number of seals (both non-pups and pups) on icebergs divided by the total survey area for each survey day. Further details on the survey methods are provided in Womble and others (2020).

2.5. Iceberg segmentation

Iceberg segmentation was conducted using the openCV Python module (https://opencv.org), building on previous methods developed by McNabb and others (2016). In order to avoid detecting glacier ice or other non-fjord features, we only used photos that were entirely contained within the fjord (i.e. photos containing the fjord sidewalls or glacier terminus were removed). Our approach is similar to that of McNabb and others (2016). The photos were first converted from RGB to HSV before applying a segmentation step based on pixel brightness values. We then applied an 11 × 11 (40 cm × 40 cm) Gaussian blur to remove small, noisy features such as whitecaps. Next, we enhanced the photo contrast by stretching the images from 0 to 255. A threshold of 170 was then applied where we considered anything above this value to be ice. Iceberg edges were detected using a Canny edge detector with a 3 × 3 (11 cm × 11 cm) kernel before the edges were dilated. Finally, we found all closed regions within the image, excluding interior contours, and calculated the area within the contours.

Figure 3 provides an example of the iceberg segmentation method. Our segmentation method worked well for most images. However, some issues arose when icebergs were closely packed, causing small icebergs or brash ice to be clumped into larger icebergs; this is a well-known issue with iceberg segmentation algorithms (e.g. Rezvanbehbahani and others, 2020). While this does not impact the total ice coverage calculations, it does affect the iceberg size distributions; however, we suggest that the impact on the shape of the distributions is relatively small.

Figure 3. Aerial photos of the fjord overlain with results from the iceberg segmentation method for (a) low and (b) high ice concentrations.

3. Results

3.1. Glacier velocities

Annual (2007–2014) and seasonal speeds (2014–2022) of Johns Hopkins Glacier are shown in Figure 4a. Overall, the velocity record indicates a gradual slowdown, especially from 2013 to 2021, during which time the velocity decreased by ~45%.

Figure 4. Time series of (a) glacier velocity at the points labeled in Figure 1a (stair plots are from ITS_LIVE annual velocities and point velocities are from ITS_LIVE-Scene-pairs Version 2), (b) glacier length relative to the confluence of the tributary glaciers, (c) change in elevation at points 3.5 and 5.5 km relative to 2000 (IfSAR and ICESat-2 data are denoted by the triangle and star, respectively). Error bars are normalized median absolute deviation (NMAD) values, (d) ice fraction and (e) seal concentration.

Strong seasonal velocity variations are evident during the time period of 2014–2022 (i.e. when the available velocity data can resolve seasonal velocities) on the order of ~1 km a⁻¹, with peak velocities a factor of 4–5 times higher than minimum velocities. Peak velocities occurred in mid-May and minimum velocities were in early October. This pattern likely reflects a strong runoff (rain and meltwater) influence where a large influx of water to the base of the glacier increases basal motion in spring and early summer, when the subglacial drainage system is poorly developed. Velocities then decrease throughout the summer as the drainage system becomes more efficient.

3.2. Terminus advance and glacier thickening

Between 2007 and 2021, the terminus of Johns Hopkins Glacier advanced ~230 m and the lower reaches of the glacier thickened by ~30 m. However, the rates of advance and thickening were not steady and were marked by varying amounts of seasonal retreat and surface lowering.

The glacier steadily advanced and thickened during the 2007–2013 time period, with ~190 m of advance and ~12 m of thickening (the elevation in 2007 was estimated using a linear interpolation between the 2000 STRM and 2010 IfSAR DEMs). The rate of terminus advance is fairly stable during this time period and seasonal retreat is limited to a maximum of ~80 m. In contrast, the 2013–2019 period is marked by more drastic seasonal cycles. The most significant seasonal retreat occurred in summer 2016 (~330 m of retreat occurred between May and September of 2016). We also found a surface lowering of ~22 m between 2014 and 2017, coincident with the retreat that occurred at that time. The terminus position appears to have (re-)stabilized during the 2019–2021 period; seasonal advance/retreat is now <~60 m, the terminus is slowly advancing, and the lower glacier has thickened by 20 m.

3.3. Changes in bathymetry and moraine growth

The bathymetry of Johns Hopkins Inlet reveals a fjord with steep side walls and a relatively flat floor (Fig. 5a). A large amount of glaciomarine sediment has accumulated in Johns Hopkins Inlet since the last deglaciation and during its current advance. Seismic campaigns indicate that the average sediment accumulation rate throughout the fjord was about 1.8 − 2.0 × 10⁷ m³a⁻¹ from 1892, when the glacier completely filled the fjord and was retreating, to 1979 (Cai and others, 1997).

Figure 5. Comparison of bathymetric surveys. (a) Fjord bathymetry in 2009. White line illustrates the centerline track from 1972 used for cross-sectional analysis in (c). (b) Sedimentation rate between 2009 and 2020. (c) Cross-section of bathymetry. Vertical lines show the position of the glacier terminus. Colors correspond to the colorbar in Figure 1c.

NOAA bathymetry surveys show that from 1972 to 2009 the fjord filled with sediment at rates of >2 m a⁻¹ in the region ~0.5 km from the glacier terminus, 1–2 m a⁻¹ in the near-terminus region and 0.5 m a⁻¹ near the mouth of the fjord (Hodson and others, 2013). Analysis of the 2009–2020 time period found a similar sedimentation rate of ~1.0 m a⁻¹ in the near-terminus region, with a total of ~10 m of sediment added over the 11-year time period (Fig. 5b). Conversely, the data suggest a substantial loss of sediment along the fjord walls of ~5–0 m a⁻¹, although this could also be due to data issues in areas of steep terrain (as discussed by Hodson and others, 2013). The end moraine continues to thicken and shoal, and first surfaced at low tide during July 2019 (Fig. 6). Analysis of Worldview imagery from 2014–2019 with acquisition timestamps close to low tide did not reveal surfacing prior to 2019. Following the moraine surfacing, the terminus position stabilized and exhibited less seasonal retreat and the glacier thickened ~15 m from 2019 to 2021.

Figure 6. Photos documenting the surfacing of the moraine in summer 2019. Figures (a)–(c) were taken during aerial surveys and (d) was taken from a kayak.

3.4. Variations in iceberg habitat

Concurrent with a reduction in glacier velocities and the shoaling surfacing of the moraine, we observe a reduction in ice fraction (iceberg area divided by the total surveyed area) and harbor seal concentration (number of harbor seals divided by the total surveyed area) (Figs 4d,e). Seal concentration appears to depend logarithmically on the ice fraction (Fig. 7), suggesting that changes in ice availability influence the number of seals in the fjord, especially when ice coverage is low.

Figure 7. Harbor seal concentration versus ice fraction for each aerial survey. Colors indicate the pupping (June) and molting (August) seasons.

Despite these changes in ice coverage, we observe relatively little variability in the iceberg size distributions. The iceberg size distributions are heavy-tailed and plot as approximately straight lines in log-log space (Fig. 8), suggesting that they are power law distributions. Some of the distributions exhibit a kink at around 50 m², which is due to the iceberg segmentation process lumping small icebergs into one. These erroneously large icebergs are very low probability and typically have little impact on the rest of the distribution (except for surveys with high ice fractions). We therefore exclude icebergs >50 m² in our analysis.

Figure 8. Empirical complementary cumulative distribution function across all aerial surveys.

In addition to setting a maximum iceberg size of a _max, power law probability density functions require a minimum iceberg size of a _min, which we set equal to 1 m². The probability density function is given by

(1)$$p( a) = \left({1-\alpha\over a_{\rm max}^{1-\alpha}-a_{\rm min}^{1-\alpha}}\right)a^{-\alpha},\; $$

where α is a constant. The cumulative distribution function (CDF), which indicates the probability that a randomly selected iceberg has an area less than or equal to a, is found by integrating Eqn (1), yielding

(2)$$P( A\leq a) = {a^{1-\alpha} - a_{\rm min}^{1-\alpha}\over a_{\rm max}^{1-\alpha}-a_{\rm min}^{1-\alpha}}.$$

We estimate α for each survey by using the maximum likelihood method in the Python powerlaw module (Alstott and others, 2014) to optimize the fit between Eqn (2) and our empirical cumulative distribution function. The complementary cumulative distribution function (CCDF), which indicates the probability that a randomly selected iceberg has an area greater than a, is given by

(3)$$P( A> a) = 1 - P( A\leq a) = {a_{\rm min}^{1-\alpha} - a^{1-\alpha}\over a_{\rm max}^{1-\alpha}-a_{\rm min}^{1-\alpha}}.$$

Note that the CCDF does not produce a straight line in log-log space because we have selected an upper bound of a _max. We find that power law distributions provide a good fit to the data (Fig. 9), with a D statistic (maximum difference between the empirical and theoretical CDFs) typically around 0.01.

Figure 9. Power law fit to iceberg size distributions. (a)–(c) Example of the best-fit power law distribution for a survey on 14 August 2013. The best-fit power law exponent for each survey is shown vs. (d) time and (e) ice fraction. The solid line indicates the mean value and the dashed lines indicate the standard deviation from the mean.

Across 89 surveys, we observe a mean exponent of $\overline \alpha = 2.33$ and a standard deviation of $\sigma _\alpha = 0.093$ (Fig. 9). These values are in agreement with those of Neuhaus and others (2019), who report values of α = 2.08–2.35 for Columbia Glacier, Alaska, which also produces small icebergs compared to those found in Greenland and Antarctica. Glaciers that produce large, tabular icebergs are expected to yield iceberg size distributions with lower values of α (Åström and others, 2021). We observe no clear long-term trends in the power law exponent. There is a tendency for the exponent to be lower when the ice coverage is larger (typically in early summer) but this is largely a consequence of the iceberg segmentation process performing poorly when icebergs and brash ice are closely packed, which we verified by visually inspecting segmented images.

Theoretical work has suggested that iceberg calving and fragmentation processes should produce power law distributions with exponents similar to what we observe (Åström and others, 2021). We cannot exclude other heavy-tailed distributions, such as lognormal distributions. Observations of icebergs that have experienced significant drift suggest that the distributions tend toward lognormal distributions after experiencing significant melt (Kirkham and others, 2017), which may explain the slight curvature seen in some of the empirical CCDFs (Fig. 8). Nonetheless, the observed distributions can be reasonably approximated as power law distributions.

Variations in α can be attributed to variations in melting, which tends to reduce the exponent because large icebergs lose area at a faster rate than small icebergs, or by variations in calving and iceberg fragmentation processes. Since we are unable to detect any clear trends in the exponent that might elucidate how changes in glacier or fjord dynamics influence iceberg size distributions we therefore take $\sigma _\alpha$ to represent the uncertainty in the exponent (and not as a reflection of changes in environmental conditions).

4. Discussion

4.1. Impact of moraine growth on glacier dynamics, ice coverage and ice habitat for seals

Johns Hopkins Glacier has been advancing since the 1920s. Our observations indicate that the glacier has continued to advance and thicken during the early part of the 21st century, although its rate of advance has declined over the past three decades (McNabb and Hock, 2014). The advance has been associated with the progradation and thickening of an end moraine.

Some tidewater glaciers experience seasonal variations in velocity, especially in their lower reaches, due to changes in stress associated with seasonal terminus advance and retreat (e.g. van der Veen, 2002), while others appear to respond strongest to meltwater input (Moon and others, 2014). Johns Hopkins Glacier falls into the latter category. Seasonal variations in velocity are large, up to about a factor of five, but (i) peak velocities occur in early summer when the terminus is in an advanced position and (ii) seasonal variations in terminus position are generally small (<50 m). Larger variations in terminus position did occur from about 2016–2019, concomitant with reductions in surface elevation. However, these changes were not sufficiently large to cause a significant glacier dynamic response upstream. We suggest that the larger seasonal variations in length occurred because the glacier began to retreat off of the moraine, but the retreat stalled when changes in glacier mass balance or sedimentation at the terminus allowed the glacier to regain footing. From 2019 onward, the advance has been steady with very little seasonality in terminus position.

Over seasonal timescales, ice coverage and seal concentrations are typically highest in June, when glacier velocities are also high (Fig. 4). However, Womble and others (2021) found a stronger correlation between seal numbers and ice coverage during the pupping season in June than molting season in August. This pattern suggests that seals may respond to changes in ice habitat differently depending upon life-history events and energetic constraints imposed by a dependent pup in June but not in August, once pups have been weaned.

Over longer timescales, we observe that ice velocities, ice coverage and seal concentration have decreased steadily since 2014. The continued growth of the end moraine (Fig. 5), and its emergence above sea level in 2019 (Fig. 6) suggests a direct link between glacier dynamics and iceberg habitat for seals. Growth of the moraine has caused Johns Hopkins Glacier to behave dynamically like a land-terminating glacier, with strong seasonality in ice flow linked to the seasonal evolution of the subglacial drainage system (e.g. Nienow and others, 2017). As the moraine grew, it provided additional flow resistance, causing the glacier to slow and thicken (see e.g. Amundson, 2016) and it reduced the surface area available for iceberg calving. Consequently, this has led to a reduction in icebergs and seals in the fjord in both June and August since the beginning of our observations in 2007. However, this reduction is most prominent in June when seals are likely more dependent upon iceberg habitat as it provides a refuge from predation for young pups and also provides a stable platform for nursing. We found a positive relationship between iceberg and seal concentrations, especially with low ice coverage, suggesting that seals may adjust their use of Johns Hopkins Inlet when shifts in ice availability occur (Fig. 7). When ice availability is reduced, seals may move to other sites (ice and/or terrestrial) or spend more time in the water which may result in behavioral changes and ultimately have fitness-level implications.

4.2. Implications for iceberg habitat

Our observations indicate relatively little variability in iceberg size distributions at Johns Hopkins Inlet during a 15-year period, despite significant variations in glacier velocities and fluxes. This suggests that iceberg size distributions are an intrinsic property of calving and fragmentation processes, as supported by the modeling work of Åström and others (2021). Furthermore, Womble and others (2021) did not find a relationship between the number of seals and iceberg size during pupping or molting periods from 2007 to 2014. Thus, as long as Johns Hopkins Glacier continues to calve small (i.e. non-tabular) icebergs, then the form of the iceberg size distributions that we observe should hold and we can use them to make predictions about the number of habitable icebergs in a fjord at any given point in time.

The number of habitable icebergs, n, in a fjord is

(4)$$n = NP( A> a_{\rm h}) ,\; $$

where N is the total number of icebergs in a fjord and P(A > a _h) is the probability of selecting an iceberg larger than the minimum habitable iceberg size a _h. McNabb and others (2016) suggest that harbor seals require a minimum iceberg size of 1.6 m².

We previously set a maximum iceberg size of a _max due to issues with the iceberg segmentation procedure that sometimes resulted in small icebergs being lumped into much larger icebergs. For the purposes of this discussion, we assume that iceberg size distributions can be reasonably approximated as power law distributions without a maximum size threshold. Thus, the probability density function and complementary cumulative distribution functions are given by

(5)$$p( a) = {\alpha -1\over a_{\rm min}}\left({a\over a_{\rm min}}\right)^{-\alpha}$$

and

(6)$$P( A> a) = \left({a\over a_{\rm min}}\right)^{1-\alpha},\; $$

respectively.

The number of icebergs in the fjord is related to the ice coverage by

(7)$$A_{\rm total} = \int_{a_{\rm min}}^\infty Np( a) a\, {\rm d}a,\; $$

where A _total represents the total ice coverage. Substituting in the probability density function (Eqn (5)) and evaluating yields

(8)$$A = N\int_{a_{\rm min}}^\infty ( \alpha-1) \left({a\over a_{\rm min}}\right)^{1-\alpha}\, {\rm d}a = Na_{\rm min}\left({\alpha-1\over \alpha-2}\right),\; $$

where α is required to be >2 (as is the case with our data).

Combining Eqns (4), (6) and (8), we find that the number of habitable icebergs for a given ice coverage is

(9)$$n = \left({\alpha-2\over \alpha-1}\right){A_{\rm total}\over a_{\rm min}}\left({a_{\rm h}\over a_{\rm min}}\right)^{1-\alpha}.$$

We estimate the uncertainty in n by using the standard propagation of uncertainty and assuming that the only uncertainty is from the exponent α, resulting in

(10)$$\sigma_n = \sqrt{\left({{\rm d}n\over {\rm d}\alpha}\right)^2\sigma_\alpha^2}.$$

Inserting Eqn (9), this becomes

(11)$$\sigma_n = n\left[{1\over ( \alpha-2) ( \alpha-1) } - \ln{\left({a_{\rm h}\over a_{\rm min}}\right)}\right]\sigma_\alpha.$$

Using values of α = 2.33, $\sigma _\alpha = 0.093$, a _min = 1 m² and a _h = 1.6 m² yields σ _n = 0.17n. Thus, for a given ice coverage, the number of habitable icebergs will be within 17% of the number predicted by Eqn (9). Since the number of habitable icebergs is approximately proportional to the ice coverage, and the ice coverage scales with ice flux (see Figs 4a,d), we can expect the number of habitable icebergs to continue to decrease with decreasing ice flux if the proglacial moraine continues to grow and impede ice flow and calving. This would not necessarily be the case if the iceberg size distributions exhibited more stochastic or secular variability than what we have observed. Many factors beyond iceberg size determine whether seals will utilize particular icebergs, such as timing of critical life events as well as seal age, sex and social behavior; future studies should attempt to link physical processes to iceberg utilization by harbor seals.

5. Conclusions

We used multiple data sources to characterize the evolution of the glacier–fjord environment at Johns Hopkins Glacier and Inlet over the past two decades. Satellite-derived velocity measurements and digital elevation models indicate the glacier has continued to advance and thicken and has slowed down in recent years. From 2007 to 2021, the glacier advanced ~160 m, thickened ~23 m and exhibited a gradual slowdown, especially from 2013 to 2021, during which time the velocity decreased by ${\sim }45\%$. Analysis of aerial photographs indicates that concurrent with the slowdown was a decrease in ice coverage and harbor seal concentrations on icebergs during the pupping period in June and the molting period in August. Moreover, satellite and aerial observations indicate the surfacing of an end moraine in 2019. We see the influence of this growing moraine over the past two decades throughout many of our datasets, such as a reduction in seasonal retreat following its surfacing and an overall decrease in iceberg discharge as the terminus became increasingly grounded. In addition, despite large changes in ice fluxes over the past two decades, we find little variability in iceberg size distributions at Johns Hopkins Inlet, thus implying that the number of habitable icebergs is proportional to overall ice coverage and establishing a direct link between glacier dynamics and seal habitat. We expect a similar pattern at other advancing tidewater glaciers where sustained terminus advance eventually leads to the formation of a large end moraine, a decrease in calving rates and a reduction in habitable icebergs for seals.