Coupled-Mode Parabolic Equations for the Modeling of Sound Propagation in a Shallow-Water Waveguide with Weak Elastic Bottom

Abstract

In this work, a mode parabolic equation method with interacting modes accounting for the weak elasticity at the bottom is developed. An important feature of the proposed method is that computations of elastic modes are avoided and that the solution is obtained in the form of expansion over acoustic modes. A numerical technique for solving resulting mode parabolic equations is developed, and the accuracy and efficiency of the resulting solution is validated by a direct comparison against source image solutions in the 3D wedge benchmark problem. Satisfactory agreement of the two solutions is achieved for sufficiently small values of shear wave speed that are typical for soft sediments of the sea bottom. The developed approach may be used for solving 3D problems of sound propagation with the elastic properties of bottom sediments taken into account.

1. Introduction

The modeling of sound propagation in three-dimensional (3D) shallow water waveguides is one of the most important areas of research in underwater acoustics [1]. Despite the impressive achievements of the past two decades, the simulation of acoustic fields in large domains with 3D inhomogeneities of various kinds remains a significant challenge, especially when broadband signals are considered. Such simulations are necessary, e.g., for the estimation of anthropogenic acoustic noise levels, which is important for the protection of marine fauna [2]. It is also important that the modeling is carried out in reasonable time (ideally, in real time). Many approaches to the simulation of 3D acoustic fields are being actively developed, including 3D parabolic equations, Gaussian beams [3], rays-theoretical techniques [4], as well as the methods based on the finite-difference or finite-element discretizations [5] of wave equations (see the nice collections of papers in recent special issues [6,7] dedicated to 3D propagation effects and techniques in underwater acoustics and references therein).

Mode parabolic equations (MPE) are a promising and convenient tool for solving 3D sound propagation problems of ocean acoustics. They have been used in adiabatic form, for instance, by Collins [8] and in the refined version by Trofimov [9]. Wide-angle and pseudodifferential MPE were also proposed in recent works [2,10] and used as a basis for the computational code AMPLE [11] that is used for computing noise levels due to broadband sources over large sea areas (e.g., during seismic survey [2]).

In the work in [12,13], this approach has been extended to the case of interacting modes and obtained very good agreement between numerical solutions by MPE and by the source image method in the wedge benchmark problem [14], which was achieved both for up-slope and cross-slope propagation.

The MPE approach involves two separate steps in calculating the acoustic field. First, we calculate the acoustic normal modes along the trace, and then we solve the system of amplitude parabolic equations, in which we obtain the amplitudes of the acoustic modes. The summation of the modes with the calculated amplitudes gives us the desired field. The first step, the computation of acoustic normal modes in waveguides, is a delicate and important task. We should use numerical methods which provide sufficient accuracy along with a high speed of calculations. The number of such methods has been developed in the last few years. For instance, in the article [15] the Chebyshev–Tau spectral method was implemented to solve acoustic normal modes with a stratified ocean [15]. In [16] The Legendre collocation method based on domain decomposition is proposed to calculate normal modes. In [17], the algorithm using the Chebyshev–Tau spectral method was proposed to solve for the horizontal wavenumbers and modes of approximately range-independent segments. In [18], a discrete PE model using the Chebyshev spectral method is derived based on the wide-angle rational approximation. In [19], a multidomain Chebyshev collocation method for the accurate computation of normal modes in open oceanic and atmospheric waveguides was devised. In [20], a three-dimensional, coupled-mode two-way model using the direct-global-matrix technique as well as Fourier synthesis was presented.

For many practical reasons, it would be very important to generalize the MPE theory to the case of elastic bottom. Although, of course, one might perform it by computing elastic normal modes, there is an alternative and a more convenient way to take weak bottom elasticity into account in sound propagation problems of underwater acoustics [21]. The weak elasticity assumption significantly simplifies the derivation of the equations, but the equations obtained in this way can be applicable in practical problems when the shear wave velocity does not exceed the value of about 300 m/s. This may be true for liquid sediments and for sand or clay.

Various parabolic equations for elastic media have been derived in some articles (see [22,23] etc.). As a rule, they have different restrictions and an insufficient amount of test examples. Previously, we have derived an adiabatic mode parabolic equation taking into account the weak elasticity of the bottom [24]. It can be used in the seismoacoustics of liquid sediments when the shear modulus is small.

In this work, we derive a system of parabolic equations with interacting modes, taking into account the small shear modulus at the bottom and the only interface between the water column and elastic bottom. The derivation is based on the so-called generalized multiscale expansions method [25]. The obtained equations are numerically solved by the Crank–Nicolson scheme with iterations and are included into the software package developed for the modeling of acoustical fields using the MPE approach [13]. Then, we compare numerical results obtained by the derived elastic MPE with the source image solutions and validate the efficiency and accuracy of the equations.

2. Basic Equations and Expansions

We consider the propagation of time-harmonic acoustic waves in the three-dimensional waveguide $\Omega =\left\{\right(x,y,z\left)\right|0\le x\le \infty ,-\infty \le y\le \infty ,-H\le z\le 0\}$ with weak elastic bottom, described by the linear elasticity equations [26]

(1)$-\rho {\omega }^2{u}_k=\sum _{j=1}^3\frac{\partial {\sigma }_{kj}}{\partial x_j},k=1,2,3,$

where $\omega$ is the angular frequency of the waves, and the stress tensor of elasticity ${\sigma }_{ij}$ is (here $u_1=u,u_2=v,u_3=w$)

${\sigma }_{ij}=\lambda \left(\sum _{k=1}^3\frac{\partial u_k}{\partial x_k}\right){\delta }_{ij}+\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right).$

Here, $\lambda$ and $\mu$ are the Lame coefficients, $\rho$ is the density of the medium, u, v, and w are displacement components in the Cartesian coordinates with the z-axis directed upwards.

We postulate the following expansions of parameters and displacements

$\begin{array}{cc}\hfill & E_{\mathrm{ef}}=\lambda +2\mu ,E_{\mathrm{ef}}=E_0(X,z)+ϵE_1(X,Y,z),\hfill \\ \hfill & \mu =ϵ{\mu }_1(X,Y,z),\frac{1}{\rho }={\gamma }_0(X,z)+ϵ{\gamma }_1(X,Y,z),\hfill \end{array}$

(2)$\begin{array}{cc}\hfill & u=[u_0(X,Y,z)+ϵu_1(X,Y,z)+\dots ]{\mathrm{e}}^{\mathrm{i}\theta \left(X\right)/ϵ}\hfill \\ \hfill & +u_S(X,Y,z,\eta ){\mathrm{e}}^{\mathrm{i}\theta \left(X\right)/ϵ},\hfill \\ \hfill & v={ϵ}^{1/2}[v_{1/2}(X,Y,z)+ϵv_{3/2}(X,Y,z)+\dots ]{\mathrm{e}}^{\mathrm{i}\theta \left(X\right)/ϵ}\hfill \\ \hfill & +v_S(X,Y,z,\eta ){\mathrm{e}}^{\mathrm{i}\theta \left(X\right)/ϵ},\hfill \\ \hfill & w=[w_0(X,Y,z)+ϵw_1(X,Y,z)+\dots ]{\mathrm{e}}^{\mathrm{i}\theta \left(X\right)/ϵ}\hfill \\ \hfill & +w_S(X,Y,z,\eta ){\mathrm{e}}^{\mathrm{i}\theta \left(X\right)/ϵ}.\hfill \end{array}$

Effective Young’s modulus is $E_{\mathrm{ef}}=\rho c^2$, where c is the velocity of compressional sound waves [27]. Shear modulus is $\mu =\rho c_s^2$, where $c_s$ is the velocity of shear waves. We introduce slow variables $X=ϵx$ and $Y={ϵ}^{1/2}y$, where small parameter $ϵ$ is the ratio of the wavelength to a typical size of horizontal inhomogeneities of the propagation medium. We also introduce fast variable $\eta =\xi (X,Y,z)/\sqrt{ϵ}$ and assume that displacements in shear waves $u_S,v_S,w_S$ also depend on this variable $w_S=w_S(X,Y,z,\eta )$, etc.

The Ansatz (2) is derived from a general form of expansion of the displacements by considering expressions at different powers of $ϵ$ in the basic equations (see Appendix A).

Using the notation introduced above, we expand the elasticity Equation (1). Then, we collect terms with the same powers of $ϵ$. For ${ϵ}^0$, we have the following equality

$\begin{array}{cc}\hfill & {\omega }^2{u}_0=-{\gamma }_0{E}_0\mathrm{i}{\theta }_X(w_{0z}+\mathrm{i}{\theta }_X{u}_0),\hfill \\ \hfill & {\omega }^2{w}_0=-{\gamma }_0{\left[E_0(w_{0z}+\mathrm{i}{\theta }_X{u}_0)\right]}_z.\hfill \end{array}$

Let us call $P=-E_{\mathrm{ef}}(u_x+v_y+w_z)$ quasi-pressure. Performing a multi-scale expansion of P (see [25]), we obtain the following expansion for this quantity $P=(P_0+ϵP_1+\dots ){\mathrm{e}}^{\mathrm{i}\theta \left(X\right)/ϵ}$, where $P_0=-E_0(w_{0z}+\mathrm{i}{\theta }_X{u}_0)$. Then,

(3)${\left({\gamma }_0{P}_{0z}\right)}_z+\left[\frac{{\omega }^2}{c_0^2}-{\left({\theta }_X\right)}^2\right]{\gamma }_0{P}_0=0,$

where $c_0(X,z)=\sqrt{{\gamma }_0{E}_0}$ is the respective zero-order approximation of the sound speed. Next, $P_0$ is sought in the form $P_0=A(X,Y)\varphi (X,z)$, and a Sturm–Liouville (SL) problem is obtained for $\varphi (X,z)$ (called acoustic spectral problem) with $k^2={\left({\theta }_X\right)}^2$ as a spectral parameter. The interface conditions for this SL problem are discussed in Section 4, and the pressure-release boundary conditions have the form

(4)${\left({\gamma }_0{\varphi }_{0z}\right)}_z+\left(\frac{{\omega }^2}{c_0^2}-k^2\right){\gamma }_0{\varphi }_0=0,$

It is known that such an SL problem has a countable set of solutions $(k_j^2,{\varphi }_j)$ indexed by $j=1,2,\dots$, where eigenfunctions ${\varphi }_j$ are orthogonal to each other and can be normalized as

(5)$⟨{\varphi }_i,{\varphi }_j⟩={\int }_{-H}^0{\gamma }_0{\varphi }_i{\varphi }_{j}dz={\delta }_{ij}.$

Each eigenfunction corresponds to a particular solution $P_{0j}=A_j(X,Y){\varphi }_j(X,z)$ of the equation for $P_0$.

Zero-order approximations for the displacements can be also expressed via $P_0$ (index j is omitted):

$u_0=\frac{{\gamma }_0}{{\omega }^2}\mathrm{i}k{P}_0,v_{1/2}=\frac{{\gamma }_0}{{\omega }^2}{P}_{0Y},w_0=\frac{{\gamma }_0}{{\omega }^2}{P}_{0z}.$

Consider now the equations at ${ϵ}^1$. The first-order approximation for the quasi-pressure has the form

$\begin{array}{cc}\hfill & P_1=-E_0(u_{0X}+v_{0Y}+w_{1z}+\mathrm{i}k{u}_1)-E_1(\mathrm{i}k{u}_0+w_{0z}),\hfill \end{array}$

Additionally, the respective equations for the displacement components $u_1$ and $w_1$ write as

$\begin{array}{cc}\hfill & {\omega }^2{u}_1=\mathrm{i}k{\gamma }_0{P}_1+\mathrm{i}k{\gamma }_1{P}_0+{\gamma }_0{P}_{0X}\hfill \\ \hfill & -\frac{\mathrm{i}k}{{\omega }^2}{\gamma }_0[2{\gamma }_0{\mu }_{1z}{P}_{0z}+{\left({\gamma }_{0z}{\mu }_1{P}_0\right)}_z],\hfill \\ \hfill & {\omega }^2{w}_1={\gamma }_0{P}_{1z}+{\gamma }_1{P}_{0z}-\frac{k^2}{{\omega }^2}{\gamma }_0{P}_0(2{\gamma }_0{\mu }_{1z}+{\gamma }_{0z}{\mu }_1).\hfill \end{array}$

Expressing all involved quantities in terms of $P_1$ and $P_0$, we obtain in the usual manner at the $O({ϵ}^1)$ the boundary value problem for $P_1$, which is not always solvable.

(6)$\begin{array}{cc}\hfill & {\left({\gamma }_0{P}_{1z}\right)}_z+\left(\frac{{\omega }^2}{c_0^2}-k^2\right){\gamma }_0{P}_1=k^2{\gamma }_1{P}_0-{\left({\gamma }_1{P}_{0z}\right)}_z\hfill \\ \hfill & -2\mathrm{i}k{\gamma }_0{P}_{0X}-{\gamma }_0{P}_{0YY}-\mathrm{i}k{\gamma }_{0X}{P}_0-\mathrm{i}{k}_X{\gamma }_0{P}_0\hfill \\ \hfill & +\frac{E_1{\omega }^2}{E_0{c}_0^2}{\gamma }_0{P}_0+\frac{k^2}{{\omega }^2}\left[2{\left({\gamma }_0^2{\mu }_{1z}\right)}_z+{\gamma }_{0z}^2{\mu }_1\right]P_0=F.\hfill \end{array}$

3. WKB Solutions for Shear Waves

It is interesting to note that expressions for the displacements in the shear waves in WKB approximation are obtained automatically in the frame of the used multi-scale method. In comparison with the usual MPE, when we use parabolic scaling with the asymptotic variables $X=ϵx$ and $Y=\sqrt{ϵ}y$ and phase $\zeta =\theta \left(X\right)/ϵ$ in the case of the weak elastic media, we should introduce one more phase variable.

Consider the wave equation for the vertically propagating shear wave in the media with constant $c_s$: $u_{zz}+({\omega }^2/c_s^2)u=0$. However, $c_s^2=\mu /\rho$ and $\mu =ϵ{\mu }_1$. So, we obtain equation $ϵu_{zz}+({\omega }^2\rho /{\mu }_1^2)u=0$. For this equation, we can obtain the approximate solution by the classical WKB method even without the assumption of constant $c_s$. In the multi-scale approach, we can remove $ϵ$ from this equation by introducing the new variable $\eta =z/\sqrt{ϵ}$.

In our more general case, when $c_s$ and other parameters are not constants, we should introduce the phase variable $\eta =\xi (X,Y,z)/\sqrt{ϵ}$. Thus, in our considerations, only the powers of $\sqrt{ϵ}$ take parts in the asymptotic variables. Then, it is reasonable to seek solutions for the displacements u, v, and w as the asymptotic series in the power of $\sqrt{ϵ}$. Additionally, dependent variables should depend on X, Y, z, $\zeta$, and $\eta$. Successive consideration of the terms of these series is performed in the Appendix A.

For the shear waves displacements, we use Expansion (A3) obtained as a result of the considerations in the Appendix A.

$\begin{array}{cc}\hfill & u_S={ϵ}^{1/2}{u}_{S1/2}(X,Y,z,\eta )+ϵu_{S1}(X,Y,z,\eta )+\dots ,\hfill \\ \hfill & v_S=ϵv_{S1}(X,Y,z,\eta )+{ϵ}^{3/2}{v}_{S3/2}(X,Y,z,\eta )+\dots ,\hfill \\ \hfill & w_S=ϵw_{S1}(X,Y,z,\eta )+{ϵ}^{3/2}{w}_{S3/2}(X,Y,z,\eta )+\dots .\hfill \end{array}$

Let us substitute them into the elasticity Equation (1) and collect terms at the same powers of $ϵ$. Note that subscript j is omitted everywhere in this section.

At ${ϵ}^{1/2}$ in the equation for $u_S$ and at ${ϵ}^0$ in the equation for $w_S$, we obtain:

(7)$\begin{array}{cc}\hfill & {\gamma }_0{\mu }_1{\xi }_z^2{u}_{S1/2,\eta \eta }+{\omega }^2{u}_{S1/2}=0,\hfill \\ \hfill & {\xi }_z{w}_{S1,\eta }+\mathrm{i}k{u}_{S1/2}=0.\hfill \end{array}$

To solve these equations, we put $u_{S1/2}=\tilde{u}(X,Y,z)exp\left(\mathrm{i}\eta \right)$ and $w_{S0}=\tilde{w}(X,Y,z)$$exp\left(\mathrm{i}\eta \right)$:

$(-{\gamma }_0{\mu }_1{\xi }_z^2+{\omega }^2)\tilde{u}=0,\mathrm{i}({\xi }_z\tilde{w}+k\tilde{u})=0.$

Thus, we obtain the expression for phase variable $\xi (X,Y,z)$:

${\xi }_z^2=\frac{{\omega }^2}{{\gamma }_0{\mu }_1}\Rightarrow \xi =\pm {\int }_{-H}^{z}q(X,Y,s)ds=\pm \left|\xi \right|,$

where $q=\omega /\sqrt{{\gamma }_0{\mu }_1}$. Moreover, $\tilde{w}=-k\tilde{u}/{\xi }_z$. From the equations at ${ϵ}^1$ and ${ϵ}^{1/2}$ for $u_S$ and $w_S$, we obtain:

(8)$(2{\xi }_z{\tilde{u}}_z+{\xi }_{zz}\tilde{u}){\mu }_1+{\xi }_z{\mu }_{1z}\tilde{u}=0.$

This equation with separable variables can be easily solved

$2\frac{{\tilde{u}}_z}{\tilde{u}}+\frac{{\xi }_{zz}}{{\xi }_z}+\frac{{\mu }_{1z}}{{\mu }_1}=0\Rightarrow \tilde{u}=\frac{C(X,Y)}{\sqrt{|{\xi }_z{\mu }_1|}}.$

Finally, we have solutions for $u_{S1/2}$ and $w_{S1}$:

$\begin{array}{cc}\hfill & u_{S1/2}=\frac{1}{\sqrt{q{\mu }_1}}\left(C_1{\mathrm{e}}^{\mathrm{i}\left|\xi \right|/\sqrt{ϵ}}+C_2{\mathrm{e}}^{-\mathrm{i}\left|\xi \right|/\sqrt{ϵ}}\right),\hfill \\ \hfill & w_{S1}=-\frac{k}{q\sqrt{q{\mu }_1}}\left(C_1{\mathrm{e}}^{\mathrm{i}\left|\xi \right|/\sqrt{ϵ}}-C_2{\mathrm{e}}^{-\mathrm{i}\left|\xi \right|/\sqrt{ϵ}}\right).\hfill \end{array}$

The expressions at ${ϵ}^{3/2}$ for $u_S$ and $v_S$ give us the equation for $\tilde{v}$ in $v_{S1}=\tilde{v}(X,Y,z)$$exp\left(\mathrm{i}\eta \right)$, which is identical to Equation (8) for $\tilde{u}$, and therefore we have

$v_{S1}=\frac{1}{\sqrt{q{\mu }_1}}\left(C_3{\mathrm{e}}^{\mathrm{i}\left|\xi \right|/\sqrt{ϵ}}+C_4{\mathrm{e}}^{-\mathrm{i}\left|\xi \right|/\sqrt{ϵ}}\right).$

We also obtained the equation for ${\tilde{u}}_1$ in $u_{S1}={\tilde{u}}_1(X,Y,z)exp\left(\mathrm{i}\eta \right)$:

$\begin{array}{cc}\hfill & (2{\xi }_z{\tilde{u}}_{1,z}+{\xi }_{zz}{\tilde{u}}_1){\mu }_1+{\xi }_z{\mu }_{1z}{\tilde{u}}_1=\hfill \\ \hfill & [k^2+{\xi }_Y^2+{\xi }_z^2({\gamma }_1/{\gamma }_0)]{\mu }_1\tilde{u}-{\left({\mu }_1{\tilde{u}}_z\right)}_z.\hfill \end{array}$

${\tilde{w}}_1$ can be found from

${\xi }_z{\tilde{w}}_1=\mathrm{i}{\tilde{w}}_z-k{\tilde{u}}_1-{\xi }_X\tilde{u}-{\xi }_Y\tilde{v},$

which also appeared when collecting terms of the order ${ϵ}^{1/2}$.

4. Interface Conditions

We consider the only interface between the water layer and the elastic bottom. At the surface described by equation z = h(x, y), we introduce the interface conditions:

(9)$\begin{array}{cc}\hfill & ({\mathbf{u}}_{+}-{\mathbf{u}}_{-})·\mathbf{n}=0,\mathbf{n}·({\mathbf{\sigma }}_{+}-{\mathbf{\sigma }}_{-})\mathbf{n}=0,\hfill \\ \hfill & {\mathbf{t}}_{\mathbf{1}}·\left({\mathbf{\sigma }}_{-}\right)\mathbf{n}=0,{\mathbf{t}}_{\mathbf{2}}·\left({\mathbf{\sigma }}_{-}\right)\mathbf{n}=0,\hfill \end{array}$

where n is a unit normal vector to the surface, and the subscript ‘+’ denotes upper limit of the respective quantities $z=z_0$ (i.e., the limit $z\to z_0+0$), while the subscript ‘−’ corresponds to the limit $z\to z_0-0$, and $\mathbf{u}$ is the particle velocity vector.

The normal vector $\mathbf{n}$ has coordinates

$\mathbf{n}=(-h_x,-h_y,1)·{(1+h_x^2+h_y^2)}^{-1/2},$

while tangent vectors ${\mathbf{t}}_{\mathbf{1}}$ and ${\mathbf{t}}_{\mathbf{2}}$ can be written as having the components:

$\begin{array}{cc}\hfill & {\mathbf{t}}_{\mathbf{1}}=(1,0,h_x)·{(1+h_x^2)}^{-1/2},\hfill \\ \hfill & {\mathbf{t}}_{\mathbf{2}}=(-h_x{h}_y,1+h_x^2,h_y)\hfill \\ \hfill & \times {(1+h_x^2)}^{-1/2}{(1+h_x^2+h_y^2)}^{-1/2}.\hfill \end{array}$

We further assume that surface function $h(x,y)$ can be represented in the form

$h=h_0\left(X\right)+ϵh_1(X,Y),$

and therefore $h_x=ϵh_{0X}+{ϵ}^2{h}_{1X}$ and $h_y={ϵ}^{3/2}{h}_{0Y}$.

The interface conditions for the stress tensor modulo terms of the order $O({ϵ}^2)$ are as follows:

(10)${\left[{\sigma }_{33}\right]}_{-}={\left[{\sigma }_{33}\right]}_{+},{\left[{\sigma }_{13}\right]}_{-}=0,{\left[{\sigma }_{23}\right]}_{-}=0.$

Using the fact that the interface perturbation $ϵh_1$ is small, we transfer the interface conditions to the unperturbed surface $z=h_0$, retaining only the terms up to the order $O(ϵ)$.

Thus, we obtain the following Taylor expansion:

$\begin{array}{cc}\hfill & {\left.\mathbf{u}(X,Y,z,\eta )\right|}_{z=h_0+ϵh_1}=\hfill \\ \hfill & {\left.\mathbf{u}(X,Y,z,\eta )\right|}_{z=h_0}+{\left.ϵh_1\mathbf{u}{(X,Y,z,\eta )}_z\right|}_{z=h_0}+\dots ,\hfill \\ \hfill & {\left.\mathbf{\sigma }(X,Y,z,\eta )\right|}_{z=h_0+ϵh_1}=\hfill \\ \hfill & {\left.\mathbf{\sigma }(X,Y,z,\eta )\right|}_{z=h_0}+{\left.ϵh_1\mathbf{\sigma }{(X,Y,z,\eta )}_z\right|}_{z=h_0}+\dots .\hfill \end{array}$

With all the mentioned assumptions, the interface Condition (10) at $z=h_0\left(X\right)$ can be rewritten in terms of $P_0$ and $P_1$ in the following form:

(11)$\begin{array}{cc}\hfill & {\left(P_0\right)}_{+}-{\left(P_0\right)}_{-}=0,\hfill \\ \hfill & {\left(P_1\right)}_{+}-{\left(P_1\right)}_{-}=-h_1{\left({\gamma }_0{P}_{0z}\right)}_{+}\left({\gamma }_{0+}^{-1}-{\gamma }_{0-}^{-1}\right)\hfill \\ \hfill & -\frac{2k^2}{{\omega }^2}{\left({\mu }_1{\gamma }_0{P}_0\right)}_{-},\hfill \end{array}$

(12)$\begin{array}{cc}\hfill & {\left({\mu }_1{\xi }_z{u}_{S0\eta }\right)}_{-}=-{\left[\frac{\mathrm{i}k{\mu }_1}{{\omega }^2}\left(2{\gamma }_0{P}_{0z}+{\gamma }_{0z}{P}_0\right)\right]}_{-},\hfill \\ \hfill & {\left({\mu }_1{\xi }_z{v}_{S0\eta }\right)}_{-}=-{\left[\frac{{\mu }_1}{{\omega }^2}{\left(2{\gamma }_0{P}_{0z}+{\gamma }_{0z}{P}_0\right)}_Y\right]}_{-}.\hfill \end{array}$

The interface condition for the normal displacement component from (9) is

$\begin{array}{cc}\hfill & {[w_0+ϵ(w_1-h_{0X}{u}_0+h_1{w}_{0z}+w_{S0})]}_{-}=\hfill \\ \hfill & {[w_0+ϵ(w_1-h_{0X}{u}_0+h_1{w}_{0z})]}_{+}.\hfill \end{array}$

Let us express this condition through $P_0$ and $P_1$:

$\begin{array}{cc}\hfill & {\left({\gamma }_0{P}_{0z}\right)}_{-}={\left({\gamma }_0{P}_{0z}\right)}_{+},\hfill \\ \hfill & \left[{\gamma }_0{P}_{1z}+{\gamma }_1{P}_{0z}-\frac{k^2}{{\omega }^2}{\gamma }_0{P}_0(2{\gamma }_0{\mu }_{1z}+{\gamma }_{0z}{\mu }_1)\right.\hfill \\ \hfill & {\left.-\mathrm{i}k{h}_{0X}{\gamma }_0{P}_0+h_1{\left({\gamma }_0{P}_{0z}\right)}_z+{\omega }^2{w}_{S0}\right]}_{-}\hfill \\ \hfill & ={\left[{\gamma }_0{P}_{1z}+{\gamma }_1{P}_{0z}-\mathrm{i}k{h}_{0X}{\gamma }_0{P}_0+h_1{\left({\gamma }_0{P}_{0z}\right)}_z\right]}_{+}.\hfill \end{array}$

From (7), we can express ${\omega }^2{w}_{S0}=\mathrm{i}k{\gamma }_0{\mu }_1{\xi }_z{u}_{S0\eta }$ and, considering the first equation in (12), we obtain ${({\omega }^2{w}_{S0})}_{-}={[{\gamma }_0{\mu }_1{k}^2(2{\gamma }_0{P}_{0z}+{\gamma }_{0z}{P}_0)/{\omega }^2]}_{-}$, and finally obtain:

(13)$\begin{array}{cc}\hfill & {\left({\gamma }_0{P}_{0z}\right)}_{+}-{\left({\gamma }_0{P}_{0z}\right)}_{-}=0,\hfill \\ \hfill & {\left({\gamma }_0{P}_{1z}\right)}_{+}-{\left({\gamma }_0{P}_{1z}\right)}_{-}=\mathrm{i}k{h}_{0X}{P}_{0+}({\gamma }_{0+}-{\gamma }_{0-})\hfill \\ \hfill & -h_1[{\left\{{\left({\gamma }_0{P}_{0z}\right)}_z\right\}}_{+}-{\left\{{\left({\gamma }_0{P}_{0z}\right)}_z\right\}}_{-}]\hfill \\ \hfill & -{\left({\gamma }_0{P}_{0z}\right)}_{+}\left[{\left(\frac{{\gamma }_1}{{\gamma }_0}\right)}_{+}-{\left(\frac{{\gamma }_1}{{\gamma }_0}\right)}_{-}\right]\hfill \\ \hfill & +\frac{2k^2}{{\omega }^2}{\left[{\gamma }_0^2({\mu }_1{P}_{0z}-{\mu }_{1z}{P}_0)\right]}_{-}.\hfill \end{array}$

This expression contains only $P_0$ and $P_1$ variables, so we can consider a pure acoustical case with the perturbations due to shear waves.

5. Boundary Conditions

At the upper boundary z = 0, we postulate $\sigma \mathbf{n}=0$, which is reduced to

${\sigma }_{33}=-P_0{\mathrm{e}}^{\mathrm{i}\theta /ϵ}-ϵP_1{\mathrm{e}}^{\mathrm{i}\theta /ϵ}=0.$

Thus, we have Dirichlet conditions or pressure release conditions $P_0=0$ and $P_1=0$.

At the lower boundary $z=-H$, we can choose $\sigma \mathbf{n}=0$, and therefore

${\sigma }_{33}=-P_0{\mathrm{e}}^{\mathrm{i}\theta /ϵ}-ϵ[P_1-2(k^2/{\omega }^2){\mu }_1{\gamma }_0{P}_0]{\mathrm{e}}^{\mathrm{i}\theta /ϵ}=0.$

Again, we have pressure release conditions for $P_0=0$ and $P_1=0$. Moreover, ${\sigma }_{13}=0$ at $z=-H$ reduces to the following equality:

$\left({\mu }_1{\xi }_z{u}_{S0\eta }\right)=-\frac{2\mathrm{i}k}{{\omega }^2}\left({\mu }_1{\gamma }_0{P}_{0z}\right).$

In a similar way, the condition ${\sigma }_{23}=0$ at $z=-H$ reduces to:

$\left({\mu }_1{\xi }_z{v}_{S0\eta }\right)=-\frac{2}{{\omega }^2}\left({\mu }_1{\gamma }_0{P}_{0Yz}\right).$

More useful types of boundary conditions at the lower boundary for comparison of our solutions with the solutions of ASA wedge wave propagation problems are transparent-like conditions at $z=-H$, where we require

$u_{S0}\sim {\mathrm{e}}^{-\mathrm{i}\left|\xi \right|/\sqrt{ϵ}},w_{S0}\sim {\mathrm{e}}^{-\mathrm{i}\left|\xi \right|/\sqrt{ϵ}},v_{S0}\sim {\mathrm{e}}^{-\mathrm{i}\left|\xi \right|/\sqrt{ϵ}},$

Which can be satisfied explicitly by assuming $C_1=0$ and $C_3=0$. We also assume $P_0=0$ and $P_1=0$ and introduce an absorbing layer in order to suppress reflections of compressional waves from bottom.

6. Elastic Mode Parabolic Equations (EMPE)

Let $\left\{{\theta }_j\right|j=M,\cdots ,N\}$ be a set of phases. We introduce multi-mode representation following [13]:

$\begin{array}{cc}\hfill & \mathcal{P}={\mathcal{P}}_0+ϵ{\mathcal{P}}_1+\dots \hfill \\ \hfill & =\sum _{j=M}^N\left[P_{0j}(X,Y,z)+ϵP_{1j}(X,Y,z)+\dots \right]{\mathrm{e}}^{\mathrm{i}{\theta }_j/ϵ}.\hfill \end{array}$

Which gives a solution for the compressional acoustic waves. Here,

$\begin{array}{cc}\hfill & P_{0j}=A_j(X,Y){\varphi }_j(z,X),\hfill \\ \hfill & P_{1j}=\sum _{l=0}^{\infty }{B}_{jl}(X,Y){\varphi }_l(z,X),\hfill \end{array}$

where $B_{jl}=⟨P_{1j},{\varphi }_l⟩$ (projections obtained via inner product). Eigenfunctions ${\varphi }_j$ with ${\theta }_j$ are obtained from the spectral Problem (3). Amplitudes $A_j$ can be found from the solvability condition for the boundary value Problem (6) for $P_{1j}$. For this following [13], we multiply the expression for ${\mathcal{P}}_1$ by ${\varphi }_l$ and then integrate the resulting equation from $-H$ to 0 by parts twice with the use of the interface Conditions (11) and (13).

$\begin{array}{cc}\hfill & \sum _{j=M}^N{\mathrm{e}}^{\mathrm{i}{\theta }_j/ϵ}{\int }_{-H}^0\left[{\left({\gamma }_0{\varphi }_{lz}\right)}_z+\left(\frac{{\omega }^2}{c_0^2}-k_j^2\right){\gamma }_0{\varphi }_l\right]P_{1j}dz\hfill \\ \hfill & =\sum _{j=M}^N{\mathrm{e}}^{\mathrm{i}{\theta }_j/ϵ}{\int }_{-H}^0{F}_j{\varphi }_{l}dz+\sum _{j=M}^N{\mathrm{e}}^{\mathrm{i}{\theta }_j/ϵ}\hfill \\ \hfill & \times \{{\varphi }_l\left[{\left({\gamma }_0{P}_{1jz}\right)}_{+}-{\left({\gamma }_0{P}_{1jz}\right)}_{-}\right]\hfill \\ \hfill & -{\gamma }_0{\varphi }_{lz}\left[{\left(P_{1j}\right)}_{+}-{\left(P_{1j}\right)}_{-}\right]+(k_j^2-k_l^2)B_{jl}\}.\hfill \end{array}$

The left part of this equality turns zero due to (3), and the right part represents the desired parabolic equation. The terms $(k_j^2-k_l^2)B_{jl}$ in these expressions can be omitted because of the resonant condition $|k_j-k_l|=O\left(ϵ\right)$ [13].

As a result, we obtain a system of parabolic wave equations for $l=M,\dots ,N$

(14)$2\mathrm{i}{k}_l{A}_{l,X}+\mathrm{i}{k}_{l,X}A+A_{l,YY}+\sum _{j=M}^N{\beta }_{lj}{A}_{j}exp\left({\theta }_{lj}\right)=0,$

where ${\theta }_{lj}=\mathrm{i}({\theta }_l-{\theta }_j)/ϵ$, and ${\beta }_{lj}$ is given by the formula

$\begin{array}{cc}\hfill & {\beta }_{lj}={\int }_{-H}^0{\gamma }_0\nu {\varphi }_j{\varphi }_{l}dz+{\int }_{-H}^0{\gamma }_1\left(n_0^2-k_j^2\right){\varphi }_j{\varphi }_{l}dz\hfill \\ \hfill & -{\int }_{-H}^0{\gamma }_1{\varphi }_{jz}{\varphi }_{lz}dz-\mathrm{i}{k}_j(C_{lj}-C_{jl})\hfill \\ \hfill & -\frac{k_j^2}{{\omega }^2}{\int }_{-H}^0[2{\left({\gamma }_0^2{\mu }_{1z}\right)}_z+{\gamma }_{0z}^2{\mu }_1]{\varphi }_j{\varphi }_{l}dz\hfill \\ \hfill & +\{h_1{\varphi }_j{\varphi }_l\left[k_j^2\left({\gamma }_{0+}-{\gamma }_{0-}\right)-{\left(n_0^2{\gamma }_0\right)}_{+}+{\left(n_0^2{\gamma }_0\right)}_{-}\right]\hfill \\ \hfill & -h_1{\gamma }_0^2{\varphi }_{jz}{\varphi }_{lz}\left({\gamma }_{0+}^{-1}-{\gamma }_{0-}^{-1}\right)\hfill \\ \hfill & {\left.-\frac{2k_j^2}{{\omega }^2}{\gamma }_0^2[{\mu }_1{\left({\phi }_l{\phi }_j\right)}_z-{\mu }_{1z}{\phi }_l{\phi }_j]\}\right|}_{z=h_0}.\hfill \end{array}$

Here, the values $\nu$ and $n_0$ are the same as in the study [24], coefficients $C_{lj}$ have been obtained in the work [13], and ${\beta }_{lj}$ differs from ${\alpha }_{lj}$ from the latter paper only by two terms containing ${\mu }_1$. Thus, in the case of low shear velocities, the problem reduces to the acoustic case with the correction in the form of these terms.

7. Initial-Boundary Value Problem for the System of Elastic Mode Parabolic Equations

We refer to the quantities and variables X, Y, $\nu$, ${\gamma }_0$${\mu }_1$, ${\theta }_j$, $u_{S1/2j}$, $v_{S1j}$, $w_{S1j}$, and $A_j$ as the asymptotic ones. Considerations of initial-boundary value problems in a partially bounded domain require the use of physical quantities and variables, which are x, y, $\overline{\nu }=ϵ\nu$, $\overline{\gamma }={\gamma }_0$$\overline{\mu }=ϵ{\mu }_1$, ${\overline{u}}_{Sj}={ϵ}^{1/2}{u}_{S1/2j}$, ${\overline{v}}_{Sj}=ϵv_{S1j}$, ${\overline{w}}_{Sj}=ϵw_{S1}$, ${\overline{A}}_j(x,y)=A_j(ϵx,\sqrt{ϵ}y)=A_j(X,Y)$, and ${\displaystyle {\overline{\theta }}_j={\int }^X\frac{1}{ϵ}{\theta }_{j,X}dX={\int }^x{k}_{j}dx}$, ${\overline{c}}_s=\sqrt{\overline{\gamma }\overline{\mu }}$. It can be easily verified that Equation (14) in physical variables are read as

(15)$2\mathrm{i}{k}_l{\overline{A}}_{l,x}+\mathrm{i}{k}_{l,x}{\overline{A}}_j+{\overline{A}}_{l,yy}+{\overline{\beta }}_{ll}{\overline{A}}_l+\sum _{j=M,j\ne l}^N{\overline{\beta }}_{lj}{\overline{A}}_{j}exp\left({\overline{\theta }}_{lj}\right)=0,$

where ${\overline{\beta }}_{lj}$ are computed by the same formulae as ${\beta }_{lj}$, with $\nu$ replaced by $\overline{\nu }$, ${\overline{\theta }}_{lj}=({\overline{\theta }}_l-{\overline{\theta }}_j)$.

For Equation (15), we consider the initial-boundary value problem in the domain of the form $\left\{(x,y)\right|0\le x<\infty ,-y_0/2\le y\le y_0/2\}$, with the initial condition

(16)${\overline{A}}_j(0,y)=g_j\left(y\right)=(g,{\varphi }_j)={\int }_{-H}^0{\gamma }_{0}g(z,y){\varphi }_j\left(z\right)dz,$

and the transparent boundary conditions at $y=-y_0/2$ and $y=y_0/2$ (i.e., requiring that only outgoing waves are present there). In practice, this condition is ensured by the use of perfectly matched layers.

Now, we formulate the interface conditions for the displacements in the case of the absorbing conditions at the bottom.

(17)${\left(u_{S0j}\right)}_{-}=\frac{A_j{k}_j}{{\omega }^3}{\left[\sqrt{{\gamma }_0{\mu }_1}\left(2{\gamma }_0{\varphi }_{jz}+{\gamma }_{0z}{\varphi }_j\right)\right]}_{-}$

(18)${\left(v_{S0j}\right)}_{-}=\frac{A_{jY}}{{\omega }^3}{\left[\sqrt{{\gamma }_0{\mu }_1}\left(2{\gamma }_0{\varphi }_{jz}+{\gamma }_{0z}{\varphi }_j\right)\right]}_{-}$

(19)${\left(w_{S0j}\right)}_{-}=\frac{A_j{k}_j^2}{{\omega }^4}{\left[{\gamma }_0{\mu }_1\left(2{\gamma }_0{\varphi }_{jz}+{\gamma }_{0z}{\varphi }_j\right)\right]}_{-}$

We have expressions for the displacements for the considered case:

$\begin{array}{cc}\hfill & u_{S0j}=C_2{\left(\frac{{\gamma }_0}{{\omega }^2{\mu }_1}\right)}^{1/4}exp(-\mathrm{i}\eta ),v_{S0j}=C_4{\left(\frac{{\gamma }_0}{{\omega }^2{\mu }_1}\right)}^{1/4}exp(-\mathrm{i}\eta ),\hfill \\ \hfill & w_{S0j}=k_j{C}_2{\left(\frac{{\gamma }_0^3{\mu }_1}{{\omega }^6}\right)}^{1/4}exp(-\mathrm{i}\eta ),\eta =\frac{1}{ϵ}{\int }_{-H}^z\frac{\omega }{\sqrt{{\gamma }_0{\mu }_1}}dz.\hfill \end{array}$

Final expressions for the displacements in the shear waves in the physical variables are

$\begin{array}{cc}\hfill & {\overline{u}}_{S0j}=\frac{{\overline{A}}_j{k}_j}{{\omega }^2}\Phi (z,h){\left(\frac{\overline{\gamma }}{{\overline{c}}_s}\right)}^{1/2}{\left[\sqrt{\overline{\gamma }}{\overline{c}}_s^{3/2}\left(2{\varphi }_{jz}+\frac{{\overline{\gamma }}_z}{\overline{\gamma }}{\varphi }_j\right)\right]}_{-},\hfill \\ \hfill & {\overline{v}}_{S0j}=\frac{{\overline{A}}_{j,y}}{{\omega }^2}\Phi (z,h){\left(\frac{\overline{\gamma }}{{\overline{c}}_s}\right)}^{1/2}{\left[\sqrt{\overline{\gamma }}{\overline{c}}_s^{3/2}\left(2{\varphi }_{jz}+\frac{{\overline{\gamma }}_z}{\overline{\gamma }}{\varphi }_j\right)\right]}_{-},\hfill \\ \hfill & {\overline{w}}_{S0j}=\frac{A_j{k}_j^2}{{\omega }^3}\Phi (z,h){\left(\overline{\gamma }{\overline{c}}_s\right)}^{1/2}{\left[\sqrt{\overline{\gamma }}{\overline{c}}_s^{3/2}\left(2{\varphi }_{jz}+\frac{{\overline{\gamma }}_z}{\overline{\gamma }}{\varphi }_j\right)\right]}_{-},\hfill \\ \hfill & \Phi (z,h)=exp\left(\mathrm{i}{\int }_z^{-h\left(x\right)}\frac{\omega }{{\overline{c}}_s}dz\right).\hfill \end{array}$

8. Numerical Examples

The developed model based on EMPE (15) can handle any waveguides where horizontal changes of bathymetry and hydrology are slow in comparison with the wavelength. Nevertheless, for the test calculations, it is convenient to choose a sufficiently simple waveguide with the elastic bottom for which the exact solutions, or solutions, obtained by essentially different methods are known. In our opinion, the most popular of such waveguides is the ASA (Acoustical Society of America) wedge benchmark, for which we can obtain solutions for the acoustic fields by the source image method, even in the case of elastic bottom.

In order to illustrate the accuracy of the solution obtained using the equations derived in this study, we perform a numerical simulation of sound propagation in the standard ASA wedge benchmark problem (Figure 1) with the bottom slope angle $\approx 2.{86}^{\circ }$ [1,14]. Moreover, in some calculations, we used a wedge with an angle $\approx 1.{14}^{\circ }$. The sound speed in the water is 1500 m/s. The sound speed at the bottom is 1700 m/s. The bottom density is 1500 kg/m³, while the water density is 1000 kg/m ${}^3$. We assume that the seawater is a lossless medium, while at the bottom the attenuation is 0.5 dB/ $\lambda$. The point source with frequency 25 Hz is located at the depth of 100 m, and the depth of the receiver is 30 m. The water depth is 200 m at the source location. For calculation purposes, we truncate the computational domain at the depth of 1500 m below the sea surface. Shear speed at the bottom in the numerical experiments considered below is varied in the range $c_s$ = 0–500 m/s. The attenuation of the shear waves at the bottom is assumed to be zero.

Figure 2 shows transmission loss in the x–z plane obtained by the EMPE method and the first 10 modes in this case at distance x = 1000 m from the source.

Figure 3 illustrates the comparison between transmission losses obtained by the EMPE method and by the source images method for $c_s=350 \mathrm{m}/\mathrm{s}$ (top) and $c_s=500 \mathrm{m}/\mathrm{s}$ (bottom). Moreover, for the case with $c_s=350 \mathrm{m}/\mathrm{s}$, we additionally assumed bottom attenuation to be $0.65 \mathrm{d}\mathrm{B}/\lambda$ to illustrate the effect found in computational experiments increasing the accuracy of calculations with the introduction of additional wave attenuation at the bottom. The curves for the case $c_s=500 \mathrm{m}/\mathrm{s}$ (bottom) illustrate how EMPE fails at large values of $c_s$. The comparison of the transmission loss obtained by EMPE with the one obtained by the source image method reveals the discrepancy of the curves at distances $x>2$ km.

In Figure 4, we estimate the dependence of the meansquare errors of calculations from the bottom shear speed, when $c_s$ = 0–500 m/s. We considered it along wedge propagation for the wedge angle $\approx 2.{86}^{\circ }$ (left) and for the wedge angle $\approx 1.{14}^{\circ }$ (right). One can see that in both cases of EMPE, taking into account the weak elasticity of the bottom provides better accuracy for $c_s<270$ m/s compared with the simple MPE without such consideration.

In Figure 5, we investigate the influence on the accuracy of calculations by an additional small attenuation at the bottom, which could be useful for future work. Moreover, we illustrate the influence of the type of the boundary condition (Dirichlet or Neumann) at the lower boundary of the area on the accuracy of calculations.

9. Discussion

In the the case of cross-wedge propagation, we discovered discrepancies at large distances between solutions obtained by the source image method and by EMPE that can be explained by the narrow-angle nature of the derived EMPE, which leads to additional errors. These errors can be significantly reduced by using wide-angle equations that take into account the elastic properties of the bottom (e.g., similar to the ones discussed in [10]). Moreover, such equations can be derived within the framework of the formalism presented in this study by considering terms arising at higher powers of a small parameter $ϵ$. This can be a matter for future work.

The calculation time for the classical ASA wedge benchmark problem with 44 propagation modes was about 426 s on the one core of the Intel ${}^{\circledR }$ Core™ i5-4690 K CPU, 3.50 GHz. Out of this time, 59 s were required to calculate all 44 modes along the 4 km trace with the horizontal step of 4 m (1000 steps). the vertical step was 1 m with the overall waveguide depth H = 1.5 km (1500 dots). For calculations of modes, we used the method of the inverse iteration with the first order Richardson extrapolation to increase the accuracy. To calculate the wavenumbers along the trace, the bisection algorithm was used. Of course mode calculations can be significantly accelerated with the use of parallel computing. The EMPE System (14) was solved numerically by the Crank–Nicholson implicit difference scheme in combination with the Gauss–Seidel iteration method. The step on variable y was 3 m, with the width of the waveguide being 3 km, so the calculation grid was $1000\times 1000$ nodes. To prevent reflections from the boundaries, we used transparent boundary conditions at $y=\pm 1500 \mathrm{m}$. As transparent boundary conditions, we used the perfectly matched layers from [28]. Reflections from the bottom were suppressed by the artificially absorbing layer beginning from the depth of 1 km. The time of calculation for the EMPE system was about 367 s. Unfortunately, at this stage, calculations are hard to parallelize.

The above examples show that the elastic mode parabolic Equation (14), despite the relatively small influence of the shear modulus, has a solution different from the purely acoustic case. The derived equations are applicable in our test cases up to a shear velocity ≈300 m/s.

In our opinion, the MPEs used for the case of a weak elastic bottom developed in this study are the most computationally efficient tool for solving 3D sound propagation problems in shallow-water acoustics. This opinion may be justified by the fact that MPEs can be considered a result of model order reduction procedure as compared with 3D parabolic equations known from the literature (see, e.g., [29,30,31] and references therein).

Although one should expect that sound propagation models based on ray theory or Gaussian beams are even more efficient [3,4], they can be somewhat less accurate for low frequencies and relatively small water depth.

10. Conclusions

• In this work, with the use of the approach in articles [9,13,24], a mode parabolic equation method for resonantly interacting modes accounting for weak elasticity at the bottom was developed.

• The proposed method is numerically validated. The test calculations carried out for the ASA wedge benchmark prove to be in excellent agreement with the source image method [14] for shear wave speeds up to $c_s\approx$ 300 m/s at the bottom and a rather good agreement of up to $c_s\approx$ 400 m/s.

• We have developed the software package [32] for the modeling of sound propagation in 3D waveguides based on the derived equations.

• This software package has been successfully used to plan and analyze the results of the acoustic experiments on the propagation of sound in a shallow sea [33].

Institutional Review Board Statement

Not applicable.

Data Availability Statement

Data sharing not applicable.

Conflicts of Interest

The author declares no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. ASA wedge with the wedge angle [Forumla omitted. See PDF.].

Figure 2. Transmission loss in the x–z plane obtained by the EMPE method for [Forumla omitted. See PDF.] m/s and bottom attenuation 0.5 dB/[Forumla omitted. See PDF.] (top). The first 10 modes in this case at distance X = 1000 m from the source (bottom).

Figure 3. Transmission loss obtained by the elastic mode parabolic equation method compared with the source images method. Shear waves velocity [Forumla omitted. See PDF.] m/s (top) and [Forumla omitted. See PDF.] m/s (bottom), attenuation is 0.5 dB/[Forumla omitted. See PDF.].

Figure 4. Dependence of the errors of calculations from the bottom shear speed, when [Forumla omitted. See PDF.] = 0–500 m/s. Along wedge propagation for the wedge angle [Forumla omitted. See PDF.] (left) and for the wedge angle [Forumla omitted. See PDF.] (right).

Figure 5. Influence on the accuracy of calculations by an additional small attenuation at the bottom for [Forumla omitted. See PDF.] m/s, zero bottom attenuation (left), and 0.5 dB/[Forumla omitted. See PDF.] (right).

Appendix A. Derivation of the Main Ansatz

Here, we derive ansatz for our problem from the general form of displacement component expansion. $$E_{\mathrm{ef}}=\lambda +2\mu ,E_{\mathrm{ef}}=E_0(X,z)+ϵE_1(X,Y,z),\mu =ϵ{\mu }_1(X,Y,z),$$ $$\frac{1}{\rho }={\gamma }_0(X,z)+ϵ{\gamma }_1(X,Y,z),$$ $$\mathbf{U}=\sum _{j=0}^{\infty }{ϵ}^{j/2}{\mathbf{U}}_{j/2}(X,Y,z,\eta ,\zeta ),$$ where $\mathbf{U}=(U,V,W)$ is the vector of displacements, $\eta =\xi (X,Y,z)/\sqrt{ϵ}$, and $\zeta =\theta (X,Y,z)/ϵ$. Now, we expand the equations of elasticity and consider terms at the different powers of $ϵ$.

At $O({ϵ}^{-2}),O({ϵ}^{-3/2})$, we have:$$\begin{array}{cc}\hfill & W:O({ϵ}^{-2}):{\gamma }_0{E}_0{W}_{0,\zeta \zeta }{\theta }_z^2=0,\hfill \\ \hfill & V:O({ϵ}^{-3/2}):{\gamma }_0{E}_0{W}_{0,\zeta \zeta }{\theta }_Y{\theta }_z=0.\hfill \end{array}$$

To satisfy the equations at ${ϵ}^{-2}$ and at ${ϵ}^{-3/2}$, we can choose ${\theta }_z=0$; thus, now $\theta =\theta (X,Y)$.

At $O({ϵ}^{-1}),O({ϵ}^{-1/2})$, we have:$$\begin{array}{cc}\hfill & V:O({ϵ}^{-1}):{\gamma }_0{E}_0(V_{0,\zeta \zeta }{\theta }_Y+W_{0,\eta \zeta }{\xi }_z){\theta }_Y=0\hfill \\ \hfill & W:O({ϵ}^{-1}):{\gamma }_0{E}_0(V_{0,\eta \zeta }{\theta }_Y+W_{0,\eta \eta }{\xi }_z){\xi }_z=0,\hfill \\ \hfill & U:O({ϵ}^{-1/2}):{\gamma }_0{E}_0(V_{0,\zeta \zeta }{\theta }_Y+W_{0,\eta \zeta }{\xi }_z){\theta }_X=0.\hfill \end{array}$$

From the equation at ${ϵ}^{-1}$ for V, we can conclude that ${\theta }_Y=0$. Consequently $\theta =\theta (X)$, and note that shear waves essentially depend on $\xi$, and ${\xi }_z\ne 0$ in all cases. At ${ϵ}^{-1/2}$ in the equation for U and at ${ϵ}^{-1}$ in equation for W, we obtain $W_{0,\eta }=0$, so $W_0=w_0(X,Y,z,\zeta )$. $$\begin{array}{cc}\hfill & W:O({ϵ}^{-1/2}):{\gamma }_0{E}_0(U_{0,\eta \zeta }{\theta }_X+W_{1/2,\eta \eta }{\xi }_z+V_{0,\eta \eta }{\xi }_Y){\xi }_z=0\hfill \\ \hfill & V:O({ϵ}^0):{\gamma }_0{E}_0(U_{0,\eta \zeta }{\theta }_X+W_{1/2,\eta \eta }{\xi }_z+V_{0,\eta \eta }{\xi }_Y){\xi }_Y\hfill \\ \hfill & +{\gamma }_0{\mu }_1{V}_{0,\eta \eta }{\xi }_z^2+{\omega }^2{V}_0=0.\hfill \end{array}$$

From the equation for W at ${ϵ}^{-1/2}$, we can conclude that $U_{0,\eta \zeta }{\theta }_X+W_{1/2,\eta \eta }{\xi }_z+V_{0,\eta \eta }{\xi }_Y=0$. Using this result, from the equation for V at ${ϵ}^0$, we have ${\gamma }_0{\mu }_1{V}_{0,\eta \eta }{\xi }_z^2+{\omega }^2{V}_0=0$, i.e., the Helmholtz equation for the elastic shear waves. However, in our problem, we generate by the point source only compressional waves in the liquid, without any independent shear waves in the bottom. So, we can postulate that $V_0=0$. Thus, we drop out the solutions with dominating shear waves. Now, we should require $$U_{0,\eta \zeta }{\theta }_X+W_{1/2,\eta \eta }{\xi }_z=0.$$

To resolve the equations at different powers of $ϵ$ with respect to $\zeta$, we can seek variables $U_j$ in the form $$U_j(X,Y,z,\eta ,\zeta )={\tilde{U}}_j(X,Y,z,\eta )exp(\mathrm{i}\zeta ),j=0,1/2,1,3/2,\dots .$$

This is also true for $V_j$ and $W_j$, and further tildes are omitted. Integrating the above equation by $\eta$, we obtain $U_0=-W_{1/2,\eta }{\xi }_z/(\mathrm{i}{\theta }_X)+u_0(X,Y,z)$, where $u_0(X,Y,z)$ is the integration constant. Denote $u_{S0}(X,Y,z,\eta )=-W_{1/2,\eta }{\xi }_z/(\mathrm{i}{\theta }_X)$ and integrate this equality by $\eta$. As a result, we have $$W_{1/2}=w_{S1/2}(X,Y,z,\eta )+w_{1/2}(X,Y,z),$$ where $w_{1/2}$ is also the integration constant and $w_{S1/2,\eta }=-\mathrm{i}{\theta }_X{u}_{S0}/{\xi }_z$.

At $O\mathrm{\left(}{ϵ}^0\mathrm{\right)},O\mathrm{\left(}{ϵ}^{1/2}\mathrm{\right)}$, we have:$$\begin{array}{cc}\hfill & U:O\mathrm{\left(}{ϵ}^0\mathrm{\right)}:\mathrm{i}{\theta }_X{\gamma }_0{E}_0\mathrm{(}{w}_{0,z}+\mathrm{i}{\theta }_X{U}_0\mathrm{)}+\mathrm{i}{\theta }_X{\gamma }_0{E}_0{W}_{1/2,\eta }{\xi }_z\hfill \\ \hfill & +{\gamma }_0{\mu }_1{U}_{0,\eta \eta }{\xi }_z^2+{\omega }^2{U}_0=0\hfill \\ \hfill & W:O({ϵ}^0):{\gamma }_0{[E_0(w_{0,z}+\mathrm{i}{\theta }_X{u}_0)]}_z\hfill \\ \hfill & +{\gamma }_0{E}_0(\mathrm{i}{\theta }_X{U}_{1/2,\eta }+W_{1,\eta \eta }{\xi }_z+V_{1/2,\eta \eta }{\xi }_Y){\xi }_z+{\omega }^2{w}_0=0\hfill \\ \hfill & V:O({ϵ}^{1/2}):{\gamma }_0{E}_0{(w_{0z}+\mathrm{i}{\theta }_X{u}_0)}_Y+{\gamma }_0{\mu }_1{V}_{1/2,\eta \eta }+{\omega }^2{V}_{1/2}\hfill \\ \hfill & +{\gamma }_0{E}_0(\mathrm{i}{\theta }_X{U}_{1/2,\eta }+W_{1,\eta \eta }{\xi }_z+V_{1/2,\eta \eta }{\xi }_Y){\xi }_Y=0.\hfill \end{array}$$

The consideration of the equations for U at ${ϵ}^0$ gives two equations:(A1)$$\mathrm{i}{\theta }_X{\gamma }_0{E}_0(w_{0z}+\mathrm{i}{\theta }_X{u}_0)+{\omega }^2{u}_0=0,{\gamma }_0{\mu }_1{u}_{S0\eta \eta }{\xi }_z^2+{\omega }^2{u}_{S0}=0.$$

The last equation is the Helmholtz equation for the SV elastic shear waves. However, in our problem, we can postulate $u_{S0}=0$ and also $w_{S1/2}=0$. Thus, we drop out solutions with dominating shear waves.

The consideration of the equations for W at ${ϵ}^0$ also gives two equations:(A2)$${\gamma }_0{[E_0(w_{0z}+\mathrm{i}{\theta }_X{u}_0)]}_z+{\omega }^2{w}_0=0,\mathrm{i}{\theta }_X{U}_{1/2,\eta }+W_{1,\eta \eta }{\xi }_z+V_{0,\eta \eta }{\xi }_Y=0.$$

The last equation here eliminates the last term in the equation for V at ${ϵ}^{1/2}$. Thus, from the equation for V at ${ϵ}^{1/2}$, we obtain two separate equations:$${\gamma }_0{E}_0{(w_{0z}+\mathrm{i}{\theta }_X{u}_0)}_Y+{\omega }^2{v}_{1/2}=0,{\gamma }_0{\mu }_1{v}_{S1/2,\eta \eta }{\xi }_z^2+{\omega }^2{v}_{S1/2}=0.$$

The equation for V at ${ϵ}^{1/2}$ thus (regarding the above results) leads to $V_{1/2}=v_{1/2}(X,Y,z)$$+v_{S1/2}(X,Y,z,\eta )$, where ${\omega }^2{v}_{1/2}=-{\gamma }_0{E}_0{(w_{0z}+\mathrm{i}{\theta }_X{u}_0)}_Y$, and ${\gamma }_0{\mu }_1{v}_{S1/2,\eta \eta }+{\omega }^2{v}_{S1/2}=0$. Again, the last equation is the Helmholtz equation for the SH elastic shear waves. However, in our problem, any independent shear waves in the main order are absent. So, we can postulate that $v_{S1/2}=0$. Moreover, from the second equation in (A2), we can conclude that $U_{1/2}=u_{1/2}(X,Y,z)+u_{S1/2}(X,Y,z,\eta )$ and $W_1=w_{S1}(X,Y,z,\eta )+w_1(X,Y,z)$, where $w_{S1,\eta }=-\mathrm{i}{\theta }_X{u}_{S1/2}/{\xi }_z$.

Introduce value $P_0=-E_0(w_{0,z}+\mathrm{i}{\theta }_X{u}_0)$. From the first equations in (A1) and (A2), we obtain the equation for $P_0$: $${({\gamma }_0{P}_{0,z})}_z+{\gamma }_0\left(\frac{{\omega }^2}{{\gamma }_0{E}_0}-{\theta }_X^2\right)P_0=0.$$

With the appropriate boundary and interface conditions for $P_0$, this equation consists of a typical acoustic Sturm–Liouville problem, with the spectral parameter ${\theta }_X^2$ depending on the slow variable X.

From the equations for U at ${ϵ}^{1/2}$, we have two equations:$$\mathrm{i}{\theta }_X{\gamma }_0{E}_0(w_{1/2,z}+\mathrm{i}{\theta }_X{u}_{1/2})+{\omega }^2{u}_{1/2}=0,{\gamma }_0{\mu }_1{u}_{S1/2,\eta \eta }{\xi }_z^2+{\omega }^2{u}_{S1/2}=0.$$

Additionally, from the equation for W at ${ϵ}^{1/2}$, we have:$${\gamma }_0{[E_0(w_{1/2,z}+\mathrm{i}{\theta }_X{u}_{1/2})]}_z+{\omega }^2{w}_{1/2}=0.$$

Combining the above equations and introducing value $P_{1/2}=-E_0(w_{1/2,z}+\mathrm{i}{\theta }_X{u}_{1/2})$, we obtain the equation for $P_{1/2}$:$${({\gamma }_0{P}_{1/2,z})}_z+{\gamma }_0\left(\frac{{\omega }^2}{{\gamma }_0{E}_0}-{\theta }_X^2\right)P_{1/2}=0.$$

This equation exactly coincides with the one for $P_0$. Without loss of generality, we can put $P_{1/2}=0$, and also $w_{1/2}=0$ and $u_{1/2}=0$. Because $w_{1/2}$, $u_{1/2}$, and, related with them, $v_1$ are absent in the equations for U, V at ${ϵ}^1$, we can introduce the following representation of displacements, which gives us the final parabolic equation. (A3)$$\begin{array}{cc}\hfill & u=[u_0(X,Y,z)+ϵu_1(X,Y,z)]{\mathrm{e}}^{\mathrm{i}\theta (X)/ϵ}\hfill \\ \hfill & +[{ϵ}^{1/2}{u}_{S1/2}(X,Y,z,\eta )+ϵu_{S1}(X,Y,z,\eta )]{\mathrm{e}}^{\mathrm{i}\theta (X)/ϵ}+\dots ,\hfill \\ \hfill & v=[{ϵ}^{1/2}{v}_{1/2}(X,Y,z)+{ϵ}^{3/2}{v}_{3/2}(X,Y,z)]{\mathrm{e}}^{\mathrm{i}\theta (X)/ϵ}\hfill \\ \hfill & +[ϵv_{S1}(X,Y,z,\eta )+{ϵ}^{3/2}{v}_{S3/2}(X,Y,z,\eta )]{\mathrm{e}}^{\mathrm{i}\theta (X)/ϵ}+\dots ,\hfill \\ \hfill & w=[w_0(X,Y,z)+ϵw_1(X,Y,z)]{\mathrm{e}}^{\mathrm{i}\theta (X)/ϵ}\hfill \\ \hfill & +[ϵw_{S1}(X,Y,z,\eta )+{ϵ}^{3/2}{w}_{S3/2}(X,Y,z,\eta ){\mathrm{e}}^{\mathrm{i}\theta (X)/ϵ}+\dots .\hfill \end{array}$$

Thus, the obtained ansatz contains two groups of members: describing compressional waves in the water and at the bottom and describing shear waves at the bottom.

References

1. Jensen, F.B.; Porter, M.B.; Kuperman, W.A.; Schmidt, H. Computational Ocean Acoustics; 2nd ed. Springer: Berlin/Heidelberg, Germany, 2011.

2. Manul’chev, D.; Tyshchenko, A.; Fershalov, M.; Petrov, P. Estimating Sound Exposure Levels Due to a Broadband Source over Large Areas of Shallow Sea. J. Mar. Sci. Eng.; 2022; 10, 82. [DOI: https://dx.doi.org/10.3390/jmse10010082]

3. Porter, M.B. Beam tracing for two-and three-dimensional problems in ocean acoustics. J. Acoust. Soc. Am.; 2019; 146, pp. 2016-2029. [DOI: https://dx.doi.org/10.1121/1.5125262] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31590567]

4. de Moraes Calazan, R.; Rodríguez, O.C. Simplex based three-dimensional eigenray search for underwater predictions. J. Acoust. Soc. Am.; 2018; 143, pp. 2059-2065. [DOI: https://dx.doi.org/10.1121/1.5030922] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29716240]

5. Bottero, A.; Cristini, P.; Komatitsch, D.; Brissaud, Q. Broadband transmission losses and time dispersion maps from time-domain numerical simulations in ocean acoustics. J. Acoust. Soc. Am.; 2018; 144, pp. EL222-EL228. [DOI: https://dx.doi.org/10.1121/1.5055787]

6. Lin, Y.T.; Porter, M.B.; Sturm, F.; Isakson, M.J.; Chiu, C.S. Introduction to the special issue on three-dimensional underwater acoustics. J. Acoust. Soc. Am.; 2019; 146, pp. 1855-1857. [DOI: https://dx.doi.org/10.1121/1.5126013]

7. Petrov, P.; Katsnelson, B.; Li, Z. Modeling Techniques for Underwater Acoustic Scattering and Propagation (Including 3D Effects). J. Mar. Sci. Eng.; 2022; 10, 1192. [DOI: https://dx.doi.org/10.3390/jmse10091192]

8. Collins, M.D. The adiabatic mode parabolic equation. J. Acoust. Soc. Am.; 1993; 94, pp. 2269-2278. [DOI: https://dx.doi.org/10.1121/1.407498]

9. Trofimov, M.Y. Narrow-angle parabolic equations of adiabatic single-mode propagation in horizontally inhomogeneous shallow sea. Acoust. Phys.; 1999; 45, pp. 575-580.

10. Petrov, P.S.; Ehrhardt, M.; Tyshchenko, A.G.; Petrov, P.N. Wide-angle mode parabolic equations for the modelling of horizontal refraction in underwater acoustics and their numerical solution on unbounded domains. J. Sound Vib.; 2020; 484, 115526. [DOI: https://dx.doi.org/10.1016/j.jsv.2020.115526]

11. Tyshchenko, A.; Zaikin, O.; Sorokin, M.; Petrov, P. A Program based on the Wide-Angle Mode Parabolic Equations Method for Computing Acoustic Fields in Shallow Water. Acoust. Phys.; 2021; 67, pp. 512-519. [DOI: https://dx.doi.org/10.1134/S1063771021050110]

12. Abawi, A.T.; Kuperman, W.A.; Collins, M.D. The coupled mode parabolic equation. J. Acoust. Soc. Am.; 1997; 102, pp. 233-238. [DOI: https://dx.doi.org/10.1121/1.419819]

13. Trofimov, M.Y.; Kozitskiy, S.B.; Zakharenko, A.D. A mode parabolic equation method in the case of the resonant mode interaction. Wave Motion; 2015; 58, pp. 42-52. [DOI: https://dx.doi.org/10.1016/j.wavemoti.2015.06.003]

14. Tang, J.; Petrov, P.; Kozitskiy, S.; Piao, S. On the source images method for sound propagation in a penetrable wedge: Some corrections and appendices. Proceedings of the IEEE International Conference Days on Diffraction 2017; St. Petersburg, Russia, 19–23 June 2017; pp. 304-310.

15. Tu, H.; Wang, Y.; Lan, Q.; Liu, W.; Xiao, W.; Ma, S. A Chebyshev–Tau spectral method for normal modes of underwater sound propagation with a layered marine environment. J. Sound Vib.; 2021; 492, 115784. [DOI: https://dx.doi.org/10.1016/j.jsv.2020.115784]

16. Tu, H.; Wang, Y.; Lan, Q.; Liu, W.; Xiao, W.; Ma, S. Applying a Legendre collocation method based on domain decomposition to calculate underwater sound propagation in a horizontally stratified environment. J. Sound Vib.; 2021; 511, 116364. [DOI: https://dx.doi.org/10.1016/j.jsv.2021.116364]

17. Tu, H.; Wang, Y.; Yang, C.; Wang, X.; Ma, S.; Xiao, W.; Liu, W. A novel algorithm to solve for an underwater line source sound field based on coupled modes and a spectral method. J. Comput. Phys.; 2022; 468, 111478. [DOI: https://dx.doi.org/10.1016/j.jcp.2022.111478]

18. Tu, H.; Wang, Y.; Ma, X.; Zhu, X. Applying the Chebyshev–Tau spectral method to solve the parabolic equation model of wide-angle rational approximation in ocean acoustics. J. Theory Comput. Acoust.; 2022; 30, 2150013. [DOI: https://dx.doi.org/10.1142/S2591728521500134]

19. Sabatini, R.; Cristini, P. A multi-domain collocation method for the accurate computation of normal modes in open oceanic and atmospheric waveguides. Acta Acust. United Acust.; 2019; 105, pp. 464-474. [DOI: https://dx.doi.org/10.3813/AAA.919328]

20. Luo, W.; Yu, X.; Yang, X.; Zhang, Z.; Zhang, R. A three-dimensional coupled-mode model for the acoustic field in a two-dimensional waveguide with perfectly reflecting boundaries. Chin. Phys. B; 2016; 25, 124309. [DOI: https://dx.doi.org/10.1088/1674-1056/25/12/124309]

21. Chapman, D.M.F.; Godin, O.A. Dispersion of interface waves in sediments with power-law shear speed profiles. 1. Exact and approximate analytical results. J. Acoust. Soc. Am.; 2001; 110, pp. 1890-1907. [DOI: https://dx.doi.org/10.1121/1.1401739]

22. Hudson, J.A. A parabolic approximation for surface waves. Geophys. J. Int.; 1981; 67, pp. 755-770. [DOI: https://dx.doi.org/10.1111/j.1365-246X.1981.tb06952.x]

23. Hudson, J.A. A parabolic approximation for elastic waves. Wave Motion; 1980; 2, pp. 207-214. [DOI: https://dx.doi.org/10.1016/0165-2125(80)90002-5]

24. Trofimov, M.Y.; Kozitskiy, S.B.; Zakharenko, A.D. Weak shear modulus in the acoustic mode parabolic equation. Proceedings of the IEEE International Conference Days on Diffraction 2016; St. Petersburg, Russia, 27 June–1 July 2016; pp. 416-421.

25. Nayfeh, A.H. Perturbation Methods; John Wiley and Sons: New York, NY, USA, 1973; 437.

26. Babich, V.M.; Kiselev, A.P. Elastic Waves: High Frequency Theory; CRC Press: Boca Raton, FL, USA, 2018.

27. Brehovskih, L.M.; Goncharov, V.V. The introduction to the Continuous Media Mechanics; Nauka: Moscow, Russia, 1982.

28. Collino, F. Perfectly matched absorbing layers for the paraxial equations. J. Comput. Phys.; 1997; 131, pp. 164-180. [DOI: https://dx.doi.org/10.1006/jcph.1996.5594]

29. Lin, Y.T.; Duda, T.F.; Newhall, A.E. Three-dimensional sound propagation models using the parabolic-equation approximation and the split-step Fourier method. J. Comput. Acoust.; 2013; 21, 1250018. [DOI: https://dx.doi.org/10.1142/S0218396X1250018X]

30. Sturm, F. Leading-order cross term correction of three-dimensional parabolic equation models. J. Acoust. Soc. Am.; 2016; 139, pp. 263-270. [DOI: https://dx.doi.org/10.1121/1.4939735] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26827022]

31. Petrov, P.S.; Sturm, F. An explicit analytical solution for sound propagation in a three-dimensional penetrable wedge with small apex angle. J. Acoust. Soc. Am.; 2016; 139, pp. 1343-1352. [DOI: https://dx.doi.org/10.1121/1.4944692]

32. Software “AmpeQ” for Simulating of the Acoustic Fields in the Sallow Water Waveguides by the Elastic Mode Parabolic Equation Method. Available online: https://disk.yandex.ru/d/Ao2lyJQBuNB5pQ (accessed on 14 September 2022).

33. Rutenko, A.N.; Kozitskii, S.B.; Manul’chev, D.S. Effect of a sloping bottom on sound propagation. Acoust. Phys.; 2015; 61, pp. 72-84. [DOI: https://dx.doi.org/10.1134/S1063771015010108]