1. Introduction
Glacier velocity and its response to both internal and external changes is a crucial parameter for understanding ice fluxes and mass balance (van der Veen, Reference van der Veen2013; Benn and others, Reference Benn, Fowler, Hewitt and Sevestre2019). Traditionally, satellite-based remote sensing methods, such as Interferometric Synthetic Aperture Radar (InSAR), speckle tracking of SAR imagery and optical feature tracking (e.g. Murray and others, Reference Murray, James, Macheret, Lavrentiev, Glazovsky and Sykes2012; Strozzi and others, Reference Strozzi, Paul, Wiesmann, Schellenberger and Kääb2017) or speckle/offset tracking (e.g. Gardner and others, Reference Gardner2018; Lei and others, Reference Lei, Gardner and Agram2021) are used to monitor glacier movement. Nevertheless, the temporal resolution, especially at higher latitudes is limited to at least a few days and, in the case of optical data, depends on cloud cover and visibility, often affecting the spatio-temporal data quality. Furthermore, there are velocity changes that occur on hourly and daily timescales that may be missed by the feature-tracking estimates (Schmid and others, Reference Schmid, Radić, Tedstone, Lea, Brough and Hermann2023; Zwally and others, Reference Zwally, Abdalati, Herring, Larson and Steffen2002).
Since 2014, the launch of Sentinel-1 and Sentinel-2 as part of the Copernicus program has enabled the continuous acquisition of glacier velocities over 6–12 day intervals in near real-time over large areas (e.g. Lei and others, Reference Lei, Gardner and Agram2021; Koch and others, Reference Koch, Seehaus, Friedl and Braun2023). Among typically employed methods, InSAR performs best over short timescales when surface decorrelation is small, while feature-tracking techniques can track the movement of features from days to weeks (Heid and Kääb, Reference Heid and Kääb2012; Sund and others, Reference Sund, Lauknes and Eiken2014). At well-established monitoring sites, detailed measurements of glacier velocity may be performed with high precision using differential GNSS (dGNSS) stations (e.g. Winberry and others, Reference Winberry, Anandakrishnan, Alley, Bindschadler and King2009; Pickell and Hawley, Reference Pickell and Hawley2024). However, running long-term observations requires maintenance visits and is more prone to data gaps (e.g. Richter, Reference Richter2013; Błaszczyk, Reference Błaszczyk2024).
Seismic stations (seismometers, accelerometers, geophones) record ground shaking with very high-temporal resolution, ranging most often from 100 or above samples per second, and are widely used for glacier monitoring (Podolskiy and Walter, Reference Podolskiy and Walter2016; Aster and Winberry, Reference Aster and Winberry2017). To guarantee accurate sampling and synchronization among all the independently operating stations, their digitisers are equipped with an internal GPS receiver for timing calibration (Ringler, Reference Ringler2021). GPS signals are acquired continuously, or intermittently, to reduce power consumption, and measurements’ spatial precision does not need to be high. Except for time calibration, the GPS data are usually discarded as other techniques, like dGNSS, are often used in conjunction to record its position which in most of the typical seismological applications is stable over time. This is, however, not true for an increasing number of cryoseismological studies (Podolskiy and Walter, Reference Podolskiy and Walter2016), where seismic stations are deployed in glaciated regions (e.g. Köhler and others, Reference Köhler, Maupin, Nuth and Van Pelt2019; Sergeant, Reference Sergeant2020; Hudson and others, Reference Hudson, Brisbourne, Kufner, Kendall and Smith2023; Nanni and others, Reference Nanni, Roux and Gimbert2024).
Here, we propose that the internal GPS data from seismic stations can be used to monitor glacier velocity despite its limited precision. This can be particularly advantageous when noother locational data are available, or the existing time-series have gaps and/or insufficient temporal resolution. We focus the analysis on a case study at Borebreen, Svalbard, where seismic stations were deployed to monitor glacial seismicity connected to a current glacier surge. By comparing the velocities inferred from the internal GPS data with satellite-derived time series, we conclude that these systems can provide reliable data for glacier dynamics at no extra cost as a by-product of their primary function.
2. Materials and methods
2.1. Study area and instruments
Borebreen, located in the Svalbard archipelago, is a marine-terminating glacier that began a surge in 2021 and continues at the time of writing, making it a suitable location for testing our hypothesis. We employed seismic digitisers DATA-CUBE3 from DiGOS manufacturer, with a GPS antenna mounted inside the unit. Three stations were installed in the central part of Borebreen (Fig. 1) for 98 days between 10 April and 16 July 2024.
Figure 1. Orthophoto map of the borebreen study area. The dots indicate the positions of seismic stations. The colour of the dots represents the total displacement in the horizontal direction measured by the internal GPS over the time of deployment. The black rectangles mark the area contributing to the Sentinel-1 velocity time series. The baseline is a Planet Labs SkySat photo from 11 September 2023. The inset map in the top right corner presents the location of Svalbard and the location of Borebreen in Svalbard (red dot).
2.2. Seismic station internal GPS data
Seismic stations need to include GPS systems for internal clock synchronization. These systems, while not intended for precise geodetic measurements, provide positional data that can potentially be repurposed for tracking glacier movement. For this study, the internal GPS data from three seismic stations, installed in the ablation zone of Borebreen, were extracted and processed to derive a time series of glacier velocities at each location. Some seismic instrument manufacturers provide an option to extract GPS positions from the raw data files. We extracted all timestamps with their latitude, longitude and elevation from raw seismic digitiser data files. In our experiment, GPS synchronization was performed in cycled mode for 3 minutes within 30-minute intervals for a reason of power saving yielding about 360 position and time samples per hour (however, the exact number of samples varies). Upon extraction, a minimalistic quality control (QC) procedure is applied. We filter out all measurements where the number of satellites used for positioning is smaller than five and then remove all measurements with latitude and longitude values falling outside of one sigma limits computed over 12 hour-long intervals (Fig. 2). Subsequently, for each interval, mean and standard deviation of latitude, longitude, elevation and number of satellites is computed, which effectively resamples the data to 12 hours. The displacement is calculated as the horizontal distance from the first averaged sample and the velocity as a gradient of displacement over time. The resulting position (with associated uncertainty) of the top station BH3 over time is presented in Fig. 3.
Figure 2. Distribution of Northing and Easting (UTM) coordinates in the first 12-hour-long analysis interval at the BH3 (top) site. (b) A scatter plot showing all (black dots) and filtered data (red dots) with a black error bar representing the mean and standard deviation recalculated for the filtered (red) dataset. Data points are plotted with 90% transparency; (a) and (c) Northing and Easting coordinates, respectively, distributions presented as histograms with grey bars for all the data and red bars for filtered data. The red curve represents the probability density function (PDF) inferred for all the data (scale is not presented) and dashed blue lines represent the limits of inferred pdf’s one, two, and three sigma intervals of the whole dataset.
Figure 3. Position of seismic station BH3 on borebreen over time mapped as (a) Easting and northing in UTM coordinates; (b) total horizontal displacement and elevation change. the cross centres mark the mean coordinates values in each 12-hour-long window, the cross’s arms’ size shows standard deviation (note different XY scales in b), and the colour shows the time elapsed since 10 April 2024. As a background in (a) a satellite photo from 11 September 2023 is used with visible surface crevasses across the flowline. Two insets a and b in (a) show zoomed-in 20 × 20 m sections of the areas marked with black rectangles.
2.3. Satellite-based glacier velocity
We compared the seismic GPS data to satellite-derived ice velocity from Sentinel-1 SAR data using standard feature-tracking methods (e.g. Strozzi and others, Reference Strozzi, Luckman, Murray, Wegmuller and Werner2002). This method detects the displacement of intensity features between consecutive images and, where phase coherence is preserved (e.g. when temperatures are below freezing), can also detect the displacement of speckle features. All 12 day repeat-pass image pairs from Sentinel-1A were utilised, the image patch size was 416 × 128 pixels (approximately 1 km in ground range), the velocity was sampled every 50 × 10 pixels (∼100 m in ground range), and the resulting velocity map was geocoded to EPSG:32633 using the TanDEM-X DEM. Finally, we used the resulting Sentinel-1 velocity time series to infer the representativeness of seismic GPS data for ice flow velocity mapping.
Alongside SAR glacier velocities from Sentinel-1, we compare the seismic station GPS data against the velocity data from NASA MEaSUREs Intermission Time Series of Land Ice Velocity and Elevation (ITS_LIVE) (Lei and others, Reference Lei, Gardner and Agram2021). We acquired a time series of ITS_LIVE velocity at the coordinates matching the initial positions of three seismic stations (Fig. 1). The time series is primarily composed of glacier velocities derived from optical data such as Sentinel-2, Landsat 8 and 9, while Sentinel-1 is unavailable for the studied period. The ITS_LIVE data processing methodology is detailed by Lei and others (Reference Lei, Gardner and Agram2022).
3. Results and discussion
3.1. Performance of seismic internal GPS
Surface sliding velocity data from seismic station BH2 internal GPS were compared against satellite-derived time series for the same area (Fig. 4a). For direct comparison, the GPS data, averaged in 24-hour windows, were downsampled to match the Sentinel-1 data by selecting only the measurements available on common dates. Both methods successfully captured the major temporal changes in glacier motion. Notably, the internal GPS at BH2 showed good agreement with Sentinel-1-derived velocities with a strong Pearson’s correlation coefficient of 0.93, indicating the validity of this approach.
Figure 4. Comparison of satellite and seismic station internal GPS time series. (a) Direct velocity comparison at BH2 area: red dots represent sentinel-1 data, black dots represent downsampled GPS data, the grey curve represents daily averaged BH2 GPS data, and green dots represent ITS_LIVE data; (b–d) cumulative displacements along the flowline at the sites BH1, BH2, and BH3, respectively: red dots represent sentinel-1 data, black dots show GPS data in 12 hour averaged windows with uncertainty, and green dashed line represents ITS_LIVE data. Grey solid and dashed lines represent displacements from sliding velocities indicated at respective panels, that were fitted to satellite data (red shading) and GPS data (grey shading), respectively. The difference between sentinel- and GPS-derived displacement is shown with blue crosses and between ITS_LIVE and GPS with blue diamonds (right-hand side axis).
From 20 April to 30 May, when the surface sliding velocity at the BH2 site was approximately 2.3 m d−1, the horizontal positioning uncertainty increased to around three m across all three stations (Fig. 5b), exceeding the daily displacement. This substantial uncertainty was likely caused by a layer of fresh or redeposited snow covering the installation sites, which significantly degraded the GPS signal-to-noise ratio. A marked improvement in positioning precision to ∼ 0.5 m occurred as air temperatures rose from 30 May onward (Fig. 5b), melting the attenuating snow layer. Over this period, the ratio of uncertainty to daily displacement decreased notably to 0.2. To provide context, the mean ratio of standard deviation to velocity in the ITS_LIVE dataset is 0.1. Improved site design could help mitigate or minimise such issues in future installations. While the limited precision of seismic station internal GPS is a drawback, the method remains viable for long-term experiments or fast-flowing glaciers with flow rates significantly exceeding the measurement uncertainty.
Figure 5. Seismic station GPS horizontal positioning precision along the flowline. (a) Histogram of the qc-ed standard deviation values in each 12 hour window for all three stations; (b) The qc-ed standard deviation over time marked with coloured scatter points. The colour of each point corresponds to the number of samples used during the averaging. The daily temperature at Borebreen measured by each BH1–BH3 station and averaged is marked by blue–red line. It may be affected by direct sunlight. Additionally, the average number of satellites used for position determination in each 12 hour window is plotted with black dots.
The elevation change initially appeared too noisy for reliable interpretation. However, following a drop in positioning uncertainty, the elevation change became clearer, declining approximately 15 m over the distance of 150 m within 40 days (Fig. 3b). This value falls within the expected range for surface inclination in the studied area.
The velocity domain is unsuitable for analysing location data with high-temporal sampling, as moderate spatial uncertainty is greatly amplified during gradient computation. Therefore, we use denser, 12 h averaged GPS data in the displacement domain. Fig. 4b–d presents comparisons of total horizontal displacement and the difference between the seismic station GPS and satellite data, which reaches a maximum of 15 m after 3 months (comparable to Sentinel-1 data pixel size of 10 m). Changes in surge dynamics are discussed in the following section.
While seismic internal GPS outperformed Sentinel-1 in terms of temporal resolution, it still has significant limitations in precision when compared to dedicated geodetic systems, e.g. low-cost cryosphere-dedicated GNSS stations with low power consumption characterised by Pickell and Hawley (Reference Pickell and Hawley2024), which can achieve precision levels of ∼ 0.01 m. In our study, after processing, the uncertainty of seismic internal GPS horizontal positioning ranged from 0.4 to 4 m (Fig. 5a). The post-processed seismic GPS data has similar temporal resolution and precision as ITS_LIVE data (Fig. 4a). However, during the polar night the ITS_LIVE catalogue must rely solely on Sentinel-1 derived velocities, with availability varying across regions and periods. Notably, at the time of writing, the Sentinel-1 data in the ITS_LIVE catalogue for Svalbard is only available up to 2022, meaning it is not accessible for the studied period. Hence, during poor visibility conditions and polar night, the ITS_LIVE time series contains gaps, potentially missing short-term variations in sliding velocity. Furthermore, the feature-tracking process inherently smooths the 2D ice velocity field by using feature-matching windows that overlap in space. This leads to the matching of similar features and therefore more gradual changes in the ITS_LIVE velocity field (Lei and others, Reference Lei, Gardner and Agram2021) compared to Sentinel-1 and GPS data. This difference is clearly visible as consistently observed fluctuations in residual values for all sites (Fig. 4b–d). Therefore, we propose that the method described here can serve as a promising alternative for obtaining higher temporal sampling than SAR Sentinel-1 and more complete temporal coverage than ITS_LIVE, especially at sites where installing more advanced sensors is not feasible due to risks of losing equipment, limited funding, or power supply constraints. Additionally, it can be used to fill observational gaps in time series from other systems when necessary, e.g. during polar night.
3.2. Surge dynamics inferred from seismic internal GPS
A significant spring speed-up of Borebreen, occurring in June, is evident in all analysed datasets (Fig. 4), with a two-fold increase in surface velocity in all three sites between 15 June and 1 July (from 2.27 to ∼ 4.7 m d−1 at BH2 site). A more detailed pattern emerges from the higher temporal resolution seismic GPS data when analysing cumulative horizontal displacements detrended by the initial sliding velocity of 2.27 m d−1 (Fig. 6). The increase in air temperature, followed by a rise in sliding velocity, suggests that meltwater reaches the glacier bed, increasing water pressure and enhancing ice flow. This effect shows a lag time of about one week at all sites. Similar velocity increases have been documented in other tidewater glaciers in Svalbard, e.g. Kronebreen and Stonebreen (Benn and others, Reference Benn, Hewitt and Luckman2022). However, Borebreen was already surging, and the observed velocity increase implies that basal energy or enthalpy levels were high. Consequently, even a small rise in air temperature could trigger a significant sliding speed-up. This aligns with the enthalpy model of Benn and others (Reference Benn, Fowler, Hewitt and Sevestre2019), which suggests that enthalpy gains from supraglacial sources can exceed enthalpy losses, thereby enhancing sliding velocities.
Figure 6. (a) Cumulative displacement detrended by a linear trend of 2.27 m d−1, as derived from satellite observations at BH2 (red shading in Figure 4c) (i.e. negative slope informs about sliding slower than 2.27 m d−1). sentinel-1 data are shown as large dots. GPS data are represented by small dots with uncertainty, and ITS_LIVE data are plotted with lines. Blue, orange and green colours indicate BH1, BH2 and BH3 sites, respectively. The opaque grid of solid and dashed grey lines indicate displacements corresponding to sliding velocities of 2.62 and 4.74 m d−1, respectively. The 4.74 m d−1 velocity was interpreted as best fit to BH2 data after speed-up (blue shading in Figure 4c), while the 2.62 m d−1 velocity was fitted to BH1 GPS data prior to speed-up in the interval marked by vertical blue dashed lines. (b) Temperature at Borebreen (blue–red line) and precipitation at Longyearbyen (black bars). The temperature was measured by each BH1–BH3 station and averaged. It may be affected by direct sunlight.
Importantly, the higher temporal resolution of the detrended displacement allows for easy identification of stable periods and changes in the sliding velocity revealing dynamics of acceleration with duration in the order of one week. Such short-term velocity variations cannot be resolved by Sentinel-1 data, which, per the Nyquist–Shannon theorem, only captures variations longer than 20 days with 10-day sampling, as for Borebreen (Shannon, Reference Shannon1949). Crucially for investigating glacier dynamics, the enhanced resolution allows more precise identification of the speed-up initiation, at the BH2 site, to the first days of June, while differences in the timing of speed-up initiation among the three stations suggest that sliding velocities propagated upstream from the terminus (Fig. 6a). Moreover, the high-resolution GPS data, which are not subject to the inherent smoothing present in feature-tracking methods, allow to infer an intermediate ∼ 2.6 m d−1 mode during the speed-up phase and determine approximate one-week duration of switching between intermediate to fast sliding velocities, which, at all three locations, eventually stabilised at ∼ 4.74 m d−1.
Additionally, while the acceleration at BH2 may appear coincidental with precipitation events, a lower-magnitude speed-up was recorded at BH1 since 20 April, coinciding with air temperatures nearing 0°C. This suggests that meltwater is the primary driver of the observed speed-up and supports the theory that high water pressures and subglacial enthalpy are maintained (Thøgersen and others, Reference Thøgersen, Gilbert, Schuler and Malthe-Sørenssen2019), enabling even minor inputs of surface meltwater to the bed to accelerate glacier flow. These findings highlight the significant influence of seasonal ice flow changes on surge velocities.
4. Conclusions
This study highlights the potential of seismic stations with internal GPS to monitor glacier velocity. Results from Borebreen, Svalbard, demonstrate a strong correlation between velocities derived from internal GPS data and those from satellite-based methods, documenting that these systems may be a valuable extra resource for glacier monitoring. Despite its lower precision (up to 0.5 m) compared to dedicated geodetic systems and potential temporal instabilities, quality-controlled and temporarily averaged data can effectively capture glacier movement at much higher temporal resolution than the Sentinel-1 data, which had a sampling of 2 and 10 days interchangeably. The ability to pinpoint the timing and duration of speed-up phases, as well as to identify intermediate velocity modes, underscores the importance of temporal resolution for capturing short-term variations in sliding velocity. The presented method could be particularly useful in remote study sites, especially on fast-flowing glaciers, where maintaining a large number of geodetic stations is logistically challenging due to difficult terrain and a lack of power sources. Given the increasingly more frequent deployment of seismic stations in polar regions, this approach could enhance existing glacier observation networks with relatively high-temporal resolution of position tracking data at no extra cost.
We provide insights into the dynamics of surging glaciers, emphasizing the role of meltwater as a key driver of speed-up events. The observed upstream propagation of sliding velocities and the stabilization at ∼ 4.7 m d−1 highlight the temporal and small-scale spatial variability of surge development. Identifying intermediate velocity phases and the ∼ 20 consecutive positive-degree days required to transition into steady fast flow underline the importance of thermal and hydrological conditions in modulating surge behaviour.
Our findings highlight the critical role of meltwater in driving seasonal ice flow dynamics during glacier surges. The observed upstream propagating spring speed-up, with velocities increasing from 2.27 to ∼ 4.7 m d−1, underscores the sensitivity of basal conditions to temperature variations. These findings align with the enthalpy model, where even minor supraglacial meltwater inputs can enhance basal enthalpy and trigger significant flow acceleration. Additionally, the lower-magnitude speed-up recorded earlier in the season suggests that sustained high water pressures and enthalpy are key to maintaining surge dynamics. This study demonstrates how seasonal meltwater variations influence glacier behaviour and highlights the usefulness of auxiliary high-temporal resolution observations in understanding surge mechanisms.
Data availability statement
The original data presented in the study and the Python codes used for data processing and figure generation are openly available at https://doi.org/10.5281/zenodo.14936324. The precipitation data were downloaded from Meteostat (https://meteostat.net/en, Svalbard Lufthavn). For extracting position data from DiGOS Data-CUBES, we used ‘gipptools’ (version 2024.170) provided by the manufacturer. The ITS_LIVE dataset was downloaded from https://itslive-dashboard.labs.nsidc.org/.
Acknowledgements
We would like to acknowledge Heidi Sevestre and Gabrielle Kleber for fieldwork assistance, Eirik Hellerud, Piotr Kupiszewski and Einar Berntson from the Norwegian Polar Institute (NPI) and CHC pilots for their help in planning and executing the fieldwork. Also, we acknowledge Shad O’Neel, Hester Jiskoot and Frank Pattyn for constructive comments which substantially improved this paper.
Author contributions
Conceptualization, W.G.; methodology, W.G., A.L.; software, W.G., A.L.; validation, W.G., A.L.; formal analysis, W.G., W.D.H.; investigation, W.D.H., W.G., D.M.P. and R.C.; resources, W.G., A.L., W.D.H.; data curation, W.G, A.L.; writing—original draft preparation, W.G.; writing—review and editing, All; visualization, W.G, project administration, W.D.H.; funding acquisition, W.D.H, R.C. All authors have read and agreed to the published version of the manuscript.
Funding statement
We acknowledge funding from the Royal Geographical Society (RGS) Walters Kundert Fellowship, a Norwegian Research Council (RCN) Arctic Field Grant, and from the Svalbard Integrated Arctic Earth Observing System (SIOS) via an access to Planet Imagery grant.
Open access
