⇦Chinese Physics Letters, 2017, Vol. 34, No. 7, Article code 074302 Coupled Perturbed Modes over Sloping Penetrable Bottom Fei-Long Zhu(朱飞龙)1,2, Eric I. Thorsos3, Feng-Hua Li(李风华)1** Affiliations 1State Key Laboratory of Acoustics, Institute of Acoustics, Chinese Academy of Sciences, Beijing 100190 2University of Chinese Academy of Sciences, Beijing 100049 3Applied Physics Laboratory, University of Washington, WA 98105, USA Received 24 February 2017 **Corresponding author. Email: lfh@mail.ioa.ac.cn Citation Text: Zhu F L, Thorsos E I and Li F H 2017 Chin. Phys. Lett. 34 074302 Abstract In environments with water depth variations, one-way modal solutions involve mode coupling. Higham and Tindle developed an accurate and fast approach using perturbation theory to locally determine the change in mode functions at steps. The method of Higham and Tindle is limited to low frequency ($\le$250 Hz). We extend the coupled perturbation method, thus it can be applied to higher frequencies. The approach is described and some examples are given. DOI:10.1088/0256-307X/34/7/074302 PACS:43.30.Bp, 43.20.Bi, 43.20.Mv © 2017 Chinese Physics Society Article Text In environments with water depth variations, normal mode solutions involve mode coupling. One approach is to model the depth variation as a series of steps, and uses a new set of modes at each step. The sound field is determined by matching the boundary conditions at each step interface.[1] Higham and Tindle (H&T)[2] developed an accurate and fast approach, which uses perturbation theory to locally determine the change in mode functions and mode coupling at steps, and this limits the required number of new sets of mode functions. The method of H&T is limited to low frequency ($\le$250 Hz), while the method described here can be extended to higher frequency. When the normal modes are known for a specific environment, perturbation theory can be used to solve for the new normal modes for a slightly different environment. This approach can be employed when normal modes are used for range-dependent environments. This method is relatively straightforward when the only range dependence is a changing sound speed profile.[3-5] For the case of a range-dependent depth, a modification of the modal equation is required to account for the change in density profile over a sloping bottom. In a horizontal stratified environment, the mode functions $U_n (z)$ are the solutions of the depth-separated modal equation[6] $$\begin{alignat}{1} \rho (z)\frac{d}{dz}\Big[{\frac{1}{\rho (z)}\frac{dU_n (z)}{dz}}\Big]+[{k^2(z)-k_n^2 }]U_n (z)=0.~~ \tag {1} \end{alignat} $$ In the work described here, a water column overlies a bottom of finite depth, and thus only discrete modes are present, denoted by the integer index $n$. In general, the densities in water and the sea floor bottom are different, and the assumptions of a constant density in the water column and piecewise constant density in the bottom are often used in acoustic models. Here the simplest assumption is made: $\rho (z)=\rho_{\rm w}$ for $z < h$, and $\rho (z)=\rho_{\rm b}$ for $h < z < H$ ($\rho_{\rm w} \ne \rho_{\rm b}$), where $H$ is the total depth, $\rho_{\rm w}$ and $\rho_{\rm b}$ are constant densities in the water and bottom, respectively. Also in Eq. (1) $k(z)$ is the acoustic wave number. The bottom is assumed to include absorption. Thus the mode functions, the mode eigenvalues $k_n^2$, and $k(z)$ for $h < z < H$ are complex quantities. For many applications using modes, continuous modes can be assumed unimportant and are ignored. However, for propagation over a sloping bottom, it is found that coupling into and out of the continuous modes is important. By taking the vertical domain to be finite, the continuous modes are approximated here by a set of closely spaced discrete modes. The gradients of the mode functions are discontinuous at the water and bottom interface, and the depth of the 'kink' in the mode functions will vary with the water depth, which is difficult to handle by the perturbation method. To remove the gradient discontinuity of the mode functions, a transformation made on the $z$ coordinate is given by $$ \tilde {z}=\begin{cases} \!\! z,&0\leqslant z < h, \\\!\! h+(\rho_{\rm b} /\rho_{\rm w})(z-h),&h < z < H. \end{cases}~~ \tag {2} $$ Transforming Eq. (1) from $z$ to $\tilde {z}$ yields $$\begin{align} \frac{d^2\tilde {U}_n (\tilde {z})}{d\tilde {z}^2}+X(\tilde {z})[{\tilde {k}^2(\tilde {z})-k_n^2 }]\tilde {U}_n (\tilde {z})=0,~~ \tag {3} \end{align} $$ where $X(\tilde {z})=1$ for $z < h$, $X(\tilde {z})=(\rho_{\rm w} /\rho_{\rm b})^2$ for $h < z < H$, $\tilde {k}(\tilde {z})=k[{z(\tilde {z})}]$, and $\tilde {\rho }(\tilde {z})=\rho [{z(\tilde {z})}]$. For an eigenvalue equation of the form of Eq. (3) when there is no attenuation in the bottom and when $\tilde {\rho }(\tilde {z})$ and $\tilde {k}(\tilde {z})$ are continuous functions of depth (not the case being considered here), it is well known[7] that the operator in Eq. (3) is self-adjoint, and the mode functions and eigenvalues are real. Also the mode functions $\tilde {U}_n (\tilde {z})$ and $\tilde {U}_m (\tilde {z})$ are orthogonal and satisfy $$\begin{align} \int\limits_0^{\tilde {z}(H)} {X(\tilde {z})\tilde {U}_n (\tilde {z})\tilde {U}_m (\tilde {z})d\tilde {z}=\delta _{nm}},~~ \tag {4} \end{align} $$ and form a complete set so that an arbitrary function $f(\tilde {z})$ on $(0,\tilde {z}(H))$ can be expanded as $$\begin{align} f(\tilde {z})=\sum\limits_m {a_m} \tilde {U}_m (\tilde {z}).~~ \tag {5} \end{align} $$ For the more general case where $\tilde {\rho }(\tilde {z})$ and $\tilde {k}(\tilde {z})$ are real but may be only piecewise continuous (e.g., with jumps at the water-bottom interface), the above statements remain true,[8] including Eqs. (4) and (5). When the absorption is included in the bottom, the mode functions are complex, yet Eq. (4) remains true as a bi-orthogonal condition[9,10] without any complex conjugation. For the examples considered here, $\tilde {\rho }(\tilde {z})$ and $\tilde {k}(\tilde {z})$ are piecewise continuous, and absorption is included in the bottom with the result that the mode functions are complex. Equation (4) has been verified numerically, and it will be assumed that Eq. (5) remains true. The final results support this assumption. The above expressions can be applied directly to a range-independent environment. When the depth changes with range, the depth change is modeled as a series of small steps in the water-bottom interface. The perturbation scheme presented here allows new complex mode wave numbers and new complex mode functions to be obtained at a series of steps, reducing the number of times that new mode functions need to be obtained with a full eigenfunction solution. As the sound propagates forward in a small range increment that crosses a step in depth, the changes in sound speed and density profiles lead to the known changes $\tilde {k}^2(\tilde {z})\to \tilde {k}^2(\tilde {z})+\Delta \tilde {k}^2(\tilde {z})$ and $X(\tilde {z})\to X(\tilde {z})+\Delta X(\tilde {z})$, then the mode functions change to $\tilde {U}_n (\tilde {z})+\Delta \tilde {U}_n (\tilde {z})$ and the horizontal wave numbers for each mode change to $k_n +\Delta k_n$ according to perturbation theory. The new values must satisfy Eq. (3) as well, namely, $$\begin{alignat}{1} \!\!\!\!\!\!\!\!&\frac{d^2(\tilde {U}_n (\tilde {z})+\Delta \tilde {U}_n (\tilde {z}))}{d\tilde {z}^2}+(X(\tilde {z})+\Delta X(\tilde {z}))[\tilde {k}^2(\tilde {z})\\ \!\!\!\!\!\!\!\!&+\Delta \tilde {k}^2(\tilde {z})-(k_n +\Delta k_n)^2](\tilde {U}_n (\tilde {z})+\Delta \tilde {U}_n (\tilde {z}))=0.~~ \tag {6} \end{alignat} $$ In solving for these changes, H&T used Eq. (5) to express the change in the mode functions as $$\begin{align} \Delta \tilde {U}_n (\tilde {z})=\sum\limits_m {a_{nm} } \tilde {U}_m (\tilde {z}),~~ \tag {7} \end{align} $$ with $a_{nm}$ becoming unknowns to solve for. H&T just kept first-order terms, which yields first-order results for $\Delta k_n$ and $a_{nm}$ in Eqs. (15) and (16), respectively.[2] However, in general $\Delta X$ is almost the same order of $X$, thus the term $\Delta X\Delta \tilde {k}^2(\tilde {z})$ should be kept as well. Therefore, better expressions than given in Ref. [2] are the following $$\begin{align} \Delta k_n =\,&\Big[\int\limits_0^{\tilde {z}(H)} (X\Delta \tilde {k}^2+\Delta X\tilde {k}^2+\Delta X\Delta \tilde {k}^2\\ &-\Delta Xk_n^2)\tilde {U}_n^2 d\tilde {z}\Big]/(2k_n),~~ \tag {8}\\ \end{align} $$ $$\begin{align} a_{nm} =\,&\frac{1}{k_n^2 -k_m^2 }\int\limits_0^{\tilde {z}(H)} [X\Delta \tilde {k}^2+\Delta X\tilde {k}^2+\Delta X\Delta \tilde {k}^2\\ &-\Delta Xk_n^2]\tilde {U}_n (\tilde {z})\tilde {U}_m (\tilde {z})d\tilde {z},~(n\ne m),~~ \tag {9} \end{align} $$ $$\begin{align} a_{nn} =\,&-\frac{1}{2}\Big[\int\limits_0^{\tilde {z}(H)} {\Delta X\tilde {U}_n^2 d\tilde {z}} +\sum\limits_{m\ne n}(a_{nm}^2\\ &+2a_{nm} \int\limits_0^{\tilde {z}(H)}\Delta X\tilde {U}_n \tilde {U}_m d\tilde {z})\Big]\Big/\Big[1\\ &+\int\limits_0^{\tilde {z}(H)} {\Delta X\tilde {U}_n^2 d\tilde {z}}\Big].~~ \tag {10} \end{align} $$ Even with these improvements the low-order approximation results are only effective at relatively low frequencies ($\leqslant 250$ Hz), where there are few trapped modes. As the frequency is increased, the results become divergent over a series of steps. The result at higher frequencies can be improved significantly by keeping all the expansion terms in Eq. (6), with the result $$\begin{alignat}{1} \Delta k_n =\,&\frac{J_n +\sum\limits_m {a_{nm} T_{nm} } }{2k_n (1+a_{nn})},~~ \tag {11} \end{alignat} $$ $$\begin{alignat}{1} a_{nm} =\,&\frac{T_{nm} +\sum\limits_p {a_{np} T_{pm} } }{({k_n +\Delta k_n })^2-k_m^2 }({n\ne m}),~~ \tag {12} \end{alignat} $$ $$\begin{alignat}{1} a_{nn} =\,&-\frac{1}{2}\Big[\int\limits_0^{\tilde {z}(H)} {\Delta X\tilde {U}_n^2 d\tilde {z}} +\sum\limits_m (a_{nm}^2 \\ &+2a_{nm} \int\limits_0^{\tilde {z}(H)} {\Delta X\tilde {U}_n \tilde {U}_m d\tilde {z}} ) \\ &+\sum\limits_p \sum\limits_q a_{np} a_{nq} \int\limits_0^{\tilde {z}(H)} \Delta X\tilde {U}_p \tilde {U}_q d\tilde {z}\Big].~~ \tag {13} \end{alignat} $$ In these equations $$\begin{align} J_n =\,&\int\limits_0^{\tilde {z}(H)} [X\Delta \tilde {k}^2+\Delta X\tilde {k}^2+\Delta X\Delta \tilde {k}^2\\ &-\Delta X(k_n +\Delta k_n)^2]\tilde {U}_n^2 d\tilde {z},~~ \tag {14} \end{align} $$ $$\begin{align} T_{nm} =\,&\int\limits_0^{\tilde {z}(H)} [X\Delta \tilde {k}^2+\Delta X\tilde {k}^2+\Delta X\Delta \tilde {k}^2\\ &-\Delta X(k_n +\Delta k_n)^2]\tilde {U}_n \tilde {U}_m d\tilde {z}.~~ \tag {15} \end{align} $$ The low-order approximation given by Eqs. (8)-(10) ignores the high-order terms, such as $\sum\limits_m {a_{nm} T_{nm}}$ and $\sum\limits_p {a_{np} T_{pm}}$, which become important as the frequency is increased, since the coupling between the nearby modes at high mode number is stronger. In Eqs. (11)-(13), terms on the left side also appear as unknowns on the right hand side. An iteration algorithm has been developed to solve these equations, and the low-order approximation results from Eqs. (8)-(10) are used as the initial values. After iteration, the perturbed approximation results have converged to the exact solutions of Eqs. (11)-(15). In the one-way coupled mode theory, the method referred to as the single scatter approximation agrees almost perfectly with the exact solution.[11] In the single scatter approximation, the mode coefficients in the $(j+1)$th step are $$\begin{align} A_n^{j+1} =\sum\limits_m {C_{nm}^j A_m^j e^{ik_n^j (r_j -r_{j-1})}},~~ \tag {16} \end{align} $$ where $$\begin{align} C_{nm}^j =\,&\frac{1}{2}\Big[\int\limits_0^H {\frac{U_n^{j+1} (z)U_m^j (z)}{\rho ^{j+1}(z)}} dz\\ &+\frac{k_m^j }{k_n^{j+1}}\int\limits_0^H {\frac{U_n^{j+1} (z)U_m^j (z)}{\rho ^j(z)}dz}\Big].~~ \tag {17} \end{align} $$ The new mode functions $U_n^{j+1} (z)$ can be easily determined by the inverse transformation of $\tilde {U}_n^{j+1} (\tilde {z})$. Thus we can determine the mode amplitudes at the $(j+1)$th step from those at the $j$th step without solving the depth separated wave equation. This will save considerable computation effort. For the initial conditions, the mode coefficients are[6] $$\begin{align} A_n^1 =\sqrt {\frac{2\pi }{k_n }} \frac{e^{i\pi /4}}{\rho (z_{\rm s})}U_n (z_{\rm s}),~~ \tag {18} \end{align} $$ where $z_{\rm s}$ is the source depth. All calculations were conducted using MATLAB code, and the mode functions in specific environments were calculated by KRAKENC.[12] To test the accuracy of this method, the coupled modes' results in the single scatter approximation are set as the reference solutions. The configuration used with up-slope and down-slope segments in the environment is shown in Fig. 1. The water column is homogeneous with $c_{\rm w}=1500$ m/s, $\rho_{\rm w}=1$ g/cm$^{3}$, bounded by a pressure release sea surface and homogeneous bottom ($c_{\rm b} ={\rm 1600}$ m/s, $\rho_{\rm b}=1.5$ g/cm$^{3}$, $\beta =0.5$ dB/$\lambda$). The water depth remains constant at $h_0 =50$ m from the source position to the range $r_0=1$ km, then decreases linearly from $h_0$ to $h_1=26$ m in the range of $r_1=3$ km, and then increases linearly from $h_1$ to $h_0$ at $r_2=5$ km. The computational region extends to a 'false floor' depth $H$ that remains constant at 100 m, which is deep enough to avoid false floor reflections. In the calculations discussed here, a total of 60 normal modes were used, including 23 trapped modes at the source position where the water depth is 50 m. The distance between steps was set as 10 m, thus the step height $\Delta h=0.12$ m. It has been found that a suitable step height is approximately $\lambda /15$. If the step height is much greater than this, the result will quickly become divergent. If the step height is much less than this, the computation burden in unnecessarily increased.
cpl-34-7-074302-fig1.png
Fig. 1. The environment configuration with up-slope and down-slope segments.
Table 1. The normal mode computation time after step 1 with 60 total modes. The computation for the method of this work and of H&T is based on starting with exact modes just before step 1 and using perturbation theory to account for step 1. KRAKENC uses the exact environment for the modes after step 1.
Frequency (Hz) Normal mode computation time (s)
This work H&T KRAKENC
900 0.110 0.101 33.23
920 0.119 0.100 33.28
940 0.109 0.099 33.74
960 0.109 0.102 33.62
980 0.126 0.100 34.08
1000 0.109 0.103 33.81
1020 0.117 0.105 33.57
1040 0.113 0.103 32.73
1060 0.121 0.099 32.78
1080 0.114 0.108 33.00
1100 0.110 0.100 32.80
cpl-34-7-074302-fig2.png
Fig. 2. The errors for the mode 15 horizontal wave number in comparison with KRAKENC. Here 60 modes are used in total.
Table 1 lists the computation time for obtaining the local normal modes after the first step in the bottom, 60 modes included. The results show that the extended perturbation method cost slightly more computation time than the H&T method, but is still faster by more than two orders of magnitude when compared with the computation time of KRAKENC. Figure 2 shows the errors for the mode 15 horizontal wave number at different frequencies after the first step. The errors for the method of H&T increase approximately linearly as the frequency increases, and becomes divergent over a series of steps, while the extended perturbation method described here has a nearly constant error at about $1.25\times 10^{-5}$ m$^{-1}$ over the frequency range shown. Thus the extended perturbation method leads to a minor increase in computation time, but improves the normal modes prediction accuracy, and increases the range of frequencies over which the method is accurate.
cpl-34-7-074302-fig3.png
Fig. 3. Transmission loss comparison between the reference solution and the adiabatic approximation solution at depth 20 m, and the frequency is 1 kHz.
A 1 kHz source is positioned at $z_{\rm s}=30$ m. A transmission loss comparison between the reference solution and the solution obtained using the adiabatic approximation is shown in Fig. 3.
cpl-34-7-074302-fig4.png
Fig. 4. Transmission loss comparison between the reference solutions and the extended perturbation solutions at depth 20 m, and the frequency is 1 kHz.
The large errors in Fig. 3 show that strong mode coupling is caused by the sloping bottom, which cannot be ignored for obtaining accurate results. The reference solution results and extended perturbation results are compared in Fig. 4. In the calculation, the perturbed mode functions are replaced by mode functions from KRAKENC every 10 steps in up-slope propagation and 5 steps in down-slope propagation. The results in Fig. 4 show that the extended perturbation solution agrees quite well with the reference solutions over the whole range, and the errors at the peaks are smaller than 0.5 dB. Thus the extended perturbation solution gives a very good approximation to the reference solutions. In Fig. 4, both computations use 400 stair steps. In the reference solution, the eigenvalue equation is solved 400 times. By contrast, the perturbation approximation results use exact normal modes only 10% of the time for the up-slope propagation, and 20% of the time for down-slope propagation compared with the reference solution, which leads to considerable saving in computation efforts. In summary, the extended perturbation coupled mode theory works well at higher frequencies, such as 1 kHz, and the method is fast and accurate. The method can include any number of normal modes, and can be applied to both far-field and near-field calculations. This method can also be easily incorporated into other mode-based sound field models, such as the couple mode-parabolic equation (CMPE).[13,14] As the slope angle is increased, the method still works if the step length is decreased. One issue not discussed here is the fact that $\tilde {z}(H)$ varies as the water depth changes, which has not been taken into account. This may contribute to the small error that does occur and lead to the requirement of using new mode sets after a modest number of steps. This issue will be addressed in future work. Also, the ultimate limitation of the frequency range for accurate results and the consistency of mode functions' update interval between up-slope propagation and down-slope propagation are still in need of further investigation. References A coupled mode solution for acoustic propagation in a waveguide with stepwise depth variations of a penetrable bottomCoupled perturbed modes over a sloping penetrable bottomCoupled mode perturbation theory of range dependenceCoupled perturbed modes and internal solitary wavesA vertical eigenfunction expansion for the propagation of sound in a downward-refracting atmosphere over a complex impedance planeThe problem of energy conservation in one-way models
[1] Evans R B 1983 J. Acoust. Soc. Am. 74 188
[2] Higham C J and Tindle C T 2003 J. Acoust. Soc. Am. 114 3119
[3] Tindle C T, O' Driscoll L M and Higham C J 2000 J. Acoust. Soc. Am. 108 76
[4] Higham C J and Tindle C T 2003 J. Acoust. Soc. Am. 113 2515
[5]Wang H Z, Wang N and Gao W 2009 Periodical of Ocean University of China 39 160 (in Chinese)
[6]Jensen F B, Kuperman W A, Porter M B and Schmidt H 2011 Computational Ocean Acoustics 2nd edn (New York: Springer)
[7]Courant R and Hilbert D 1953 Methods of Mathematical Physics (New York: Interscience) vol I chap VI
[8]Wilcox C H 1984 Sound Propagation in Stratified Fluids (New York: Springer) chap 2 sec 1
Wilcox C H 1984 Sound Propagation in Stratified Fluids (New York: Springer) chap 3 sec 9
[9]Friedman B 1956 Principles and Techniques of Applied Mathematics (New York: Wiley) chap 4
[10] Waxler R 2002 J. Acoust. Soc. Am. 112 2540
[11] Porter M B, Jensen F B and Ferla C M 1991 J. Acoust. Soc. Am. 89 1058
[12]Porter M B The KRAKEN Normal Mode Program (DRAFT) (SACLANT Undersea Research Centre)
[13]Abawi A T, Kuperman W A and Collins M D 1997 J. Acoust. Soc. Am. 102 233
[14]Peng Z H and li F H 2001 Sci. Chin. A 31 165 (in Chinese)