Acoustic scattering from layered seafloors exhibits dependence on both the mean geoacoustic layering, as well as the roughness properties of each layer. Several theoretical treatments of this environment exist, including the small roughness perturbation approximation, the Kirchhoff approximation, and three different versions of the small slope approximation. All of these models give different results for the scattering cross section and coherent reflection coefficient, and there is currently no way to distinguish which model is the most correct. In this work, an integral equation for scattering from a layered seafloor with rough interfaces is presented, and compared with small roughness perturbation method, and two of the small slope approximations. It is found that the most recent small slope approximation by Jackson and Olson [J. Acoust. Soc. Am. 147(1), 56–73 (2020)] is the most accurate when the root-mean-square (rms) roughness is large, and some models are in close agreement with each other when the rms roughness is small.

The ocean floor contains variations in both its roughness and layering structure. At low frequencies, sound can interact with sub-bottom layers but the effects of roughness may be relatively small (at least for modest roughness). At high frequencies, the attenuation in the ocean bottom is higher, and the interaction with sub-bottom layers is reduced, but the effect of scattering may be more important. At intermediate frequencies, acoustic waves interact with both the sub-bottom layering and roughness. It is often of practical interest to remotely sense properties of both sub-bottom layers and the rough interfaces that separate them. Several models have been previously developed to solve the forward problem using a point source.1–5 

Models used for these purposes are limited in their applicability. The Kirchhoff approximation (KA) is restricted to angles close to the specular direction, and the small roughness perturbation method (SPM, or perturbation theory) performs best away from specular. The small-slope approximation was introduced by Voronovich,6 and applies to the entire angular range for certain parameters of the rough interface.7 Its original incarnation was for Dirichlet boundary conditions, but it has been applied to fluid, elastic,8,9 and poroelastic10 halfspaces.

Recently, the small-slope approximation has been expanded to encompass layered media, but there are three competing models. One small-slope approximation for layered media was developed by Jackson,11 but is not explored here due to its strange behavior for slow sediment layers. Another small-slope approximation was developed by Gragg and Wurmser12 in the acoustics literature, and later by Berrouk et al.13 in the electromagnetics literature. It is denoted SSL2 in this work. The last and most complicated small slope approximation was developed by Jackson and Olson,14 and is denoted SSL3. The SSLn convention comes from Jackson and Olson14 and is retained here. As shown by Jackson and Olson,14 all of these models disagree for certain roughness and layer geoacoustic properties. This ambiguity is troubling. Although SSL3 has the most physically relevant motivation, it is not clear which approximation should be used in a given situation.

In this work, we remedy this ambiguity by providing a comparison between SSL2, SSL3, SPM, and the exact solution using integral equations. Two geoacoustic environments with two sets of roughness parameters each are used. The integral equation method is based on Monte-Carlo averaging, so individual realizations must be produced. A recent application of the Kirchhoff approximation3 treats the seafloor layering more faithfully than previous work but is specialized to the point-source, point-receiver geometry, not plane waves. Since formally-averaged quantities such as the scattering cross section are not available, we make no comparisons to this model in this work. A comparison to these models is certainly a fruitful area for future work. We find that SSL3 provides the best match with exact results for the scattering strength and coherent reflection coefficient. We do not present a systematic study of the region of validity for these models, although that is also a productive area for future work.

In Sec. II we present the geometry and environment. The basic concepts for the models used here are presented in Sec. III. The integral equation method is detailed in Sec. IV. Comparisons are made to theory in Sec. V. Discussion and conclusions are presented in Sec. VI.

The geometry of the problem is shown in Fig. 1. Although arbitrary fluid layering is treated in theoretical work,14 we limit the problem here to an overlying water column (a half-space), a fluid layer, and an underlying fluid half-space (which we refer to as the basement). These domains are denoted as Ω0,Ω1, and Ω2, respectively. Each domain, Ωn is bounded by one of two boundaries. Γ1 bounds Ω0 from Ω1 and is the water-sediment interface. Γ2 bounds Ω1 from Ω2, and is the interface between the sediment layer and the sediment basement. The boundary of a domain Ωn is denoted Ωn, with Ω0=Γ1,Ω1=Γ1Γ2, and Ω2=Γ2. The normal vectors associated with each of these boundaries are shown in Fig. 1. Note that both normal vectors point into Ω1. This property is important for derivation of the boundary integral equations.

FIG. 1.

(Color online) Layered environment and geometry. Although only the upper interface is depicted as rough, the integral equations defined here can be used with two rough boundaries. Γn denotes each interface, and Ωn denotes the medium immediately above Γn. The arrows show the direction of integration used in the integral equations developed in Sec. IV.

FIG. 1.

(Color online) Layered environment and geometry. Although only the upper interface is depicted as rough, the integral equations defined here can be used with two rough boundaries. Γn denotes each interface, and Ωn denotes the medium immediately above Γn. The arrows show the direction of integration used in the integral equations developed in Sec. IV.

Close modal

Each domain, Ωn, is characterized by a phase speed cn, density ρn, and dimensionless loss parameter, δn. The complex sound speed in each domain can be written as

c̃n=cn1+iδn.
(1)

The wavenumber in each domain is related to the complex sound speed through kn=ω/c̃n, where ω is the acoustic angular frequency with units of radians per second. Dimensionless ratios are defined as ac1=c1/c0, and ac2=c2/c0 for sound speed, and aρ1=ρ1/ρ0 and aρ2=ρ2/ρ0.

The incident acoustic wave vector is specified by

ki=k0k̂i,
(2)

where k̂i is the incident acoustic unit wave vector and is given by

k̂i=cosθix̂sinθiẑ,
(3)

where x̂ is the unit vector in the x (horizontal) direction, and ẑ is the unit vector in the z (vertical) direction. The scattered wave vector into Ω0, back into the water column, is similarly given by

ks=k0k̂s,
(4)

where k̂s is the scattered acoustic unit wave vector and is given by

k̂s=cosθsx̂+sinθsẑ.
(5)

The incident grazing angle θi, and scattered grazing angle θs, are both measured from the horizontal axis.

The rough interfaces are described in terms of their power spectra. Let Wn be the power spectrum of the rough interface constituting Γn. The rough interface Γn is specified by the function fn(xn). The Fourier transform of fn(x) is denoted Fn(kx) with wavenumber argument kx. The power spectrum is defined by Wn(kx1)δ(kx2kx1)=Fn(kx1)Fn(kx2), where the angle brackets denote ensemble averaging. The truncated power law roughness spectrum known as the “von Kármán” spectrum is used here and is specified by

Wn(Kx)=w1n(K0n2+kx2)γ1n/2,
(6)

where w1n is the one-dimensional (1D) spectral strength for interface n with units of m 3γ1n,γ1n is the dimensionless 1D spectral exponent, and K0n is the spectral cutoff for interface n with units of rad/m. The mean square height for interface n is denoted hn2, and is equal to the integral of Wn over the real line. For the von Kármán spectrum,

h1n2=w1nπΓ((γ1n1)/2)K0nγ1n1Γ(γ1n/2).
(7)

Values of γ1n greater than unity result in a finite h1n2. The 1D correlation function for the nth interface, Cn(x), is defined as

Cn(x)=h1n2fn(x)fn(x+x).
(8)

For the von Karman spectrum, the correlation function is

Cn(x)=21νnΓ(νn)(K0nx)νnKνn(K0nx),
(9)

where νn=(γ1n1)/2. Kν(x) is the modified Bessel function of the second kind with argument x and order ν. We assume the rough interfaces do not intersect.

The three scattering models compared in this work are the SPM, SSL2, and SSL3. Since a complete description of these models is quite lengthy, only the elements will be provided here. The reader is referred to Jackson and Olson14 where these models are presented in complete form. We focus on two quantities, the coherent reflection coefficient, and the scattering cross section, both of which are defined in terms of the T-matrix, T(ksx,kix). The T-matrix is a transfer function between an incident plane wave with horizontal wave vector kix, and a scattered plane wave with horizontal wave vector ksx. The scattering cross-section due to 1D roughness, σ(ksx,kix), is defined as15 

σ(ksx,kix)=ksz2k0C(ksx,kix),
(10)
C(ksx,kix)δ(kixkix)=T(ksx,kix)T(ksx,kix)T(ksx,kix)T(ksx,kix),
(11)

where C(ksx,kix) is the incoherent second-moment of the T-matrix. The coherent reflection coefficient, |R(kix)| is defined as

|R(kix)|δ(ksxkix)=|T(ksx,kix)|,
(12)

where the delta function must be included on the left hand side since it is always present in the average T-matrix for stationary roughness. The quantity |R| is actually the magnitude of the complex coherent reflection coefficient, R. However, for brevity, we refer to |R| as the coherent reflection coefficient. The above definitions assume stationary roughness, and we will further assume Gaussian statistics for this random process. The coherent reflection coefficient is frequently used with an angular argument, |R(θi)|, instead of the horizontal component of the wave vector.

For the scattering cross section, all models here use the factor An(ksx,kix) for the nth interface, which is defined as14 

An(ksx,kix)=1ac(n1)2aρ(n1)×An1(ksx)An1(kix)Ãn(ksx,kix),
(13)

where An1(kx) is the amplitude of the downgoing plane wave coefficient in medium n – 1 (just above interface n) due to a plane wave incident from medium Ω0, and

Ãn(ksx,kix)=12{an[1+Vn(kix)][1+Vn(ksx)]bn[1Vn(kix)][1Vn(ksx)]},
(14)

where

an=(1aρ(n1)aρn)ksxkixkn121+ac(n1)2aρ(n1)acn2aρn,
(15)
bn=(aρnaρ(n1)1)βn1(kix)βn1(ksx).
(16)

Vn(kix) is the flat-interface reflection coefficient of interface Γn assuming an overlying infinite halfspace in medium Ωn1. The sine of the angle in Ωn is βn(kx)=1kx2/kn2. For the upper interface, V1(kx) is

V1(kx)=V1H(kx)+V2H(kx)e2ik1β1(kx)D1+V1H(kx)V2H(kx)e2ik1β1(kx)D,
(17)

where D is the mean thickness of Ω1, and VnH(kx) is the reflection coefficient of the nth layer assuming both sides consist of halfspaces defined as

VnH(kx)=Zn1Zn+1,
(18)
Zn=aρnacnβn1(kx)aρ(n1)ac(n1)βn(kx).
(19)

For perturbation theory, the two-dimensional (2D) T-matrix for interface n is

TnSPM(ksx,kix)=ik0β0(ksx)An(ksx,kix)Fn(ksxkix).
(20)

For SSL2, the 2D T-matrix for the layered, rough seafloor is

TSSL2(ksx,kix)=k02πβ0(ksx)Δkz×n=1NAn(ksx,kix)ei(ksxkix)xiΔkzfn(x)dx,
(21)

where Δkz=kszkiz is the difference between the vertical component of the scattered and incident wavenumbers.

SSL3 requires a version of An(ksx,kix) where interface n has been displaced by an amount fn, which is denoted An(ksx,kix,fn). This expression is rather complicated, and the full version is presented in Eqs. (37) and (82) of Jackson and Olson.14 The SSL3 T-matrix for interface n in 2D geometry is

TnSSL3(ksx,kix)=ik02πβ0(ksx)ei(ksxkix)x×0fn(x)An(ksx,kix,f)dfdx.
(22)

The main difference between SSL2 and SSL3 is that in SSL3, the factor An(ksx,ksi,fn) depends on the height of layer n, and the integral over space includes variations in the sediment layering due to roughness fn(x), whereas neither are true for SSL2. Because of this dependence on fn, SSL3 takes into account changes due to roughness in the interference pattern produced by the layered seafloor, whereas SSL2 assumes that the interference pattern is unchanged by roughness. In this way, SSL2 may be thought of as a hybrid between SPM and a true small-slope approximation.

For the scattering cross section and coherent reflection coefficient, we refer the reader to Jackson and Olson.14 Scattering strength for SPM is easy to compute using the formulas provided there. The coherent reflection coefficient for SSL2 is computed from Eq. (74), and scattering strength from Eq. (79) of that reference. For SSL3, the coherent reflection coefficient can be found in Eqs. (83)–(87), and scattering strength in Eq. (90)–(96), and (A1)–(A35). These formulae are omitted due to the large amount of space required to express these approximations and all of their definitions, and interpretations of the models will rely on expressions for the T-matrices presented above.

The Jackson and Olson14 analysis focused on 2D roughness, whereas we consider 1D rough interfaces in this work. These differences are minor for SPM, and formally averaging the SPM result is simple. The 1D version of SSL2 can be found by simply replacing the Kirchhoff integral, Eq. (63) in Jackson and Olson,14 with

In1D(η)=2eη2hn20cos(ΔKxuk0)[eη2hn2Cn(u/k0)1]du,
(23)

where ΔKx=ksxkix is the difference between the horizontal component of the the outgoing and incoming wave vectors. Similarly, the small slope integral, Eq. (89) in Jackson and Olson,14 should be replaced by

Issln1D(ηa,ηb)=e(1/2)(ηaηb)2hn2In1D(ηaηb)
(24)

to compute SSL3.

Integral equations (IE) provide a method to produce the exact scattered pressure due to rough surfaces. Methods for pressure release16 and fluid-fluid17 boundary conditions have been previously presented in the underwater acoustics literature. A numerical method for layered media was presented by Tang and Hefner,18 but its derivation was not based on an integral equation for that environment. Rather, in that reference, a discretized matrix equation is derived from the integral equation for a single interface, and a discretized matrix equation is given for the layered case via physical intuition. The method presented here is based on matching boundary conditions for three different integral equations. This method can be shown to be equivalent to the method of Tang and Hefner after correcting a few errors on the diagonal terms and rearranging the density ratio factors.

The pressure, pn, in any domain n must follow the Helmholtz equation in each domain,

2pn+kn2pn=0,
(25)

where 2 is the Laplacian operator. The Green's function is the solution to the Helmholtz equation with a point source on the right-hand side. In two dimensions, the free space solution (i.e., without boundaries) in domain n using a point source is

2Gn(r,r0)+kn2Gn(r,r0)=δ(xx0)δ(zz0),
(26)
Gn(r,r0)=i4H0(1)(kn|rr0|),
(27)

where i=1 is the imaginary unit, δ(x) is the Dirac delta function, and H0(1)(x) is the Hankel function of the first kind of order zero, with argument x. The position vectors are defined as r=rxx̂+rzẑ and r0=rx0x̂+rz0ẑ. We denote the position vector restricted to Γn as rn, and the normal vector as n̂n.

Within each domain, the pressure field satisfying a Helmholtz equation can be solved using the Helmholtz integral formula, also known as the Helmholtz-Kirchhoff integral equation (HKIE).19 The pressure on the boundary can be expressed, in the absence of an incident pressure field, and with an outward-pointing normal vector (corresponding to the “exterior” boundary value problem), as

α(rl)pn(rl)=Vl,mnpn(rm)nmKl,mnpn(rm),
(28)

where the integral operators are defined as

Vl,mn[ϕ(rl)]=ΓmGn(rl,rm)ϕ(rm)dSm,
(29)
Kl,mn[ϕ(rl)]=ΓmGn(rl,rm)nmϕ(rm)dSm.
(30)

The function ϕ(r) is arbitrary and square-integrable, dSm indicates that the integration is carried out over the boundary with respect to the subscript variable, and /nm=n̂m·m is the normal derivative with respect to the m argument (as opposed to l). The subscript of l, m on the integral operators indicates that integration is carried out along Γm, and the operator output is a function defined on Γl. The parameter α is equal to β/(2π), where β is the angle subtended by the tangent lines on each side of a given point. For a smooth surface, α=1/2 at all points. In this work, we form the integral equation along a piecewise continuous, non-smooth surface, and must calculate α at each point. The operator V is commonly referred to as the single-layer potential operator, and K as the double-layer potential operator. In this work, the exterior form of the Helmholtz integral equation is used with domains having a single boundary—only Ω0 and Ω2.

We may also form the equivalent integral equation for a domain with inward-pointing normal vector (the “interior” boundary value problem). In this case, the interior integral equation is defined for Ω1 only, which is bounded by Γ1 and Γ2. Here, we write the integral operators on each boundary separately, giving

(1α(r1))p1(r1)=V1,11p1(r1)n1+K1,11p1(r1)V1,21p1(r2)n2+K1,21p2(r2),
(31)
(1α(r2))p1(r2)=V2,21p1(r2)n2+K2,21p1(r2)V2,11p1(r2)n2+K2,11p(r1).
(32)

Although written separately, these equations should be thought of as being a single integral equation since the first specifies the pressure on the upper boundary, and the second specifies the pressure on the lower boundary. Both are required for a solution of the boundary value problem.

To form an integral equation for the union of all domains, continuity conditions for pressure and normal velocity are enforced between all domains, keeping track of the normal vector direction. The incident pressure from domain Ω0 is added to the right hand side of the integral equation for that domain, and terms are rearranged to give the integral equations in terms of the following matrix of operators

[α(r1)I+K1,10V1,10(1α(r1))IK1,11aρ1V1,11K1,21aρ1V1,21K2,11aρ1V2,11(1α(r2))IK2,21aρ1V2,21α(r2)I+K2,22aρ2V2,22][p0(r1)p0(r1)n1p1(r2)aρ11p1(r2)n2]=[pi(r1)],
(33)

where a blank spot in a matrix denotes either a zero operator or a variable that is identically zero, and I denotes the identity operator (which maps a function onto itself). The right hand side of this system of integral equations is the incident pressure on Γ1 from Ω0 in the first row, and is zero for all other rows. The unknown variables consist of the pressure in Ω0 and Ω1, as well as their normal derivatives. Note that the normal derivative for Ω1 has the factor aρ11, which is due to the boundary conditions for the continuity of the normal velocity across Γ2.

In this equation, the direction of integration determines the direction of the unit normal vector. The convention followed here is that the normal vector points to the right of the integration direction along each boundary, Γn. In Fig. 1, the direction of integration along each boundary and in each domain has been specified. In Ω1, which has both boundaries, the integration can be thought of as being in the clockwise direction, to the right on the top, and to the left on the bottom. Formally, the integral should be closed in Ω1 between Γ1 and Γ2, but this part of the integral may be neglected if the pressure field decays to zero, which we assume here.

Extensions of this method to multiple layers can be made by formulating the HKIE in each domain and matching boundary conditions. A systematic method to perform this type of calculation was presented by von Petersdorff and Leis,20 although their analysis uses the operators in Eqs. (29) and (30) and their normal derivatives (the adjoint double layer, and hypersingular potential operators, respectively—both of which are not used here). Although the method of von Petersdorff and Leis has superior stability and numerical conditioning than the method used here, it is more complicated due to the hypersingular operator, which is difficult to implement numerically. The numerical condition number for the method detailed in this work has been found to be adequate for our purposes (on the order of 106 or 107).

The incident field used here is an approximation to a plane wave developed by Thorsos.16 This field is incident from Ω0 onto interface Γ1 and takes the form (for our time convention)

pi(r1,f)=pieiki·r1(1+w(r1))(x1z1cotθi)2/g2,
(34)
w(r1)=2(x1z1cotθi)2/g21(k0gsinθi)2,
(35)

where g is a parameter controlling the width of the incident field, and pi is the complex pressure amplitude. The 3 dB angular width of this beam is

Δθ=22log(2)k0gsinθi.
(36)

As the product k0g grows large, the incident field better approximates a plane wave, and it is valid at lower grazing angles. The angular width increases as θi decreases, so small grazing angles are more computationally demanding for the numerical solution of scattering problems.16 

These integral operators can be discretized using standard techniques, such as the boundary element method (BEM).16,21 In this work, these operators were discretized using the collocation method with linear basis functions to approximate the pressure and normal derivative, resulting in a square matrix for each of the integral operators. The matrices were assembled into a fully discrete block matrix according to Eq. (33).

After the pressure and pressure normal derivative on each boundary is found, it is propagated to the far-field using the HKIE. If the pressure in Ω0 is sought, then this becomes

p0(rf)=Vf,10p0(r1)n1Kf,10p0(r1),
(37)

where the subscript f denotes the field pressure point locations. The pressure in other domains can be found from the integral equations for that domain. Once the field pressure is found, the scattering cross section can be estimated by16 

σ=rL|p0(rf)|2|pi|2,
(38)

where r is the distance from the center of the top mean interface to the field point, and

L=gπ2[10.5(1+2cot2θi)(kgsinθi)2]
(39)

is the effective ensonified length of the rough interface. The angle brackets denote ensemble averaging.

The coherent reflection coefficient is a bit more difficult to estimate. Instead of an analytic formulation, we follow the method used by Thorsos.22 We compare the scattered pressure due to the rough layered environment to the scattered pressure in Ω0 due to a flat, rigid boundary of the same length, p0flat. Namely,

|R|=|p0(rf)p0flat(rf)|.
(40)

These calculations use the same tapered incident field.

We present results for two different geoacoustic environments. The first environment has a layer with a greater sound speed than that of water, where ac1=1.05,aρ1=1.8,ac2=1.8,aρ2=2.5. The attenuation parameters are δ1=0.01 and δ2=0.02. We call this environment the “fast layer.” The second environment is a slow mud layer overlying a fast basement with ac1=0.99,aρ1=1.4,ac2=1.8,aρ2=2.5. The attenuation parameters here are set to δ1=0.0005 and δ2=0.02 since softer sediments typically have smaller attenuation coefficient values. This environment is called the “slow layer.” Geoacoustic parameters are summarized in Table I. These geoacoustic properties correspond to the second and third geoacoustic environments presented in Jackson and Olson.14 The acoustic frequency was set to 2 kHz, with ω12.6×103 rad/s.

TABLE I.

Geoacoustic parameters used in numerical examples. All computations use water sound speed c0= 1500 m/s and density ρ0=1000 kg/m3, although only the ratios are important.

CaseDomainThickness (m)Sound speed ratioDensity ratioLoss parameter
Fast 1.05 1.8 0.02 
Layer  1.8 2.5 0.01 
Slow 0.99 1.4 0.005 
Layer  1.8 2.5 0.01 
CaseDomainThickness (m)Sound speed ratioDensity ratioLoss parameter
Fast 1.05 1.8 0.02 
Layer  1.8 2.5 0.01 
Slow 0.99 1.4 0.005 
Layer  1.8 2.5 0.01 

In all results, interface 2 is smooth and interface 1 is rough (except for the integral equation test case). Two sets of roughness parameters for each environment are used. One set has small k0h1, and the other has larger k0h1. These two parameter sets are presented in Table II. Note that w11 and γ11 for the large roughness case are the 1D equivalent to the parameters in the examples presented in Jackson and Olson.14 Formulas found in Appendix D and the errata list from Jackson and Richardson23 were used to perform this conversion. In this table, the root-mean-square (rms) height of each interface multiplied by k0 is shown, as is the rms height divided by the average layer thickness, D, set to 1 m.

TABLE II.

Roughness spectrum parameters used in numerical examples.

Casew11 (m3γ11)γ11K01 (rad/m)k0h1h1/D
Large k0h1 2 ×103 0.66 0.079 
Small k0h1 2 ×104 0.21 0.025 
Casew11 (m3γ11)γ11K01 (rad/m)k0h1h1/D
Large k0h1 2 ×103 0.66 0.079 
Small k0h1 2 ×104 0.21 0.025 

The integral equation results used 48 independent roughness realizations. The incident field width parameter, g, was set to 40λ, which limited the range of grazing angles over which the integral equation results are valid. At 25° grazing, the incident field relative angular width is about 5%, and is 10% at 18° grazing. Therefore, conservatively, results should be trusted above 25 degrees, but plots are shown down to 18 degrees. The total surface length of the realizations was set to L=5g, so that multiple reflections between the interfaces could be captured accurately. This value was chosen by gradually decreasing the value of L until noticeable effects were seen (starting at L=16g). The surface was sampled at Δx=λ/16, which is a rather small sampling interval, but was chosen because coarser sampling did not converge within 1 dB. The rough surfaces are generated using the spectral method of Thorsos16 with the specified sampling interval. However, the power spectrum was low-pass filtered so that slopes at very small scales did not cause numerical issues with the discretized integral equations. The power spectra at wavenumbers between 6k0 and 6.5k0 were smoothly tapered to zero using a raised cosine function (inspired by LePage and Schmidt24), and were set to zero between 6.5k0 and 8k0 (the Nyquist wavenumber). This transition region corresponded to about 100 points of the sampled wavenunber domain, with a total of 3480 points.

With these parameters, SPM was the fastest model to compute, as it required only evaluating functions at the angular and wavenumber arguments. SSL2 was the next fastest, but required evaluating the Kirchhoff integral for each incident and scattered angle. The integration was performed using the trapezoidal rule, and the time required was 0.5 s for each scattering strength figure presented below. SSL3 required computing the Kirchhoff integral many times for each incident and scattered direction, and required 28 min for each scattering strength figure presented below. The integral equation calculations for these parameters took about 12 h for all 48 realizations. Times listed here were for an equivalent single-core processor, but we used a quad-core processor at 4.2 GHz and parallel for-loops to speed up the calculations.

The first result is a validation of the integral equation. The fast layer geoacoustic properties are used with both interfaces flat. The numerically calculated reflection coefficient is compared to the theoretical plane wave reflection coefficient, Eq. (17). These two curves are compared in Fig. 2. The lower limit of the grazing angles is set to 18 degrees since the angular width of the incident field and finite length of the surfaces that enter into the integral equations can cause discrepancies at low grazing angles. The model-IE comparison is quite good, although the IE result is slightly less than the model at small grazing angles by about 1%.

FIG. 2.

(Color online) Flat-interface integral equation result compared with the plane wave reflection coefficient model.

FIG. 2.

(Color online) Flat-interface integral equation result compared with the plane wave reflection coefficient model.

Close modal

A comparison for the coherent reflection coefficient between the theoretical models and integral equations for the fast layer with small roughness is presented in Fig. 3. The small rms roughness causes the integral equation coherent reflection coefficient to depart only slightly from the flat-interface cases. SSL3 agrees with the integral equation result in this case, but SSL2 shows noticeable departures from all other models and integral equations. This figure serves to show that for small values of the rms roughness, defined here as k0h10.2, SSL3 is accurate for the coherent reflection coefficient, but SSL2 is less accurate, although not by much.

FIG. 3.

(Color online) Coherent reflection coefficient for fast layer, small roughness case.

FIG. 3.

(Color online) Coherent reflection coefficient for fast layer, small roughness case.

Close modal

Results for the coherent reflection coefficient for large roughness are presented in Fig. 4. The integral equation departs significantly from the flat interface case, especially near 80 degrees grazing angle, and near the peaks. SSL3 is the best model presented here, although it has some small errors near the peaks. SSL2 follows the integral equation less closely. A notable difference is that SSL2 has a different local minimum near 45 degrees grazing angle than both the integral equation, flat interface, and SSL3. We may conclude that when the roughness is increased, SSL3 is a more physically realistic model, since it matches the integral equation results. Physically, this improved accuracy is due to the fact that SSL3 accounts for changes to the interference pattern when roughness is present in a layered seafloor.

FIG. 4.

(Color online) Coherent reflection coefficient for fast layer, large-roughness case.

FIG. 4.

(Color online) Coherent reflection coefficient for fast layer, large-roughness case.

Close modal

Scattering strength results from the fast layer case with small rms roughness presented in Fig. 5. The integral equation result is shown, along with SPM, and both SSL2 and SSL3. For these values of the geoacoustic and roughness parameters, all three models agree quite well with each other, and the models appear to fall within the uncertainty of the Monte-Carlo simulations. For these values of the roughness and geoacoustic parameters, we conclude that all models examined perform adequately for the scattering cross section. This figure gives confidence that SSL2 and SSL3 agree for small rms roughness, which they should in the limit as k0h10. The agreement of all the models with the integral equation results gives confidence that the implementations of both the theoretical models and integral equations are sound. It is interesting to note that for these roughness and geoacoustic parameters, SSL2 and SSL3 agree for scattering strength, but not for the coherent reflection coefficient.

FIG. 5.

(Color online) Scattering strength results for the fast layer, small roughness case.

FIG. 5.

(Color online) Scattering strength results for the fast layer, small roughness case.

Close modal

When the roughness is increased in Fig. 6, all the models depart from one another, but only slightly. Perturbation theory becomes inaccurate near the specular direction, which is expected, but also contains some small errors at moderate angles. SSL2 performs better than perturbation theory near the specular direction, but still contains moderate errors compared to the integral equation at moderate angles. SSL3 performs better than SSL2, notably near 55 degrees grazing. Another notable difference between SSL2 and SSL3 are the changes in shape near 43 degrees grazing. Close examination of this region shows that SSL3 has a much different shape than both SSL2 and SPM. It appears that SSL3 is predicting an alteration of the interference pattern compared to SPM. The uncertainty of the integral equation results is too large to make a determination about which model is correct, but it appears that SSL3 follows the IE curve more closely between about 45 and 50 degrees grazing.

FIG. 6.

(Color online) Backscattering strength comparison for fast layer, large roughness case.

FIG. 6.

(Color online) Backscattering strength comparison for fast layer, large roughness case.

Close modal

The coherent reflection coefficient for the slow layer geoacoustic environment with small roughness properties is shown in Fig. 7. It is compared with SPM, SSL2, and SSL3. The oscillations in the coefficient are small compared to that of the fast layer. For this case, the integral equation, flat interface, and SSL3 agree quite well with each other. SSL2 contains small discrepancies compared to the integral equation, although less than the fast layer case with small roughness. All the models presented here perform an adequate job for this small roughness case—even the assumption of a flat interface.

FIG. 7.

(Color online) Coherent reflection coefficient for slow layer, small roughness case.

FIG. 7.

(Color online) Coherent reflection coefficient for slow layer, small roughness case.

Close modal

The coherent reflection coefficient for the slow layer with large roughness is plotted in Fig. 8. SSL3 is the best model here with an error of less than 1%. SSL2 departs significantly from the integral equation result, especially near 60–70 degrees grazing and around 35 degrees grazing. SSL2 and the flat-interface model all have errors of less than 10%. Based on these small errors, we conclude that the roughness present for this case has a small effect on the reflection coefficient for a slow layer, but is best modeled by SSL3. However, this conclusion cannot be extended to the scattered field in general, as will be seen in the next example.

FIG. 8.

(Color online) Coherent reflection coefficient for slow layer, large roughness case.

FIG. 8.

(Color online) Coherent reflection coefficient for slow layer, large roughness case.

Close modal

Scattering strength from the slow layer with small roughness is presented in Fig. 9, and the IE results are compared to the same models. The shape of the IE curve is much different from the fast layer case, and there are some deep nulls at several angles. SSL3 follows the IE curve the best. SSL2 and SPM agree with the integral equation result over most of the angular range, except for the three local minima present, near 35, 40, and 60 degrees grazing. Here, SPM and SSL2 underestimate the scattering cross section. SSL3 provides the correct fit near these grazing angles. If the attenuation of Ω1 increased slightly, then SSL2 and SPM would match the IE result nearly as well as SSL3.

FIG. 9.

(Color online) Backscattering strength comparison for slow layer, small roughness case.

FIG. 9.

(Color online) Backscattering strength comparison for slow layer, small roughness case.

Close modal

Scattering strength results from the slow layer case with larger rms roughness are presented in Fig. 10. The integral equation result is shown along with SPM and both SSL2 and SSL3. For these values of the geoacoustic and roughness parameters, SSL2 and SPM agree over all but the largest grazing angles. SSL3 disagrees with both SSL2 and SPM close to specular, and near the local minima away from the specular direction. The integral equation result agrees quite well with SSL3, but not SSL2 or SPM. The large differences between SSL2 and SSL3 are surprising given the similar results presented for the coherent reflection coefficient in Fig. 8. However, we note that the scattering cross section depends on the factors [1±V1(θi)] and [1±V1(θs)], as seen in Eq. (13). Since the magnitude of V1(θ) is close to unity for the slow layer, small changes in V1(θ) can lead to large relative changes in [1±V1(θ)], depending on the sign of V1(θ). We can conclude from this example that, for a slow layer, increasing the roughness causes large discrepancies between SSL2 and SSL3, and also SSL3 is the most accurate model for these parameters to within about 0.3 dB.

FIG. 10.

(Color online) Backscattering strength comparison for slow layer, large roughness case.

FIG. 10.

(Color online) Backscattering strength comparison for slow layer, large roughness case.

Close modal

In Sec. V, we have seen that for the larger values of spectral strength, SSL3 departs from SSL2, which was shown previously in Jackson and Olson.14 We have also seen that the Monte-Carlo integral equation method agrees very well with SSL3 in the large-roughness cases, and agrees with all the models in the small roughness case. In Jackson and Olson,14 it was postulated that the disagreement between SSL2 and SSL3 was due to the fact that SSL3 accounts for changes to the interference pattern due to changes in layer thickness caused by the rough interfaces. The agreement between SSL3 and the integral equation supports the conclusion that this effect is indeed present in scattering from layered surfaces. It remains to be seen whether these differences can be seen in field experiments.

To investigate what is causing these large differences, another numerical experiment was performed with the same sound-speed and density as the fast layer, but with the attenuation of the slow layer and the large roughness parameters. These results are not shown, but the decreased attenuation of the fast layer showed deep nulls in SSL2 and SPM that were not apparent in SSL3 or the IE results. When large roughness is present, its effect is more pronounced if the attenuation is small for both a fast and a slow layer. We can conclude that the presence of roughness changes the interference pattern caused by the layering structure, and this effect is more pronounced with small values of the attenuation coefficient in Ω1.

The confirmation of the differences between SSL3 and SPM has implications for geoacoustic inversion of layered rough seafloors. The alteration of the interference pattern due to roughness could cause inversion schemes that use the reflection coefficient25,26 or scattering strength27 to provide incorrect results if an inappropriate model is used, such as the flat-interface assumption for the reflection coefficient, or the SPM for the scattering cross section. SSL3 provides a promising model to use in such cases.

This work was focused on deciding between competing models for layered rough interfaces for a few sets of geoacoustic and roughness parameters. A systematic study of the validity of each of these scattering models was not performed but would be a valuable avenue for future work. The integral equation methods presented here could be used for such a study.

Funding was provided by the U.S. Office of Naval Research.

1.
D.
Tang
and
D.
Jackson
, “
A time-domain model for seafloor scattering
,”
J. Acoust. Soc. Am.
142
(
5
),
2968
2978
(
2017
).
2.
D. R.
Olson
and
C. W.
Holland
, “
Fast computation of time-domain scattering by an inhomogeneous stratified seafloor
,”
J. Acoust. Soc. Am.
147
(
1
),
191
204
(
2020
).
3.
S.
Pinson
,
J.
Cordioli
, and
L.
Guillon
, “
Spherical wave reflection in layered media with rough interfaces: Three-dimensional modeling
,”
J. Acoust. Soc. Am.
140
(
2
),
1108
1115
(
2016
).
4.
S.
Pinson
, “
Spherical wave scattering from rough surfaces and array processing: Application to sound-speed profile measurement uncertainty analysis
,”
J. Acoust. Soc. Am.
142
(
3
),
1189
1198
(
2017
).
5.
C. W.
Holland
,
S.
Pinson
,
C. M.
Smith
,
P. C.
Hines
,
D. R.
Olson
,
S. E.
Dosso
, and
J.
Dettmer
, “
Seabed structure inferences from TREX13 reflection measurements
,”
IEEE J. Oceanic Eng.
42
(
2
),
268
288
(
2017
).
6.
A. G.
Voronovich
, “
Small-slope approximation in wave scattering by rough surfaces
,”
Sov. Phys. JETP
62
,
65
70
(
1985
).
7.
S. L.
Broschat
and
E. I.
Thorsos
, “
An investigation of the small slope approximation for scattering from rough surfaces. Part II. numerical studies
,”
J. Acoust. Soc. Am.
101
,
2615
2625
(
1997
).
8.
T.
Yang
and
S. L.
Broschat
, “
Acoustic scattering from a fluid–elastic-solid interface using the small slope approximation
,”
J. Acoust. Soc. Am.
96
(
3
),
1796
1804
(
1994
).
9.
R. F.
Gragg
,
D.
Wurmser
, and
R. C.
Gauss
, “
Small-slope scattering from rough elastic ocean floors: General theory and computational algorithm
,”
J. Acoust. Soc. Am.
110
(
6
),
2878
2901
(
2001
).
10.
T.
Yang
,
S.
Broschat
, and
C.
Galea
, “
A comparison of perturbation theory and the small-slope approximation for acoustic scattering from a rough interface for a Biot medium
,”
IEEE J. Oceanic Eng.
27
(
3
),
403
412
(
2002
).
11.
D.
Jackson
, “
The small-slope approximation for layered seabeds
,”
Proc. Mtgs. Acoust.
19
(
1
),
070001
(
2013
).
12.
R. F.
Gragg
and
D.
Wurmser
, “
Scattering from a rough seafloor with stratification
,” in
Proceedings 2005 Conference on Boundary Influences in High Frequency, Shallow Water Acoustics
, Bath UK (September 5–9,
2005
).
13.
A.
Berrouk
,
R.
Dusséaux
, and
S.
Afifi
, “
Electromagnetic wave scattering from rough layered interfaces: Analysis with the small perturbation method and the small slope approximation
,”
Prog. Electromagn. Res. B
57
,
177
190
(
2014
).
14.
D.
Jackson
and
D. R.
Olson
, “
The small-slope approximation for layered, fluid seafloors
,”
J. Acoust. Soc. Am.
147
(
1
),
56
73
(
2020
).
15.
E. I.
Thorsos
and
D. R.
Jackson
, “
The validity of the perturbation approximation for rough surface scattering using a Gaussian roughness spectrum
,”
J. Acoust. Soc. Am.
86
,
261
277
(
1989
).
16.
E. I.
Thorsos
, “
The validity of the Kirchhoff approximation for rough surface scattering using a Gaussian roughness spectrum
,”
J. Acoust. Soc. Am.
83
,
78
92
(
1988
).
17.
E. I.
Thorsos
,
D. R.
Jackson
, and
K. L.
Williams
, “
Modeling of subcritical penetration into sediments due to interface roughness
,”
J. Acoust. Soc. Am.
107
,
263
277
(
2000
).
18.
D.
Tang
and
B. T.
Hefner
, “
Modeling interface roughness scattering in a layered seabed for normal-incident chirp sonar signals
,”
J. Acoust. Soc. Am.
131
(
4
),
EL302
EL308
(
2012
).
19.
A. D.
Pierce
,
Acoustics: An Introduction to its Physical Principles and Applications
(
The Acoustical Society of America
,
Melville, NY
,
1994
).
20.
T.
von Petersdorff
and
R.
Leis
, “
Boundary integral equations for mixed Dirichlet, Neumann and transmission problems
,”
Math. Methods Appl. Sci.
11
(
2
),
185
213
(
1989
).
21.
T. W.
Wu
,
Boundary Element Acoustics: Fundamentals and Computer Codes
(
WIT Press
,
Southampton, UK
,
2000
).
22.
E. I.
Thorsos
, “
Acoustic scattering from a ‘Pierson-Moskowitz’ sea surface
,”
J. Acoust. Soc. Am.
88
(
1
),
335
349
(
1990
).
23.
D. R.
Jackson
and
M. D.
Richardson
,
High-Frequency Seafloor Acoustics
, 1st ed. (
Springer
,
New York
,
2007
).
24.
K. D.
LePage
and
H.
Schmidt
, “
Spectral integral representations of monostatic backscattering from three-dimensional distributions of sediment volume inhomogeneities
,”
J. Acoust. Soc. Am.
113
(
2
),
789
799
(
2003
).
25.
J.
Dettmer
,
S. E.
Dosso
, and
C. W.
Holland
, “
Trans-dimensional geoacoustic inversion
,”
J. Acoust. Soc. Am.
128
(
6
),
3393
3405
(
2010
).
26.
C. W.
Holland
,
P. L.
Nielsen
,
J.
Dettmer
, and
S.
Dosso
, “
Resolving meso-scale seabed variability using reflection measurements from an autonomous underwater vehicle
,”
J. Acoust. Soc. Am.
131
(
2
),
1066
1078
(
2012
).
27.
G.
Steininger
,
C. W.
Holland
,
S. E.
Dosso
, and
J.
Dettmer
, “
Seabed roughness parameters from joint backscatter and reflection inversion at the Malta plateau
,”
J. Acoust. Soc. Am.
134
(
3
),
1833
1842
(
2013
).