Random Matrix Theory for Sound Propagation in a Shallow-Water Acoustic Waveguide with Sea Bottom Roughness

: The problem of sound propagation in a shallow sea with a rough sea bottom is considered. A random matrix approach for studying sound scattering by the water–bottom interface inhomogeneities is developed. This approach is based on the construction of a statistical ensemble of the propagator matrices that describe the evolution of the waveﬁeld in the basis of normal modes. A formula for the coupling term corresponding to inter-mode transitions due to scattering by the sea bottom is derived. The Weisskopf–Wigner approximation is utilized for the coupling between waterborne and sediment modes. A model of a waveguide with the bottom roughness described by the stochastic Ornstein–Uhlenbeck process is considered as an example. Range dependencies of mode energies, modal cross coherences and scintillation indices are computed using Monte Carlo simulations. It is shown that decreasing the roughness correlation length enhances mode coupling and facilitates sound scattering.


Introduction
Sound scattering by random inhomogeneities of the marine environment is one of the fundamental problems in ocean acoustics.Usually, multiple scattering of acoustic signals causes phase and amplitude fluctuations and complicates interference pattern and dispersion characteristics [1].This process acts as an obstacle preventing the application of acoustic methods to solving various practical problems in marine science, e.g., underwater communication and navigation [2][3][4][5], or hydroacoustical tomography [6][7][8][9][10].
The stochastic nature of inhomogeneities anticipates the usage of statistical methods for their description.Typically, longitudinal inhomogeneities of the marine environment are relatively weak and can be fairly described as perturbations superimposed onto a range-independent waveguide, in which the acoustic field is conveniently described in terms of normal modes.This explains why mode-based approaches taking the effect of scattering into account via the mode coupling are widely used in underwater acoustics community [11].
One of the most popular approaches is the derivation of coupled kinetic equations for ensemble-averaged modal intensities.It was originally proposed in the seminal work of Tappert of Dozier [12], and further developed in a series of papers [13][14][15][16].The ray-based theory of mode coupling was developed by Virovlyansky [17].In the present paper, we consider another promising approach based on the construction of a wavefield propagator based on the random matrix theory [18,19].By definition, a propagator is an operator that governs wavefield evolution in the course of its propagation along the waveguide [20].The matrix representation of the propagator can be obtained using the basis of normal modes, and the matrix elements are complex-valued random amplitudes of intermode transitions (the propagator matrices considered here are not to be confused with the ones emerging in the vector-matrix formalism [21] for the representation of the solutions for wave equations).The main advantage of this approach is the possibility to take into account interference effects within a full-wave picture.We conducted early studies of a wavefield propagator for idealized range-periodic waveguides in the context of the wave chaos problem [22][23][24].The propagator properties for waveguides with random inhomogeneities caused by linear internal waves were considered in [25][26][27].The case of nonlinear internal waves was considered in [28][29][30].In [31], this approach was generalized onto the case of waveguides with adiabatic longitudinal variations.The propagator taking into account temporal variability of the medium was considered in [5].
In the majority of works concerning the wavefield propagators for randomly inhomogeneous waveguides, deep-water propagation scenarios with volumetric inhomogeneities were investigated.At the same time, the propagators describing sound scattering on random sea bottom roughnesshave not been thoroughly studied until now.Indeed, scattering from rough surfaces is a fairly extensive field of wave physics (see, for example, a brief review [32], or another paper from the present special issue [33]).However, multiple scattering from roughness in waveguides has some specific features and still is a predominantly open problem.Sound scattering of this kind has a lot in common with the general problem of wave transport in disordered quantum mesoscopic waveguides [34,35], and some of the ideas of the matrix propagator theory are drawn from there (this issue is discussed in [18]).In fact, the problem of scattering from bottom roughness is of great importance for shallowsea acoustics, mainly due to the fact that in real-world applications, the information on the bathymetry is often insufficiently accurate.To our knowledge, the first paper studying acoustic mode interaction in a waveguide with range-dependent bathymetry is [36].A kind of a wavefield propagator for a shallow sea was studied in [37], where a waveguide with randomly perturbed interfaces was considered, and the staircase approximation was used for the interface representation.In addition, we can mention the paper [38] devoted to the effect of scattering from the rough sea bottom on modal energies, and the paper [39] with the application of the cross-section method to a waveguide with random bathymetry.
In the present paper, we develop a theoretical approach that goes beyond the staircase approximation used in [37] and allows one to take into account the actual form of the sea bottom roughness.In addition, the particular case of the roughness described by the stochastic Ornstein-Uhlenbeck process is considered.
The paper is organized as follows.In the next section, we represent the general theory of the non-unitary matrix propagator for a shallow sea.Section 3 is devoted to the construction of the matrix propagator in the presence of bottom roughness.The particular case of the roughness described by the Ornstein-Uhlenbeck process is analyzed in Section 4. Section 5 presents the result of the numerical simulation of a model waveguide with a near-bottom sound channel.In the Discussion setion, we summarize the main results and outline the ways of further development of the approach.

General Theory
We start with the Helmholtz equation in cylindrical coordinates r, θ, z.Under the assumption that azimuthal coupling is weak (i.e., the derivative with respect to θ can be neglected), it can be written as where r is the range, z is the depth, P(r, z) is the acoustic pressure, ρ is the density, k 0 = 2π f /c 0 is a reference wavenumber, c 0 is a reference sound speed, and n(r, z) is the refractive index.Invoking the far-field approximation and neglecting back-scattering, we can transform Equation (1) to a one-way equation for the wavefield envelope function where the prefactor √ k 0 r takes into account cylindrical sound spreading, and the operator Q is given by the expression [40] The solution of Equation ( 2) can be formally written in terms of the propagator Ĝ defined as In order to facilitate the calculation of Ĝ, it is reasonable to use the basis of normal modes, i.e., solutions of the Sturm-Liouville problem [40] with properly imposed boundary and interface conditions [40].A wavefield can be represented as an expansion over normal modes: The amplitude of the m-th mode is determined as Hereafter, we consider the case of weak range dependence, assuming that the operator Q can be expressed as a sum of the range-independent part Q0 and a small perturbation V, i.e., ∂Ψ ∂r In the rest of this study we assume that modes in the Sturm-Liouville problem (5) are computed for Q = Q0 , i.e., for the reference range-independent waveguide.Substituting the ansatz (6) into Equation (8), and taking into account that Ψ m satisfies (5) for Q = Q0 , we obtain a system of coupled equations for the modal amplitudes a m (r): where are the matrix elements of the operator V. Equation ( 9) is similar to the perturbation theory for the Schrödinger equation in the way it is derived in [41], where the wave function of the perturbed (time-dependent) quantum system is represented in the form of expansion over the eigenstates of its counterpart with the time-independent Hamiltonian.Other approaches to the derivation of the mode coupling equations can be found in [11,40,[42][43][44].
Using the substitution we can transform Equation (9) to the following form: It is formally equivalent to a Volterra-type integral equation The respective expression for mode amplitudes a m reads as Basically, there are two kinds of modes in a shallow-water waveguide: waterborne modes with weak attenuation, and modes predominantly propagating inside the bottom sediment layer and experiencing strong losses.In practice, the latter ones are poorly known (due to the lack of data about precise bottom structure), and it is reasonable to exploit the fact that their contribution to the field is relatively small.This can be performed in the spirit of the Weisskopf-Wigner approach that is well known in the theory of open quantum systems [45].Without loss of generality, we can assume that waterborne modes correspond to m ≤ M, and modes with m > M have maxima in the bottom.
The expansion ( 6) can be rewritten as where M max is the total number of discrete-spectrum modes taken into account.Separating the contributions of water and sedimental modes in (12), we have Waterborne and sedimental modes overlap weakly; therefore, transitions between them can be considered a second-order effect.Except for very low frequencies, the sedimental modes have very strong attenuation; therefore, we can neglect the scattering of acoustic energy from sedimental to water modes.Since it is natural to assume that only the water modes are excited at r = r 0 , it means that solutions for the sedimental modes can be written as Substituting ( 17) into ( 16), we obtain equations for modes with m ≤ M If r is sufficiently large, only the terms with n = m efficiently contribute to this double sum, and we rewrite the latter equality as where The second term describes the mode attenuation caused by coupling to sedimental modes.Here, we do not imply any symmetries for the perturbation matrix, assuming that in general V mn = V nm .Below, we consider only the case of an operator V describing some stationary random process.It means that the corresponding matrix elements V mn (r) are stationary random functions as well, and Equation ( 19) is stochastic with non-Markovian kernels [15,46].The presence of V mn within the kernels anticipates fluctuations of attenuation.However, as long as perturbation is assumed to be weak, matrix elements V mn are small.Therefore, we can ignore fluctuations of attenuation and replace the product V mn V nm by its statistical average.Also, we can take into account that sedimental modes have much higher damping rates α m .This implies that the function rapidly decays with increasing |r − r |.This allows us to invoke the Markov approximation and ignore range non-locality.Thus, we obtain Here the damping coefficients Γ m are evaluated by means of the Wiener-Khintchin theorem, where Cmn (k) is Fourier-transformed C mn (r).Invoking the perturbation theory for Equation (21) and taking into account (11), we obtain the first-order approximation for the waterborne modes in the following form: where γ m = α m + Γ m .This equation can be rewritten in the matrix form as where I is the identity matrix, and Λ is a diagonal matrix corresponding to the unperturbed propagator and A is the perturbation matrix with elements Note that within this study bold letters denote matrices, e.g., A, while their respective elements are written as A mn .
The key idea of the approach based on random matrix theory is to represent matrix elements A mn as random Gaussian variables whose variance is determined by variance of the integral (26).This approach is well developed in the case of deep-water acoustic waveguides with volume inhomogeneities, when the propagator is nearly unitary [18,19,27,31].
The validity of the first-order perturbation theory demands the propagator step to be sufficiently small.In order to perform long-range propagation modeling, we can partition the waveguide into short segments of length ∆r, and represent the resulting propagator as a product Ĝ(r F , r 0 ) = Ĝ(r F , r F − ∆r) Ĝ(r F − ∆r, r F − 2∆r)... Ĝ(r 0 + ∆r, r 0 ), (27) or, in the matrix representation, Now the segment propagators can be obtained from Equation (24).

Random Sea Bottom Roughness
Let us now consider the case where the water depth randomly fluctuates with the range r, i.e., h(r where δh(r) is some random function to be specified later.To derive the expression for the perturbation operator V, it is convenient to introduce the rescaled depth variable where Similar transformations of variables are often used in quantum-mechanical problems with moving boundaries [47][48][49].We assume that variations of the water depth are weak, i.e., that η 1. Transformation of the depth variable (30) leads to the transformation of derivatives in (2) and (3): ∂ ∂r Substituting ( 32) and ( 33) into ( 2) and (3), we obtain As long as η is small compared to 1, we can separate the range-independent and range-dependent parts of n 2 (r, Z) and ρ(r, Z) using the Taylor expansion.Leaving only the zero-order and first-order terms, we have where For the sake of simplicity, we neglect the terms involving ∂ρ/∂Z.Also, we can take into account that 1 h Now Equation ( 34) can be rewritten as where differential operators Da and Db are given by expressions In order to find out the resulting perturbation operator V (see ( 8)), we have to extract the perturbation term containing η(r) from the square-root operator.This problem can be solved by means of the Taylor approximation (which is widely used in similar studies, see, for example, [50]): After such transformation, the perturbation operator reads as Matrix elements of the operator V are expressed as where and normal modes Ψ n (z) correspond to the unperturbed problem with V = 0. Elements of the perturbation matrix A are given by the formula A mn (∆r) = k 0 (υ a,mn D a,mn + υ b,mn D b,mn ) , where Here, υ a,mn and υ b,mn are stochastic integrals which depend on statistical properties of η(r).The presence of nonzero imaginary parts of wavenumbers k rn destroys symmetries of the matrix A. However, the attenuation of water modes is usually weak and obeys inequality k rn ∆r −1 .Therefore, we can neglect the impact of imaginary parts in Equations ( 47) and ( 48) and assume k rn k n .Then, Equations ( 47) and ( 48) read as

Ornstein-Uhlenbeck Sea Bottom Roughness
In this section, we consider the sea bottom roughness described by the formula where η(r) is the normalized stochastic Ornstein-Uhlenbeck process being the solution of the Langevin equation In this formula, r c is the roughness correlation length, and ξ is Gaussian white noise, δ(r) is the delta function, and angular brackets denote ensemble averaging.Roughness autocorrelation function is given by [51] Taking into account (52), we find where ξ a,mn is a Gaussian complex random variable with zero mean and unit variance.These Gaussian variables obey the relation In order to find the corresponding expression for υ b,mn , we can use spectral representation of the Ornstein-Uhlenbeck process: where φ(k) is a random phase function obeying and S(k) is the power spectrum of η(r) given by the formula Substituting ( 57) into (48), we obtain where ξ a,mn is another Gaussian random variable with zero mean and unit variance, also obeying and the corresponding variance is determined via the formula The resulting expression for perturbation matrix elements can be written as where In order to evaluate the coefficients of damping caused by coupling to sedimental modes, one has to compute the correlator C mn (r , r ) (20).Let us denote η = dη/dr and recall that ξ(r)η(r − s) = 0 .
Then, according to (52), we have where δ(s) is the delta function.The correlator can therefore be written as Applying the Wiener-Khintchin theorem, we find that where

Numerical Simulation
For the numerical simulation, we use a model of a shallow-water waveguide with the sound-speed profile described by the formula [52][53][54]  The magnitude of the bottom roughness is σ η = 0.01; this corresponds to the r.m.s.displacement of 1 m.The correlation length of roughness r c is varied from 100 to 1000 m.The sediment attenuation coefficient is taken to be 0.5 dB per wavelength.
Calculation of the propagator described above is based on the mode perturbation theory that provides accurate prediction only for relatively low signal frequencies.We carried out the statistical modeling of propagator matrices for two values of frequency, namely 100 and 200 Hz.It was found that the number of waterborne modes in the considered model of a waveguide is related to the acoustic frequency via the formula f 20 M. That is, we have 5 waterborne modes for 100 Hz, and 10 waterborne modes for 200 Hz.The propagator step ∆r is taken to be equal to the correlation length of the roughness ∆r = r c .This ensures the applicability of the perturbation theory to the evaluation of the matrix elements of the propagator, and, on the other hand, meets the condition of statistical independence of the propagators for neighboring waveguide segments.This condition allows us to neglect the correlations between the segment propagators.For each set of the parameters, a statistical analysis is carried out using 1000 realizations of propagator matrices.
Firsly, let us examine range variations of modal energies.In the case of a single-mode source, normalized mean energies are determined by the formula The angular brackets in this formula denote averaging over the ensemble of roughness realizations.Decreasing the roughness correlation length enhances the impact of scattering and transitions of mode energy to higher-order modes, exhibiting faster attenuation by the bottom.Therefore, in this case, the scattering increases energy losses.It is confirmed by the data presented in Figures 2 and 3. Notably, in the case of the acoustic frequency of 100 Hz, the decay is exponential (this is revealed by the linear decreasing of the energy logarithm).
In the case of f = 200 Hz, higher-order modes decay slower than the exponential rate due to the partial re-pumping from much more energetic lower-order modes.The presence of energy flux from lower to higher modes due to mode coupling is illustrated in Figure 3a,b, where the range dependence of the first mode energy is depicted.It is shown that decreasing r c , anticipating stronger coupling, leads to faster attenuation.
- Decay of the first mode energy consists of two components: energy transfer to the sedimental modes with high loss, and energy transfer to the waterborne modes with relatively low loss.To eliminate the contribution of the latter component, we can estimate the roughness-induced enhancement of losses.In particular, we can consider the depthintegrated transmission loss for a wavefield created by a single-mode source.Using the propagator, it can be calculated using the formula where r ref = 1 m.The corresponding data for n = 1 are presented in Figure 3c,d.It turns out that the coupling of the first mode to lossy sedimental modes is fairly weak; therefore, the bottom roughness weakly influences transmission loss.Another characteristics that can be used for the description of mode coupling is modal cross-coherence defined as where range dependencies of modal amplitudes are calculated for the initial conditions of the following form: We focus our attention on the cross-coherence of the first mode with other waterborne modes.Data for frequency of 100 Hz are presented in Figure 4.The data suggest that the first three modes maintain a high degree of coherence throughout the entire distance under consideration.On the other hand, the fifth mode experiencing the strongest scattering among the waterborne modes undergoes the fastest decoherence.
In the case of f = 200 Hz (see Figure 5), we also consider the first five modes.However, as the frequency is increased, these modes correspond to flatter rays [55] that are less in contact with the sea bottom roughness.Therefore, they expose much weaker scattering: considerable decoherence is observed only in the case of short-range roughness.In the cases of r c = 500 and 1000 m, the cross-coherence decay is very weak.Fluctuations of modal amplitudes can be estimated by means of the corresponding scintillation index defined as The data presented in Figure 6 demonstrates the growth of fluctuations as the unperturbed (i.e., calculated in the unperturbed waveguide without roughness) modal attenuation rate increases.However, with a distance of 20 km, only the highest waterborne modes manage to achieve statistical saturation with SI 1 [56].Moreover, in the case of f = 100 Hz (see Figure 4a,b), even the highest five modes do not achieve plateau, despite the high values of the scintillation index.It means that saturation for low-frequency propagation should occur at longer distances.According to Figure 7, the first mode exposes very weak fluctuations for all values of the roughness correlation length, anticipating nearly adiabatic propagation.
Comparing Figure 7a,b, a twofold increase in frequency leads to a twofold increase in the scintillation index, implying that the fluctuation strength is proportional to the frequency.

Discussion
In the present study, we outline a novel approach to the modeling of sound scattering by random sea bottom roughness.This approach is based on the construction of random propagator matrices.It is a very computationally efficient tool for wavefield modeling that has many other attractive features as well.In particular, the propagator takes into account almost all information about the acoustic properties of the waveguide segment under consideration.Therefore, it can provide comprehensive information about the scattering physics, including statistics of losses and fluctuations, and estimates of crossmodal coherence, as as the behavior of many other wavefield parameters.
We derive explicit analytical expressions for the propagator matrix elements for roughness described by the stochastic Ornstein-Uhlenbeck process.It is very important from the viewpoint of practical implementation, as the Ornstein-Uhlenbeck process represents a typical example of roughness with exponentially decaying autocorrelation function that can be a good model of realistic water depth variations in a shallow sea.Statistical analysis of sound propagation in a shallow-sea with the Ornstein-Uhlenbeck bottom roughness shows that the low-order modes are weakly affected by the scattering and exhibit highly unsaturated statistics of fluctuations.Nevertheless, mode coupling due to scattering results in significant increasing of attenuation rates of the lowest-order modes.Decreasing the roughness correlation length leads to drastic enhancement of the scattering, especially for the highest waterborne modes.It is reflected as the fast decay of cross-mode coherence functions, as well as the sharp increase in the modal scintillation indices.Such enhancement of scattering was reported in [39].
Further development of the propagator approach based on the random matrix theory, in our opinion, should involve generalization to 3D propagation scenarios [57][58][59].This can be accomplished, for example, in the framework of pseudodifferential mode parabolic equations [60,61].Indeed, the effect of the horizontal refraction on the field at a reception point can be very significant in a shallow sea.Some phenomena similar to the branching of light beams in optics can be expected in this case [62].Another promising direction of generalization of the results presented here is the development of the random matrix theory for broadband signals that would accurately handle inter-frequency correlations.Finally, the results of the present work can be extended to the case of an elastic bottom by using the approach proposed in [63].Arguably, bottom elasticity would not significantly change the results of the numerical experiments reported here, as one can expect its contribution to manifest primarily in somewhat greater attenuation coefficients of the bottom modes.

where c 0 =Figure 1 .
Figure 1.Sound-speed profile for the waveguide model used in numerical simulation .

Figure 2 .
Figure 2. Dependence of mode energies on range.Values of acoustic frequency: 100 Hz for panels (a,b), and 200 Hz for panels (c,d).Values of roughness correlation length r c : 100 m for panels (a,c), and 500 m for panels (b,d).

Figure 3 .
Figure 3. Range dependence of the first mode energy (panels (a,b)) and depth-integrated transmission loss of a wavefield with the initial condition in the form of the first mode.Values of acoustic frequency: 100 Hz (a,c), and 200 Hz (b,d).Values of roughness correlation length are indicated in the figure.

Figure 6 .Figure 7 .
Figure 6.Dependence of mode scintillation indices on range.Values of acoustic frequency: 100 Hz for panels (a,b), and 200 Hz for panels (c,d).Values of roughness correlation length r c : 100 m for panels (a,c), and 500 m for panels (b,d).