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

: 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 efﬁciency 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 sufﬁciently 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.


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

Basic Equations and Expansions
We consider the propagation of time-harmonic acoustic waves in the three-dimensional where ω is the angular frequency of the waves, and the stress tensor of elasticity σ ij is (here Here, λ and µ are the Lame coefficients, ρ 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 Effective Young's modulus is E ef = ρc 2 , where c is the velocity of compressional sound waves [27]. Shear modulus is µ = ρc 2 s , 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 η = ξ(X, Y, z)/ √ 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, η), 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 Let us call P = −E 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 + . . .)e iθ(X)/ , where P 0 = −E 0 (w 0z + iθ X u 0 ). Then, where c 0 (X, z) = √ γ 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)φ(X, z), and a Sturm-Liouville (SL) problem is obtained for φ(X, z) (called acoustic spectral problem) with k 2 = (θ X ) 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 It is known that such an SL problem has a countable set of solutions (k 2 j , φ j ) indexed by j = 1, 2, . . ., where eigenfunctions φ j are orthogonal to each other and can be normalized as Each eigenfunction corresponds to a particular solution P 0j = A j (X, Y)φ 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): . Consider now the equations at 1 . The first-order approximation for the quasi-pressure has the form Additionally, the respective equations for the displacement components u 1 and w 1 write as 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.

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 = √ y and phase ζ = θ(X)/ 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 + (ω 2 /c 2 s )u = 0. However, c 2 s = µ/ρ and µ = µ 1 . So, we obtain equation u zz + (ω 2 ρ/µ 2 1 )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 η = z/ √ . In our more general case, when c s and other parameters are not constants, we should introduce the phase variable η = ξ(X, Y, z)/ √ . Thus, in our considerations, only the powers of √ 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 √ . Additionally, dependent variables should depend on X, Y, z, ζ, and η. 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.
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: ξ z w S1,η + iku S1/2 = 0 .
To solve these equations, we put u S1/2 =ũ(X, Y, z) exp(iη) and w S0 =w(X, Y, z) exp(iη): (−γ 0 µ 1 ξ 2 z + ω 2 )ũ = 0 , i(ξ zw + kũ) = 0 . Thus, we obtain the expression for phase variable ξ(X, Y, z): From the equations at 1 and 1/2 for u S and w S , we obtain: This equation with separable variables can be easily solved Finally, we have solutions for u S1/2 and w S1 : The expressions at 3/2 for u S and v S give us the equation forṽ in v S1 =ṽ(X, Y, z) exp(iη), which is identical to Equation (8) forũ, and therefore we have We also obtained the equation forũ 1 in u S1 =ũ 1 (X, Y, z) exp(iη): w 1 can be found from which also appeared when collecting terms of the order 1/2 .

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: 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 → z 0 + 0), while the subscript '−' corresponds to the limit z → z 0 − 0, and u is the particle velocity vector.
The normal vector n has coordinates while tangent vectors t 1 and t 2 can be written as having the components: We further assume that surface function h(x, y) can be represented in the form The interface conditions for the stress tensor modulo terms of the order O( 2 ) are as follows: [ 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: With all the mentioned assumptions, the interface Condition (10) at z = h 0 (X) can be rewritten in terms of P 0 and P 1 in the following form: The interface condition for the normal displacement component from (9) is Let us express this condition through P 0 and P 1 : From (7), we can express ω 2 w S0 = ikγ 0 µ 1 ξ z u S0η and, considering the first equation in (12), we obtain (ω 2 w S0 ) − = [γ 0 µ 1 k 2 (2γ 0 P 0z + γ 0z P 0 )/ω 2 ] − , and finally obtain: 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.

Boundary Conditions
At the upper boundary z = 0, we postulate σn = 0, which is reduced to 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 σn = 0, and therefore Again, we have pressure release conditions for P 0 = 0 and P 1 = 0. Moreover, σ 13 = 0 at z = −H reduces to the following equality: In a similar way, the condition σ 23 = 0 at z = −H reduces to: 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 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.

Elastic Mode Parabolic Equations (EMPE)
Let {θ j |j = M, . . . , N} be a set of phases. We introduce multi-mode representation following [13]: Which gives a solution for the compressional acoustic waves. Here, where B jl = P 1j , φ l (projections obtained via inner product). Eigenfunctions φ j with θ 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 P 1 by φ l and then integrate the resulting equation from −H to 0 by parts twice with the use of the interface Conditions (11) and (13).
The left part of this equality turns zero due to (3), and the right part represents the desired parabolic equation. The terms (k 2 j − k 2 l )B jl in these expressions can be omitted because of the resonant condition |k j − k l | = O( ) [13].
As a result, we obtain a system of parabolic wave equations for l = M, . . . , N where θ lj = i(θ l − θ j )/ , and β lj is given by the formula Here, the values ν and n 0 are the same as in the study [24], coefficients C lj have been obtained in the work [13], and β lj differs from α lj from the latter paper only by two terms containing µ 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.

Initial-Boundary Value Problem for the System of Elastic Mode Parabolic Equations
We refer to the quantities and variables X, Y, ν,γ 0 µ 1 , θ 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,ν = ν,γ = γ 0 µ = µ 1 ,ū Sj = 1/2 u S1/2j ,v Sj = v S1j ,w Sj = w S1 ,Ā j (x, y) = A j ( x, It can be easily verified that Equation (14) in physical variables are read as 2ik lĀl,x + ik l,xĀj +Ā l,yy +β llĀl + whereβ lj are computed by the same formulae as β lj , with ν replaced byν,θ lj = (θ l −θ j ). For Equation (15), we consider the initial-boundary value problem in the domain of the form {(x, y)|0 ≤ x < ∞ , −y 0 /2 ≤ y ≤ y 0 /2}, with the initial condition 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.
We have expressions for the displacements for the considered case: Final expressions for the displacements in the shear waves in the physical variables arē z ω c s dz .

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 ≈ 2.86 • [1,14]. Moreover, in some calculations, we used a wedge with an angle ≈ 1.14 • . 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 3 , 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.   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 ≈ 2.86 • (left) and for the wedge angle ≈ 1.14 • (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.

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 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 × 1000 nodes. To prevent reflections from the boundaries, we used transparent boundary conditions at y = ±1500 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.

•
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 ≈ 300 m/s at the bottom and a rather good agreement of up to c s ≈ 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.

Appendix A. Derivation of the Main Ansatz
Here, we derive ansatz for our problem from the general form of displacement component expansion.
where U = (U, V, W) is the vector of displacements, η = ξ(X, Y, z)/ √ , and ζ = θ(X, Y, z)/ . Now, we expand the equations of elasticity and consider terms at the different powers of .
From the equation for W at −1/2 , we can conclude that U 0,ηζ θ X + W 1/2,ηη ξ z + V 0,ηη ξ Y = 0. Using this result, from the equation for V at 0 , we have γ 0 µ 1 V 0,ηη ξ 2 z + ω 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,ηζ θ X + W 1/2,ηη ξ z = 0 .
At O( 0 ), O( 1/2 ), we have: The consideration of the equations for U at 0 gives two equations: 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.
Introduce value P 0 = −E 0 (w 0,z + iθ X u 0 ). From the first equations in (A1) and (A2), we obtain the equation for P 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 θ 2 X depending on the slow variable X.
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.