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.

## I. INTRODUCTION

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 poroelastic^{10} 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 Wurmser^{12} 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 SSL*n* convention comes from Jackson and Olson^{14} 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 approximation^{3} 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.

## II. GEOMETRY AND ENVIRONMENT

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 $\Omega 0,\u2009\Omega 1$, and $\Omega 2$, respectively. Each domain, Ω_{n} is bounded by one of two boundaries. $\Gamma 1$ bounds $\Omega 0$ from $\Omega 1$ and is the water-sediment interface. $\Gamma 2$ bounds $\Omega 1$ from $\Omega 2$, and is the interface between the sediment layer and the sediment basement. The boundary of a domain Ω_{n} is denoted $\u2202\Omega n$, with $\u2202\Omega 0=\Gamma 1,\u2009\u2202\Omega 1=\Gamma 1\u222a\Gamma 2$, and $\u2202\Omega 2=\Gamma 2$. The normal vectors associated with each of these boundaries are shown in Fig. 1. Note that both normal vectors point into $\Omega 1$. This property is important for derivation of the boundary integral equations.

Each domain, Ω_{n}, is characterized by a phase speed *c _{n}*, density

*ρ*, and dimensionless loss parameter,

_{n}*δ*. The complex sound speed in each domain can be written as

_{n}The wavenumber in each domain is related to the complex sound speed through $kn=\omega /c\u0303n$, 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\rho 1=\rho 1/\rho 0$ and $a\rho 2=\rho 2/\rho 0$.

The incident acoustic wave vector is specified by

where $k\u0302i$ is the incident acoustic unit wave vector and is given by

where $x\u0302$ is the unit vector in the *x* (horizontal) direction, and $z\u0302$ is the unit vector in the *z* (vertical) direction. The scattered wave vector into $\Omega 0$, back into the water column, is similarly given by

where $k\u0302s$ is the scattered acoustic unit wave vector and is given by

The incident grazing angle *θ _{i}*, and scattered grazing angle

*θ*, are both measured from the horizontal axis.

_{s}The rough interfaces are described in terms of their power spectra. Let *W _{n}* be the power spectrum of the rough interface constituting $\Gamma n$. The rough interface $\Gamma n$ is specified by the function $fn(xn)$. The Fourier transform of $fn(x)$ is denoted $Fn(kx)$ with wavenumber argument

*k*. The power spectrum is defined by $Wn(kx1)\delta (kx2\u2212kx1)=\u27e8Fn(kx1)Fn(kx2)\u2217\u27e9$, 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

_{x}where $w1n$ is the one-dimensional (1D) spectral strength for interface *n* with units of $m 3\u2212\gamma 1n,\u2009\gamma 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 *W _{n}* over the real line. For the von Kármán spectrum,

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

For the von Karman spectrum, the correlation function is

where $\nu n=(\gamma 1n\u22121)/2$. $K\nu (x)$ is the modified Bessel function of the second kind with argument *x* and order *ν*. We assume the rough interfaces do not intersect.

## III. MODELS FOR SCATTERING FROM 1D ROUGHNESS

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 Olson^{14} 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 *k _{ix}*, and a scattered plane wave with horizontal wave vector

*k*. The scattering cross-section due to 1D roughness, $\sigma (ksx,kix)$, is defined as

_{sx}^{15}

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

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(\theta 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 *n*th interface, which is defined as^{14}

where $An\u22121(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 $\Omega 0$, and

where

$Vn(kix)$ is the flat-interface reflection coefficient of interface $\Gamma n$ assuming an overlying infinite halfspace in medium $\Omega n\u22121$. The sine of the angle in Ω_{n} is $\beta n(kx)=1\u2212kx2/kn2$. For the upper interface, $V1(kx)$ is

where *D* is the mean thickness of $\Omega 1$, and $VnH(kx)$ is the reflection coefficient of the *n*th layer assuming both sides consist of halfspaces defined as

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

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

where $\Delta kz=ksz\u2212kiz$ 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 *f _{n}*, 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

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 *f _{n}*, 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 Olson^{14} 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

where $\Delta Kx=ksx\u2212kix$ 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

to compute SSL3.

## IV. INTEGRAL EQUATIONS

Integral equations (IE) provide a method to produce the exact scattered pressure due to rough surfaces. Methods for pressure release^{16} and fluid-fluid^{17} 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, *p _{n}*, in any domain

*n*must follow the Helmholtz equation in each domain,

where $\u22072$ 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

where $i=\u22121$ is the imaginary unit, $\delta (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\u0302+rzz\u0302$ and $r0=rx0x\u0302+rz0z\u0302$. We denote the position vector restricted to $\Gamma n$ as $rn$, and the normal vector as $n\u0302n$.

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

where the integral operators are defined as

The function $\varphi (r)$ is arbitrary and square-integrable, *dS _{m}* indicates that the integration is carried out over the boundary with respect to the subscript variable, and $\u2202/\u2202nm=n\u0302m\xb7\u2207m$ 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 $\Gamma m$, and the operator output is a function defined on $\Gamma l$. The parameter

*α*is equal to $\beta /(2\pi )$, where

*β*is the angle subtended by the tangent lines on each side of a given point. For a smooth surface, $\alpha =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 $\Omega 0$ and $\Omega 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 $\Omega 1$ only, which is bounded by $\Gamma 1$ and $\Gamma 2$. Here, we write the integral operators on each boundary separately, giving

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 $\Omega 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

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 $\Gamma 1$ from $\Omega 0$ in the first row, and is zero for all other rows. The unknown variables consist of the pressure in $\Omega 0$ and $\Omega 1$, as well as their normal derivatives. Note that the normal derivative for $\Omega 1$ has the factor $a\rho 1\u22121$, which is due to the boundary conditions for the continuity of the normal velocity across $\Gamma 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, $\Gamma n$. In Fig. 1, the direction of integration along each boundary and in each domain has been specified. In $\Omega 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 $\Omega 1$ between $\Gamma 1$ and $\Gamma 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 10^{6} or 10^{7}).

The incident field used here is an approximation to a plane wave developed by Thorsos.^{16} This field is incident from $\Omega 0$ onto interface $\Gamma 1$ and takes the form (for our time convention)

where *g* is a parameter controlling the width of the incident field, and *p _{i}* is the complex pressure amplitude. The 3 dB angular width of this beam is

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 $\Omega 0$ is sought, then this becomes

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 by^{16}

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

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 $\Omega 0$ due to a flat, rigid boundary of the same length, $p0flat$. Namely,

These calculations use the same tapered incident field.

## V. RESULTS

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,\u2009a\rho 1=1.8,\u2009ac2=1.8,\u2009a\rho 2=2.5$. The attenuation parameters are $\delta 1=0.01$ and $\delta 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,\u2009a\rho 1=1.4,\u2009ac2=1.8,a\rho 2=2.5$. The attenuation parameters here are set to $\delta 1=0.0005$ and $\delta 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 $\omega \u224812.6\xd7103$ rad/s.

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 *w*_{11} 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 Richardson^{23} were used to perform this conversion. In this table, the root-mean-square (rms) height of each interface multiplied by *k*_{0} is shown, as is the rms height divided by the average layer thickness, *D*, set to 1 m.

The integral equation results used 48 independent roughness realizations. The incident field width parameter, *g*, was set to $40\lambda $, 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 $\Delta x=\lambda /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 Thorsos^{16} 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 Schmidt^{24}), 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.

### A. Integral equation validation

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%.

### B. Fast layer

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 $k0h1\u22480.2$, SSL3 is accurate for the coherent reflection coefficient, but SSL2 is less accurate, although not by much.

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.

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 $k0h1\u21920$. 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.

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.

### C. Slow layer

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.

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.

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 $\Omega 1$ increased slightly, then SSL2 and SPM would match the IE result nearly as well as SSL3.

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\xb1V1(\theta i)]$ and $[1\xb1V1(\theta s)]$, as seen in Eq. (13). Since the magnitude of $V1(\theta )$ is close to unity for the slow layer, small changes in $V1(\theta )$ can lead to large relative changes in $[1\xb1V1(\theta )]$, depending on the sign of $V1(\theta )$. 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.

## VI. DISCUSSION AND CONCLUSION

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 $\Omega 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 coefficient^{25,26} or scattering strength^{27} 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.

## ACKNOWLEDGMENTS

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