Improving the Performance of Mode-Based Sound Propagation Models by Using Perturbation Formulae for Eigenvalues and Eigenfunctions

: Numerous sound propagation models in underwater acoustics are based on the representation of a sound ﬁeld in the form of a decomposition over normal modes. In the framework of such models, the calculation of the ﬁeld in a range-dependent waveguide (as well as in the case of 3D problems) requires the computation of normal modes for every point within the area of interest (that is, for each pair of horizontal coordinates x,y). This procedure is often responsible for the lion’s share of total computational cost of the ﬁeld simulation. In this study, we present formulae for perturbation of eigenvalues and eigenfunctions of normal modes under the water depth variations in a shallow-water waveguide. These formulae can reduce the total number of mode computation instances required for a ﬁeld calculation by a factor of 5–10. We also discuss how these formulae can be used in a combination with a wide-angle mode parabolic equation. The accuracy of such combined model is validated in a series of numerical examples.


Introduction
Computational technique based on the normal modes theory is widely used in underwater acoustics and its applications that cover a large area of marine sciences [1][2][3]. It is important that normal modes provide both quantitative and qualitative understanding of many physical effects related to sound propagation in the sea. Investigation of the dependence of modal wavenumbers k j and eigenfunctions φ j (z) on certain environment factors and parameters is often sufficient for estimating the influence of these factors on the structure of acoustical field (in particular, on interference patterns, shadow zones, arrival times, and energy distribution [4,5]).
Indeed, sound field in a 3D oceanic waveguide formed by the sea surface z = 0 and sea bottom z = h(x, y) (here z stands for depth, while x, y are the horizontal coordinates, see in [1] for detailed discussion of waveguides in underwater acoustics) can be represented in the form of the expansion over local normal modes φ j (z, x, y) (i.e., the eigenfunctions φ j (z) computed for the given values of x, y). Under the adiabatic assumption, mode amplitudes A j (x, y) satisfy the so-called horizontal refraction equation (HRE) [1,5,6] and therefore the dependence of eigenvalue k 2 j = (k j (x, y)) 2 on the variables x, y is what actually determines the propagation of sound waves associated with each mode in the horizontal plane.
In this case, the distribution of k j (x, y) plays the role of effective refractive index for horizontal rays [5,7]. In many similar sound propagation models (e.g., in the mode parabolic equation theory [8][9][10]), variation of media parameters in the horizontal plane also influences acoustical field via the corresponding variability of horizontal wavenumbers in Equation (2) and eigenfunctions in Equation (1).
It is therefore desirable to have explicit formulae describing the dependence of eigenvalues k j and modal functions φ j (z) on the horizontal variables x, y.
The variability of modes in x and y can result from the sound speed dependence on x, y. In this case, the standard perturbation theory from quantum mechanics [1,11] can be applied to take this dependence into account. For example, it is routinely used to the compute imaginary parts of eigenvalues and modal group velocities [1].
However, eigenvalues and mode functions can also depend on the horizontal coordinates due to variations of bottom relief. In shallow-water acoustics, this factor is of primary importance, especially for low-frequency sound that is less sensitive to volume inhomogeneities of the sound speed in the water column. This problem is obviously connected with differentiation of eigenvalues and mode functions with respect to the water depth h. Explicit formulae for the derivatives of eigenvalues with respect to range in the presence of water depth variations h = h(x, y) were first derived by Brekhovskikh and Godin [12]. Later, Godin [13] also obtained expressions for the derivatives of mode functions from the so-called generalized orthogonality relationships [12]. Obviously, the problem of differentiation of eigenvalues and eigenfunctions with respect to the range r can be reduced to the differentiation with respect to water depth h, provided that environmental gradient, i.e., dh/dr, is known. Identical formulae were derived independently by Trofimov [14] who used multi-scale approach. To our knowledge, however, neither a second-order perturbation theory nor explicit formulae for second-order derivatives of eigenvalues and mode functions in the case of water depth variations h = h(x, y) = h 0 + ∆h(x, y) (see Figure 1) have been presented in the literature. It appears that the derivation techniques from in [12][13][14] cannot be extended to the second-order case. We also observed that the convergence of Godin's series for the derivatives of mode functions is slow, and the respective approximation error decreases in a non-monotonic way. Clearly, this makes the formulae for the derivatives of modal functions from in [13,14] ill-suited for practical computations. Very recently, explicit formulae for the first and second derivatives of eigenvalues and modal functions in the case of water depth variations were outlined by Petrov et al. [15]. The approach used in the latter paper is different from those of Godin [13] and Trofimov [14]. For the first-order derivatives of eigenvalues, the resulting formulae are the same as derived by Trofimov and Godin [13,14] (modulo some simple transformations). At the same time, the numerical examples presented in this study indicate that our expressions for the derivatives of mode functions appear to be more practical and computationally robust.
The main goal of this study is to provide a detailed derivation of the formulae from in [15], including the generalization to the third-order perturbations of horizontal wavenumbers that were not given in the mentioned paper, and to show how this perturbation theory can be used in realistic problems of sound propagation. Note that in [15] only the very idea of the derivation of the perturbative formulae is given, while the calculations are largely omitted.
The paper is organized as follows. We start with the formulation of the Sturm-Liouville problem from which horizontal wavenumbers and mode functions are determined (Section 2). In Section 3, we introduce change of variables that reduces it to a form for which the standard results of the perturbation theory for linear operators [16] can be applied. Next, we derive perturbation formulae for the solutions of the latter problem. In Section 4, we use the derived perturbative formulae for the wavenumbers and mode functions to compute acoustic field by the formula (1), where the mode amplitudes A j are computed using a wide-angle mode parabolic approximation for Equation (2). We show that even a single call of the spectral problem solver together with our perturbation theory allows to perform accurate simulation of sound propagation in a very large area.

Acoustical Spectral Problem
The main goal of the present study is to improve the computational efficiency of sound propagation models based on the normal modes theory that rely on representation of acoustic field in the form (1). In the course of the field computation, such models usually make numerous calls of a solver of the following acoustic spectral problem [1,2]: on an interval z ∈ [0, H], where F| z=h ± denotes one-sided limits of some function F(z) (possibly discontinuous) as z approaches the point z = h from below or above, respectively. The function c(z) in (3) is assumed to be continuously differentiable on the intervals I 1 = [0, h) and I 2 = (h, H] (it can therefore have a finite jump discontinuity at z = h), and the density is a piecewise-constant function: The second and the third equalities in Equation (3) express pressure-release boundary conditions at the surface z = 0 and at some sub-bottom z = H (the sub-bottom is introduced to avoid complications associated with continuous spectrum of the halfspace). The value of H is chosen to be sufficiently large in order to ensure the convergence of the solution of the propagation problem of interest. The third and the fourth equalities are continuity conditions for sound pressure and particle velocity at the interface z = h between the water and the upper sediment layer of the bottom.
Solutions of the problem (3) are pairs (k j , φ j (z)), where k j is called horizontal (modal) wavenumber, and φ j (z) is the respective (vertical) mode function. The Dirichlet boundary conditions at the endpoints of the interval [0, H] in Equation (3) ensure that the problem has purely discrete spectrum consisting of a countable set of real eigenvalues k 2 j , j = 1, 2, . . . . We assume that they are ordered in such a way that k 2 j ≥ k 2 j+1 for all j. Note that the mode functions form an orthonormal basis with respect to the scalar product i.e., (φ i , φ j ) = δ ij (see [1,2]). Hereafter, we are mainly interested in waterborne modes, i.e., such modes that their eigenfunctions φ j (z) have maxima in the water column, and their eigenvalues k j belong to the interval [ω/c b , ω/c min ] (c b is the sound speed in the upper layer of the bottom, and c min is the minimum value of the sound speed in the water column). Away from the source the field is mostly formed by waterborne modes, and therefore in practical problems of underwater acoustics they are of primary importance.
In the literature, horizontal wavenumbers k j are often called eigenvalues of the problem (3) [1], although strictly speaking this term should be used only for their squares k 2 j . According to this tradition, hereafter we also call k j eigenvalues.
Note that the problem (3) can be considered a regular counterpart to the singular spectral problem for the Pekeris operator [17,18] (where the condition at z = H is replaced by the boundedness condition at z → ∞). Within this work, we restrict our attention to the self-adjoint case, i.e., to the case of lossless bottom (with negligible attenuation). The regular case considered below is arguably more important from the practical point of view.
In this case, k j and φ j can be considered functions of h (while the latter is in turn a function of horizontal coordinates x, y). An illustration of this dependence for a typical shallow-water waveguide with penetrable bottom is presented in Figure 1.
Observe that for the field computation by the formula (1) it is necessary to know k j (x, y) and φ j (z, x, y) for all values of x, y (or for all values of h(x, y) along in the area of interest). This requires the solution of the problem (3) for all values of h with certain sufficiently small step. On the other hand, we can consider the variations of water depth as a perturbation with respect to certain average value h 0 . Expressing the bathymetry variations h = h 0 + ∆h(x, y) as fluctuations of water depth ∆h(x, y) around some average value h 0 we can use Taylor series expansion for horizontal wavenumbers and mode functions where prime denotes the derivative of k j or φ j (z) (for a fixed value of z = z r ) with respect to the parameter h.

Perturbation of Normal Modes by the Bathymetry Variations
The problem for the bathymetry perturbation (i.e., interface perturbation) can be reduced to a well-studied problem of the potential perturbation in the stationary Schrödinger equation [11]. Let us introduce new variable t = zh 0 /h and denote T = Hh 0 /h. After substitution into Equation (3), we obtain the following spectral problem: where k j are the same wavenumbers as those obtained by solving (3), and eigenfunctions ψ j (t) are related to the eigenfunctions of the original problem (3) by the relation where Q = ω 2 c 2 , and superscripts in parentheses stand for derivatives with respect to . Now, we follow the standard scheme of the perturbation theory. Substituting expansions (9) into (8) and separating terms with identical powers of epsilon , we obtain for 0 the spectral problem (8) for k In order to obtain k where dots denote the terms containing ψ j and their derivatives that vanish according to the following identities: Separating in (12) terms of different powers of , we obtain a family of boundary j that read as Now, let us compute the scalar product of Equation (11) with the function ψ using the following standard orthogonality identities (see, e.g., in [11]) Integrating the first term on the right-hand side of Equation (15) by parts and using (16) and (14), we obtain the following expression for k where we introduced the following notation: for the terms associated with the sub-bottom stretching in new coordinate t (obviously, these terms can be neglected for trapped/waterborne modes).
Multiplying Equation (11) by the function ψ (0) i and repeating the steps described above, we arrive at the following formula for the scalar product of ψ (1) j and ψ (0) i : As {ψ i } form an orthogonal basis, the function ψ (1) j can be expressed in the form of the following series: Functions ψ (1) by their construction are the derivatives ∂ψ ∂ = h 0 ∂ψ ∂h . Consequently, the derivatives ∂φ ∂h of eigenfunctions φ(z) of the original problem (3) can be found as Let us now proceed with the terms of the order 2 that arise when substituting (9) into (8). We obtain the following equation: Multiplying it by ψ (0) j (in the sense of the scalar product), and using the orthogonality relations (16) and boundary conditions (12), we obtain an expression for k The scalar product of the equality (23) and ψ (0) i leads to the expression for ψ while the second derivative of the eigenfunction of the original problem (3) can be found as Let us conclude this section with the formula for k (3) j , that can be obtained in by reiterating the calculations above:

Numerical Example: the 3D Coastal Wedge Problem
Consider a coastal wedge shown in Figure 2 with the opening angle α = 2.86 • and water depth near the source h = 200 m. Assume that the sound speed and the density in the water column are c 1 = 1500 m/s and ρ 1 = 1 g/cm 3 , respectively, and the values of these parameters for the bottom are c 2 = 1700 m/s and ρ 2 = 1.5 g/cm 3 . Sound propagation in the wedge with the parameters specified above is considered a standard benchmark problem for 3D models of sound propagation in underwater acoustics (see, e.g., in [1,20]). A Cartesian coordinate system in the waveguide is chosen in such a way that x axis is aligned along the isobath, while the direction of the y axis coincides with the water depth gradient, and a point source of the frequency f = 25 Hz, is deployed at the depth z s = 100 m (its horizontal coordinates are x s = 0, y s = 0). The problem of sound propagation in the coastal wedge with the parameters specified above is usually used for the validation of various mathematical models in shallow-water acoustics. It was recently shown that an accurate solution of this problem in the cross-slope direction can be obtained under the adiabaticity assumption [21]. In our opinion, this result is remarkable, as it indicates that many propagation problems can be successfully tackled by a computationally lightweight adiabatic models.
The main idea of the present study is that the efficiency of the respective approaches can be further improved if we reduce the number of calls of the Sturm-Liouville problems solver. The first thing to do in this direction is to answer the question on how far can we go from the source with just one call of the normal modes computation. Arguably, the wedge problem can be considered a worst case scenario for the application of our perturbation theory, as depth variations in this environment are by no means small. Indeed, an accurate solution of the wedge problem reported in [21] required the computation of modes for the values of water depth in the range 50-350 m, i.e., and therefore we have h 0 = 200 m, and ∆h = −150 · · · + 150 m. Figure 3 shows the dependencies of the wavenumbers k j on the water depth h for all three waterborne modes excited by the source. The red lines correspond to the "exact" eigenvalue dependencies on h. Approximations of these dependencies obtained by the perturbation formulae of the first, second, and third order are represented by magenta, blue, and black lines, respectively. It can be seen that the 3rd-order approximation is almost perfect for ∆h ∈ [−100 m, 100 m] for the first and second modes, and also works pretty well for ∆h ∈ [−50 m, 100 m] for the third mode. The 2nd-order perturbation theory provides accurate approximations of the wavenumber variations for |∆h| ≤ 50 m, while the linear (1st-order) approximation of k j (h) ensures reasonable quality only for |∆h| ≤ 10-20 m. Note that all approximations are clearly invalid beyond the cut-off depth value. Yet, highorder approximations offer considerable improvement over the accuracy of the linear one.
It is also important that the calculations of the derivatives k Although the accuracy of approximation of the eigenvalues k j dependence on h can be of interest per se, it is more important to investigate the accuracy of the field computation if our formulae are used to compute φ j (z, x, y) in the expansion (1) and k j (x, y) in Equation (2). In this study, we compute the solution of the wedge problem using the perturbation theory for the normal modes in combination with the field representation (1) and the pseudodifferential parabolic equations for mode amplitudes A j (PDMPEs). The latter equations can be used to compute the solution of Equation (2) neglecting the back-scattering (which is known to be small in the case of the penetrable wedge), see in [21,22] for the details.
The solution of the wedge problem obtained by the scheme outlined here is obtained under several assumptions and approximation. In order to validate its accuracy, we compare it against the solution computed using source images technique [20]. The comparison along the x axis for y = 0 and z = z r = 30 m is presented in Figure 4. From Figure 4d, it is clear that the solution based on the PDMPE theory and adiabatic approximation (we actually use a wide-angle parabolic approximation for Equation (2) and the expansion (1) according to the theory from in [21,22]). exhibits excellent agreement with the reference solution if the spectral problem is solver for each value of h (for all x). If we solve it only for h = h 0 and use first-order perturbative formulae, then the field can be computed accurately only for x < 6-7 km (see Figure 4a). A significant improvement can be achieved by employing 2nd-order formula (Figure 4b), and the agreement with the reference solution in this case is very good up to x = 15 km, which is already sufficient for many practical problems arising in underwater acoustics, where the speed of the field simulation is of primary importance. Note that in this case the 3rd-order formula extends the correct solution by a couple of kilometers in range, as shown in Figure 4c. As for h < 100 m (i.e., for |∆h| > 100 m), even the third-order approximation for the wavenumbers k j (h) fails to reproduce actual dependence (see Figure 3), the solution for x > 17 km cannot be accurate. Indeed, it is known that interference pattern along this interval is formed by the field of the first mode that reaches this area by "direct" and "refracted-by-the-wedge" paths [23,24], and the turning point of the latter is located very close to the wedge apex, where perturbation theory does not provide sufficient accuracy for the wavenumbers.   The slice of acoustical field P(x, y, z) by the horizontal plane z = z r = 30 m for the normal modes computed by the 3rd-order perturbation theory and the respective reference solution are shown in Figure 5. It can be seen that within the rectangle |y| ≤ 2 km, x ≤ 15 km that the solution computed using k j and φ j approximated by the formulae (6), (7) cannot be distinguished from the reference solution. At the same time, the perturbation theory allows to reduce the computational cost by a factor of 10 in this case.
The simulation results confirm that the perturbation theory for normal modes developed in our study can be successfully used to improve the efficiency of the sound propagation model based representation of acoustical field in the form (1), and that sometimes it is sufficient to solve the spectral problem (3) only once for some average water depth value h 0 .
We also observed that even if the accuracy provided by the perturbation theory and a single call of mode solver is not sufficient, one can use solve the spectral problem for maximal and minimal values of the depth, compute the derivatives of wavenumbers and mode functions by our formulae, and then use a clamped spline approximation for the entire area of interest. For example, in the wedge case this results in a solution that perfectly coincides with the one obtained by using k j computed for all values of the water depth. This two-point clamped spline approach can be considered a good compromise between the efficiency and the accuracy.

Conclusions
In this study, formulae for the first-and second-order derivatives of wavenumbers and eigenfunctions of normal modes in a shallow-water waveguide with respect to the water depth are obtained (we stress again that main formulae were outlined in [15], but the detailed derivation was not presented there). They are derived by using vertical coordinate transformation that leads to the reformulation of the problem in such a way that the interface perturbation turns into the potential perturbation in a stationary Schrödinger equation. Clearly, our results can be generalized to the interface perturbation theory of arbitrary order. On the practical side, however, most problems can be covered by the second-and third-order formulae presented here.
Note that first-order derivatives of wavenumbers and mode functions with respect to h were computed by a different method in [9,12,13], however unfortunately the latter approach cannot be generalized to obtain higher-order formulae. It can be shown that our expression for the first derivative (20) can be reduced to the respective formula from [9]. However, the two expressions are absolutely different from the computational point of view. Indeed, in our case the main contribution to the first derivative of the eigenfunction with respect to the h parameter is made by the derivatives of the latter with respect to z, and the series over other unperturbed eigenfunctions can be considered a small correction to the latter. By contrast, the expression for the eigenfunction perturbation derived by Trofimov is a pure expansion over unperturbed eigenfunctions. As such, it is much less robust, as the convergence of the series is quite slow.
In this study, we presented an example where the acoustical field is computed in a 3D wedge benchmark problem. Even for this idealized scenario the computational time can be reduced by a factor of 10 using our perturbative formulae. In more realistic problems, the reduction in the computational cost can be even more significant, especially when it is necessary to take mode coupling effects into account.
Author Contributions: Methodology, P.P. and M.T.; formal analysis, A.Z. and P.P.; investigation, A.Z. and P.P. All authors have read and agreed to the published version of the manuscript.
Funding: This study was supported by POI FEB RAS Program "Modeling of various-scale dynamical processes in the ocean" (project No. 121021700341-2).