1. Introduction
In recent decades, it has become possible to deploy active and passive seismic instrumentation in the remote regions of Antarctica and Greenland. Since then, it has been demonstrated that seismic activity in polar glacial environments is largely controlled by strain and stress fields due to the dynamics of the ice sheet, ice caps and large outlet glaciers (Podolskiy and Walter, Reference Podolskiy and Walter2016 and references therein). When huge masses of ice are in motion, the glacier dynamics strongly affect or even control the seismicity within the ice layer and also at the interface with the bedrock (Podolskiy and Walter, Reference Podolskiy and Walter2016) involving a broad range of frequencies of excitation of seismic waves.
It has been observed that episodes of crevassing and ice cracking (Walter and others, Reference Walter, Clinton, Deichmann, Dreger, Minson and Funk2009; Colgan and others, Reference Colgan, Rajaram, Abdalati, McCutchan, Mottram, Moussavi and Grigsby2016; Hudson and others, Reference Hudson, Brisbourne, White, Kendall, Arthern and Smith2020), stick-slip and basal motion (Anandakrishnan and Alley, Reference Anandakrishnan and Alley1994; Danesi and others, Reference Danesi, Bannister and Morelli2007; Creyts and Schoof, Reference Creyts and Schoof2009; Smith and others, Reference Smith, Smith, White, Brisbourne and Pritchard2015; Röösli and others, Reference Röösli, Helmstetter, Walter and Kissling2016), iceberg calving (O’Neel and others, Reference O’Neel, Larsen, Rupert and Hansen2010; Bartholomaus and others, Reference Bartholomaus, Larsen, O’Neel and West2012; Walter and others, Reference Walter, Olivieri and Clinton2013), glacier noise and tremor (Röösli and others, Reference Röösli, Walter, Husen, Andrews, Lüthi, Catania and Kissling2014), hydraulic fractures (Hudson and others, Reference Hudson, Brisbourne, White, Kendall, Arthern and Smith2020), transient subglacial drainage (Winberry and others, Reference Winberry, Anandakrishnan and Alley2009), tidal forcing (Zoet and others, Reference Zoet, Anandakrishnan, Alley, Nyblade and Wiens2012), are some of the main processes capable of raising the tensile and shear stress values until the exceeding of thresholds that trigger seismicity. A wide range of frequencies can be excited in the rupture processes, from icequake episodes (with characteristic frequencies even more than 100 Hz) to tidal modulation (with characteristic frequencies in the order of 10−3 Hz).
The David Glacier is the imposing outlet glacier of the northern Victoria Land region, draining the inner part of the plateau flowing from eastern Dome C, southern Talos Dome and northern Taylor Dome (Frezzotti and others, Reference Frezzotti, Tabacco and Zirizzotti2000). It collects a catchment basin of ∼4% of the East Antarctic Ice Sheet (Rignot, Reference Rignot2002) and feeds the ∼100 km long Drygalski Ice Tongue that floats into the Ross Sea. When the glacier crosses the Transantarctic Mountains, its basal topography has a steep slope change which generates a remarkable icefall (∼400 m downhill) and terminates in the so-called David Cauldron, at the grounding line level (Bindschadler and others, Reference Bindschadler2011, Fig. 1 and Fig. S1 in supplementary material).
Figure 1. a) Map of the seismic networks in operation between 2003 and 2016 around the David Glacier. Triangles and crosses are used for seismic station sites, coloured in agreement with legend. The meteorological station Sofiab and the MZS tide gauge (installed at Mario Zucchelli station) are plotted with yellow diamonds. The map is developed using Quantarctica (Matsuoka and others, Reference Matsuoka, Skoglund and Roth2018) environment while the DEM is the landsat image mosaic of Antarctica (Lima 15 m; Bindschadler and others, Reference Bindschadler2011). Subglacial lakes are extracted from the compilation of Smith and others (Reference Smith, Fricker, Joughin and Tulaczyk2009) and named according to the same nomenclature (D1–D5); the subglacial water flux is from Lebrocq and others, 2013; the grounding (blue) and hydrostatic (green) lines are extracted from ASAID (Bindschadler and others, Reference Bindschadler2011). b) The availability of the waveforms included in the database. DY: Italian network, ZL: New Zealand network. STAR is a semi-permanent Italian seismic station that has been working since 2003.
For the David Glacier area, the ice speed derived from satellite radar interferometry (MEaSUREs collection, Rignot and others, Reference Rignot, Mouginot and Scheuchl2017) varies in a range between a few tens of m/y and ∼500 m/y from the plateau towards the proximity of the grounding line, respectively, with strong acceleration in correspondence of the steepest topographic slopes. For its dimensions, its massive ice flow, and its impact on the dynamics of the front polynya, the David Glacier and Drygalski Ice Tongue have been studied at length since the early 90s (Frezzotti, Reference Frezzotti1993) with important outcomes concerning pioneer and revised estimates of the local mass balance (Frezzotti, Reference Frezzotti1997; Rignot and others, Reference Rignot, Mouginot and Scheuchl2011), the ice thickness, surface speed and flux (Frezzotti and others, Reference Frezzotti, Tabacco and Zirizzotti2000; Wuite and others, Reference Wuite, Jezek, Wu, Farness and Carande2009; Le Brocq and others, Reference Le Brocq, Ross, Griggs, Bingham, Corr, Ferraccioli, Jenkins, Jordan, Payne, Rippin and Siegert2013; Rignot and others, Reference Rignot, Mouginot and Scheuchl2017; Moon and others, Reference Moon, Cho and Hoonyol2021), the morphology of the bottom (Tabacco and others, Reference Tabacco, Bianchi, Chiappini, Zirizzotti and Zuccheretti2000), the grounding line definition (Bindschadler and others, Reference Bindschadler2011; Stutz and others, Reference Stutz2021), the subglacial lake system (Smith and others, Reference Smith, Fricker, Joughin and Tulaczyk2009; Lindzey and others, Reference Lindzey, Beem, Young, Quartini, Blankenship, Lee, Lee, Lee and Lee2020), the stability of the Drygalski Ice Tongue (Indrigo and others, Reference Indrigo, Dow, Greenbaum and Morlighem2021), the glacier thinning and retreat in recent geological eras (Stutz and others, Reference Stutz2021).
Based on the analyses of NASA’s Ice, Cloud and Land Elevation Satellite (ICESat) laser altimeter data, Smith and others (Reference Smith, Fricker, Joughin and Tulaczyk2009) identified the presence of six subglacial, hydrologically linked, lakes in the region of the David Glacier catchment (labelled as D1–D5 in Fig. 1). Although the extent of lakes may be intended to be qualitative and ultimately variable over time due to potential measurement uncertainties (Smith and others, Reference Smith, Fricker, Joughin and Tulaczyk2009), Pattyn (Reference Pattyn2008) recognises that lakes would be responsible for detectable ice mass variations due to transient episodes of drainage and refilling through the subglacial hydrographic network. More recently, Indrigo and others (Reference Indrigo, Dow, Greenbaum and Morlighem2021) have demonstrated the hydrological continuity between the floating ice channelization below the David Cauldron and the upstream drainage subglacial network, which would extend for about 400 km within the David Glacier drainage basin (Le Brocq and others, Reference Le Brocq, Ross, Griggs, Bingham, Corr, Ferraccioli, Jenkins, Jordan, Payne, Rippin and Siegert2013). A similar speculation was inferred by Moon and others (Reference Moon, Cho and Hoonyol2021), who measured the David Glacier velocity changes between 2016 and 2020 by applying the offset tracking technique to Sentinel-1A SAR images. These authors suggest that fluctuations and transient increases in the ice stream flow speed can be attributed to the injection of bottom water through a diffuse supraglacial hydrological system. Miles and others (Reference Miles, Stokes, Jamieson, Jordan, Gudmundsson and Jenkins2022) pointed out that the variability of the David ice discharge between 2005 and 2018 was associated with changes in ice shelf buttressing and the modulating effect of local glacier geometry. The drainage and refilling cycles of subglacial lakes, as well as the discharge of subglacial water, cause basal water pressure and basal sliding to vary, which potentially leads to seismic activity.
Temporary seismological experiments have evidenced the occurrence of low-energy seismicity in the David Glacier area (Bannister and Kennett, Reference Bannister and Kennett2002; Danesi and others, Reference Danesi, Bannister and Morelli2007; Zoet and others, Reference Zoet, Anandakrishnan, Alley, Nyblade and Wiens2012) providing interesting inferences about possible trigger mechanisms. Bannister and Kennett (Reference Bannister and Kennett2002) advanced the hypothesis that a deep tectonic lineament or possible fractures in the ice layer were activated. Subsequently, using data from a dedicated local seismic network deployed around the glacier—partially re-analysed in the present study—Danesi and others (Reference Danesi, Bannister and Morelli2007) suggested that the coupling processes between the ice flow and the bedrock activated the repeating brittle failure of one or more asperities at the ice/bedrock interface with a stick-slip mechanism. Further supporting this, Zoet and others (Reference Zoet, Anandakrishnan, Alley, Nyblade and Wiens2012), analysing data from the remote Transantarctic Mountains Seismic Experiment (TAMSEIS) network, observed a significant correlation between seismic activity and ocean tides, indicating that tidal modulation likely influences the recurrence of these seismic events.
The dynamic process that drives glacier motion and arrest, and the subsequent occurrence of seismicity for basal resistance, is controlled by the competition between sliding velocity—which induces a frictional loading on basal asperities—and the effective normal stress at the ice-rock interface. By studying seismicity at the local scale on an Alpine glacier, Gräff and Walter (Reference Gräff and Walter2021) have demonstrated the key role of the subglacial hydraulic system in the instability of basal asperities.
At larger scales, Zoet and others (2018) have analysed seismicity in the Whillans Ice Stream (WIS, Antarctica), in the light of laboratory experiments, to investigate how subglacial water regelation controls the stick-slip process, in the presence of till. The authors show that a regelation process, and therefore a regelation time, are necessary elements for triggering repeating seismicity, in addition to influencing glacier movement by changing the contact area and deforming the till bed.
A revised analysis of the seismicity in Antarctica from 2000 to 2021 was developed using the Machine Learning approach (Pena Castro and others, Reference Peña Castro, Schmandt, Nakai, Aster and Chaput2024), finding both repeating and non-repeating earthquakes in the area of the David Glacier. To better discriminate between sources in the area, we analyse the evolution of the seismicity observed in the area of David Glacier during three measurement campaigns in a 14-year-long period, between November 2003 and February 2016. We observe that the seismicity shows highly clustered spatial distributions and indicates evidence of a repeated source which, on the other hand, does not reveal characteristics of continuity or seasonality over the long time that could be unequivocally correlated with regular recurrent forcing (temperature cycles, oceanic tides). Based on the spatio-temporal distribution of the earthquakes and their deviation from the magnitude-frequency scaling laws characteristic of tectonic seismicity, we infer that these events are likely induced by glacier dynamics, and we aim to provide a qualitative assessment of potential environmental processes that may have triggered these seismic events. After examining several environmental parameters, we find that clustered and repeated seismic events may be correlated with transient episodes of ice mass discharge correlated to the subglacial hydrographic system (observed by remote sensing analysis), which could lubricate the bedrock interface.
2. Methodology
2.1. Seismic data collection
We have collected all available seismic data registered in the study area during three Italian temporary experiments, funded by the Italian National Program for Research in Antarctica (PNRA, Programma Nazionale di Ricerca in Antartide) and run along 2003–04, 2005–06 and 2015–16 austral summer campaigns (November-January) around the David Glacier. During the first experiment—between November 2003 and February 2004 and jointly conducted by an Italian/New Zealand team (Danesi and others, Reference Danesi, Bannister and Morelli2007)—9 temporary broad-band seismic stations were installed on rock outcrops around the ice stream. The Italian team repeated temporary observations in summer campaigns 2005–06 and 2015–16 by reoccupying only the sites with the best signal-to-noise ratio (green triangles and yellow stars in Fig. 1). For these two last experiments, the seismic network was composed of seven and six stations, respectively.
New Zealand data and metadata are available through IRIS Data Management Center with the ZL network code (stations equipped with a CMG-40 seismic sensor and an Orion broad-band digitizer). All Italian stations (network code DY) were equipped with a Trillium T40 seismic sensor and a Reftek130-01 broad-band digitizer, sampling data at 100 and 125 Hz, and powered by photovoltaic systems.
One of the Italian stations was located at Starr Nunatak (STAR) in November 2003, more than 50 km away from the glacier, and during the years, it has been equipped to become a semi-permanent seismic station, able to record for more than 9 months per year continuously. Due to its stability and performance, the STAR station has so far provided a quasi-continuous database of daily data, including most of the winter seasons.
All the data recovered from the experiments were organized in a SeisComp3 structure (Helmholtz Centre Potsdam GFZ German Research Centre for Geosciences and gempa GmbH 2008), which allowed us to perform a systematic analysis of the whole database; the availability and data volume for each seismic station involved in one or more experiments between 2003 and 2016, marked by network code, is shown in the right panel of Figure 1.
2.2. New 3D seismic velocity model of the glacier based on radio echo sounding data
To build a more accurate 3D glacier structural model we used all the Radio Echo Sounding (RES) datasets recorded during the campaigns conducted by the PNRA program. The datasets cover a period from 1995 to the most recent RES measurement campaign (austral summer 2015–16), and all were focused on the David Glacier area (coloured lines in Fig. 2a) using a low-frequency version of the INGV ‘GlacioRadar’ instrument operating at 12 MHz. This particular configuration was used to integrate the lower number of bedrock reflections available in the Cauldron area in datasets collected at 60 MHz and 150 MHz (Fig. 2b).
Figure 2. PNRA RES data coverage in the David Glacier area. (a) flight tracks used and divided by years of sampling and (b) example of radargram at 12 mhz frequency obtained upstream the cauldron area along the track sketched by black circles on a panel A.
Reflection travel times were converted into depth using a constant electromagnetic (EM) wave velocity of 168 m/μs while cross-point analysis showed a difference in ice thickness of less than 20 m in 86% of cases. Unfortunately, RES data were not specifically designed to perform wet/dry reflection analysis impeding an effective identification of the hydrographic network directly from the radar data.
The 3D bedrock model of the David Glacier has been obtained by subtracting the ice thickness values from the corresponding RAMPDEM2 (Liu and others, Reference Liu, Jezek, Li and Zhao2015) elevation points for all available RES reflection points since 1995. All the elevation data (surface and bedrock) were gridded using the ‘IDW’ (Inverse to Distance Weighting) method. Figure 3 shows the overlay between the ice surface (grey) and the bedrock topography as colour-based relief with isolines at 100 m. Finally, we converted the model into a regular 3D matrix for seismological analyses: we first parameterized a volume of 124 × 100 × 23 km3 over a regular-cell grid (each cell of 0.5 × 0.5 × 0.5 km3) and then we associated ice or bedrock (or both) to each cell, as a result of the difference between surface and bedrock heights in the corresponding point of the ice thickness model. For each cell, we accordingly associated the values of P-wave and S-wave velocity.
Figure 3. The two datasets used to build up the 3D model of the David Glacier area used for the absolute location are represented in this figure. In grey, the topographic surface elevation taken from RAMP2 elevation model (Liu and others, Reference Liu, Jezek, Li and Zhao2015) is overlaid to the bedrock elevation (colour scaled; isolines at 100 m) derived by PNRA RES datasets (vertical exaggeration = 16). All maps are reported in UTM58S projection (WGS84).
For ice layers, we used the compressional-wave velocity v p = 3.8 km/s, and the shear-wave velocity v s = 2.0 km/s (Röthlisberger, 1972). For the bedrock, we used vp = 6.1 km/s from 0 to 9.5 km depth, 6.3 km/s from 10 to 17.5 km, 6.6 km/s for deeper layers, and vp/v s ratio of 1.73 in agreement with the model used by previous work in the area (Danesi and others, Reference Danesi, Bannister and Morelli2007).
Where the grid exceeded the area of the RES survey and direct observations were not available, we merged mean velocity values extracted from a 1D-layered velocity model (Danesi and others, Reference Danesi, Bannister and Morelli2007). Given the geographical extension of the spatial distribution of nodes, the final 3D velocity structure counts 2,352,303 values.
2.3. Seismic data analysis
The seismic waveforms available for the three campaigns were used to characterize the evolution of seismicity around David Glacier with only one of the stations, STAR, located at Starr Nunatak, operating continuously since its installation in 2003 (Fig. 1).
Icequakes (episodes of fracturing of the ice layer) and basal events (events occurring at the ice-bedrock interface) are the most typical signals recorded in glaciated regions (Podolskiy and Walter, Reference Podolskiy and Walter2016).
The seismic network recorded thousands of these events per day, which we classified and distinguished mainly based on frequency content and duration (Supplementary Fig. S2). In particular, (i) the seismic spectrum of icequakes has maximum amplitudes at frequencies higher than 5–10 Hz, while the earthquake spectrum is shifted to lower frequencies; (ii) the high-frequency icequake signal can be recorded only at one or two nearby stations because it decays within 2–5 s, while basal (ice bedrock interface) events can be detected even at stations tens of km away; (iii) icequakes have a signal duration of a few seconds, while basal events generally have longer waveform durations (up to tens of seconds). Mainly due to the very large station spacing (20–100 km), this work focuses on the analysis of basal earthquakes, while all signals classified as icequakes are not considered.
To obtain a catalogue as rich as possible of basal earthquake locations, we followed a four-step scheme: (1) absolute location using a seismic 1D velocity structure; (2) absolute location using the new local 3D velocity structure; (3) relative earthquake relocation with a double-difference approach and definition of clusters of seismicity and (4) seismic catalogue enrichment using a phase-matched filter detection algorithm.
2.3.1. Absolute location—1D velocity model
For austral summer seasons, when the full network data were available, we performed a standard absolute earthquake location by using a simple layered 1D Earth velocity structure as already proposed in Danesi and others (Reference Danesi, Bannister and Morelli2007), which counts one layer of ice (1.5 km deep) overlapping four crustal layers down to 33 km depth, and a half-space. For each station, we performed the standard STA/LTA trigger algorithm (1 and 30 s intervals, respectively) to detect each variation in amplitude exceeding the signal-to-noise ratio (SNR) equal to 3, after a bandpass filter in the range 0.4–4 Hz. We used the seismological software SeisComp3 for manual review of earthquake locations, provided that at least five arrival times (P and S) had been obtained at three stations. Finally, we excluded locations with RMS greater than 0.5 s.
2.3.2. Absolute location—3D velocity model
After the manual location, the events were then relocated following a probabilistic approach with the NonLinLoc software (Lomax and others, Reference Lomax, Virieux, Volant, Berge-Thierry, Thurber and Rabinowitz2000) using the 3D velocity model which was obtained from the GlacioRadar survey (Section 2.2). The 3D-absolute manual location catalogue included 350 weak events most located along the steep slope at the entrance of the David Cauldron, where the glacier topography indicates its maximum declivity and ice flux (Figs 3–5). Black circles in Figure 4 shows the horizontal uncertainty related to each epicentral estimation. In Table S1 and Figure S3 (both in the supplementary material), all of the absolute locations obtained with the 3D model and global statistics of the results are listed and shown.
Figure 4. Map of epicentres (in red) and uncertainties (circles) after relocation with the nonlinloc locator (Lomax and others, Reference Lomax, Virieux, Volant, Berge-Thierry, Thurber and Rabinowitz2000) and the 3D radio-glacial model. See Table S1 of supplementary material for a complete list of the 3D absolute relocated events.
2.3.3. Relative relocation—double-difference approach
Absolute locations obtained by the NonLinLoc inversion were used as input catalogue for a relative double-difference location with the HypoDD algorithm (Waldhauser and Ellsworth, Reference Waldhauser and Ellsworth2000; Waldhauser, Reference Waldhauser2001). We intended to apply a coherent inversion scheme for all data available between 2003 and 2016. After the pre-processing, more than 1700 P arrival times and more than 500 S arrival times were selected. We defined an inversion weighting scheme with a weight equal to 1 for P-phases and 0.5 for S-phases, a residual threshold of 0.5 s, a maximum distance allowed between linked pairs of 1 km, and damping equal to 80. After relocation, median residual times are lower than 0.1 s and median spatial errors lower than 95 m, 110 m, and 195 m, respectively, for the three coordinates x, y and z. After this relative relocation step, we obtained 139 events distributed on 10 clusters of seismicity (Table 1).
Table 1. Geographic coordinates of the centroids of seismic clusters and number of annual occurrences
2.3.4. Phase-matched filter detection analysis and catalogue enrichment
Once the relative earthquake location catalogue was completed, we applied a phase-matched filter detection technique (Chamberlain and others, Reference Chamberlain2017) to retrieve, for each cluster, further possible detections and correlated repeating signals that may have escaped the previous detection search.
We initially considered only the records from the 2003–04 campaign. We selected the seismic trace of the hypoDD master event for each cluster and applied the phase-matched filter technique to find correlated signals recorded at TRIO, the closest station to the David Cauldron (∼40 km away), which provided continuous data for the search time interval and signal-to noise ratio > 5 for the master events. We applied the Python-based EQcorrscan package (Chamberlain and others, Reference Chamberlain2017) on the continuous seismic signals to perform the multi-parallel, matched-filter detection, and we were able to detect all the P and S chunks and extract the waveforms significantly similar (cross-correlation coefficient ≥ 0.7) to the 10 selected master events representing each cluster; we then applied the Median Absolute Deviation (MAD) method to correlate the waveforms, providing more than 1500 cross-correlated seismic signals recorded at the TRIO station for the period November 2003–February 2004.
To estimate the temporal variation of these events, we replicated the analysis for the detection of repeating earthquakes using master events recorded at station STAR, and applying the phase-matched filter technique on the continuous data available over 14 years, extending therefore the time interval to winters and full years when the temporary seismic network was not installed.
3. Results and comparison with environmental parameters
3.1. Seismicity characterization
We achieved a catalogue of repeating events (cross-correlation coefficient between waveforms ≥ 0.7) grouped in 10 clusters.
Figure 5 shows the larger groups of repeating events with sizes of yellow circles proportional to the number of events while the complete number of occurrences inside each cluster during time is listed in Table 1. As expected and in accordance with previous works (Danesi and others, Reference Danesi, Bannister and Morelli2007; Zoet and others, Reference Zoet, Anandakrishnan, Alley, Nyblade and Wiens2012), the centroids of the two largest clusters (Cl_1 and Cl_4, Fig. 5) were located in the David Cauldron at the foot of the icefall (section CC’ in Fig. 5). A significant number of basal earthquakes, collected in the remaining eight clusters, occurred on top of the icefall, mainly after 2005.
Figure 5. Distribution of the clustered seismicity obtained in this work after the relative relocation of the catalog of 3D absolute location. The clusters of repeating earthquakes are represented with orange circles, plotted at their cluster centroid coordinates (see Table 1), and scaled with the cumulative occurrences (when larger than 50) over 14 years. The map is developed using Quantarctica (Matsuoka and others, Reference Matsuoka, Skoglund and Roth2018) environment, the DEM is Lima (15 m; Bindschadler and others, Reference Bindschadler2011), the subglacial water flux is from Lebrocq and others, 2013 and the grounding (blue) and hydrostatic (green) lines are extracted from ASAID (Bindschadler and others, Reference Bindschadler2011). The ice speed map is taken from measures collection (Rignot and others, Reference Rignot, Mouginot and Scheuchl2017). The three transects sketched in the main figure (white dots), depict the topographic variation of the bedrock and the ice thickness along transversal and longitudinal directions of the David Cauldron. A simplified location of the clustered seismicity is represented by the black dots, most of them located at the base of the ice and in correspondence with the bedrock topographic variation.
In Figure 6a, the vertical component of raw seismic signals for the 1588 events, filtered in the frequency band 0.4–4 Hz, recorded by station TRIO between November 2003 and February 2004 with cross-correlation coefficient > 0.95 (i.e. 50 waveforms in Fig. 6b).
The frequency content of most of these signals ranges between 0 and 5 Hz with the most energetic signal reaching 15 Hz (Fig. 6c). The similarity among the waveforms, signal duration, and frequency content suggests that these seismic signals may have a common source.
Figure 6. Examples of cross-correlated events and their frequency distribution. (a) Superimposition of the vertical components of 1588 correlated events, recorded at station TRIO between November 2003 and January 2004. (b) Fifty vertical components of raw seismic signals recorded at station TRIO on day 324 of 2003 (11 Nov 2003), filtered in the frequency band 0.4–4 hz. Signals have cross-correlation coefficients greater than 0.95. (c) A 3-h long waveform for TRIO station and related spectrogram where frequency content of each event is visible.
The distribution of seismicity in space and time (Fig. 5 and Table 1) indicates that David Cauldron hosted the events with the highest similarity threshold (cross-correlation coefficient > 0.8) and a very high number of repetitions (Cl_01 and Cl_04, nearly 1900). For these two event clusters, we adopt the source parameter estimates already obtained by Danesi and others (Reference Danesi, Bannister and Morelli2007) through waveform spectral analysis. In particular, the authors give a mean value for the seismic moment of M0 ∼ 3.8・1011 Nm, and a peak slip velocity of about 60 mm/s (assuming fracture in the ice).
However, the seismic activity in this area had a limited duration (November 2003–December 2003) and ceased definitively after 2004. In the same period, a very large cluster was active on the top of the icefall (Cl_03, 1188 events) and, again, it ceased activity after 2004. Two significant clusters of seismicity were continuously active between 2003 and 2016, although with fewer events (Cl_06 and Cl_10, respectively, 318 and 225 events). Both were located upstream of the icefall, along the main branch of the David ice stream. The remaining clusters (Cl_02, Cl_05, Cl_07 and Cl_09) had few events and were excluded from further interpretation. These results however raised at least two main questions concerning (i) the location of the seismicity in a few specific areas and (ii) the reduction in the rate of seismicity after 2004.
In the weeks following the network installation in 2003, from Julian day 320 to 353, the repeating seismic activity not only shows similar waveforms but also a regular time interval between events, which is about 24.8 ± 4.7 min (Fig. 7a and 7c), in agreement with previous observations (Zoet and others, Reference Zoet, Anandakrishnan, Alley, Nyblade and Wiens2012). Since previous studies have confirmed that the rate of seismicity is highly correlated with the local tidal modulation (Casula and others, Reference Casula, Danesi, Dubbini and Vittuari2007; Zoet and others, Reference Zoet, Anandakrishnan, Alley, Nyblade and Wiens2012), we extracted the available sea level height measurements from the permanent tide gauge observatory located about 70 km north of the Drygalski Ice Tongue, offshore from the Italian base Mario Zucchelli Station (MZS in Fig. 1). The correlation between seismicity and sea level is shown in Figure 7.
Figure 7. Examples of correlation between sea level height and parameters statistics for seismicity. (a) distribution of inter-event time with the sea level height (horizontal axis is time expressed in days after 2003/01/01). In blue MZS sea level height as measured by tide gauge (black vertical axis), in red inter-event time spacing between consecutive events in minutes (red vertical axis). (b) Distribution of the number of events per day (in black) with the mean inter-event time (in red). (c) and (d) Histograms show the distribution of inter-event time spacings for the two periods 320–353 and 354–390 days after 2003/01/01, respectively. The corresponding density functions are superimposed in cyan and the red vertical lines give the inter-event value corresponding to the maximum of density. In general, the density values d(xi) satisfy the following relation ∑id(xi)(bi +1−bi) = 1, where (bi +1−bi) is the interval between bins.
The evolution of seismicity exhibits an abrupt change after day 354 (11 Dec, 2003) when the number of events per day significantly drops (Fig. 7b) and the inter-event time doubles up to ∼50 min (Fig. 7d). Seismicity located at the top of the David Glacier icefall persisted over the following years; conversely, further repeating occurrences have not been observed in the David Cauldron area after 2004, suggesting that the activation of the downstream events was likely triggered by a transient driver.
The decrease in the number of events per day after Julian day 354 may be due to the effect of the katabatic wind, which would be responsible for the increase in seismic background noise and the subsequent decrease in detection efficiency, especially when the wind speed exceeds 12 m/s (6 knots), as observed by Frankinet and others (Reference Frankinet, Lecocq and Camelbeeck2021).
For the 2003–04 campaign, therefore, we extracted the wind speed measurements for the SofiaB weather station, located 35 km upstream from the TRIO station (Fig. 1), from the open-access database ‘Antarctic Meteo-Climatological Observatory at MZS and Victoria Land’ (http://www.climantartide.it).Comparing the wind speed measurements recorded by SofiaB with our seismic signal at TRIO station (Fig. 8), we confirm the effect of the wind speed on the detection efficiency (the number of events is lower when the wind speed is high). On the other hand, wind gusts exceeding 10 m/s do not appear to significantly impact the root mean squares (RMS) values in the frequency band 0.1–4 Hz, which is the characteristic frequency band of the events under study. The number of seismic occurrences in the second half of the 2003–2004 summer campaign is consistently lower than the number of detections in the first half, regardless of the wind speed registered (Fig. 8).
Figure 8. Effect of the wind on the detectability of the events. From top to bottom: (a) daily noise at TRIO station in the two main frequencies of seismic events (0.1–4 hz) and wind (5–15 hz); (b) the number of picks (in green) and its 3-days running mean (in red) obtained by the triggering detection of station TRIO; (c) mean hourly wind speed recovered from meteorological station Sofiab.
3.2. Wavelet cross-correlation with tide height, wind speed and meteorological parameters
The correlations between the occurrence of seismicity and some environmental parameters, such as tide height and wind speed, were verified also using the cross-wavelet transform approach which enables the examination of signal coherence over time and the identification of any common characteristic frequencies. The WaveletComp 1.1 software, integrated with the R package (Rösch and Schmidbauer, Reference Rösch and Schmidbauer2018), was employed for this purpose. Specifically, the inter-event time spacings of the 2003–04 seismic clusters were compared to the tide height (Fig. 9a).
Figure 9. (a) Cross-wavelet power spectrum between the sea level heights and the inter-event time spacing of the seismic events as a function of time (horizontal axis is in days after 2003/01/01). (b–e) cross-wavelet power spectrum between meteorological parameters and the inter-event time spacing of the seismic occurrences. In all plots, the coloured scale indicates the cross-wavelet power level at each period, that is, a dimensional because the data sets have been normalised: the yellow colour indicates low correlation, while the blue colour indicates high correlation. Black arrows indicate the phase shift between the two time-series: when arrows point to the right the time-series are in phase, when arrows point to the left they are in counter phase. Note the logarithmic vertical scale. The white lines delineate the statistically significant areas, at 10% significance level against a white noise null.
As the tide heights data have one sample per hour whereas the inter-event time spacings are irregularly distributed over time, we chose to compute the upper and lower envelope of the tide time series, corresponding to the combination of spring and neap, according to Figure 7. These envelopes were used to obtain the tide heights corresponding to the epoch of each event of the seismic clusters by interpolation. Interestingly, the resulting cross-wavelet power spectrum (Fig. 9a) implies that the inter-event time during the 5-week clustered seismicity is not correlated with the tide period, because the values are close to zero (the cross-wavelet power spectrum can be interpreted as the local covariance between two time-series). Conversely, a strong correlation between the two time-series can be observed just after day 355 (12 Dec 2003), with a period of 14 days and the smallest phase difference (Fig. 9a). This result suggests that the tidal modulation is the most probable forcing of the seismicity after day 355, possibly controlling the clusters located between the grounding and the floating lines, while the source of the clustered activity before that date should be attributed to a different cause, limited in time and superimposed over the tidal forcing itself.
Furthermore, the main meteorological parameters recorded at the SofiaB weather station (atmospheric pressure, air temperature, relative humidity, wind speed) for the summer campaign of 2003–04 were extracted and the wavelet cross-correlations between inter-event time spacing and each meteorological parameter were calculated (Fig. 9). A notable correlation between these meteorological parameters can be observed between Julian days 005 and 020 (January 2004), with a characteristic period of approximately 3 days. However, no correlation is evident with the seismic clusters recorded at the end of 2003.
Finally, in the presence of large crevasses along the steepest stretches upstream of the icefall, surface meltwater can infiltrate to the bed of the glacier, lubricating the bedrock interface and changing the seismic signals. This process of hydrofracturing has been studied both theoretically and observationally for its potential effects on Antarctic ice shelf weakening (Lai and others, Reference Lai2020). However, since the surface temperature recorded at the SofiaB weather station never exceeded −5°C between November 2003 and January 2004, with mean value −16°C (Fig. S4 in Supplementary Materials), we tend to rule out hydrofracturing as a driving process of instability in this area. Additionally, seismic events resulting from surface fractures in polar regions typically exhibit high-frequency characteristics, generally exceeding 10 Hz (Lombardi and others, Reference Lombardi, Gorodetskaya, Barruol and Camelbeeck2019). These frequencies are associated with near-surface brittle deformation processes, such as crevasse formation, and can vary with environmental conditions and ice dynamics. Given that the observed seismic signals fall outside this typical frequency range, it is unlikely that surface fractures are their source.
3.3. Ice mass variation from GRACE measurements
The Gravity Recovery And Climate Experiment (GRACE; 2002-2017) and its Follow-On mission (GRACE-FO; 2018-still working) (Dahle and Murboeck, Reference Dahle and Murboeck2019) give important information related to the Earth’s gravity field changes due to mass variations. For the entire hydrological Antarctic Ice Sheet_315 basin (AIS_315) that includes the David Glacier and the area under investigation, an increase in the mass is recorded from mid-November 2003 to mid-February 2004 (red box Fig. 10a), followed by a rapid decrease up to April 2004. The Gravity Information Service (GravIS) of the German Research Centre for Geosciences (GFZ) makes the products derived from the satellite missions GRACE and GRACE-FO (Sasgen and others, Reference Sasgen, Groh and Horwath2019, Reference Sasgen, Groh and Horwath2020) available. In particular, the mass changes of the Antarctic Ice Sheet are provided in terms of a time series of gridded ice-mass changes per surface area with a spatial resolution of 50 × 50 km. The David Glacier region was identified in these grids as a small portion (approximately 7 pixels, yellow box Fig. 10b) within the AIS_315 basin. These pixels span a length of 350 km, extending from the glacier catchment area towards the coast. GravIS database made available the surface mass density and boundary grids for the basin, which allowed the calculation of the corresponding variations in glacial mass pixel by pixel, with a resolution higher than the average value on the entire basin (AIS_315) which can be directly downloaded from the GravIS website. A noteworthy variation of ice mass, up to 0.4 Gt (yellow colours in Fig. 10b), was recorded between November and December 2003, in correspondence with the 5 weeks under investigation, when the highly correlated seismicity was active. Pixels close to the grounding line of the David glacier have registered an increase in mass, whereas the backward pixels are characterized by mass reduction.
Figure 10. Ice-mass variation in the David Glacier area observed by GRACE. (a) ice mass variation inside the whole AIS_315 basin from March 2003 and October 2004. Blue and orange lines are the values by GFZ and COST-G (international combination service for time-variable gravity field; Meyer and others, Reference Meyer2020) estimates. The red box indicates the period of interest of the seismicity, when both the models give an ice mass increment in the basin. (b) map of the AIS_315 basin with the overlay of the 50 × 350 km2 area (yellow box), where the local ice mass discharge is calculated monthly (lower rectangles). The position of 31dpt is marked with a black cross. The ice mass variation is given in GTON scale and is calculated for each cell and each year between 2003 and 2016 as recorded by GRACE and GRACE-FO.
The mass variation observed by GRACE could be due to snow accumulation. However, Frezzotti and others (Reference Frezzotti, Urbini, Proposito, Scarchilli and Gandolfi2007) observed that during 2003 the snow accumulation at 31Dpt (74 S, 156E; Fig. 10) was higher than 2002 and 2004 and similar to average between 1966 and 1998 and during the period 1998–2004. This observation allows us to reject the hypothesis that the mass variation is due to change in snow accumulation.
Smith and others (Reference Smith, Fricker, Joughin and Tulaczyk2009) proposed the occurrence of a transient ice mass reduction of approximately 1.1 km3 between November 2003 and March 2008, potentially linked to the discharge of water from the subglacial lakes D2 and D3 and the subsequent refilling of D1 (Fig. 1). The maximum volume displacement rate was concentrated in the initial six months of the period under consideration.
Despite the paucity of direct observations of changes in the glacier basal conditions over the 5 weeks under study, both the GRACE measurements and the evidence provided by Smith and others (Reference Smith, Fricker, Joughin and Tulaczyk2009) lead us to believe that a significant circulation of water may have been injected into the subglacial drainage system below the David Glacier during that period.
The clustered seismic signals that we found in this work show variable duration between 4 and 30 s, and frequency content concentrated in the 1–10 Hz band.
Due to the geometry and interdistances of the seismic network, the possibility that the observed low frequencies are caused by the combined effect of surface wave propagation in the surface ice layer -thereby acting as a wave guide- and attenuation cannot be excluded. However, from a physical point of view, this frequency content is compatible with the resonance of subglacial fluids and hydraulic transients (Podolskiy and Walter, Reference Podolskiy and Walter2016).
This analysis therefore suggests that transient injection of subglacial fluids from upstream can reasonably be considered as the main trigger for the seismic events observed at the floating zone level.
However, the seismicity in the area is most affected by tidal modulation (especially after day 355), being periodically decoupled from the bedrock and rich in water-saturated sediments (till). The periodic injection of the water stream could therefore favor the reduction of the basal shear stress and the acceleration of the glacier.
4. Discussion
We have analysed the evolution of seismicity over a non-continuous period of 14 years around the David Glacier; it should be noted that the geometry of the network, especially the large inter-station distances (up to 100 km), together with the assumptions concerning the crustal propagation model, generate substantial systematic uncertainties in the determination of the depth of hypocentres. Errors in depth determination necessarily affect the estimation of magnitude. Previous studies (Danesi and others, Reference Danesi, Bannister and Morelli2007) inferred the magnitude Mw for the 2003–04 correlated events using a standard Brune approach and found that the magnitude values followed a delta-like distribution peaking around the value 1.3. In this paper, together with the absolute manual location, the Ml values were calculated from the waveform amplitudes and compared with the above-mentioned previous results (Fig. S5 of supplementary material). We obtain varying Ml in the range − 1 ≤ Ml ≤ 4 with the significant limitations described above.
Nevertheless, despite these limitations and given the favourable azimuthal distribution of the stations, the epicentral location of the seismicity is well resolved. Based on the frequency content of the waveforms, we infer that the origin of the seismicity is at the bottom of the ice layer, at the interface with the bedrock, including the till where present (section AA’, BB’ and CC’ of Fig. 5).
Most of the seismicity is gathered in two clusters of low-energy basal earthquakes (Cl_01 and Cl_04 in Fig. 5 and Table 1) characterized by highly correlated waveforms (cross-correlation coefficient ≥ 0.95, Fig. 6), concentrated in a space of ∽2 km2 and in a time-range of 5 weeks between November and December 2003. The high similarity of the waveforms suggests a common source, supporting the recurrent activation of an asperity at the ice-rock interface. No further seismic activity of a similar nature was recorded in the following months or years, up to the last available 2015–16 austral summer campaign. However, significant seismicity has occurred continuously over the 14 years along the largest ice stream that feeds the David Cauldron and channels one of the main branches of the hydrographic network (Cl_03, Cl_06, Cl_9 and Cl_10 in Fig. 5).
We attribute the seismicity of clusters Cl_03, Cl_06, Cl_09 and Cl_10 to elastic loading due to sliding velocity over local asperities. In fact, this is the section of David Glacier where the ice stream changes its slope and direction (sections in Fig. 5) upstream of the icefall. In the following, however, we will focus the discussion on the origin of Cl_01 and Cl_04. The latter events are located in the area of the glacier grounding line. This seismicity may be triggered by the injection of ocean water at the interface between ice and bedrock. The rising ocean tide may force the decoupling at the base of the glacier and cause both fracturing and stick-slip motion (Lucas and others, Reference Lucas2023). However, given the characteristics of the recorded seismic activity, including the number of occurrences, inter-event spacing time, and location, we suggest that the events occurring between November and mid-December 2003 may have originated from a different source. If it is reasonable that the tidal modulation represents the primary force acting throughout the entire observation period (see Figs 7 and 9), the origin of clusters Cl_01 and Cl_04 is more likely to be found as a transient, local, non-seasonal forcing superimposed on the tide. The effect of weather conditions on seismicity can be excluded, at least in November 2003 (Figs 8 and 9), as there is no evidence that they affect the occurrence of seismicity during the three austral campaigns taken into account (Fig. S6 in Supplementary Material). Instead, the frequency content and the duration seem compatible with a hydraulic transient.
For the period 2003–04 under study, the GRACE measurements of ice mass variation show a trend compatible with the injection of a sizable amount of water into the hydrographic network beneath the David Glacier. The associated discharge would account for the increment in basal lubrication, the possible glacier acceleration, and the reduction of the basal shear stress which could favour the stick-slip seismic mechanism. In fact, studies around the coupling between glaciers and the rigid and deformable beds suggest that the dominant processes by which ice advances and moves past asperities on the bedrock are regelation and plastic flow until the ice speed velocity reaches a threshold depending on environmental factors such as ice rheology and roughness (Hooke, Reference Hooke2005; Zoet and others, Reference Zoet2020) which activates the seismicity. Zoet and others (Reference Zoet2020) have estimated that basal stick-slip seismicity can be generally triggered when the sliding speed is between 500 and 2000 m/year, which is compatible with the David Glacier flow velocity (Vittuari and others, Reference Vittuari, Zanutta, Lugli, Martelli and Dubbini2023; Moon and others, Reference Moon, Cho and Hoonyol2021 and references therein).
However, aero-geophysical surveys of the David Glacier catchment area (Lindzey and others, Reference Lindzey, Beem, Young, Quartini, Blankenship, Lee, Lee, Lee and Lee2020) suggest the revision of the subglacial lake locations and even cast doubts about the presence of real lakes, rather preferring the interpretation of a distributed system of subglacial drainage. On the other hand, significant fluctuations in the David Glacier flow velocity, characterized by sudden transient increases in the ice speed (about 5–10%) with no regularity over time, were observed by Moon and others (Reference Moon, Cho and Hoonyol2021) and Miles and others (Reference Miles, Stokes, Jamieson, Jordan, Gudmundsson and Jenkins2022). Moon and others (Reference Moon, Cho and Hoonyol2021) pointed out that the observed glacier velocities during the Antarctic summer increased for at least three years after 2016, probably due to the extension of the summer melt in the ice surface and thus the increase in basal sliding. In addition, Miles and others (Reference Miles, Stokes, Jamieson, Jordan, Gudmundsson and Jenkins2022) found that the change in ice discharge could be driven by both external forcing (e.g. katabatic wind) and internal ice-sheet processes. The area of velocity increase observed by Moon and others (Reference Moon, Cho and Hoonyol2021) is upstream of the David Cauldron and occurs at an elevation of more than 600 m. This area does not show any significant melting, as indicated by AWS (Automated Weather Station) data (e.g. Fig. S4), surface, and satellite image surveys. Therefore, our interpretation is that the glacier velocity increase observed by Moon and others (Reference Moon, Cho and Hoonyol2021) should be more likely attributed to bottom sliding rather than induced by surface melting. Unfortunately, the frequency of satellite image acquisition during the 2003–04 period was insufficient to permit the surveying of any surface velocity change.
5. Conclusions
Monthly GRACE data from the austral summer of 2003 show a mass variation along the catchment of David Glacier towards the coast, which can be related to the emptying and refilling of the regional subglacial hydrographic network and the subsequent acceleration of glacier flow. The basal lubrication conditions correlate with seismic events that we have identified in relation to the main flow towards the David Cauldron. Notably, there is a decrease in the number of seismic events detected during the 2005–06 and 2015–16 field campaigns. This could be due to a reduced number of recording stations, or it could indicate a reduced inflow of lubricating water through the subglacial hydraulic network.
Recent Synthetic Aperture Radar measurements (Moon and others, Reference Moon, Cho and Hoonyol2021) confirm that the dynamics of David Glacier are influenced by transient fluctuations in flow velocity. Although the seismic and satellite observations were not made simultaneously, both seem to indicate recurrent (at least seasonal) behaviour in the glacier. This suggests a correlation between seismic events and periodic injection of liquid water into the subglacial hydraulic system.
Regardless of the origin or nature of the subglacial water transfer, our study shows that seismic activity appears to be triggered by basal lubrication, as suggested by Zoet and others (Reference Zoet2020). Large-scale observations emphasize that the dynamics of outlet glaciers are highly sensitive to significant changes in the hydrological system (Maier and others, Reference Maier, Andersen, Mouginot, Gimbert and Gagliardini2023), which can affect their velocity and geometry over time. Therefore, our results suggest that continuous monitoring of seismic activity at the ice-bedrock interface could be a valuable tool for tracking ice sheet dynamics in the context of global warning.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2025.10062.
Data availability statement
All raw data can be provided by the corresponding authors upon request.
Acknowledgements
We are very grateful to both reviewers for their fruitful and encouraging comments. Maps were produced in the QGIS Geographic Information System environment, using the free database Qantarctica provided by the Norwegian Polar Institute. We used R and Rstudio IDE (Posit Team 2025) for data analysis and representation. We thank DISTART University of Bologna (IT) and DIMEC University of Modena e Reggio (IT) for providing local tide gauge data. Meteorological parameters were obtained from ‘Meteo-Climatological Observatory at MZS and Victoria Land’ of PNRA - http://www.climantartide.it. The Gravity Information Service (GravIS, we thank Ingo Sansgen) of the German Research Centre for Geosciences (GFZ) provided observational data derived by the satellite missions.
Author contributions
S.D. designed the experiments, carried them out, analysed the seismic data, and conceived the manuscript. S.S. carried out the experiments and analysed seismic data. A.B. analysed non-seismic time series and GRACE data. S.U. and A.Z. performed the glacio-radar survey and analysed radar data. M.F. reviewed all available data and led the building of the conceptual model. S.D. prepared the manuscript with contributions from all co-authors.
Funding statement
This work was supported by the Progetto Nazionale di Ricerca in Antartide (PNRA 2010/A2.09 2013/AZ2.09) funded by the Italian Ministry of Research (MUR), and by the MACMAP Project funded by Istituto Nazionale di Geofisica e Vulcanologia (Environment Department).
Competing interests
The authors declare that they have no conflict of interest.
Open access
