Long-term underwater noise measurements were conducted in the East Siberian Shelf over the course of approximately one year, spanning from August 22, 2019 to August 13, 2020. An autonomous passive acoustic recorder (AURAL-M2, Multi-Electronique Inc.) was moored 13 m above the seafloor at location 74° 37.327’N, 174° 56.397’E in waters with a depth of approximately 70 m ( Figure 1A ). Acoustic data were recorded for 10 minutes every hour at a sampling rate of 32,768 Hz. From September 2 to 10, 2019, the R/V Araon operated the airgun (GI-SOURCE 355, Sercel) for underwater geological survey in the Chukchi-East Siberian Continental Margin (Jin, 2020). Low-frequency impulsive airgun sounds were emitted by the airgun shots at approximately 11-s intervals while maintaining a constant firing depth of ~6 m. The waveform of the airgun shot, measured from a near-field hydrophone (AGH-7100-C, Geophysical Products Inc.) positioned approximately 1 m above the airgun, exhibited a spike-shaped pattern with a peak frequency of ~27 Hz and dominant energy concentrated below 300 Hz (Han et al., 2023).

During the seismic survey, the R/V Araon approached the receiver system from the northeast direction, coming closest at 15:47 on September 9, after which the vessel moved away in the southwest direction ( Figure 1B ). Notably, bathymetry in the northeast direction from the receiver exhibited range-dependent variations, while bathymetry in the southwest direction displayed minimal changes in water depth. Our current analysis focused on cases where the source moved away from the receiver in the southwest direction within a range-independent environment. It is important to note that when the source-receiver distance was less than 18.6 km, the energy of received signal exceeded the upper dynamic limit of the hydrophone, leading to signal truncation. Moreover, airgun data received at distances exceeding 121.5 km exhibited poor signal-to-noise ratio (SNR). Therefore, our analysis was confined to the range encompassed between these two distances.

Acoustic waves propagating in a waveguide exhibit dispersion as normal modes at different frequency propagate at different group velocities. Therefore, several dispersive modes with different frequency-dependent group velocities can be observed in the spectrogram of the received signal after propagating at least several kilometers in shallow water, and the group velocity v m ( f ) is calculated by Equation 1 (Frisk, 1994; Jensen, 2011; Bonnel et al., 2020):

where, k m ( f ) is the horizontal wavenumber of mode m at frequency f. Therefore, the modal travel time t_m as a function of frequency over propagation range r is given by

Figure 2A shows the time series of the airgun sound received at a source-receiver distance of 18.6 km. Interestingly, arrivals that appear to be precursor arrivals were received before the water-borne arrivals. The precursor arrival is a signal that propagates primarily through the sediment layer and arrives prior to any water-borne arrival (Dahl and Choi, 2006; Choi and Dahl, 2007), and the observation of precursor arrivals in our measurements means that the sound speed in lower sediment layer might be faster than that in the water column. Figure 2B shows the spectrogram for water-borne arrivals corresponding to the red box in Figure 2A . The first two modes were clearly observable in the spectrogram. However, higher modes were not clearly visible.

In this study, the warping method was applied to extract each dispersion curve from the spectrogram. In normal mode theory, for an impulsive source signal propagating in a Pekris waveguide with a rigid bottom, the received signal as a function of time t can be expressed as (Jensen, 2011; Bonnel et al., 2020)

where M is the mode number of propagating modes, a_m is the amplitude of mode m, c_w is the water sound speed, and f_cm represents the cutoff frequency of the m-th mode which can be calculated by ( 2 m − 1 ) c 4 D , where D is the water depth. Note that the phase term of (Equation 3) is a non-linear function for time t. To linearize the non-linearity, we used t 2 + ( r c w ) 2 as the warping function h(t) under the assumption that the impulsive source signal propagates in an isovelocity waveguide. The warped signal y_w (t) is calculated by

where h(t) is the warping function, h ′ ( t ) is the time derivative of h(t) (Bonnel and Chapman, 2011; Bonnel et al., 2020; Liu et al., 2020). Figure 2C shows the result of the warping transformation, in which the first two modes are well separated in the warped frequency.

Now that the two modes have distinct warped frequency bands, we extracted the first and second modes using a bandpass filter, as shown in Figures 3A, B , respectively. Each warped modes were then inversely warped by replacing h(t) with h ( t ) − 1 = t 2 − ( r c w ) 2 in (Equation 4). As a next step, wavelet synchrosqueezing transform was applied to sharpen the time-frequency resolution of the dispersion curves of each mode (Daubechies et al., 2011; Thakur et al., 2013), and then time-frequency ridge tracking algorithm (Meignen et al., 2015; Iatsenko et al., 2016) was used to extract the maximum-energy ridges from wavelet synchrosqueezing transform results, and the results are shown in Figures 3C, D . Finally, the extracted dispersion curves were smoothed out using a 5-point moving average filter, with the results shown as red dashed lines in Figures 3C, D .

The seafloor of the experimental site was flat with a water depth difference of only ~5 m from the receiver position to a point 18.6 km southwest. The sediment structure was reported to be composed of a thin layer of unconsolidated mud overlying the high-density sediment created by the grounding events of ice masses through the repeated advance and retreat of glaciers (Niessen et al., 2013; Dove et al., 2014; Han et al., 2023). The sound speed profile in the water column was measured using a conductivity-temperature-depth (CTD) cast on 14 September at points ~5 and ~56 km from the receiver (see the black and red open circles in Figure 1B ), and the measurements showed that the sound speed profiles at two points were similar and distributed within the range of 1,438 and 1,448 m/s ( Figure 4A ).

In this study, a two-layer bottom model for geoacoustic parameter inversion was constructed based on the previous survey results (Han et al., 2023), as shown in Figure 4B . It was assumed that the ocean environment is range-independent, and the bottom consists of a homogeneous fluid sediment layer overlying a homogeneous fluid half-space. The genetic algorithm was used to estimate geoacoustic parameters by matching the extracted dispersion curves from the acoustic data with the replicas predicted by the KRAKEN acoustic propagation model. As mentioned in Section 2.1, the shortest range from the receiver position where the received signals were not truncated was 18.6 km, and at which point the water depth difference between the receiver and airgun positions was only ~ 5 m. Therefore, we used the data received at this point for geoacoustic parameter inversion.

Based on the two-layer bottom model, we estimated six unknown parameters including the layer thickness H, sound speed c _s and density ρ s in the surficial sediment layer and sound speed c _b and density ρ b of the basement, and the water depth D. The sound speed profile measured at the point ~5 km from the receiver position was used as model input to run the KRAKEN acoustic propagation model. The search spaces for each parameter are shown in Table 1 , which were set to sufficiently wide ranges based on previous survey results (Han et al., 2023). To create model replicas, the KRAKEN model was run in the frequency range between 16 and 130 Hz for mode 1 and 36 and 170 Hz for mode 2 at a 1-Hz interval. Then, the group velocity predictions as a function of frequency were converted to modal travel time for the source-receiver distance of 18.6 km by (Equation 2). To compare the modal travel times between the model replicas and those extracted dispersion curves from the acoustic data, each result was time-aligned by adjusting the arrival time for 130 Hz of the first mode to 0 on the time axis. Then, the least square method was used, in which the objective function to be minimized is given by

where M is the total number of modes used for the geoacoustic inversion, and M was set to 2 in this study. The vector X = { H , c s , ρ s , c b , ρ b , D } . t m ( f n ) and t ^ m ( f n ,  X ) is the measured and modeled modal travel times, respectively. N_m is the total number of frequency segments used for the mode m.

To determine the optimal values of geoacoustic parameters, a global search was conducted using a genetic algorithm over the parameter search spaces. The genetic algorithm parameters were set as follows: a population size of 64, a crossover fraction of 0.8, and a mutation probability of 0.05 to prevent local minima. The algorithm terminated either when the number of generations reached 300 or when the objective function no longer decreased within an additional 60 generations from the generation that achieved the best result.

The best-fit values and their uncertainties for each parameter within the search space are presented in Table 1 . The inversion results reveal a water depth of approximately 68.2 m, with a surficial sediment layer approximately 4.1 m in thickness overlaying a high-velocity basement. The surficial sediment exhibits an approximate sound speed of 1422.4 m/s and a density of about 1.58 g/cm³. These results suggest that the estimated geoacoustic parameters of the surficial sediment layer closely align with our previous survey results (Han et al., 2023). Additionally, the underlying basement displays a sound speed of approximately 1733.6 m/s and a density of approximately 1.84 g/cm³. Figure 5 shows the sensitivities of J ( X ) to parameter variations around the optimal values of the six environmental parameters. A sensitivity plot for each parameter is created by calculating J ( X ) within the search spaces while keeping other parameters at their optimal values. As expected, J ( X ) is most sensitive to the sound speed of the lower sediment layer and water depth D, while the other parameters are relatively less sensitive in the inversion process. In the sensitivity results for the surficial layer thickness, 0 m represents a scenario where the sediment is not structured with two sediment layers but solely with a half-space. Figure 6 compares the measured modal curves with model replicas predicted using the inversion results for the first four modes. Although the mode-3 and 4 are not clearly visible in the measured spectrogram, the predictions of mode-1 and 2 are in good agreement with the corresponding measured modal curves.

As a subsequent step, the Bayesian approach was applied to estimate the uncertainties in geoacoustic parameter estimates derived from the inversion process. Let X and d denote the given vectors representing the six geoacoustic parameters and the measured data, respectively. Following Bayes’ rule, the posterior probability density P ( X| d ) can be expressed as Equation 6 (Gerstoft and Mecklenbräuker, 1998; Dosso, 2002; Dosso and Dettmer, 2011)

where P ( X ) is the prior distribution, and P ( d |X ) is the conditional probability density of given vector X for the measured data d. P ( d ) is the probability density of measured data d, acting as a normalizing factor. To calculate P ( X| d ) , the prior distribution is assumed to be uniform within the search bounds of each parameter. The conditional probability density P ( d |X ) , interpreted as the likelihood function L(X), is given by Equation 7

where J ( X ) represents the data misfit function (considered in section 2.3). Consequently, the posterior probability density can be written as Equation 8

In this study, the Markov-chain Monte Carlo method utilizing Metropolis sampling (Metropolis et al., 1953; Hastings, 1970) was employed to estimate the marginal probability densities for each parameter. The inversion results were then utilized as initial values for the estimation. The sampling process involved 10,000 iterations using a proposal distribution in the form of a normal distribution, centered on the sampled results. Considering the parameter scale, the standard deviation of the proposal distribution for sound speed parameters was set to 1, while for the remaining parameters, it was set to 0.1. Subsequent to estimating the marginal probability densities for each parameter, parameter uncertainties were quantified using 95% Highest Probability Density (HPD) credible intervals (Bonnel et al., 2013; Gelman et al., 2013), as depicted in the last right column of Table 1 .

Up to this point, we have derived estimates for geoacoustic parameters and water depth by comparing measured modal curves with model replicas for a source-receiver distance of 18.6 km. Conversely, if ocean environmental parameters, including geoacoustic parameters, are assumed to be known, the distance between the acoustic source and receiver can be estimated using the same method as described above. As indicated in section 2.1, reliable signals were received within the range of 18.6 and 121.5 km. Therefore, in this session, we estimate source-receiver distances for this range. Initially, assuming a range-independent environment, we executed the KRAKEN acoustic propagation model with input parameters comprising the estimated geoacoustic parameters and a water depth of 68.2 m. Subsequently, model replicas were generated for source-receiver distances ranging from 1 to 140 km at 1-m intervals. These replicas were compared with the measured modal curves extracted using the same method outlined in section 2.2. Figures 7A, B illustrate the comparison between estimated source-receiver distances and those measured by GPS, along with their respective distance errors. Notably, distances estimated under the assumption of a range-independent environment (depicted by black open circles) exhibited increasing errors as the distance extended, in contrast to the actual distances measured by GPS. This discrepancy appears to be attributed to the gradual decrease in water depth from approximately 70 to 47 m as the source-receiver distance increases from 18.6 to 121.5 km. Consequently, the group velocity, as a function of frequency, decreases accordingly.

To mitigate the discrepancy in the range-independent environment, we employ an adiabatic approximation for model propagation in the range-dependent environment (Jensen, 2011; Bonnel et al., 2022). The modeled modal travel times, denoted as t ^ m ( f n ,  X ) in (Equation 5), are calculated by dividing the source-receiver distance into range-independent segments based on a 1-m change in water depth for the bottom bathymetry depicted in Figure 1C . The modal travel times for each segment are then summed using Equation 9

Here N r is the total number of range segments, Δ r i represents the range in the ith range segment, and v m ( r i ,   f n ,  X ) is the group velocity predictions of mode m as a function of frequency for the ith range segment.

Subsequently, the inversion process was reiterated for the source-receiver distance of 18.6 km, and the inversion results are presented in Table 2 . The estimated geoacoustic parameters in the range-dependent environment exhibited consistency with those in the range-independent environment, falling within 95% HPD credible intervals, except for the thickness of the surficial sediment, estimated to be 5.3 m. the revised results indicate a surficial sediment layer with a thickness of approximately 5.3 m. The surficial sediment layer has an approximate sound speed of 1432.6 m/s and a density of about 1.41 g/cm³. The underlying basement exhibits a sound speed of approximately 1737.6 m/s and a density of approximately 1.87 g/cm³. Based on these optimal values, we now re-estimate the source-receiver distances to 121.5 km using the optimal values for the range-dependent environment. Consequently, as the source-receiver distance increased, the distance error−up to 30% in the range-independent environment was reduced to within 10%. It was assumed that the estimated geoacoustic parameters were independent of changes in water depth with distance, possibly explaining why the error could not be further reduced.