An application of a seismic nodal system with seismic ambient noise near Kunlun Station, Antarctica: estimating ice thickness and firn structure

1.   Introduction

The Antarctica ice sheet, covering over 98% of this continent, stands as the largest ice mass on Earth. Its thickness emerges as a pivotal parameter in the modeling of ice sheet dynamics and contributions to sea level changes (Budd, 1991; Alley et al., 2005; Hanna et al., 2013). The upper structure of the ice sheet, generally within 150 m depth of the snow surface, is equally vital. It is characterized by a depth-increasing density profile that results from the densification of snow into glacial ice, in which the firn layer is defined as the intermediate stage between fresh snow and underlying ice (Alley, 1987; Cuffey and Paterson, 2010). This depth-density profile is governed mainly by temperature and snow accumulation rates, which are known to vary significantly across the climatically diverse Antarctic continent (van den Broeke et al., 2008). Detailed knowledge of the snow and firn layer is critical for conducting accurate ice mass balance studies and paleoclimate reconstruction (Parrenin et al., 2012; Shepherd et al., 2012).

Traditionally, the investigation of ice sheet thickness has been approached by radio-echo sounding (i.e., ground-penetrating radar investigation) alongside some contributions from other methods, such as active seismic techniques (Drewry et al., 1982; Plewes and Hubbard, 2001; Bingham and Siegert, 2007; Phạm and Tkalčić, 2018). The densification profile of the upper ice sheet can be directly measured by ice core drilling and borehole logging, or indirectly using seismic refraction or radar measurements (Godio and Rege, 2015; Morris et al., 2017; Schlegel et al., 2019; Hollmann et al., 2021).

Rapid developments of seismic ambient noise methods have made it practical to investigate ice sheet thickness and subsurface structures by using ambient noise recorded by seismic sensors, at lower logistical cost and greater deployment versatility compared to traditional methods (Shapiro et al., 2005; Podolskiy and Walter, 2016; Picotti et al., 2017). Specifically, the horizontal-to-vertical spectral ratio (H/V) method and the ambient noise cross-correlation method (seismic interferometry) have been applied in different glacier environments. In the H/V method, the peak frequency of the H/V curve is related to S-wave velocity and thickness of the subsurface medium (Lundin et al., 2017). Lévêque et al. (2010) used the H/V method in Dome C, Antarctica, to estimate the thickness of the snow layer. Picotti et al. (2017) validated the reliability of the H/V method in detecting glacier thickness at different glaciers. Yan P et al. (2018) estimated ice sheet thickness across Antarctica by applying the H/V method to single station ambient noise records. Berbellini et al. (2019), and Jones et al. (2021) obtained the seismic velocity profile of ice sheet and uppermost crust in polar regions using Rayleigh wave ellipticity (a branch of the H/V method). For the ambient noise cross-correlation method, surface wave dispersion curves can be extracted from noise cross-correlation functions, and then inverted for S-wave velocity profile estimations (Xia JH et al., 2003). Zhan ZW et al. (2014) used noise cross-correlations and auto-correlations to estimate water column thickness under the Amery Ice Shelf. Walter et al. (2015) retrieved Greenland Ice Sheet phase velocity dispersion measurements for thickness estimation. Diez et al. (2016) inverted the seismic structure of the Ross Ice Shelf by using dispersion curves obtained from ambient noise cross-correlations. Chaput et al. (2022), and Zhang ZD et al. (2022) deployed small dense seismic arrays in Antarctica to study the upper structure of the ice sheet using the ambient noise cross-correlation method. Sergeant et al. (2020) investigated Rayleigh waves by employing various noise sources in glacier sites in Greenland and Apls. Zhou W et al. (2022) performed distributed acoustic sensing (DAS) measurements to investigate the snow and firn layers on the Rutford Ice Stream, Antarctica. Moreover, Fu L et al. (2022, 2023) deployed seismic arrays near the Dalk Glacier, East Antarctica and then used the ambient noise H/V method and cross-correlation method to investigate glacier structures and crustal structures, but their results from the two methods were interpreted individually.

Joint inversion of the H/V curve and dispersion curve from ambient noise cross-correlation has been proven to be a powerful and relative simple approach that helps to resolve trade-off issues within the individual methods (Piña-Flores et al., 2017, 2020; Sivaram et al., 2018; García-Jerez et al., 2019). The dispersion curves are sensitive to the absolute average S-wave velocity; however, they suffer from a non-uniqueness problem due to the wide range of possible depths (Spica et al., 2018; Petrescu et al., 2023). In contrast, the H/V curves are most sensitive to S-wave impedance contrast (i.e., density × velocity) with the trade-off between depth and average velocity of the sedimentary layers (Picozzi et al., 2005). Combining these two complementary methods in a joint inversion scheme produces results that have been shown to be more reliable (Giancarlo, 2010; Piña-Flores et al., 2017).

To verify the usefulness and accuracy of our joint inversion scheme in glacier studies, it was important to have independent investigation results for comparison. Thus, we chose Kunlun Station, Dome A, located at the highest summit of the Antarctic ice sheet, as our study site. The ice sheet at Kunlun Station is generally flat and stable; an ice divide separates the opposing directions of ice-flows passing through Kunlun Station. Its location makes Kunlun Station an international research hotspot for searching the oldest ice (Severinghaus et al., 2010) and understanding the ice sheet system (Sun B et al., 2009). A variety of studies have been carried out near Kunlun Station, such as ice core drilling and ground-penetrating radar investigations (Cui XB et al., 2010, 2016; Jiang S et al., 2012; Yang WX et al., 2021). A number of investigations regarding subglacial topography and ice flow velocities have also been made along the inland traverse route from the coastal Zhongshan Station to Kunlun Station since 2016 by techniques such as ice cores and snow pits (Xiao CD et al., 2008; Zhang SK et al., 2008; Ding MH et al., 2015; Tang XY et al., 2015, 2016; Yan P et al., 2018). All these previous investigations provide independent constraints to verify our results, estimated from ambient noise methods.

To obtain the seismic structure of the ice sheet near Kunlun Station, Dome A, Antarctica, and to enhance our understanding of the usefulness of a seismic nodal system in glacier environments, we deployed a dense linear seismic array and a single seismic node to record seismic ambient noise during China’s 39th expedition to Antarctica. Below, we first present the study area and describe the seismic ambient noise experiment in the field. Then we detail the results of the H/V method from a single seismic node and the ambient noise cross-correlation method from the seismic array. We invert the ice sheet thickness using the lower frequency portion of the H/V curve. We then apply joint inversion of the dispersion curve from noise cross-correlations and the higher frequency portion of the H/V curve to describe the upper densification structure of the ice sheet. Our results are found to be consistent with the BedMachine ice thickness dataset (Morlighem et al., 2020; Morlighem, 2022) and with a density profile from radar investigations, which in turn has been validated by an ice core (Yang WX et al., 2021).

2.   Study Area and Data

Kunlun Station is located at the inland Dome A central area (Figure 1a); its surface elevation of more than 4000 m makes it the highest region of the Antarctica ice sheet. Horizontal stability at such an ice divide makes this an excellent place for ice studies; the mean surface velocity of the ice sheet over the entire Dome A area is ~8.8 cm a⁻¹; the measured velocity at Kunlun Station is even lower, at ~2.9 cm a⁻¹ (Yang YD et al., 2018). The Gamburtsev Subglacial Mountains range, the largest subglacial mountains in Antarctica, lies beneath the ice sheet in this region. Thus, multiple studies have been performed there. A 109 m deep ice core drilling at Dome A was conducted to reconstruct the history of volcanism (Jiang S et al., 2012). Ground-penetrating radar investigations of the topography of the Gamburtsev Subglacial Mountains have been carried out around Kunlun Station (Cui XB et al., 2010, 2016; Yang WX et al., 2021).

Figure  1.  (a) Location of Kunlun Station in Antarctica. (b) A Smartsolo seismic node installed on the surface of ice sheet. (c) BedMachine (Morlighem, 2022) ice thickness map showing the location of the linear dense seismic array (white circles) deployed near Kunlun Station (yellow star). The seismic node where the H/V method is applied is represented by the blue diamond. The inset illustrates the nearby area around Kunlun Station. The coordinates are projected in EPSG:3031 Antarctic Polar Stereographic projection and the units are in meters.

DownLoad: Full-Size Img PowerPoint

During China’s 39th expedition to Antarctica, we conducted an ambient noise seismic experiment at close proximity to Kunlun Station. Between 5 and 10 January 2023, 20 Smartsolo three-component 5 Hz seismic nodes were deployed in a dense linear pattern near Kunlun Station (Figure 1b).

The first seismic node was located at the far end from Kunlun Station and the distance between each node was approximately 500 m, which resulted in a total length of 11.1 km. In addition, one seismic node was set up for H/V measurement on the array line. All seismic nodes were installed on the surface of the ice sheet with central spikes (Figure 1c). The sampling frequency was 250 Hz for all nodes. The continuously recorded seismic data were divided into one-day segments, and the preprocessing steps included instrument response and trend removal. Unfortunately, the horizontal-component of the seismic array data has proved to be unusable due to possible instrument installation issues (Figure S1).

3.   Methods

3.1   H/V Method

The ambient noise horizontal-to-vertical spectral ratio (H/V) method computes the ratio between the average of the horizontal and vertical Fourier spectrum. This method, popularized by Nakamura (1989), is widely used to assess the predominant frequency of a site when high impedance contrast (i.e., density × velocity) exists (Bard, 1999) and to detect sediment thickness (Sylvette et al., 2006). Over the years, methods have been developed for modeling the H/V curve to be used in geophysical investigations. Reviews of these methods have been made by Lunedei and Malischewsky (2015), and Molnar et al. (2022).

In this study, we employed a more recently developed H/V forward calculation and inversion method based on the diffuse field assumption (DFA) and the contribution of body and surface waves (Sánchez-Sesma et al., 2011; García-Jerez et al., 2016). Under DFA, the H/V curve can be expressed in terms of the imaginary parts of the Fourier-transformed Green’s functions for coinciding source and receiver. This relationship can be expressed as:

The imaginary component of the Green’s function, $I\left[{G}_{ii}\left(x,x,\omega \right)\right]$, is the response at a point in the $i$ direction resulting from a unit force applied at that same location. Here, the indices 1 and 2 denote the horizontal orientations, while 3 indicates the vertical orientation, and ω stands for circular frequency. The Green's function components are associated with a specific geometry and certain material properties in a horizontally layered structure. Consequently, Equation (1) computes the H/V curve theoretically and allows inversion. In the inversion process, model spaces are set for P- and S-wave velocities, density, Poisson ratio, and thickness of each layer. Monte Carlo sampling inversion method was applied to find the optimal parameters in the model space according to the misfit function (García-Jerez et al., 2016), which is defined as:

where ${{H}{V}}^{\mathrm{o}\mathrm{b}\mathrm{s}}$ and ${{H}{V}}^{\mathrm{t}\mathrm{h}}$ refer to the observed H/V curve and the H/V curve obtained through forward calculation, m represents a model, subscript $j$ denotes the $j$th point of the curve, and $\sigma$ denotes the standard deviation associated with the observed H/V curve, ${HV}^{\mathrm{o}\mathrm{b}\mathrm{s}}$ .

We employed the python package “hvsrpy” (Vantassel, 2021) to compute the H/V curve in the 0.1–10 Hz frequency band from ambient noise recorded by one seismic node near Kunlun Station. During the relatively short observation period, we selected one continuous hour of seismic data that were recorded under calm conditions and without transient signals. The H/V curve was calculated using 60-second time windows without overlapping. Fourier amplitude spectra were smoothed through a Konno and Ohmachi (1998) window with a coefficient of 40. The horizontal east (E) and north (N) components were averaged by their geometric mean when divided by the vertical-component (Albarello and Lunedei, 2013). We used a frequency-based window-rejection algorithm that implemented in “hvsrpy” to increase the quality of the H/V curve (Cox et al., 2020).

After obtaining the H/V curve, the 0.1–2 Hz frequency band of the curve was used to invert the ice sheet thickness. The H/V inversion procedure was conducted in the widely used “HV-inv” (v2.5) software (García-Jerez et al., 2016). To speed up the inversion procedure, the observed H/V curve was downsampled to 50 points. Two models were used to perform ice sheet thickness inversion; we followed the model space set up by Yan P et al. (2018) with BedMachine ice thickness as a reference (as listed in Table 1). Model A is a simple model: homogeneous ice overlying the bedrock. Model B is a two-layer ice sheet structure: a low S-wave-velocity bottom-ice layer is between the surface-ice layer and the bedrock. After the inversion, two optimal solutions were best fitted to the observed H/V curve.

Table  1.  Parameter space for model A and model B for ice sheet thickness inversion. The range of Poisson ratio is 0.1–0.4 for all layers.

Model	Layer	${V}_{{\mathrm{p}}}$ (m s−1)	${V}_{{\mathrm{s}}}$ (m s−1)	$\rho$ (kg m−3)	Thickness (m)
A	Ice	3800–4000	1800–2000	917	1820–3380
	Rock	4300	2480	2500	Half space
B	Ice I	3750–4000	1800–2000	917	1560–1950
	Ice II	3500–3750	1400–1600	917	650–1040
	Rock	4300	2480	2500	Half space

 | Show Table

DownLoad: CSV

3.2   Ambient Noise Cross-correlation Method

Ambient Noise Cross-correlation Functions (NCFs) between pairs of seismic stations can be used to generate the surface wave dispersion curve and to obtain the corresponding shear-wave velocity structure (Shapiro et al., 2005). This method relies on the cross-correlation of ambient noise recorded at seismic arrays to approximate the Green’s function between station pairs (Shapiro and Campillo, 2004). We employed the “CC-FJ” (Li ZB et al., 2021) software to calculate the vertical-component cross-correlation functions from the linear seismic node array. We followed the cross-correlation processing steps of Bensen et al. (2007). Specifically, we downsampled noise records to 50 Hz, applied a 1–10 Hz band-pass filter, and adopted spectral whitening during cross-correlation.

We applied dispersion analysis to the NCFs using the phase-weighted slant-stacking technique (Cheng F et al., 2021), which offers higher spectral resolution than the widely used phase shift method (Park et al., 1998). From the NCFs derived from vertical-component seismic data, the dispersion curve of the fundamental mode Rayleigh wave can be manually picked from the peak values of the spectra while avoiding artifacts in the dispersion spectra.

3.3   Joint Inversion of H/V and Dispersion Curves

To determine the upper structure of the ice sheet, the dispersion curve and the higher frequency band (i.e., 2–9 Hz) H/V curve were processed together by a joint inversion scheme. Non-uniqueness of the inversion problem can be reduced by joint inversion of these two measurements and by exploiting their complementary sensitivity to the structure parameters (Piña-Flores et al., 2017; Spica et al., 2018). The H/V curve carries information associated with S-wave impedance contrasts of a layering structure, but is weakly sensitive to its absolute velocity. The dispersion curve, on the other hand, is sensitive to the absolute velocity value, but only weakly sensitive to the details of layering (Spica et al., 2018). The “HV-inv” software and the Monte Carlo sampling inversion method were used in this final stage in order to characterize the seismic velocity structure. Under DFA, the synthetic dispersion curve can be computed at the same time as the synthetic H/V curve (García-Jerez et al., 2016; Piña-Flores et al., 2017). A joint misfit is calculated as the sum of the misfits from each data type:

where $C$ represents a dispersion curve and $k$ denotes the $k$th point of the dispersion curve. Note that in this study we assigned equal weights to the contributions of H/V curves and dispersion curves; we also assigned an error of 100 m s⁻¹ to the picked dispersion curve when calculating the joint misfit.

In our study, a three-layer model was set up to represent the seismic structures of the upper ice sheet. The model parameter space used for inversion is listed in Table 2; the bottom half-space of the model was set to glacier ice. In joint inversion, the range for seismic properties was wide, and the thicknesses for layers were referenced from a firn densification model (Ligtenberg et al., 2011). In the inversion process, we normalized the amplitude of the synthetic H/V curves by dividing the whole frequency band by three to fit the peak frequency of the observed H/V curve, even if the amplitudes were less matched (details are elaborated in Section 4 below).

Table  2.  Parameter space for joint inversion for upper structure of the ice sheet. Poisson ratio is 0.1–0.4 for all layers.

Layer	${V}_{{\mathrm{p}}}$ (m s−1)	${V}_{{\mathrm{s}}}$ (m s−1)	$\rho$ (kg m−3)	Thickness (m)
Snow	600–4000	200–2000	200–917	10–30
Firn	600–3750	200–2000	200–917	30–120
Ice	3900	1900	917	Half space

 | Show Table

DownLoad: CSV

4.   Results and Discussion

The ambient noise H/V curve was obtained from data collected by a seismic node very close to Kunlun Station (Figure 1b). The node was often disturbed by human and vehicle activities; therefore, for our H/V curve computation we selected a one-hour long seismic record collected during a calm period. In the practice of similar H/V measurements, a one-hour long record is considered adequate to obtain a stable H/V curve (Yan P et al., 2018; Zhang YH et al., 2023). The seismic nodal system adopted in this study has been demonstrated to be able to record seismic phases down to 0.1 Hz, well beyond its natural frequency, thus empowering it in a wide range of seismic applications (Zeckra et al., 2022). Available in the supplementary materials (Figure S2–S4) are an example of the ambient seismic record, improvement achieved after applying the frequency-based window-rejection algorithm, and a comparison of H/V curves calculated from two different record lengths.

The H/V curve we constructed from the seismic node data shows three clear peaks in the 0.1–10 Hz frequency band (Figure 2). Specifically, there is a peak at ~0.17 Hz and two other peaks at ~3 Hz and ~6 Hz. The possibility that the peak at ~6 Hz is the higher-mode of the ~3 Hz peak can be ruled out; it is not three times the frequency of the ~3 Hz peak, as would be expected from the theory pointed out by Carcione et al. (2017). Our preliminary interpretation of the H/V peaks, consistent with existing H/V theory, is that the peak frequencies are associated with different layer interfaces (Picotti et al., 2017; Yan P et al., 2018; Molnar et al., 2022). We estimate the corresponding interface depths using the relationship equation $f={V}_{\mathrm{s}}/4h$, where $f$ is the peak frequency, ${V}_{\mathrm{s}}$ is the S-wave velocity (1900 m s⁻¹), and $h$ is the thickness of the medium layer of interest (ice, here). This simple and quick estimation reveals that the major peak at ~0.17 Hz corresponds to a 2794 m layer, suggesting that this peak represents the impedance contrast between ice and bedrock. Meanwhile the two peaks in the higher frequency bands correspond to features of the upper structure of the ice sheet.

Figure  2.  (a) H/V curve calculated using seismic ambient noise data recorded by the H/V seismic node. The solid line represents the average H/V curve; the dashed line represents the standard deviation. (b) Zoomed-in area of the H/V curve showing the two high-frequency peaks.

DownLoad: Full-Size Img PowerPoint

Figure 3 presents optimal S-wave velocity models, A and B, derived from the H/V inversion in the 0.1–2 Hz frequency band. The synthetic H/V curves of the optimal results for model A and model B both fit the observed H/V curve at the peak frequency (i.e., 0.17 Hz, Figure 3a). Model B estimates the bottom-ice layer thickness to be 759 m, and the total ice thickness to be 2577 m, which is consistent with the BedMachine ice thickness (i.e., 2600 ± 40 m at Kunlun Station). The model A estimation of total ice thickness (2480 m), however, deviates 120 m from the BedMachine reference thickness (Figure 3b). The two-layer ice sheet structure (model B) is therefore in better agreement with the BedMachine data, and is consistent with the results given by Yan P et al. (2018), which computed the ice sheet thickness across Antarctica using the H/V method. Their findings support the argument originally stated in Wittlinger and Farra (2012), which suggests that the Antarctic ice sheet is stratified and consists of a bottom-ice layer that makes up around 1/3 of the entire ice column and exhibits low S-wave velocity. The S-wave velocity drop in the bottom-ice layer may be caused by the preferred orientation of ice crystals and fine layering of soft and hard ice layers (Wittlinger and Farra, 2012).

Figure  3.  (a) The synthetic H/V curves and the observed H/V curve in the 0.1–1 Hz frequency band. The synthetic H/V curves are modelled using the optimal inverted S-wave velocity profiles for both model A and model B. (b) The optimal inversion S-wave velocity profiles for model A and model B. The horizontal dashed line indicates the reference BedMachine ice thickness; shaded area shows the uncertainty (Morlighem, 2022).

DownLoad: Full-Size Img PowerPoint

The subglacial topography and ice thickness in the central 30 km × 30 km region at Dome A have been studied in ground-penetrating investigations (Cui XB et al., 2010, 2016), which conclude that the subglacial topography under Kunlun Station is a U-shaped glaciated main trough valley and tributaries, and that the bedrock elevation is undulating with a minimal ice thickness of 1618 m and a maximal ice thickness of 3139 m (Cui XB et al., 2010). It should be noted that ice thickness obtained by H/V inversion reflects the average thickness of an area (spatial resolution) beneath the seismic station; the radius of this area is $\lambda /4$, where $\lambda$ is the S-wave seismic wavelength (determined by peak frequency and S-wave velocity of ice) (Picotti et al., 2017; Yan P et al., 2018). In this study, the spatial resolution is around 5.6 km. Therefore, the ice thickness result of H/V inversion coincides with the interpolated nature of BedMachine data (500 m resolution); however, BedMachine data may have neglected small-scale subglacial topography, a difference that should be taken into consideration when comparing H/V inversion results to model ice sheet dynamics where the bedrock undulates intensely (Cui XB et al., 2016).

The NCFs obtained by correlating vertical-component ambient noise of seismic node pairs with clear Rayleigh waves retrieved from the NCFs are shown in Figure 4a (The envelope of the NCFs is shown in Figure S5). We note that the NCFs were asymmetric and that the Rayleigh wave signal on the acausal part (i.e., negative times) of the NCFs was much stronger than that on the causal part (positive times). Such asymmetry suggests that the ambient noise sources were not uniformly distributed azimuthally (Sabra et al., 2005; Zeng XF et al., 2017; Sergeant et al., 2020). Given that anthropic activities concentrated around Kunlun Station (Figure 1b) are the mostly likely sources of noise, the acausal part of the NCFs is probably clear energy, a conclusion supported by the fact that the noise is concentrated in the high-frequency (5–25 Hz) band (Figure S6), which is associated with anthropogenic origin and has been utilized in clarifying sources in seismic ambient noise studies in Antarctica (Chaput et al., 2022). Although the non-uniform distribution of noise sources may bias the retrieved Rayleigh waves, the effects of directional noise sources are found to be small (Yang YJ and Ritzwoller, 2008; Zeng XF et al., 2017; Pearce et al., 2024).

Figure  4.  (a) Noise cross-correlation functions (NCFs) filtered in the 1–10 Hz band using the vertical-component noise records. The dashed red lines denote moveouts with propagating velocities of 400 m s⁻¹ and 1700 m s⁻¹, respectively. (b) Dispersion spectra extracted from NCFs by use of the phase-weighted slant-stacking technique. The black circles and solid line represent the picked fundamental mode dispersion curve.

DownLoad: Full-Size Img PowerPoint

We calculated the dispersion spectra using a phase-weighted slant-stacking technique (Cheng F et al., 2021) (Figure 4b). The extracted fundamental mode dispersion curve falls within 3–10 Hz. At 4 Hz, the velocity indicated by our dispersion curve is ~1700 m s⁻¹; at a higher frequency (10 Hz), the velocity is ~900 m s⁻¹. In addition to the fundamental mode, artifacts of different kinds can be seen in the dispersion spectra. These artifacts are probably spatial aliasing artifacts related to the linear array geometry, rather than higher modes, with which they might be incorrectly identified (Cheng F et al., 2023). However, in this study, which utilized a small-length seismic array with short interstation distances and a short observation period to constrain the near-surface velocity structure (Spica et al., 2018), those artifacts do not interfere significantly with the fundamental mode dispersion curve. The effective investigation depth determined from the dispersion curve can be estimated using the empirical relationship between depth and the wavelength $\lambda$ of the fundamental mode Rayleigh wave (i.e., the depth equals $\lambda /3\lambda /2$). The wavelength $\lambda$ can be estimated from the velocity associated with the minimal frequency (${f}_{\mathrm{m}\mathrm{i}\mathrm{n}}$) of the dispersion curve as ${\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}}={C}_{R}\left({f}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)/{f}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ (Park et al., 1999; García-Jerez et al., 2019). Under this assumption, we estimate that the maximum depths are between 141 m and 212 m, which is deep enough to characterize the upper structure of the ice sheet related to firn densification.

The dispersion curve extracted from NCFs and the 2–9 Hz frequency band H/V curve were processed together in a joint inversion scheme under DFA. We observed double peaks in the 2–9 Hz frequency band H/V curve (Figure 2). Such a double peak phenomenon hints at the presence of two impedance contrast interfaces, that is, at a three-layer structure (Bard, 1999; Molnar et al., 2022). Therefore, we used a three-layer model in the inversion. Additionally, we noticed that the double peaks of the H/V curve show a relatively lower peak amplitude (i.e., ~1.2) compared to the peak (i.e., ~3) at ~0.17 Hz. Despite their significantly lower amplitudes, the double peaks are distinctly discernible, and the standard deviation of the H/V curve is minimal. In a similar ice sheet environment (at Dome C, Concordia Station) Lévêque et al. (2010) reported H/V peaks with an amplitude range of 1–2, occurring within the 1–10 Hz frequency band, corresponding to a 23 m snow layer. The interpretation of H/V amplitudes actually remains a subject of ongoing debate, as the H/V amplitude can vary depending on many factors, such as the energy content of the recorded signal, the damping of the material, or even the coupling conditions between soil and the sensor (Albarello and Lunedei, 2011; Bignardi et al., 2016; Zhang YH et al., 2023). However, the main focus of the H/V method in this study was to characterize the depth of impedance contrast layers by constraining the peak frequencies, because they can be robustly determined even when the amplitudes are not well matched (Bignardi et al., 2016); thus, we normalized the amplitudes of the synthetic H/V curves under DFA in the inversion to match the peak frequencies.

The results of the joint inversion are displayed in Figure 5. It shows that the joint inversion technique can retrieve a three-layer structure. Above the glacier ice half-space are a shallow layer at 27 m depth and a deeper layer at 115 m depth. These results suggest that the snow/firn/ice composition within the upper ice sheet column reflects a three-stage densification process, whose resulting stages can be separated by critical density values (Herron and Langway, 1980; Hörhold et al., 2011; Fujita et al., 2014). Specifically, the first stage (snow layer) is between the glacier surface and the depth associated with density of 550 kg m⁻³. This stage is characterized by rapid compression of low-density snow. The second stage (firn layer) is located between densities of 550 kg m⁻³ and 830 kg m⁻³, where the main mechanism of densification has been recrystallization of ice crystals, and pores of atmospheric air have closed off. And the third stage is the region with densities from 830 kg m⁻³ slowly increasing to 917 kg m⁻³, as trapped air becomes more and more compressed (Lundin et al., 2017). Consequently, strong impedance contrasts exist at the two interfaces separating the three stages, which we associate with the presence of the two clear H/V peaks observed in this study. In this context, we suggest that the 27 m depth shallow layer from the joint inversion results is the snow layer, and that the second layer between 27 m and 115 m is the firn layer.

Figure  5.  Results of the joint inversion of H/V and dispersion curves. In (a) and (b), the synthetic curves are coloured based on misfit values. The observed H/V and the dispersion curve are represented by solid black lines with error bars. The error of dispersion curve remains constant at 100 m s⁻¹. In (c), (d), and (e), the P-wave velocity, S-wave velocity, and density models are coloured corresponding to the synthetic curves. In each panel, the best inverted model is marked in red. Note that the amplitudes of all of the synthetic H/V curves are normalized by dividing by a factor of three the entire frequency band in panel (a).

DownLoad: Full-Size Img PowerPoint

To verify the estimated depth of the snow and firn layers, we compared our results to a density profile obtained from ground-penetrating radar at Kunlun Station, which is further validated by a density profile derived from a deep ice core drilling (DA2005) (Jiang S et al., 2012; Yang WX et al., 2021). The radar investigation yielded a depth of 23 m for the snow layer (corresponding to the critical density value of 550 kg m⁻³) and a depth of 98 m for the firn layer (critical density value of 830 kg m⁻³). Thus, results from our joint inversion for the upper densification structure of the ice sheet are consistent with results of these independent surveys at the same site. It also should be noted that the inverted depths correspond to layer interfaces with significant impedance contrasts, not necessarily the same depths indicated by the critical density values (Hörhold et al., 2011). In addition, although we suggest that the joint inversion alleviates the non-uniqueness problem, we acknowledge that inversion results’ uncertainty cannot be eliminated, because the current misfit function under DFA fits the curve, rather than giving more weight to fitting the peak frequency. In the future, we will explore ways to fit the the peak frequency better.

At Kunlun Station, a short-time (1-hour) ambient noise record from a seismic node was used to extract a clear broad frequency band H/V curve in which three clear peaks can be observed. The H/V curve carries information on both shallow (<120 m) and deeper (~2600 m) impedance contrasts (Piña-Flores et al., 2017; Spica et al., 2018). The lower (~0.17 Hz) peak corresponds to the ice-bedrock interface; the two higher peaks (~3 Hz and ~6 Hz) correspond to features of the upper structure of the ice sheet. The H/V method can be used to decipher deep ice structures from data collected by a single seismic nodal system. A small-length array with a short inter-station distance, using the same kind of sensors, is capable of simultaneously collecting sufficient data to allow construction of a dispersion curve that can reveal details of the upper structures. Joint inversion of the Rayleigh wave dispersion curve and the H/V curve were carried out to reduce the non-uniqueness problem of independent inversion and to enhance the reliability of the model results (García-Jerez et al., 2019).

5.   Conclusions

During China’s 39th expedition to Antarctica, a dense linear seismic array comprised of a portable seismic nodal system was deployed near Kunlun Station, Dome A, Antarctica, for ambient noise cross-correlation and dispersion curve measurements. In addition, a seismic node that recorded three-component ambient noise was adopted for H/V measurement. Continuous seismic ambient noise was recorded for 5 days between 5 and 10 January 2023 using the nodal systems. The H/V curve was calculated; three clear peaks were observed — a low frequency peak, at ~0.17 Hz, and two high frequency peaks, at ~3 Hz and ~6 Hz. We conducted H/V inversion under DFA in the 0.1–2 Hz frequency band to determine the ice sheet thickness, which we estimate to be 2577 m, consistent with the BedMachine ice thickness dataset. Furthermore, we found that a stratified ice sheet model can fit the observed H/V curve better than a single layer model. A dispersion curve in the 3–10 Hz frequency band was extracted from the NCFs using a phase-weighted slant-stacking technique. The dispersion curve and the H/V curve in the 2–9 Hz band were jointly inverted to determine the upper structure of the ice sheet; joint inversion was employed to reduce non-uniqueness in inversion and to improve model reliability. The inversion result revealed that the upper densification structure of the ice sheet consists of a surface snow layer of 27 m thickness above a firn layer extending down to a depth of 115 m. These findings align with density profiles derived from independent surveys at the same site. Our study demonstrates the versatility and broad band capability of the seismic nodal system to collect high quality data in the extreme glacier environment, as well as the potential of joint inversion of dispersion and H/V curves, in ice sheet structure investigations, due to the low cost of this approach, and the complementary advantages of the two curves.

Supplementary Materials

H/V Calculation

Figure S2 shows an example of seismic ambient noise recorded by the H/V node, where prolonged interfering signals can be found in the data. Figure S3 illustrates the H/V curves before and after the employment of the frequency-based window-rejection algorithm (Cox et al., 2020). After the rejection, the resulting H/V curve shows a smaller standard deviation of the mean frequency value, which goes from 0.06 to 0.02, of the lower peak at ~0.17 Hz. Although at ~0.17 Hz, the standard deviation curve of the mean H/V curve shows a larger amplitude range (2–6), the mean H/V peak frequency value is stable. The comparison of H/V results obtained from 9-hour and 1-hour ambient noise records shows no significant differences in the H/V curve results and mean peak frequency (Figure S3). This suggests a one-hour long record is adequate to obtain a stable H/V curve (Yan P et al., 2018; Zhang YH et al., 2023).

  S1.  The amplitude spectra of seismic ambient noise recorded by linear array node 01. The two horizontal-components have strong noise signal interference, making their amplitude larger than that of the vertical component. There is an order of magnitude difference in the root mean square (RMS) values of the three-component amplitude of the seismic signal. We believe that the issue is likely caused by the poor coupling between the instrument and the snow surface. This makes the H/V method not applicable.

DownLoad: Full-Size Img PowerPoint

  S2.  Example of seismic ambient noise recorded by the H/V node on 8 January 2023.

DownLoad: Full-Size Img PowerPoint

  S3.  The H/V curve calculated from 1-hour long seismic record between 19:00 to 20:00 on 8 January 2023. Panel (b) and (d) depict the H/V results before and after the employment of the frequency-based window-rejection algorithm, respectively. Vertical dashed line denotes the mean frequency, and red shaded area denotes standard deviation of mean frequency. Solid and dashed lines illustrate the mean H/V curve and its standard deviation curves. The white circles indicate the peak frequency of all the individual 60-second windows. And the green diamond indicates the peak frequency of the mean H/V curve. Panel (a), (c), and (e) are three-components seismic signals. After the rejection, the mean frequency value of the lower frequency peak is 0.17, and its standard deviation is 0.02. Before the rejection, mean frequency is 0.17, while its standard deviation is 0.06.

DownLoad: Full-Size Img PowerPoint

  S4.  Comparison of H/V results obtained from 9-hour and 1-hour ambient noise records. (a) The H/V curve calculated from 9-hour long seismic record between 13:00 to 24:00 on 8 January 2023. (b) The H/V curve calculated from 1-hour long seismic record between 9:00 to 10:00 on 8 January 2023.

DownLoad: Full-Size Img PowerPoint

  S5.  Noise cross-correlation functions (NCFs) for 1–10 Hz band vertical-component signals in the time-domain. The red lines denote upper envelope of the NCFs. The dashed lines denote velocities of 400 m s^–1 and 1700 m s^–1.

DownLoad: Full-Size Img PowerPoint

  S6.  (a) Vertical-component seismic ambient noise record from node 10 of the linear array on 5 January 2023. (b) Spectrogram of the seismic ambient noise in the 0.1–25 Hz frequency band.

DownLoad: Full-Size Img PowerPoint