## ABSTRACT

The small-slope approximation (SSA) for rough-interface scattering is most commonly applied to the upper boundary of either impenetrable media or uniform half-space media, but has been recently developed for layered media in the acoustic and electromagnetic cases. The present work gives an overview of three forms of the SSA for layered media. The first has been previously presented in the acoustics literature. The second is from the electromagnetics literature and in the present work is converted to the fluid-sediment problem. A missing proof is supplied of a key consistency condition demanded of the small-slope ansatz. As is usual, these small-slope results are expressed in *k*-space. A third SSA for layered seafloors follows from conversion of the usual half-space formulation from *k*-space to coordinate space. This form turns out to be useful for reverberation simulations. The three different approaches are compared with respect to scattering strength and the coherent reflection coefficient, but an assessment of their relative merits will require comparison with exact calculations.

## I. INTRODUCTION

The small-slope approximation (SSA) for rough-interface scattering^{1} has been applied to unlayered media, whether fluid, elastic,^{2,3} or poroelastic.^{4} The intent of this article is to study the SSA for layered fluid media. There is a disconcerting lack of uniqueness in the problem, and three different SSAs will be examined. Gragg and Wurmser^{5} and Jackson^{6} have developed methods applicable to the case in which only the water–sediment interface is rough, with a fluid layer of finite thickness overlying a semi-infinite fluid layer. The SSA has been developed for the general layered case in the electromagnetics literature,^{7–9} and this method is applied in this article to the acoustic problem. A third approach is developed in this article as a possible improvement over other methods, and is a result of a translation of the usual small-slope ansatz to coordinate space.

Any approach to the SSA for layered media must satisfy several conditions: (1) as interface roughness becomes small, the small-slope result should agree with small-roughness perturbation theory; (2) for vanishing roughness, it should yield the usual flat-interface reflection coefficient; (3) it must satisfy reciprocity; and (4) it must revert to the usual half-space approach when there is no layering. All of the approaches considered in this article satisfy all of these conditions. The SSA is usually expressed in *k*-space, but a new small-slope ansatz is developed in coordinate-space. This ansatz is more difficult to implement than that of Refs. 7–9, but appears to be a more physically realistic treatment of phase shifts appearing in the small-slope integrand. Use of the coordinate-space formulation has another advantage, as it is consistent with reverberation models that employ this approach.^{10}

The primary goals of this article are to define and implement three different approaches to the application of the SSA to layered seafloors, and then to show by examples that they yield different results. The question as to which approach is best is not addressed and will require use of accurate numerical methods to provide ground-truth solutions. The Monte Carlo method of Thorsos^{11} could be applied to a discretized version of the integral equations appropriate for the multilayer problem. A general formulation of multi-domain integral equations is given in Ref. 12, and would be appropriate in this case.

The problem is discussed in the following sequence. Section II defines the geoacoustic model to be used and presents basic expressions used in treating the half-space small-slope problem in *k*-space. Results for the small-roughness perturbation theory in layered media, needed for all three small-slope approaches, are given in Sec. III. Next, the three different SSAs for layered media are presented in chronological order in Secs. IV–VI. Section VII treats the coherent reflection coefficient and the bistatic scattering strength. Numerical illustrations are given in Sec. VIII, and conclusions and recommendations are presented in Sec. IX.

## II. BACKGROUND

Figure 1 shows the type of layered seafloor to be discussed in this article. Although two-dimensional roughness is shown, this article develops models for three-dimensional roughness. The acoustic properties of each layer are assumed to be those of a fluid with sediment-water sound-speed ratios $\nu pn$, density ratios $b\rho n$, and loss parameters *δ _{n}*. The complex sound speed ratios incorporating loss are $apn=\nu pn/(1+i\delta n)$. The thickness of the

*n*th layer is denoted

*H*, and the mean

_{n}*z*-coordinate of the

*n*th interface is

*z*with positive

_{n}*z*pointing up. The mean

*z*-coordinate of the sediment–water interface (interface number 1) is taken to be zero.

The discussion of this section will center on the *T*-matrix, which can be viewed as a generalized reflection coefficient capable of characterizing acoustic scattering by rough and/or heterogeneous media. Our goal in this section is to build up an understanding of the fundamentals of the SSA for half spaces in a way that motivates the development of the layered case. Consider a plane wave of unit amplitude incident upon a medium whose average upper boundary is the plane *z* = 0. Assuming $exp\u2009(\u2212i\omega t)$ time dependence, the plane wave can be characterized by the wave vector

where $ex,\u2009ey$, and $ez$ are unit vectors, $k0=\omega /c0$ is the wavenumber in water, and

with $K=Kx2+Ky2$. The subscript *n* is used to designate different seafloor layers, and *n* = 0 corresponds to the overlying water. The *k _{n}* are complex for

*n*> 0, with $kn=k0/apn$.

The scattered pressure field due to an incident plane wave with wave vector having horizontal component $Ki$ is

Note $R=exx+eyy$. Expression (3) can be regarded as the definition of the *T*-matrix. This expression and others in this article place the symbols denoting integration together rather than using them to bracket the integrand. This practice is common in the physics literature and aids in keeping track of complicated integrands. In the interest of brevity, a single integral symbol is used for the multiple integral, and the limits (minus infinity to plus infinity) are not shown.

One of several properties of the *T*-matrix that must be preserved in any small-slope ansatz is reciprocity, which takes the form of Eq. (J.5) in Ref. 13

Most published work on the SSA does not consider layered media. For example, Dirichlet boundaries have been considered,^{14,15} as well as boundaries of homogeneous half spaces (e.g., water above fluid, elastic,^{2,3} or poroelastic^{4} media). The small-slope formalism is a systematic expansion,^{14} but most work, including the present one, employs the lowest-order term in this expansion (lowest order in slope). Some rather subtle questions regarding the small-slope expansion have been addressed by Thorsos and Broschat.^{14} These authors show that the expansion is not actually in powers of true interface slope, as sometimes claimed, but rather in powers of what they term “generalized slope.” Practically speaking, it is less important to acquire deeper understanding of the expansion than it is to assess its accuracy through numerical calculations, and this is the primary effort in Ref. 15.

The *T*-matrix for uniform half-space media is of the general form

In this equation, SSA stands for “small-slope approximation,” and

is the change in the vertical component of the wave vector upon scattering.

The rough interface is defined by the equation $z=f(R)$. The function $A(Ks,Ki)$ is independent of the roughness and only depends on the boundary conditions. It is found from the first-order perturbation theory, for which the *T*-matrix up to first order in the interface relief $f(R)$ is

Here, SPM stands for “small-perturbation method,” and $F(K)$ is the Fourier transform of the relief function

As the first-order perturbation theory obeys reciprocity, $A(Ks,Ki)=A(\u2212Ki,\u2212Ks)$. This can be used to show that the ansatz (5) satisfies the reciprocity condition (4). In Eq. (7), the first term has the flat-interface reflection coefficient as a factor and is the zeroth-order [in $f(R$)] solution. The second term is the first-order correction to the zeroth-order solution. A key requirement is that the small-slope ansatz (5) yields Eq. (7) in the appropriate limits. Expanding Eq. (5), the first-order term agrees with that in Eq. (7). The zeroth-order terms match provided

This relation is absolute; that is, it must be satisfied by the first-order perturbation theory, and no changes or adjustments in the formalism are possible to insure this outcome. Equation (9) requires evaluation of $A(Ks,Ki)$ in the specular direction ($Ks=Ki$). For this, Voronovich (Ref. 16, Chap. 5) provides the following perturbative method relating the first-order function to the reflection coefficient. The reflected pressure with flat interface and a unit-amplitude incident plane wave is

If the interface is perturbed by being moved upward a small distance *f*, the first-order change in pressure is

where $V(K,f)$ is the reflection coefficient (referred to *z* = 0) when the interface is moved upward by *f*. Comparing with Eq. (3), the first-order *T*-matrix must be

As the perturbation in interface relief is independent of position, Eq. (8) gives $F(K)=f\delta (K)$, from which Eq. (7) gives

For the half space, the reflection coefficient when the interface is at some general *z* = *f* is

where the reflection coefficient is still referred to the plane *z* = 0, and the exponential factor gives the phase change resulting from movement of the actual interface. Inserting Eq. (14) into Eq. (13), one finds that the consistency condition [Eq. (9)] is satisfied.

Three important points can now be made. First, the discussion of consistency made no reference to the properties of the half-space medium. It could be fluid, elastic, or poroelastic. Second, the ansatz lends itself to formal averaging if the interface relief is assumed to be a stationary random process. In this case, the expected value of the exponentiated random variable (multiplied by $\u22121=i$) involves the exponentiated second moment. Formal averages are discussed in Sec. VII. Finally, it should be noted that the ansatz cannot be used for layered media, as Eq. (14) does not apply in this case. Thus, Eq. (5) must be replaced by a more general ansatz if layered media are to be treated. There is no unique choice for a more general ansatz, and three different candidates are discussed below.

## III. SMALL-ROUGHNESS PERTURBATION FOR LAYERED MEDIA

Any small-slope ansatz must agree in the appropriate limits with the small-roughness perturbation theory, and this gives rise to a consistency condition such as Eq. (9). This relation applies only to the case of homogeneous half spaces. Berrouk *et al.*^{7} give a consistency condition for layered media. This condition is reviewed, and a missing proof that it is satisfied for any number of layers is supplied. As there is such a strong connection between the SSA and perturbation theory, a summary of the latter as applied to layered media is also given here. Published work on the small-roughness perturbation theory for layered media, for example,^{7,17} is not presented in forms convenient for the immediate problem. Consequently, a very brief derivation is given here.

The starting point is a general expression for the first-order *T*-matrix (from perturbation theory) for scattering by the *n*th interface of a layered fluid medium

Here, $Fn(K)$ is the Fourier transform of the *n*th relief function $fn(R)$

Expression (15) defines $An(Ks,Ki)$. For scattering by roughness of the sediment–water interface $An(Ks,Ki)$ must be identical to that in Ref. 18 and Chap. 13 of Ref. 13

The factors *a*_{1} and *b*_{1} will be defined shortly. This expression accounts for layering below the sediment–water interface through its effect on the sediment-water reflection coefficient, *V*(*K*). All that is necessary to generalize to buried rough interfaces (*n* > 1) is to replace this reflection coefficient by that at the interface of interest and to introduce transmission factors that account for plane-wave propagation through layers above that interface

The amplitudes of down-going waves just above the *n*th interface are denoted $An\u22121(K)$ with the unit-amplitude plane wave defined by **K** incident on the water–sediment interface. The function $A\u0303n(Ks,Ki)$ is similar to $A1(Ks,Ki)$ for the case in which a unit-amplitude plane wave is incident on interface *n* with the layer above (*n* – 1) taken to be semi-infinite. While in $A1(Ks,Ki)$ the upper medium is assumed to be water, here, it is assigned the properties of layer *n* – 1.

and

In these equations, $Vn(K)=Bn\u22121(K)/An\u22121(K)$ is the plane-wave reflection coefficient at the *n*th interface, taking into account all layering below this interface. The down-going wave amplitude $An(K)$ was introduced previously, and $Bn\u22121(K)$ is the amplitude of the up-going wave immediately above the *n*th interface. Note $b\rho 0=ap0=1$. Expressions (18)–(21) complete the construction of $An(Ks,Ki)$. Reciprocity is obeyed and takes the form $An(Ki,Ki)=An(\u2212Ks,\u2212Ks)$. These functions are employed in all the small-slope methods to be discussed, and appear in the consistency condition given by Berrouk *et al.*^{7} for layered media

Here, *V*(*K*) is the reflection coefficient of the entire seafloor when all interfaces are flat, referred to the water–sediment interface. For a medium with several layers, the terms on the left- and right-hand sides of Eq. (22) are very complicated, making direct verification difficult. The authors of Ref. 7 state that they have verified Eq. (22) for $N=$ 1 and 2 interfaces, but indicate that numerical calculations are required to verify it for larger *N*. In the following, a general proof of the validity of Eq. (22) for all *N* is given, avoiding use of analytical expressions for the individual terms. Considerable simplification is possible because Eq. (22) requires evaluation of $An(Ks,Ki)$ only in the specular direction ($Ks=Ki$). For this, Voronovich's perturbative method [Eq. (13)] can be used to obtain

The consistency condition [Eq. (22)] can now be expressed as

The needed treatment of reflection by layered seafloors departs from typical published approaches. One such approach is to iterate the expression for the reflection coefficient for a layer of finite thickness placed above a layered medium of a known reflection coefficient, starting from the lowest interface and working upward, for example, Chap. 1 of Ref. 19. Another is to iterate the impedance in similar fashion, for example, Chap. 3 of Ref. 20, and use this to obtain the reflection coefficient. The propagator matrix method of Gilbert^{21} proceeds in the opposite direction from the uppermost interface to the lowest. A matrix approach is used here as well, but it will proceed from the lowest interface upward to the water–sediment interface. The reflection coefficient can be expressed in terms of a matrix product as follows:

The column vectors represent the amplitudes of up- and down-going waves, with the column vector on the right expressing the fact that there is only a down-going wave in the semi-infinite *N*th layer with amplitude *τ* corresponding to the transmission coefficient $\tau (K)/A(K)$. The amplitude of the plane wave incident from the water onto the sediment is denoted *A*(*K*), and is not equal to unity, as $\tau (K)$ is held fixed as the various interfaces are moved in forming derivatives. The column vector on the left expresses the fact that a down-going wave in the water of amplitude *A*(*K*) produces an up-going wave of amplitude $A(K)V(K)$. The square matrices *P _{n}* are

for *n* > 0, where *H _{n}* is the thickness of the

*n*th layer. The matrices

do not depend upon layer thicknesses. They are conversion matrices giving the amplitudes immediately above the *n*th interface (in layer *n* – 1) in terms of the amplitudes immediately below (in layer *n*), where the impedance ratio *Z _{n}* associated with the

*n*th interface is

The matrices *P _{n}* include conversion into (the $Cn+1$ matrix) and propagation through (the square matrix involving the exponentials) the

*n*th layer. The matrix associated with the sediment–water interface is a special case, and is

Next, consider the derivatives in the sum of Eq. (24) starting with the last.

This expression follows from the fact that variations in the depth of the *N*th interface cause corresponding (but opposite) changes in the thickness of the $(N\u22121)$th layer. The next-to-last term in Eq. (24) is

The two derivative terms follow from the decrease in thickness in layer *N* – 2 and increase in thickness of layer *N* – 1 as interface *N* – 1 is moved upward. Note that when Eqs. (30) and (31) are added, as required in Eq. (24), the derivative terms with respect to the thickness of the $(N\u22121)$th layer cancel. Proceeding to other terms in the sum, this cancellation continues up through the *n* = 2 term giving

The uppermost interface must be treated separately with derivative

The first term on the right-hand side of Eq. (33) cancels Eq. (32). Evaluating the second term, the sum over all derivatives is

As a final step, one can use Eq. (34) and

to obtain the consistency condition in the form of Eq. (24).

## IV. LAYERED MEDIA: ACOUSTICS LITERATURE

Gragg and Wurmser^{5} and Jackson^{6} have given small-slope approaches to layered seafloors in which a fluid layer of finite thickness sits above a semi-infinite fluid “basement.” The water–sediment interface is rough, and the layer–basement interface is flat. As the approach of Ref. 5 is a special case of that in Ref. 7, this section will only consider the method of Ref. 6. The *T*-matrix is

where “SSL1” indicates that this is the first of three SSAs for layered media to be considered, and $V(K,f)$ is the reflection coefficient for the layered medium referred to *z* = *f*, with the water-sediment boundary moved upward by *f*, while all other layer interfaces are unchanged in the *z*-coordinate,

Here, $V1(K)$ is the reflection coefficient that would be observed if the upper layer were extended downward to be semi-infinite,

with *Z*_{1} given by Eq. (28). In Eq. (37), $V2(K)$ is the reflection coefficient (referred to *z* = 0) that would be seen if the upper layer extended to all positive *z*, but the layer-basement interface is not moved,

with *Z*_{2} given by Eq. (28). The function $\Psi (Ks,Ki)$ is obtained by forcing agreement with first-order perturbation theory and is

In this equation, the second argument of $V(K,f)$ is suppressed. It is set to zero, including the terms involving the derivative of the reflection coefficient. The derivative of the reflection coefficient is

with all arguments on the right-hand side suppressed but equal to *K*. Note that this derivative is not the same as that in Eq. (11) because of the difference in reference planes. It can be shown that $\Psi (Ki,Ki)=1$. Setting $f(R)=0$, it follows that the ansatz gives the correct zeroth-order result, the first term in Eq. (7). By expanding Eq. (36), it can be shown to be consistent with the first-order perturbation theory, thanks to the construction embodied in Eq. (40). It also satisfies reciprocity and, of course, obeys the consistency condition [Eq. (22)], as it must, agreeing in zeroth- and first-order perturbations.

One unattractive feature of the present ansatz is that the *T*-matrix is a complicated function of interface relief, so the simple averaging method used in the half-space problem does not apply. Another drawback is that it is not obvious how this approach might be extended to multiple rough interfaces, but the two following approaches do not have this problem.

## V. LAYERED MEDIA: ELECTROMAGNETICS LITERATURE

A small-slope *T*-matrix for the layered-fluid case is readily adapted from the electromagnetic case treated in Refs. 7–9, and is

As noted earlier, the function $fn(R)$ is the relief of the *n*th interface, and the $An(Ks,Ki)$ are given in Sec. III. This ansatz is similar to that of Gragg and Wurmser,^{5} but their version, intended for scattering in the case of a finite-thickness layer over a semi-infinite basement, does not include the *n* = 2 term in the sum. This term does not appear in the scattering cross section, so its absence does not affect the results of Ref. 5. The *n* = 2 term is needed, however, in obtaining the coherent field even if buried interfaces are flat.

Setting all the relief functions to zero in Eq. (42) gives

This is to be compared with the *T*-matrix for a medium with flat interfaces and reflection coefficient *V*(*K*)

One rather non-intuitive consequence of the consistency condition is that it involves *all* interfaces between layers, even those that are flat, not rough. It can be seen that the layered small-slope approach of Ref. 5 would be equivalent to Ref. 7 if a term for the buried, flat interface had been added to the *T*-matrix. As presented in Ref. 5, the scattered field agrees with that given by Ref. 7, but the buried interface term is needed in order to obtain the coherent field.

This SSA will be denoted SSL2 and has some attractive properties. The formalism is not much more complicated than the half-space case, and formal averaging proceeds in the same fashion as will be seen in Sec. VII. One disturbing note is that the phase term inside the Kirchhoff integral containing the interface relief involves the wavenumber in water but not the wavenumbers for the various layers. This may be an oversimplification, since small changes in each layer interface result in complex non-linear changes in the phase and amplitude of the scattered field, rather than the linear change in phase implied by Eq. (42). Although SSL2 satisfies the consistency condition [Eq. (22)], it does not follow that the consistency condition gives the correct change in the scattered field due to uniform perturbations of an interface but only for small perturbations. This can be confirmed by comparing Eqs. (42) and (11), which describes the first-order change in the scattered field due to a perturbation in a layer height. These two are clearly related, which seems to imply that SSL2 is restricted to small perturbations of the layer interface. Any larger uniform perturbations would change the background field, requiring orders beyond first in height. This contradicts the requirement that any small-slope ansatz must be valid for all orders in height.^{1,14} In the SSL1 approach,^{6} buried rough interfaces cannot be treated, but phase depends on the sound speeds of both the water and the uppermost layer. In Sec. VI, a small-slope ansatz is developed that can treat buried rough interfaces with the phase associated with a given boundary dependent on the sound speeds of the media on either side of the boundary.

## VI. NEW SSA FOR LAYERED MEDIA

A new small-slope ansatz (denoted SSL3) for layered fluid media is developed in this section. The ansatz is most simply expressed in coordinate space in contrast to SSL1 and SSL2, which are expressed in *k*-space. SSL3 overcomes the drawback of SSL2's phase term discussed in Sec. V.

The developments to follow will be summarized briefly. First, the scattered field in coordinate space using SPM is reviewed. Next, the half-space small-slope ansatz, Eq. (5), is used to express the scattered field in this approximation in coordinate space, and several identities involving Green's functions are used to bring it into a similar form to the coordinate space representation of SPM. This form has an obvious generalization to buried rough interfaces. By comparing the small slope and SPM expressions for the scattered field, the angular factor $An(Ks,Ki)$ can be identified, and the identities required by any small-slope ansatz can be verified. Since the coherent reflection coefficient and scattering cross section are wavenumber-domain quantities, the small-slope scattered field is converted back into **K** space (i.e., the *T*-matrix), again through the use of Green's function identities.

Because of the close connection between the SSA and the SPM, the discussion begins with a coordinate-space SPM expression for the half-space case. This expression serves as a model for the sort of small-slope expression that is desired. The scattered field at **r** due to a unit amplitude point source at $r0$ to first order is^{10}

The Green's functions and their derivatives are evaluated an infinitesimal distance above the interface. The following unperturbed Green's function for the field above the interface appears in Eq. (45)

The first term inside the square brackets corresponds to the direct path from the source to the receiver, and the second to the reflected path from the source to the seafloor and then to the receiver. The vertical coordinate *z* increases in the upward direction and is zero at the sediment–water interface. The source height *z*_{0} is positive, and the Green's function is normalized such that its value when the field point **r** is very close to the source point $r0$ (i.e., the separation is much less than a wavelength) is $1/|r\u2212r0|$.

A means of converting *k*-space methods to coordinate space is provided by expression (J.4) from Ref. 13. The scattered field at position **r** due to a unit point source in the water at position $r0$ is given in terms of the *T*-matrix as

Inserting the half-space small-slope *T*-matrix [Eq. (5)], the scattered field in the SSA is

For the half-space case, $A(Ks,Ki)$ can be specialized to

where

In order to bring the small-slope result [Eq. (48)] into a form resembling the SPM result [Eq. (45)], the Green's functions in the latter must be identified as arising from the integrals over $Ki$ and $Ks$ in Eq. (48). This requires that the integrand in Eq. (48) can be factorized into separate functions of $Ki$ and $Ks$. A barrier to this factorization is seen in the factor $1/\Delta kz=1/{k0[\beta 0(Ks)+\beta 0(Ki)]}$. This is not a simple product, and so it prevents the desired factorization. Factorization can be accomplished by using a “trick” in which the offending part of the integrand is expressed as a factorizable integral. Note that

Rearranging a bit, a factor in Eq. (48) can be expressed as

Considering Eq. (6), it can be seen that the integrand in Eq. (52) can be factorized into separate functions of $Ki$ and $Ks$. To proceed further, careful inspection of Eq. (48) supplemented by Eq. (52) shows that it contains the following unperturbed Green's function:

A key point is that this is the half-space Green's function with the interface shifted vertically by *f*. It is evaluated at an infinitesimal distance above the interface, at $z=f+$, and the point $r0$ is in the water ($z0>0$). Examination of Eqs. (49) and (50) shows that the term containing the factor $Ks\u22c5Ki$ corresponds to a product of transverse derivatives of the Green's functions. Similarly, the term containing the factor $\beta 1(Ki)\beta 1(Ks)$ corresponds to a product of vertical derivatives of the Green's functions *below* the interface. These can be converted to derivatives above the interface by using the fact that the vertical derivative divided by the density is continuous across the interface. The other terms correspond to products not involving derivatives. The last term on the right-hand side of Eq. (52) cannot be factorized but can be recognized as giving the zeroth-order field.

Thus, Eq. (48) finally can be expressed as

This expression is of the form used in reverberation simulations.^{10} One bonus of the coordinate-space approach is that the two small-slope constraints can be verified with little effort. First, when $f(R\u2032)=0$, Eq. (54) trivially gives the desired zeroth-order result, as the first term vanishes in this case. Second, when $f(R\u2032)$ is small, the Green's functions can be expanded in a Taylor series around *f* = 0, and the integral over *f* in the first term of Eq. (54) simply gives a factor $f(R\u2032)$ and the desired first-order result.

The integral over *f* may seem unnatural in the context of small-slope expressions for the scattered field, but it should be recognized that it arose from a mathematical identity rather than an *ad hoc* ansatz. This integral appears in Berman's coordinate-space small-slope treatment of the pressure-release interface.^{22} The present result [Eq. (54)] is equivalent to Eq. (21) in Ref. 22 after specialization to the Dirichlet boundary, conversion to the time domain, and allowing for differences in normalization. Berman shows that the integral over *f* can be performed analytically in this case, yielding a very efficient means for numerical computation. It seems doubtful that this can be accomplished in the fluid-fluid case.

The advantages of the coordinate-space approach are (1) the phase involving the rough interface explicitly takes into account the field at the rough surface due to the layering (through the Green's functions), (2) uniform perturbations to layer (const *f*) change the field in the sediment exactly due to the use of the Green's functions, and (3) the generalization to layered media is simple. The scattered field due to roughness of the *n*th interface is taken to be

The index *n* takes on the value of unity for the water–seafloor interface and larger integral values for buried interfaces. The unperturbed coordinate of the *n*th interface is *z _{n}*, and it is negative. The Green's function $G(0)(r\u2032,r,f)$ is that for the layered medium with the

*n*th interface displaced vertically by

*f*. The sound-speed and density ratios $apn$ and $b\rho n$ are ratios with respect to water, not to the layer above. As noted earlier $ap0=1$ and $b\rho 0=1$. The total scattered pressure is the sum of Eq. (55) over all

*n*plus the zeroth-order reflected pressure. Reciprocity is manifest, the zeroth-order constraint is built in, and comparison of Eq. (55) as $fn(R)$ becomes small with expression (13) in Ref. 17 shows that the first-order constraint is satisfied.

Expression (55) yields the following computational “recipe” for reverberation simulations:

Compute Green's functions and their derivatives at all ranges $R\u2032$ just above interface

*n*with this interface displaced by a range of heights ($\u2212fmax<f<fmax$), where $fmax$ is the largest anticipated value of relief.For each $R\u2032$ evaluate the usual perturbation-theory integrand [except for the factor $fn(R\u2032)$] over all

*f*from zero to $f=fn(R\u2032)$.Integrate over

*f*from 0 to $fn(R\u2032)$.Integrate over $R\u2032$.

Sum over all interfaces.

Add flat-interface pressure.

Section VII develops formal averages for the coherent reflection coefficient and the scattering strength. These are plane-wave concepts and require that a *k*-space formulation be available. For SSL3, the initial steps leading from *k*-space to coordinate space can be followed in reverse. First, the *k*-space representation of the Green's function [Eq. (46)] is substituted in Eq. (55). Next, comparison of the altered version of Eq. (55) with Eq. (47) allows identification of the *T*-matrix for scattering from the *n*th layer

These contributions must be evaluated and summed for all rough interfaces, and the complete *T*-matrix consists of this sum plus the flat-interface *T*-matrix, $V(Ki)\delta (Ks\u2212Ki)$. While this last term may seem peculiar, it must be realized that parts of the sum over *n* cancel part of the flat-interface result, leading to an alteration in the coherent reflection coefficient compared to the flat-interface case. In the spatial domain version of the half-space ansatz, the zeroth-order pressure must be added as well. The notation of Eq. (56) obscures the steps required to implement this approximation. The function $An(Ks,K,f)$ is the function $An(Ks,K)$ defined in Eq. (18) evaluated for a vertical shift of the *n*th interface by an amount *f*. Such a shift will change both the reflection coefficient $Vn(K)$ and the wave amplitude $An\u22121(K)$. These matters will be made more concrete in Sec. VII.

## VII. FORMAL AVERAGING

This section employs first and second moments of the *T*-matrix to determine the coherent reflection coefficients and the scattering strengths for the three different SSAs. The coherent reflection coefficient $Vcoher(K)$ can be found from the expression

Following Ref. 16, Chap. 2, Ref. 23, and using the notation of Ref. 13, the scattering cross section per unit area per unit solid angle, $\sigma (Ks,Ki)$, can be found using

For brevity $\sigma (Ks,Ki)$ will be referred to as the “scattering cross section,” even though this implies dimension length squared when it is actually dimensionless. In most applications one deals with $10\u2009\u2009log10[\sigma (Ks,Ki)]$, the bistatic scattering strength.

It is necessary to assume that the relief functions $fn(R)$ are stationary Gaussian processes. The second moments of the relief functions, the “covariances,” provide a complete statistical characterization owing to the assumption of Gaussian behavior.

where $hn2$ is the mean-square relief of the *n*th interface, $\rho n(R\u2212R\u2032)$ is the spatial correlation function [with $\rho n(0)=1$], and $\delta nn\u2032$ is the Kronecker delta, equal to unity if $n=n\u2032$ and zero otherwise. While the assumption of stationarity implies a seafloor of infinite extent, this is an approximation that does not interfere with practical application, just as time series analysis often is based on stationarity in the time domain. It is assumed that the roughness is uncorrelated between the different interfaces, but this assumption could be avoided by employing cross-covariances as in Ref. 24. Specific choices for mean-square relief and spatial correlation are made in Sec. VIII (Numerical examples). A key expression used in formal averaging in SSAs is

where *x* is a Gaussian random variable having zero mean, and where “$\u27e8\u2009\u27e9$” denotes averaging over a theoretical infinite ensemble. As the discussion progresses, it will emerge that the covariances enter through formal averages of the following general form:

All three SSAs employ the following integral:

This integral is commonly used in the Kirchhoff approximation. The argument *η* is different for each approximation and will be defined as each is discussed. If the relief statistics are isotropic, that is, if the spatial correlation $\rho n(R)$ depends only on the magnitude *R* of the vector **R**, the integral (62) can be written as

where $J0(x)$ is the zeroth-order Bessel function of the first kind, and $\Delta K=|\Delta K|$ with $\Delta K=Ks\u2212Ki$. The integrals $In(\eta )$ are usually sharply peaked in the specular direction where $\Delta K=0$.

### A. Formal averaging for SSL1

The results of this subsection are limited to the relatively simple case of a rough fluid layer over a flat, semi-infinite fluid basement. The *T*-matrix for SSL1 given in Eq. (36) is difficult to average formally because of the square roots of the reflection coefficients, which themselves contain exponential factors in the denominator. Exponential factors can be brought to the numerator if the denominator of the reflection coefficient is expanded in a geometric series

This expression is also used in formal averaging for SSL3. Taking the first moment of Eq. (36), the assumption of stationarity for $f(R)$ causes the integrand to be independent of **R** with the exception of the factor $exp[i(Ks\u2212Ki)\u22c5R]$, which yields a delta function. Comparison with Eq. (57) gives

In application, the sum is truncated at some finite *m* = *M*. No rule is offered for determining *M*, rather *M* is found by increasing the number of terms until the result shows no significant change.

It is more difficult to obtain the second moment of Eq. (36) as needed for the scattering cross section. As a first step, the square root of the reflection coefficient appearing in Eq. (36) is expanded in a power series in $exp\u2009(2ik1\beta 1f)$,

where the first few expansion coefficients are

The *T*-matrix is then expressed as a sum over *m*_{1} and *m*_{2} with corresponding powers of $exp\u2009(2ik1\beta 1f)$. The second moment of the *T*-matrix can be found using Eq. (60), from which the scattering cross section follows:

Here,

and the argument of the integral (63) is

In applying SSL1 to numerical examples to be presented later, a pitfall was discovered in performing the multiple sums. It was found that convergence of the sums was greatly improved if only terms for which $m1+m2\u2264M$ and $m3+m4\u2264M$ were retained. When this is done, the factors multiplying $exp\u2009(2mik1\beta 1f)$ are complete for $0<m<M$; that is, they will not change as higher values of *m*_{1}, *m*_{2}, *m*_{3}, and *m*_{4} are included. If, on the other hand, a term such as $m1=1,\u2009m2=M$ is included, this will bring in a term in $exp\u2009[2(M+1)ik1\beta 1f]$. This term is not complete because, for example, the excluded term $n1=0,\u2009n2=M+1$ would also contribute to the sum. There is evidently a large degree of cancellation in the expansion coefficients such that incomplete factors are inaccurate and detrimental to convergence. This conclusion has been reinforced by considering the backscatter case for which *K _{i}* =

*K*so that the square root in Eq. (36) disappears and the simpler geometric series (64) can be used. The backscattering cross section obtained in this manner agrees with the result obtained by restricting the sums.

_{s}### B. Formal averaging for SSL2

The second SSA under consideration, SSL2, can be applied to the general case in which an arbitrary number of layers is allowed, and all interfaces may be rough. Averaging is straightforward, being no more complicated than that used for the half-space small-slope. The coherent reflection coefficient is found by taking the first moment of Eq. (42) and comparing the result to Eq. (57). The Gaussian random process of interest for substitution in Eq. (60) is

Then

As $\u27e8x2\u27e9$ does not depend upon the integration coordinate **R** in Eq. (42), the integration gives a delta function. Comparing the result to Eq. (57) gives the following general expression for the coherent reflection coefficient:

In the numerical examples to be presented later, there are two layers (*N* = 2), and the following evaluations of Eq. (23) are needed.

Here, $\theta =cos\u22121(K/k0)$ is the grazing angle in water. These functions are also used in evaluating the coherent reflection coefficient for SSL3.

Turning to the scattering cross section per unit area per unit solid angle, it is necessary to evaluate

where

Omitting steps that are encountered in many derivations of the SSA, the bistatic scattering cross section for the general case of an arbitrary number of layers is

### C. Formal averaging for SSL3

Formal averaging for the new small-slope ansatz, SSL3, is more difficult than for SSL2. As for SSL1 and SSL2, numerical examples to be given later consider only the simplest example: a rough fluid layer of finite thickness over a flat, semi-infinite fluid basement. For the moment, however, consider the general case in which all interfaces may be rough. Then, Eq. (56) with the zeroth-order reflection included gives the *T*-matrix

where the $An(Ks,Ki,f)$ are found from Eqs. (17) and (18) by displacing the *n*th interface by the amount *f _{n}*.

Taking the first moment of Eq. (80), the coherent reflection coefficient for the general case is

The remainder of this discussion of formal averaging for SSL3 is restricted to the special case in which a fluid layer with a rough upper interface lies above a flat fluid basement. The following function is needed for derivation of both the coherent reflection coefficient and the scattering cross section:

In this equation, *a* is *a*_{1} in Eq. (20), *b* is *b*_{1} in Eq. (21), $V(K,f)$ is given by Eq. (64), and $\Delta kz$ is given by Eq. (6). The factor $exp\u2009(\u2212i\Delta kzf)$ arises because the *T*-matrix is referred to the plane *z* = 0. Using this reference does not limit the validity of SSL3 in the same way that SSL2 is limited by only using the water sound speed in the phase terms. The final result for the coherent reflection coefficient is

where

The first of the $V\xaf(k)$ is

and, using the expansion (64),

As for SSL1, the sums will be truncated in any application.

The double sums in Eq. (87) appear because the *T*-matrix (80) contains a product of the reflection coefficient with itself. In deriving the scattering cross section, the *T*-matrix is squared, so quadruple sums will appear. As with the coherent reflection coefficient, it is assumed that only the sediment–water interface is rough. As the derivation involves extensive and tedious algebra and a large number of terms, it is merely sketched. Although the cross section involves the second central moment of the *T*-matrix (58), it is sufficient to evaluate the ordinary second moment, the first term on the left-hand side of Eq. (58). This is because the subtraction of the second term on the left-hand side removes all terms containing a product of two delta functions that would otherwise appear. These terms represent the coherently reflected field. In taking the second moment, however, a single delta function $\delta (Ks\u2212Ks\u2032)$ appears, as demanded by Eq. (58). Out of the plethora of terms that arise in the second moment, all except those of the form

can be ignored. Integrals of this form are the first term in the Kirchhoff integral (62). The missing second term will always be present and can be added with confidence, as it removes a delta function that is part of the coherent field. The following factor appears in each term in the sums:

The arguments *γ _{a}* and

*γ*are defined later.

_{b}The scattering cross section can be written in the form

In order to give the flavor of the cross section in SSL3, the function *q*_{1} is given here. The functions $q2,\u2009\u2009q3,$ and *q*_{4} are more complicated and are given in the Appendix. For *χ* = 1

and for $\chi =\u22121$

Here,

The variable *χ* has been introduced merely to keep track of conjugation of the reflection coefficient and is always +1 when there is no conjugation and −1 otherwise. This notation is not needed in the quadruple-sum term *q*_{4}, which is used with only a single set of arguments. In numerical examples to follow, the multiple sums are truncated after only a few terms and usually do not impose an excessive computational burden.

## VIII. NUMERICAL EXAMPLES

The simplest case that exercises the features of the three SSAs for layered seafloors is illustrated in Fig. 2. A sediment layer of finite thickness $H1=1$ m overlies a semi-infinite basement, and only the water–sediment interface is rough. The acoustic frequency is taken to be 2 kHz, and the water sound speed is 1500 m/s. Both the finite and semi-infinite sediments are fluids characterized by density ratio $b\rho n$, sound-speed ratio *ν _{n}*, and the “loss parameter”

*δ*, related to attenuation in dB/m as follows:

_{n}Three illustrative cases will be presented with parameters given in Table I. The “low-contrast case” has relatively little difference in the parameters for the finite and semi-infinite layers. Thus, the reflection coefficient at interface two is rather small. As a result, the three approximations (1)–(3) give similar results. The reflection coefficient in the high-contrast case is considerably larger, and more difference is seen in the approximations. The greatest difference occurs for the mud-layer case, where basement reflection is very important.

The roughness of the water–sediment interface is described by the “von Kármán” spectrum

This is an isotropic spectrum, dependent only on the magnitude of the two-dimensional wave vector **K**. The covariance of the relief function is the inverse Fourier transform of the spectrum.

The mean-square roughness is

and the spatial correlation is

where $\Gamma (x)$ is the gamma function, and $K\nu (x)$ is the modified Bessel function of the second kind of order $\nu =\gamma 2/2\u22121$.

In all computations to follow, $w2=0.001\u2009m4\u2212\gamma 2,\gamma 2=3.0$, and $K0=1.0\u2009\u2009m\u22121$. Given the acoustic frequency of 2 kHz, $k0h1=0.664$, which indicates the roughness may be large enough to pose a problem for the SPM. This is the normal expectation for pressure-release surfaces, where lowest-order perturbation theory is expected to be accurate when this product is much less than unity. For fluid–fluid interfaces, however, Thorsos *et al.*^{25} show that accurate results can be obtained even when this product is equal to one.

The three SSAs will be compared with respect to coherent reflection coefficient and backscattering strength. It should be noted that, in the cases considered, the scattering strengths given by SSL2 are the same as would be obtained using the formalism of Gragg and Wurmser.^{5}

### A. Low-contrast case

In this section, the parameters for the low-contrast case in Table I are used. Using Eqs. (66), (74), and (83), the coherent reflection coefficient can be computed for all three approximations. The transverse component of the wave vector appearing in these expressions is $K=k0\u2009cos\u2009\u2009\theta $, where *θ* is the grazing angle.

Figure 3 compares the magnitude of the coherent reflection coefficient for SSL1 with an upper summation limit *M* = 2 with that for the flat-interface coefficient $V(k0\u2009cos\u2009\u2009\theta )$ and for the “Eckart” formula. The latter reflection coefficient arises in the Kirchhoff approximation and also in the SSA for the half-space case. Although there is no theoretical justification for its use in the layered case, it is of interest as a possible simple alternative to more complicated approximations and is given by the expression

The oscillations in the reflection coefficient magnitudes are due to interference caused by reflection from the buried (flat) interface. The coherent reflection coefficient of SSL1 is smaller than the flat-interface reflection coefficient, as expected, with the deficit being due to incoherent scattering quantified by the scattering cross section. The Eckart curve gives a fair approximation to the SSL1 curve. Figure 4 compares the magnitudes of the coherent reflection coefficients in the three approximations using an upper summation limit *M* = 2 for SSL1 and SSL3. Increasing this limit produces no noticeable change in the results. The coherent reflection coefficient magnitudes match quite well for all three approximations. As the coherent reflection coefficient is determined by the amount of incoherently scattered power that is subtracted from the flat-interface reflection, this comparison suggests that the scattering strengths, which measure incoherent scattering, should match as well.

As noted earlier, the scattering strength is the dB version of $\sigma (Ks,Ki)$, the scattering cross section per unit area per unit solid angle. Figure 5 compares the backscattering ($Ki=\u2212Ks=exk0\u2009cos\u2009\u2009\theta $) strength for the water–sediment interface as computed using the SSA (SSL3), the SPM, and the Kirchhoff approximation. As expected, the SSA and perturbation theory agree at smaller grazing angles. Also as expected, the Kirchhoff approximation fails at small grazing angles. It is interesting that the small-slope and Kirchhoff approximations do not agree near the specular direction, which is vertical incidence in this backscattering example. This contrasts to the half-space case in which they agree. This may have implications in the use of the Kirchhoff approximation for layered seafloors (see Ref. 26, Chap. 6 and Refs. 27–30). Agreement between the Kirchhoff and SSAs for half-spaces in the specular direction can be proven analytically. A similar demonstration for the layered case cannot be carried out, to the authors' knowledge, and this is most likely because the Kirchhoff approximation does not properly account for phase and amplitude changes as the scattering interface is displaced.

The SPM scattering strength rises above that for the SSA for grazing angles approaching 90°. Assuming this shows inaccuracy in the SPM, it is at odds with the results of Ref. 30, which show SPM to be accurate in the roughness regime employed here. That previous work, however, was restricted to uniform half-spaces and may not be applicable here. The large scattering strength near normal incidence may be indicative of the inability of the SPM to enforce energy conservation. This problem arises because the scattered intensity in SPM rises linearly with mean-square roughness. As a result, there is no upper bound to scattered intensity even though energy conservation demands one.

The three different SSAs are compared for backscattering by the water–sediment interface in Fig. 6. These computations used Eq. (69) for SSL1, Eq. (79) for SSL2, and Eq. (90) for SSL3. The summation limit *M* = 2 was used for SSL1 and SSL3. An upper summation limit of *M* = 1 for SSL3 produced essentially identical results, while the limit *M* = 0 gave noticeable error. The three approximations give backscattering strength curves that are practically indistinguishable. This accords with the match seen in the coherent reflection coefficients.

### B. High-contrast case

While the low-contrast geoacoustic parameters in Table I produce nearly identical results in all three SSAs, the high-contrast parameters lead to noticeable differences. For example, in Fig. 7, the coherent reflection coefficient for SSL2 differs significantly than that for SSL1 and SSL3 (both computed with upper summation limit *M* = 2). The latter two approximations, however, are essentially indistinguishable. Interestingly, the flat-interface reflection coefficient is *smaller* than the coherent reflection coefficients near normal incidence for all three approximations. The flat-interface coefficient is suppressed by wave interference to such a degree that scattering can yield an increase in coherent reflection. In support of this statement, note that in Fig. 7, the flat-interface coefficient essentially vanishes for a grazing angle of 80°. In light of this, it is not surprising that the roughness scattering contribution to the coherent reflection coefficient would partially fill in the region near 80°. This interpretation is borne out by the fact that, when the layer thickness is changed from 1 m to 0.75 m, the flat-interface reflection is stronger than the coherent reflection.

Figure 8 compares backscattering strengths computed using SSL1, SSL2, and SSL3 for the high-contrast case. The upper summation limit *M* = 3 was necessary for SSL1, but *M* = 2 was sufficient for SSL3. Figure 8 also gives the SPM result for backscattering, and it can be seen that, as the grazing angle increases, the small-slope backscattering strengths fall below the SPM curve. The curves for the three SSAs are similar except for the region between 40° and 60°. This region is magnified in the lower panel of Fig. 8.

### C. Mud-layer case

The magnitudes of the coherent reflection coefficients for the three SSAs are shown for the mud-layer case of Table I in Fig. 9. The curves for SSL1 and SSL3 are similar and do not depart greatly from the flat-interface curve. The curve for SSL2 is quite different and shows a peculiar upturn at small grazing angles. The mud-layer parameters yield large differences in backscattering between the three SSAs as seen in Fig. 10. For small-to-medium grazing angles, the SSL2 curve follows the SPM curve, including the deep nulls at grazing angles below 45°. Unlike the other two SSAs, SSL1 shows suspicious behavior, departing from SPM even at rather small angles, and exhibiting narrow peaks not seen in any of the other approximations. Given that the coherent reflection coefficients for SSL1 and SSL3 are similar, it is surprising that their backscattering strengths are quite different. The nulls for SSL3 are less pronounced than those for SSL1 and SSL2, and this difference, in addition to those noted above, makes the mud-layer case a good test for the superiority of one of the three approximations if they are compared to accurate numerical calculations.

Case . | Layer number . | Thickness (m) . | Sound speed ratio . | Density ratio . | Loss parameter . |
---|---|---|---|---|---|

Low contrast | 1 | 1 | 1.1 | 1.8 | 0.02 |

2 | $\u221e$ | 1.2 | 2.0 | 0.01 | |

High contrast | 1 | 1 | 1.05 | 1.8 | 0.02 |

2 | $\u221e$ | 1.8 | 2.5 | 0.01 | |

Mud layer | 1 | 1 | 0.99 | 1.4 | 0.005 |

2 | $\u221e$ | 1.8 | 2.5 | 0.01 |

Case . | Layer number . | Thickness (m) . | Sound speed ratio . | Density ratio . | Loss parameter . |
---|---|---|---|---|---|

Low contrast | 1 | 1 | 1.1 | 1.8 | 0.02 |

2 | $\u221e$ | 1.2 | 2.0 | 0.01 | |

High contrast | 1 | 1 | 1.05 | 1.8 | 0.02 |

2 | $\u221e$ | 1.8 | 2.5 | 0.01 | |

Mud layer | 1 | 1 | 0.99 | 1.4 | 0.005 |

2 | $\u221e$ | 1.8 | 2.5 | 0.01 |

## IX. DISCUSSION

The SSA applied to layered seafloors shows significant differences compared to the unlayered case. First, there seems to be no unique ansatz. Three different approximations are considered in this article, all reduce to the usual half-space result when layering is absent, and all meet the usual requirements imposed on any small-slope ansatz. Unlike the half-space case, for layered seafloors the SSA does not agree with the Kirchhoff approximation near the specular direction with implications for practical codes that employ the Kirchhoff approximation in this angular regime. The SSA of Ref. 6, here labeled SSL1, shows suspicious behavior in an example with a mud layer overlying a harder basement. Another negative aspect is that it is not presently known how to extend this approximation to multiple layers. The small-slope approach of Berrouk *et al.*^{7} can treat multiple layers and is the most tractable of the three but may oversimplify the problem. In the low-contrast case, where internal reflections are relatively weak, the three approximations give nearly identical results, but this is not true in cases with strong internal reflections.

The differences between SSL2 and SSL3 may result from the way they treat perturbations of the layer interfaces. If a uniform perturbation is made to a flat water-sediment interface, the scattered field will change, due to alterations in the wave interference from the sediment layering. In SSL3, perturbations of this type are treated exactly, since the coordinate space Green's functions are used, and are explicit functions of the layer position, *f*. In SSL2, this type of perturbation results in simple linear changes to the phase of the *T*-matrix. Therefore, SSL3 takes into account changes to wave interference in the sediment due to roughness, whereas SSL2 holds the interference pattern (background plane wave field) constant. This interpretation would result in changes to the peaks and nulls present in the backscattering cross section and coherent reflection coefficient, which are observed in the model comparisons presented in this work. To confirm these speculations, comparison with exact numerical calculations (e.g., using the Monte Carlo approach of Ref. 11 combined with a discretized version of the integral equations presented in Ref. 12) and the various SSAs are needed.

Although two of the approximations can treat scattering by buried interfaces, this article does not provide illustrations of such scenarios. Future work in this direction may highlight differences in the approximations. The primary conclusion of this article is that the two most promising approximations yield different results, and exact calculations are needed to decide which approximation is best.

## ACKNOWLEDGMENTS

This work was supported by the U.S. Office of Naval Research. D.J. benefited from discussions with Dr. D. Tang and Dr. A. Ivakin.

### APPENDIX: FUNCTIONS NEEDED FOR SSL3

This appendix furnishes the functions *q*_{2}, *q*_{3}, and *q*_{4} needed in Eq. (90).

For $\chi 1=1$ and $\chi 2=1$

The following parameters appear in Eq. (A1). For the sake of economy, the same symbols are used as for the parameters of *q*_{1}, even though the two parameter sets are not the same. This practice is followed for *q*_{3} and *q*_{4} as well.

If $\chi 1=\u22121,\u2009V1(K1)$ and $V2(K1)$ are replaced by $V1*(K1)$ and $V2*(K1)$. Likewise, if $\chi 2=\u22121$, $V1(K2)$ and $V2(K2)$ are replaced by $V1*(K2)$ and $V2*(K2)$.

For $\chi 1=1,\u2009\chi 2=1$, and $\chi 3=1$,

The following parameters appear in Eq. (A10):

Conjugation is handled in the same fashion as for *q*_{2}. For example, if $\chi 2=\u22121,\u2009V1(K2)$ and $V2(K2)$ are replaced by $V1*(K2)$ and $V2*(K2)$.

The last function needed is *q*_{4}, given by

The following parameters appear in Eq. (A27):