Influence of seabed on very low frequency sound recorded during passage of merchant ships on the New England shelf

Since the early 1980s, the acoustic intensity striations observed in time-frequency spectrograms have been extensively studied because they contain information about the acoustic propagation in the ocean environment. This information content is useful for both environmental characterization and source localization. In particular, the waveguide invariant β has been defined as a quantity that correlates directly to the environment.^1–9 This work identifies a feature in the striation pattern from surface ships over the very low frequency (VLF) portion of the spectrum that may have gone unnoticed in previous studies from shallow water environments with multilayered seabeds. Specifically, this VLF feature is observed on merchant ship spectrograms measured during the 2017 Seabed Characterization Experiment¹⁰ in an area called the New England Mudpatch.¹¹ Here, we correlate the characteristics of this feature to the geoacoustic structure of a multilayered seabed using a normal mode methodology. Evidence is given that this VLF feature occurs at frequencies where the group velocities of pairs of modes that contribute to a striation become equal, which coincides where $β→±∞$⁠.

The waveguide invariant has been defined as the slope of the interference striation on a log-log plot of a range-frequency spectrogram. As explained in Ref. 3, by setting a first-order Taylor's series expansion of the acoustic intensity to zero, the waveguide invariant corresponding to the interference between modes m and n is the ratio of the difference in the phase slowness (inverse of the modal phase speed) to the difference in the group slowness (inverse of the model group velocities). This direct dependence on the modal properties of the ocean environment is potentially what makes β a useful tool in environmental characterization of an ocean waveguide.

The remainder of this paper is organized as follows. The feature in ship spectrograms observed on the New England Mudpatch is described in Sec. II. The theoretical formulation that explains the observed feature is discussed in Sec. III. The approach to the inverse problem appears in Sec. IV and Sec. V shows the results of the analyses. The conclusion appears in Sec. VI.

Ship radiated noise data were recorded on two Marine Physical Laboratory (MPL) Scripps vertical line arrays (VLAs) on the New England continental shelf in 2017.¹⁰ The received level versus frequency and time records were processed for the passage of many merchant ships. About twenty ship passage events possessed an adequate signal to noise ratio to allow direct observation of the VLF feature in the spectrograms. An example is shown in Fig. 1. These low-frequency spectrograms come from the container ship KALAMATA, measured on March 24, 2017 on hydrophone 8 (about 40 m deep) on VLA 1 and VLA 2. Spectrograms are shown for the 20–100 Hz band and the smaller 20–40 Hz band. From automatic identification system (AIS) data, the closest point of approach (CPA) is about 3.3 km and 6.0 km for VLA 2 and VLA 1, respectively, and the ship speed was about 20 kn.

To first order, these spectrograms exhibit a nested set of hyperbolas (also known as striations) with the vertices occurring at CPA. The striations are approximately symmetric about the CPA time over a wide frequency band. The symmetry is not exact due in part to a stern/bow aspect dependence of the broadband source level and the fact that the sound propagated upslope to the receivers. In addition to this temporal symmetry, the VLF portion of the spectrograms shows a partial symmetry axis at $f=F$ (shown by the white dashed line in Fig. 1) that separate two hyperbola branches. The main observation is that the slope of all striations, with respect to frequency, vanishes at $f=F$⁠. In Sec. 3 it will be shown that $f=F$ is the frequency where β_mn becomes infinite for a mode pair m and n, or equivalently the modal group velocities for modes m and n become equal.

Generally large merchant ships radiate a portion of their sound at discrete frequencies below 50 Hz. One indeed observes in Fig. 1 that superimposed on the broadband striations are a groups of tonals in the 20–40 Hz band. The tonals appear to be associated with the KALAMATA because one observes the same set of tonals as the ship passes both VLA 1 and VLA 2. If a ship is undergoing uniform motion, there will be a Doppler shift to the lower frequencies. However, as the ship passes both CPAs, the tonals have a positive Doppler shift which is consistent with AIS data that had the KALAMATA undergoing a course and speed (increase) change as it approached both VLAs from the east. A detailed analysis of the effect of the motion of the ships on the received spectra is beyond the scope of the present analysis, and no attempt is made to include motion effects or discrete tonals in the acoustic modeling of the spectrograms.

Table I shows the frequencies at which $f=F$ for twenty merchant ship data samples. $f=F$ varied from about 22.7 to about 26.2 Hz, with a mean value of 24.5 Hz and a standard deviation of about 1.0 Hz. The observed variability may be ascribed to the variability of the characteristics of the deeper sediment layers in the proximity of the New England Mudpatch. The basis of this statement will be discussed in Sec. VI.

This section presents a physical explanation for why $β=±∞$ at $f=F$⁠. The acoustic pressure in an ocean waveguide can be expressed as a linear response

where $G$ is the complex Green's function solution to the Helmholtz equation for a unit point source and S(f) is the source function. $G$ can be expressed as a normal mode expansion

where k_m and $ϕm$ are the horizontal wavenumber eigenvalues and depth dependent eigenfunctions, respectively; r is the source-receiver range; z_s and z are the source and receiver depths, respectively; and $ρ(zs)$ is the ratio of the density of the fluid medium relative to the density of water at the source depth. When the magnitude of the argument of the Hankel function, $H01(kmr)$⁠, is greater than 5, its asymptotic approximation can be used and Eq. (2) becomes

To understand the characteristics of the observed striations in Fig. 1 it is useful to examine the acoustic intensity

where $|Pmn|2$ is the acoustic intensity that results from the interaction of two modes m and n:

where

and

From Eq. (5), the range at which modal interference of modes m and n causes a striation is

Thus, the condition that the slope of a striation vanishes at $f=F$ is

or

which is equivalent to

For $kn≠km$ this condition leads to

where the group velocity for mode n is defined as

The waveguide invariant β_nm for modes n and m can be defined as³ 

where

and the modal phase speed of the nth mode is defined as

Thus, one observes that the condition for the vanishing of the slope of a striation [Eq. (9)] for the non-degenerate case $km≠kn$ is $β=±∞$⁠, which corresponds to $Vm(2πF)=Vn(2πF)$⁠. Finally, $βnm(2πF)=+∞$ for C_m > C_n and $βnm(2πF)=−∞$ for C_m < C_n.

The VLF feature in the ship spectrograms at $f=F$ occurs because the group velocities of two modes are the same. Thus, it is necessary to identify the mode pairs that contribute to the observed striation patterns in the frequency band of interest. The need to identify mode pairs is analogous to the need to identify multipath pairs to estimate relative time delays in localization methods such as that reported by Westwood and Knobles.¹² In this study, preliminary simulations of spectrograms (such as those in Fig. 1) with a normal mode model^13,14 showed that modes 1 and 2 were responsible for the striation pattern in the 20–40 Hz band. As the frequency increases, higher-order modes begin to contribute. Thus, an inversion approach based on β_mn over a larger bandwidth would need to a priori identify multiple pairs of modes that are responsible for the striation patterns. This approach would generally result in a spectrum of $F$ values. Such an involved study is beyond the scope of the present paper. To clearly understand the cause of the VLF feature showing $β=∞$⁠, we concentrate on the 20–40 Hz band where the striations are dominated by modes 1 and 2.

The inversion idea adopted in this study was a matched-feature approach to find optimal values for $F$⁠. In a feature-based inversion, a solution in parameter space H is found that gives the lowest value for an error function E(H). In this work,

where $Fmod(H)$ is the modeled value of the frequency where the group velocities of a selected mode pair become equal.

A critical step of this inversion approach is to select a frequency resolution with which to extract the observed $F$ value from spectograms. The frequency spacing is determined from a choice of an integration time of $Tint=(1/Δf)=2N/fs$⁠, where the sampling rate $fs=$ 25 kHz and N is an integer and $Tint$ is much less than the full observation period of the spectrogram. Generally increasing $Tint$ provides finer resolution of $F$⁠. However, a practical limit on $Tint$ exists because the ship is moving relative to the receiver, and for each instant of time, the propagation path from the source to the receiver changes. Thus, spatial variability places limitations on the size of $Tint$⁠. In such cases, taken to the extreme too large of $Tint$ would smear the spectrograms. Fluctuations in the water column that affect the amount of temporal coherence can impose limitations on $Tint$⁠. Fortunately, during SBCEXP, the spatial and temporal variability in the water column were rather small.

The choice of N = 21 which gives $Tint≈84$ s appeared to be a reasonable choice for the spectrogram data analyzed in this work. This choice gives a frequency resolution of $Δf=0.0119$ Hz and yields 1681 frequency bins in the 20–40 Hz band. Using this $Δf$⁠, the $F$ values that were extracted for the spectrogram data samples shown in Fig. 1 for the merchant ship KALAMATA were 24.9982 Hz and 23.3531 Hz for VLA 2 and VLA 1, respectively. For the application of feature inversion we only utilize these values of $F$ in separate optimizations of selected model geoacoustic parameter values.

The next step of the inversion process was to select a parameter space H. The goal of this work was to understand with the simplest waveguide model possible the physical properties of the seabed that could explain the observed values of $F$ in the New England Mudpatch. The approach adopted was to create an effective horizontally stratified representation of the seabed. For this purpose, a geoacoustic profile was adopted that has four sediment layers over an acoustic basement, as shown in Table II. The model for the first three layers is based in large part on the sub-bottom layering survey and analyses performed by Goff et al.¹⁵ The geoacoustic properties in the first layer are based on an analysis by Buckingham,¹⁶ whose grain shearing model provided a good fit to the reported sound speed ratio (inferred and direct measurements from about 15 separate measurements and analyses).¹⁰ The geoacoustic values for the 2 m transition layer followed by another layer of about 7.5 m of sand are based on a transdimensional inversion analyses of reflection data by Belcourt et al.¹⁷ The transition layer has a very large sound speed gradient, and has been previously explained by Goff.¹⁵ The total thickness of the mud and transition layer is about 12.3 m, and was established using the two-way travel time bottom layering data from Ref. 15 combined with an analysis of Signal, Underwater Sound (SUS) waveform data.¹⁸ The sound speed in the basement was previously determined independently by Wan¹⁹ and Potty and Miller²⁰ in modeling broadband data from SUS where the Airy phase frequency for mode 2 was modeled with a basement sound speed of about 1812 m/s. It is of interest to note that the reported Airy phase frequency for mode 2 lies between 26 and 28 Hz which is larger than the observed $F$⁠. The Airy phase frequency is interpreted as the minimum group velocity for a given mode.

Pre-modeling of the measured striations demonstrated that variations in the sediment structure of the first three layers could not explain the observed $F$ values. It was hypothesized that the observed value of $F$ could be ascribed to the existence of a deep fourth layer not reported by Goff.¹⁵ Assuming nominal parameter values for the density and attenuation of this layer, the inversion parameter space H, thus, consists of two parameters: the sound speed with upper and a lower bounds of 1751 and 1812 m/s, respectively, and the thickness with a lower and upper bound of 100 and 300 m. The lower and upper bounds of the sediment thickness were determined by a sensitivity study.

The inversion procedure adopted in this study is to use Monte Carlo sampling over the parameter space H. The total number of Monte Carlo samples $NMC$ was selected to be 5000. For each parameter sample, the normal mode eigenvalues $kn(2πf)$ were computed for 1681 frequencies in the 20–40 Hz band for a pre-selected number of modes. The criteria of total mode number selection was made based on the manner in which the ORCA normal mode algorithm finds modes in the complex horizontal wavenumber plane.^13,14

For each sampling of parameter values in H, $Fmod$ is found such that $Vm(2πFmod)=Vn(2πFmod)$⁠, and $Fmod$ is one of the 1681 processing frequencies in the 20–40 Hz band. Generally, the result of the optimization calculation yields multiple combinations of the sediment thicknesses and the sound speeds that give an error function value E = 0. For larger $Tint$⁠, a smaller number of solutions result in E = 0. The selection of $Tint≈$ 80 s, which gave 1681 frequencies in the 20–40 Hz band, generally resulted in about 5–10 solutions with E = 0. In this work, we report the average value and standard deviations of the sediment thickness and sound speed in the fourth layer for the E = 0 solutions.

The feature–based inversion scheme was employed with $F=24.9982$ Hz and $F=23.3531$ Hz for the KALAMATA spectograms measured on VLA 2 and VLA 1, respectively. The optimum sound speeds obtained for the fourth sediment layer were 1810.6 ± 3 m/s and 1753.3 m/s ± 3 m/s for VLA 1 and VLA 2, respectively. The optimum sediment thickness obtained for the fourth sediment layer were 270 m ± 30 m and 173 m ± 30 m for VLA 1 and VLA 2, respectively. The uncertainty arises only from the standard deviation of values that produced a value of E = 0. If one were to include values of E > 0 the uncertainty would increase.

Figure 2 shows the group velocity dispersion curves for the optimized geoacoustic profiles for the VLA 1 and VLA 2 values for $F$⁠. The group velocity versus frequency was constructed with a five-point numerical derivative²² of $kn(2πf)$⁠. The critical observation in Fig. 2 is denoted by the blue dashed vertical lines indicating when $V1(2πF)=V2(2πF)$⁠. Finally, one observes that mode 2 has a broad minimum in its group velocity from about 25–28 Hz, which is consistent with the reporting of the Airy phase frequency in Refs. 18–20 that utilized SUS data to infer seabed properties of the New England Mudpatch.

Figures 3 and 4 compare the optimized model spectrograms to the measured spectrograms for the merchant ship KALAMATA on hydrophone 8 (about 42 m depth) of VLA 1 and VLA 2, respectively. The modeled spectrograms are computed using Eqs. (1) and (2) with the optimized sound speed and sediment thickness in the fourth sediment layer. In the mode sum in Eq. (2), an assumed effective source depth of 8.3 m was used, which is the reported draft of the KALAMATA. The CPA ranges and ship speeds were slightly adjusted to provide a qualitatively good fit between the measured and modeled spectrogram. The adjusted CPA ranges for VLA 1 and VLA 2 were 6000 m and 3250 m, respectively. The adjusted ship speed that was used in the computations was 20.9 kn. The modeled and measured spectrograms are in qualitative agreement.

Merchant ship spectograms measured on the New England Mudpatch in 2017 revealed a low-frequency feature attributed to the group velocities of modes 1 and 2 being equal at a frequency $F$⁠. This observations is equivalent to the statement that the waveguide invariant for modes 1 and 2 is infinite at $F$⁠; i.e., $β12(2πF)=+∞$⁠. The observed average value for $F$ for about 20 merchant ship data samples was about 24 Hz. The Airy phase frequency for mode 2 is slightly greater than $F$⁠; the mode 2 group velocity possesses a broad minimum from about 26 to 28 Hz.

It was shown that variations in $F$ on the order of about 1.5 Hz could be ascribed to changes in the sound speed (50 m/s) and thickness of an effective deep sediment layer. This inferred effective layer appears to be consistent with the results of a seismic survey of the New England shelf performed by Siegel et al.²³ in 2014. The seismic study revealed zones of low sediment sound speeds about 250 m below the water sediment interface in an area just east of the Mudpach that extends down to the continental slope. Specifically, Fig. 2 in in Ref. 23 shows a reflector at about 250 m beneath the seafloor and is labeled as a glacial unconformity formed by an ice stream. The two estimates of a deep sediment lie above and below the deepest layer of sediments identified from the survey as overburden pressure. This is consistent with the sound speed only being 1812 m/s. Thus, the existence of a $βnm(2πF)=+∞$ regime may indicate a means for detecting deep sediment properties.

This research was supported by the Office of Naval Research.