# Coupled Mode Sound Propagation in Inhomogeneous Stratified Waveguides
## 3 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, \(\rho_{1}=(1000+0.1z+0.1x)\) kg/m\({}^{3}\), \(c_{2}=(1700-0.2z+0.1x)\) m/s, and \(\rho_{2}=(1500+0.2z+0.1x)\) kg/m\({}^{3}\). In the transmission region (\(x\geq 400\) m), parameters are \(c_{1}=(1540-0.1z)\) m/s, \(\rho_{1}=(1040+0.1z)\) kg/m\({}^{3}\), \(c_{2}=(1740-0.2z)\) m/s, and \(\rho_{2}=(1540+0.2z)\) kg/m\({}^{3}\). The scatterer has 20m radius, (\(x_{s}=350,z_{s}=25\)) m center, and is fluid-filled with \(c_{s}=1800\) m/s and \(\rho_{s}=1600\) kg/m\({}^{3}\). The shapes of the interface and bottom are formulated by Equations (33) and (34), respectively,  
\[h_{1}(x)= \begin{cases}40+10(1-\cos(\pi x/400)),&0<x<400,\\ 60.&x\geq 400.\end{cases}  \tag{33}\]  
\[h_{2}(x)= \begin{cases}100-10(1-\cos(\pi x/400)),&0<x<400,\\ 80,&x\geq 400.\end{cases}  \tag{34}\]  
The waveguide is excited by a distributed source \(p_{i}(z)=\sin(1.5\pi 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, M, R, and T are available. Then, acoustic pressure can be computed using Equation (27) with the source condition \(\mathbf{P}_{i}=(0,1/\sqrt{0.02},0,\ldots,0)^{\mathrm{T}}\). Figure 2a,b, respectively, show the acoustic-pressure moduli calculated by the present two-way coupled-mode method (CMM) and COMSOL. The results were computed on a quad-core 3.6 GHz 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\times 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 3 shows wave focusing at \((450,20)\) m in the transmission region by injecting the incident wave of Figure 3, which is computed from Equation (32). Figure 3 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.