4.1 Stationary-phase zone data selection

Near-surface imaging based on traffic noise is a passive source approach, the core principle of which is to obtain virtual shot gathers through seismic interferometry and to subsequently invert surface-wave dispersion curves to constrain the shear-wave velocity structure. According to the theory of seismic interferometry, the noise sources must be located within the stationary-phase zone of the geophone array (Wapenaar and others, Reference Wapenaar, Draganov, Snieder, Campman and Verdel2010); otherwise, strong energy sources outside this zone can lead to an overestimation of apparent velocities in the virtual shot gathers (Park and Miller, Reference Park and Miller2008; Chamarczuk and others, Reference Chamarczuk, Draganov, Quiros and Malinowski2022). The true velocity ${V_{{\text{true}}}}$ can be expressed as the product of the apparent velocity $V$ and $\cos \theta $, where $\theta $ denotes the angle between the source direction and the survey line (Morton and others, Reference Morton, Ivanov, Peterie, Miller and Livers-Douglas2021). In this study, because the roadway and the survey line intersect at an oblique angle, only traffic noise signals within the stationary-phase zone were selected for interferometric processing. To achieve this, the convoy length and speed were first estimated to determine the specific time periods during which the convoy entered and exited the stationary-phase zone.

In the recorded traffic noise data, the time of the maximum absolute amplitude ${T_{{\text{Peak}}}}$ at each geophone corresponds to the moment when the midpoint of the convoy is closest to that geophone, plus the seismic wave propagation time. Accordingly, the arrival time of the convoy midpoint at each geophone can be estimated as follows:

(1)\begin{equation}{T_n} = {T_{n{\text{Peak}}}} - \frac{{{h_n}\cos \theta }}{{{v_c}}},\end{equation}

where n denotes the geophone number, ${h_n}$ is the distance from geophone n to the roadway and $\theta $ represents the angle between the geophone array and the roadway. ${v_{{c}}}$ is the phase velocity of high-frequency surface waves, whose value can be roughly estimated by dispersion imaging of field data. The angle $\theta $ can be calculated using the slopes of the survey line ${m_A}$ and the roadway ${m_B}$ derived from the geophone coordinate positions, as given by Eqn (2). In this study, the angle between the survey line and the roadway is ∼8.23°:

(2)\begin{equation}\theta = \arctan \left( {\left| {\frac{{{m_A} - {m_B}}}{{1 + {m_A}{m_B}}}} \right|} \right).\end{equation}

The apparent average velocity of the convoy is estimated from the difference between the peak-amplitude times recorded at two selected geophones and their projected separation along the roadway direction, as expressed in Eqn (3). This projection accounts for the intersection angle between the survey line and the roadway.

(3)\begin{equation}v = \frac{{\left( {{n_2} - {n_1}} \right) \cdot x \cdot \cos \theta }}{{{T_{n2{\text{Peak}}}} - {T_{n1{\text{Peak}}}}}},\end{equation}

where x is the receiver spacing, ${n_2}$ and ${n_1}$represent the channel numbers of two receivers along the survey line and ${T_{n2{\text{Peak}}}}$ and ${T_{n1{\text{Peak}}}}$ denote the times of peak amplitude in the traffic noise recorded at ${n_2}$ and ${n_1}$. Using this method, the average velocity of the convoy in the dataset is calculated to be ∼2.636 m s⁻¹.

If ${T_{{\text{start}}}}$ and ${T_{{\text{end}}}}$ represent the times when the convoy enters and leaves the receiver array, respectively, the relative motion between the convoy and the receiver array can be expressed as follows:

(4)\begin{equation}\frac{{\left( {N - 1} \right) \cdot x \cdot \cos \theta + L}}{v} \approx {T_{{\text{end}}}} - {T_{{\text{start}}}},\end{equation}

where $\left( {N - 1} \right) \cdot x$ is the total array length and L is the convoy length, as illustrated in Fig. 4.

The duration ${T_{{\text{end}}}} - {T_{{\text{start}}}}$ is determined from the onset and termination times of the convoy-passage signal, following Rezaeifar and others (Reference Rezaeifar2023). Specifically, the absolute amplitudes of all receiver traces are summed at each time sample to form an array-stacked amplitude trace. This stack is used only to enhance the temporal envelope of the convoy passage. The array-stacked amplitude trace is then smoothed using a median filter and processed with the short-term average/long-term average (STA/LTA) algorithm, which is commonly applied in earthquake seismology to determine first arrivals, in order to identify ${T_{{{start}}}}$ and ${T_{{{end}}}}$, as shown in Fig. 5.

The convoy length L was obtained by substituting ${T_{{start}}}$ and ${T_{{end}}}$ into Eqn (4), yielding a value of ∼159.46 m.

The convoy speed and length were first derived from the recorded traffic-noise signals. The geometric relationship among the convoy, roadway, survey line and stationary-phase zone is illustrated schematically in Fig. 6. The arc-shaped boundary in Fig. 6 is used only to indicate the conceptual extent of the stationary-phase zone and is not drawn as a quantitatively scaled contour. The actual entry and exit times of the convoy through the stationary-phase zone were calculated from Eqns (5)–(10), and the corresponding time interval was then identified in the seismic records, as illustrated in Fig. 8 (You and others, Reference You, Mi, Guan, Zhang and Liu2023).

(5)\begin{equation}{T_1} = T - \frac{{{x_2} + L/2}}{v} + \frac{h}{{{v_{\text{c}}}\sin (\theta - {\theta _{\text{F}}})}}\end{equation}

(6)\begin{equation}{T_2} = T - \frac{{{x_1} - L/2}}{v} + \frac{h}{{{v_{\text{c}}}\sin (\theta + {\theta _{\text{F}}})}},\end{equation}

where

(7)\begin{equation}{x_1} \approx \frac{h}{{\tan \left( {\theta + {\theta _{\text{F}}}} \right)}}\end{equation}

(8)\begin{equation}{x_2} \approx \frac{h}{{\tan (\theta - {\theta _{\text{F}}})}}\end{equation}

(9)\begin{equation}T = {T_{{\text{peak}}}} - \frac{{h \cdot \cos \theta }}{{{v_{\text{c}}}}}\end{equation}

(10)\begin{equation}{\theta _{\text{F}}} = \sqrt {\frac{\lambda }{{2r}}}\end{equation}

where h represents the distance from the nearest geophone to the roadway; ${\theta _{{F}}}$ is the angular width of the first Fresnel zone (Sadeghisorkhani and others, Reference Sadeghisorkhani, Gudmundsson, Roberts and Tryggvason2017); $\theta $ is the intersection angle between the roadway and the geophone array, which was calculated as 8.23°; ${v_{{}}}$ is the convoy velocity, estimated to be ∼2.636 m s⁻¹; r is the total length of the array; and $\lambda $ is the characteristic wavelength, which was calculated using the velocity corresponding to the highest frequency.

Since the roadway intersects the geophone array, the geophones are distributed on both sides of the road. Therefore, the array was divided into two subarrays, labeled as array ① and array ②, as shown in Fig. 4. The time interval from ${T_1}$ to ${T_2}$ corresponds to the period during which the convoy traveled within the stationary-phase zone of the geophone array. The values of ${T_1}$ and ${T_2}$ were calculated separately for arrays ① and ②, and the seismic data within each respective ${T_1}$– ${T_2}$ interval were extracted for subsequent seismic interferometry processing.

In Fig. 7, the orange sections of the seismic records represent the data generated while the convoy traveled through the stationary-phase zones associated with geophones No. 30 and No. 60. These time windows were determined using Eqns (5) and (6). The stationary-phase data from all geophone channels in both arrays were then extracted and used for subsequent cross-coherence computations.

4.2 Virtual shot gather obtained by cross-coherence

Cross-coherence (Eqn 11) was applied to both the complete traffic-noise dataset and the data filtered by the stationary-phase zone (the orange sections in Fig. 7) to generate the virtual shot gathers.

(11)\begin{equation}C\left( \omega \right) = \frac{{u\left( {{X_A},\omega } \right){u^ * }\left( {{X_B},\omega } \right)}}{{\left| {u\left( {{X_A},\omega } \right)} \right|\left| {u\left( {{X_B},\omega } \right)} \right|}}\end{equation}

(12)\begin{equation}C\left( \omega \right) = \frac{{G\left( {{X_A},\omega } \right){G^ * }\left( {{X_B},\omega } \right)}}{{\left| {G\left( {{X_A},\omega } \right)} \right|\left| {G\left( {{X_B},\omega } \right)} \right|}},\end{equation}

where $u\left( {{X_A},\omega } \right)$ and $u\left( {{X_B},\omega } \right){\text{ }}$denote the Fourier transforms of the seismic signals recorded at points ${X_A}$ and ${X_B}$, while $G\left( {{X_A},\omega } \right)$ and $G\left( {{X_B},\omega } \right)$ represent the corresponding Green’s functions. Because Eqn (11) does not depend on the source term, this method suppresses source-related effects and enhances the signal-to-noise ratio (Aki, Reference Aki1957).

Figures 8 and 9 present the virtual shot gathers of arrays ① and ②, in which the surface waves constitute the dominant energy component. Linear coherent phases with apparent velocities significantly higher than those of the surface waves can be observed. These features are consistent with the P-guided waves observed in active-source seismic records from ice-covered regions (Li and others, 2021; Nakata and others, Reference Nakata2023; Zhang and others, Reference Zhang2024; Shi and others, Reference Shi, Li and Chen2025). This observation indicates that traffic noise provides stronger seismic energy than conventional passive source approaches. Comparisons between Figs 8a,d and 9a,d show that the apparent velocities of the surface waves (indicated by red dashed lines) in the full-data reconstructed virtual shot gathers are higher than those derived from the stationary-phase data. This agrees with the theoretical expectation that sources located outside the stationary-phase zone result in overestimated apparent velocities.

4.3 Dispersion spectrum calculated by Frequency-Bessel Transform

The Frequency-Bessel Transform was used to calculate the dispersion spectra (Wang and others, Reference Wang, Wu and Chen2019), and the computation employed the code developed in Li and others (Reference Li2021). The results are presented in Fig. 10. Figure 10a,b display the dispersion spectra of arrays ① and ② using the complete traffic-noise data, while Fig. 10c,d show those obtained from the stationary-phase data. Comparison indicates that the spectra derived from the full traffic-noise data are strongly affected by noise, leading to low signal-to-noise ratios and discontinuous P-guided wave dispersion. In contrast, the stationary-phase data produce more continuous and distinct surface and P-guided wave dispersion energy bands with minimal noise interference and improved signal-to-noise ratios.

4.4 Monte Carlo inversion of dispersion curves

This study concentrates on retrieving subsurface information from surface-wave dispersion. Accordingly, only the surface-wave energy in the dispersion spectra is analyzed. For each virtual shot gather, dispersion spectra were calculated, and the fundamental-mode surface-wave dispersion curves were then extracted through manual picking, as illustrated in Fig. 11. The high-frequency portion of the dispersion curves reaches up to 27 Hz, indicating that traffic noise provides higher-frequency content and thus allows for higher-resolution imaging of shallow structures compared with traditional passive seismic methods.

The Monte Carlo algorithm was employed to invert the fundamental-mode surface-wave dispersion curves. This algorithm is a global optimization method based on probabilistic simulation, which performs random sampling of the parameter space and evaluates the misfit between forward-modeled and observed data to achieve nonlinear inversion of the dispersion curves (Sambridge and Mosegaard, Reference Sambridge and Mosegaard2002; Shapiro and Ritzwoller, Reference Shapiro and Ritzwoller2002; Maraschini and Foti, Reference Maraschini and Foti2010; Aleardi and Stucchi, Reference Aleardi and Stucchi2021). The initial model parameters and inversion settings are summarized in Table 1.

As shown in Fig. 12, the S-wave velocity model that best fits the theoretical and observed dispersion curves was selected. By inverting each dispersion curve, a 1-D S-wave velocity structure at each imaging point was obtained as illustrated in Fig. 12b. All 1-D S-wave velocity models derived from the inversion were then spatially interpolated using the Local Polynomial interpolation method to construct a two-dimensional S-wave velocity profile of the survey area (Fig. 13).

Figure 12.

(a) Dispersion-curve fitting results from Monte Carlo inversion. Observed dispersion curve (DC) represents the DC derived from the field data, and Best DC denotes the DC corresponding to the optimal S-wave velocity model obtained from the inversion. (b) The red solid line indicates the optimal inverted S-wave velocity model, and the black dashed lines show the range of the initial inversion model.

Figure 12 Long description

The image contains two panels labeled a and b. Panel a shows dispersion-curve fitting results from Monte Carlo inversion. The horizontal axis is Frequency (Hz), ranging from 0 to about 25 Hz, and the vertical axis is Phase velocity (m/s), ranging from about 1300 to 2000 m/s. The observed dispersion curve derived from the field data and the best-fit dispersion curve corresponding to the optimal S-wave velocity model are plotted over multiple modelled dispersion curves. The modelled curves are colored according to the misfit value shown by the vertical color scale on the right. The dispersion curves generally decrease in phase velocity with increasing frequency, from about 1900 m/s at low frequency to about 1350 m/s at high frequency. Panel b shows the inverted S-wave velocity model as a stepped depth profile. The horizontal axis is shear-wave velocity in m/s, and the vertical axis is Depth (m), ranging from 0 to 120 m. The red solid line indicates the optimal inverted S-wave velocity model, while the black dashed lines indicate the minimum and maximum bounds of the initial inversion model. The optimal model generally increases in shear-wave velocity with depth.

Figure 13.

Two-dimensional S-wave velocity profile constructed by interpolating all one-dimensional inverted models using the Local Polynomial interpolation method. The left part of the profile corresponds to array ② (geophones No. 1–44) located on the left side of the roadway, and the right part corresponds to array ① (geophones No. 45–100) located on the right side of the roadway.

Figure 13 Long description

The image shows a two-dimensional S-wave velocity profile constructed by interpolating one-dimensional inverted models using the Local Polynomial interpolation method. The profile is divided into two adjacent parts separated by a central label, Road, with dashed lines indicating the roadway position. The left part corresponds to array ②, consisting of geophones No. 1–44 located on the left side of the roadway. Its horizontal axis is Inline (m), with tick marks from about 100 to 700 m. The right part corresponds to array ①, consisting of geophones No. 45–100 located on the right side of the roadway. Its horizontal axis is Inline (m), with tick marks from about 1200 to 1900 m. Both parts share the same vertical depth scale, labeled Depth (m), ranging from 0 to 120 m. The color scale below the profiles is labeled Vs (m/s), with values from about 1350 to 2550 m/s. In both parts, the S-wave velocity structure varies mainly with depth, with relatively lower velocities in the shallow layer and higher velocities at greater depth. Minor lateral variations are visible along the inline direction.

To investigate the density variation of the overlying snow and firn layers with depth, this study adopts the empirical relationship between shear-wave velocity and density proposed by Diez and others (Reference Diez2014) based on in situ measurements from the Colle Gnifetti Glacier. The relationship was derived by applying the Herglotz–Wiechert inversion to diving-wave travel-time data combined with density measurements from the KCI ice core, yielding a functional relationship between ${v_s}\left( z \right)$ and ${\rho _{\left( z \right)}}$:

(13)\begin{equation}\rho \left( z \right) = \frac{{{\rho _{{\text{ice}}}}}}{{1 + {{\left[ {\left( {{v_{{\text{s}},{\text{ice}}}} - {v_{\text{s}}}\left( z \right)} \right)/950} \right]}^{1.17}}}},\end{equation}

where ${\rho _{\left( z \right)}}$ represents the density at depth z (kg m⁻³), ${\rho _{{{ice}}}}$ is the density of ice and ${v_s}\left( z \right)$ denotes the shear-wave velocity at that depth. This empirical model effectively characterizes the density distribution of the snow and firn layers and is particularly applicable in cold glacier environments, where it captures significant density gradients within the firn.

As shown in Fig. 14a, the density-depth curve (red dashed line) was derived from the shear-wave velocity curve (blue solid line) using this relationship. To evaluate the applicability of Eqn (13) to this study area, the calculated velocity–density–depth curve (Fig. 14a) was compared with the borehole stratigraphic data from the same site (Fig. 14b). The comparison shows strong consistency in both trend and key stratigraphic boundaries. In the surface snow layer, the modeled density of ∼550 kg m⁻³ agrees precisely with borehole measurements. This relatively high near-surface density suggests that the shallow snow has already undergone substantial compaction and may also reflect local firn modification associated with wind packing, meltwater refreezing or near-surface mechanical disturbance, rather than representing freshly deposited snow. Within the firn layer (6–51 m), the theoretical density smoothly increases from 550 to ∼830 kg m⁻³. The density of ∼830 kg m⁻³ corresponds to the firn–ice transition, at which the interconnected air-filled pore spaces are largely closed off and firn is transformed into glacier ice. Between 51 and 90 m, the modeled density remains above the pore-close-off threshold and gradually approaches the density of pure glacier ice. This interval therefore represents the post-close-off transition zone, where the firn has already transformed into glacier ice but residual air bubbles and remaining porosity continue to be compressed with increasing depth. Glaciologically, this section records the final stage of snow-to-ice transformation rather than ordinary firn densification. Below ∼90 m, within the compact glacial ice layer, the modeled density (∼917 kg m⁻³) is very close to the measured mean of 901 kg m⁻³, with a relative error of only ∼1.8%. Overall, the model captures the main increasing density trend with depth and provides a continuous estimate of the snow–firn–ice transition structure. Rather than relying only on discrete borehole observations, this velocity-derived density profile offers a way to trace the progressive densification process and to infer the approximate depth range over which pore close-off and the transition toward glacier ice occur. These results further suggest that, when calibrated against local stratigraphic constraints, the Diez et al. relationship can provide a useful seismic proxy for interpreting firn densification and the transition toward compact glacier ice in polar ice-sheet settings.

Figure 14.

(a) Shear-wave velocity profile with depth (blue solid line) and density profile derived from empirical formula (red dashed line); (b) Borehole log of the snow-ice-bedrock sequence at the study site. Source: Modified from Talalay and others, Reference Talalay2025.

Figure 14 Long description

The image contains two panels labeled a and b. Panel a shows shear-wave velocity and density profiles with depth. The vertical axis is Depth (m), ranging from 0 to 120 m. The lower horizontal axis is labeled Shear Wave Velocity Vs (m/s), with values from about 1000 to 2500 m/s, and the upper horizontal axis is labeled Density ρ (kg m−3), with values from about 500 to 2000 kg m−3. The blue solid curve represents the shear-wave velocity profile, and the red dashed curve represents the density profile derived from the empirical formula. Additional average curves are shown in the legend as Vs avg and ρ avg. The profiles generally increase with depth, and the plot is divided into snow, firn and glacial ice zones. A horizontal dashed line near 90 m marks the approximate transition to compact glacial ice. Panel b shows a borehole log of the snow–ice–bedrock sequence at the study site. The diagram is titled Stratigraphy and density and includes labeled layers of snow, firn and glacial ice, with boundary depths marked at 6 m, 51 m and 89 m. Density values are shown as 550 kg m−3, 830 kg m−3 and 901 kg m−3. Near the bottom of the borehole log, basal ice and bedrock are labeled with depths around 540.93 m, 541.12 m and 541.60 m. Overall, the figure compares the velocity-derived density profile with the borehole stratigraphic information from the study site.