1. Introduction
Glaciers and ice sheets have been rapidly melting owing to rising temperatures, contributing a significant component (60%) to global sea level rise (Oppenheimer and others, Reference Oppenheimer, Portner, Roberts, Masson-Delmotte, Zhai, Tignor, Poloczanska, Mintenbeck, Alegria, Nicolai, Okem, Petzold, Rama and Weyer2019; IPCC, Reference Masson-Delmotte, Zhai, Pirani, Connors, Péan, Berger, Caud, Chen, Goldfarb, Gomis, Huang, Leitzell, Lonnoy, Matthews, Maycock, Waterfield, Yelekçi and Zhou2021). At tidewater glaciers, a significant component of the ice loss is attributed to submarine melting at the glacier–ocean interface, with estimates ranging from 6% up to 57% at different glaciers (Motyka and others, Reference Motyka, Hunter, Echelmeyer and Connor2003; Wagner and others, Reference Wagner2019). The submarine melting also drives some amount of calving, which forms the other significant part of the ice loss at the glacier termini (Motyka and others, Reference Motyka, Hunter, Echelmeyer and Connor2003; Fried and others, Reference Fried2015; Wagner and others, Reference Wagner2019). Thus, submarine melting is one of the leading processes driving tidewater glacier retreat (Slater and others, Reference Slater2019), and a key phenomenon that needs to be accounted for in climate-change assessments and projections. For example, in Greenland, where net mass loss from the ice sheet doubled over the last decade, accelerated ice retreat began at tidewater glaciers (Straneo and others, Reference Straneo2012) driven by rising submarine melting at the glacier termini (Straneo and Cenedese, Reference Straneo and Cenedese2015), suggested to be due to oceanic warming and atmospheric-driven ice-sheet melting (Slater and Straneo, Reference Slater and Straneo2022). Though significant work has gone into quantifying the submarine melt rate, it is still an open problem with many uncertainties, and its estimates are poorly bounded (Motyka and others, Reference Motyka, Hunter, Echelmeyer and Connor2003; Benn and others, Reference Benn, Warren and Mottram2007; O’Leary and Christoffersen, Reference O’Leary and Christoffersen2012). For example, in a recent study, estimates of the submarine melt contributions at Saqqarliup Sermia glacier in Greenland were noted to be at odds with the volume of icebergs and floating ice observed (Wagner and others, Reference Wagner2019), suggesting that this problem is far from fully understood. Furthermore, observations of meltwater intrusions in LeConte glacier in Alaska indicated that the submarine melting was two orders of magnitude higher than that predicted using existing theory (Jackson and others, Reference Jackson2020), and direct measurements of the melting at LeConte glacier in Alaska showed similar order of discrepancies with theoretical predictions (Sutherland and others, Reference Sutherland2019).
Ice loss mechanisms at glaciers emit distinct acoustic signatures, presenting an opportunity to use sound to monitor these events and tackle the uncertainty on their contributions (Schulz and others, Reference Schulz, Berger and Jansen2008; Pettit and others, Reference Pettit, Nystuen and O’Neel2012; Deane and others, Reference Deane, Glowacki, Dale Stokes and Pettit2019), which has spurred a burgeoning interest in this field of glacier cryoacoustics. Considerable research has been undertaken on the acoustics in polar regions, including in sea-ice covered (Milne and Ganton, Reference Milne and Ganton1964; Pritchard, Reference Pritchard1990; Kinda and others, Reference Kinda, Simard, Gervaise, Mars and Fortier2013; Dziak and others, Reference Dziak2015) and glacier environments (Pettit and others, Reference Pettit, Nystuen and O’Neel2012; Pettit, Reference Pettit2012; Glowacki and others, Reference Glowacki, Deane, Moskalik, Blondel, Tegowski and Blaszczyk2015; Deane and others, Reference Deane, Glowacki, Dale Stokes and Pettit2019) and that of icebergs (Urick, Reference Urick1971; Glowacki and others, Reference Glowacki, Deane and Moskalik2018). Passive-acoustic monitoring is a promising tool to understand the climate-change related mechanisms in glacial bays because it can provide (i) long-term continuous monitoring data with good temporal resolution on the order of a few seconds, (ii) information on ice–ocean interactions that is otherwise difficult to obtain, (iii) a large monitoring area-of-coverage, (iv) it can distinguish between the different mechanisms transpiring in glacial bays (Vishnu and others, Reference Vishnu, Deane, Chitre, Glowacki, Stokes and Moskalik2020) and (v) sound data are relatively simple to acquire compared to other data. Direct measurements at or near the termini of tidewater glaciers are dangerous due to calving activity and have been done using active acoustics (Fried and others, Reference Fried2015; Sutherland and others, Reference Sutherland2019) and marine robots (Kimball and others, Reference Kimball2014; Carlson and others, Reference Carlson2019; Howe and others, Reference Howe2019; Wagner and others, Reference Wagner2019; Bruzzone and others, Reference Bruzzone, Odetti, Caccia and Ferretti2020) which may need heavy investment for field campaigns and can be used for temporally sparse measurements. Large scale ice loss monitoring has been achieved so far using remote-sensing techniques such as satellite (Slater and others, Reference Slater2021), terrestrial or aerial photogrammetry (Kappas, Reference Kappas, Singh, Singh and Haritashya2011; Brun and others, Reference Brun2016) and seismology (Podolskiy and Walter, Reference Podolskiy and Walter2016) techniques, often coupled with model-based approaches using buoyant plume theory (Jenkins, Reference Jenkins2011) or oceanographic flux-gate methods (Motyka and others, Reference Motyka, Hunter, Echelmeyer and Connor2003; Jackson and Straneo, Reference Jackson and Straneo2016). Acoustic remote sensing can complement these existing methods and address some of their disadvantages, which may include limited time-resolution or spatial coverage, high costs, inability to operate in low-light/cloudy conditions or inability to directly monitor underwater ice loss mechanisms (Slater and others, Reference Slater2021).
Acoustics has already been proven successful in monitoring and quantifying calving fluxes and styles (Glowacki and others, Reference Glowacki, Deane, Moskalik, Blondel, Tegowski and Blaszczyk2015; Glowacki and Deane, Reference Glowacki and Deane2020; Glowacki, Reference Glowacki2022; Tegowski and others, Reference Tegowski2023). The other significant but less-understood component, submarine melting, leads to the release of ice-trapped air bubbles which, in turn, may enhance the melt activity (Wengrove and others, Reference Wengrove, Pettit, Nash, Jackson and Skyllingstad2023). This release generates a loud and distinctly observable crackling sound underwater similar to ‘frying bacon’ (Urick, Reference Urick1971; Pettit and others, Reference Pettit, Lee, Brann, Nystuen, Wilson and Neel2015), opening up the possibility of using passive acoustics to quantify submarine melting. The challenge herein is to unravel the complex relationship between the sound level and the submarine melt rate. This complexity arises due to the fact that tidewater glacial bays are time-varying environments with varying ice mélange cover and complex underwater acoustic channels whose thermohaline structure varies across seasons, thus altering the sound speed profile and introducing curvature in the acoustic paths, impacting sound propagation (Deane and others, Reference Deane, Glowacki, Tegowski, Moskalik and Blondel2014; Gł owacki and others, Reference Głowacki, Deane, Moskalik, Tȩgowski and Blondel2015; Pettit and others, Reference Pettit, Lee, Brann, Nystuen, Wilson and Neel2015; Glowacki and others, Reference Glowacki, Moskalik and Deane2016; Sanjana and others, Reference Sanjana, Latha, Thirunavukkarasu and Venkatesan2018; Reference Sanjana, Latha and Thirunavukkarasu2024). If this link between sound and melt activity can be deciphered, it can facilitate the development of an acoustic method to monitor submarine melting. Acoustics can then concurrently measure both calving and melting components of ice loss which excite different frequency bands in the acoustic spectrum. Therefore, this development can be used to form a suite of acoustic-based methods to fully decompose and monitor glacial ice loss, or moving forward, even quantify it.
Towards understanding the link between the sound level and melt activity, a series of studies on submarine melt-induced sound were undertaken focusing on Svalbard (Tegowski and others, Reference Tegowski, Deane, Lisimenka and Blondel2011; Deane and others, Reference Deane, Glowacki, Tegowski, Moskalik and Blondel2014) and Alaska (Pettit and others, Reference Pettit, Lee, Brann, Nystuen, Wilson and Neel2015). Svalbard has been warming six times faster than the global average (Wawrzyniak and Osuch, Reference Wawrzyniak and Osuch2020), and its temperature change is among the largest on Earth (IPCC, 2019), highlighting its increased sensitivity to climate change. Studies by Schuler and others Reference Schuler(2020) estimate the climatic mass balance to be 
$-7\pm4$ Gt a−1 across all Svalbard glaciers, with 
$-2\pm7$ Gt a−1 attributed to frontal ablation. A number of studies have concentrated on tidewater glaciers in Hornsund fjord, whose average retreat rate has been faster than the average retreat of tidewater glaciers across Svalbard (Blaszczyk and others, Reference Blaszczyk, Jania and Kolondra2013; Reference Blaszczyk2019). For example, the cumulative glaciological mass balance at Hansbreen in Hornsund was estimated to be higher than many other large area (
$ \gt 50\,\mathrm{km}^2$) glaciers (Schuler and others, Reference Schuler2020). This region has good research infrastructure and is well studied in terms of oceanography, hydrology and calving acoustics, complementing investigations into submarine melting acoustics and making it a good test-bed for additional studies.
In order to study the complex link between the underwater sound and submarine melt, we made acoustic recordings at four glaciers in Hornsund fjord, Spitsbergen. Here, we show that the recorded sound level exhibits a clear link with the melt activity. This link is robust, but complicated, due to the spatial variability in the sound level arising from several environmental factors. To remove the spatial variability, we develop a forward model describing the dependence of the sound level in the glacial bays on the location of measurement, thus capturing the effect of the environmental factors into a single mathematical framework. This allows us to decipher the spatial variability of the recorded sound level to a certain degree and interpret it in terms of the generated sound at the glacier terminus. The generated sound is shown to be correlated with the ablation rate at the terminus and the water temperature in the bay, establishing its link to the submarine melting. This analysis shows that the sound radiated by melting glacier ice provides a robust indicator of the melt rate at the glacier terminus and can be a promising signal to assess glacial melt rate variations.
2. Data and methods
2.1. Study area and experiment methodology
We study the underwater acoustic field in the bays of four glaciers in Hornsund fjord in June–July 2019: Hansbreen, Paierlbreen, Muhlbacherbreen and Samarinbreen (location shown in Fig. 1(a), photos in Fig. 1(b) and (c)). These are polythermal glaciers which have a mixed basal thermal regime (Glasser, Reference Glasser, Singh, Singh and Haritashya2011). Underwater sound was recorded from a boat using a vertical hydrophone array. It consisted of four HTI-96-MIN (High-Tech Inc., 2023) and two ITC6050C (ITC, 2023) hydrophones. These were mounted on a line and nonuniformly separated with an average separation of 2.63 m. The array was kept vertical by attaching a weight to the end of the line. Despite this, the array still exhibited a small amount of tilt due to the water flow drag. Depth sensors were mounted on the array line at two positions to detect and estimate this tilt—at 0.34 m below the surface mark and 0.34 m below the bottom-most hydrophone. The attachment of the hydrophones to the line was designed with sufficient isolation and mechanical damping using foam pieces, in order to reduce the effect of flow-induced strum noise on the array recordings. An illustration of the recording setup is shown in Figure 1(d).
We used a DR-680 Tascam digital recorder to capture audio at 96 kHz. During the trials, we made detailed logs of events occurring in the bay such as passing growlers or bergy bits that were close enough to contribute significantly to the recordings. Additionally, we logged the boat’s GPS location, measured the water temperature along the transect using a logger from RBR (RBR, 2022), took photographs of the bay supplemented with compass data and measured the range to the glacier using a laser-based device. Ranges were also measured to growlers and bergy bits that came close to the array. Conductivity-temperature-depth (CTD) casts were used to make estimates of the sound speed profile in the four bays (Fig. 2). These show the presence of a near-surface low sound speed layer due to the intrusion of freshwater from submarine melting and subglacial-discharge (illustrated in Fig. 1(d)). This occupies a surface layer reaching down to about 15–20 m below the surface and affects the acoustic propagation in the bay (Glowacki and others, Reference Glowacki, Moskalik and Deane2016; Vishnu and others, Reference Vishnu, Deane, Chitre, Glowacki, Stokes and Moskalik2020).
Figure 1. (a) The experiment location in Hornsund fjord, Svalbard, where recordings were made in the bays of four glaciers, (b) and (c) photos from field deployments showing the glacier cliff, floating ice mélange and growlers in the bay, and (d) illustration of the recording setup. The sound was recorded using a vertical hydrophone array deployed off a boat. The boat is allowed to drift outward from an initial starting point while recordings are made. At the glacier terminus, melt water resulting from submarine glacier melting and subglacial freshwater discharge gives rise to a glacially modified water layer on top. This water is colder and less saline than that at the bottom and has a different sound speed, which distorts the paths of melt-induced sound rays traveling from the glacier to the hydrophones.
Figure 2. (a) Temperature and (b) salinity, as determined on the Practical Salinity Scale, measured using a CTD sensor in the glacial bays on the dates nearest to the acoustic deployment dates. The warmest parts of the warm saline water layer extend from depths of at least 15 m upto 30–70 m across the four locations for the deployment dates considered.
We aim to assess the link of the submarine melting activity with a metric of the sound generated at the glacier termini, which is independent of the recorder’s location and depends only on the ice properties and the average melt rate. To estimate this quantity, we made recordings using the hydrophone array along linear transects in the glacial bays to observe the spatial variation of sound and obtain spatial diversity in the acoustic information. Within a glacial fjord, the starting points of transect were chosen to be different each time, as much as possible within the safety and accessibility conditions within the fjord, in order to get as much spatial diversity as possible in the acoustic data. During each of these transects, the boat was allowed to drift outward from the starting point while recordings are made (Fig. 1(d)), and the direction of the drift often varied across transects. Transect durations were planned to be roughly within 20–40 min, though this was sometimes varied depending on the safety and accessibility conditions within the fjord, and the transects were occasionally cut short when the array was surrounded by too much ice mélange. The transects are named as T1, T2, etc., where the number indicates the chronological order of their acquisition in each bay. The dates and glaciers at which the transect data were recorded are specified in Table 1.
Table 1. Dates in 2019 when acoustic measurement transects were undertaken in glacial bays, and start times and durations of each of the takes
2.2. Beamforming and acoustic summary
We first visualise the acoustic data recorded at the vertical array using a processing method known as beamforming (Johnson and Dudgeon, Reference Johnson and Dudgeon1992). Beamforming visualises the distribution of acoustic energy arriving at various elevation angles, in the form of a vertical angular profile of the recorded sound and how it evolves with time. This allows us to focus our attention on sound arriving from specific angles of interest. The distribution across elevation angles is computed based on the apparent velocity of the incoming sound waves at the different hydrophones of the array, which draws from an assumption of the sound speed based on measurements.
We use a frequency-domain Bartlett beamformer (Johnson and Dudgeon, Reference Johnson and Dudgeon1992) to assess the acoustic energy arriving from different elevation angles, assuming the waves have a planar wavefront. The beamformer output is computed within the 1–3 kHz spectral band because the submarine melt sound is highest in this band and tapers off at higher frequencies (Pettit and others, Reference Pettit, Lee, Brann, Nystuen, Wilson and Neel2015; Vishnu and others, Reference Vishnu, Deane, Chitre, Glowacki, Stokes and Moskalik2020). Frequency-domain beamforming allows us to focus on this narrow frequency band with low computational complexity.
The nonuniform spacing of the hydrophone elements in the line array ensured the array’s plane-wave Bartlett beamformer response yielded lower sidelobe levels than a uniformly spaced array, without incurring a significant increase in mainlobe width. The beamformer response (beam-pattern) at non-broadside angles, shown in Figure 3, is roughly 8 dB lower than at broadside, showing the rejection of sound from other directions during beamforming. When focusing on the energy incident horizontally at the array, the rejection of energy from non-horizontal directions is not perfect since the array aperture is short (Johnson and Dudgeon, Reference Johnson and Dudgeon1992), allowing some contributions from the growlers to leak into this horizontal band. Nevertheless, the 8 dB reduction of the energy from growlers allows us to better observe the sound from the glacier terminus and ice mélange.
Figure 3. Normalised plane-wave Bartlett beamformer response (beam-pattern) to a plane wave with frequencies spanning 1–3 kHz. This is computed for the nonuniformly spaced array deployed by us and compared against what would have been obtained if a uniformly spaced array with the same aperture (18.4 m) was deployed. The decrease in sidelobe level due to use of the nonuniform configuration is highlighted.
2.3. Ablation rate estimation
To compare the acoustic activity at the different glaciers against the melt activity, the ablation rates for the four glaciers are estimated using satellite photogrammetry-based observations of retreat. The rates are estimated by combining calving front retreat estimates and glacier velocity estimates obtained from Błaszczyk and others (Reference Błaszczyk2023) for Paierlbreen, Muhlbacherbreen and Samarinbreen and calving front retreat estimates from Li and others (Reference Li, Heidler, Mou, Igneczi, Zhu and Bamber2023); (Reference Li, Heidler, Mou, Igneczi, Zhu and Bamber2024) and velocity estimates from Blaszczyk and others (Reference Blaszczyk2024) for Hansbreen. The estimates in Błaszczyk and others (Reference Błaszczyk2023) were derived based on monthly averaged velocity mosaics from Friedl and others (Reference Friedl, Seehaus and Braun2021) and thus have a resolution of 1 month. The estimates in Li and others Reference Li, Heidler, Mou, Igneczi, Zhu and Bamber(2023) have an average temporal resolution of 4 days. Blaszczyk and others (Reference Blaszczyk2024) provide glacier velocities based on measurements from 11 survey stakes placed at different points along Hansbreen glacier, measured with resolutions of between a week to a month. The velocity estimate for Hansbreen is obtained using the velocity of stake 4 presented in this work (based on daily measurements), as it was closer to the calving front and located towards its centre. This velocity is multiplied by a factor of 3.6 to obtain the average calving front velocity, based on observations for 2015 detailed in Blaszczyk and others (Reference Blaszczyk2021), which uses a subset of the same stake velocity data. This estimate is also verified to be comparable to velocity estimates obtained from Friedl and others (Reference Friedl, Seehaus and Braun2021) which also provides calving front velocity estimates, albeit at lower resolution.
3. Results
3.1. Beamformer output
The beamformer outputs for all transects (of which one example is plotted in Fig. 4) reveal that there is a persistent arrival of sound energy at the array from a narrow angular band around the horizontal (0∘) (Vishnu and others, Reference Vishnu, Deane, Chitre, Glowacki, Stokes and Moskalik2020). These arrivals come from two types of sources: (i) submarine melting at the glacier terminus and (ii) melting of the ice mélange, a mixture of icebergs and smaller glacier ice pieces, scattered across the bay (Fig. 1(b)). Most of the submarine melt-induced sound production at the glacier is restricted to a narrow strip of limited depth (Vishnu and others, Reference Vishnu, Deane, Chitre, Glowacki, Stokes and Moskalik2020). This is because the expulsion of bubbles and consequent release of acoustic energy occurs due to the difference between the internal gas pressure in the ice and the hydrostatic pressure in the corresponding near-terminus water region, which quickly reduces with depth below the water surface. The depth limit of significant acoustic activity was estimated to be about 12.8 m in controlled experiments with glacier ice in Svalbard (Vishnu and others, Reference Vishnu2023). This sound from the submarine melting at the glacier and the ice mélange propagates underwater and arrives at the array at nearly horizontal angles, which results in the band in the beamformer output centred around 0∘ (Vishnu and others, Reference Vishnu, Deane, Chitre, Glowacki, Stokes and Moskalik2020).
Figure 4. Beamformer output plotted against trial time and elevation angles at which sound arrives at the array, for data recorded at Hans glacier, T7. The output is expressed in decibels with reference to µPa
$^2/^\circ$, which shows the amount of acoustic power per degree at each time point. Trial time is measured from the start of the recording. The narrow band of angles around the horizontal (
$-5^\circ$ to 5∘) mostly contains energy from the direct and surface-reflected paths from the glacier terminus and ice mélange in the bay. The green lines represent limits of the band within which the beamformer output is integrated to form the acoustic summary of this transect. The output also reveals the presence of a growler floating nearby during the recording, leading to an energy band in the output with time-varying elevation angles. Since it was close to the array at the beginning, the melt-induced sound from the growler arrives at the array at large positive elevation angles (band of energy at 15∘–45∘ in the plot). As the recording progresses, the growler moves away from the array, resulting in the gradual decrease in the amplitude and elevation angle of its energy contributions.
Apart from the melt signal, there may also be contributions to the recording from growlers or bergy bits that float near the array, as seen in the example output in Figure 4. Growlers are not always present near the array during recordings, but if they are, their melting contributes significantly to the recorded sound level (Vishnu and others, Reference Vishnu, Deane, Chitre, Glowacki, Stokes and Moskalik2020), potentially masking the signal from the glacier terminus and complicating its analysis. Here, we can use the capability of beamforming with a vertical array of hydrophones to disentangle this sound from nearby growlers. Since this sound arrives at the array at steeper angles than the sound from the glacier terminus, beamforming allows us to exploit this to extract the sound contributions from the terminus alone and reject that from nearby growlers. It also allows us to account for the effect of the unknown and often-complicated sea-floor topography on the recorded sound, because sound reflected off the sea-floor arrives at steeper elevation angles. Accordingly, we use beamforming to reduce the relative contribution from growlers and bottom-reflected sound and focus on the glacier terminus melt-induced sound arriving at near-horizontal directions, and observe its spatial variation clearly. This horizontal band of energy is our best estimate of the signal from submarine glacier melting observed at a given recording location, henceforth referred to as the melt signal.
3.2. Melt signal
After beamformer-based processing, we compute the melt signal as a function of time by numerically integrating the beamformer output over elevation angles within a horizontal band using the trapezoidal rule. By doing so, we obtain an acoustic summary plot of the recording made during the transect. The angular extent of the band depends on the glacial bay geometry and the array’s location. The melt signal level is calibrated such that it matches the average power across the hydrophones in transects when there were minimal contributions from growlers, as shown in Figure 5.
The variation of the melt signal against the range of the array from the closest point on the glacier is shown in Figure 6(a) for all transects. The melt signal is directly dependent on the sound generated at the ice–water interfaces in the bay, which, in turn, is dependent on
1. The ice properties at each glacier that influence the acoustics (Grossman and others, Reference Grossman, Johnson, Stokes and Deane2024) and
2. The average melt rate in the bay.
Figure 5. Scatter plot of calibrated melt signals against the average power across hydrophones. The beamformer output is calibrated such that the melt signal (integrated beamformer output within the horizontal band) matches the average recorded power across hydrophones for transects where there is little or no contribution from growlers and bottom-reflected energy. The calibration is done by multiplying the output by a constant factor that ensures the melt signal matches the average recorded power in such cases. From the plot, the calibrated melt signal roughly matches the average hydrophone power for two such transects where growler contributions and bottom-reflected energy are minimal, namely Hans T1 and Paierl T1. For two other transects during which growlers were present, namely Paierl T4 and Samarin T2, the melt signal is lower than the average hydrophone power because the beamformer suppresses contributions from the growlers using the array’s directionality.
Figure 6. (a) Summary plots of the melt signal versus range from the glacier for all transects. Each line represents one transect. The three colours used for Hans represent different times—yellow at the beginning of the campaign, cyan in the middle of the campaign and green at the end of the campaign, when the temperature was higher. The melt signal levels in different transects are clustered depending on the glacier, due to the combined influence of several glacier-specific factors. (b) The MSI plotted for all the transects during the field campaign (in dB normalised with respect to the median of the Muhlbacher T3 curve, plotted with a 15 dB span on y-axis similar to (a)). The vertical spread in these curves is reduced compared to (a), showing that the modelling is able to account for some of the variability.
To understand whether the spread in melt signals seen in Figure 6(a) at the different glaciers can be attributed to the melt activity, we assess the ablation rates at the different glaciers. The ablation rates are found to vary by a factor of at most ∼5 within the timeframe of interest. Recall that the dependence of the melt signal on the melt rate is expected to be linear (Glowacki and others, Reference Glowacki, Moskalik and Deane2016)—i.e. a doubling of the melt rate should double the melt signal (increase it by 3 dB in decibel scale)—thus, variation in melt rate alone does not fully explain the spread observed in the melt signals at different glaciers, which differ by up to ∼22 times (13.4 dB). Furthermore, as the array drifts away from the glacier, Figure 6(a) shows that the melt signal varies but does not decrease uniformly with the range to the glacier. These observations are attributed to many additional factors other than the melt rate and geometrical spreading from the closest point on the glacier which considerably modulate the variation of the melt signal resulting from the generated sound, including:
3. The geometry and span of the glacier front, because the recorded sound contains contributions from across the submerged glacier face.
4. The presence and coverage of melting ice mélange.
5. The array location with respect to the sound sources, namely, the glacier terminus, growlers and the ice mélange.
6. The sound speed profile in the underwater acoustic channel between the source and array, which is determined primarily by its thermohaline properties. During summer, the channel typically has an upward-refracting profile due to a near-surface cold freshwater layer composed of melt water and subglacial-discharge of freshwater resulting from atmospheric melting (Vishnu and others, Reference Vishnu, Deane, Chitre, Glowacki, Stokes and Moskalik2020; Slater and Straneo, Reference Slater and Straneo2022). This freshwater layer acts like a lens, ducting sound through it (Glowacki and others, Reference Glowacki, Moskalik and Deane2016), and overlies a warm saline water layer, as shown by our thermohaline measurements (Fig. 2). The warmest parts of this layer were observed to extend from depths of at least 15 m down to 30–70 m in the water column.
All these factors are taken into account in the modelling (next section) to explain the variability in the melt signal.
3.3. Modelling melt signal variation
To interpret the melt signal (Fig. 6(a)) in terms of melt rate variations, we first need to understand and separate out the variations caused by the location-dependent factors mentioned previously. We do this by starting from a canonical formulation fundamentally used in ambient noise oceanography (Deane, Reference Deane1997) (Eqn (A7) in Appendix A), which describes the contribution of a radiating elemental surface area 
$\mathrm{d}A$ to the power spectrum of the recorded acoustic pressure at the hydrophone at a location r0, as
\begin{equation}
\mathrm{d}P(\omega, \textbf{r}, \textbf{r}_0) = M\left(\textbf{r}, \kappa \right) \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}A,
\end{equation} where 
$\omega = 2\pi f$ is the angular frequency. 
$M(\textbf{r}, \kappa)$ is a term describing the sound generated by the source at a location r and depends on source parameters κ. In the current context, these source parameters include the source depth, ice properties and the melt rate (which, in turn, depends on water temperature T). 
$G(\omega, \textbf{r}, \textbf{r}_0)$ is the frequency-domain Green’s function that describes the effects of underwater sound propagation from the source to the hydrophone and its dependence on acoustic channel properties such as the sound speed profile.
Equation (1) is sufficiently general that it can be adapted to the geometry and propagation factors relevant for the glacial melt problem. This is done in Appendix A—we start from Eqn (1) and develop a forward model that mathematically describes the melt signal variation in terms of all the controlling factors discussed in the previous section. The mathematics of the problem allows us to express the final model for the melt signal spectrum in terms of two separate terms. One is the ‘source term’ describing the sound generated due to melting, which we are interested in, and the other encompasses the location-dependent variations in the signal. The model for the melt signal recorded at array location ra is given by Eqn (A23), reproduced here:
\begin{equation}
P_a(\omega, \textbf{r}_a) = U(\omega, T) G_a(\omega, \textbf{r}_a),
\end{equation} where 
$U(\omega, T)$ is the source term, which depends on the melt rate—specifically, the source power spectrum per unit area at 1 m from the glacier surface recorded at a reference depth of 0 m. We define the melt-source intensity (MSI) in Eqn (A27) as an integral of this quantity within 1–3 kHz. The second term, referred to as the melt area-integrated Green’s function (MGF), models the melt signal’s spatial variations within each glacial bay due to propagation effects and its dependence on the spatial distribution of melting ice, the array location, and the sound speed profile. In the context of melt assessment, this is a ‘nuisance term’ capturing the many location-dependent factors that add complexity to the problem and confound the melt signal analysis, which we need to mitigate eventually. Thus, it is paramount that we carefully model this term so that it incorporates the physical effects of each of these factors as accurately as possible in order to aid our analysis moving forward. This rigorous exercise is undertaken (in Appendix A), where we incorporate the complex co-interacting effects of the multiple environmental factors into a single mathematical quantity describing the location-dependent variation as
\begin{align}
G_a(\omega, \textbf{r}_a) &= \underbrace{\frac{1}{M}\sum_{j=1}^{M}}_\text{average over sensors}\left[ \underbrace{\int_{A_1} \overbrace{\alpha(z)}^\text{depth-dependence} \overbrace{\left|G(\omega, \textbf{r}, \textbf{r}_j)\right|^2}^\text{Green's function} \mathrm{d}\textbf{r}}_\text{glacier terminus term}\right.\nonumber\\ &\quad + \underbrace{\int_{A_2} \alpha(z) \left|G(\omega, \textbf{r}, \textbf{r}_j) \right|^2 \mathrm{d}\textbf{r}}_\text{growler term} \nonumber\\ &\quad +\left. \underbrace{\int_{A_3} \alpha(z) \left|G(\omega, \textbf{r}, \textbf{r}_j) \right|^2 \mathrm{d}\textbf{r}}_\text{ice m\'{e}lange term} \right].
\end{align} The MGF accounts for the melt signal’s variation due to three types of sound sources - the glacier terminus (incorporating the effect of its geometry), growlers and ice mélange, with submerged melting areas A 1, A 2 and A 3, respectively. 
$\alpha(z)$ describes the depth-dependence of the sound generated. Equations (2) and (3) together describe the effect of all the factors discussed in the previous section, combined into one single model.
Evaluating Eqns (2) and (3) involves substantial groundwork in order to obtain measurements or estimates of the several environmental parameters describing the experiment at each location. To do this, we obtained CTD measurements in the glacial bays (Fig. 2) to evaluate the sound speed profile of the underwater channel, aerial photographs and experimental logs to evaluate ice mélange cover and the presence of growlers, satellite-aided estimates of the glacier geometry, and GPS in conjunction with a rangefinder to estimate the array location within the bay. These measurements are used in Eqns (2) and (3), and an acoustic propagation model (Chitre, Reference Chitre2023) is used to evaluate the Green’s function for each transect of data recorded (details in Appendix B). Finally, the MGF-predicted variations for all transects are plotted together in Figure 7.
Figure 7. Melt-area integrated Green’s function versus range from the glacier. The plot’s y-axis span is the same as that of the melt signal plot in Figure 6(a). The MGF captures some of the small-scale variability in the melt signal, though the relative levels between transects may not be the same because of differing melt rates. For example for Muhlbacher T1, we observe a steep decrease in the beginning till 290 m, followed by a gentle decrease. For Hans T1, there is little or no variation until 340 m, followed by a slight increase, and then it flattens out again. For Hans T8, there is a sharp increase in the beginning due to a growler’s presence, followed by a sharp decrease when the growler recedes from the array, and the decrease becomes gentler when the growler is far away. These variations are visible in the MGF too.
In Figure 7, there is a 5 dB (factor of 3) spread in the MGF across different transects in the four bays. This shows that there is a significant variation in melt signal arising only from location-specific factors. The modelled MGF is able to capture some of the small-scale variability in Figure 6(a) and explains the nonuniform variation observed in melt signal with increasing distance. In other words, it is able to characterise the effect of the different co-interacting environmental factors that together determine the complex spatial variation. Some features of the variability are not captured by the MGF, which could be due to (i) noise contributions from sources other than melting, (ii) not fully incorporating the presence of the ice mélange or growlers, (iii) incomplete representation of the acoustic medium or (iv) inaccurate assumptions. For example, the two spikes in the sound levels for Muhlbacher T1 at ranges of 275 m and 280 m are due to ice hitting the boat. The spike in Samarin T2 at 390 m is attributed to a calving event. While it is not straightforward to fully account for the uncertainty arising due to such aspects, it is partly characterised and quantified in Appendix C and incorporated into the forthcoming discussion of results.
3.4. Estimating MSI from melt signal
We now assess the MSI, an acoustic metric representing the sound generated by the melt activity. This quantity is independent of the recorder’s location and depends only on the ice properties and the average melt rate. Having modelled the MGF, we can estimate the MSI by subtracting (in decibels) the model-predicted variation in the melt signal—this is equivalent to shifting the term Ga in Eqn (2) to the left side to estimate U. The estimated MSIs for all the recorded transects are plotted in Figure 6(b) (the prequalifier ‘estimate’ is dropped henceforth for brevity). We see that removing the location-dependent variations reduces the vertical spread in the melt signal levels by a factor of 2, from 13.4 dB in Figure 6(a) to 10.3 dB in Figure 6(b).
We now summarise Figure 6 into Figure 8 which plots the median MSI and melt signal at 275 m from the glacier against the (a) water temperature recorded in each transect at 26 m depth, (b) maximum temperature across the depth of the water column recorded in the bay and (c) the calving front ablation rate. Deployments T1 and T2–T5 are grouped together as they happened at the same location on successive days. We also plot the best linear fits for the MSI plots. The temperatures and ablation rate are used here as proxies for the melt activity—the temperature is one of the main drivers of melt activity, and the ablation encompasses the net effect of melting and calving. The maximum water temperature across the depth is indicative of the temperature of the warm saline water layer. The uncertainty in each of the points is also represented in Figure 7(b).
Figure 8. Plots showing correlation of median melt-source intensity and melt-signal with water temperature and ablation rate. Circles show median MSI, and squares show the melt signal at a range of 275 m from the glacier, or the range closest to 275 m. This reference range value allows us to obtain data points from most transects. The statistical significance of the linear regression coefficient is tested using a F test utilising the Python ‘statsmodels’ package. A P-value <0.05 is considered statistically significant. (a) The median MSI (in dB with respect to the minimum point) across each transect plotted versus median water temperature along the transect, (b) average of the transect-wise median MSI for each deployment day plotted versus the maximum temperature across depth obtained via CTD measurement in the corresponding bay, (c) average of the median MSI versus the ablation rate. Water temperature measurements shown in (a) were made at 26 m depth, which lies within the warmest band of the water column comprising of incoming waters. Note that the temperatures in (a) and (b) were measured at different depths and locations in the bay. Both the MSI and melt signal are plotted with a 13 dB range on the y-axis. The whiskers indicate the uncertainty in terms of ±1 standard deviation. The overall spread in MSI around the linear fit (5.4 dB in (a) and 3.2 dB in (c)) is smaller than that in the melt signal (7.6 dB in (a) and 4.9 dB in (c)), showing the improvement from removing the location-dependent variations.
4. Discussion
Figure 8 shows that the MSI/melt signal are correlated with the water temperature and ablation rate, which provides a two-pronged validation of the acoustic technique for monitoring melt activity. There are, however, deviations in these points from a monotonically increasing curve, which arise from effects that can broadly be classified into three categories. Before we further discuss the correlation of the MSI/melt signal, we discuss these categories to provide necessary subtext to understand the plots:
1. The first category encompasses effects arising from the fact that we use the mean ablation rate per day and water temperature in the bay, as proxies for the melt rate. We use these to assess the effectiveness of the acoustic method because they are currently our best available options. However, while the ablation rate includes the effect of submarine melting, it also includes contributions from calving and subglacial discharge. There are two limitations in using water temperature as a proxy for melting. (i) The values we measured are not fully representative of the water temperature fields at the ice–water interfaces which drive the melting. There is likely a complex nondeterministic relationship between the two, and the temperatures at the interfaces are likely to be highly variable locally. (ii) Submarine melting at the glacier is affected by many variable local co-interacting factors beyond the temperatures alone, including glacier shape/morphology, water circulation (Jenkins, Reference Jenkins2011), salinity and bubbles in the water column (Wengrove and others, Reference Wengrove, Pettit, Nash, Jackson and Skyllingstad2023). The resultant submarine melting is usually nonuniform across the glacier face—for example, it is highest in the vicinity of subglacial discharge plumes (Motyka and others, Reference Motyka, Hunter, Echelmeyer and Connor2003; Fried and others, Reference Fried2015; Wagner and others, Reference Wagner2019), and the difference between regions with and without subglacial discharge plumes may differ by up to an order of magnitude or more as shown at Greenland glaciers (Fried and others, Reference Fried2015; Wagner and others, Reference Wagner2019). Thus, the melting is not driven just by a single water temperature value, but rather a mix of several factors, which are not considered in the plot.
2. The second category stems from uncertainty in the MSI estimation technique and the parameters used in it. This could be attributed to (i) uncertainty in the data and (ii) uncertainty in the model. Uncertainty in the data could arise, for example, due to noise contributions from sound sources that are unaccounted for. Uncertainty in the model arises due to incomplete representation of the underwater medium parameters and their variability. These deviations could also partly be attributed to the difference in ice properties across glaciers. The vertical whiskers in Figure 8 quantify some of the uncertainty effects in this category and are thus partly indicative of the magnitude of this uncertainty (see Appendix C).
3. The third category stems from the fact that the MSI is indicative only of the melt activity occurring down to a limited depth at the glacier terminus which generates the sound, and it is not certain to what degree this is indicative of the melt across the full terminus depth.
The horizontal whiskers in Figure 8 depict the variability in water temperature within each of the transects where the acoustic recordings were made. Clearly, even this within-transect variability is considerable. While these do not adequately characterise the first category of variability effects owing to the complex link between the measured water temperatures and the submarine melt rate, they depict significant local variability in the factors that drive melting.
The results in Figure 8 show a distinct robust correlation between the melt signal and (a) the median water temperature along the array transect (regression coefficient 
$R_\text{melt signal}$= 0.55), (b) maximum water temperature across depth (
$R_\text{melt signal}$ = 0.79) and (c) ablation rate (
$R_\text{melt signal}$ = 0.75), despite the limitations and sources of variability discussed above. Furthermore, the correlations of the MSI are generally higher than that for the melt signal—the MSI versus transect temperature regression coefficient is R MSI = 0.70, which is statistically significant (P = 0.0018), and a coefficient of R MSI = 0.85 is obtained for the MSI versus ablation rate variation (P = 0.03). The improved regression coefficients for MSI show that the incorporation of environmental factors has reduced variability in the melt signal, suggesting that the MGF modelling is successful.
5. Conclusion
In this study, we detail the measurement and analysis of the sound radiated due to submarine melting at glacier termini, made using a vertical hydrophone array in Hornsund fjord, Svalbard. We develop a physics-based model to explain and interpret the sound level variation in the bay and estimate the generated sound at the glacier terminus using the model developed and the data obtained. The analysis provides a two-pronged validation of the acoustic technique in terms of the correlation of the recorded sound level and generated sound with two physical parameters—(1) the water temperature which is one of the main drivers of submarine melting, and (2) the net frontal ablation rate which comprises the melt rate as one of its significant components.
The melt signal’s correlation with the melt-activity proxies implies that it can be used to monitor melting activity, if the limitations of the method are taken into account. Going a step further, the higher correlations of the MSI with the melt activity proxies show that the MGF is a valuable tool when interpreting the melt signal, and that the melt signal combined with MGF shows promise for monitoring submarine melting based on passive acoustic records. Overall, these observations support the conclusion that the sound radiated by melting submerged glacier ice at the terminus provides a robust indicator of the melt rate and is a promising signal to assess the melt rate provided its spatial variations can be satisfactorily accounted for. The study presented here also provides a necessary framework for using acoustic data to estimate submarine melt rates in the future, facilitating long-term, large area-of-coverage monitoring crucial for climate-change studies. The acoustic technique for monitoring is uniquely positioned because the acoustic signal is a direct outcome of melting and does not require an understanding or measurement of the many complex and locally varying interacting parameters that lead to or result from melting. This makes a strong case for using acoustic monitoring to complement other approaches.
The next step forward in the evolution of this technique is to develop practical long-term monitoring systems with simpler recording setups, which will form a part of our future work. While a single-hydrophone recording setup has obvious advantages in terms of simplicity and ease of deployment for long-term monitoring, an array with vertical aperture is advantageous in disregarding interfering noise from growlers and the sea floor, as we have showcased in this work. A limitation in the vertical-array acoustic data is that it can only discriminate sounds in the vertical direction and not the horizontal direction. Thus, neither can it distinguish between sound coming from the glacier versus that from the opposite direction of the glacier, nor can it capture the variability in melt horizontally across the glacier face, observed in several earlier studies (Motyka and others, Reference Motyka, Hunter, Echelmeyer and Connor2003; Fried and others, Reference Fried2015; Wagner and others, Reference Wagner2019). Since there is no way for us to distinguish this aspect of spatial variability, we only provide an estimate of the melt-induced sound aggregated across the melting ice–water interfaces using these data, with the assumption that the melt-noise is the dominating component in the data, based on our observations. This may not necessarily hold at other glaciers and fjord settings and would need to be accounted for. This limitation can be tackled in the future by using an array that has both horizontal and vertical aperture in order to distinguish both the horizontally arriving sound as well as the spatial variability along the glacier front.
Other steps that can contribute towards building a robust and effective monitoring system include
1. Unraveling the dependence of sound level on glacial ice properties—e.g. the bubble density ρ and average spectrum of a melting bubble β discussed in Eqn (A27). These parameters could vary significantly between glaciers, or even at different locations on the same glacier, potentially complicating our understanding of the melt signal and application of the acoustic technique and modelling. Hence, understanding these would aid in generalising the acoustic technique to different glaciers. Additional controlled experiments and physics-based modelling towards understanding these aspects were undertaken in the campaign, with the objective of an improved understanding that can help in this regard (Grossman and others, Reference Grossman, Johnson, Stokes and Deane2024).
2. Better characterisation of the environmental parameters at the monitored glacial bays, which will also allow us to refine the sound-propagation modelling and improve the accuracy of the analysis. This would also make the method more robust to differences in the environment.
Apart from monitoring the temporal trends in melt activity, the acoustic technique could potentially be used to estimate the melt rate from the MSI, which would also require an accurate independent verification of this estimate. This estimation is not possible with the current dataset and state of knowledge of the physics of melting at the glacier terminus but will be pursued in future work in this area.
The acoustic technique can be used for monitoring autonomously over a long term in the harsh glacial environments with good time-resolution and can shed light on the temporal variability within or across seasons. Passive acoustics is already in use widely for long-term monitoring in several applications (Thode and others, Reference Thode, Skinner, Scott, Roswell, Straley and Folkert2010; Caruso and others, Reference Caruso2020; Lamont and others, Reference Lamont2022) including polar monitoring (Kinda and others, Reference Kinda, Simard, Gervaise, Mars and Fortier2013; Reference Kinda, Simard, Gervaise, Mars and Fortier2015; Glowacki and others, Reference Glowacki, Deane, Moskalik, Blondel, Tegowski and Blaszczyk2015; Matsumoto and others, Reference Matsumoto2019; Mahanty and others, Reference Mahanty2020), which means that the technology for achieving long-term recordings for glacial melt monitoring is not impractical to achieve. Recent breakthroughs in passive acoustics via techniques such as distributed acoustic sensing using fiber-optic cables (Bouffaut and others, Reference Bouffaut2022) have the potential to further enable long-term monitoring over a large area in glacial regions (Walter and others, Reference Walter2020) with more spatial information that can enhance the melt-activity monitoring. The observations from such long-term recordings could serve as valuable inputs to climate-change assessments and projections, complementing existing data and improving the accuracy of these assessments.
Acknowledgements
This work was possible due to the support of grants from the National Research Foundation, Prime Minister’s Office, Singapore under its Marine Science Research and Development Program MSRDP P-42, US Office of Naval research grant ONR-N00014-24-1-2394, National Science Foundation grant number OPP-1748265, the National Science Centre, Poland (grant 2021/43/D/ST10/00616), the Institute of Geophysics, Polish Academy of Sciences Grant 2a/IGF PAN/2017, and a subsidy for the Institute of Geophysics, Polish Academy of Sciences provided by the Ministry of Science and Higher Education of Poland. We thank the staff at the Polish Polar Station, Hornsund, for all the support they have provided and for maintaining the oceanographic monitoring with conductivity-temperature-depth profiles.
Competing interests
The authors declare no competing interests.
Appendix A. Modelling of sound from submarine melting—theoretical framework
We begin based on the approach in Glowacki and others Reference Glowacki, Moskalik and Deane(2016). Figure 1 shows an illustration to understand the model development. Consider a small volume element of ice located at r containing escaping bubbles which are the sound sources, and the sound is recorded by a hydrophone located at r0. At a time tj, the volume-element radiates a sound signature 
$s(t - t_j, z, a_j, \gamma_j)$ due to bubble j, which propagates through the underwater channel. The sound generated depends on the depth z of the volume-element, radius aj of bubble j, and γj which denotes properties of the glacier ice in the volume-element that affect the sound signature. We do not detail these properties but note that the crystalline structure of ice may be one such property that affects the release timescale of bubbles from the ice and thus its sound signature (Grossman and others, Reference Grossman, Johnson, Stokes and Deane2024).
The effects of geometrical spreading, attenuation, refraction and boundary scattering on the sound wave as it propagates from r to r0 are characterised by an impulse response function 
$g(t,\textbf{r},\textbf{r}_0)$. Assuming the ocean channel can be treated as a linear, time-invariant system, the acoustic pressure at the hydrophone due to bubble j is given by Glowacki and others Reference Glowacki, Moskalik and Deane(2016)
\begin{equation}
\zeta_j(t,\textbf{r}, \textbf{r}_0) = \int_{-\infty}^{t} s(\tau - t_j, z, a_j, \gamma_j) g(t - \tau, \textbf{r}, \textbf{r}_0) \mathrm{d}\tau.
\end{equation} If there are 
$N(\textbf{r})$ bubbles emitting noise at various times tj within the volume-element during some observation interval τc, then the total contribution to the noise field due to bubbles released from the volume-element is given by the sum
\begin{equation}
\mathrm{d}p(t,\textbf{r}, \textbf{r}_0) = \sum_{j = 1}^{N(\textbf{r})}\int_{-\infty}^{t} s(\tau - t_j, z, a_j, \gamma_j) g(t - \tau, \textbf{r}, \textbf{r}_0) \mathrm{d}\tau,
\end{equation} where τc is assumed to be long enough that many events have occurred but not so long that the geometry of the problem has changed significantly within this time span. We are concerned with sound-emitting bubbles on the glacier wall or from the melting surfaces of floating ice. Denote by 
$v_t(\textbf{r}, T)$ the glacier terminus melt rate defined as the velocity at which the calving front recedes through melting, which depends on the in situ water temperature T. vt also depends on other factors not detailed in this analysis for the sake of simplicity, which potentially includes local water circulation, plume freshwater flux, effects of bubbles in the water column, and glacier shape and morphology (Jenkins, Reference Jenkins2011). Let 
$\rho(\textbf{r}, a, \gamma)$ denote the bubble density (number of bubbles per unit volume) which can be described using the product of marginal probability density distributions ρa and ργ that describe bubble radius and ice properties respectively:
\begin{equation}
\rho(\textbf{r}, a, \gamma) = \rho(\textbf{r}) \rho_a(a) \rho_\gamma(\gamma),
\end{equation}where
\begin{equation}
\int_{V}^{} \rho(\textbf{r}) \mathrm{d} \textbf{r} = n(\textbf{r}),\\
\int_{B}^{} \rho(a) \mathrm{d}a = 1, \text{and}\\
\int_{\Gamma}^{} \rho(\gamma) \mathrm{d}\gamma = 1,
\end{equation}where 
$n(\textbf{r})$ is the number of bubbles per unit volume at r, B and Γ denote the bubble radii and glacier ice properties spans of interest and V represents a unit volume.
Let the volume-element correspond to an elemental area 
$\mathrm{d}A$ on the glacier terminus that is melting at 
$v_t(\textbf{r}, T)$ m s−1 over the observation time τc. Then
\begin{equation}
N(\textbf{r}) = \rho(\textbf{r}) v_t(\textbf{r}, T) \tau_c \mathrm{d}A \rho(\gamma) \mathrm{d}\gamma \rho(a) \mathrm{d}a,
\end{equation} where 
$\mathrm{d}A$ is an elemental surface area on the ice exposed to water. The form of the dependence of vt on T is not explored here, but we note that previous works have suggested an algebraic dependency on the difference between the temperature and the freezing point (Slater and Straneo, Reference Slater and Straneo2022) or even an approximately linear dependence within the plume (Jenkins, Reference Jenkins2011). Substituting (A5) into (A2),
\begin{align}
\mathrm{d}p(t,\textbf{r}, \textbf{r}_0) = \sum_{j = 1}^{\rho(\textbf{r}) v_t(\textbf{r}, T) \tau_c \mathrm{d}A \rho(\gamma) \mathrm{d}\gamma \rho(a) \mathrm{d}a}\int_{-\infty}^{t} s(\tau - t_j, z, a_j, \gamma_j) g(t - \tau, \textbf{r}, \textbf{r}_0) \mathrm{d}\tau .\nonumber\\
\end{align} If the bubble release events are assumed to be occurring randomly and independently during the acoustically active period, so that the event times tj are Poisson-distributed (Papoulis, Reference Papoulis1984), then by Carson’s theorem (Rice, Reference Rice1944) the power spectrum of 
$\mathrm{d}p$ is
\begin{equation}
\mathrm{d}P(\omega, \textbf{r}, \textbf{r}_0) = M\left(\textbf{r}, \kappa \right) \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}A,
\end{equation} where 
$\omega = 2\pi f$ is the angular frequency, 
$G(\omega, \textbf{r}, \textbf{r}_0)$ is the frequency-domain Green’s function and M describes the sound generated, which depends on a set of source parameters κ including source depth z, ice properties γ, bubble radius a and melt rate. M is further decomposed as
\begin{equation}
M\left(\textbf{r}, T, \kappa \right) = \rho(\textbf{r}) v_t(\textbf{r}, T) S(\omega, z, a, \gamma) \rho(\gamma) \mathrm{d}\gamma \rho(a) \mathrm{d}a
\end{equation} where 
$S(\omega, y, a, \gamma)$ is the power spectrum of bubble releases. The acoustic power spectrum recorded at r0 can be expressed as
\begin{equation}
P(\omega, \textbf{r}_0) = P_1(\omega, \textbf{r}_0) + P_2(\omega, \textbf{r}_0) + P_3(\omega, \textbf{r}_0),
\end{equation}where P 1, P 2 and P 3 correspond, respectively, to contributions from submarine melting at the glacier terminus, large growlers floating in the bay and ice mélange. Together, these account for most of the sound power in the 1–3 kHz band (see illustration in Fig. 1). Contributions due to calving in this spectral band are minimal. P 1 can be expressed as an area integral of (A7) over the melting face of the glacier:
\begin{equation}
P_1(\omega, \textbf{r}_0) = \int_{A_1}^{} \rho(\textbf{r}) v_t(\textbf{r}, T) \bar{S}(\omega, z) \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}\textbf{r}
\end{equation} where 
$\bar{S}(\omega, z) = \int_{\Gamma}^{} \int_{B}^{} S(\omega, z, a, \gamma) \mathrm{d}a \mathrm{d}\gamma$ is the source spectrum average over the bubble radii and glacier properties of interest. Let 
$\textbf{r} = (x, y, z)$ denote the position in terms of easting x, northing y and depth z (with respect to the water surface), respectively, where x and y are measured with respect to some arbitrary origin (without loss of generality). Assuming the glacier face is vertically flat in z and can be represented by a curve W in x–y coordinates, and H is the depth of the glacier front face below the sea surface, we have
\begin{equation}
P_1(\omega, \textbf{r}_0) = \int_{W}^{} \int_{0}^{H} \rho(\textbf{r}) v_t(\textbf{r}, T) \bar{S}(\omega, z) \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}\textbf{r},
\end{equation}where the first integral is the line integral over the curve W. Previous work (Vishnu and others, Reference Vishnu2023) has shown that the power spectrum P of sound generated by bubbles escaping from melting ice can be expressed as the product of separate functions that depend on ω, water temperature T and source depth z. This allows us to decompose the source spectrum in (A11) into a product of depth- and frequency-dependent terms, as
\begin{equation}
\bar{S}(\omega, z) = \alpha(z) \beta(\omega),
\end{equation}where 
$\alpha(z)$ indicates the depth-dependence of the spectrum normalised with respect to a certain depth (set as 0 m here) and β is the average source spectrum of bubble releases at a depth of 0 m.
Though the bubble density and melt rate may vary across the glacier face (Wagner and others, Reference Wagner2019), there is no way for us to determine these variations using the vertical array used, and we are only interested in an average value of these parameters that lead to the resultant melt signal. The second mean-value theorem of integrals guarantees the existence of an average spectrum-weighted bubble release rate per unit area given by
\begin{equation}
\eta_1(T) = \frac{\int_{W}^{} \int_{0}^{H} \rho(\textbf{r}) v_t(\textbf{r}, T) \bar{S}(\omega, z) \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}z~\mathrm{d}x~\mathrm{d}y}{\int_{W}^{} \int_{0}^{H} \bar{S}(\omega, z) \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}\textbf{r}}.
\end{equation}It is worth noting that η 1 is weighted over the horizontal span and depth of acoustic activity at the terminus. In practise, contributions to this term come only from the region of the glacier termini that produces sound, and thus, its direct interpretation is limited to the melt activity at shallow depths, unless further efforts to expand the region of interpretability are taken into account. Substituting (A13) and (A12) in (A11),
\begin{equation}
P_1(\omega, \textbf{r}_0) = \eta_1(T)\beta(\omega) \int_{W}^{} \int_{0}^{H}\alpha(z) \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}\textbf{r}.
\end{equation}For melting growlers, we assume the form of the melt-spectrum is similar to that of the ice at the glacier terminus. The sound power contribution from this component for the ith growler can similarly be expressed as an area integral
\begin{equation}
P_{2,i}(\omega, \textbf{r}_0) = \eta_{2,i}(T_{2,i}) \beta(\omega) \int_{A_{2,i}} \alpha(z) \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d} \textbf{r},
\end{equation} where 
$\eta_{2,i}(T_{2,i})$ denotes the average spectrum-weighted bubble release rate per unit area for the ith growler dependant on the effective water temperature at the growler 
$T_{2,i}$, and the integral is over the submerged surface area of the growler. The total contribution from all growlers, then, is expressed as
\begin{equation}
P_{2}(\omega, \textbf{r}_0) = \beta(\omega) \sum_{i = 1}^{N_2} \eta_{2,i}(T_{2,i}) \int_{A_{2,i}} \alpha(z) \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}\textbf{r},
\end{equation}where N 2 is the number of growlers.
Likewise, for melting ice mélange, if we assume the form of the average melt-spectrum is similar to that at the glacier terminus, the noise power contribution from this component can be expressed as
\begin{equation}
P_3(\omega, \textbf{r}_0) = \eta_3(T_{3}) \beta(\omega) \int_{A_3} \alpha(z) \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}A,
\end{equation}where 
$\eta_3(T_3)$ and T 3 denote the average bubble release rate per unit area and effective water temperature for the ice mélange, and the integral is over the submerged surface area of the ice mélange. We assume the ice mélange can be treated as quasi-2-D (Burton and others, Reference Burton, Amundson, Cassotto, Kuo and Dennin2018), with the exposed submerged area sheet being at a constant depth of y 3 (which can be set depending on the nominal thickness of the surface ice).
Thus, the total power spectrum of the recorded data at the hydrophone can be expressed as a sum of (A14), (A16) and (A17), as
\begin{align}
& P(\omega, \textbf{r}_0) = \beta(\omega) \left[ \eta_1(T) \int_{W}^{} \int_{0}^{H} \alpha(z) \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}\textbf{r} \right.\\
& + \sum_{i = 1}^{N_2}\eta_2(T_{2,i}) \int_{A_{2,i}} \alpha(z) \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 d\textbf{r} + \eta_3(T_{3}) \alpha(z_3)\nonumber\\ & \quad\times\left. \int_{A_3} \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}\textbf{r} \right].\nonumber
\end{align} To evaluate the contributions from the different factors, detailed measurements of the growler sizes and shapes and ice mélange coverage areas would be needed. However, these parameters are not all accurately known in the current dataset collected. For the sake of throwing light on the relative contributions of the three factors, we consider some simplifying assumptions. For the contribution from the ice mélange, we first denote a horizontal bounding box within which the ice mélange exists, with its northwestern corner at 
$(x_3, y_3)$ and with easting width W 3 and northing length L 3. Then the contribution from the ice mélange can be expanded as
\begin{align}
&P_3(\omega, \textbf{r}_0) = \eta_3(T_{3}) \beta(\omega) \alpha(z_3) \int_{y_3}^{y_3+L_{3}} \int_{x_3}^{x_3+W_{3}} I(x, y) \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}x~\mathrm{d}y ,
\end{align} where 
$I(x, y)$ is an indicator which is 1 if there is ice mélange at northing-easting location (x, y) and 0 otherwise. It is not easy to obtain 
$I(x, y)$ over the glacial bay, but we may estimate a statistical average of the coverage from photographs of the bay during deployments. From the second mean-value theorem of integrals, we know that there exists a weighted indicator value 
$\bar{I}$ such that
\begin{align}
&\bar{I} = \frac{\int_{y_3}^{y_3+L_{3}} \int_{x_3}^{x_3+W_{3}} I(x, y) \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}x~\mathrm{d}y}{\int_{y_3}^{y_3+L_{3}} \int_{x_3}^{x_3+W_{3}} \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}x~\mathrm{d}y} .
\end{align} Clearly, 
$0\le \bar{I} \lt 1$, and a smaller coverage of ice mélange within the bounding box should lead to a smaller value of 
$\bar{I}$. Substituting (A19) and (A20) into (A18), we get
\begin{align}
& P(\omega, \textbf{r}_0) = \beta(\omega) \eta_1(T) \left[ \int_{W}^{} \int_{0}^{H} \alpha(z) \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}\textbf{r} + \sum_{i = 1}^{N_2}k_{2,i} \right.\nonumber\\ & \int_{A_{2,i}} \alpha(z) \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}\textbf{r} \\
& \left. + k_3 \bar{I} \alpha(z_3) \int_{y_3}^{y_3+L_{3}} \int_{x_3}^{x_3+W_{3}} \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}x~\mathrm{d}y \right],\nonumber
\end{align} with 
$k_{2,i} = \frac{\eta_{2,i}(T_{2,i})}{\eta_1(T)}$ and 
$k_3 = \frac{\eta_3(T_{3})}{\eta_1(T)}$. We are not in a position to exactly determine 
$k_{2,i}$ and k 3, but we do have some intuition about what these constants mean. For example, if the area-normalised bubble release rate of all growlers and ice mélange (which depends on their melt rates) is equal to that at the glacier terminus, 
$k_{2,i}= 1$ and 
$k_3 = 1$. The average acoustic power recorded across the J hydrophones of the array as a function of its 2-D location ra in the glacial bay can be modelled as
\begin{align}
\tilde{P_a}(\omega, \textbf{r}_a) = \frac{1}{J} \sum_{j= 1}^{J} P(\omega, \textbf{r}_j)
\end{align}where rj is the 3-D position vector of the jth hydrophone. Recall that the melt signal is calibrated to be equal to the acoustic-power averaged across hydrophones when the sound is incident at the array only at horizontal angles (Fig. 5). Thus, (A22) also models the melt signal when the latter is true. However, to model the melt signal when a nontrivial percent of the recorded acoustic energy may be incident at non-horizontal angles, further modifications will now be incorporated into this equation.
The horizontal band in the beamformer output contains contributions from only direct and surface-reflected raypaths arising from submarine melting at the glacier terminus and the ice mélange. Thus, the modelling henceforth will focus on only these raypaths, and the Green’s function will be taken to represent only this subset of possible sound raypaths (and discard bottom-reflected raypaths). The acoustic contribution from nearby growlers is incident at angles mostly outside the horizontal band. However, a small amount of this energy leaks into the horizontal band even after spatial filtering, as quantified by the beam-pattern (Fig. 3). For the 
$i\text{th}$ growler with location ri and height 
$H_{2,i}$, the arrival angle of sound at the hydrophone array can be estimated from the geometry. The growler sound’s leakage into the horizontal band is represented by the beam-pattern sidelobe level for this arrival-angle, denoted as 
$B(\textbf{r}_i,\textbf{r}_a, H_{2,i})$. Using this, a model for the melt signal modified from (A22) by incorporating beamforming leakage is
\begin{equation}
P_a(\omega, \textbf{r}_a) = U(\omega, T) G_a(\omega, \textbf{r}_a),
\end{equation}where 
$U(\omega, T) = \beta(\omega) \eta_1(T)$, and 
$G_a(\omega, \textbf{r}_a)$ is summarised as
\begin{align}
G_a(\omega, \textbf{r}_a) &= \frac{1}{J}\sum_{j=1}^{J} \left[\underbrace{\int_{A_1} \alpha(z)\left|G(\omega, \textbf{r}, \textbf{r}_j)\right|^2 \mathrm{d}\textbf{r}}_\text{glacier terminus term} + \underbrace{\int_{A_2} \alpha(z) \left|G(\omega, \textbf{r}, \textbf{r}_j) \right|^2 d\textbf{r}}_\text{growler term} \right. \nonumber\\ & \quad +\left. \underbrace{\int_{A_3} \alpha(z) \left|G(\omega, \textbf{r}, \textbf{r}_j) \right|^2 \mathrm{d}\textbf{r}}_\text{ice m\'{e}lange term} \right], \text{where}
\end{align}
\begin{align*}
& \int_{A_1} \alpha(z)\left|G(\omega, \textbf{r}, \textbf{r}_j)\right|^2 \mathrm{d}\textbf{r} = \int_{W}^{} \int_{0}^{H} \alpha(z) \left|G(\omega, \textbf{r}, \textbf{r}_j) \right|^2 \mathrm{d}\textbf{r},
\end{align*}
\begin{align*}
& \int_{A_2} \alpha(z) \left|G(\omega, \textbf{r}, \textbf{r}_j) \right|^2 \mathrm{d}\textbf{r} = \sum_{i = 1}^{N_2} B(\textbf{r}_i,\textbf{r}_a,H_{2,i}) \int_{A_{2,i}}^{} \alpha(z) \left|G(\omega, \textbf{r}, \textbf{r}_j) \right|^2 \mathrm{d}\textbf{r}, \\
& \int_{A_3} \alpha(z) \left|G(\omega, \textbf{r}, \textbf{r}_j) \right|^2 \mathrm{d}\textbf{r} = \bar{I} \alpha(z_3) \int_{y_3}^{y_3+L_{3}} \int_{x_3}^{x_3+W_{3}} \left|G(\omega, \textbf{r}, \textbf{r}_0) \right|^2 \mathrm{d}x~\mathrm{d}y.
\end{align*} To summarise the acoustic metric within the frequency band of interest, we consider the integral of 
$P_a(\omega, \textbf{r}_a)$ within the melt-spectrum between 1 and 3 kHz given by
\begin{equation}
\int_{2\pi\times1000}^{2\pi\times3000} P_a(\omega, \textbf{r}_a) \mathrm{d}\omega = \int_{2\pi\times1000}^{2\pi\times3000} \beta(\omega) \eta_1(T) G_a(\omega, \textbf{r}_a) \mathrm{d}\omega.
\end{equation} 
$G_a(\omega, \textbf{r}_a)$ is found to be approximately flat within the 1–3 kHz band, allowing us to drop its dependence on ω and simplify the above equation as
\begin{equation}
\int_{2\pi\times1000}^{2\pi\times3000} P_a(\omega, \textbf{r}_a) \mathrm{d}\omega = U(T) G_a(\textbf{r}_a),
\end{equation}where
\begin{align}
U(T) = \eta_1(T) \int_{2\pi\times1000}^{2\pi\times3000}\beta(\omega) \mathrm{d}\omega,
\end{align}is the sound power spectrum per unit area at reference 1 m from the melting ice surface at a depth of 0 m integrated within the melt frequency band, which we refer to as the MSI. 
$G_a(\textbf{r}_a)$ is the MGF and models the spatially varying component of the melt signal due to the combined effect of the channel properties, depth-dependence of source level, glacier front geometry and location of ice in the bay. Moving 
$G_a(\textbf{r}_a)$ to the left hand side of (A23) translates to removing the location-dependant variations in the melt signal to estimate the MSI U(T), which was undertaken in Section 3.
As per this model, the MSI depends on the 
$\eta_1(T)$ which, in turn, depends on the water temperature. It also depends on 
$\beta(\omega)$ which represents the average power spectrum of sound generated by bubbles in a glacial fjord, which, in turn, depends on additional parameters that may be specific to the glacier such as the distribution of bubble radii and ice properties.
Appendix B. Modelling parameters
We compute the Green’s functions using a ray-tracing propagation code (Chitre, Reference Chitre2020; Reference Chitre2023), which accounts for the effect of sound speed profiles on the propagation. We consider only direct, surface-reflected and ice-reflected rays in the modelling, since the data processing considers only the horizontal band and does not include any bottom-reflected raypaths. The Green’s function is computed for a frequency of 2 kHz which is the centre of the spectral band where the melt spectrum is strongest (Pettit and others, Reference Pettit, Lee, Brann, Nystuen, Wilson and Neel2015; Vishnu and others, Reference Vishnu, Deane, Chitre, Glowacki, Stokes and Moskalik2020). The reflection-coefficient of the water surface is set to 0.91, estimated as the specularly reflected sound component from a randomly rough surface with a nominal root-mean-square roughness of 5 cm (Jensen and others, Reference Jensen, Kuperman, Porter and Schmidt2011). The reflection-coefficient from the ice wall is computed according to the relations for a liquid–solid interface (Jensen and others, Reference Jensen, Kuperman, Porter and Schmidt2011) based on ice acoustic properties detailed in Glowacki and Deane Reference Glowacki and Deane(2020). The geometry of the glacier wall is computed based on digitised satellite photographs obtained from Sentinel-2A via the USGS EROS Archive (Earth Resources Observation and Science (EROS) Center, 2020).
The values of 
$k_{2,i}$ and k 3 are unknown. We assume the melt rate of growlers is approximately equal to that at the glacier terminus and set 
$k_{2,i} = 1$ for all i. Nominal values of 
$A_{2,i}$ are chosen based on knowledge of existing growler dimensions, aerial photos taken at the field, and logs describing the ice growlers, and the acoustic signatures of these growlers recorded at the array. We treat k 3 as a free variable with small values within 
$[0,1)$ and set it based on the logs and photos taken at the trials.
Appendix C. Quantifying uncertainty in results
Here, we discuss how the uncertainty in Figure 8 is characterised incorporating as many factors as possible. Firstly, there is uncertainty in the water temperature, which varies across the bay. We do not assess what the water temperature may be at the glacier front itself since this region is inaccessible and may show large local fluctuations (Motyka and others, Reference Motyka, Hunter, Echelmeyer and Connor2003). The temperature values shown in Figure 8 are the median of the measurements obtained from the temperature sensor along each of the transects when mounted at a constant depth. The uncertainty shown in temperature values (horizontal whiskers) for points corresponding to each glacial bay is the standard deviation of temperature measurements over all the transects in each bay. This is representative of the spatial variability in the water temperature in regions in the bay at least 150 m away from the glacier front, which delivers heat to the terminus and is one of the main factors driving the melting.
Uncertainty in the MSI estimate (Eqn (A27)) may arise from its two components—the melt signal and the MGF. Melt signal uncertainty may arise from sound contributions from sources other than those considered, namely submarine glacier melting, and melting of growlers and ice mélange. There is no reliable method to measure this using the data available, so we estimate this uncertainty as the standard deviation in the estimated MSI from its mean value along each transect (Fig. 6(b)). The intuition here is that if the MGF had fully accounted for all the sound sources leading to the melt signal, it would have successfully removed all location dependent variations, and thus each of the lines in Figure 6(b) would have been flat. Thus, any residual deviation in these lines from a flat line is attributed to unaccounted-for sources leading to variations in the recorded sound.
Uncertainty in the MGF (Eqn (A24)) may arise from uncertainty in the Green’s function term Ga and in the depth-dependence term 
$\alpha(z)$. The depth-dependence was measured in a previously conducted controlled experiment (Vishnu and others, Reference Vishnu2023) which used different trials with glacier ice blocks. The uncertainty in 
$\alpha(z)$ is modelled by assuming it is Gaussian distributed with a standard deviation 
$\sigma_d(z)$, which is set equal to the standard deviation in the measured depth-dependences across all the trials in the study (Vishnu and others, Reference Vishnu2023). Uncertainty in Ga arises due to the possible deviation in environmental parameters from those assumed in the modelling, such as (1) sound speed profile, (2) reflection coefficients from water surface and (3) contributions from ice mélange and growlers and their sizes. We characterise this in three steps:
1. A Monte-Carlo method is used to characterise the effect of the uncertainty in sound speed profile. We run 500 different realisations of the ray-tracing code (Chitre, Reference Chitre2023) for the source-receiver ranges considered, with the sound speed profile in each realisation generated randomly via a Gaussian process kernel (Rasmussen, Reference Rasmussen, Bousquet, von Luxburg and Rätsch2004), varying from the one measured in the glacial bay by a certain standard deviation σc. σc is estimated as the standard deviation between 12 sound speed profiles measured from CTD casts at different locations in Hans bay on
$27\text{th}$ June. This value is assumed to be representative of the horizontal variation in sound speed profiles within the bay and turns out to be
$\sigma_c = 0.668$ m s−1. This procedure yields the uncertainty in amplitudes of each sound raypath attributed to uncertainty in sound speed.2. Uncertainties in growler and ice mélange sizes are modelled by assuming these values have a nominal standard deviation of 10% around the value estimated based on photographs and logs.
3. Uncertainty in the reflection coefficients may arise due to variability in the roughness of the water–air or water–ice boundaries which leads to sound raypaths being reflected by varying amounts. This uncertainty is modelled by assuming that the root-mean-square roughness is Gaussian distributed with a mean and standard deviation of 5 cm. The reflection coefficients are then estimated as the specularly reflected component of sound incident at a randomly rough boundary (Jensen and others, Reference Jensen, Kuperman, Porter and Schmidt2011), with the boundary roughness distributed as mentioned above.
Finally, we combine all the above components together with the melt signal uncertainty by using a linear approximation of error propagation theory (Ku, Reference Ku1966; Clifford, Reference Clifford1973) via the ‘uncertainties’ package (Lebigot, Reference Lebigot2023) to compute the uncertainty in the MSI. Uncertainty in the ablation rates is not plotted as it is not characterised in the sources based on which it is estimated.
Open access
