Coupled Mode Sound Propagation in Inhomogeneous Stratiﬁed Waveguides

: An efﬁcient coupled mode method for modeling sound propagation in horizontally stratiﬁed inhomogeneous waveguides, in which the seabed is modeled as a (layered) acoustic medium, is presented. The method is based on Fawcett’s coupled mode method and the multimodal admittance method. The acoustic ﬁeld is expanded onto the unusual local eigenfunctions composed by normal modes in the corresponding one-layer homogeneous waveguides with constant depth equal to the local total depth of the multilayered waveguide. A set of energy-conserving ﬁrst-order differential equations governing the modal amplitudes of acoustic ﬁelds is derived. The admittance method is employed to solve the differential equations in a numerically stable manna. The coupled mode method considers the backscattering effect of inhomogeneities and full coupling between local modes, and offers improvement from the viewpoint of efﬁciency and computational cost. The acoustic ﬁelds predicted by the method agree well with those computed by the commercial ﬁnite element software COMSOL Multiphysics. The method can be extended to further establish fast and accurate 3D sound propagation models in complex shallow water environments.


Introduction
Sound propagation in shallow water in the presence of inhomogeneities has been studied for several decades, conventionally modeled by the nonseparable Helmholtz equation in horizontally stratified inhomogeneous waveguides [1][2][3][4][5]. The step-wise coupled mode method is commonly used for sound propagation in waveguides in which a step-wise variation of the properties with range occurs [3,4]. However, for the sound propagation problem in waveguides with continuously varying properties, the stepwise approximation of the properties renders the numerical solution procedure cumbersome [6]. The continuous coupled mode method is a classical approach to sound propagation problems in range-dependent waveguides, originally presented by Pierce [1] and Milder [2]. They expanded the sound field in terms of usual local modes with range-dependent modal amplitudes and derived a system of coupled mode equations containing two coupled coefficients. Rutherford and Hawker indicated that the Pierce-Milder coupled mode equations do not conserve the energy among modes because the usual local modes satisfy different boundary conditions than the acoustic pressure [7]. Fawcett derived a coupled system of differential equations containing two coupling matrices and two interface matrices, and treated this coupled mode system by a finite difference scheme [8]. It is noted that, the coupling matrices and interface matrices, which require a knowledge of derivatives of usual local modes, can only be approximately computed [9]. Fawcett's approach considers the backscattering effect from range-dependent environments and full coupling between modes, and conserves the energy among modes. However, the complexity involved in computing coupling matrices and interface matrices along with the complex numerical scheme for solving differential equations makes this approach impractical, especially in long range and high frequency applications, where the computational requirements increase substantially. In order to reduce computational complexity, several asymptotic coupled mode methods have been developed. The adiabatic coupled mode method neglects mode coupling among modes and provides accurate predictions when the range dependence of a medium is sufficiently weak [5,10]. The one-way approximation coupled mode method is efficient for sound propagation in waveguides where the backscattering effect can be ignored [5,11]. However, asymptotic coupled mode models can provide very poor solutions when the range dependence of a medium is significant. Thus, it is still very important to study an efficient two-way coupled mode method for sound propagation in complex underwater waveguides.
The multimodal admittance method presented by Pagneux [12,13] is a two-way coupled mode method for wave propagation in waveguides with range-dependent crosssections. This method basically consists in rewriting the Helmholtz equation into a set of first-order evolution equations and then projecting onto the usual local modes to obtain a set of first-order differential equations governing the modal amplitudes of the wave fields. In the case where these differential equations are to be integrated numerically, the problem of numerical divergence owing to the presence of evanescent modes appears. The admittance matrix, namely, the Dirichlet-to-Neumann operator in modal domain, is introduced to avoid numerical problems. Using the admittance matrix, first-order differential equations are reformulated into the Riccati equation governing the admittance matrix and the first-order evolution equation governing the modal amplitudes of acoustic pressure. Both the Riccati equation and the evolution equation can be computed with classical numerical schemes [12][13][14][15]. This method is efficient and numerically stable. Liu and Li [16] derived a coupled system of first-order differential equations for sound propagation in simple one-layer waveguides with range-dependent environments based on Fawcett's coupled mode method and the multimodal admittance method. We extend the method here to present an efficient coupled mode scheme for sound propagation and wave scattering problems in complex inhomogeneous stratified waveguides. This paper focuses on the two-way solutions in horizontally stratified waveguides with both bathymetric and volumetric variations. In particular, waveguides with penetrable scatterers are also inspected, modeled by horizontally stratified waveguides with intersecting interfaces. We note that the paper deals with acoustic waveguides in which the seabed is also modeled as a (layered) acoustic medium. In many applications, the modeling of seabed needs to account for the shear rigidity (seabed elasticity); see, e.g., [17] and the references cited therein. In an attempt to reduce the complexity involved in the computation of usual local modes and their derivatives for complex multilayered environments, the local modes for expansion of acoustic pressure in this paper are defined by the transverse eigenfunctions in a one-layer waveguide with constant depth equal to the local total depth of the multilayered waveguide and constant physical parameters. We note that such local modes and their derivatives can be analytically computed. A set of second-order coupled mode equations governing the modal amplitudes is derived by projecting the local mode onto the Helmholtz equation. The energy conservation is guaranteed by introducing interface matrices. Remarkably, we introduce a set of unknown quantities with respect to the modal amplitudes and their horizontal derivatives. Using this unknown vector, the second-order coupled mode equations are reduced into a set of first-order differential equations with respect to the modal amplitudes and the unknown quantities. This first-order coupled system preserves the energy conservation property, but no interface matrix is contained. The required coupling matrices can be computed with numerical integral methods since their integrands are analytical. Moreover, it is straightforward to implement the admittance method for numerically stable solutions. The use of the unknown vector substantially increases the efficiency of the method. Since a more realistic representation of the ocean bottom is a penetrable infinite half space, the application of our coupled mode method to waveguides with a penetrable infinite half space is studied on basis of the elimination of branch cut and branch line integrals [18]. Numerical results for sound propagation in inhomogeneous multilayered waveguides predicted by commercially available finite element software COMSOL multiphysics are also presented to validate the efficiency and accuracy of this method.
The remainder of the paper is organized as follows: Section 2 presents the basic formulation of two-way and one-way asymptotic coupled mode methods for sound propagation in horizontally stratified inhomogeneous waveguides. The sound field expression generated by an arbitrary incident wave in such an environment is provided. By constructing mapping between incident waves and wave fields, optimal incident waves resulting in sound self-focusing at any position [19] in the multilayered waveguide are studied. Section 3 illustrates the sound fields predicted by our coupled mode method and COM-SOL for a complex multilayered waveguide, and the sound self-focusing pattern in this numerical example is shown. Section 4 presents the extension of the present coupled mode formalism to waveguides with the effect of internal waves and a penetrable bottom. Section 5 gives the conclusions.

Two-Dimensional Two-Way Coupled Mode Equations
Let us start with sound propagation in the multilayered waveguide sketched in Figure 1 with interfaces h j (x), j = 1, 2, . . . , J − 1, separating the layers. The upper water column is confined between a flat pressure-release surface and underlying sediment layers of any shape. The Neumann boundary condition is imposed at the bottom of the lower layer. The whole region is divided into three subregions. Subregion 0 ≤ x ≤ x R corresponds to the scattering region where both volumetric and bathymetric variations exist. In subregions x < 0 (not shown) and x > x R , mass density and sound speed are assumed to be only depth-dependent, and boundaries and interfaces are horizontal. This assumption allows for consistently formulating the source and radiation conditions by means of normalmode expansions.
supplemented by boundary conditions and continuous conditions where j = 1, 2, . . . , J − 1, c(x, z) is sound speed, and ρ(x, z) is density. c(x, z) and ρ(x, z) could exhibit sharp discontinuities at the interfaces. According to the multimodal method, acoustic pressure can be expressed as a sum of transverse eigenfunctions ψ n (z; x): where N is the truncation number, P n (x) are modal amplitudes, and ψ n (z; x) are defined to be normal modes of a one-layer homogeneous waveguide with constant depth equal to the local depth of the multilayered waveguide, satisfying the following equations: Eigenfunctions ψ n (z; x) in the series representation have analytical expressions, i.e., In contrast, the usual local modes that traditional coupled mode method used in the series representation for acoustic pressure, namely, transverse modes at (x, z) of the complex multilayered waveguide, can only be numerically computed. The use of eigenfunctions ψ n (z; x) simplifies computations in this sense.
Then, a projection of the Helmholtz equation onto ψ m (z; x) gives By substituting Equation (4) into Equation (7) and applying Fawcett's approach [8], one can obtain second-order coupled-mode equations where A mn , K mn , and L mn are correction coefficients resulting from the use of unusual eigenfunctions ψ n (z; x); B mn , C mn , F mn , and G mn are mode coupling coefficients; and D mn and E mn are interface coefficients given by, respectively, Since ψ n (z; x) satisfy different boundary and interface conditions than acoustic pressure when the bottom boundary and interfaces are range dependent, interface coefficients D mn and E mn are introduced in the derivation to guarantee energy conservation among modes by properly applying the boundary and interface conditions. Details of the derivation of D mn and E mn are presented in Appendix A. The interface coefficients can be analytically computed. Other coefficients can be numerically computed by numerical integration method. In this paper, we use the Clenshaw-Curtis method [20] to calculate the correction and coupling coefficients.
Although the fully two-way coupled system Equation (22) conserves the energy among modes, as any classical coupled mode method, the local modal series representation Equation (4) has a slow rate of convergence with the modal amplitudes following the order O(n −2 ), where n is the mode number. The poor convergence is attributed to the differences in the interface and boundary conditions satisfied by ψ n (z; x) and acoustic pressure. There are several ways that allow for accelerating the convergence of the local mode series [6,[21][22][23][24]. One of them [22,23] uses the idea of boundary modes. It is based on an enriched modal series representation in terms of supplementary boundary modes, which are orthogonal to the local eigenfunctions but satisfy different boundary and interface conditions. Boundary modes have the beneficial effect of acting as evanescent modes and do not change the form of the coupled mode equations. However, the numerical implementation of the boundary modes for multilayered waveguides is complicated and generally requires an amount of computing power. We therefore only consider the nonimproved version of coupled mode equations for sound propagation analysis in multilayered environments.
In order to solve the second-order coupled mode equations, we define a set of unknown quantities S m , m = 0, 1, . . . , N − 1, by The x-derivative of S m gives where By substituting Equations (18) and (19) into Equation (8), and rewriting the resulting equation into matrix form, we have where P = (P n ), S = (S m ), and the elements in matrices A, C, L, and K are A mn , C mn , L mn and K mn , respectively. In Equation (20), we use the identical relation Substituting Equation (21) into the expressions of C and H yields C = AQ and H = Q T AQ. Thus, one can obtain H = C T A −1 C.
From Equations (18) and (20), we obtain the first-order coupled mode equations d dx We note that, to calculate acoustic fields, only four of the nine matrices in Equation (8) remain to be computed. The first-order coupled system Equation (22) preserves the energy conservation property but does not contain any interface matrix. For details of energy conservation of Equation (22), see Appendix B. It is the proper application of S that gives the coupled mode equations a simple form and substantially reduces computing power, especially in long-range application.
The first-order differential equations account for source and radiation conditions. Directly integrating Equation (22) starting from the radiation condition may encounter the problem of numerical divergence owing to the exponential growth of evanescent modes. In order to avoid numerical problems, we introduce admittance matrix Y and propagator matrix M by S(x) = Y(x)P(x) and P(x) = M(x)P(0), respectively. Y satisfies a Riccati equation and M a first-order differential equation coupled to Y: Y can be computed by integrating from right to left with the initial value Y( are horizontal wavenumbers in region x > x R (for more details, see Appendix C). M(x) can be solved from left to right with the initial value M(0) = I, where I is the identity matrix. In this paper, numerical integration of Equations (23) and (24) is performed using the Magnus method [25,26]. Further, the reflection and transmission matrices that characterize the scattering properties of the scattering region can be computed. Reflection matrix R is defined by P r (0) = RP i (0), and transmission matrix T is defined by P t (x R ) = TP i (0), where P i , P r , and P t denote the modal amplitudes of the incident, reflected and transmitted waves, respectively. We obtain and where Y(0) is the admittance matrix at x = 0, and Y 0 = A(0) L(0) − K(0) is the local admittance matrix at x = 0. Eventually, one can write the acoustic pressure in the multilayered waveguide as where Ψ is a vector with elements being ψ n . To summarize, once the geometry and medium parameters of the waveguide are identified, admittance matrix Y and propagator matrix M are available for all x. Then, reflection matrix R and transmission matrix T can be deduced. Finally, the sound field generated by any incident wave can be obtained. Although the derivation of the coupled mode system in the paper is on basis of multilayered waveguides in which the properties are continuously varying with range in each layer, the method can be naturally extended to solve sound propagation problems for waveguides in which sharp discontinuities of the properties with range exist; see, e.g., [27][28][29].

Two-Dimensional One-Way Coupled Mode Equations
For waveguides where the range dependence of medium is weak, one-way approximation coupled-mode formulation provides a speed-up way to solve acoustic pressure by neglecting backscattering energy [30]. This means that elements in reflection matrix R are relatively small in this case. From Equation (25), it follows directly that the admittance matrix at any position approximates the local admittance matrix, and the iterative computation of admittance matrix is no longer necessary. Thus, the admittance matrix takes the following form:Ỹ The one-way approximation coupled-mode equation reads with initial condition P i (0) = P + (0) and P + (0) representing the modal amplitudes of forward propagating wave. Similarly, we introduce propagator matrixM(x) by P(x) =M(x)P(0).M(x) satisfies the equationM = (−A −1 C + A −1Ỹ )M with the initial valueM(0) = I. Then, acoustic pressure under one-way approximation can be written as

Optimal Incident Waves for Sound Focusing
As a byproduct of our coupled mode algorithms, the optimal incident wave that leads to sound focusing at a certain position in a multilayered waveguide is investigated.
From Equation (27), the acoustic pressure at any position (x 0 , z 0 ) can be rewritten as where < a 1 , a 2 >= a † 1 a 2 denotes the scalar product of vectors in Hilbert space, and ' † ' denotes the conjugate transpose. The energy flux of the incident wave is assumed to be unity. Then, acoustic pressure modulus at (x 0 , z 0 ) can reach its maximum value, max p i |p(x 0 , z 0 )|, when the modal amplitudes of the incident wave take the following form: where Λ is fixed by the normalization of the sound energy flux of incident wave E pi = Im(

Numerical Results
In this section, we validate the present coupled mode method by comparing the results for a numerical example with those obtained using the acoustics module of the commercially available finite element software COMSOL Multiphysics and show the sound focusing pattern by injecting optimal incident waves.
Let us consider sound propagation in a two-layer waveguide with a penetrable scatterer being located in the upper water column. We model such an inhomogeneous waveguide as a virtual four-layer waveguide with intersecting interfaces. In scattering region (0 < x < 400 m), sound speed and mass density are c 1 = (1500 − 0.1z + 0.1x) m/s, ρ 1 = (1000 + 0.1z + 0.1x) kg/m 3 , c 2 = (1700 − 0.2z + 0.1x) m/s, and ρ 2 = (1500 + 0.2z + 0.1x) kg/m 3 h 2 (x) = 100 − 10(1 − cos(πx/400)), 0 < x < 400, 80, The waveguide is excited by a distributed source p i (z) = sin(1.5πz/100) from the left at frequency f = 100 Hz. There are 10 propagating modes in the left lead and 12 propagating modes in the right lead. From radiation condition Y(x = 400), Y, M, R, and T are available. Then, acoustic pressure can be computed using Equation (27) with the source condition P i = (0, 1/ √ 0.02, 0, . . . , 0) T . Figure 2a,b, respectively, show the acousticpressure moduli calculated by the present two-way coupled-mode method (CMM) and COMSOL. The results were computed on a quad-core 3.6GHz processor running the 64-bit Windows operating system. The truncation number used for the local series representation is N = 20. The average computation time per output grid [31] is 1.4 × 10 −4 s. Clearly, the sound field wildly oscillates in the scattering region, which is expected for the case where the backscattering energy is large and the mode coupling effect is significant. Figure 2c illustrates the acoustic pressure distributions in the horizontal direction at a depth z = 20 m predicted by the two-way CMM, one-way CMM and COMSOL. The spatial discretization and truncation numbers used in the two-and one-way CMM are identical. We see that the calculations by the two-way CMM and COMSOL are in excellent agreement. The slight differences are attributed to the truncation of the infinite modal series Equation (4), and the interpolation difference between COMSOL and our code. It can be minified, of course, by increasing truncation number N and improving mesh quality. However, the results predicted by the one-way CMM are remarkably different, meaning that backscattering energy cannot be neglected in this environment. Since the modal method only necessitates the trace of the the modal amplitudes of the first N local basis functions, the present coupled mode algorithm is much more efficient than COMSOL is, which is a direct Helmholtz solver, for sound propagation over long-range distances. Next, the acoustic wave self-focusing pattern in a multilayered waveguide is presented. The geometry and medium parameters of the waveguide are kept the same as those in Figure 2; the frequency is f = 100 Hz. Figure 3a shows wave focusing at (450, 20) m in the transmission region by injecting the incident wave of Figure 3b, which is computed from Equation (32). Figure 3c plots the corresponding distribution of acoustic pressure along the z direction at x = 450 m. We see that energy in the vicinity of the focus occupies the vast majority of energy along the z direction. This method allows one to enhance the energy density of any position in a multilayered waveguide.

Application to Waveguides with Penetrable Infinite Half Space and Internal Wave
In this section, we apply the coupled mode formalism to sound propagation problems in a multilayered waveguide, considering the effect of a penetrable infinite half space and internal waves, which is a common scenario in real ocean environments.
In previous sections, we assumed that the bottom is acoustically rigid. However, a more realistic representation of the ocean bottom is an infinite half space that allows for energy penetrating into the half space. The rigid bottom assumption neglects one or more components of sound fields depending on the choice of branch cut, for instance, leaky modes plus an integral along the Pekeris cut for Pekeris branch cut, or an integral of the continuous spectrum along the EJP cut for EJP branch cut. Difficulties in the application of bounded coupled mode models to a Pekeris waveguide consist of the discrete approximation of continuous spectrum. An available algorithm eliminates the branch cut and branch line integrals by introducing a small complex sound speed gradient in the half space [18]. The branch cut contribution is represented by discrete modes. This method has the advantage that leaky modes do not increase in amplitude with increasing depth. Since branch line modes make a minor contribution to the sound field in the water column, it is possible to ignore the branch line modes in numerical implementation [32]. Therefore, it provides a straightforward way to incorporate with the couple mode algorithms. We extend this approach here to coupled mode sound propagation in multilayered waveguides.
We consider the waveguide shown in Figure 4. A two-layer system of sound speed in the water column is used. The upper mixed layer has sound speed c ml and water density ρ 1 . The lower water layer has sound speed c 1 and the same density ρ 1 . The shape of the internal solitary wave is formulated by where ∆h is the amplitude of the internal wave, x wave denotes the wave center, and L is the width parameter [33]. The ocean bottom is characterized by a penetrable infinite half space with sound speed c 2 and density ρ 2 . To apply the coupled mode formalism, the penetrable infinite half space is approximately reproduced by a bounded layer. The corresponding sound speedc 2 has a small gradient,c where the attenuation α is in Gaussian form, α = α 0 e −(z−h 2 (x)) 2 /(2λ 2 2 ) , α 0 is in the unit of dB/λ 2 , λ 2 = c 2 / f , and α c = 40πlog 10 e. The thickness of the lower bounded layer, h 2 − h 1 (x), is meant to be larger than 3λ 2 , such that attenuation increases gently in the beginning, then grows rapidly with increasing depth, and finally smoothly reaches its maximum at z = h 2 . We use such a form of α in an attempt to reduce the spurious reflection that can return to the water column and provide accurate solutions in the water column using an acceptable truncation depth.
where x R = 16 km. The water-sediment interface is sloping, In numerical implementation, a false rigid boundary is placed at h 2 = 200 + 5λ 2 beneath the water column, and the sound speed in the bottom layer is replaced byc 2 (Equation (36)) with α 0 = 10 dB/λ 2 . The point source is located in (x s , z s ) = (0, 30) m at a frequency of 50 Hz. The source condition, which is represented by the modal amplitudes of the Green's function in the waveguide of Figure 4, is [20] where Ψ is a vector with elements being functions ψ n , and Y x s = A(x s ) L(x s ) − K(x s ) .
The acoustic pressure can be computed using Equation (27) with the source condition and radiation conditions. In addition, the eigenvalues of which fall in the interval ω/c 2 , max(ω/c 1 (z)) , correspond to the horizontal wavenumbers of the trapped modes. At the chosen frequency, there are six trapped modes at x = 0 and four trapped modes in the transmission region. Figure 5a presents the computed sound field using the present two-way coupled mode method. The truncation number is N = 30, which has been proved to be enough for numerical convergence. Figure 5b compares the transmission loss predicted by the two-way CMM, one-way CMM, and COMSOL. The numerical model used for COMSOL is a waveguide with a PML beneath the water column. The receiver depth is z r = 60 m. The transmission loss is defined by We see that the agreement between the two-way CMM and COMSOL is excellent. Differences are attributed to the deviations in spatial discretization. The results of the one-way CMM are accurate for the range internal from 0 m to about 6 km, implying that the backscattered energy is small in this region. As sound propagates through the internal wave, waves interact with the sloping interface and internal waves, and differences in the results between the one-and two-way CMM increase gradually. Since the computation model used in COMSOL mimics the real waveguide with a penetrable infinite half space, one can conclude that acoustic pressure in the water column of a waveguide with a penetrable infinite half space and internal waves predicted by our two-way coupled mode method is accurate.

Conclusions
An efficient coupled mode method is presented for sound propagation through multilayered inhomogeneous waveguides based on Fawcett's coupled mode method and the multimodal admittance method. The inhomogeneous environments include rangedependent physical parameters and range-dependent interfaces (boundaries). The method uses the unusual local transverse eigenfuctions in the series representation of the acoustic pressure, which are defined to be the transverse modes in a one-layer homogeneous waveguide with constant depth equal to the local total depth of the multilayered waveguide and constant physical parameters. The unusual local modes have analytical solutions. Projecting the Helmholtz equation onto the unusual local modes yields a set of second-order differential equations with respect to the modal amplitudes of acoustic pressure. Energy conservation among modes is guaranteed by the introduction of two interface matrices in the case of a lossless acoustic waveguide. By properly applying a set of unknown quantities S in terms of the modal amplitudes of acoustic pressure and their horizontal derivatives, the second-order coupled mode equations are reduced into a set of first-order differential equations. This first-order coupled system preserves the energy conservation property but requires no knowledge of interface matrices. The admittance method was employed to solve the first-order differential equations for avoiding numerical problems. The reflection and transmission matrices that characterize the scattering properties are also given. Our coupled mode method considers the full coupling between modes and the backscattering effect and substantially reduces computing power compared with fully numerical methods, especially in long-range and high-frequency applications. Numerical results show that our two-way coupled mode method provides us with accurate solutions for multilayered inhomogeneous waveguides with acoustically rigid bottom or a penetrable infinite half space. Our coupled system is appropriate for two-way solutions in long-range/high-frequency propagation and scattering problems and can be naturally extended to treat sound propagation through real ocean environments.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Interface Matrices D and E
Here, we present the derivation of D mn and E mn . Integrating by part for the projection of the term ∂/∂z(ρ −1 ∂p/∂z) in Equation (1), and using the following boundary and interface conditions, where j = 1, 2, . . . , J − 1, one can obtain where w i,n are the modal amplitudes, and ψ n (z) are transverse normal modes in a homogeneous waveguide with a constant depth equal to h J , i.e., ψ n (z) = 2/h J sin (n + 1/2)πz/h J . Applying the operator where the elements in W are w i,n , and K x is a diagonal matrix with diagonal elements being k xi . Therefore, the eigenvalue problem governed by Equation (A10) is transformed into a matrix eigenvalue problem of Equation (A12). The axial wavenumbers k xi and the modal amplitudes w i,n can be obtained by eigen-decomposing the matrix A −1 (K − L).
Alternatively, for environments of which the properties are invariant in the horizontal, the coupling among modes is lost. The second-order coupled mode equation is then reduced to P = A −1 (L − K)P.
Since S = YP and P = A −1 S, the admittance matrix in the transmission region takes the following form Y = A(L − K). (A14)