Glacial fronts are important summer habitat for narwhals ( Monodon monoceros ); however, no studies have quantified which glacial properties attract whales. We investigated the importance of glacial habitats using telemetry data from n = 15 whales tagged in September of 1993, 1994, 2006 and 2007 in Melville Bay, West Greenland. For 41 marine-terminating glaciers, we estimated (i) narwhal presence/absence, (ii) number of 24 h periods spent at glaciers and (iii) the fraction of narwhals that visited each glacier (at 5, 7 and 10 km) in autumn. We also compiled data on glacier width, ice thickness, ice velocity, front advance/retreat, area and extent of iceberg discharge, bathymetry, subglacial freshwater run-off and sediment flux. Narwhal use of glacial habitats expanded in the 2000s probably due to reduced summer fast ice and later autumn freeze-up. Using a generalized multivariate framework, glacier ice front thickness (vertical height in the water column) was a significant covariate in all models. A negative relationship with glacier velocity was included in several models and glacier front width was a significant predictor in the 2000s. Results suggest narwhals prefer glaciers with potential for higher ambient freshwater melt over glaciers with silt-laden discharge. This may represent a preference for summer freshwater habitat, similar to other Arctic monodontids.

1. Introduction

Arctic and subarctic glacial fjords are characterized by high rates of productivity that lead to rich marine ecosystems, including high densities of seabirds, marine mammals and fishes [1]. High productivity in Greenland's glacial fjords and their downstream regions has been attributed to glacial meltwater, with a strong correlation between the presence of meltwater nutrients and phytoplankton blooms [2]. These plumes may aggregate plankton or stun plankton via freshwater osmotic shock [3], making them easy prey for larger surface-feeding predators and multiple trophic levels. Nutrient fluxes at the glacier fronts are also used for post-bloom plankton production, lengthening overall feeding opportunities in summer. In some areas of the Arctic where the permanent multi-year sea ice has vanished, glacial fjords are replacing sea ice habitat for ice-breeding species [3].

The West Greenland narwhal (Monondon monoceros) subpopulation, with a mean subpopulation abundance estimated at approximately 6000 animals in 2007 [4], occurs in Melville Bay and frequents glacial fronts in summer and autumn [5,6]. It is unknown why narwhals have an affinity for glaciers; physical properties of fjords may offer enhanced feeding opportunities, though to date there has been little evidence of summertime feeding. Narwhals may also be attracted to subsurface freshwater melt at the glacier face which may resemble estuarine habitat used by other Arctic odontocetes, e.g. beluga (Delphinapterus leucas) [7].

Using satellite remote sensing data collected over two decades, we examined the suite of glaciers visited by narwhals in Melville Bay. We developed quantitative covariates to describe individual glaciers and used proximity analyses in statistical models to examine relationships between narwhal occurrence and glacier fjord covariates. We examined narwhal use of the ‘near-glacier’ region (within approx. 7 km of the ice front); however, we refer to this generally as the ‘glacier front’ for lack of pre-existing terminology. This study sheds light on what glacier features may be selected by narwhals and improves our understanding of how future changes in freshwater melt [8] may influence narwhal habitat.

2. Methods

(a) Narwhal data

Narwhals were captured and instrumented with satellite-linked time–depth–temperature recorders in September of 1993, 1994, 2006 and 2007 in Melville Bay, West Greenland [6,9–11]. We included locations with argos classes of less than or equal to 1.5 km accuracy and positions between September and November, including the start of the southbound migration [5]. Locations were removed using speed (greater than or equal to 1.8 m s−1) and angular (default) filters in R v. 2.13.2 [12] using the package ‘argosfilter’ [13]. Resulting whale locations were reduced to a single position per whale per day during peak of satellite passage to decrease autocorrelation bias, standardize temporal sampling and address the effects of different duty cycles. We used a correlated random walk model to estimate locations based on observed filtered locations and associated argos error (‘crawl’ package [14]). The result was a dataset of predicted and observed locations across all days of the study containing 815 predicted and 763 observed locations (figure 1). We also confirmed narwhal presence at each glacier by manually checking each individual's locations. Figure 1. Map of study area in Greenland showing observed narwhal and glacier locations. Inset: glacier front detail north of Kullorsuaq (red lines are front locations) with yellow lines representing 7 km buffer zones. Imagery from https://nsidc.org/data/NSIDC-0633.

(b) Glacier and fjord data

We focused on 41 northwest Greenland glaciers (figure 1), covering the coastal region of narwhal activity. Table 1 provides information on the suite of glacier variables. For each glacier, we created a single ‘glacier point’ by computing the terminus line centroids and spatially joined narwhal presence or visits to each glacier point and glacier covariates ([15,16], table 1, figure 2). Bathymetry was from Oceans Melting Greenland (http://dx.doi.org/10.5067/OMGEV-BTYSS) and gravity data [20]. Figure 2. (a) Glacier schematic with covariate labels, adapted from [23], (b) modelled relationships (and 95% CI) among covariates from top models in the 2000s at 7 km, explaining the number of visits, proportion of narwhals and probability of visits by narwhals at glaciers.

Table 1.List of glacier covariates. Collapse covariate data time period description use of covariates reference glacier front width (m)a autumn to winter of: 1992 and 2000 linear distance (m) across glacier terminus ice front. Owing to change in ice front positions, terminus width can be difficult to define. We use a conservative measurement that provides minimum width value and captures relative width across all study glaciers — 1992 widths used for 1990s narwhal data — 2000 widths used for 2000s and combined 1990s/2000s narwhal data [15] glacier front change (km)a 1992–2005, 1992–2000, 2000–2005 change in glacier ice front position (either advance or retreat) during the specified time period — 1992–2000 used for 1990s narwhal data — 2000–2005 used for 2000s narwhal data — 1992–2005 used for combined 1990s and 2000s narwhal data [15] glacier ice velocity (m yr−1)a autumn to winter of: 2000 and 2005 surface ice velocity measured approximately 1/2 glacier width up the glacier from the ice front. Measurement made at roughly glacier centreline (area of fastest flow). Individual glacier velocities are highly correlated (r = 0.9 or higher) for each glacier and relative glacier velocities across all glaciers are stable — 2000 velocities used for combined 1990s and 2000s narwhal data — 2005 velocities used with 2000s narwhal data [16] glacier ice thickness (m)a 2006–present glacier ice depth from Operation IceBridge MCoRDS L3 gridded ice thickness data. Value is 2 km-zonal average calculated from the glacier point (see text for detail on glacier point) single thickness measurement used with all narwhal data [17] iceberg discharge: sea surface coverage (classification based on % coverage)a Sep 1–Oct 31 of years 2006 and 2007 estimated per cent (categorized in: 0–25%, 25–50%, 50–75% and 75–100%) ice coverage within a 3 km buffer proximity zone of the glacier front. Based on 61 RGB Moderate Resolution Imaging Spectroradiometer (MODIS) images (250 m resolution) during the narwhal visitation period georeferenced with a common set of control points and a first-degree polynomial transformation using ArcMap (ESRI) single discharge average value used for all data produced for this study iceberg discharge: spatial extent (km) Sep 1–Oct 31 of years 2006 and 2007 classified by per cent coverage in a 3 km zone at the glacier. If ice discharge sea surface coverage was more than 50%, iceberg discharge extent was measured. Measurements of ice discharge distance started at the midpoint of the glacier terminus line and ended as far away as the discharge extent could be reasonably attributed to the individual glacier, with median distance selected for multiple possible endpoints. Discharge extent was assigned to all glaciers sharing the same discharge region single discharge average value used for all data produced for this study subglacial water discharge (km3 yr−1)a yearly mean calculated from model data spanning 1999–2013 discharge estimates are derived from the Racmo2.3 surface mass balance climate model, partitioned into individual meltwater outlets using a D8 flow routing algorithm [18,19] subglacial sediment flux (Mt yr−1) derived from several datasets covering 1999–2013 sediment flux is estimated from an algorithm that uses ice velocity and ice thickness data to estimate the sediment concentration of the water leaving each hydrologic catchment. This concentration is multiplied by the yearly discharge to get sediment flux not used in models because correlated to water discharge [19] bathymetry: max. depth (m)a OMG from 2015 and gravity data 2010–2012 classified as the maximum depth in a 3 km buffer offshore from each glacier front. OMG data were used for all available glaciers. Bathymetry data from Operation IceBridge L4 inversions of gravity anomalies were used where OMG data were not available. Survey data have a line spacing of 5 km and along-track resolution of 4 km. Average bathymetry values were computed from both OMG and gravity anomalies using a 3 km zonal average at each glacier point OMG bathymetric measurements were used where available, with bathymetry from gravity anomalies used for the remaining glaciers [20–22]

(c) Narwhals and glaciers proximity analysis

We created three proximity buffers around each glacier centroid point (at 5, 7 and 10 km to serve as a sensitivity analysis) and quantified whale visits within these regions on decadal (1990s, 2000s) and combined (1990s plus 2000s) time periods. We estimated (i) whale presence or absence, (ii) the total number of visits by whales in a 24 h period, and (iii) the fraction of whales that visited a glacier in each decade.

(d) Statistical methods

We assessed collinearity among predictors by calculating Pearson's correlation coefficients, resulting in a reduced set of variables (with pairwise correlations ≤0.6) to estimate the relationship between narwhals and glacial predictors. We used generalized linear models (GLMs) to identify physical predictors associated with narwhal attendance for each proximity buffer and time period. The total number of visits by whales was modelled with a Poisson's error structure, while the fraction of whales visiting glaciers and probability that a whale visited a glacier were modelled as binomial GLMs. We used stepwise model selection based on the lowest Akaike's Information Criteria value [24].

3. Results and discussion

We tracked the autumn movements of 15 adult narwhals over 4 years in Melville Bay (1993 and 1994: n = 8, 3M : 5F; 2006 and 2007: n = 7, 2M : 5F). Tracking durations between September and November ranged from 21 to 88 days in the 1990s and 48 to 90 days in the 2000s. Across both decades there were two clusters of common glaciers that were visited by whales: north of Cape Seddon and the Fisher Islands (figure 1). Additionally, in the 2000s, whales visited glaciers north of Kullorsuaq. At the largest distances (10 km), whales visited twice as many unique glaciers in the 2000s as the 1990s (21 glaciers and 10 glaciers, respectively), yet differences between decades declined with declining distance radii. In the 2000s, whales visited a larger numbers of glaciers owing to the loss of autumn fast ice and increased availability of habitat as demonstrated by an expanded range along the coast (figure 1). Changes in the timing of sea ice advance and retreat have been profound in Baffin Bay and Melville Bay [1]. Sea ice freezes up 3.5 weeks later than in 1979 [25] and fast ice at the glacier fronts in summer is now rarely present.

It is unknown at what distances glacier fjords attract narwhals; thus sensitivity analysis examined predictors at multiple scales. The set of significant predictor variables was consistent across all scales (5, 7 and 10 km) and for the three visitation metrics (table 2). Sensitivity analyses were important because observation and modelling studies at Greenland outlet glaciers have demonstrated notable spatial differences in fjord water properties across scales used in this study (5–10 km), including variations in salinity, temperature and sediment from subglacial water plumes [26,27]. We did not include subsistence hunting pressure in our analyses because it was difficult to quantify. Glaciers close to Savissivik and Kullorsuaq have higher hunting pressure than glaciers inside Melville Bay, where no hunting is supposed to occur because the area is protected. There may be an avoidance response around these communities owing to hunting pressure regardless of glacial features and this may impact habitat selection. Finally, some models at 5 km in the 1990s did not converge owing to low sample sizes.

Table 2.Final GLM results for response metrics in all years pooled at three radii from glacial fronts. Bold numbers indicate p < 0.05. GLMs at 10 km were estimated with quasi-likelihood model structures to account for overdispersion for sum total visits (all time periods) and proportion of narwhals visiting glaciers. Dashes indicate covariate was not included in the final model; blanks indicate models fail to converge at the smaller radii. Collapse narwhal response covariate 10 km radius 7 km radius 5 km radius estimate s.e. p-value estimate s.e. p-value estimate s.e. p-value sum total visits glacier ice velocity −0.001 0.000 0.016 −0.001 0.000 0.005 −0.001 0.001 0.036 glacier ice thickness 0.013 0.002 <0.001 0.015 0.003 <0.001 0.015 0.005 0.001 proportion visiting glaciers glacier ice velocity — — −0.001 0.000 0.042 −0.001 0.000 0.126 front width (2000s) 0.001 0.000 0.030 — — — — — — glacier ice thickness 0.008 0.002 <0.001 0.014 0.003 <0.001 0.012 0.004 0.003 probability of visiting a glacier glacier ice velocity — — — −0.003 0.002 0.067 −0.002 0.001 0.088 glacier ice thickness 0.025 0.009 0.004 0.084 0.031 0.008 0.022 0.009 0.013 subglacial water discharge −1.780 0.967 0.065 −6.966 4.247 0.101 — — — iceberg discharge extent — — — 2.611 1.352 0.053 — — — bathymetry: max. depth — — — −0.012 0.007 0.095 — — —

Ice front thickness, or vertical glacier height from the seafloor, was a significant covariate in all models with narwhals consistently visiting thicker glaciers (figure 2). Most glaciers are at approximately 90% flotation owing to the density of glacier ice, so this metric provides an estimated height of the submerged ice front face. In the 2000s, the front width also entered the models as a significant variable, with narwhals using wider (longer) ice fronts. The consistent use of thicker fronts and, when significant, wider ice faces may represent an attraction to ambient freshwater melt across the wall of underwater ice, with narwhals choosing maximal freshwater areas.

Surprisingly run-off, though included in some models, was never significant. When included, the relationship was negative, indicating narwhals prefer low subglacial run-off glaciers. Combined with the preference for thick fronts, the data suggest narwhals prefer glaciers with higher ambient melt from freshwater ice over glaciers with silt-laden discharge. Although research suggests subglacial discharge rises in buoyant plumes and increases glacier ice melt along the plume path [28], the subglacial discharge plumes may change water properties so they are not as attractive to narwhals as ambient melt.

Finally, a negative relationship with glacier velocity was included in several models but was often not significant. When included, narwhals used slower moving glaciers (low velocity). Glacier velocity represents both speed and iceberg calving activity (assuming a stable front location). Thus, use of lower velocity glaciers may suggest a preference for glaciers with less calving activity. Given use of thick glaciers, the preference for lower velocities was surprising. Thicker glaciers generally have higher velocities (e.g. [29]). High glacier velocities are also often associated with larger drainage basins, with larger subglacial discharge and fjord sediment flux, elements our models suggest narwhals select against. Our data suggest there may be unique glacier fjords preferred by narwhals—those with sufficiently thick ice fronts but low to moderate calving activity.

Ethics

Narwhal tagging was conducted under permits provided by the Greenland Government and IACUC protocol (no. 4155-01) from the University of Washington.

Data accessibility

The data used for this study have been deposited in Dryad as http://dx.doi.org/10.5061/dryad.ms812 [22].

Authors' contributions

K.L.L., M.P.H.-J. and R.D. conducted the fieldwork, K.L.L, T.M., D.D.W.H., and R.M. developed the database and analysed data, K.L.L. drafted the manuscript and integrated comments from all authors. All authors are accountable for this work and approved the final version for publication.

Competing interests

We have no competing interests.

Funding

Narwhal tagging was funded by Greenland Institute of Natural Resources and Greenland Environmental Research Institute in 1993 and 1994, NOAA Ocean exploration and DANCEA in 2006 and 2007. K.L.L. was funded by ONR grant no. N00014-11-1-0201.

Acknowledgements Thanks to local hunters in NW Greenland for assisting in narwhal captures. Thanks also to J. Willis, I. Fenty and K. Tinto for bathymetry data. Three reviewers improved the paper.

Footnotes

One contribution to the special feature ‘Effects of sea ice on Arctic biota’. Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3518946.