DAS to discharge: using distributed acoustic sensing (DAS) to infer glacier runoff

Abstract

Observations of glacier melt and runoff are of fundamental interest in the study of glaciers and their interactions with their environment. Considerable recent interest has developed around distributed acoustic sensing (DAS), a sensing technique which utilizes Rayleigh backscatter in fiber optic cables to measure the seismo-acoustic wavefield in high spatial and temporal resolution. Here, we present data from a month-long, 9 km DAS deployment extending through the ablation and accumulation zones on Rhonegletscher, Switzerland, during the 2020 melt season. While testing several types of machine learning (ML) models, we establish a regression problem, using the DAS data as the dependent variable, to infer the glacier discharge observed at a proglacial stream gauge. We also compare two predictive models that only depend on meteorological station data. We find that the seismo-acoustic wavefield recorded by DAS can be utilized to infer proglacial discharge. Models using DAS data outperform the two models trained on meteorological data with mean absolute errors of 0.64, 2.25 and 2.72 m³ s⁻¹, respectively. This study demonstrates the ability of in situ glacier DAS to be used for quantifying proglacial discharge and points the way to a new approach to measuring glacier runoff.

1. Introduction

Glaciers are an important yet diminishing reservoir of freshwater for communities and ecosystems (Casassa and others, 2009). In the European Alps, for example, modeled future trends indicate a large reduction or disappearance of glaciers on decadal timescales due to climate change (Haeberli and others, 2007; Linsbauer and others, 2013; Zekollari and others, 2019). Glaciated catchments provide a river discharge buffering mechanism, particularly important during the dry season. This mechanism will likely be disrupted if alpine glaciers continue to retreat and to disappear (Mark and Seltzer, 2003) with immediate effects on the downstream ecology, which is particularly susceptible to changes in glacier-sourced freshwater input to proglacial streams (Cauvy-Fraunié and others, 2016). In addition, hydroelectric power production is expected to decrease within the century as a substantial part of the current hydroelectic power is produced by unsustainable glacier mass loss caused by the warming climate (Schaefli and others, 2019). As infrastructure grows and glaciers retreat, it will become increasingly important to measure or infer glacier melt runoff, to help to accurately predict its contribution to the catchment's freshwater resources on seasonal and diurnal timescales.

Glacier surface melt is the primary contributor to the mid latitude glacier hydrological system (Shreve, 1972). However, it remains difficult to observe the dominant processes that drive surface melt with sufficient spatial and temporal resolution (Landmann, 2022). Conventional in situ methods for measuring glacier surface ablation include ablation stakes (Fountain and Vecchia, 1999; Pratap and others, 2015; Landmann and others, 2021) and the use of meteorological data to calculate energy fluxes that result in glacier surface melt (Braithwaite, 1995; Hanna and others, 2005; Lenaerts and others, 2019). Although ablation stake measurements and reconstruction from meteorological station data are foundational methods, they come with the significant disadvantage of being labor intensive and therefore difficult to implement glacier-wide, long-term studies. Satellite remote sensing, in contrast, offers the only feasible way to monitor glacial melt at a global scale. A wide variety of remote-sensing methods have been used to infer glacier surface melt indirectly through observed changes in glacier elevation (Markus and others, 2017; Sutterley and others, 2018), mass (Wouters and others, 2008) or surface backscatter (Ridley, 1993; Trusel and others, 2013; Bevan and others, 2018). Although satellite remote sensing may offer true global coverage, it oftentimes lacks the spatial or temporal resolution required to resolve rapid, local variations in surface melt (Yang and Smith, 2013; Yang and Li, 2014; Wille and others, 2019). More fundamentally, even when remote sensing of glacier surface melt is able to attain a desired spatial and temporal resolution (Trusel and others, 2013; Bevan and others, 2018; Sutterley and others, 2018), such platforms nevertheless benefit from – and in many cases require – in situ observations for calibration and validation. Advances in satellite remote sensing of glacier melt therefore motivate the need for improved in situ observations of glacier surface melt.

The familiar variety of sounds associated with flowing water attests to the ubiquity of flow-induced acoustics. A correspondingly large number of previous studies have examined the seismo-acoustic wavefield generated by water flow. Basic physical processes implicated in the generation of sound from flowing water include wave breaking (Manasseh and others, 2006), hydraulic jump formation (Ronan and others, 2017), low-frequency fluid pulsing in conduits (Podolskiy, 2020) and the entrainment and collapse of air bubbles in turbulent flows (Prosperetti, 1988; Morse and others, 2007). In terrestrial rivers, both discharge and bedload transport contribute to the seismic wavefield (Burtin and others, 2008, 2011; Gimbert and others, 2016; Roth and others, 2016, 2017; Cook and others, 2018), as do roughness elements such as boulders (and resulting rapids) and engineered blocks and weirs (Schmandt and others, 2013; Osborne and others, 2021, 2022). In glaciers, flow in subglacial conduits is constrained by conduit size with an observable impact on the seismic wavefield (Bartholomaus and others, 2015; Nanni and others, 2020).

Here, we utilize distributed acoustic sensing (DAS) to record the seismo-acoustic wavefield originating from turbulent supraglacial water flow. The sensing component of DAS is a single mode optical fiber cable deployed on the surface of the glacier. The principle of DAS is that the phase shift of Rayleigh-backscattered light in an optical fiber is used to infer the fiber axial strain rate with spatial resolution on the order of several tens of centimeters and at frequencies, dependent on cable length, of millihertz to several kilohertz (Shatalin and others, 2021), therefore enabling observation of seismo-acoustic wavefields (Lindsey and Martin, 2021; Douglass and others, 2023). Fluid flow velocities within pipes have been estimated using regression of DAS data (Vahabi and others, 2020; Titov and others, 2022). Several studies have previously described glacier surface (Walter and others, 2020; Hudson and others, 2021) and borehole DAS deployments (Booth and others, 2023) for investigating the en- and subglacial environment. Here, we leverage DAS observations from a 9 km long optical fiber deployed along the flow line of an alpine glacier to examine the relationship between glacier melt and the in situ glacier surface seismo-acoustic wavefield.

2. Field site and data

2.1 Rhonegletscher

Our measurements were conducted at Rhonegletscher, a temperate mountain glacier located in the central Swiss Alps, in the summer of 2020 (Fig. 1a). The glacier covers a total area of 15.5 km² and ranges from 3600 m above sea level (a.s.l.) to 2200 m a.s.l. at its terminus with a length of about 8 km (GLAMOS, Bauder and others, 2020). During the field study, the surface of Rhonegletscher in the accumulation zone primarily consisted of firn (Fig. 1b). The ablation zone was characterized by bare ice, crevasses and distributed supraglacial meltwater streams (Fig. 1c).

Figure 1. (a) Map of the study site. Approximate path of the fiber optic cable deployment and location of the distributed acoustic sensing (DAS) interrogator including outline of Rhonegletscher (Consortium, 2005). Orthophoto provided from the Swiss Federal Office of Topography. (b) Photo of the glacier surface and deployed cable in the accumulation zone (credit: Małgorzata Chmiel), consisting mostly of firn at the time of deployment (July 2020). (c) Photo of the glacier surface and deployed cable in the ablation zone (credit: Sara Klaasen), consisting primarily of bare ice with areas of crevassing, meltwater surface streams, meltwater pools and glacier moulins.

2.2 Distributed acoustic sensing (DAS) deployment

A Silixa iDAS$^{\rm TM}$ interrogator was deployed in a tent west of the terminus of Rhonegletscher from 4 July 2020 to 4 August 2020. A 9 km single-mode fiber optic cable was laid out on the surface of the glacier approximately along the glacier flow line spanning across ablation and accumulation zones. During the first portion of the experiment, interrogator recording settings such as channel spacing and sampling rate were varied for instrument and sensitivity testing. Starting on 13 July, settings remained constant for the remainder of the experiment. To avoid complexities with different instrument settings, in this study we only use the data from 13 July to 4 August 2020. During this time, data were recorded continuously at 1 kHz sampling frequency, 4 m channel spacing and 10 m gauge length over 2496 channels. At this sampling frequency, cable length and gauge length, the iDAS$^{\rm TM}$ is sensitive to 2 picostrain per square root Hertz. The last 188 channels contain instrument noise only, because the actual fiber optic cable length was shorter than the length set in the interrogator settings. Thus, we only use the first 2308 channels for our analysis. For most of our analysis, we high-pass filtered the data above 50 Hz. In later analysis, we investigate the unfiltered DAS data to determine the influence of the broad band spectrum on discharge prediction. The high-pass filter also mitigates the effects of thermal expansion with a diurnal period (Klaasen and others, 2021), shading from transient and local cloud cover, and from other anthropogenic sources (Huynh and others, 2022) such as nearby hydropower production causing narrow-banded seismic energy at 16.7 and 50 Hz. For each channel, we calculated the root mean square (RMS) of the fiber strain-rate for each 30 s window of each channel in the DAS data (Fig. 2a).

Figure 2. (a) DAS time series over analysis period. Data are high-pass filtered above 50 Hz and normalized to peak RMS strain rate over all channels per time step. Low channel numbers are located closest to the terminus down glacier (i.e. closer to the interrogator) and higher channel numbers are located progressively up glacier according to the plotted cable layout in Figure 1a. The dashed line denotes roughly the transition from the ablation zone down glacier and the accumulation zone up glacier. (b) Rhône river discharge recorded about 3 km downstream of the proglacial lake. During the final 2 d of the experiment, a standing wave formed in the proglacial stream in the location of the discharge measurement resulting in the three crest pattern that is evident. (c) Hourly temperature and precipitation data from 10 min recordings at Grimsel Hospiz meteo station (Swiss Federal Office of Meteorology and Climatology MeteoSwiss).

2.3 Discharge measurements

During summer, meltwater from Rhonegletscher is the primary contributor to the highest reaches of the Rhône river near Oberwald, Switzerland. A radar-based discharge gauge (Swiss Federal Office for the Environment, station ID number 2268) located in Gletsch about 3 km downstream of Rhonegletscher's proglacial lake recorded hourly averaged discharge of the Rhône river throughout the duration of DAS data collection. Discharge data (Fig. 2b) were linearly interpolated to 30 s to match the 30 s RMS time steps calculated from the raw DAS data.

2.4 Meteorological measurements

We used meteorological data from the station Grimsel Hospiz (Swiss Federal Office of Meteorology and Climatology MeteoSwiss) located 5–8 km southwest of Rhonegletscher behind a mountain ridge. Temperature data were collected at 10 min intervals and precipitation data were recorded as the sum over the 10 min period (Fig. 2c).

3. Machine learning models

3.1 Architectures: linear, neural network, long short-term memory

In order to quantify the relationship between glacier melt and the recorded glacier surface seismo-acoustic wavefield, we employ three separate machine learning (ML) models using Keras TensorFlow (Martin and others, 2015) and assess their relative performance. We first implement a linear model with a single dense layer with linear activation. This model mostly serves as a baseline point of comparison with two more flexible models. Second, we implement a Neural Network (NN) model with two dense layers containing 32 units and a rectified linear unit activation function each, a flattening layer and a dense layer with one unit. Finally, we implement a Long Short-Term Memory (LSTM) model with a single LSTM layer containing 32 units and a dense layer with one unit. The features (independent variables) in our analysis consist of the multivariate time series of DAS strain rate data. The labels (dependent variables) in our analysis consist of the measured discharge values from the downstream discharge gauge. These models are each associated with learning rate, batch size and data input window size hyperparameters; we choose these hyperparameters based on the results of 90 experiments per model (see Fig. S1). As a result of the analysis, we choose a learning rate of 0.001, a batch size of 32 feature-target pairs, a window size of 200 time steps as these parameters produced stable and robust results. The Supplemental Information further describes hyperparameter tuning.

3.2 Cross-validation scheme

Previous studies of changes in supraglacial hydrology through space and time (Nicholson and others, 2021, e.g.) demonstrate that the surfaces of glaciers are inherently non-stationary over the timescale of several weeks during the melt season. Supraglacial stream geometry changes throughout the melt season and responds to change in water flow (Germain and Moorman, 2019). For this reason, we randomly shuffled the time series windows used for inputs prior to data separation into training, validation and test sets. We therefore ensure that all possible glacier surface melt regimes occurring during the observation period are captured in the model training dataset. In addition to shuffling, we use standard cross-validation (CV) techniques (Bishop and Nasrabadi, 2006, Chapter 14.2) wherein we perform 100 model trainings, each with a uniquely seeded test/training split. CV allows us to quantify model sensitivity to input data and estimate the non-stationary effect of the glacier surface on model performance.

3.3 Meteo-LSTM model

We consider an intermediate complexity, ‘Meteo-LSTM’ model that uses an LSTM model architecture with temperature and precipitation data as features and discharge as labels. The goal of this model is to understand the impact of model complexity versus the underlying usefulness of different datasets by testing a model which has similar complexity to the DAS-LSTM model but only relies on the meteo station data.

3.4 Positive degree-day (PDD) model

PDD models are widely used to infer glacier melt from limited meteorological observations (Braithwaite, 1984). We implement a PDD model following Hock (Hock, 2005). We carry out a minimization analysis to select the melt rate factor and lapse rate value that resulted in the lowest absolute error in discharge. Temperatures as collected at Grimsel Hospiz were corrected over elevation bands of 100 m. Then the discharge prediction at each elevation band was summed to get the final predicted discharge,

(1)$$D = \sum_{z = 2.3\ {\rm km}}^{3.6\ {\rm km}} \left\{\matrix{ \left\{\left[( T + \gamma \;( z - z_0) \right]\,f + P\right\}A \hfill & T > 0 \hfill \cr PA \hfill & T\leq 0 \hfill \cr \hfill }\right.$$

where z is the altitude, z ₀ is the terminus altitude, D is the total predicted discharge, T is temperature, γ is the calibrated lapse rate, f is the calibrated melt factor, P is the precipitation rate and A is the area of the glacier within each step in the summation. The glacier area is given as an idealized rectangle with the glacier area, width and elevation range as found in GLAMOS (GLAMOS, Bauder and others, 2020). The PDD model results were interpolated to match the times of discharge measurements used as LSTM model targets. In order to compensate for meltwater transport from the proglacial lake to the discharge gauge downstream, which is evident from the phase lag between a basic PDD model and measured discharge curves, the PDD model results were shifted based on the phase of maximal cross-correlation between modeled and observed discharge.

4. Results

The results of our analysis are listed in Table 1. For all of our analyses, we present results in terms of the mean absolute error (MAE) of the residuals and the standard deviation of the residuals between model outputs and discharge gauge measurements. All of these performance statistics are reported for the test dataset in order to quantify model performance when evaluated on data that were not used for parameter estimation. Overall, the best performing models use an LSTM architecture with input DAS data. These models perform about 40% better than the NN model in terms of MAE. The LSTM models also result in a more than 200 times reduction in MAE compared to a linear model.

Table 1. Model types and mean absolute error (MAE) for test dataset

We plot the estimated discharge time series and residuals from our model (Figs 3a–f, respectively). Examination of these time series confirms that the DAS-LSTM model is able to capture the phase of discharge (Fig. 3a). In contrast, the Meto-LSTM and PDD models suffer from both poor amplitude and phase response (Figs 3b,c).

Figure 3. (a) DAS-LSTM model ensemble mean (red dashed) line and confidence interval (grey region) from cross-validation (CV). (b) same as (a), but with the meteo-LSTM model. (c) Positive degree-day (PDD) model results. (d–f) Residuals for the DAS-LSTM, Meteo-LSTM and PDD models, respectively.

Model residuals for the DAS-LSTM model show no systematic relationship with increasing discharge (Fig. 3d). The Meteo-LSTM model, in contrast, shows both poor amplitude and phase response (Fig. 3d) which is likely due to the poor correlation between temperature and precipitation amplitude and phase. The PDD model estimates reasonable amplitudes with a phase shift. We therefore calculate PDD residuals using a best fit time shift. Residuals for the PDD model are uncorrelated with increasing discharge and an order of magnitude larger than the residuals from the DAS-LSTM model.

4.1 Ablation zone versus accumulation zone

Models trained on ablation zone data performed better and have less variance than models that only used accumulation zone data. Models trained on data from the ablation zone have a mean MAE of 0.64 m³ s⁻¹ and standard deviation of 0.1 m³ s⁻¹ whereas models trained on accumulation zone data have a mean MAE of 1.07 m³ s⁻¹ and standard deviation of 0.24  m³ s⁻¹. This can also be seen in the sensitivity analysis discussed in the Discussion section and shown in Figure 4 where particular sectors in the ablation zone generally show higher sensitivity to discharge than areas in the accumulation zone.

Figure 4. Channel sensitivity analysis from applying a uniform in time Gaussian pulse with a width of 50 channels. A new discharge prediction is made each time the Gaussian pulse is centered on the next channel. The mean prediction is calculated from the predicted discharge of the 100 LSTM models produced. Predictions are given in values of a normalized discharge. A spatial trend in discharge sensitivity arises at four locations highlighted in red: three sectors in the ablation zone and one sector in the accumulation zone. At these locations, a given increase in normalized strain rate results in higher predicted normalized discharge values than would be expected at other locations along the cable. The dashed line denotes the approximate location of the transition from the ablation zone to the accumulation zone as determined by the drop in correlation of strain rate RMS with wind speed which reflects the cable melting into snow. This point had moved roughly a kilometer up glacier over the course of the experiment and may explain the significant peak in predicted discharge near the transition line.

4.2 Meteo-LSTM and PDD results

The results of the LSTM model run with temperature and precipitation as inputs are shown in Figure 3b. The PDD predictions were shifted according to highest correlation coefficient, corresponding to 5.4 h, before the residual was calculated to account for meltwater transport to the discharge gauge. The MAE of the residuals of the predictions on the test sets of data are 2.29 and 2.27 m³ s⁻¹ for the Meteo-LSTM and the PDD models, respectively.

4.3 Low-frequency versus high-frequency

Models trained on low frequency (<50 Hz) filtered DAS data perform slightly worse with an MAE of 0.68 m³ s⁻¹ compared to 0.67 m³ s⁻¹ of the high-frequency trained models while also having a larger residual standard deviation of 1.25 m³ s⁻¹ compared to 1.18 m³ s⁻¹ of high-pass filtered models. An analysis of 100 LSTM models trained on unfiltered DAS data was also done which performs slightly better than both filtering methods with an MAE and standard deviation of 0.64 and 1.15 m³ s⁻¹, respectively, which may be explained by the broadband nature of the surficial hydrological soundscape (Podolskiy and others, 2023).

5. Discussion

Our study demonstrates the potential for DAS-based glacio-hydrological sensing to be a robust technique for potential in situ measurements of glacier runoff. We find good agreement with 0.64 m³ s⁻¹ MAE between DAS-LSTM-inferred and stream gauge-measured discharge values. We begin this section by discussing why the seismo-acoustic wavefield carries so much correlation with glacier discharge.

5.1 The physical basis relating discharge to the seismo-acoustic wavefield

As described in the Introduction, a wide variety of processes contribute to the glacier seismo-acoustic wavefield. A key result that allows us to decipher the origin of our wavefield–discharge relationship is that our regression analysis performs equally well or slightly better in the range 50–500 Hz as compared to the range 0–50 Hz. This high-frequency band eliminates the possibility that the dominant signal in our analysis has its origin in subglacial processes such as conduit flow (Bartholomaus and others, 2015), gurgling crevasses (Podolskiy, 2020) and bedload transport (Roth and others, 2016, 2017), all of which are thought to create signal below 50 Hz. Furthermore, crevassing and basal stick-slip sliding is expected to generate seismic signals above 50 Hz (Podolskiy and Walter, 2016) in addition to anthropogenic activity and wind (Podolskiy and others, 2023) will also cause increased RMS. However, we infer that the sound generated from supraglacial streams is the dominant contributor to our discharge regression analysis due to its persistent existence during our melt-season measurement. Our basis for this inference is by comparison with previous studies that have examined the same acoustic frequency range in the context of terrestrial rivers (Bolghasi and others, 2017; Osborne and others, 2021, 2022; Podolskiy and others, 2023). Additionally, Figure S3 compares wind from the nearby meteo station to daily means of DAS strain rate and variance RMS observations and shows little correlation throughout the experiment which suggests that supraglacial turbulent flow to be the dominant signal. However, we emphasize that supraglacial turbulent flow is not the sole contributor to the seismo-acoustic wavefield and on glaciers with less prominent supraglacial runoff, it may become more important to disentangle seismic signals generated from other processes from the signal generated from the supraglacial hydrological system.

5.2 DAS offers a stable observation platform on melting glacier surfaces

Ensuring the stability of instrumentation on the surface of glaciers is notoriously challenging (Carmichael, 2019). As a result, most melt season seismic deployments in the ablation zone of glaciers, for example, only cover spatial apertures on the order of 1 km (e.g. Röösli and others, 2014). Studies that have employed dense glacier surface arrays have generally avoided the melt season due to melt-induced tilt and toppling of the instruments (Gimbert and others, 2021). Stream discharge in terrestrial rivers is usually measured by establishing a relationship, called a rating curve, that empirically relates stream height (also called stage) to discharge (Kennedy, 1984). In order to bypass logistical complexities associated with this approach, recent studies have elected to pursue passive acoustic observation of river height (Osborne and others, 2021; Podolskiy and others, 2023). The motivation to use seismo-acoustic observations to study surficial glacier hydrology is even stronger given that seasonal variations in stream morphology (Knighton, 1981; Marston, 1983; Karlstrom and others, 2013) would be expected to result in a strongly time-dependent rating curve. For our study, the deployment of the cable along the glacier flow line allows for sensitivity to source mechanisms in a wide area encompassing both the ablation and accumulation zones. A particular benefit of fiber optic sensing over other methods is that the fiber optic cable can be deployed strategically and is not limited to the specific instrumentation requirements such as the availability of electrical power at the sensing location that hinder many other types of seismic sensing equipment or other in situ instrumentation. In this study, the cable transects many features typical of mountain glaciers: crevasses, supraglacial streams, rock debris, firn, snow, etc.

5.3 Sensitivity analysis

All model iterations show a spatial sensitivity to predicting discharge. In Figure S2, we investigated prediction performance relative to different parts of the cable by isolating the observed acoustic noise in these locations. The DAS data were sectioned in three different ways and used as model input to predict discharge: the whole cable, only channels within the ablation zone and only channels within the accumulation zone. We find model improvement when data within the ablation zone, where we expect the most pervasive surface hydrology to exist, are used for training and prediction. When only the data from the accumulation zone are used for training and prediction, the models perform markedly worse, 1.03 m³ s⁻¹ mean MAE as compared to 0.63  m³ s⁻¹ mean MAE for the ablation data. In addition, the standard deviation of the residuals is three times higher than that of the models using ablation data alone. In the following subsections, we discuss possible mechanisms by which changes in the meltwater flow within supraglacial streams as a result of temporal variation in discharge cause fluctuations in acoustic noise power as observed by DAS.

Figure 4 shows a model sensitivity analysis where we generate a synthetic strain rate Gaussian pulse with a width of 50 channels and uniform in time. The pulse is then centered on each channel before making a discharge prediction. We iterate this procedure for each LSTM model trained on the whole cable DAS data. Increased values of predicted normalized discharge for a given channel in Figure 4 indicate that an increase of measured DAS strain rate or acoustic noise results in an increase of predicted discharge. Three sectors of cable in the ablation zone centered around channel 150, 650 and 1400 are shown to be of more importance to predicting discharge from DAS strain rate. Interestingly, a sector around channel 2250 in the accumulation zone near the glacier headwall also imparts some sensitivity to predicted discharge. The most sensitive portion of the cable is the sector around channel 1400 where the snow line is located during the cable deployment time and the ice fall of Rhonegletscher is located. The melting of snow around the snow line during the observation period caused the snow line to recede and exposed more bare ice to the fiber. Surface crevassing, newly formed meltwater streams and audible drainage from within exposed crevasses may have all contributed to the high RMS strain rate signal in this area. This provides a first step into the potential of forming a spatially distributed, rather than integrated, inference of glacier surface melt.

5.4 Scaling up to other glaciers, longer time spans and the potential for monitoring

We have demonstrated that DAS can be used to infer glacier runoff on Rhonegletscher. The model that we have trained is not expected to be immediately portable to other glaciers, however, for the simple reason that changing the layout of the cable would result in different channel weightings. Glaciers that have differing contributions of runoff and glacier melt to total discharge may require further independent discharge measurements, at least initially, to validate the model inference. Despite these complications, the acoustic noise–discharge relationship does appear to persist with a variety of flow regimes (Podolskiy and others, 2023) which we expect to be the case for supraglacial DAS observations as well.

We have demonstrated discharge inference over a 1-month time period. We do not expect our model to perform well when extrapolated to more diverse glacier surface conditions than those encountered during our deployment. Additional observations would likely be necessary to capture the wide variety of surface energy balance regimes that occur over long time periods. We expect this to be true, for example, when comparing end-member summertime and wintertime conditions (e.g. Chapter. 5, Cuffey and Paterson, 2010). However, it may also be the case that smaller term variations, for example, due to supraglacial stream rearrangement during the summer melt season (Pitcher and Smith, 2019), may also require more detailed training data in order to match the performance attained by our model. Over these long time periods, a number of environmental factors could change and thereby limit the performance of our model. These factors include the coupling of the fiber to the surface of the glacier and, at longer time scales, even the geometry of the cable. During our deployment, the coupling of the fiber in certain areas did change (e.g. melting of snow around the snow line), the effect of which is accounted for in our CV scheme and shows no significant reduction in performance through time (see Fig. 3a). For this reason, we suggest that minimal model retraining would be required on account of changes in fiber coupling alone.

Our data analysis was retrospective, however, there is no fundamental reason why a trained model such as ours could not be used for near-real-time discharge estimation. Once a model is trained on DAS data and tested for accuracy, it can be applied to previously unseen DAS data recorded from the same fiber array to infer discharge. Carrying out prediction with our trained model is orders of magnitude faster than the training and testing steps. Depending on the location of a discharge gauge downstream, this time limitation may be sufficiently faster than a reading from a proglacial discharge gauge and indicates an ability to supplement traditional discharge monitoring networks.

Our approach could provide particularly useful information in several glaciological settings. The most salient example is for glaciers that terminate in the ocean or in lakes. For these glaciers, it is not possible to deploy traditional stream gauges and discharge is generally estimated through melt stakes or surface energy balance calculations (Jackson and others, 2022). In glaciated catchments that have complex networks of proglacial streams or multiple points of runoff, traditional discharge instrumentation would become logistically burdensome.

It is not known at the present time whether it would be possible to train a model on sufficiently many different glaciers or cable layouts so as to arrive at a general model that would capture the relationship between DAS and discharge for an arbitrary new glacier where discharge measurements have yet to be made.

6. Conclusion

In situ measurements of glacier runoff have previously been logistically difficult to obtain, particularly in areas with geographically complicated catchments or glaciers with distributed surface hydrological regimes. We demonstrate a correlation between the in situ seismo-acoustic wavefield measured from the surface of a glacier and proglacial discharge measured by a radar-based gauge. Our ML model that relates these quantities identifies spatial variability and coherence in discharge sensitivity to acoustics. The ability to quantify glacier runoff using turbulent flow-generated seismo-acoustics as observed by DAS opens the door to gaining insights into these regions. Discharge inferences produced by DAS and ML could one day be ingested in glacier mass-balance models that have typically been limited by a lack of in situ glacier runoff validation (Lenaerts and others, 2019). In addition, seasonality of accumulation, ablation and runoff may be characterized by changes in acoustic signals that we observe here during the melt season; however, this will need to be investigated in subsequent studies. Here, we have demonstrated the first of its kind application of DAS for inferring glacier runoff informed by radar-based discharge and other observations.