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.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⁻¹, 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⁻¹ 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⁻¹ (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⁻¹, 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⁻¹). 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⁻¹, respectively. The 4.74 m d⁻¹ velocity was interpreted as best fit to BH2 data after speed-up (blue shading in Figure 4c), while the 2.62 m d⁻¹ 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⁻¹ 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⁻¹.

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.