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

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.

Tau spectral method is far faster than the Legendre-Galerkin spectral method and slightly slower than the classic finite difference method.

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].

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 where ω = 2πf is the angular frequency and f is the frequency of the sound source. By using the technique of the separation of variables, the timeindependent sound pressure p(r, z) can be decomposed into where R(r) satisfies the Hankel function, and ψ(z) satisfies the following modal equation: where k 2 r is a separation constant and k r is the horizontal wavenumber. The function ρ(z) is the density, k = (1 + iηα)ω/c(z) is the complex wavenumber, α is the attenuation in dB/λ (λ is the wavelength of the sound source), and η = (40π log 10 e) −1 . In this paper, a two-layer media problem consisting of a layer of water (z ∈ [0, h]) and a layer of sedimentation (z ∈ [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: The boundary conditions at the ends of the interval 0 ≤ z ≤ H and the interface conditions at z = h, where h ≤ H, are imposed. The upper boundary condition is The lower boundary condition is: Perfectly free or rigid boundary Perfectly free boundary Water Sediment Figure 1: Sound speed, density and attenuation profiles of the layered marine environment.
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 ψ at the interface: where the superscripts − and + indicate limits from above and below, respectively. Eq. (3) has a set of solutions (k r,m , ψ m ), m = 1, 2, . . . when supplemented by the above boundary conditions, where k r,m and ψ 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 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 where z s is the depth of the sound source, and H (1) 0 (·) 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: 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

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 Any smooth differentiable function u(x) defined in the interval [−1, 1] can be expanded according to this set of basis functions as follows: We use the hat symbol to represent the Chebyshev spectral expansion coefficient of the corresponding function. For example,û k and (uv) k represent the k-th spectral expansion coefficients of functions u and uv, respectively. The expansion coefficientû k iŝ 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]: 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 u k of Eq. (15) can be approximately calculated aŝ and the discretization version of Eq. (14) is thus Similar to expanding the function u(x) in Eq. (14), the first derivative function of u(x) can also be expanded approximately as The expansion coefficients satisfy the following relationship: where c k is defined in Eq. (15). Equivalently, we can also express Eq. (20) in the form of a matrix-vector product whereû ′ = [û ′ 0 , . . . ,û ′ N ] T ,û = [û 0 , . . . ,û N ] T , and 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]: its equivalent vector form is where (uv) = [ (uv) 0 , . . . , (uv) N ] T and C v is a (N + 1)-order square matrix, whose elements depend on the values ofv 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]:

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 Both ψ w (z) and ψ b (z) satisfy Eq. (3) Here, we use the fact that dx/dz = −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): 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 By following these conventions, we apply Eqs. (21) and (23) to Eq. (28) twice, and we can obtain As a result, the matrix form of Eq. (3) in the Chebyshev-Tau method can be finally written as Applying this result to the water column with |I| = h and the bottom sediment with |I| = H − h, we can obtain the following: Both matrices A and B are square matrices of order (N + 1), and we split matrix A, for example, into four blocks: where block matrices A 11 and A 22 are square matrices of order (N − 1) and 2, respectively. In a similar way, we can also obtain four block matrices of B to be B 11 , B 12 , B 21 and 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ψ w andψ 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 Thus, the boundary conditions and interface conditions of Eqs. (5) to (9) can be transformed and expressed as We divide the four row vectors into blocks, and we name the row vectors consisting of the first (N − 1) elements as s 1 , t 1 , p 1 , and q 1 and the row vectors consisting of the last two elements as s 2 , t 2 , p 2 , and q 2 , respectively. Similarly, let the first (N − 1) elements of the column vectorsψ w andψ b bê ψ w1 andψ 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 According to the block structure of the matrix and/or vectors in the above formula, Eq. (38) can be abbreviated as where L 11 is a square matrix of order (2N − 2), and L 22 is a square matrix of order 4. The (2N − 2)-dimensional column vectorψ 1 is obtained by stacking the (N − 1)-dimensional column vectorsψ w1 andψ b1 and the 4-dimensional The remaining matrices satisfy the relevant compatibility conditions. The resulting system of Eq. (39) can be solved by the following two steps: Therefore, it is necessary only to solve the (2N −2)-order algebraic eigenvalue problem in Eq. (40). For each set of eigenvalues / eigenvectors (k 2 r ,ψ), a set of eigensolutions (k r , ψ) 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ψ w andψ b to the depth at the same resolution. The vectors ψ w and ψ b should be stacked into a single column vector to form ψ, 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.

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 α 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 MAT-LAB 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. The pseudolinear profile [25] is shown in Fig. 2(a) and given by: 2. The surface duct profile [26] c 2 (z) is shown in Fig. 2(b).
3. The linear sound speed profile in the bottom layer is shown in Fig. 2(c) and given by 4. The exponentially increasing sound speed profile in the bottom layer is shown in Fig. 2(d) and given by 5. The segmented sound speed profile in the water column is shown in Fig. 2(e) and given by 6. The other linear increasing sound speed profile in the bottom sediment is given by c 6 (z) = 0.5z + 1500 The plot of c 6 (z) is similar to that of c 3 (z) and is not shown in Fig. 2. 7. The linear increasing attenuation profile in the bottom layer is shown in Fig. 2(f) and given by

Single-layer waveguide
Example 1 of Table 1 is a single-layer waveguide problem. The density of seawater is constant at 1 g/cm 3 , 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.

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    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.    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.  Example 4 adds bottom attenuation 0.5 dB/λ based on example 3. Table  0  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. 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.
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. 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.
From the analysis of the above examples, we can see that in general, NM-   Table 9: Comparison of k r,m of example 7. 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.

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.