Generalized Form of the Invariant Imbedding Method and Its Application to the Study of Back-Scattering in Shallow-Water Acoustics

: A generalized form of the matrix-invariant imbedding method was developed to solve boundary-value problems for coupled systems of Helmholtz-type equations. Within this approach, a boundary-value problem solution can be obtained by solving evolutionary ﬁrst-order imbedding equations for a matrix-valued function. The proposed method is applied to the solution of coupled equations for mode amplitudes describing the propagation of acoustic waves in a range-dependent shallow-water waveguide. The back-scattering of modes by bathymetry features is investigated, and the coefﬁcients of the modal expansion of the wave reﬂected by an inhomogeneity in the bottom relief are computed. It is demonstrated that back-scattering is strongly connected with the modal interactions and that the back-scattered ﬁeld consists of modes with numbers different from the number of the incident mode.


Introduction
The solution of sound propagation problems in practical applications of underwater acoustics often results in very intense computations, since areas of interest often extend hundreds or even thousands of kilometers along horizontal directions. In order to handle such large-scale problems, one usually has to neglect certain physical effects (which are considered to be of secondary importance) when simplifying propagation equations. For example, in all models based on the parabolic equations, the back-scattering of sound is not taken into account for the sake of the simplification of the problem and in exchange for the possibility of using marching computational schemes with large step sizes in range [1][2][3]. Another common approach is to neglect the mode-coupling effect when using computational models based on the theory of normal modes [1,[4][5][6]. This simplification allows one to uncouple the equations for mode amplitudes, thus greatly simplifying their solution (while retaining their capability of handling many practical problems).
When performing the aforementioned simplifications, it is always necessary to investigate the inaccuracies in the resulting models of sound propagation that arise due to the neglect of certain terms/effects. Such an investigation can be carried out by using numerically exact sound propagation techniques (e.g., [7,8]). Some forty years ago, a method that allows a reformulation of a boundary-value problem (BVP) for the Helmholtz equation (in particular, related to wave propagation and scattering) in the form of an initial-value problem for the first-order Riccati equation was proposed by Babkin and Klyatskin [9]. The latter equation can be solved with efficient marching numerical algorithms, and yet the reformulation allows one to fully take back-scattering effects into account. This technique is known in the literature as the invariant imbedding method.
In this study, we propose a generalization of this method so that it can be applied to a system of coupled Helmholtz equations (or, in the same way, to a Helmholtz equation for an unknown vector function). Within the proposed approach, the BVP for a system of coupled Helmholtz equations is replaced by a matrix Riccati equation, where the unknown quantity is the matrix of reflection coefficients as a function of the right boundary of an inhomogeneous layer. The generalized invariant imbedding method is applied to the investigation of the back-scattering of sound on a bottom-relief inhomogeneity in a shallow-water waveguide (see Figure 1). In the framework of the normal mode theory, the solution of this problem consists in the computation of mode amplitudes satisfying a coupled system of 1D elliptic equations [1]. The proposed method allows us to investigate the interplay between the back-scattering and the mode-coupling effects in a shallowwater environment. The computation of the matrix of reflection coefficients for different incident modes reveals the fact that the back-scattering is always accompanied by strong mode coupling. It is interesting to note that matrix Riccati equations for reflection coefficients can also be found in the literature in the context of elastic wave propagation [10][11][12]. Although their derivation in the mentioned works follows a different scheme, it appears that such equations can be applied to the simulation of any sort of waves described by a vectorvalued function.

Problem Statement
The problem of sound propagation in the 2D waveguide of a shallow sea Ω = {(x, z) | z ≥ 0} that is depicted in Figure 1 is considered. The waveguide consists of a water layer 0 ≤ z ≤ h(x) and a liquid-penetrable bottom. These layers are separated by the interface z = h(x) (the function h(x) is assumed to be continuous). It is supposed that the water depth h outside the interval x ∈ [L 0 , L] is constant (h = h 0 ). Thus, the bottom-relief inhomogeneity is localized on this interval. Sound propagation in the waveguide Ω is described by the 2D Helmholtz equation: where P(x, z) is the acoustic pressure, c is the sound speed, and ω is the angular frequency. It is assumed that a normal mode exp(−ik 0,j x)φ 0,j (z) is incident at the bottom-relief inhomogeneity from the right. Here, φ 0,j is an eigenfunction of the acoustic spectral problem [1]: and k 0,j is the corresponding horizontal wavenumber, ρ denotes the density (note that there is usually a density jump at the water-bottom interface), and H max is a sufficiently large value of depth at which the computational domain is truncated (so that this value does not actually influence the simulation results). Throughout this paper, we assume that the mode functions are normalized, i.e., that The solution of Equation (1) in the waveguide can be represented [1] in the form of a truncated modal expansion (hereinafter, only the first n modes are taken into account) where φ j (x, z) are eigenfunctions obtained from (2) for a fixed value of x, and A j (x) are mode amplitudes satisfying the following coupled system of equations: where ρ dz are elements of the square n × n matrices U and V , respectively.
The set of n scalar equations (3) can be replaced by a single equation for the unknown vector function a(x) = (A 1 (x), A 2 (x), . . . , A n (x)) T : where is a diagonal matrix with eigenvalues k j (x) on the main diagonal).
In order to formulate a BVP for Equation (4), we derive the boundary conditions at x = L 0 and x = L in the case of incidence of a single normal mode at the right boundary of the inhomogeneous section of the waveguide Ω. In what follows, we write a l (x; L) instead of a(x), emphasizing the dependence of the solution on the position of the right boundary (in accordance with the notation used in [9]). The subscript l corresponds to the number of the incident mode. In addition, it is supposed that time dependence is suppressed in such a way that e ikx corresponds to propagation in the positive direction of the X axis [1].
Since the field for x ≥ L is a superposition of the l-th incident normal mode and back-scattered wave consisting of all modes taken into account, the vector a(x; L) can be represented in the following form: . . .
where R p,q are reflection coefficients to be computed.
For x ≤ L 0 , the vector a(x; L) can be written in the form where T p,q are transmission coefficients. It is easy to see that a(x; L) given by (5) and (6) satisfies the following boundary conditions at x = L and x = L 0 : where K 0 ≡ K(L) = K(L 0 ), and the elements of the vector on the right-hand side of the second condition are all zeros, except for the l-th one. The boundary conditions (7) and the vector Helmholtz Equation (4) can be considered as a vector boundary-value problem. All n column vectors a l (x; L) can be combined into a square matrix A(x; L) that satisfies the matrix boundary-value problem where B = ∂A(x;L) ∂x , D 2 = −D 0 , and D 0 = iK 0 (i.e., both D 0 and D 2 are diagonal matrices). In order to simplify the derivation, the solution of the boundary-value problem with the transformed right part of the second boundary condition will be sought.
where E is in the n × n identity matrix. The solution of the original boundary-value problem (8) can be found using the formula Hereinafter, the tilde symbols are suppressed. In order to obtain the imbedding equations for the BVP (9), we differentiate it with respect to the parameter L (following the same derivation scheme as that used by Babkin and Klyatskin [9] in the scalar case).
Using the system (9) and the second boundary condition, we obtain where, by definition, A L ≡ A(L, L) (similarly, B L ≡ B(L, L)). If we compare the last BVP with the problem (9), we can see that their solutions are connected by the following relation: which can be considered as an ordinary differential equation with independent variable L. Its initial condition can be written as The last thing that we need is the equation for the function A(L; L) ≡ A L . Let us differentiate this matrix-valued function: The last equation can be transformed into a new one by using (12) and (13).
The equation can be considered as a matrix Riccati equation. The initial condition for the equation can be obtained from the boundary conditions (12).
From (5) and (10), a formula for the reflection coefficients can be derived: Note that for n = 1 (that is, when the mode coupling is neglected), the derived imbedding Equations (13) and (15) reduce to their scalar counterparts proposed by Babkin and Klyatskin [9].
Both imbedding Equations (13) and (15) can be easily solved, e.g., with a Runge-Kutta method (the fact that the unknown function is matrix-valued does not lead to any additional complications). Matrices U and V can be easily computed by using formulae for the derivatives of the eigenfunctions with respect to the media parameters [13,14]. The formulae for the first-order derivatives from the earlier work of Godin [13] are more general and allow one to handle multiple interfaces, as well as range-dependent sound speed and density profiles. The main advantage of the approach from the recent study [14] is the possibility to obtain high-order derivatives with respect to h (the formulae for the first-order derivatives from [14] are equivalent to their counterparts from [13]).

Modal Decomposition of the Back-Scattering Field: Numerical Results
An example of a waveguide with the following parameters is considered. The water depth h 0 outside the inhomogeneity is equal to 90 m. The sound speed in the water is c 1 = 1500 m/s, and the sound speed in the bottom is c 2 = 2000 m/s. It is assumed that the density of the water is ρ 1 = 1 g/cm 3 and the density of the bottom is ρ 2 = 2 g/cm 3 . The frequency of the source is f = 50 Hz. The lower boundary of the computational domain was set to H max = 200 m. For the specified problem parameters, the source excites four waterborne modes for h = h 0 . We first take into account 20 modes excited by the source, and this number is sufficient to ensure the convergence of the obtained solution.
The depth change associated with the inhomogeneity is varied from 1 to 10 m. The profile of the bottom-relief inhomogeneity is described by the formula where L 0 = 4600 m and L = 5200 m. This is shown in Figure 2.    It is easy to see that the values of non-adiabatic reflection coefficient are much larger than the respective quantities obtained in the adiabatic approximation. In addition, the incident modes excite modes with different numbers in the reflected field (the reflected amplitude of the incident mode is relatively small). Thus, we conclude that the back-scattering of acoustic waves in range-dependent environments cannot be considered separately from the mode-coupling effect.
Note that, in fact, all coefficients have relatively small values, since the back-scattering in shallow-water environments with smooth bathymetry variations is not very pronounced. In Figure 3b (which corresponds to the incidence of the fourth mode), one can see that the relations between the values of the reflection coefficients are similar.
The second example of a waveguide is shown in Figure 4. The waveguide has the same parameters, but the depth of the inhomogeneity is varied from −1 to −10 m. Basically, this is a small underwater mountain.

Conclusions
In this study, we generalized the invariant imbedding method to the case of a BVP for a vector function. In this case, the derivation results in imbedding equations for a matrixvalued function. In particular, we obtained a matrix Riccati equation for the reflection coefficients. The solution of this equation can be easily obtained with standard ODE integration tools, e.g., with numerical Runge-Kutta schemes. Note that, in the present work, we restricted our attention to the computation of the reflection coefficients from the solution of Equation (15), whereas the computation of the mode amplitudes also requires the solution of Equation (13) (together, they are called imbedding equations by the analogy with the respective scalar equations from the pioneering work [9]). Note that the former is a nonlinear first-order ODE for a matrix-valued function, while the latter equation is linear and, therefore, much easier to solve. We do not present its solution here, as the magnitude of the back-scattered field would anyway be too small to see against the forward-propagating wave (it would be interesting, however, to compute the field, e.g., in 3D environments with rotational symmetry in the horizontal plane by using our approach).
The generalized invariant imbedding method proposed here was applied to the investigation of the modal content of a back-scattered acoustic field in a waveguide with a bottom-relief inhomogeneity. More precisely, we obtained the imbedding equations describing coupled mode amplitudes and computed the modal decomposition coefficients for the reflected field in the case in which the incident wave consisted of a single normal mode. Our results indicate that back-scattering is inevitably accompanied by strong mode coupling, and the reflected field is mostly formed by modes that are different from the incident one. In particular, the latter fact implies that one can always opt for a one-way propagation model once the adiabatic approximation is invoked.
Although in the idealized example considered here, back-scattering occurs only due to the variations in the water depth, the proposed technique can be also used in the case of more realistic environments with a lossy bottom, multiple interfaces (i.e., several layers of the bottom with finite jumps in sound speed and density at the respective boundaries), and range-dependent sound speed profiles. Such cases will be studied in future work. It would be also interesting to use the generalized invariant imbedding equations to solve 3D sound propagation problems in shallow-water waveguides.
Author Contributions: Methodology, P.P. and K.K.; formal analysis, M.K. and P.P.; investigation, M.K. and P.P. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.