High‐Resolution Characterization of the Firn Layer Near the West Antarctic Ice Sheet Divide Camp With Active and Passive Seismic Data

We construct a high‐resolution shear‐wave velocity (VS) model for the uppermost 100 m using ambient noise tomography near the West Antarctic Ice Sheet Divide camp. This is achieved via joint inversion of Rayleigh wave phase velocity and H/V ratio, whose signal‐to‐noise ratios are boosted by three‐station interferometry and phase‐matched filtering, respectively. The VS shows a steep increase (0.04–0.9 km/s) in the top 5 m, with sharp interfaces at ∼8–12 m, followed by a gradual increase (1.2–1.8 km/s) between 10 and 45 m depth, and to 2 km/s at ∼65 m. The compressional‐wave velocity and empirically‐obtained density profile compares well with the results from Herglotz–Wiechert inversion of diving waves in active‐source shot experiments and ice core analysis. Our approach offers a tool to characterize high‐resolution properties of the firn and shallow ice column, which helps to infer the physical properties of deeper ice sheets, thereby contributes to improved understanding of Earth's cryosphere.


Introduction
Antarctic ice sheets preserve important information about past environmental conditions, climate change impacts, and necessary information for predicting future sea-level rise (e.g., Fegyveresi et al., 2011;Schmidt et al., 2023).Of particular note is the West Antarctic Ice Sheet (WAIS), which has long been considered sensitive to warming of the ocean waters (Joughin et al., 2014;Rignot et al., 2014).Projections estimate a contribution to sea-level rise of up to 3 m if WAIS were to collapse (Holland et al., 2019).Overlying the ice sheet is the firn layer, an intermediate phase between fresh snow and glacial ice (Cuffey & Paterson, 2010).Characterized by its substantial air content, the near-surface firn layer thermally and mechanically isolates the underlying glacial ice from external environmental influences.Firn properties are very sensitive to surface climate change, such as temperature variation, surface accumulation, and wind patterns (Herron & Langway, 1980;Ligtenberg et al., 2011;Maeno & Ebinuma, 1983;Wilkinson, 1988), and its properties (e.g., depth, density) are essential components for estimating ice sheet mass and understanding ice-atmosphere interactions (e.g., Helsen et al., 2008;van Wessem et al., 2017).
Seismological tools, including active-source surveys (e.g., Agnew et al., 2023;Horgan et al., 2011;Karplus et al., 2024;Pearce et al., 2023) and ambient noise tomography (Li et al., 2021;e.g., Nakata et al., 2019;Pyle et al., 2010;Shapiro et al., 2005), have been valuable for probing the internal structures of Antarctic ice sheets (e. g., Aster & Winberry, 2017;Podolskiy & Walter, 2016).However, challenges persist in accurately characterizing firn properties (Maupin, 2017;Tanimoto et al., 2013).For instance, Zhang et al. (2022) constructed a shear-wave velocity (V S ) model in the upper 300 m through ambient noise tomography using data from a seismic array close to the WAIS Divide camp (Figure 1a), acquired as part of the Thwaites Interdisciplinary Margin Evolution (TIME) project (hereafter termed the "TIME array").However, the model of Zhang et al. (2022) offers little constraint of the uppermost firn layer given limited data quality and sensitivity to shallow structures, and the predicted Rayleigh wave phase velocities (V ph ) consistently exceed observations at frequencies >25 Hz.This also reduces the value of the method for further modeling studies.
To improve the resolution of the uppermost firn (the so-called firn-air layer), we now use the Rayleigh wave H/V ratio, defined as the ratio of horizontal and vertical amplitude of the elliptical particle motion.This method has much higher sensitivity to shallow structures than V ph alone (e.g., Berbellini et al., 2019;Jones, Kulessa, et al., 2023;Jones, Ferreira, et al., 2023).We therefore undertake a joint inversion of V ph and H/V ratios of Rayleigh waves, extracted from ambient noise cross-correlation (ANC), to enhance the resolution of V S in the uppermost 20 m.Additionally, we employ advanced denoising techniques (Herrin & Goforth, 1977;Qiu et al., 2021) on ANCs before measuring Rayleigh wave V ph and H/V ratios to mitigate measurement errors, which fortifies the reliability of subsequent inversion results.

Data
The TIME array operated during Jan 5-16th, 2019, and consisted of two parts: 24 seismic nodes arranged at 15°i ntervals in a circle of radius of 500 m, and 76 nodes installed at 30 m intervals along a 1-D linear array of ∼2.2 km length (Figure 1a).At each station, a Magseis Fairfield 5-Hz 3-component seismic node is installed (Ringler et al., 2018), programmed with 1,000 Hz sampling rate, pre-amp gain of 12 dB, linear-phase Nyquist filter, and DC offset removal.In this study, only the 1-D linear array is analyzed.
We calculate ANCs between 0.5 and 40 Hz and follow the methodology of Berg et al. (2018) to preserve relative amplitude between components (Text S1 in Supporting Information S1).In Figure 1b, surface waves traveling at ∼2 km/s are observed with larger amplitudes at negative time lag (i.e., acausal part), while high-frequency surface waves appear only in the positive time lag between stations less than 1.5 km apart.This asymmetry of surface wave amplitudes has been widely observed in previous studies (Shapiro et al., 2005;Wang et al., 2019), likely due to the uneven distribution of noise sources (Chaput, Aster, Karplus, & Nakata, 2022, Chaput, Aster, Karplus, Nakata, et al., 2022).Given the complexity of signals with low signal-to-noise ratio (SNR) in the positive time lag of ANCs, we use the acausal part for a more robust analysis.
The ANCs show multiple phases traveling at different speeds (labels 1-3 in Figures 1c and 1d).Phase 1 is the most prominent signal in the 3-30 Hz range and represents the fundamental-mode Rayleigh wave with an array-mean V ph dispersion consistent with that of Zhang et al. (2022).The two other phases are less dispersive and apparent at higher frequencies (Phase 2: 20-30 Hz with V ph ≈1.8 km/s; Phase 3: ∼20 Hz with V ph ≈3.8 km/s).Phase 3 is absent from the TT component.Polarization analysis (Text S1 in Supporting Information S1) suggests Phase 2 may correspond to a higher mode Rayleigh wave, and Phase 3 is likely a P wave or leaky mode.

Ambient Noise Tomography
Here we briefly describe the denoising techniques (Texts S2.1-S2.2 in Supporting Information S1).Leveraging the 1-D linear configuration of the array, we first employ three-station interferometry (Qiu et al., 2021) to enhance fundamental-mode Rayleigh waves.Based on ray theory, for surface waves traveling between three stations i, k and j (i < j) along a straight line, the travel times should satisfy but the relation does not hold for diving body and scattered waves.Therefore, surface waves between stations i and j can be constructed by interferograms between stations ik and jk, and stacking the interferograms from all available k gives the denoised ANC between i and j.The SNR of fundamental-mode Rayleigh waves markedly  (d).Multichannel analysis diagrams, where waveforms are slant-stacked with different velocities in the frequency domain (Qiu et al., 2021).
improves in the denoised wavefield (Figure 2a).The other two phases, visible in the raw ANCs (Figure 1c, RRcomponent), can no longer be seen after denoising since this approach does not preserve amplitude information and only retain the dominant surface wave signal -the fundamental-mode Rayleigh wave (Qiu et al., 2021).Thereafter, we adopt an automatic frequency-time analysis to measure dispersion curves for all station pairs and apply eikonal tomography (Lin et al., 2008;Qiu et al., 2019) to derive phase velocity maps at frequencies between 3 and 30 Hz.The resulting V ph model (Figure 3a), consistent with the findings of Zhang et al. (2022), lacks lateral heterogeneity and gradually decreases from 1.8 km/s to 1.2 km/s as frequency increases.
For H/V ratios, we employ a different noise reduction strategy, which involves phase-matched filtering (Herrin & Goforth, 1977) to improve measurement qualities.The strategy aligns fundamental-mode Rayleigh waves at different frequencies to zero lag using the obtained phase dispersion measurements.Figure 2b shows the phasematch-filtered ANCs for station pairs associated with a specific virtual receiver at RR and RZ components, respectively.We taper the phase-match-filtered ANCs within the time range of 0 ± 0.1 s, as other signals are expected to deviate from zero lag (Herrin & Goforth, 1977).Given the localized sensitivity of H/V ratio to shallow structures right beneath the receiver (e.g., Maupin, 2017), H/V ratios measured from ANCs associated with the same receiver but different virtual sources (Figure 2c) appear mutually consistent.Thus, we measure H/V ratio from the waveform stacked over phase-match-filtered ANCs of all virtual sources for each target receiver (Figure 2c).The obtained H/V ratio is linearly interpolated onto a regular 2-D frequency-space grid.The station spacing (30 m) is used as the spatial grid size.Unlike the V ph model, the H/V ratio demonstrates a more pronounced lateral variation (Figure 3a).
Finally, we employ the improved Neighborhood Algorithm (Wathelet, 2008) to jointly invert V ph and H/V ratios for piecewise 1-D V S profile (Text S2.3 in Supporting Information S1). Figure 3a shows the optimal 2-D V S model, comprising 1-D V S profiles from all grids: this shows prominent vertical variations and subtle lateral heterogeneities.

1-D P-Wave Travel-Time Inversion
Considering the limited sensitivity of V ph and H/V ratio to the P-wave velocity (V P ), we validate the model inverted from Section 2.2 with an independent V P model from diving waves in active-source shot experiments, acquired along the linear section of the TIME array (Text S3 in Supporting Information S1).First-break traveltime picks, as a function of offset (Tx) (Figure 3b), are supplied to a Herglotz-Wiechert inversion (HWI).HWI places stringent constraints on the smoothness and continuity of Tx data and their derivatives (Aki & Richards, 1980).We follow, in part, a HWI workflow widely-applied in cryospheric studies (e.g., Brisbourne et al., 2020;King & Jarvis, 2007) but deviate from the way that Tx picks are smoothed (specifically, using a Levenberg-Marquardt damped least-squares fit to match an analytic Tx form ("KB90"; Kirchner & Bentley, 1990)).Here, Tx is smoothed using a more generalized Monte Carlo (MC) approach of McKeague (2023), to optimally satisfy HWI criteria for Tx but preferentially avoiding KB90's prescribed form for Tx (and, consequently, for the final velocity-depth model).The output curves from KB90 and MC fits are shown in Figure 3b.
The MC analysis converged on a solution after 10 5 iterations, and improves the root-mean-square (RMS) fit to the observed travel-time data by ∼25% compared to KB90 (Figure 3b).The 1-D velocity model output from HWI is presented in Figure 4a.

Results
The 1-D V S profile (Figure 4a), averaged over all grids, illustrates notable velocity variations within the top ∼10 m.V S increases from 0.04 km/s to 0.9 km/s in the uppermost ∼5 m, followed by a sharp velocity increase to ∼1.2 km/s at ∼8 m depth.Below this, V S gradually increases, reaching 1.8 km/s by ∼45 m depth and 2 km/s by ∼65 m.V S is virtually constant in the deepest section, corresponding to glacier ice.Compared to the results of Zhang et al. ( 2022) (Figure 4a), our Vs model reveals more intricate details of shallow firn structures within the uppermost ∼20 m.The enhanced resolution is attributed to the inclusion of H/V ratios in the tomography (Text S2.3 in Supporting Information S1), and the improved SNR of measurements achieved from advanced denoising techniques.Although the V P profile from the inversion (Figure 4a) has inherent uncertainties due to limited sensitivity of V ph and H/V ratio to V P , it still provides a reliable approximation of the subsurface structure beneath the TIME array.Figure 4a compares the average V P profile from surface wave inversion (orange solid line) with the V P profile derived from active-source shot experiments (orange dashed lines).Both V P models are similar, implying a rapid increase from 0.5 km/s to 2.2 km/s in the top 8 m, a further increase to 3.1 km/s at ∼45 m depth, and reach 3.4 km/s at ∼60 m depth.This consistency suggests that the inversion has converged to a desirable status, affirming the reliability of the obtained imaging results.
The density of the firn column can be evaluated from V S using the empirical relationship of Diez et al. (2014).
where ρ ice and V S.ice are the density and shear-wave velocity of glacier ice, respectively.We set ρ ice as 917 kg/m 3 , a typical value for glacier ice, and V S. ice as 2.0 km/s, the mean value of V S below 70 m.The resulting density profile (Figure 4a) exhibits large gradients from ∼320 kg/m 3 to ∼500 kg/m 3 in the top ∼8 m followed by a gradual increase in the deeper section.
Applying the effective medium approximation, we estimate the porosity (P) as where ρ air = 1.25 kg/m 3 is the density of air.In Figure 4a, the porosity approaches 70% near the surface, consistent with an expected 1:2 volume ratio of ice and air in the firn (Arthern et al., 2010).Then the porosity decreases to 45% within the top ∼8 m, corresponding to the initial stage of firn densification when ρ ≤ 506 kg/m 3 , mainly controlled by grain settling and packing.Further compression of air from the ice continues until reaching the pore close-off depth Z PCOD , where no air connection to the surface remains and the density increases to pore close-off density ρ PCOD .ρ PCOD is site dependent and usually between 800 and 830 kg/m 3 (Cuffey & Paterson, 2010).This process is often referred to as second stage of firn densification (Herron & Langway, 1980), and multiple mechanisms may have contributed to the depletion of air content.
Analysis of ice-core samples at the study site suggests ρ PCOD and Z PCOD are ∼823 kg/m 3 and ∼57 m, respectively (Fegyveresi et al., 2011).Assuming ρ PCOD = 823 kg/m 3 , we determine Z PCOD at ∼49 m using density profile in Figure 4a.Similarly, Diez et al. (2016) estimate Z PCOD at ∼47 m, and Zhang et al. (2022) at ∼62 m, with ρ PCOD of 830 kg/m 3 based on V S imaging results from Rayleigh wave phase velocity.In comparison, model calculations at the study site suggest Z PCOD at ∼40-60 m (Battle et al., 2011;Kaspers et al., 2004).We note the pore close-off process is not instantaneous, thus these estimations are overall consistent with each other, but we estimate the depths with better resolution.
Density variation over depth due to change in air space can be expressed as (Herron & Langway, 1980) where C is a site-dependent constant and estimated by linear least-squares regression of the ln ρ ρ ice ρ curves for the two stages of densification (Figure 4a).Then we calculate the ice age using (Text S4 in Supporting Information S1) ( where ρ 0 = 327 kg/m 3 is the zero-depth density, and A = 0.22 m • year 1 (Kaspari et al., 2004) is the accumulation rate.Here k is the rate constant, and can be estimated via assuming an Arrhenius type temperature sensitivity for compaction (Herron & Langway, 1980).R = 8.314 J • K 1 mol 1 is the gas constant, and T = 30°C is the annual mean temperature at WAIS (Orsi et al., 2012;Steig et al., 2013).The constant a can be determined as The estimated ice age is shown as color of the porosity line in Figure 4a, in years prior to the data collection.The age stops at 142 years at ∼47 m depth, where the density exceeds 800 kg/m 3 .We also calculate age of ice based on the mean annual layer thickness (∼0.24 m/year) from ice core analysis at WAIS (Sigl et al., 2016), with differences never exceeding 20 years (inset in Figure 4a).These differences could result from the error of densityvelocity empirical relation or the densification model (e.g., Arthern et al., 2010;Thompson-Munson et al., 2023),

Geophysical Research Letters
10.1029/2024GL108933 but they nonetheless remain within the measurement uncertainty of ice core analysis (Sigl et al., 2016).This consistency further validates the reliability of our velocity models.

Discussion
Using only one week of passive seismic data recorded at surface nodes, our study improves the resolution within the firn layer, particularly in the upper ∼20 m, compared to previous seismic studies.The model compares well with the P-wave velocity profile from active-source shot experiments.Figure 4b presents a schematic illustration of the ice structure, highlighting the snow compaction process within the firn.This achievement is attributed to the joint inversion of high-quality V ph and H/V ratio after applying advanced denoising procedures to ANCs.
Furthermore, the Vs model can help capture and interpret critical aspects of firn dynamics with a robust empirical velocity-density relation.
One interesting finding of this study is the lateral heterogeneities of the V S model.Given the laterally homogeneous nature of V ph (Figure 3a), the lateral variability in the 2-D V S model primarily stems from the H/V ratio variation.This observation raises questions about whether the lateral variability in the V S model reflects true firn structure or artifacts induced by uncorrelated noise.To address this question, we conduct forward modeling of surface waves (Text S5 in Supporting Information S1).This demonstrates that the lateral variations in H/V ratio obtained from a laterally homogeneous V S model are negligible compared to those observed in field data, assuming a Gaussian-type normal distribution for uncorrelated noise.Lateral V S variations are therefore unlikely to be entirely due to artifacts, indicating a degree of lateral variability in the shallow firn layer.These variations could be induced by various erosion and deposition patterns from localized surface processes, such as winddriven snow drifting.Notably, variable firn depth and density have been documented within a range of few meters (Weinhart et al., 2020;Wever et al., 2022).
Although ice-cores will always provide valuable ground-truth data, the accuracy of our results potentially allows seismic observations to replace them in data-poor areas.Our estimates of density, for example, can help calibrate mass balance estimations from satellite altimetry, where variations in surface height are used to infer alterations in ice mass or firn density changes (e.g., Holland et al., 2011).Furthermore, inferences from a reliable density profile could reliably include age estimates.The key consideration is with a deployment of passively-recording seismic instrumentation, we have a means of long-term ice-sheet monitoring that combines high-spatial and hightemporal resolution of ice sheet processes, requiring only the field logistics to deploy and service sensors.
In this study, higher mode surface waves are seen in ANCs at frequencies above 20 Hz (Figures 1c and 1d).Performing a multi-modal surface wave inversion can further improve the velocity model of the ice structure beneath the TIME array.Measurements of high-quality phase and amplitude information of higher mode surface waves require the fundamental mode signals being removed from the ANCs (Qiu et al., 2021).By subtracting fundamental mode surface waves given by the velocity model from this study, the joint inversion of multi-modal surface wave phase velocity and H/V dispersion can image the ice structure beneath the TIME array with higher spatial resolution.This will be a subject of future study.

Conclusion
In this study, we obtain high-resolution V S model for ice structures in the uppermost 100 m near the WAIS Divide camp using 1-week ambient seismic data collected from a dense linear 1-D nodal array.To boost the energy of the fundamental-mode Rayleigh wave and enhance the precision of the resulting Vs model, we employ three-station interferometry and phase-matched filtering on ANCs to mitigate measurement errors in V ph and H/V ratios, respectively.The obtained high-resolution V S model reveals prominent vertical variations and subtle lateral heterogeneities in firn structure, the latter shown to be at least partially indicative of firn structure through forward modeling.The V P model inferred from the ambient noise tomography, although with large uncertainties, is consistent with results derived from P-wave travel-time inversion of active-source data.The density profile, empirically converted from the V S model, provides valuable insights into the firn densification process.It suggests a pore close-off depth at ∼49 m, consistent with previous studies.Assuming steady state compaction, we derive the ice age up to 142 years in the top 49 m, compatible with measurements from ice core analysis.The method can be applied to other regions in Antarctica with dense array deployments and easily implemented to long-term ice sheet monitoring.By enhancing our understanding of firn layer properties, the tomographic tool developed in this study can significantly benefit studies related to dynamic modeling and broader seismic investigations of glaciers.

Figure 1 .
Figure 1.Array information and ANCs.(a).The Thwaites interdisciplinary margin evolution array layout.The 1-D array stations are plotted in red and named according to their distance in meters from the center of the circular array, ranging from 500 to 2750 with a 30-m increment.The inset provides an overview of the array location (red triangle) on a large-scale map.(b).ANCs for RR (radial-radial), RZ (radial-vertical) and TT (transverse-transverse) components.Potential reflections in ANCs are indicated by black arrows.(c).ANCs at 20 Hz after narrow-bandpass filtering, annotated with phases 1, 2, and 3.(d).Multichannel analysis diagrams, where waveforms are slant-stacked with different velocities in the frequency domain(Qiu et al., 2021).

Figure 2 .
Figure 2. Denoising surface waves.(a).Narrow-bandpass filtered RR-component ANCs at 20 Hz after three-station interferometry.(b).Phase-match-filtered RR-and RZ-component ANCs at station 1220.The black curves show the stacked ANCs, and blue vertical lines represent the tapering time window (i.e.0 ± 0.1 s).Stations less than 0.2 km away from 1220 are removed (white area).Negative distance means stations are closer to the origin than 1220.(c).H/V ratio at station 1220, measured from individual ANCs (gray curves) and stacked ambient noise cross-correlation (red curve) after phase-match filtering.

Figure 3 .
Figure 3. (a).Tomography results: (i) V ph , (ii) H/V ratio and (iii) inverted V S model.The black dotted lines in (iii) depict V S contour lines of 1.2 km/s and 1.8 km/s.(b).Travel-time analysis of active-source data: (i).Shot record with first-break picks, Tx (red annotation).(ii).Comparative fit to Tx data of KB90 (black) and Monte Carlo (MC) (orange) smoothing functions.For clarity, the time axis is shown as reduced time, assuming an average velocity of 3,200 m/s.(iii).Fitting error of KB90 (black) and MC (orange) curves.Horizontal lines show RMS errors of KB90 and MC (0.08 and 0.06 ms, respectively).

Figure 4 .
Figure 4. 1-D profiles averaged over all grids.(a).(i): Average V S and V P profiles from this study (red and orange solid curves), and V S from Zhang et al. (2022) (red dashed curve), V P from active-source Herglotz-Wiechert inversion (dashed orange curve).The blue solid line represents the empirically obtained density.The small inset presents enlarged velocity profiles in the top 20 m. (ii): Depth profile of ln ρ ρ ice ρ (black dotted line) with best-fitting linear trends (blue and green lines) and porosity (solid line) colored with ice age for ρ ≤ 800 kg/m 3 .The inset shows the difference in ice ages ("error") between estimations made here and from ice core analysis.(b).Schematic cartoon of the ice column structure.The vertical axis is not to scale.The depth of bedrock is obtained from Zhang et al. (2022).