# A Chebyshev-Tau spectral method for normal modes of underwater sound
propagation with a layered marine environment

Houwang Tu  Yongxian Wang<EMAIL_ADDRESS>Qiang Lan  Wei Liu  Wenbin Xiao
Shuqing Ma  College of Meteorology and Oceanography, National University of
Defense Technology, Changsha, China

###### Abstract

The normal mode model is one of the most popular approaches for solving
underwater sound propagation problems. Among other methods, the finite
difference method is widely used in classic normal mode programs. In many
recent studies, the spectral method has been used for discretization. It is
generally more accurate than the finite difference method. However, the
spectral method requires that the variables to be solved are continuous in
space, and the traditional spectral method is powerless for a layered marine
environment. A Chebyshev-Tau spectral method based on domain decomposition is
applied to the construction of underwater acoustic normal modes in this paper.
In this method, the differential equation is projected onto spectral space
from the original physical space with the help of an orthogonal basis of
Chebyshev polynomials. A complex matrix eigenvalue / eigenvector problem is
thus formed, from which the solution of horizontal wavenumbers and modal
functions can be solved. The validity of the acoustic field calculation is
tested in comparison with classic programs. The results of analysis and tests
show that compared with the classic finite difference method, the proposed
Chebyshev-Tau spectral method has the advantage of high computational
accuracy. In addition, in terms of running time, our method is faster than the
Legendre-Galerkin spectral method.

###### keywords:

Chebyshev-Tau method , spectral method , normal modes , underwater sound
propagation , computational ocean acoustics

††journal: Journal of Sound and Vibration

A Chebyshev-Tau spectral method for normal modes of underwater sound
propagation with a layered marine environment

Research highlight 1 : This paper proposes a Chebyshev-Tau spectral method for
calculating acoustic propagation in a layered marine environment.

Research highlight 2 : The accuracy of the proposed method is higher than that
of the classic finite difference method.

Research highlight 3 : Regarding the running time, our Chebyshev-Tau spectral
method is far faster than the Legendre-Galerkin spectral method and slightly
slower than the classic finite difference method.

##  1 Introduction

Sound waves are the main means of transmitting information remotely
underwater. The study of the laws of underwater acoustic propagation helps
investigators develop and utilize the ocean efficiently. The propagation of
underwater sound waves can be modeled by the basic wave equation. Solving the
wave equation directly is very difficult, but by approximating the wave
equation from different perspectives, many mature underwater acoustic
propagation models have gradually been developed. The depth separated wave
equation for normal modes is common in underwater acoustics. One of the
earliest papers was published in 1948 by Pekeris [1], who developed a theory
for a simple two-layer model (ocean and sediment) with a constant sound speed
in each layer. Progress in the development of normal mode methods is presented
in an excellent summary by Williams [2] published in 1970. Today, there are
many models available that are based on normal modes [3, 4]. In terms of
numerical discretization, the finite difference method is the most common
method of discretizing the normal mode model, although some studies have used
other methods, such as the finite element method, to discretize the normal
mode model [5, 6, 7]. As one of the widely used programs of the normal mode
model in the acoustics community, Kraken [8] was developed based on the finite
difference method and has the advantages of robustness, accuracy and high
efficiency.

In engineering numerical simulations, the spectral method is a numerical
discrete method with higher accuracy than the finite difference method [9].
Since the 1970s, with the discovery of the fast Fourier transform and the
development of high-performance computing technology, the shortcoming of the
large amount of calculations of the spectral methods is gradually being
overcome. Due to the advantages of high accuracy and fast convergence [10],
the spectral method has been widely used in the research studies of many
disciplines, such as meteorology [11], oceanography [10], and physical and
chemical measurements [12]. A number of researchers have tried to introduce
spectral methods in the study of acoustic problems. Adamou et al. [13]
proposed an efficient, accurate, and flexible numerical scheme based on
spectral methods to determine the dispersion curves and displacement/stress
profiles for modes in elastic guiding structures, possibly curved, layered,
damped, inhomogeneous, or anisotropic. Quintanilla et al. [14] attempted to
model guided elastic waves in generally anisotropic media using a spectral
collocation method. Colbrook et al. [15] presented a new approach for solving
acoustic scattering problems. Wise et al. [16] presented an arbitrary acoustic
source and sensor distributions in a Fourier collocation method. Wang et al.
[17] presented an improved expansion scheme for the acoustical wave
propagator. Evans [18] proposed a Legendre-Galerkin technique for differential
eigenvalue problems with complex and discontinuous coefficients in underwater
acoustics. This method requires that the lower boundary be the pressure
release boundary, and the computational speed is slow. Subsequently, Evans et
al. [19] studied the Legendre-Galerkin spectral method for constructing
atmospheric acoustic normal modes. Tu et al. [20] showed normal mode and
parabolic equation models using the Chebyshev spectral method to process a
single layer of a body of water with constant density and no attenuation.

In this article, we introduce a Chebyshev spectral method that can address the
discontinuity of the sound speed, density and attenuation profiles based on
the normal mode model. This method projects the original underwater acoustic
propagation problem onto a finite-dimensional subspace with Chebyshev
polynomials as the basis functions [21]. After solving the problem in the
spectral space, the expansion coefficients are inversely transformed back to
the real physical space. The Tau method is used in this study to perform the
spectral expansion of the boundary conditions and explicitly apply the
constraints of the boundary conditions to the solution of the spectral
coefficients. In the proposed method, it is necessary to consider the
discontinuity of the parameters at the interface of the water column and
bottom sediment; thus, domain decomposition is used in the depth direction to
address this issue [22].

##  2 Normal mode model in a layered marine environment

Considering a cylindrical coordinate system where the sound source is a simple
harmonic point source and the marine environment is a cylindrically symmetric
two-dimensional sound propagation case, the Helmholtz equation can be written
as

$\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial p}{\partial
r}\right)+\rho(z)\frac{\partial}{\partial
z}\left(\frac{1}{\rho(z)}\frac{\partial p}{\partial
z}\right)+\frac{\omega^{2}}{c(z)^{2}}p=0$ (1)

where $\omega=2\pi f$ is the angular frequency and $f$ is the frequency of the
sound source. By using the technique of the separation of variables, the time-
independent sound pressure $p(r,z)$ can be decomposed into

$p(r,z)=\psi(z)R(r)$ (2)

where $R(r)$ satisfies the Hankel function, and $\psi(z)$ satisfies the
following modal equation:

$\rho(z)\frac{\mathrm{d}}{\mathrm{d}z}\left(\frac{1}{\rho(z)}\frac{\mathrm{d}\psi(z)}{\mathrm{d}z}\right)+k^{2}\psi(z)=k_{r}^{2}\psi(z)$
(3)

where $k_{r}^{2}$ is a separation constant and $k_{r}$ is the horizontal
wavenumber. The function $\rho(z)$ is the density,
$k=(1+i\eta\alpha)\omega/c(z)$ is the complex wavenumber, $\alpha$ is the
attenuation in dB$/\lambda$ ($\lambda$ is the wavelength of the sound source),
and $\eta=(40\pi\log_{10}{e})^{-1}$. In this paper, a two-layer media problem
consisting of a layer of water ($z\in[0,h]$) and a layer of sedimentation
($z\in[h,H]$) is considered. As shown in Fig. 1, the marine environmental
parameters in Eq. (3) are defined in the water and bottom separately, and they
can be discontinuous at the interface $z=h$:

$\begin{split}c(z)&=\begin{cases}c_{w}(z),&0\leq z\leq h\\\ c_{b}(z),&h\leq
z\leq H\end{cases},\cr\rho(z)&=\begin{cases}\rho_{w}(z),&0\leq z\leq h\\\
\rho_{b}(z),&h\leq z\leq
H\end{cases},\cr\alpha(z)&=\begin{cases}\alpha_{w}(z),&0\leq z\leq h\\\
\alpha_{b}(z),&h\leq z\leq H\end{cases}\end{split}$ (4)

Figure 1: Sound speed, density and attenuation profiles of the layered marine
environment.

The boundary conditions at the ends of the interval $0\leq z\leq H$ and the
interface conditions at $z=h$, where $h\leq H$, are imposed. The upper
boundary condition is

$\psi(z=0)=0$ (5)

The lower boundary condition is:

$\psi(z=H)=0$ (6)

or

$\frac{\mathrm{d}\psi(z=H)}{\mathrm{d}z}=0$ (7)

which correspond to a perfectly free surface at $z=0$ and a perfectly free or
rigid bottom at $z=H$. At the interface, the sound pressure and the normal
particle velocity of sound pressure must be continuous; these conditions are
imposed by allowing for a unique value $\psi$ at the interface:

$\psi(z=h^{-})=\psi(z=h^{+})$ (8)
$\frac{1}{\rho(z=h^{-})}\frac{\mathrm{d}\psi(z=h^{-})}{\mathrm{d}z}=\frac{1}{\rho(z=h^{+})}\frac{\mathrm{d}\psi(z=h^{+})}{\mathrm{d}z}$
(9)

where the superscripts $-$ and $+$ indicate limits from above and below,
respectively. Eq. (3) has a set of solutions $(k_{r,m},\psi_{m}),m=1,2,\dots$
when supplemented by the above boundary conditions, where $k_{r,m}$ and
$\psi_{m}$ are also called the $m$-th horizontal wavenumber and eigenmode,
respectively. The eigenmodes of Eq. (3) are arbitrary to a non-zero scaling
constant. We assume that the modes are scaled (normalized) so that

$\int_{0}^{H}\frac{{[\psi_{m}(z)}]^{2}}{\rho(z)}\mathrm{d}z=1,\quad
m=1,2,\dots,M$ (10)

where $m$ is the order of modes and $M$ is the number of modes. Finally, the
fundamental solution of the 2D Helmholtz equation can be approximated by

$p(r,z)\approx\frac{i}{4\rho(z_{s})}\sum_{m=1}^{M}\psi_{m}(z_{s})\psi_{m}(z)H_{0}^{(1)}(k_{r,m}r)$
(11)

where $z_{s}$ is the depth of the sound source, and $H_{0}^{(1)}(\cdot)$ is
the Hankel function, which corresponds to $R(r)$.

To demonstrate the sound field results, the transmission loss (TL) of the
sound pressure is defined as follows:

$\text{TL}=-20\log_{10}\left(\frac{|p|}{|p_{0}|}\right)$ (12)

The unit of TL is the decibel (dB), where $p_{0}$ is the sound pressure at a
distance of 1 m from the sound source. In an actual sound field display, TL is
usually used.

##  3 Chebyshev-Tau spectral method for the normal mode model

###  3.1 Chebyshev spectral method

The spectral method is based on using finite-order function expansion and
truncation to approximate the original function. The continuous smooth
function $u(x)$ is expanded according to a set of basis functions in a
weighted sum form. This set of basis functions, also known as trial functions,
constitutes a set of orthogonal bases in a continuous function space. If the
basis functions are Chebyshev orthogonal polynomials, the Chebyshev spectral
method is formed. In this case, the trial function $T_{k}(x)$ is

$T_{k}(x)=\cos{(k\cos^{-1}(x))},\quad k=0,1,2,\cdots$ (13)

Any smooth differentiable function $u(x)$ defined in the interval $[-1,1]$ can
be expanded according to this set of basis functions as follows:

$u(x)=\sum_{k=0}^{\infty}\hat{u}_{k}T_{k}(x)$ (14)

We use the hat symbol to represent the Chebyshev spectral expansion
coefficient of the corresponding function. For example, $\hat{u}_{k}$ and
$\widehat{(uv)}_{k}$ represent the $k$-th spectral expansion coefficients of
functions $u$ and $uv$, respectively. The expansion coefficient $\hat{u}_{k}$
is

$\hat{u}_{k}=\frac{2}{\pi
c_{k}}\int_{-1}^{1}\frac{u(x)T_{k}(x)}{\sqrt{1-x^{2}}}\mathrm{d}x,\quad
c_{k}=\begin{cases}2,&k=0\\\ 1,&k>0\end{cases}$ (15)

Eqs. (14) and (15) are called the inverse Chebyshev transform and forward
Chebyshev transform, respectively.

When numerically calculating the integral in Eq. (15), it is necessary to
introduce numerical discretization for the variable $x$ in the interval
$[-1,1]$. Considering the convenience of the processing boundary conditions
and the accuracy of the numerical calculations, the Gauss-Lobatto nodes are
often used [10]:

$x_{j}=\cos{\left(\frac{j\pi}{N}\right)},\quad j=0,1,2,\dots,N$ (16)

where $N$ is the number of discretization points, which is equal to the
truncation order in the spectral method. In this case, the expansion
coefficient $\hat{u}_{k}$ of Eq. (15) can be approximately calculated as

$\begin{split}\hat{u}_{k}\approx\frac{1}{d_{k}}\sum_{j=0}^{N}u(x_{j})T_{k}(x_{j})\omega_{j},\quad
k=0,1,2,\dots,N\\\ \omega_{j}=\begin{cases}\frac{\pi}{2N},\quad j=0,N\\\
\frac{\pi}{N},\quad\text{otherwise}\end{cases},d_{k}=\begin{cases}\pi,\quad
k=0,N\\\ \frac{\pi}{2},\quad\text{otherwise}\end{cases}\end{split}$ (17)

and the discretization version of Eq. (14) is thus

$u(x)\approx\sum_{k=0}^{N}\hat{u}_{k}T_{k}(x)$ (18)

Similar to expanding the function $u(x)$ in Eq. (14), the first derivative
function of $u(x)$ can also be expanded approximately as

$u^{\prime}(x)\approx\sum_{k=0}^{N}\hat{u}^{\prime}_{k}T_{k}(x)$ (19)

The expansion coefficients satisfy the following relationship:

$\hat{u}^{\prime}_{k}\approx\frac{2}{c_{k}}\sum_{\begin{subarray}{c}j=k+1,\\\
j+k=\text{odd}\end{subarray}}^{N}j\hat{u}_{j}$ (20)

where $c_{k}$ is defined in Eq. (15). Equivalently, we can also express Eq.
(20) in the form of a matrix-vector product

$\hat{\mathbf{u}}^{\prime}\approx\mathbf{D}\hat{\mathbf{u}}$ (21)

where
$\hat{\mathbf{u}}^{\prime}=[\hat{u}^{\prime}_{0},\dots,\hat{u}^{\prime}_{N}]^{T}$,
$\hat{\mathbf{u}}=[\hat{u}_{0},\dots,\hat{u}_{N}]^{T}$, and $\mathbf{D}$ is a
square matrix of order $(N+1)$. The superscript $T$ denotes the transpose of a
vector.

The spectral coefficients of the product of two known functions can also be
represented by the respective spectral coefficients [10]:

$\widehat{(uv)}_{k}\approx\frac{1}{2}\sum_{m+n=k}^{N}\hat{u}_{m}\hat{v}_{n}+\frac{1}{2}\sum_{|m-n|=k}^{N}\hat{u}_{m}\hat{v}_{n}$
(22)

its equivalent vector form is

$\widehat{\mathbf{(uv)}}\approx\mathbf{C}_{v}\hat{\mathbf{u}}$ (23)

where
$\widehat{\mathbf{(uv)}}=[\widehat{(uv)}_{0},\dots,\widehat{(uv)}_{N}]^{T}$
and $\mathbf{C}_{v}$ is a $(N+1)$-order square matrix, whose elements depend
on the values of $\hat{v}_{n}$.

To normalize each mode by using Eq. (10), we must evaluate the integration of
a known function $u(x)$ in interval $[-1,1]$. Here, we give only the following
formula briefly [23]:

$\int_{-1}^{1}u(x)\text{d}x\approx-2\sum_{\begin{subarray}{c}n=0,\\\
n=\text{even}\end{subarray}}^{N}\frac{\hat{u}_{n}}{n^{2}-1}$ (24)

###  3.2 Discrete normal mode model using the Chebyshev-Tau spectral method

The construction of the Chebyshev polynomial spanning starts with a choice of
the maximum order $N$ of the polynomials, which depends on the number of modes
needed, as discussed below. A single set of Chebyshev polynomials cannot span
both the water column and bottom sediment since it will not have the required
derivative discontinuity at the water-sediment interface. It is advisable to
use the domain decomposition method of Min and Gottleib [22]. The domain
decomposition takes the form

$\psi(z)=\begin{cases}\psi_{w}(z)\approx\sum_{k=0}^{N}\hat{\psi}_{w_{k}}T_{k}(x_{w}),\quad
x_{w}=-\frac{2}{h}z+1,&0\leq z\leq h\\\
\psi_{b}(z)\approx\sum_{k=0}^{N}\hat{\psi}_{b_{k}}T_{k}(x_{b}),\quad
x_{b}=\frac{2}{h-H}z+\frac{H+h}{H-h},&h\leq z\leq H\end{cases}$ (25)

Both $\psi_{w}(z)$ and $\psi_{b}(z)$ satisfy Eq. (3) on $[0,h]$ and $[h,H]$,
respectively.

To obtain the discretization equation of Eq. (3), as the first step, we
introduce the coordinate transformation of $z\in I\mapsto x\in[-1,1]$. Here,
the interval is $I=[0,h]$ for the water column layer and $I=[h,H]$ for the
bottom sediment layer. By applying the coordinate transformation to Eq. (3),
we obtain

$\frac{4}{|I|^{2}}\rho(x)\frac{\mathrm{d}}{\mathrm{d}x}\left(\frac{1}{\rho(x)}\frac{\mathrm{d}\psi(x)}{\mathrm{d}x}\right)+k^{2}\psi(x)=k_{r}^{2}\psi(x),\quad
x\in[-1,1]$ (26)

Here, we use the fact that ${\mathrm{d}x}/{\mathrm{d}z}=-2/|I|$, and $|I|$
denotes the length of the interval $I$. In the next step, we discretize the
problem into a system of linear equations. Constructing linear equations by
the spectral method requires the residual error to be orthogonal to a set of
finite dimensional bases. Classically, the Tau, Galerkin, collocation and
pseudospectral methods are commonly used for this purpose. We discuss the Tau
method in detail in this paper, which is essentially a method of mean weighted
residuals. When applied to solve Eq. (3), the Chebyshev-Tau method aims to
obtain the solution of the weak form of Eq. (3) (see Eq. (3.3.15) of [10] for
further details):

$\begin{split}\int_{-1}^{1}\left[\frac{4}{|I|^{2}}\rho(x)\frac{\mathrm{d}}{\mathrm{d}x}\left(\frac{1}{\rho(x)}\frac{\mathrm{d}\psi(x)}{\mathrm{d}x}\right)+k^{2}\psi(x)-k_{r}^{2}\psi(x)\right]\frac{T_{k}(x)}{\sqrt{1-x^{2}}}\mathrm{d}x=0\\\
x\in(-1,1),\quad k=0,1,\dots,N-2\end{split}$ (27)

The orthogonality of Chebyshev polynomials $\\{T_{k}(x)\\}$ can lead to a
system of $(N-1)$ linear equations, and the two endpoints of the interval $I$
are constrained by boundary conditions. For brevity in the following text, we
introduce some temporary functions $v(x),w(x),s(x)$ on $[-1,1]$ defined as

$v(x)=\frac{1}{\rho(x)},\quad w(x)=v(x)\psi^{\prime}(x),\quad
s(x)=\rho(x)w^{\prime}(x)$ (28)

By following these conventions, we apply Eqs. (21) and (23) to Eq. (28) twice,
and we can obtain

$\hat{\mathbf{w}}=-\frac{2}{|I|}\mathbf{C}_{v}(\mathbf{D}\hat{\bm{\psi})},\quad\hat{\mathbf{s}}=-\frac{2}{|I|}\mathbf{C}_{\rho}(\mathbf{D}\hat{\mathbf{w}})=\frac{4}{|I|^{2}}\mathbf{C}_{\rho}(\mathbf{D}(\mathbf{C}_{v}(\mathbf{D}\hat{\bm{\psi}})))$
(29)

As a result, the matrix form of Eq. (3) in the Chebyshev-Tau method can be
finally written as

$\left(\frac{4}{|I|^{2}}\mathbf{C}_{\rho}\mathbf{D}\mathbf{C}_{v}\mathbf{D}+\mathbf{C}_{k^{2}}\right)\hat{\bm{\psi}}=k_{r}^{2}\hat{\bm{\psi}}$
(30)

Applying this result to the water column with $|I|=h$ and the bottom sediment
with $|I|=H-h$, we can obtain the following:

$\displaystyle\mathbf{A}\hat{\bm{\psi}}_{w}$
$\displaystyle=k_{r}^{2}\hat{\bm{\psi}}_{w},$ $\displaystyle\mathbf{A}$
$\displaystyle=\frac{4}{h^{2}}\mathbf{C}_{\rho}\mathbf{D}\mathbf{C}_{v}\mathbf{D}+\mathbf{C}_{k^{2}},$
(31) $\displaystyle\mathbf{B}\hat{\bm{\psi}}_{b}$
$\displaystyle=k_{r}^{2}\hat{\bm{\psi}}_{b},$ $\displaystyle\mathbf{B}$
$\displaystyle=\frac{4}{(H-h)^{2}}\mathbf{C}_{\rho}\mathbf{D}\mathbf{C}_{v}\mathbf{D}+\mathbf{C}_{k^{2}}$

Both matrices $\mathbf{A}$ and $\mathbf{B}$ are square matrices of order
$(N+1)$, and we split matrix $\mathbf{A}$, for example, into four blocks:

$\mathbf{A}=\left[\begin{array}[]{cc}\mathbf{A}_{11}&\mathbf{A}_{12}\\\
\mathbf{A}_{21}&\mathbf{A}_{22}\\\ \end{array}\right]$ (32)

where block matrices $\mathbf{A}_{11}$ and $\mathbf{A}_{22}$ are square
matrices of order $(N-1)$ and 2, respectively. In a similar way, we can also
obtain four block matrices of $\mathbf{B}$ to be $\mathbf{B}_{11}$,
$\mathbf{B}_{12}$, $\mathbf{B}_{21}$ and $\mathbf{B}_{22}$.

As the weak form of Eq.(3) shows, the Tau method does not require each trial
function to satisfy the boundary conditions, which are directly enforced by
the equations of the system. Instead, the boundary conditions and interface
conditions Eqs. (5) to (9) are expanded with Eq. (17) and transformed into
linear equations about $\hat{\bm{\psi}}_{w}$ and $\hat{\bm{\psi}}_{b}$. In the
question we studied, the boundary conditions and the interface conditions
yield four linear equations that determine four of the coefficients of the two
highest-order polynomials. For convenience of description, we introduce the
intermediate row vectors

$\mathbf{s}=[s_{0},s_{1},s_{2},\dots,s_{N}],\quad\mathbf{t}=[t_{0},t_{1},t_{2},\dots,t_{N}]$

where $s_{k}=T_{k}(-1)=(-1)^{k}$, $t_{k}=T_{k}(+1)=1$. Furthermore, let

$\mathbf{p}=\frac{1}{\rho_{w}(h)}\mathbf{s}\mathbf{D},\quad\mathbf{q}=\frac{1}{\rho_{b}(h)}\mathbf{t}\mathbf{D}.$

Thus, the boundary conditions and interface conditions of Eqs. (5) to (9) can
be transformed and expressed as

$\displaystyle\mathbf{t}\hat{\bm{\psi}}_{w}=0$ (33)
$\displaystyle\mathbf{s}\hat{\bm{\psi}}_{b}=0$ (34)
$\displaystyle\mathbf{s}\mathbf{D}\hat{\bm{\psi}}_{b}=0$ (35)
$\displaystyle\mathbf{s}\hat{\bm{\psi}}_{w}-\mathbf{t}\hat{\bm{\psi}}_{b}=0$
(36)
$\displaystyle\mathbf{p}\hat{\bm{\psi}}_{w}-\mathbf{q}\hat{\bm{\psi}}_{b}=0$
(37)

We divide the four row vectors into blocks, and we name the row vectors
consisting of the first $(N-1)$ elements as $\mathbf{s}_{1}$,
$\mathbf{t}_{1}$, $\mathbf{p}_{1}$, and $\mathbf{q}_{1}$ and the row vectors
consisting of the last two elements as $\mathbf{s}_{2}$, $\mathbf{t}_{2}$,
$\mathbf{p}_{2}$, and $\mathbf{q}_{2}$, respectively. Similarly, let the first
$(N-1)$ elements of the column vectors $\hat{\bm{\psi}}_{w}$ and
$\hat{\bm{\psi}}_{b}$ be $\hat{\bm{\psi}}_{w1}$ and $\hat{\bm{\psi}}_{b1}$,
respectively. Since the interface conditions must consider both the water and
the bottom sediment, Eq. (34) is taken as the lower boundary condition as an
example. The matrix form of Eq. (31) becomes

$\left[\begin{array}[]{cc|cc}\mathbf{A}_{11}&\mathbf{0}&\mathbf{A}_{12}&\mathbf{0}\\\
\mathbf{0}&\mathbf{B}_{11}&\mathbf{0}&\mathbf{B}_{12}\\\
\hline\cr\mathbf{t}_{1}&\mathbf{0}&\mathbf{t}_{2}&\mathbf{0}\\\
\mathbf{s}_{1}&-\mathbf{t}_{1}&\mathbf{s}_{2}&-\mathbf{t}_{2}\\\
\mathbf{p}_{1}&-\mathbf{q}_{1}&\mathbf{p}_{2}&-\mathbf{q}_{2}\\\
\mathbf{0}&\mathbf{s}_{1}&\mathbf{0}&\mathbf{s}_{2}\end{array}\right]\left[\begin{array}[]{c}\hat{\bm{\psi}}_{w1}\\\
\hat{\bm{\psi}}_{b1}\\\ \hline\cr\hat{\bm{\psi}}_{w,N-1}\\\
\hat{\bm{\psi}}_{w,N}\\\ \hat{\bm{\psi}}_{b,N-1}\\\ \hat{\bm{\psi}}_{b,N}\\\
\end{array}\right]=k_{r}^{2}\left[\begin{array}[]{c}\hat{\bm{\psi}}_{w1}\\\
\hat{\bm{\psi}}_{b1}\\\ \hline\cr 0\\\ 0\\\ 0\\\ 0\\\ \end{array}\right]$ (38)

According to the block structure of the matrix and/or vectors in the above
formula, Eq. (38) can be abbreviated as

$\left[\begin{array}[]{cc}\mathbf{L}_{11}&\mathbf{L}_{12}\\\
\mathbf{L}_{21}&\mathbf{L}_{22}\\\
\end{array}\right]\left[\begin{array}[]{c}\hat{\bm{\psi}}_{1}\\\
\hat{\bm{\psi}}_{2}\end{array}\right]=k_{r}^{2}\left[\begin{array}[]{c}\hat{\bm{\psi}}_{1}\\\
\mathbf{0}\end{array}\right]$ (39)

where $\mathbf{L}_{11}$ is a square matrix of order $(2N-2)$, and
$\mathbf{L}_{22}$ is a square matrix of order 4. The $(2N-2)$-dimensional
column vector $\hat{\bm{\psi}}_{1}$ is obtained by stacking the
$(N-1)$-dimensional column vectors $\hat{\bm{\psi}}_{w1}$ and
$\hat{\bm{\psi}}_{b1}$ and the 4-dimensional column vector
$\hat{\bm{\psi}}_{2}=[\hat{\bm{\psi}}_{w,N-1},\hat{\bm{\psi}}_{w,N},\hat{\bm{\psi}}_{b,N-1},\hat{\bm{\psi}}_{b,N}]^{T}$.
The remaining matrices satisfy the relevant compatibility conditions. The
resulting system of Eq. (39) can be solved by the following two steps:

$\hat{\bm{\psi}}_{2}=-\mathbf{L}_{22}^{-1}\mathbf{L}_{21}\hat{\bm{\psi}}_{1},\quad(\mathbf{L}_{11}-\mathbf{L}_{12}\mathbf{L}_{22}^{-1}\mathbf{L}_{21})\hat{\bm{\psi}}_{1}=k_{r}^{2}\hat{\bm{\psi}}_{1}$
(40)

Therefore, it is necessary only to solve the $(2N-2)$-order algebraic
eigenvalue problem in Eq. (40). For each set of eigenvalues / eigenvectors
$(k_{r}^{2},\hat{\bm{\psi}})$, a set of eigensolutions $(k_{r},\bm{\psi})$ of
Eq. (3) can be obtained by the discrete inverse Chebyshev transform Eq. (18).
In this process, it is necessary to inversely transform the eigenvectors
$\hat{\bm{\psi}}_{w}$ and $\hat{\bm{\psi}}_{b}$ to the depth at the same
resolution. The vectors $\bm{\psi}_{w}$ and $\bm{\psi}_{b}$ should be stacked
into a single column vector to form $\bm{\psi}$, and each mode should be
normalized by Eq. (10). Finally, the sound pressure field is obtained by
applying Eq. (11) to every mode chosen by the favorable phase speed interval.
A systematic theoretical study of the convergence of the Chebyshev polynomials
and Tau method can be found in Boyd [21], and the asymptotic rate of
convergence in the Chebyshev-Tau spectral method has been previously reported
[24]; no more details are provided here.

##  4 Test examples and analysis of results

To verify the validity of the Chebyshev-Tau spectral method in solving the
normal mode model, the following tests and analysis are performed through
seven examples listed in Table 1. To make these examples closer to the real
ocean, the attenuation of seawater in the marine environment is set to
$\alpha_{w}=0.0$. For comparison, we introduce the widely used Kraken program
based on finite difference discretization and the rimLG program [18] based on
the Legendre-Galerkin spectral method. To facilitate the description, the
normal mode model program based on the Chebyshev-Tau spectral method developed
by us is abbreviated as the NM-CT program. Since the Kraken program is based
on Fortran language and the rimLG program is based on MATLAB language, we
developed both a Fortran code version and a MATLAB code version of NM-CT. The
numerical error of the results of the two versions is less than $10^{-14}$.
The results listed in the following text retain only ten decimal places; thus,
the results reported by the two implementations should be the same. It must be
mentioned that the size of $N$ in NM-CT depends on the number of modes
required for the final synthesized sound field, and the number of modes is
usually related to factors such as depth and sound source frequency.
Therefore, the choice of $N$ is empirical. We estimate the number of modes of
the final synthesized sound field through the depth, sound source frequency,
and modal phase velocity limits, similar to what Kraken did; the size of $N$
can usually be 2 to 10 times the number of modes. For a large number of
measured profile data, the size of $N$ is recommended to be 6 to 20 times the
number of modes.

The various profiles used for the numerical examples are

  1. 1.

The pseudolinear profile [25] is shown in Fig. 2(a) and given by:

$c_{1}(z)=\frac{1}{\sqrt{az+b}},a=5.94\times
10^{-10}\text{s}^{2}/\text{m}^{3},b=4.16\times
10^{-7}\text{s}^{2}/\text{m}^{3}$ (41)

  2. 2.

The surface duct profile [26] $c_{2}(z)$ is shown in Fig. 2(b).

  3. 3.

The linear sound speed profile in the bottom layer is shown in Fig. 2(c) and
given by

$c_{3}(z)=0.4z+1950$ (42)

  4. 4.

The exponentially increasing sound speed profile in the bottom layer is shown
in Fig. 2(d) and given by

$c_{4}(z)=2000-50\exp{\left(-\frac{z-100}{10}\right)}$ (43)

  5. 5.

The segmented sound speed profile in the water column is shown in Fig. 2(e)
and given by

$c_{5}(z)=\begin{cases}1520,\quad 0\leq z\leq 20\\\ 1600-4z,\quad 20\leq z\leq
30\\\ 1480,\quad 30\leq z\leq 200\end{cases}$ (44)

  6. 6.

The other linear increasing sound speed profile in the bottom sediment is
given by

$c_{6}(z)=0.5z+1500$ (45)

The plot of $c_{6}(z)$ is similar to that of $c_{3}(z)$ and is not shown in
Fig. 2.

  7. 7.

The linear increasing attenuation profile in the bottom layer is shown in Fig.
2(f) and given by

$\alpha_{1}(z)=0.0001875z+0.0125$ (46)

Figure 2: Acoustic profiles in the examples. Number of | ${\rho_{w}}$ | ${c_{w}}$ | ${\alpha_{w}}$ | ${h}$ | ${\rho_{b}}$ | ${c_{b}}$ | ${\alpha_{b}}$ | ${H}$ | Bottom
---|---|---|---|---|---|---|---|---|---
Examples |  g/cm3 | m/s |  dB/$\lambda$ | m |  g/cm3 | m/s |  dB/$\lambda$ | m | free/rigid
Example 1 | 1.0 | $c_{1}$ | 0.0 | 50 | 1.0 | $c_{1}$ | 0.0 | 100 |  Eq. (7)
Example 2 | 1.0 | 1500 | 0.0 | 50 | 1.5 | 1800 | 1.5 | 100 |  Eq. (6)
Example 3 | 1.0 | $c_{2}$ | 0.0 | 100 | 2.0 | $c_{3}$ | 0.0 | 150 |  Eq. (6)
Example 4 | 1.0 | $c_{2}$ | 0.0 | 100 | 2.0 | $c_{3}$ | 0.5 | 150 |  Eq. (6)
Example 5 | 1.0 | $c_{1}$ | 0.0 | 100 | 2.0 | $c_{4}$ | 0.5 | 200 |  Eq. (7)
Example 6 | 1.0 | $c_{1}$ | 0.0 | 100 | 2.0 | $c_{4}$ | 0.5 | 200 |  Eq. (6)
Example 7 | 1.0 | $c_{5}$ | 0.0 | 200 | 1.5 | $c_{6}$ | $\alpha_{1}$ | 1000 |  Eq. (6)
Table 1:  List of examples.

###  4.1 Single-layer waveguide

Example 1 of Table 1 is a single-layer waveguide problem. The density of
seawater is constant at 1 g/cm3, and there is no attenuation in the water. For
this type of single-layer marine environment, the NM-CT program developed
based on the Chebyshev-Tau spectral method is competent and needs only to set
the acoustic parameters in the two layers as a continuous layer in the input
file. In example 1, the pseudolinear sound speed waveguide has an analytical
solution given by the Airy functions [25]. The reliability of the NM-CT
program is proven by comparison with analytical solutions. Three source
frequencies, $f$ = 100 Hz, 1000 Hz and 7500 Hz, are considered. Because rimLG
can handle only the perfectly free lower boundary, it does not appear in
example 1. Table 2 shows a few horizontal wavenumbers of example 1 calculated
by the analytical solution, Kraken and NM-CT program for three frequencies.
NM-CT uses a spectral truncation order of $N$=100 at 100 Hz, $N$=400 at 1000
Hz and $N$=1000 at 7500 Hz, and 400, 2000 and 6000 equidistant discrete points
are used in the three frequencies for the Kraken program. It can be seen in
the table that the wavenumbers calculated by the three programs are very
similar; the Kraken result matches the analytical solutions with seven
significant digits for the three frequencies. The NM-CT matches the analytical
solutions with eight, seven and seven significant digits for $f$ = 100 Hz,
1000 Hz and 7500 Hz, respectively. Fig. 3 shows the errors of the wavenumbers
calculated by the two programs compared with the analytical solution. In the
figure, it can be seen more intuitively that regardless of the frequency, the
error of the wavenumber calculated by Kraken and NM-CT is approximately an
order of magnitude in this example. However, Kraken uses far more discrete
points in the vertical direction than NM-CT, which means that NM-CT needs
fewer discrete points to achieve higher accuracy.

Freq. | $m$ | Analytic | Kraken | NM-CT
---|---|---|---|---
$f$=100 Hz | 1 | 0.428726 2353 | 0.428726 2349 | 0.428726 2349
| 2 | 0.418720 6763 | 0.418720 6759 | 0.418720 6760
| 3 | 0.411401 3864 | 0.411401 3859 | 0.411401 3859
| 6 | 0.381781 7434 | 0.381781 7430 | 0.381781 7430
| 7 | 0.366073 2255 | 0.366073 2252 | 0.366073 2250
| 8 | 0.346775 4845 | 0.346775 4845 | 0.346775 4840
| 11 | 0.258936 7165 | 0.258936 7224 | 0.258936 7158
| 12 | 0.212957 4128 | 0.212957 4258 | 0.212957 4119
| 13 | 0.147227 5253 | 0.147227 5574 | 0.147227 5241
$f$=1000 Hz | 1 | 4.32256 5484 | 4.32256 5479 | 4.32256 5480
| 5 | 4.26193 0803 | 4.26193 0798 | 4.26193 0798
| 9 | 4.22231 1533 | 4.22231 1528 | 4.22231 1528
| 61 | 3.73510 9820 | 3.73510 9815 | 3.73510 9815
| 65 | 3.66898 6745 | 3.66898 6740 | 3.66898 6740
| 69 | 3.59716 8254 | 3.59716 8249 | 3.59716 8249
| 121 | 1.80461 2920 | 1.80461 2910 | 1.80461 2909
| 125 | 1.51324 5610 | 1.51324 5598 | 1.51324 5598
| 129 | 1.13647 1493 | 1.13647 1476 | 1.13647 1476
$f$=7500 Hz | 1 | 32.47270 963 | 32.47270 960 | 32.47270 960
| 32 | 31.97126 325 | 31.97126 322 | 31.97126 322
| 63 | 31.66388 669 | 31.66388 666 | 31.66388 666
| 156 | 30.95233 922 | 30.95233 918 | 30.95233 918
| 218 | 30.55369 016 | 30.55369 012 | 30.55369 012
| 311 | 29.84310 998 | 29.84310 994 | 29.84310 994
| 404 | 28.75361 660 | 28.75361 656 | 28.75361 656
| 528 | 26.71692 999 | 26.71692 995 | 26.71692 995
| 621 | 24.67372 146 | 24.67372 141 | 24.67372 141
Table 2: Comparison of $k_{r,m}$ of example 1.

Figure 3: The errors of the wavenumbers calculated by the Kraken (small blue
circles) and NM-CT (red plus sign) programs compared with the analytical
solutions.

###  4.2 Two-layer waveguide

Examples 2 to 7 of Table 1 are layered waveguide problems. When the bottom
sediment has attenuation, the modal wavenumbers are complex numbers. KrakenC
is a version of Kraken that finds the wavenumbers in the complex plane, which
computes the complex eigenvalues exactly. Therefore, we also tested the
results of KrakenC under the same configuration as Kraken for all examples
with attenuation.

Example 2 is the simplest layered waveguide with attenuated sediment, in which
the sound speed, density and attenuation are uniform constants. Table 3 shows
the horizontal wavenumbers calculated by the Kraken, KrakenC, rimLG and NM-CT
programs at frequencies of 20 Hz and 50 Hz, respectively. Kraken and KrakenC
use 1000 discrete points in the vertical direction, and the orders of NM-CT
and rimLG are only 20 and 40 at 20 Hz and 50 Hz, respectively. It can be seen
in the table that the wavenumbers calculated by the rimLG and NM-CT programs
are completely consistent at the two source frequencies, and the wavenumbers
calculated by the four programs are basically similar. Fig. 4(a)-(c) shows the
full-field TLs calculated by Kraken, rimLG and NM-CT at 20 Hz. Fig. 4(d) shows
the TLs at a depth of 10 m calculated by the three programs. The TL field
calculated by NM-CT is almost identical to Kraken. The difference between the
TLs obtained from Kraken and NM-CT has been found to be indistinguishable at
this plotting accuracy. Table 4 lists the running times of the three programs
in this example. The time listed in the table is the average of 10 test times.
In the test, a single thread of a single process is run on MATLAB 2019a on a
HUAWEI MatebookX PRO laptop computer with an Intel i7-8550U processor, and the
Fortran compiler used in the test is GFortran 7.4.0. It can be seen in the
table that the running time of NM-CT is almost equal to Kraken but much
shorter than rimLG. Although they both belong to the class of spectral
methods, our NM-CT program is faster than the rimLG program. As one of the
possible reasons, rimLG consumes too much time in computing a complex
symmetric matrix and a positive definite symmetric matrix to form the
generalized eigenvalue problem (see Eq. (4) in [18]), while our NM-CT program
aims to solve a simpler eigenvalue problem shown by Eq. (38).

Freq. | $m$ | Kraken | KrakenC | rimLG | NM-CT
---|---|---|---|---|---
$f$=20 Hz | 1 | 0.07350 70395+ | 0.07349 26194+ | 0.07350 28581+ | 0.07350 28581+
|  | 0.37693 38382e-3i | 0.37466 68873e-3i | 0.37597 26294e-3i | 0.37597 26294e-3i
| 2 | 0.04033 76782+ | 0.04034 46000+ | 0.04040 98897+ | 0.04040 98898+
|  | 0.23700 44780e-2i | 0.23763 40719e-2i | 0.23757 23752e-2i | 0.23757 23752e-2i
$f$=50 Hz | 1 | 0.20330 26836+ | 0.20329 22496+ | 0.20329 61543+ | 0.20329 61543+
|  | 0.14684 64650e-3i | 0.14474 24293e-3i | 0.14552 80250e-3i | 0.14552 80250e-3i
| 2 | 0.18325 56952+ | 0.18318 27296+ | 0.18320 16596+ | 0.18320 16596+
|  | 0.72305 34643e-3i | 0.71311 18231e-3i | 0.71805 23083e-3i | 0.71805 23083e-3i
| 3 | 0.16334 43629+ | 0.16336 23823+ | 0.16348 65836+ | 0.16348 65836+
|  | 0.44671 05031e-2i | 0.44867 81267e-2i | 0.44892 27771e-2i | 0.44892 27771e-2i
| 4 | 0.14201 67602+ | 0.14188 87957+ | 0.14195 94443+ | 0.14195 94443+
|  | 0.25964 36347e-2i | 0.26003 91729e-2i | 0.26101 78399e-2i | 0.26101 78399e-2i
| 5 | 0.11339 56065+ | 0.11358 52045+ | 0.11371 57329+ | 0.11371 57329+
|  | 0.47551 81270e-2i | 0.47349 61331e-2i | 0.47261 24780e-2i | 0.47261 24780e-2i
Table 3: Comparison of $k_{r,m}$ of example 2.

Figure 4: Full-field TLs calculated by Kraken (a), rimLG (b) and NM-CT (c); the black dotted line indicates the interface; TL curves at a depth of 10 m calculated by the three programs in example 2 ($f$=20 Hz). Freq. | Kraken | NM-CT(Fortran) | rimLG | NM-CT(MATLAB)
---|---|---|---|---
$f$=20 Hz | 0.29 s | 0.22 s | 15.16 s | 0.12 s
$f$=50 Hz | 0.35 s | 0.24 s | 15.42 s | 0.16 s
Table 4: Comparison of the running time of example 2.

In example 3, a surface duct sound speed profile is used in the water column.
This sound speed profile is not smooth enough; thus, it requires more spectral
truncation order. In this example, there is no attenuation in the bottom
sediment, so the results of KrakenC are not shown. Table 5 shows the
horizontal wavenumbers calculated by the Kraken, rimLG and NM-CT programs at
frequencies of 50 Hz and 100 Hz. The Kraken program has 1500 discrete points
in the vertical direction, and the orders of NM-CT and rimLG are $N$=500 at 50
Hz and $N$=1000 at 100 Hz, respectively. NM-CT matches Kraken with five and
six significant digits for $f$ = 50 Hz and 100 Hz, respectively. Fig. 5 shows
the six modes calculated by NM-CT at 50 Hz; in Fig. 5, the modes obtained from
Kraken and rimLG are also shown. The difference between the mode shapes
obtained from the three programs is found to be indistinguishable for this
plotting accuracy. This finding is true for all problems without bottom
attenuation. This result shows that the modes calculated by NM-CT are in good
agreement with the Kraken results. Table 6 lists the running times of the
three programs in this example. The conclusion is the same as in example 1 and
example 2; NM-CT is much faster than rimLG and slightly slower than Kraken.

Freq. | $m$ | Kraken | rimLG | NM-CT
---|---|---|---|---
$f$=50Hz | 1 | 0.20456 66998 | 0.20456 67640 | 0.20456 67092
| 2 | 0.19825 45322 | 0.19825 46184 | 0.19825 45356
| 3 | 0.18790 98058 | 0.18790 99137 | 0.18790 98062
| 4 | 0.17209 24614 | 0.17209 26321 | 0.17209 24594
| 5 | 0.15078 56837 | 0.15078 71664 | 0.15078 56682
| 6 | 1418862 39679 | 0.14189 02497 | 0.14188 62352
$f$=100Hz | 1 | 0.41201 50812 | 0.41201 51634 | 0.41201 50943
| 2 | 0.40748 98115 | 0.40748 99355 | 0.40748 98294
| 3 | 0.40194 64244 | 0.40194 65268 | 0.40194 64373
| 6 | 0.37130 65141 | 0.37130 66596 | 0.37130 65199
| 7 | 0.35529 03522 | 0.35529 05021 | 0.35529 03512
| 8 | 0.33589 86072 | 0.33589 88120 | 0.33589 86089
| 11 | 0.29020 30496 | 0.29020 93565 | 0.29020 30549
| 12 | 0.28087 06224 | 0.28087 37829 | 0.28087 06077
| 13 | 0.25576 83615 | 0.25577 45190 | 0.25576 84904
Table 5: Comparison of $k_{r,m}$ of example 3. Figure 5: The shapes of the first six modes of example 3, were calculated by Kraken (blue dashed line), rimLG (green dashed line) and NM-CT (red solid line) at 50 Hz. Freq. | Kraken | NM-CT(Fortran) | rimLG | NM-CT(MATLAB)
---|---|---|---|---
$f$=50 Hz | 0.31 s | 8.37 s | 1348.85 s | 4.23 s
$f$=100 Hz | 0.36 s | 46.76 s | 63826.37 s | 51.42 s
Table 6: Comparison of the running time of example 3.

Example 4 adds bottom attenuation 0.5 dB/$\lambda$ based on example 3. Table 7
shows the horizontal wavenumbers calculated by the Kraken, KrakenC, rimLG and
NM-CT programs at frequencies of 50 Hz and 100 Hz. It can be seen in the table
that the wavenumbers calculated by the four programs are still similar, which
shows that NM-CT can obtain credible results regardless of whether the bottom
is attenuated.

Freq. | $m$ | Kraken | KrakenC | rimLG | NM-CT
---|---|---|---|---|---
$f$=50 Hz | 1 | 0.20456 65851+ | 0.20456 66180+ | 0.20456 65949+ | 0.20456 65399+
|  | 0.84774 86844e-5i | 0.83525 56810e-5i | 0.83563 28501e-5i | 0.83558 18884e-5i
| 2 | 0.19825 42502+ | 0.19825 39711+ | 0.19825 42260+ | 0.19825 41431+
|  | 0.20928 28166e-4i | 0.20592 49724e-4i | 0.20601 43841e-4i | 0.20599 96787e-4i
| 3 | 0.18790 92721+ | 0.18790 89290+ | 0.18790 91494+ | 0.18790 90428+
|  | 0.39680 46830e-4i | 0.39122 93174e-4i | 0.39142 56820e-4i | 0.39139 51987e-4i
| 4 | 0.17209 14240+ | 0.17208 97671+ | 0.17209 05414+ | 0.17209 03689+
|  | 0.77265 86281e-4i | 0.76662 97284e-4i | 0.76718 40570e-4i | 0.76707 12343e-4i
| 5 | 0.15077 98663+ | 0.15073 87951+ | 0.15074 45128+ | 0.15074 30467+
|  | 0.47798 91588e-3i | 0.47546 02209e-3i | 0.47667 11983e-3i | 0.47635 77140e-3i
| 6 | 0.14187 53523+ | 0.14191 73773+ | 0.14193 12614+ | 0.14192 72116+
|  | 0.11089 56529e-2i | 0.11106 71974e-2i | 0.11096 78864e-2i | 0.11099 31242e-2i
$f$=100 Hz | 1 | 0.41201 50087+ | 0.41201 49345+ | 0.41201 50486+ | 0.41201 49733+
|  | 0.53367 85874e-5i | 0.52718 84909e-5i | 0.52746 05158e-5i | 0.52744 25387e-5i
| 2 | 0.40748 96727+ | 0.40748 95194+ | 0.40748 97232+ | 0.40748 96199+
|  | 0.10237 50490e-4i | 0.10098 91480e-4i | 0.10104 23904e-4i | 0.10103 67151e-4i
| 3 | 0.40194 62152+ | 0.40194 59869+ | 0.40194 62242+ | 0.40194 61277+
|  | 0.15473 27718e-4i | 0.15240 03585e-4i | 0.15247 38530e-4i | 0.15246 90493e-4i
| 6 | 0.37130 58929+ | 0.37130 51598+ | 0.37130 57092+ | 0.37130 55634+
|  | 0.45967 88595e-4i | 0.45379 25969e-4i | 0.45402 98311e-4i | 0.45400 48023e-4i
| 7 | 0.35528 94757+ | 0.35528 81718+ | 0.35528 89150+ | 0.35528 87638+
|  | 0.64790 73341e-4i | 0.64149 75073e-4i | 0.64188 59391e-4i | 0.64184 34819e-4i
| 8 | 0.33589 72138+ | 0.33589 37153+ | 0.33589 48592+ | 0.33589 46561+
|  | 0.10299 38733e-3i | 0.10229 27681e-3i | 0.10238 88776e-3i | 0.10237 32380e-3i
| 11 | 0.29018 11398+ | 0.29003 57030+ | 0.29006 30258+ | 0.29005 65311+
|  | 0.22085 62615e-2i | 0.22281 09393e-2i | 0.22319 59465e-2i | 0.22309 85787e-2i
| 12 | 0.28085 83769+ | 0.28100 43200+ | 0.28101 60710+ | 0.28101 30935+
|  | 0.10257 56310e-2i | 0.10047 80875e-2i | 0.10013 84950e-2i | 0.10022 50558e-2i
| 13 | 0.25574 78256+ | 0.25562 82278+ | 0.25565 45944+ | 0.25564 84910+
|  | 0.21995 85067e-2i | 0.21927 77341e-2i | 0.21959 49109e-2i | 0.21951 82731e-2i
Table 7: Comparison of $k_{r,m}$ of example 4.

In example 5, a perfectly rigid boundary condition is used, and a pseudolinear
profile and exponentially increasing sound profile are used in the water
column and bottom sediment, respectively. Kraken uses 2000 discrete points in
the vertical direction, and the orders of NM-CT are $N$=50 and $N$=100 at 50
Hz and 100 Hz, respectively. To intuitively see the wavenumbers calculated by
the three programs, Fig. 6 shows the wavenumbers in the complex plane. It can
be seen in the figure that the wavenumbers calculated by NM-CT are closer to
KrakenC than Kraken; in general, the results of the three programs are highly
similar in this example. Fig. 7 shows the random four modes in the example of
50 Hz. From the shapes of the modes, we can see that there are subtle and
imperceptible differences in the last modes.

Figure 6: Horizontal wavenumbers $k_{r}$ calculated by Kraken and NM-CT at 50
Hz (a) and 100 Hz (b) of example 5. Figure 7: Mode shapes for example 5
calculated by Kraken (blue dashed line) and NM-CT (red solid line) at 50 Hz.

Example 6 is very similar to Example 5, except that the lower boundary
condition is changed to perfectly free. Table 8 shows the horizontal
wavenumbers calculated by the Kraken, rimLG and NM-CT programs at frequencies
of 50 Hz and 100 Hz. It can be seen in the table that the wavenumbers
calculated by the three programs are still similar, which shows that NM-CT can
obtain convincing results regardless of whether the bottom is perfectly rigid
or free.

Freq. | $m$ | Kraken | KrakenC | rimLG | NM-CT
---|---|---|---|---|---
$f$=50 Hz | 1 | 0.21001 77364+ | 0.21001 75130+ | 0.21001 76897+ | 0.21001 76606+
|  | 0.16131 56264e-4i | 0.16127 40702e-4i | 0.16135 42837e-4i | 0.16134 19472e-4i
| 2 | 0.20257 20089+ | 0.20257 17620+ | 0.20257 19786+ | 0.20257 19436+
|  | 0.19856 05398e-4i | 0.19851 49397e-4i | 0.19860 68575e-4i | 0.19859 33106e-4i
| 3 | 0.19226 25150+ | 0.19226 20600+ | 0.19226 24339+ | 0.19226 23760+
|  | 0.34538 89758e-4i | 0.34529 76202e-4i | 0.34546 06351e-4i | 0.34543 83693e-4i
| 4 | 0.17666 02717+ | 0.17665 90008+ | 0.17665 96941+ | 0.17665 95967+
|  | 0.65210 59254e-4i | 0.65171 19663e-4i | 0.65209 40806e-4i | 0.65205 17663e-4i
| 5 | 0.15567 70622+ | 0.15553 29067+ | 0.15553 92553+ | 0.15553 90159+
|  | 0.70465 84115e-3i | 0.69470 12846e-3i | 0.69754 71573e-3i | 0.69755 14104e-3i
| 6 | 0.15289 03336+ | 0.15302 28270+ | 0.15303 20732+ | 0.15303 20250+
|  | 0.97022 58068e-3i | 0.97893 20248e-3i | 0.97650 59872e-3i | 0.97648 88614e-3i
| 7 | 0.14377 55359+ | 0.14378 63608+ | 0.14379 95849+ | 0.14379 94008+
|  | 0.14193 21581e-2i | 0.14205 58800e-2i | 0.14204 63379e-2i | 0.14204 72405e-2i
| 8 | 0.12835 94066+ | 0.12832 98501+ | 0.12834 25492+ | 0.12834 22181+
|  | 0.13425 76929e-2i | 0.13453 64528e-2i | 0.13460 03304e-2i | 0.13460 02038e-2i
$f$=100 Hz | 1 | 0.42460 86614+ | 0.42460 84054+ | 0.42460 86003+ | 0.42460 85640+
|  | 0.17336 02035e-4i | 0.17329 20134e-4i | 0.17338 44497e-4i | 0.17336 78968e-4i
| 2 | 0.41660 99568+ | 0.41660 97405+ | 0.41660 99159+ | 0.41660 98832+
|  | 0.15615 31545e-4i | 0.15609 31461e-4i | 0.15617 30488e-4i | 0.15615 89033e-4i
| 3 | 0.40973 02081+ | 0.40972 99945+ | 0.40973 01752+ | 0.40973 01418+
|  | 0.16097 66605e-4i | 0.16091 58499e-4i | 0.16099 59230e-4i | 0.16098 18988e-4i
| 7 | 0.36408 99387+ | 0.36408 89665+ | 0.36408 95771+ | 0.36408 94714+
|  | 0.55276 64987e-4i | 0.55247 10681e-4i | 0.55279 06003e-4i | 0.55274 10256e-4i
| 8 | 0.34512 12700+ | 0.34511 90586+ | 0.34511 99730+ | 0.34511 98231+
|  | 0.83894 99567e-4i | 0.83799 71712e-4i | 0.83862 68939e-4i | 0.83854 02070e-4i
| 9 | 0.32253 76329+ | 0.32251 97554+ | 0.32252 17833+ | 0.32252 15118+
|  | 0.20056 89504e-3i | 0.19780 03179e-3i | 0.19818 31732e-3i | 0.19815 01265e-3i
| 13 | 0.29477 93267+ | 0.29474 78357+ | 0.29475 63509+ | 0.29475 60875+
|  | 0.10859 36417e-2i | 0.10237 39337e-2i | 0.10246 11119e-2i | 0.10245 10168e-2i
| 14 | 0.28749 10525+ | 0.28761 77672+ | 0.28764 27674+ | 0.28764 23007+
|  | 0.25902 70973e-2i | 0.26176 44484e-2i | 0.26155 55357e-2i | 0.26156 25849e-2i
| 15 | 0.27344 09188+ | 0.27341 42528+ | 0.27344 18466+ | 0.27344 11802+
|  | 0.29129 68544e-2i | 0.29222 85503e-2i | 0.29229 56121e-2i | 0.29229 64710e-2i
Table 8: Comparison of $k_{r,m}$ of example 6.

Example 7 is an example with the discontinuous second derivative of the sound
speed profile from the literature [18], which is used only to test the
robustness of NM-CT in calculating the sound speed profile with a
discontinuous derivative. This is especially important in cases when a large
number of measured sound speeds are used. In this example, $f$=250 Hz, the
number of discrete points used in Kraken and KrakenC are automatically chosen
by the code, and the orders of both rimLG and NM-CT are $N$=500. Table 9 shows
the horizontal wavenumbers, and Fig. 8 shows the TL fields calculated by the
four programs. From the table and figure, we can draw the conclusion that NM-
CT can still obtain satisfactory results, which means that although the
Chebyshev-Tau method proposed in this article theoretically requires the
acoustic profiles to be sufficiently smooth, for the measured data that are
not so smooth, NM-CT can also obtain valuable results for reference by using a
large $N$.

$m$ | Kraken | KrakenC | rimLG | NM-CT
---|---|---|---|---
1 | 1.06120 73981+ | 1.06120 73980+ | 1.06120 74892+ | 1.06120 73934+
| 0.31452 71656e-7i | 0.31343 84747e-7i | 0.31311 49600e-7i | 0.31309 30200e-7i
4 | 1.05908 53065+ | 1.05908 53030+ | 1.05908 68080+ | 1.05908 52520+
| 0.49312 57091e-6i | 0.49134 81574e-6i | 0.49082 92710e-6i | 0.49079 67580e-6i
8 | 1.05230 17898+ | 1.05230 17200+ | 1.05230 80953+ | 1.05230 15655+
| 0.18612 79366e-5i | 0.18540 63045e-5i | 0.18519 88927e-5i | 0.18518 55544e-5i
16 | 1.02726 88405+ | 1.02726 86270+ | 1.02727 57016+ | 1.02726 88998+
| 0.46107 53397e-5i | 0.45934 31907e-5i | 0.46093 56507e-5i | 0.45871 06521e-5i
24 | 0.99057 13001+ | 0.99057 08911+ | 0.99058 83179+ | 0.99057 09466+
| 0.18246 84646e-4i | 0.18275 41277e-4i | 0.18588 97603e-4i | 0.18221 72908e-4i
40 | 0.93895 77115+ | 0.93895 69464+ | 0.94402 33694+ | 0.93895 88486+
| 0.11022 01392e-2i | 0.11109 67196e-2i | 0.82054 60202e-3i | 0.10524 36440e-2i
48 | 0.92154 89408+ | 0.92153 90152+ | 0.92896 10252+ | 0.92154 23293+
| 0.12494 56899e-2i | 0.12630 46083e-2i | 0.10808 99895e-2i | 0.11858 48928e-2i
50 | 0.91755 92536+ | 0.91755 79958+ | 0.92582 79956+ | 0.91755 89811+
| 0.88280 41759e-3i | 0.82583 86769e-3i | 0.82816 76867e-3i | 0.78282 21102e-3i
Table 9: Comparison of $k_{r,m}$ of example 7.

Figure 8: Full-field TLs calculated by Kraken (a), Kraken (b), rimLG (c) and
NM-CT (d) of example 7; the black dotted line indicates the interface.

From the analysis of the above examples, we can see that in general, NM-CT
uses fewer discrete points to achieve or exceed Kraken’s accuracy, especially
when the acoustical profiles are smooth. The small number of discrete points
in the vertical direction means that the order of spectral truncation is
relatively low. As a result, the calculation amount of the spectral method is
small, and the computational speed is not much longer than Kraken. When the
acoustical profiles are not smooth enough, the spectral methods (NM-CT and
rimLG) require a higher truncation order, which causes an increase in the
computational time. For the case of large $N$, the discrete Chebyshev forward
transform Eq. (17) and inverse transform Eq. (18) can be accelerated in
$(2.5N\log_{2}{N}+4N)$ float operations by using the FFT (fast Fourier
transforms) shown in Appendix B of [10]. In general, compared with Kraken, the
NM-CT program has a slower computational speed but higher accuracy. Compared
with rimLG, although they are both spectral methods, the computational speed
of NM-CT is much faster than that of rimLG when the accuracy is equivalent. To
make the comparison more convenient, NM-CT uses the same truncation orders in
both the water column and bottom sediment. In fact, we have tried to use
different $N$ in the water column and bottom sediment; in actual operation,
the truncation orders in the two layers can be freely specified in the input
file of the NM-CT according to the complexity of the profiles of the two
layers. The results show that NM-CT can still obtain sufficiently accurate
horizontal wavenumbers and sound fields in this configuration.

##  5 Conclusion and future work

The Chebyshev-Tau spectral method provides a simple matrix procedure for
finding two-layer underwater acoustic normal modes when the sea surface,
bottom and interface are range independent. This method first interpolates the
acquired data of the sound speed, density and attenuation profiles to the
Gauss-Lobatto points. After a permutation on the rows and columns in the
coefficient matrix, a complex eigensystem for solving wavenumbers and modal
spectral coefficients is formed that can be solved by applying numerical
libraries and algorithms. The validity of this Chebyshev-Tau spectral method
(NM-CT program) is demonstrated in comparison with the finite difference
method (Kraken program and its variations) and Legendre-Galerkin spectral
method (rimLG program). From the perspective of the computational speed, NM-CT
is superior to rimLG. From the perspective of computational accuracy, NM-CT is
also sufficiently high compared with Kraken and rimLG. In fact, this method
can also address more complex situations, such as irregular changes in the
density and attenuation coefficient with depth, but real ocean data are too
difficult to obtain. Therefore, this article does not show particularly
complicated numerical examples. Notably, in the case where the sound speed,
density and attenuation profiles are not smooth enough, the Chebyshev-Tau
spectral method can use a higher truncation order to obtain convincing
results.

This article simply divides the actual ocean into two layers, water and
sediment. In the actual ocean, the sediment may be composed of silt or sand
with completely different densities and attenuations. In this case, the
sediment should be divided into more intermittent layers. It can be seen from
the derivation in Sect. 3 that our method has the potential to work in
multilayered environments. Regarding the issue that the computational time is
longer than Kraken, we will try to speed up the running time through parallel
computing technology in the next work.

## Acknowledgments

This work was supported by the National Key Research and Development Program
of China [grant number 2016YFC1401800]; the National Natural Science
Foundation of China [grant numbers 61972406, 51709267]; and the “Double First-
Class” Construction Guidance Project of National University of Defense
Technology [grant number 4345161111L].

## References

  * [1] C. L. Pekeris, Theory of propagation of explosive sound in shallow water, Geological Society of America Memoirs 27 (1948) 1–116.  doi:10.1130/mem27-2-p1.
  * [2] A. O. Williams, Normal mode methods in propagation of underwater sound, Wiley-Interscience, New York, 1970.
  * [3] D. C. Stickler, Normal‐mode program with both the discrete and branch line contributions, The Journal of the Acoustical Society of America 57 (1975) 856–861.  doi:10.1121/1.380525.
  * [4] M. Porter, E. L. Reiss, A numerical method for ocean‐acoustic normal modes, The Journal of the Acoustical Society of America 76 (1984) 244–252.  doi:10.1121/1.391101.
  * [5] M. J. Isakson, R. A. Yarbrough, Finite‐element modeling of long range, range‐dependent acoustic propagation in shallow water, The Journal of the Acoustical Society of America 122 (2007) 3074.  doi:10.1121/1.2942977.
  * [6] C. P. Vendhana, G. C. Diwan, S. K. Bhattacharyya, Finite-element modeling of depth and range dependent acoustic propagation in oceanic waveguides, The Journal of the Acoustical Society of America 127 (2010) 3319.  doi:10.1121/1.3392440.
  * [7] R. Kirby, W. Duan, Modelling sound propagation in the ocean: a normal mode approach using finite elements, Proceedings of ACOUSTICS (2018) 1–10.
  * [8] M. B. Porter, The Kraken Normal Mode Program, SACLANT Undersea Research Centre, 2001.
  * [9] D. Gottlieb, S. A. Orszag, Numerical Analysis of Spectral Methods, Society for Industrial and Applied Mathematics, Philadelphia, 1977.  doi:10.1137/1.9781611970425.
  * [10] C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Zang, Spectral Methods Fundamentals in Single Domains, Spring-Verlag, Berlin, 2006.  doi:10.1007/978-3-540-30726-6.
  * [11] Q. L. Zheng, Experiments on the 30-day long-range numerical weather prediction in a seven-level spectral model, Journal of Academy of Meteorological Science 10 (1989) 234–246.
  * [12] N. P. Muravskaya, Spectral methods in physical and chemical measurements, Journal of Physics: Conference Series 1420(C1) (2019) 012005.
  * [13] A. T. I. Adamou, R. V. Craster, Spectral methods for modelling guided waves in elastic media, The Journal of the Acoustical Society of America 116 (2004) 1524–1535.  doi:10.1121/1.1777871.
  * [14] F. H. Quintanilla, M. J. S. Lowe, R. V. Craster, Modeling guided elastic waves in generally anisotropic media using a spectral collocation method, The Journal of the Acoustical Society of America 137 (2015) 1180–1194.  doi:10.1121/1.4913777.
  * [15] M. J. Colbrook, L. J. Ayton, A spectral collocation method for acoustic scattering by multiple elastic plates, Journal of Sound and Vibration 461 (2019) 114904.  doi:10.1016/j.jsv.2019.114904.
  * [16] E. S. Wise, B. T. Cox, J. Jaros, B. E. Treeby, Representing arbitrary acoustic source and sensor distributions in fourier collocation methods, The Journal of the Acoustical Society of America 146 (2019) 278–288.  doi:10.1121/1.5116132.
  * [17] J. B. Wang, J. Pan, Acoustical wave propagator with modified Chebyshev expansion, Computer Physics Communications 174 (2006) 187–190.  doi:10.1016/j.cpc.2005.10.002.
  * [18] R. B. Evans, A Legendre-Galerkin technique for differential eigenvalue problems with complex and discontinuous coefficients, arising in underwater acoustics (December 2016).
URL https://oalib-acoustics.org/Modes/rimLG/(2020)

  * [19] R. B. Evans, X. Di, K. E. Gilbert, A legendre-galerkin spectral method for constructing atmospheric acoustic normal modes, The Journal of the Acoustical Society of America 143 (2018) 3595–3601.  doi:10.1121/1.5040481.
  * [20] H. W. Tu, Y. X. Wang, W. Liu, X. Ma, W. B. Xiao, Q. Lan, A chebyshev spectral method for normal mode and parabolic equation models in underwater acoustics, Mathematical Problems in Engineering 2020 (2020) 1–12.  doi:10.1155/2020/7461314.
  * [21] J. P. Boyd, Chebyshev and Fourier Spectral Methods, Second Edition, Dover, New York, 2001.
  * [22] M. S. Min, D. Gottlieb, Domain decomposition spectral approximations for an eigenvalue problem with a piecewise constant coefficient, Siam Journal on Numerical Analysis 43 (2005) 502–520.  doi:10.1137/s0036142903423836.
  * [23] J. C. Mason, D. C. Handscomb, Chebyshev Polynomials, CRC Press, New York, 2003.  doi:10.1201/9781420036114.
  * [24] J. P. Coleman, The Lanczos Tau method, Journal of Institute of Mathematics and Applications 17 (1976) 85–97.  doi:10.1093/imamat/17.1.85.
  * [25] F. B. Jensen, W. A. Kuperman, M. B. Porter, H. Schmidt, Computational Ocean Acoustics, Springer-Verlag, New York, 2011.  doi:10.1007/978-1-4419-8678-8.
  * [26] A. D. Chowdhury, C. P. Vendhan, S. K. Bhattacharyya, S. Mudaliar, A Rayleigh-Ritz model for the depth eigenproblem of heterogeneous pekeris waveguides, ACTA ACUSTICA UNITED WITH ACUSTICA 104 (2018) 597–610.  doi:10.3813/aaa.919200.