Physics-based characterization of soft marine sediments using vector sensors^a) Free

Theoretical considerations,^1,2 laboratory measurements,³ and results of numerous field experiments^4–16 indicate that shear wave speed in granular materials and, in particular, in unconsolidated marine sediments increases with depth z below the seafloor and is approximately proportional to a certain power z^ν of the depth as long as the composition of the materials remains unchanged. The power-law exponent ν is probably controlled by the shape and roughness of the grains. The gradient of the shear wave speed (or shear speed, for brevity) is very large at small z, and the shear speed experiences large relative changes over several meters or 10s of meters below the seafloor. Relative changes in density and compressional wave speed are much smaller, and these geoacoustic parameters can be modeled as depth-independent in a surficial layer of constant composition. Then, power-law depth-dependence of shear speed corresponds to the same power-law dependence on overburden pressure. Surficial unconsolidated sediments are “soft” in the sense that their shear rigidity and shear speed are small compared to the bulk modulus and compressional speed, respectively. For a more detailed discussion of the power-law depth-dependence of shear rigidity and additional references, see Refs. 2, 3, 17, and 18.

Soft sediments with power-law shear velocity profiles support horizontally (SH) and vertically polarized (SV) interface waves, which propagate along the seafloor with phase and group speeds on the order of the shear speed.^10,17 These interface waves are slow in the sense that their phase and group speeds are small compared to the sound speed in water and compressional speed in the bottom. The vertically polarized seismo-acoustic interface waves are usually referred to as Scholte waves.^19–22 The dispersion and polarization properties of slow Scholte waves supported by soft sediments, shape functions of these waves, and wave energy distribution between water and the seabed are all quite different from those of the Scholte waves that are supported by the interface of homogeneous fluid and solid half-spaces.²³ Moreover, dispersion properties of the vertically and horizontally polarized slow interface waves prove to be very similar,^10,17 making vector sensors indispensable for identifying the wave types. The distinctive feature of the slow interface waves, which is readily recognized in their measured dispersion curves, is a power-law dependence of their phase and group speeds on frequency. There is a one-to-one correspondence between the exponents of the power laws for the frequency dependence of phase or group speeds and the depth-dependence of the shear speed.^10,17 Observations of the interface waves are of considerable interest because their dispersion allows one to characterize geotechnical and geoacoustic parameters of surficial sediments that are difficult to measure by other means.^{7,18–20,24,25}

Vector sensors are increasingly employed in underwater acoustics to characterize seabed properties.^26–28 A rich dataset on wave propagation in the seabed^29,30 was obtained in 2007 in the course of shear wave surveying of the Gjøa oil/gas condensate field in the North Sea off Norway, where a seabed-coupled mechanical vibrator generated probing signals in the frequency band from a few to 60 Hz. A long, densely populated linear array of three-component vector sensors was employed, which helped to separate vertically and horizontally polarized waves, identify a number of interface waves, and measure their phase speeds (Fig. 1). Measured dispersion curves of the interface waves have been inverted to retrieve the shear speed profile in the upper 45–50 m of the seabed.^22,30

There are two striking features of the vector sensor data, which have not been previously explored. First, the vertical and radial components of the measured particle velocity have sharp peaks around 38 Hz [Fig. 1(a)], which suggest some kind of a seismo-acoustic resonance.^9,31,32 Second, when plotted on the log-log scale, the dispersion curves of the interface waves exhibit two distinct slopes at large and small phase speeds [Fig. 1(b)], which suggests that the seabed contains layers with two different power-law profiles of the shear wave speed.^7,10,17 In this paper, we re-examine the experimental results reported by Dong et al.²² with the goal of developing a simple, parsimonious geoacoustic model that qualitatively explains and quantitatively reproduces the key features of the observations. Our approach to the inverse problem is guided by a theoretical analysis of seismo-acoustic resonances and interface wave dispersion in soft sediments containing two or more layers of different composition.

The remainder of the paper is organized as follows. The experimental data underlying this work is described in Sec. II. Approximate analytic dispersion relations of interface waves supported by the seabed, which consists of two continuously stratified soft sediment layers overlaying a solid, homogeneous subbottom, are derived in Secs. III A and III B. The Wentzel–Kramers–Brillouin (WKB) approximation is employed in the derivation. The analytic dispersion relations are used in Sec. III C to find a simple geoacoustic model consistent with the interface wave observations. A physical mechanism of resonant reflection of compressional (P) waves by the seabed and geoacoustic implications of the observed resonant reflection are investigated in Sec. IV. The resulting geoacoustic model is compared to alternative models in Sec. V. Section VI summarizes our findings.

The data analyzed in this paper were acquired in a 2007 shear-wave survey^29,30 of the Gjøa field located in the Norwegian Channel in the northern North Sea off the southern coast of Norway. The water depth at the experiment site was 364 m, and the main geological interfaces at the site are flat. Surficial sediment layers are composed of soft Holocene clays deposited on glacial and glacio-marine sediments.^29,30 A massive seabed-coupled vibrator generated the seismo-acoustic wave field. The wave source was developed by the Norwegian Geotechnical Institute to efficiently generate low-frequency shear waves of different polarizations; limited P waves were also radiated by the source.^29,30 The frequency content of the probing signals generated by the source was approximately from 2 to 60 Hz with a broad maximum around 37 Hz and width of about 20 Hz at half-power level (see Fig. 5 in Ref. 29).

The signals were received on a 1-km-long ocean-bottom cable (OBC), which was deployed partially in water and partially on the seafloor in a radial direction from the source. The OBC contained 42 three-component accelerometers with 25-m spacing. To improve the resolution of short waves, a 600-m-long synthetic aperture with a much shorter 2.5-m receiver spacing was created by dragging the cable in 2.5-m steps.²⁹ Orientations of the three orthogonal receiver components were determined using airgun signals and used to represent the data in terms of the vertical and in-line (radial) and cross-range (tangential) horizontal components. This proved critical for proper discrimination and identification of various arrivals within the complex full field data.^22,29,30 Assuming a horizontally stratified seabed, the cross-range particle velocity is due to SH shear waves, while radial and vertical components of the particle velocity are due to SV shear waves and P waves. Detected arrivals included head waves, multiply reflected shear waves, and at least 10 modes of horizontally and vertically polarized interface, or surface, waves.^22,29,30

Interface waves were observed at frequencies from about 2 to 20 Hz. Dispersion curves of the SH interface waves have been extracted from the cross-range components of particle acceleration measured on the synthesized aperture horizontal array, while dispersion curves of the SV interface waves have been measured using the vertical and radial components of the acceleration.^22,30 The dispersion curves are illustrated in Fig. 1(b). The interface wave dispersion curves have been previously inverted by Socco et al.³⁰ and Dong et al.²² to retrieve the depth dependence of the shear wave speed in the top 40–50 m of the seabed. The seabed was modeled as a stack of homogeneous layers in these inversions.

Because of limitations on access to proprietary raw data, this paper focuses on re-analysis of the previously published^22,29,30 information on interface wave dispersion and power spectra of signals recorded by the three-component vector sensors. Available data consists of the frequency dependence of the phase speed of various interface waves [Fig. 1(b)], as retrieved by Dong et al.,²² and power spectra of the vertical, radial, and cross-range components of the full field. The power spectra²² averaged over multiple receivers and repeatedly emitted probing signals are shown in Fig. 1(a). For each of the vertical, radial, and cross-range components of particle velocity, the average power spectra are normalized by their respective maxima.

The main maxima of the power spectrum of the cross-range component of the field are at frequencies below 20 Hz [Fig. 1(a)]. In addition to broad low-frequency peaks below 10 Hz, which are associated with SV interface waves, the power spectra of the vertical and in-line components have significantly larger, narrow peaks around 38 Hz. (A much smaller peak at a similar frequency in the spectrum of the cross-range component is probably due to imperfect separation of the measured acceleration into the vertical, radial, and cross-range components resulting from uncertainties in the measurements of spatial orientation of individual vector sensors.) These sharp peaks are suggestive of a resonance phenomenon occurring in either the experimental equipment or the environment. In particular, as already mentioned, the source spectrum is maximum at about 37 Hz. However, the bandwidth of the source spectrum at half-power is at least 20 times larger than the sub-1-Hz width of the spectral peaks of the wave field [Fig. 1(a)]. We interpret the sharp spectral peaks around 38 Hz as a seismo-acoustic resonance originating from wave propagation conditions at the experimental site. It is shown in Secs. IV and V that such an interpretation is consistent with available geological information and results of inversion of the interface wave data.

Consider a model of soft marine sediments (Fig. 2), which consists of two layers with power-law shear velocity profiles:

The layers are located between the water column at z < 0 and a homogeneous solid half-space (subbottom) at z > H. Here h is the thickness of the upper sediment layer, and H is the vertical extent of the soft sediments. Physical considerations and available observations indicate that $0≤ν1, 2<1$⁠.^9,10,17 Shear and P wave speeds and density in the subbottom are c_sb, c_lb, and ρ_b, respectively. Sound speed and density of water near the seafloor are c_w and ρ_w; P wave speeds and densities in respective sediment layers are c_l1, ρ₁ and c_l2, ρ₂. For simplicity, we assume that variations of the sediment density and P wave speed are negligible within each soft sediment layer. We will also assume that shear speed increases steadily with depth, which implies $a1hν1≤a2(h+z0)ν2$ and $a2(H+z0)ν2≤csb.$

The increase of the shear speed c_s with depth below the seafloor creates a waveguide for shear waves. SH interface waves are normal modes of this waveguide. Despite the simplicity of the geoacoustic model, the wave equation cannot be solved analytically in terms of known mathematical functions for arbitrary values of exponents ν₁ and ν₂.^9,17 We will use a WKB-based asymptotic approach to derive the dispersion relation of the interface waves. Disregarding reflection at the interface z = h, the normal mode dispersion equation can be written as follows in the WKB approximation:²³ 

Here ω stands for wave frequency, V₁ and V₂ are plane-wave reflection coefficients at the upper, z = 0, and lower, z = z_lb, boundaries of the waveguide. The lower boundary is either the turning point z = z_t, where shear speed equals the phase speed u of the normal mode: c_s(z_t) = u, or the lower boundary z = H of the soft sediment, if there are no turning points. Note that the phase integral steadily increases with u.

Introducing a new integration variable, $w=u2/cs2−1$⁠, reduces the phase integral φ(z) in any layer with a power-law dependence of c_s to

This is a standard integral [see, e.g., Eq. (1.2.4.3) in Ref. 33], which can be expressed in terms of a hypergeometric function³⁴ for arbitrary integration limits but simplifies when one of the limits is either w = 0 or infinity. Note that w = 0 at the turning point z = z_t, and w → +∞ when z → 0.

All normal modes are evanescent waves in the subbottom and have phase speeds u < c_sb. When $0<u<a1hν1,$ the turning point z = z_t of the wave is located in the upper sediment layer at $zt=(u/a1)1/ν1.$ Then, integration is over the semi-infinite interval 0 < w < +∞ in the phase integral, and we obtain

in agreement with Ref. 17. Here Γ(⋅) is Gamma function (see Chap. 6 in Ref. 34).

When $a1hν1≤u≤a2(h+z0)ν2,$ integration in the phase integral is from z = 0 to z = h. The latter corresponds to a finite value of w. Using Eq. (1.2.4.3) in Ref. 33, we find

Here F(A, B; C; D) is the hypergeometric function, also known as the Gauss hypergeometric series or ₂F₁(A, B; C; D) hypergeometric function (see Chap. 15 in Ref. 34).

When $a2(h+z0)ν2<u<a2(H+z0)ν2,$ the wave has a turning point at $zt=(u/a2)1/ν2−z0$ within the lower sediment layer. Then, the phase integral is a sum of the integral in the upper sediment layer, which is given by Eq. (5), and an integral over h < z < z_t in the lower sediment layer. Similar to derivation of Eq. (5), we obtain

Finally, when $a2(H+z0)ν2≤u<cbs,$ there are no turning points, and the phase integral is given by

In the WKB approximation, the reflection coefficient from the turning point equals $V2=exp (−iπ/2).$²³ The reflection coefficient from the boundary z = 0, where c_s vanishes and the shear speed gradient becomes infinite, has been found in Refs. 9 and 17 and equals

for SH waves. Using these reflection coefficients V₁ and V₂ from the dispersion Eq. (3), we find the frequency of the SH interface wave with a turning point in one of the sediment layers:

Here n =1, 2, … is the order of the interface wave. Dependence of the interface wave frequency on the phase speed enters Eq. (9) via φ(z_t). Higher-order interface waves (normal modes) have higher frequencies at the same value of the phase velocity u and higher phase speeds at the same value of frequency. Explicit expressions for the phase integral in Eq. (9) are given by Eqs. (4) and (6) when the turning point is located in the upper or lower sediment layer, respectively.

When there are no turning points and $a2(H+z0)ν2≤u<cbs,$ the wave is reflected from the boundary z = H. The plane wave reflection coefficient of SH waves²³ at this boundary is

From the dispersion Eq. (3), we find

where the phase integral is given by Eq. (7). Finally, when $a1hν1≤u≤a2(h+z0)ν2,$ reflection occurs at z = h. The result is similar to Eq. (11) and differs by replacement of φ(H) with φ(h), Eq. (5). In addition, in the expression for the phase of the reflection coefficient in Eq. (10), one should use elastic parameters in the vicinity of the boundary z = h and replace ρ₂ with ρ₁, ρ_b with ρ₂, c_sb with $a2(h+z0)ν2,$ and $a2(H+z0)ν2$ with $a1hν1$(see Fig. 2).

In the special case, where a₁ = a₂, ν₁ = ν₂, z₀ = 0, and ν₁ → 0 in Eqs. (1) and (2), we have a homogeneous solid layer with the shear speed c_s = a₁ that is located between homogeneous fluid (z < 0) and solid (z > H) half-spaces. In this limit, our problem reduces to the textbook setting for Love interface waves.³⁵ The resulting waveguide for SH waves is also equivalent to the acoustic waveguide in a homogeneous fluid layer between a rigid boundary at z = 0 and a homogeneous fluid half-space z > H.²³ In the limit ν₁ → 0, Eq. (8) gives the correct result V₁ = 1 for the reflection coefficient of SH waves at the solid-fluid interface,²³ and Eq. (3) gives $φ(H)=Ha1−2−u−2$ for the phase integral. An inspection shows that the interface wave frequencies f_n that are predicted by Eq. (11) with ν₁ = 0 agree with the textbook result³⁵ for the Love wave dispersion in this special case.

Equations (4) and (9) show that frequency f_n of n-th interface wave is proportional to $u1−1/ν1$ when the turning point is located in the upper sediment layer. On the logarithmic scale, the slope of the dispersion curve, $d(ln fn)/d(ln u)=1−ν1−1$⁠, depends only on the shear-speed power-law exponent in Eq. (1).

When the phase speed u is much larger than the shear speed around z = h, the turning point is located deep in the lower sediment layer, and the vicinity of the turning point gives the main contribution into the phase integral in Eq. (3). Indeed, it follows from Eqs. (5), (6), and the equation³⁴ 

that under these conditions φ(z_t) is given approximately by Eq. (4) with a₁ and ν₁ replaced with a₂ and ν₂, respectively. Then, the slope of the dispersion curves $d(ln fn)/d(ln u)=1−ν2−1$ is controlled by the shear-speed power-law exponent in Eq. (2).

The dispersion equations, which are derived for SH interface waves in this section and for P–SV waves in Sec. II B, describe a gradual transition between the limiting cases of the constant slope of the dispersion curves.

Unlike SH shear waves, SV shear waves are coupled to P waves by the shear-speed gradients. In the case of the power-law shear velocity profile, the coupling is particularly strong near the seafloor z = 0.¹⁷ P–SV coupling leads to the appearance of two types of slow interface waves that are supported by soft marine sediments, the fundamental mode and the main sequence modes.^10,17 The main sequence modes are uncoupled from the water column, just like SH interface waves. In the WKB approximation, the dispersion Eq. (3) of the main sequence modes differs from that for SH waves by having a different reflection coefficient¹⁷ V₁ from the boundary z = 0 [cf. Eq. (8)]:

SV reflection coefficient at interfaces, where parameters of the solid are discontinuous, is also different from the reflection coefficient Eq. (10) of SH waves. In particular, the SV reflection coefficient from the boundary z = H can be written as $V2=exp (−2iΦSV),$ where

C = c_sb, M = ρ_b/ρ₂, and $N=a2−1(H+z0)−ν2csb.$C, N, and M have the meaning of the shear speed below the boundary, the ratio of the shear speeds just below and just above the boundary, and the ratio of densities above and below the boundary, respectively. Equation (14) has been obtained from the general equation for the plane wave reflection coefficient of SV waves at solid-solid interface [see, e.g., Eq. (4.2.9) in Ref. 23] in the limit when $cs/cl→0$ in both solids.

Solving the dispersion Eq. (3) for the main sequence modes with appropriate reflection coefficients V₁ and V₂, we obtain

for the waves with a turning point in one of the sediment layers. Here, as in Eq. (9) for SH modes, the phase integral is given by Eq. (4), when $0<u<a1hν1,$ and by Eq. (6), when $a2(h+z0)ν2<u<a2(H+z0)ν2$⁠. When $a2(H+z0)ν2≤u<cbs,$ waves are reflected from the boundary z = H, and we find

from Eqs. (3), (13), and (14). The phase integral in Eq. (16) is given by Eq. (7). Finally, when $a1hν1≤u≤a2(h+z0)ν2,$ waves are reflected at z = h. The result in this case differs from Eq. (16) by substitution of φ(h), Eq. (5), for φ(H). In addition, $C=a2(h+z0)ν2,$M = ρ₂/ρ₁, and $N=a1−1h−ν1a2(h+z0)ν2$ in Eq. (14) for this boundary.

The accuracy of the WKB-based asymptotic dispersion equations increases with increasing mode order,¹⁷ and the results may not be reliable at n = 1. In addition, the WKB approximation gives discontinuous results and is not accurate when turning points approach and cross interfaces, where elastic parameters are discontinuous, i.e., in the vicinity of $u=a1hν1,$ $u=a2(h+z0)ν2,$ and $u=a2(H+z0)ν2.$

An alternative approach to approximating the dispersion equation, which is particularly useful for low-order modes, was developed in Ref. 17. The approach takes advantage of the availability of an exact solution, when the power-law exponent ν₁ = 0.5, and builds a perturbation theory with respect to the parameter |ν₁ − 0.5| that is assumed to be small compared to unity. In marine sediments, |ν₁ − 0.5| < 0.5 and is often rather small. When the shear speed in soft sediments follows the power law, by neglecting terms of second and higher order in |ν₁ − 0.5|, the dispersion equation of the main sequence of P–SV interface waves can be written as¹⁷ 

for arbitrary n = 1, 2, …. Under the same assumptions, the dispersion equation of the fundamental mode is¹⁷ 

Here γ = 0.577 21… is the Euler's constant, and ψ stands for the Digamma function.³⁴ The counterpart of Eq. (17) for SH waves is¹⁷ 

As discussed in Ref. 17, Eqs. (17) and (19) can be used for interface waves in the case of a multi-layered seabed, provided the turning point is located in the upper soft sediment layer with a power-law shear speed profile. Equations (17) and (19) complement the asymptotic dispersion Eqs. (9) and (15) for the low-order, low-speed modes, for which the WKB-based results are either unavailable or not reliable.

We employ the analytical dispersion relations obtained in Secs. III A and III B as the forward model to match the measured values (Sec. II) of phase speeds of SH and SV interface waves. A nonlinear least-squares method is used to fit all the data for both wave types, simultaneously. Data from the fundamental P–SV mode and the lowest order (n = 1) SH mode are fit to the dispersion curve for the one-layer model, i.e., Eqs. (18) and (19), respectively. Data for the higher-order modes are fit to the asymptotic dispersion relations, Eqs. (9) and (15), for the two-layer model. Simultaneously fitting the data for all interface waves to multiple theoretical dispersion curves reduces the goodness of fit for any one dispersion curve, but it ensures consistency between sediment parameters estimated across all the curves.

It is assumed in the inversion that z₀ = 0 in Eq. (2) and that all modes have turning points above the bottom z = H of the second sediment layer. Then, the geoacoustic model contains six unknown parameters: depth h of the boundary between sediment layers, the density ratio ρ_w /ρ₁, and the power-law parameters a₁, ν₁, a₂, ν₂ in Eqs. (1) and (2). Results of the inversion, including 95% confidence bounds of the estimated parameters, are shown in Table I. The estimated value of the density ratio ρ_w /ρ₁= 0.537 in Table I corresponds to the density ρ₁ = 1910 kg/m³ in the top 5.6 m of the seabed.

These parameters are used to generate a dispersion curve for each P–SV and SH mode, which are drawn as solid lines in Figs. 3(a) and 3(b) for comparison with the experimental data. A dotted line marks the maximum phase speed with turning points in the first layer, $u=a1hν1,$ and a dashed line marks the minimum phase speed with turning points in the second layer, $u=a2hν2.$ All but one of the data points for the fundamental (n = 0) P–SV and the first SH modes lie below these lines, justifying the use of the single-layer model for them. The dispersion curve for the mode n = 1 in the main sequence of P–SV modes is matched with larger errors than the other modes ostensibly because the WKB approximation becomes more accurate as mode number n increases.

Line 1 in Fig. 3(c) shows the shear speed profile as a function of depth using the parameters from Table I and Eqs. (1) and (2). Line 2 is the multi-layer model from Dong et al.²² As noted in that paper, a single power-law profile is not a good fit for the data. Our two-layer model is a better fit for the data and is in reasonable agreement with the multi-layer inversion result, as discussed in more detail in Sec. V. The maximum phase speed in the data set, 350 m/s, produces the greatest turning depth, 42.5 m. These data cannot be used to estimate shear speeds at depths greater than this.

In this section we investigate the hypothesis that the strong, narrow peaks in the observed power spectra of vertical and radial components of particle velocity [Fig. 1(a)] result from the propagation conditions of P–SV waves at the site of the experiment. We offer a physical interpretation of these observations as resulting from resonantly enhanced reflection from the layered seabed, relate the resonance to the shear speed inversion results, and discuss the geoacoustic information contained in the peak frequency f_p = 38 Hz.

Seismo-acoustic resonances are often observed when surficial marine sediments have low shear speeds, but at much lower frequencies between about 0.3 and 7.5 Hz (see, e.g., Refs. 9, 31, and 32). Those resonances arise due to reflection of shear waves and, unlike the results illustrated in Fig. 1(a), are characterized by a large ratio of horizontal-to-vertical particle velocity amplitudes and do not exhibit a large difference between amplitudes of two orthogonal horizontal components of the particle acceleration.⁹ In the North Sea experiment discussed in this paper, the peak occurs at the frequency that is considerably larger than the frequencies of observed surface waves and is, therefore, likely to be caused by P waves. The travel time 1/f corresponding to the peak frequency is smaller than acoustic travel time from the source on the seafloor to the ocean surface. Thus, any interference phenomena or resonances responsible for the observed peak should be explained in terms of the ocean bottom properties.

Geoacoustic inversion of the measured dispersion curves of interface waves (Sec. III C) reveals a boundary between sediment layers at h ≈ 5.6 m below the seafloor. Shear speeds just above and just below the boundary are approximately 76 and 83 m/s, respectively, which are much smaller than the P wave speeds c_l in the sediments. Surficial sediments at the experimental site are described as soft Holocene clays.^29,30 For such sediments, c_l is expected to be somewhat less than the sound speed in water near the bottom, c_w, and increase with the depth below seafloor.^7,18,36

We will show that the power spectrum peak can be explained by the interference of P waves reflected from the seafloor and the boundary z = h within sediments. Consider first a simplified geoacoustic model, where shear rigidity is neglected at z < h, i.e., the top layer of the bottom is approximated by a homogeneous fluid with sound speed c_l1 [Fig. 4(a)]. The ocean bottom at z > h is modeled as a homogeneous solid half-space with P wave speed c_l2 and shear wave speed c_s2. The reflection coefficient of a plane acoustic wave incident from water on the seafloor will be infinite when the following condition²³ is met:

Equation (20) is similar to Eq. (3) but refers to P waves, and reflection coefficients V₁ and V₂ have a different meaning. Here V₁ and V₂ are plane-wave reflection coefficients at z = 0 and z = h for sound waves in the layer 0 < z < h. As in Eq. (3), V₁ and V₂ are the reflection coefficients for incidence from below and from above, respectively. In Eq. (20) u has the meaning of the phase speed of the trace of sound waves on a horizontal plane; in terms of u and wave frequency ω, the horizontal component of the wave vector ξ = ω/u. Equation (20) coincides with the dispersion equation of acoustic normal modes with phase speed u in the waveguide formed by the layer 0 < z < h.

For propagating (as opposed to evanescent) plane waves in the layer, the absolute values of reflection coefficients V₁ and V₂ do not exceed unity. For the condition in Eq. (20) to be met, |V₁| and |V₂| should be equal to 1 simultaneously. The reflection coefficient of a plane sound wave in fluid from a solid half-space is²³ 

Here Z_s and Z_l are impedances of shear and P waves at z > h; Z is the impedance of P waves at 0 < z < h; and θ_s is the angle that wave vector of the shear wave, below the interface, makes with the normal to the interface z = h:

For a propagating P wave incident on a solid half-space with a shear speed smaller than compressional speed c_l1, θ_s and impedances Z and Z_s are real and positive according to Eq. (22). Then, it follows from Eq. (21) that |V₂| < 1 unless u = c_l2. When u = c_l2, impedance Z_l is infinite, and V₂ = 1. This property of the reflection coefficient has a simple physical meaning. Acoustic waves cannot be totally reflected from the solid half-space because a part of the incident energy is carried away from the boundary by shear waves in the solid. The only exception occurs when the impedance of the refracted P wave in the solid becomes infinite at u = c_l2, and the amplitude of the shear wave vanishes.²³ 

The condition |V₁| = 1 will be satisfied at u = c_l2 provided

This inequality ensures that the plane wave is totally reflected at the fluid-fluid interface z = 0. The reflection coefficient from the top boundary of the layer, for incidence from below, is

at total internal reflection.²³ Hence, the resonance condition in Eq. (20) will be met at frequencies f_l,j that satisfy the following equation:

The above derivation of the resonance conditions in Eqs. (23) and (25) extends an earlier discussion by Duncan et al.³⁷ of frequencies with sharply reduced transmission losses in an underwater waveguide with a homogeneous solid bottom, when the sound speed in water is larger than the shear wave speed and smaller than the P wave speed in the bottom. The fluid-fluid boundary at z = 0 in our problem reduces to a pressure release boundary in the limit ρ_w → 0. In this limiting case, the arctangent in Eq. (25) is replaced with π/2, and our result reduces to that of Ref. 37. When ρ_w → 0, |V₁| = 1 at all incidence angles and for any c_w, and the requirement c_l2 < c_w in Eq. (23) does not apply.

The lowest-frequency P wave resonance corresponds to j = 0 in Eq. (25) and occurs at the frequency

Subsequent resonances are equally spaced in frequency with the spacing

Note that the frequency difference f_{l, j+ 1} – f_{l, j} > c_l1/2h. Under the conditions of the North Sea experiment, where h ≈ 5.6 m, the frequency spacing exceeds 85 Hz for all reasonable values of c_l1 > 1000 m/s, and—in agreement with the observations²²—only one resonance, f_l,0, is observed within the 2–60 Hz frequency band of the source.

With the resonance frequency f_l,0, layer thickness h, and sound speed in water known, Eq. (26) relates three geoacoustic parameters: P wave speeds c_l1 and c_l2 in two sediment layers and the ratio ρ_w /ρ₁ of water and sediment layer densities [Fig. 4(b)]. The value ρ_w /ρ₁ = 0.537 has been obtained from the interface wave data (Table I). If c_l2 were retrieved from, say, measured travel times of compressional head wave data,^38,39 c_l1 could be unambiguously determined from Eq. (26), and vice versa. In the North Sea experiment, the nondimensional parameter f_l,0 h/c_w ≈ 0.14 is small. Then, Eq. (26) provides a strong constraint on deviations of the ratios c₁₁/c_w and especially c_l2/c_w from unity [Fig. 4(b)]. The findings that c₁₁ and c₁₂ are smaller than but close to the sound speed in water are consistent with the available geologic information about surficial sediments³⁰ and expectations for P wave speeds in soft clays.^7,18,36

In the above discussion, we modeled the top sediment layer 0 < z < h as a fluid. To justify the application of the fluid-solid model to the interface z = h between sediment layers, it should be noted first that the layer thickness h = 5.57 m is small compared to the P wave wavelength c_l1/f_p ∼ 40 m. For P waves, the upper layer will act as a homogeneous layer with some effective (averaged) parameters. Given the very fast relative variations of the shear rigidity with depth and that shear rigidity is extremely small in the upper part of the layer, the effective shear speed will be much smaller than the 73 m/s shear speed just above the boundary z = h. Similarly, the shear modulus increases by the factor of ∼20 over the first 40 m below the boundary (see Table I). In a homogeneous half-space model of the sediments at z > h, the effective shear speed should be considerably larger than the 85 m/s value just below the interface as given by the geoacoustic inversion of the interface wave data. Hence, reflection of P waves from the boundary z = h should be treated as reflection at a solid-solid interface with a large contrast in shear speeds.

Figure 4(c) illustrates the angular dependence of the reflection coefficient V₂ of a plane P wave from the interface of two homogeneous solids with a large contrast between shear speeds (c_s2 $≫$ c_s1). The wave is incident from the solid with a smaller shear and compressional speeds (c_l2 > c_l1). Incidence angle θ_l of the wave is related to the trace velocity u by the equation sinθ_l = c_l1/u. The reflection coefficient is calculated using Eqs. (4.2.8), (4.2.13)–(4.2.15) in Ref. 23. The equations are exact but cumbersome and will not be reproduced here. V₂ is real-valued at 0 ≤ u ≤ c_l2 and positive at u = c_l2. Note that |V₂| is relatively small at steep and moderate incidence angles and, just like reflection coefficient Eq. (21) from a fluid-solid interface, has a sharp maximum at u = c_l2 [Figs. 4(c) and 4(d)]. The value of |V₂(u = c_l2)| is close to unity, and the sharp local maximum of |V₂| leads to resonance reflection of P waves from the layer 0 < z < h at the frequencies satisfying Eq. (25) as in the case of reflection from a fluid layer between fluid and solid half-spaces. In this model, the sharpness of the observed resonance peaks [see Fig. 1(a)] is related to the sharpness of the angular dependence of the reflection coefficient around its local maximum at u = c_l2 in Fig. 4(d).

When the layer 0 < z < h has small but finite shear rigidity, the reflection coefficient V₁ from the upper boundary z = 0 of the layer deviates from the reflection coefficient Eq. (24) at a fluid-fluid interface. The reflection coefficient of P waves in a solid at the solid-fluid interface is

(see, e.g., Eq. (4.2.37) in Ref. 23). The reflection coefficient is similar to that of the plane wave incident on the interface from the fluid side, Eq. (21). At boundary z = 0,

in Eq. (28). When shear speed c_s1 is small, θ_s and Z_s are proportional to the small parameter c_s1/u$≪$ 1. When c_l1 ≤ u ≤ c_w, Z is purely imaginary, and it follows from Eq. (28) that |V₁| = 1 up to terms of the third order in c_s1/u; the phase of the reflection coefficient differs from its value in Eq. (24) (i.e., at c_s1 = 0) by terms O((c_s1/u)²). Thus, deviations of V₁ from Eq. (24) are negligible.

Together with the above analysis of V₂ [Figs. 4(c) and 4(d)], these findings justify application of the resonance conditions Eqs. (23) and (25) in our problem.

Identification of the fundamental mode of P–SV interface waves as the only mode that is sensitive to sediment density has allowed us to retrieve an estimate ρ_w/ρ₁ = 0.537 of the density contrast between water and the sediment layer 0 < z < h. In previous geoacoustic inversions^22,30 of the same data set, density was not retrieved. In the two density models postulated in Ref. 30 on the basis of the available geologic information at the experimental site, the density ratio ρ_w /ρ₁ = 0.574–0.583, if the average of density in the upper 6 m of the sediments is taken for ρ₁. These values are close to the value retrieved in Sec. II C and are within the uncertainty interval of that estimate, see Table I.

Similarly, depth-independent P wave speed c_l = c_w in the seabed was postulated in Ref. 30. In Ref. 22, interface wave dispersion curves were found to be insensitive to the compressional speed, which was also assumed to be depth-independent. The relatively small deviations of c_l1 and c_l2 from c_w that are derived in Sec. IV from the measured frequency of the P wave resonance, are consistent with the rough depth-independent models.^22,30 Furthermore, the power spectrum data provides strong constraints on variations of the P wave speed across the seafloor and within top sediment layers [Fig. 4(b)].

Inversion of the interface wave dispersion data is accomplished in Sec. II C by representing the upper 40–50 m of the seabed by two layers with power-law profiles of the shear speed. The model is motivated by the observation of two distinct slopes in log-log representation of the dispersion curves [Fig. 1(b)]. To assess this shear-speed model, it is compared here to several alternative geoacoustic models of soft marine sediments. We have considered three additional models of the shear speed depth dependence: single power-law layer, three power-law layers, and two power-law layers on top of a homogeneous half-space. In the last two models, $cs(z)=a2zν2$ at h < z < H. Below the bottom of the second layer, at z > H, $cs(z)=a3zν3$ in the three-layer model; in the two-layer plus half-space model, the shear speed and density in the half-space are $csb=Na2Hν2$ and ρ_b = Mρ₂. Parameters M and N have the same meaning as in Eq. (14).

In the single-layer model, we have used the more accurate theoretical dispersion Eqs. (17)–(19) for all modes. In conjunction with the other models, Eqs. (18) and (19) have been used for the fundamental (n = 0) P–SV mode and SH mode 1, implying that those modes interact only with the uppermost layer; the WKB approximation has been used for all other modes. The P–SV mode 1 data is not well described by the WKB approximation and therefore does not have a good fit for any model. It might have been useful to exclude that data from the fit, but that has not been attempted.

Results of the interface wave data inversion in the alternative geoacoustic models are summarized in Table II and illustrated in Fig. 5. Ninety-five percent confidence bounds are given in Table II for parameters of the retrieved power-law dependencies. The two-layer model [Figs. 3(a) and 3(b)] shows major improvement over the one-layer model [Figs. 5(a) and 5(b)] in fitting the data. This is reflected in the R² values for the inversions, which increase from 0.966 for the one-layer model to 0.980 for the two-layer model. The difference in the R² values represents a decrease of the model-data misfit variance by the factor of 1.7 in the two-layer model. Comparison of Figs. 5(a)–5(b) and 3(a)–3(b) demonstrates that the one-layer model fails to fit the data at phase speeds below 75–80 m/s. The data-model mismatch is so big [Figs. 5(a)–(b)] that R² values calculated for the fundamental P-SV mode, –1.10, and the first SH mode, –3.07, prove to be negative. In contrast, the two-layer model adequately approximates the low-order mode data, with R² of 0.966 and 0.926 for the fundamental P-SV mode and the first SH mode, respectively.

The physics behind the difficulties that the one-layer model has with low-order modes can be traced back to the fact that dispersion of a slow interface wave is most sensitive to the shear speeds at depths around the turning point (Sec. III A). Parameters of the optimum one-layer model are primarily controlled by properties of the second layer (z > h), where turning points are located for most modes in the dataset. At phase speeds below 76 m/s, the turning points are located in the top layer, 0 < z < h, and the mismatch between the data and one-layer model reflects the difference between the parameters of the two sediment layers.

The two-layer plus half-space model had the same R² and produced identical estimated values of parameters of the layers and extremely close confidence intervals of these parameters (Table II) as the two-layer model (Table I), suggesting that the data does not contain the wave frequencies and mode orders that interacted with the seabed below the bottom of the second power-law layer. Despite an increase in the number of degrees of freedom, the three-layer model does not noticeably improve the dispersion data fit (R² = 0.981) and shows very low sensitivity to parameters of the deepest layer, as reflected in the confidence intervals for H, ν₃, and especially a₃. We conclude that the two-layer model is in the best agreement with available dispersion data.

We have also considered a more general two-layer model, where non-zero values of the parameter z₀ in Eq. (2) are allowed, and z₀ is considered as an additional unknown geoacoustic parameter. Despite an increase in the number of degrees of freedom, no noticeable improvement in the model-data fit was found compared to the original two-layer geoacoustic model in Table I.

A Bayesian multi-layer shear-speed inversion of the interface wave dispersion data was developed by Dong et al.²² and considered as an approximation to the linear shear speed profile in a layer overlying a homogeneous half-space. The multi-layer model²² ensures an excellent fit to the measured dispersion curves but its interpretation as an approximation to a linear profile is questionable. Sediments with linear (ν₁ = 1) profile, unlike power-law profiles with 0 < ν₁ < 1, support neither SH nor slow P–SV interface waves.^9,17 This can be traced back to the fact that, when ν₁ ≥ 1, shear speed decreases so fast near z = 0 that shear wave travel time to the seafloor becomes infinite, the waves experience extraordinary attenuation and never reach the seafloor.

The results of the multi-layer shear-speed inversion²² are compared to results of various simple, power-law-based inversions in Figs. 3(c) and 5(c). (Inversion results are extended to the 60-m depth below the seafloor, as in Ref. 22, although these may be only supported by data up to about the 45-m depth.) The results of power-law inversions, except the single-layer inversion, do not deviate far from the multi-layer geoacoustic model in the top 50 m of the seabed. The two-layer, two-layer plus half-space, and three-layer models are all well within the 95% confidence intervals²² of the Bayesian multi-layer inversions for SH and P–SV waves. Thus, the three simple models and particularly the physics-guided, parsimonious two-layer inversion provide a shear-speed depth dependence, which is arguably as consistent with the data as the much more sophisticated and computationally intensive multi-parameter, multi-layer Bayesian inversion.

Soft surficial sediments support a rich set of slow interface waves, which can account for the bulk of seismo-acoustic energy near the seafloor at low frequencies (between about 1 and a few 10s of Hertz) and are sensitive to the magnitude and depth-dependence of shear rigidity. Hydrophone measurements miss most of the interface waves. Vector sensors, such as tri-axial, bottom coupled accelerometers, are necessary to capture, separate different polarizations, and identify various interface wave modes and other components of the full wave field.

The linear dependence between logarithms of the phase (or group) speeds of the interface waves and their frequency was proposed by Chapman and Godin^10,17 as means to identify a seabed with a power-law shear-speed profile and determine its parameters. In this paper, that simple, physics-based approach to geoacoustic inversions is extended to seabeds containing several layers of soft sediments of different composition. In application to interface wave dispersion data obtained in the North Sea off Norway, the approach leads to a low-parameter model of the shear speed profile as power-law dependences in two layers. The model provides a good fit to the data and agrees with the results of a much more elaborate Bayesian inversion.²² In addition, a boundary between soft sediment layers is detected and sediment density is evaluated, with the result being consistent with available geologic information.

We identified a physical mechanism, which can lead to P wave resonances in stratified soft sediments, and demonstrated that the proposed mechanism can explain sharp peaks of the observed power spectra of the vertical and radial components of the particle velocity. The P wave resonance with a high quality factor is made possible by the fact that amplitudes of converted shear waves, which would otherwise take energy from and attenuate the P wave at reflection from a fluid-solid or solid-solid interface, are strongly suppressed at a particular incidence angle.

The simple, physics-guided approach presented in this paper results in a geoacoustic model that offers a consistent interpretation and a quantitative description of various salient features of the available data of the 2007 shear-wave experiment in the North Sea.

Statoil acquired and provided the data used in this research. This work was supported in part by the Office of Naval Research, Award Nos. N00014-17WX00773 and N00014-20WX01312. Helpful discussions with S. E. Dosso, N. R. Chapman, and A. D. Pierce are gratefully acknowledged.