Ray-based description of mode coupling by sound speed fluctuations in the ocean
aa r X i v : . [ phy s i c s . a o - ph ] J u l Ray-based description of mode coupling by sound speedfluctuations in the ocean
A.L. VirovlyanskyInstitute of Applied Physics, Russian Academy of Science,46, Ulyanova St., Nizhny Novgorod, Russia, [email protected] 5, 2018
Abstract
A traditional approach to the analysis of mode coupling in a fluctuating underwaterwaveguide is based on solving the system of coupled equations for the second statisti-cal moments of mode amplitudes derived in the Markov approximation [D.B. Creamer,J.Acoust. Soc. Am. 99, 2825–2838 (1996)]. In the present work an alternative approachis considered. It is based on an analytic solution of the mode coupling equation de-rived in the high frequency approximation [A.L. Virovlyanskii and A.G. Kosterin, Sov.Phys. Acoust. 35, 138-142 (1987)]. This solution representing the mode amplitude asa sum of contributions from two geometrical rays is convenient for statistical averag-ing. It allows one to easily derive analytical expressions for any statistical momentsof mode amplitudes. The applicability of this approach is demonstrated by comparingits predictions for a deep water acoustic waveguide with results of full wave numericalsimulation carried out using the method of wide angle parabolic equation.
PACS numbers: Introduction
An environmental model representing an unperturbed range-independent waveguide withweak sound speed fluctuations is widely used in solving different problems of ocean acoustics[1, 2]. The sound field in this model can be decomposed into a sum of normal modes ofthe unperturbed waveguide with complex amplitudes a m being random functions of range r . Mode coupling by the sound speed fluctuations is quantitatively described by statisticalmoments of mode amplitudes h a m a ∗ n i where the asterisk denotes complex conjugation andthe angular brackets denote statistical averaging. A traditional approach to the analysis ofthe modal structure of the sound field in a fluctuating waveguide is based on solving thetransport equations for moments h a m a ∗ n i [3–8].Complete system of transport equations derived in the Markov approximation includes M linear differential equations, where M is the number of propagating modes [6, 8]. Atruncated version of this system including only M equations for the mean mode intensities (cid:10) | a m | (cid:11) – they are called the master equations – was derived in Ref. [4] under assumptionthat the cross-mode coherences, that is, the moments h a m a ∗ n i with m = n , are negligible.In Ref. [8] all the statistical moments h a m a ∗ n i were calculated numerically for a deepwater acoustic waveguide with sound speed fluctuations induced by random internal waves.It turned out that the values of |h a m a ∗ n i| for close, but not equal, m and n are of thesame order of magnitude as (cid:10) | a m | (cid:11) / (cid:10) | a n | (cid:11) / and the mean intensities and cross-modecoherences ’evolve over similar range scales’. This contradicts the assumption about thesmallness of the cross-mode coherences made when deriving the master equations. In spiteof this fact, the numerical simulation demonstrated that the solutions of master equationsprovide good estimates of mean mode intensities (cid:10) | a m ( r ) | (cid:11) . However, it should be noticedthat these solutions are smooth functions of range and they do not predict small oscillationsof (cid:10) | a m ( r ) | (cid:11) with range found by solving the complete system of all the M equations.Similar results were obtained in Refs. [9, 10] for a shallow water waveguide.In the present paper we consider an alternative approach to examining the mode coupling.It was derived in Ref. [11] (see also Refs. [12–15]) in the study of ray-mode relations in awaveguide with weak sound-speed fluctuations. This approach is based on a surprisinglysimple analytical solution of the mode coupling equation obtained in the high frequencyapproximation. It expresses the mode amplitude through parameters of two geometrical rayswhich we call the mode rays. This estimate of a m is an analog of the well-known formula2f geometrical optics describing the variation of complex ray amplitude in the presence ofweak inhomogeneities of refractive index. The ray-based estimate of mode amplitude isconvenient for statistical averaging. Analytical expressions for any statistical moments ofmode amplitudes, including joint moments of amplitudes at different frequencies, are readilyfollow from this formula.We also consider a simplified expression for the mean mode intensity obtained by averag-ing the parameters of the ray-based estimate of (cid:10) | a m | (cid:11) over the ray cycle length. It is shownthat in the high frequency limit this expression, describing a smoothed range-dependenceof the mean mode intensity, satisfies the master equations. This fact agrees with numericalresult of Refs. [8–10].A weak point of our analytical approach is insufficient knowledge about the limits ofits applicability. So far, its predictions have never been compared to results of full wavesimulation. In the present paper such a comparison is made for an environmental modelsimilar to that used in Ref. [8]. Sound field excited by a points source at a frequency of100 Hz was calculated numerically using the method of wide angle parabolic equation in360 realizations of the fluctuating waveguide. Mode amplitudes were found by projectingcomputed wave fields onto eigenfuntions of the unperturbed waveguide. Parabolic equationbased (pe-based) estimates obtained in this way are compared with the ray-based estimatesfound by evaluating mode amplitudes in the same realizations of the waveguide using ourapproximate analytical solution of the mode coupling equation. Estimates of statisticalmoments were computed by the Monte Carlo method, that is, by averaging the productsof mode amplitudes over all the realizations of random waveguide. The comparison hasdemonstrated a good agreement between the pe-based and ray-based results.Numerical simulation have shown that the ray-based approach properly describes not onlythe smoothed range-dependencies of statistical moments of mode amplitudes but the smalloscillations of these moments (”missed” by the master equations), as well. These oscillationswere predicted in Ref. [12] where it was shown that the jump-like variations of the modeamplitude and its statistical moments occur in the neighborhood of the upper turning pointsof the mode rays.As in Refs. [4, 6, 8–10], we neglect the horizontal refraction of sound waves and consider atwo-dimensional environmental model. Out-of-plane wave scattering was taken into accountin Refs. [16, 17].The organization of this paper is as follows.3nalytical relation expressing the mode amplitude through parameters of two ray pathsis presented in Sec. II. Section III gives analytical expressions for a few statistical momentsof mode amplitudes derived using this relation. In Sec. IV, it is shown that in the limit ofhigh frequency the ray-based estimate of the mean mode intensity (cid:10) | a m | (cid:11) whose parametersare smoothed over the ray cycle length satisfies the master equations. Section V presentsresults of numerical simulation in a deep water waveguide with sound speed fluctuationsinduced by random internal waves. It is demonstrated that the predictions of our ray-basedapproach agree with the results of simulation carried out using the method of wide angleparabolic equation. In Sec. VI, the results of this work are summarized. II Analytical description of mode amplitudes in the pres-ence of weak sound speed fluctuations
In this section we present a simple analytical approach derived in Refs. [11, 14, 15] for aray-based description of mode amplitudes in a waveguide with weak large scale sound speedfluctuations. It is assumed that the wave field is excited by a point source.
A Mode representation of the wave field
Consider a two dimensional model of underwater sound channel with the sound speed field c ( r, z ) = ¯ c ( z ) + δc ( r, z ) , where r is the distance, z is the depth, ¯ c ( z ) is the unperturbed soundspeed profile, and δc ( r, z ) is the weak range-dependent perturbation. The refractive index is ν ( r, z ) = c /c ( r, z ) , where c is the reference sound speed. We assume that | c ( r, z ) − c | ≪ c . Due to the weakness of perturbation ν ( r, z ) = ¯ ν ( z ) + δν ( r, z ) , where ¯ ν ( z ) = c ¯ c ( z ) , δν ( r, z ) = − c ¯ c ( z ) δc ( r, z ) . We assume that the perturbation δc ( r, z ) is a zero mean Gaussian random field with thecorrelation function h δc ( r , z ) δc ( r , z ) i = K ( ξ, ζ , Z ) , (1)4here ξ = r − r , ζ = z − z , Z = ( z + z ) / . Note that K ( ξ, ζ , Z ) = K ( − ξ, ζ , Z ) . (2)Characteristic scales of correlation function K along the coordinates ξ and ζ denote ∆ ξ and ∆ ζ , respectively.The acoustic pressure field u ( r, z ) at a carrier frequency f can be expressed as u ( r, z ) = M X m =1 r πik m r a m ( r ) ϕ m ( z ) e ik m r , (3)where k m and ϕ m ( z ) are eigenvalues and eigenfunctions of the Sturm-Liouville problemin the unperturbed (range-independent) waveguide, respectively [2, 18]. For simplicity, itis assumed that the sum (3) includes only those modes whose turning points are locatedwithin the water bulk. This assumption will simplify the use of the WKB approximation fordescription of k m and ϕ m ( z ) .In what follows we will consider the wave field excited by a point source set at r = 0 and z = z . In this case, a m (0) = ϕ m ( z ) .In the WKB approximation the eigenvalues can be presented as k m = kh m , where k =2 πf /c is a reference wavenumber and h m is determined by the quantization rule [2] k Z z max z min dz p ¯ ν ( z ) − h m = π ( m − / (4)with z min and z max being the mode turning depths. In this approximation the m -th modeis associated with a ray path whose grazing angle at depth z , θ m ( z ) , is determined by therelation h m = ¯ ν ( z ) cos θ m ( z ) .The cycle length (period) of this ray path is given by D m = 2 h m Z z max z min dz p ¯ ν ( z ) − h m = 2 Z z max z min dz tan ( θ m ( z )) . (5)A ’differential’ form of the quantization rule follows from Eqs. (4) and (5) as dh m dm = − πkD m . (6)5he eigenfunction ϕ m ( z ) in the WKB approximation [2] can be presented in the form ϕ m ( z ) = ϕ + m ( z ) + ϕ − m ( z ) , (7)where ϕ ± m ( z ) = q m ( z ) e ± i [ kg m ( z ) − π/ , (8) q m ( z ) = h / m [¯ ν ( z ) − h m ] / D / m = 1[ D m tan θ m ( z )] / , (9) g m ( z ) = Z zz min dz p ¯ ν ( z ) − h m . (10)Functions ϕ ± m ( z ) represent two quasi-plane waves called the Brillouin waves.Note a useful relation following from Eqs. (4) and (10) k ∂g m ( z ) ∂m = 2 πD m Z zz min dz tan ( θ m ( z )) . (11) B Geometrical optics for modes
Within the framework of standard geometrical optics, the influence of a weak sound speedperturbation δc with spatial scales significantly exceeding the wavelength can be accountedfor using a well-known approximate formula. If in the unperturbed medium the contributionof a sound ray to the total field u is A exp ( ikS ) , where A and S are the ray amplitude andeikonal, respectively, then in the presence of perturbation its contribution becomes [1, 2] u = Ae ik ( S + X ) , (12)where X = Z Γ δν ds, (13) ds is the arc length and the integration goes over the unperturbed ray path Γ . Although thisformula is valid only at relatively short ranges it is widely used in the ocean acoustics [1]. Inparticular, it is used in solving the inverse problem in the classical scheme of ocean acoustictomography [19].In Ref. [11] (see also Refs. [12–15]) it is shown that there exists a close analog of Eq. (12)for normal modes. The point is that the m -th mode constructively interferes (adds in phase)with neighboring modes along the trajectories of two unperturbed rays leaving the source6t launch angles ± θ m ( z ) which are equal to grazing angles of the Brillouin waves ϕ ± m ( z ) at the source depth. As in Refs. [11–15], we shall call these rays the mode rays and denotetheir trajectories z ± m ( r ) . Note that for a given m the angle θ m ( z ) is a function of the carrierfrequency f . This makes the trajectories z ± m ( r ) frequency dependent. Denote the grazingangles of the ray paths z ± m ( r ) at range r by χ ± m ( r ) , so that dz ± m ( r ) /dr = tan χ ± m ( r ) . Bothmode rays have the same cycle length given by Eq. (5). Examples of mode rays are shown inFig. 1. It graphs trajectories z ± m ( r ) for the 36-th mode in the canonical sound speed profileat a carrier frequency of 100 Hz.In the presence of perturbation δc ( r, z ) the mode amplitude is expressed by the approx-imate formula a m ( r ) = ϕ + m ( z ) e ikX + m ( r ) + ϕ − m ( z ) e ikX − m ( r ) , (14)with X ± m = Z Γ ± m ds δν, (15)where the integration goes along the trajectories of mode rays Γ ± m [11, 15]. Equation (15)can be written in the form X ± m ( r ) = Z r dr ′ cos χ ± m δν (cid:0) r, z ± m ( r ) (cid:1) . (16)Formula (14) is derived under the same assumptions as its prototype for the ray amplitude(12). The simplest derivation of Eq. (14) consists in projecting the ray representation of thewave field onto eigenfunctions of the unperturbed waveguide with the evaluation of arisingintegrals using the stationary phase technique [14, 15]. Therefore Eq. (14) should haveapproximately the same range of applicability as Eq. (12).A more accurate and general expression for the mode amplitude can be derived proceed-ing from the ray representation of the wave field in a range-dependent waveguide [20–22].Besides, Eq. (14) can be generalized in a different direction. In Ref. [11, 14, 15], a moregeneral version of this formula was derived which can be used for description of wave diffrac-tion by sound speed fluctuations. In these works, the notion of Fresnel zones for modes isintroduced which is analogous to the usual Fresnel zones introduced for rays.In the present paper the indicated generalizations are not considered. All our subsequentanalysis is based on formula (14). 7 .5 1.52 1.540123 z ( k m ) c (km/s) r (km) D m z (r) z (r) χ Figure 1: Left panel: canonical sound speed profile. Right panel: trajectories of the moderays of the 36-th mode at a carrier frequency of 100 Hz.
III Statistical moments of mode amplitudes
Since X ± m ( r ) are zero mean Gaussian random functions, an analytical expression for anystatistical moment of mode amplitudes is readily derived from Eq. (14) using the well-known formula h e iα i = − e − h α i / for a zero mean Gaussian random variable α . The meanvalue (coherent component) of the mode amplitude a m is h a m i = Φ + m e − k D ( X + m ) E + Φ − m e − k D ( X − m ) E , (17)where Φ ± m = ϕ ± m ( z ) .The cross-mode coherence is given by h a m a ∗ n i = Φ + m Φ − n e − k D ( X + m − X + n ) E + Φ + m Φ + n e − k D ( X + m − X − n ) E + Φ − m Φ − n e − k D ( X − m − X + n ) E + Φ − m Φ + n e − k D ( X − m − X − n ) E . (18)In the particular case m = n , Eq. (18) gives an expression for the mean mode intensity (cid:10) | a m | (cid:11) = 2 Q m (cid:20) e − k D ( X + m − X − m ) E sin(2 kG m ) (cid:21) , (19)where Q m = q m ( z ) = 1[ D m tan θ m ( z )] / , m = g m ( z ) = 2 πD m Z z z min dz tan ( θ m ( z )) . An expression for the mean squared intensity (the fourth moment of mode amplitude) is (cid:10) | a m | (cid:11) = Q m (cid:20) kG m ) e − k D ( X + m − X − m ) E − kG m ) e − k D ( X + m − X − m ) E (cid:21) . (20)It is clear that similar formulas are readily derived for the joint statistical moments ofmode amplitudes at different frequencies.In the scope of our ray-based approach, all the moments of mode amplitudes are ex-pressed through h X ± m X ∓ n i and h X ± m X ± n i . Evaluation of these quantities is simplified underthe assumption that the horizontal correlation scale of sound speed fluctuations, ∆ r , is sub-stantially less than the cycle length D m . Then, at ranges r ≫ ∆ r the mode rays crossuncorrelated inhomogeneities, the quantities X + m and X − m become statistically independent,and h X ± m X ∓ n i = 0 .Explicit expression for the dispersions of X ± m is given by D(cid:0) X ± m (cid:1) E = Z r Z r dr ′ dr ′′ g m ( r ′ ) g m ( r ′′ ) × K (cid:18) r ′ − r ′′ , z ± m ( r ′ ) − z ± m ( r ′′ ) , (cid:0) z ± m ( r ′ ) + z ± m ( r ′′ ) (cid:1)(cid:19) , (21)where g m ( r ) = ¯ c ( z ± m ( r )) cos ( χ ± m ( r )) . Let us change the variables of integration from ( r ′ , r ′′ ) to ( r ′ , ξ ) , where ξ = r ′ − r ′′ . The main contribution to the integral comes from the interval | ξ | < ∆ ζ . Since ∆ ζ is small compared to D m , we can use the approximations z ± m ( r ′ ) − z ± m ( r ′′ ) ≃ ξ tan χ ± m ( r ′ ) and ( z ± m ( r ′ ) + z ± m ( r ′′ )) / ≃ z ± m ( r ′ ) . At ranges r ≫ ∆ ζ we can formally extendthe limits of integration over ξ to infinity. Then Eq. (21) translates to D(cid:0) X ± m (cid:1) E = Z r dr ′ ¯ c ( z ± m ( r ′ )) cos χ ± m ( r ′ ) Z ∞−∞ dξ K (cid:0) ξ, ξ tan χ ± m ( r ′ ) , z ± m ( r ′ ) (cid:1) . (22)Since the trajectories z + m ( r ) and z − m ( r ) differ only by a shift along the r -axis, the contributionto the integral over r ′ from any interval of length D m is the same for both mode rays. It iseasy to see that for an arbitrary function F ( z ± m ( r ′ ) , | χ ± m ( r ′ ) | ) we have the relation D m Z r + D m r dr ′ F ( z ± m ( r ′ ) , (cid:12)(cid:12) χ ± m ( r ′ ) (cid:12)(cid:12) ) = 2 D m Z z max z min dz tan θ m ( z ) F ( z, θ m ( z )) .
9t ranges r = N D m with N being an integer D(cid:0) X + m (cid:1) E = D(cid:0) X − m (cid:1) E = (cid:10) X m (cid:11) , (23)where (cid:10) X m (cid:11) = r D m c k k m Z z max z min dz ¯ c ( z ) tan θ m ( z ) Z ∞−∞ dξ K ( ξ, ξ tan θ m ( z ) , z ) . (24)In deriving this formula we have taken into account that due to Snell’s law cos χ ± m =( k m /k )(¯ c ( z ± m ) /c ) . At an arbitrary range r , not necessary multiple of D m , Eq. (24) gives asmoothed estimate of D ( X + m ) E and D ( X − m ) E .Let us slightly simplify formulas (17) – (19) neglecting the differences between D ( X + m ) E and D ( X − m ) E and between D ( X + m − X + n ) E and D ( X − m − X − n ) E . Assuming that h X ± m X ∓ m i = 0 and replacing D ( X ± m ) E by (cid:10) ( X m ) (cid:11) and D ( X ± m − X ± n ) E by (cid:10) ( X m − X n ) (cid:11) , we find h a m i = ϕ m ( z ) e − k h X m i , (25) h a m a ∗ n i = 2 Q m Q n h e − k h ( X m − X n ) i cos ( k ( G m − G n ))+ e − k ( h X m i + h X n i ) sin( k ( G m + G n )) i , (26) (cid:10) | a m | (cid:11) = 2 Q m h e − k h X m i sin(2 kG m ) i , (27)and (cid:10) | a m | (cid:11) = Q m h kG m ) e − k h X m i− kG m ) e − k h X m i i . (28)The value of h X m i is given by Eq. (24). We do not present an explicit expression for (cid:10) ( X m − X n ) (cid:11) . When using Eq. (26), exp h − k (cid:10) ( X m − X n ) (cid:11)i should be replaced by anyof two close functions exp h − k D ( X ± m − X ± n ) Ei . IV Ray-based approach and master equations
The equations for statistical moments h a m a ∗ n i are derived in the Markov approximationproceeding from the mode coupling equation [4, 6, 8]10 a m dr = i M X n =1 ρ mn a n e ik nm r , (29)where k nm = k n − k m , ρ mn ( r ) = k √ k m k n Z dz ϕ m ( z ) µ ( r, z ) ϕ n ( z ) , (30) µ ( r, z ) = − c ¯ c ( z ) δc ( r, z ) . (31)Note that our main formula (14) is an approximate solution of Eq. (29) [11].In Refs. [8–10], it was shown that the numerical solution of the complete system of equa-tions for all the moments h a m a ∗ n i give practically the same result as the evaluation of thesemoments in the Monte Carlo simulation based on numerical solving the mode coupling equa-tion (29) for different realizations of random perturbation δc ( r, z ) . This was the expectedresult. An unexpected result was that even though the cross-mode coherences were notsmall, the master equations properly predicted the smoothed mode intensities. In this sec-tion we will show that this result follows from our ray-based estimates of statistical moments.Namely, it will be shown that Eq. (27) obtained by smoothing the range-dependent param-eters of Eq. (19), in the limit of high frequency gives a solution to the master equations.In the notation of Refs. [8–10] the master equations have the form d (cid:10) | a m | (cid:11) dr = 2 M X n =0 Re ( I mn,nm ) (cid:0)(cid:10) | a n | (cid:11) − (cid:10) | a m | (cid:11)(cid:1) , (32)where I mn,nm are the elements of the scattering matrix defined by the relations I mn,qp = Z ∞ dξ ∆ mn,qp ( ξ ) e ik pq ξ (33)and ∆ mn,pq ( r − r ′ ) = h ρ mn ( r ) ρ pq ( r ′ ) i . (34)According to Eqs. (30), (31), and (34), ∆ mn,nm ( ξ ) = k c k n k m Z dzdz ′ ϕ n ( z ) ϕ n ( z ′ ) K ( ξ, z − z ′ , ( z + z ′ ) / c ( z ) ¯ c ( z ′ ) ϕ m ( z ) ϕ m ( z ′ ) . (35)At high frequencies, where the wavelength is small compared to the spatial scales ofperturbation δc , we deal with the small-angle forward scattering of sound waves, and eachmode couples mainly into modes with close numbers. This means that the main contributionto the sum in the right hand side of Eq. (29) comes from terms with n close to m . UsingEq. (7), we present the product of two eigenfunctions with close m and n in the form ϕ n ( z ) ϕ m ( z ) ≃ q m ( z ) (cid:8) e ik [ g n ( z ) − g m ( z )] + c.c. (cid:9) , (36)where c.c. denotes the complex conjugate of the preceding term. In the right hand side ofEq. (36) we have omitted rapidly oscillating terms whose contributions to the integral inEq. (35) are negligible. The integrand on the right of Eq. (35) is non-negligible only for | z − z ′ | = O (∆ ζ ) . We assume that the vertical scale of perturbation ∆ ζ is small comparedto the depth interval between turning points of the m -th mode. Then the Brillouin waveswithin the depth interval of width ∆ ζ can be approximated by plane waves. Substitute Eq.(36) in Eq. (35) and drop the rapidly oscillating terms. Using Eqs. (6) and (11), the phasesof the remaining terms can be represented as k [ g n ( z ) − g m ( z ) − g n ( z ′ ) + g m ( z ′ )] ≃ k ∂ g m ( z ) ∂z∂m ( n − m ) ( z ′ − z )= 2 π ( n − m ) ( z − z ′ ) D m cot θ m ( z ) . In the resulting expression, we approximately replace n and z ′ in the pre-exponential factorsby m and z , respectively. This yields ∆ nm,mn ( ξ ) = k c k m Z z max z min dz Z ∞−∞ dζ × K ( ξ, ζ , z )¯ c ( z ) D m tan θ m ( z ) h e πi ( n − m ) Dm ζ cot θ m ( z ) + c.c. i . (37)Let us plug Eq. (37) into Eq. (33) and use the relation k nm = 2 π ( m − n ) /D m which followsfrom Eq. (6). At high frequencies, the number of propagating mode becomes very large andthe sum P n . . . in Eq. (32) can be approximately replaced by the integral R dn . . . . For m satisfying the condition ≪ m ≪ M , we will use the approximate relation Z M dn exp (cid:20) πi ( m − n ) D m ( ξ ± ζ cot θ m ) (cid:21) = D m δ ( ξ ± ζ cot θ m ) . M X n =0 Re I mn,nm = Z ∞ dξ Z M dn ∆ mn,nm ( ξ ) cos (cid:20) πi ( n − m ) D m ξ (cid:21) = k c k m D m Z dz ¯ c ( z ) tan θ m ( z ) Z ∞ dξ [ K ( ξ, ξ tan θ m , z ) + K ( ξ, − ξ tan θ m , z )] . (38)From the comparison of this expression with Eq. (24) (taking into account Eq. (2)) we find M X n =0 Re ( I mn,nm ) = k d h X m i dr . (39)Plugging Eq. (27) into the left-hand side (l.h.s.) and right-hand side (r.h.s.) of Eq. (32)yields: l.h.s. = 2 k Q m sin (2 kG m ) e − k h X m i d h X m i dr , (40)and r.h.s. = A + B + C, (41)where A = 4 Q m e − k h X m i sin(2 kG m ) M X n =0 Re ( I mn,nm ) , (42) B = 4 M X n =0 Re ( I mn,nm ) (cid:0) Q n − Q m (cid:1) , (43) C = − M X n =0 Re ( I mn,nm ) Q n sin(2 kG n ) e − k h X n i . (44)According to Eq. (39), l.h.s. = A . This means that l.h.s. = r.h.s. if term A dominates insum (41).According to Eq. (27), the mean intensity (cid:10) | a m | (cid:11) varies at ranges where k h X m i = O (1) .It can be shown that at these ranges and at sufficiently large k the term A dominates inthe sum (41). The smallness of term B is caused by the fact that Q m is a smooth functionof the mode number m . Analysis of Eqs. (33) and (37) shows that the values of ∆ mn,nm and I mn,nm descrease with increasing | n − m | and the main contributions to sums (42) and(43) come from terms with n belonging to some interval | n − m | < ∆ m . Consider Brillouinwaves with grazing angles close to some fixed value. It is easy to show that the numbers m
13f corresponding modes grow with frequency but the values of Q m and ∆ m for these modeswill be approximately constant. It means that (cid:12)(cid:12) Q m +∆ m − Q m (cid:12)(cid:12) /Q m = O (1 /k ) and thereforein the high frequency limit the ratio B/A tends to zero.The smallness of C compared to A is caused by the presence in Eq. (44) a rapidlyoscillating factor sin(2 kG n ) . This can be shown by transforming sum (44) in the samemanner as it has been done for sum P n Re I mn,nm . V Numerical example
In this Section we present results of numerical simulation demonstrating the applicability ofEq. (14) and estimates of statistical moments obtained using this formula. As in Ref. [8] weconsider a deep-water waveguide with the canonical sound speed profile and perturbation δc ( r, z ) induced by random internal waves. A Environmental model and numerical simulation
In numerical simulations presented below we use an environmental model with an unper-turbed sound speed profile representing the canonical (or Munk) profile [1, 2] ¯ c ( z ) = c r [1 + ε ( e η − η − , η = 2( z a − z ) /B (45)with parameters c r = 1 . km/s, ε = 0 . , B = 1 km, and z a = 1 km. This ¯ c ( z ) is shownin the left panel of Fig. 1. The bottom was set at a depth of 5 km.It is assumed that the weak perturbation δc ( r, z ) is caused by random internal waves withstatistics determined by the empirical Garrett-Munk spectrum [1]. To generate realizationsof a random field δc ( r, z ) we apply a numerical technique developed by J. Colosi and M.Brown [23]. In their model the perturbation has the form δc ( r, z ) = c r µg N ζ ( r, z ) , (46)where g = 9 . m/s is the gravitational acceleration, µ = 24 . is a dimensionless constant, N ( z ) = N exp( z/L ) is a buoyancy frequency profile, N = 2 π/ (12 min ) = 0 . L = 1 km. The random function ζ ( r, z ) presentsinternal-wave-induced vertical displacements of a fluid parcel. Its realizations have beencomputed using Eq. (19) from Ref. [23]. We consider an internal wave field formed by14 normal modes and assume its horizontal isotropy. Components of wave number vectorsin the horizontal plane belong to the interval from π/ km − to π/ km − . An rmsamplitude of the perturbation scales in depth like exp(3 z/ L ) and its surface-extrapolatedvalue in our model is about . m/s.All the calculations were carried out at a carrier frequency of 100 Hz. The point sourceexciting the wave field was set at the sound channel axis z = z a . The complex amplitudeof the wave field was computed using the method of wide angle parabolic equation for 360realizations of random perturbation δc ( r, z ) . Parabolic equation was solved by applying theCrank–Nicolson finite-difference scheme [18]. For each realization of δc ( r, z ) , the complexamplitude of the sound field was computed up to 500 km. Starting field at r = 0 wasgenerated using the modal starter [18]. Values of a m at 501 range points uniformly samplingthe interval from 0 to 500 km were found by projecting the computed sound field ontoeigenfunctions ϕ m ( z ) . Functions a m ( r ) obtained this way we call the pe-based estimates ofmode amplitudes.Our attention was restricted to amplitudes of the first 66 modes which describe soundwaves propagating at grazing angle | χ | < . ◦ . Turning points of these modes are locatedwithin the water bulk and far enough from the boundaries for the applicability of quantizationrule (4). Starting intensities | a m (0) | of some of these modes are shown by circles in Fig. 2.Numbers of modes whose statistical moments will be shown on the plots presented below,are indicated next to the corresponding circles.Functions X + m ( r ) and X − m ( r ) determined by Eq. (16) were computed for the same 360realizations of δc ( r, z ) at the same 501 range points. Then, substituting these functions inEq. (14) we obtained the ray-based estimates of a m ( r ) .Thus, in each of 360 realizations of δc ( r, z ) for each of the first 66 modes we computeda pe-based and ray-based estimate of mode amplitude a m ( r ) . For most modes these twoestimates of a m ( r ) are in reasonable agreement. Figure 3 present typical examples of | a m ( r ) | computed these two ways for the same realization of perturbation. B Statistical moments of a m In the remaing part of this paper we will compare the pe-based and ray-based estimates ofstatistical moments calculated using the Monte Carlo method. In what follows, the angularbrackets h . . . i denote the averaging over the 360 realizations of perturbation δc .15 m | a m ( ) | ( k m − )
28 38 56 5758 16 4227 262521 24 52
Figure 2: Starting intensities (cid:10) | a m (0) | (cid:11) of normal modes with m = 16 , ..., at a carrierfrequency of 100 Hz excited by a point source set at a depth of 1 km.The ray-based estimates of statistical moments can be obtained in two ways. First, wecan substitute functions X ± m ( r ) computed for different realizations of δc into formula (14)and average a product of mode amplitudes over all the realizations. Second, we can find thesecond moments of X ± m ( r ) by averaging over the realizations and substitute these momentsin Eqs. (17)–(20). Both methods give close results. Therefore below we present only theestimates by the first method.Figure 4 shows the mean values (coherent components) of complex mode amplitudes |h a m i| as functions of range r . Comparison with similar dependencies for non-averaged modeamplitudes presented in Fig. 3 show that the averaging makes the pe-based and ray-basedresults more close.As is seen in Fig. 2, the starting mode intensity | a m (0) | is a rapidly oscillating function ofthe mode number m . For m > the values of | a m (0) | in our example are well approximatedby the WKB relation | a m (0) | = 2 Q m [1 − sin(2 kG m )] . (47)According to this formula, the oscillations are caused by term sin(2 kG m ) . In Eqs. (19) and16 | a | ( k m − / ) | a | ( k m − / ) r (km) rayperaype Figure 3: Amplitudes | a m | of the 16-th (upper panel) and 38-th (lower panel) modes asfunctions of range in a single realization of a fluctuating waveguide. The ray-based estimatepredicted by Eq. (14) (solid line) is compared to result obtained by modal decomposition ofthe sound field computed by the method of wide-angle parabolic equation (dashed line).17
100 200 300 400 5000.20.30.40.50.60.70.80.9 r (km) | < a m > | ( k m − / ) raype | 〈 a 〉 || 〈 a 〉 || 〈 a 〉 | Figure 4: Ray-based (solid lines) and pe-based (dashed lines) estimates of mean mode am-plutudes |h a m i| for m = 16 (upper lines), 38 (middle lines), and 43 (lower lines).18
100 200 300 400 50000.10.20.30.4 r (km) 〈 | a m | 〉 ( k m − ) raypemaster 〈 |a | 〉〈 |a | 〉 〈 |a | 〉 Figure 5: Mean mode intensities (cid:10) | a m | (cid:11) as functions of range for m = 56, 57, and 58. Theray-based results, pe-based results, and solutions of the master equations are shown by thicksolid, thick dashed, and thin solid lines, respectively.(27) the mode coupling manifests itself in the appearance of a weight factor at sin(2 kG m ) which monotonically decreases with range. It means that mean intensities of modes withthe numbers close to m monotonically approaches to Q m . The equalization of mean modeintensities is clearly seen in Fig. 5 where the range dependencies of mean intensities formodes 56, 57, and 58 are shown. Note that the pe-based (thick solid) and ray-based (thickdash) simulations give results close to each other and to the solution of master equations(32) (thin solid). In accord with results of Ref. [8–10] and our result derived in Sec. IV, thesolution of Eqs. (32) gives only a smoothed range dependence of the mean mode intensity (cid:10) | a m | (cid:11) and does not describe its small oscillations.The presence of these oscillations was predicted in Ref. [12]. From the viewpoint ofour ray-based approach they are caused by the fact that the strength of the sound speedfluctuations decreases with depth and therefore the main contributions to integrals (15)come from inhomogeneities located in the vicinities of upper turning points of the moderays. Near-step-like jumps of X ± m ( r ) occur at the mode rays’ upper turning points. The19 .50.6 | < a > | ( k m − / ) < | a | > ( k m − ) z ( k m ) r (km) rayperaypemaster Figure 6: Upper panel: mean amplitude of the 52-nd mode as a function of range. Middlepanel: mean intensity of the 52-nd mode as a function of range. Lower panel: mode rays ofthe 52-nd modes. Thick solid and thick dashed lines in the upper and middle panels showthe ray-based and pe-based results. Thin solid line in the middle panel presents the meanintensity (cid:10) | a | (cid:11) found by solving master equations (32).same is true for the random increment X ( r ) of eikonal of any geometrical ray described byEq. (13). This fact is well known and it underlies the so-called apex approximation [1]. Inour example this effect is most pronounced for steep enough rays with grazing angles at thesound channel axis exceeding 5 ◦ . These rays form modes with numbers m ≥ .According to Eqs. (17) – (19) the jump-like variations of X ± m ( r ) at mode rays’ upperturning points cause jump-like variations of statistical moments at the corresponding dis-tances. This phenomenon is illustrated in Fig. 6 where the range dependencies of the meanamplitude (upper panel) and mean intensity (middle panel) of the 52-nd mode are shown.Dashed vertical lines indicate ranges corresponding to upper turning points of mode rays ofthe 52-nd mode depicted in the lower panel.Formulas (25)–(28) are derived using the approximation of D ( X ± m ( r )) E by the smooth20
100 200 300 400 50000.040.080.120.16 r (km) 〈 | a m | 〉 raype 〈 |a | 〉〈 |a | 〉 〈 |a | 〉 Figure 7: Squared mode intensities for the same normal modes as in Fig. 5. The ray-basedand pe-based results are shown by solid and dashed lines, respectively.function h X m ( r ) i which has no jumps at the upper turning points of mode rays. This explainswhy the master equations (32), whose soulutions (in the high frequency approximation) aregiven by Eq. (27), do not predict the small oscillations of mean mode intensities (cid:10) | a m ( r ) | (cid:11) .Figure 7 presents the mean squared intensities of the same modes as in Fig. 5. It is seenthat the agreement between the ray-based and pe-based estimates for the fourth momentsof the mode amplitudes is less good than for the second moments. C Cross-correlations of normal modes
According to Eqs. (18), the decorrelation of modes m and n is determined by functions Y ± m,n ( r ) = exp (cid:18) − k D(cid:2) X ± m ( r ) − X ± n ( r ) (cid:3) E(cid:19) . (48)We will call Y ± m,n the correlation functions of mode rays. Figure 8 presents the values of Y + mn at 125 (upper panel), 250 (middle panel), and 500 km (lower panel). Functions Y − mn haveclose values (not shown). Let us assume that modes m and n are correlated if Y + m,n > . .21 m m m
10 20 30 40 50 60204060 0.20.40.60.81
125 km250 km500 km
Figure 8: Correlation functions of mode rays Y + mn (Eq. (48)) at ranges 125 km (upper panel),250 km (middle panel), 500 km (lower panel).Then in Fig. 8, we see that at 125 km a typical mode correlates with 20-40 neighboringmodes, at 250 km the number of correlated modes reduces to 10-20, and at 500 km to 3-5.Correlation functions of mode rays Y ± m,n ( r ) monotonically decrease with increasing r and | m − n | . But, according to Eqs. (18) and (26), the dependencies of joint statisticalmoments h a m a ∗ n i on r and | m − n | can be more complicated. In the analysis of the cross-mode coherence, as in Refs. [8–10], we will consider the normalized joint moments of modeamplitudes h b m b ∗ n i , where b m ( r ) = a m ( r ) (cid:10) | a m ( r ) | (cid:11) / . In the upper panel of Fig. 9 we present the range dependencies of cross-mode correlations22or modes with very different starting amplitudes at r = 0 (cf. Fig. 2). Joint moments |h b m b ∗ n i| of such modes within some range intervals may grow with range r and increase withincreasing | m − n | . In our example we see that |h b b ∗ i| grows in the interval from 100 to400 km, and |h b b ∗ i| exceeds |h b b ∗ i| and |h b b ∗ i| .Generally, the joint moment h b m b ∗ n i approaches Y ± m,n ( r ) only at long enough ranges wherethe factors exp (cid:16) − k D ( X ± m ) E(cid:17) become small. But Eq. (26) suggests that there are modeswhose moments h b m b ∗ n i are close to Y ± m,n ( r ) at any r . These are the modes with kG m ≃ π + 2 πJ, (49)where J is an integer. Indeed, for modes m and n satisfying this condition, sin ( k ( G m + G n )) ≃ and cos ( k ( G m − G n )) ≃ . Then from Eqs. (27) and (47), it follows that h b m b ∗ n i ≃ e − k h ( X m − X n ) i ≃ Y ± m,n ( r ) . (50)The lower panel of Fig. 9 present normalized joint moments of such modes. It is seen thatthese moments |h b m b ∗ n i| are reasonably well approximated by the corresponding correlationfunctions Y ± mn ( r ) .Notice that in the case of adiabatic perturbation, formula (50) is applicable for all themodes. Indeed, the adiabaticity of δc ( r, z ) requires that ∆ ζ ≫ D m . If this condition is met (inour theory and numeric example we deal with the inverse inequality), then X + m ( r ) = X − m ( r ) and Eq. (14) translates to a m ( r ) = ϕ m ( z ) e ikX m ( r ) . (51)Equation (50) follows immediately from this formula. Thus, it turns out that even thoughthe adiabatic approximation in our example is not applicable, the normalized cross-modecorrelations h b m b ∗ n i for modes satisfying condition (49) are properly described using a simpleadiabatic formula. It is worthwhile to note, that the starting intensities of modes satisfyingcondition (49) are | a m (0) | ≃ Q m and their mean intensities weakly vary with range (seethe preceding subsection), that is, they behave like the intensities of adiabatic modes. VI Conclusion
In this paper the predictions of our ray-based analytic approach are compared with resultsof full wave numerical simulation. Figures 3–7 and 9 present results of this comparison23 | 〈 b m b n ∗ 〉 | r (km) | 〈 b m b n ∗ 〉 | rayperaypeY +mn | 〈 b b ∗ 〉 || 〈 b b ∗ 〉 || 〈 b b ∗ 〉 || 〈 b b ∗ 〉 || 〈 b b ∗ 〉 | Figure 9: Cross-mode coherences |h b m b ∗ n i| as functions of range. Coherences for mode pairs(25,26), (25,27), and (25,28) are shown in the upper panel; for mode pairs (21,24) and (21,27)are shown in the lower panel. In both panels, the ray-based and pe-based results are shownby thin solid and thin dash lines, respectively. In the lower panel, the correlation functionsof mode rays Y + mn (Eq. (48)) are shown by thick dot-dashed lines.24hich are typical for most modes with numbers m ≤ . The comparison has confirmed theapplicability of our approach for the analysis of mode coupling in a deep water waveguidewith sound speed fluctuations induced by random internal waves. At a frequency of 100 Hzit can be used at ranges of a few hundred kilometers. However, for some modes, especiallyfor those which are weakly excited by the source and have small initial amplitudes, thecoincidence between the ray-based and pe-based estmates may be much worse (not shown).Since our approach is based on the WKB approximation, even at high frequencies itcannot be used for those modes whose turning points are located in the vicinity of thesource depth or in the water bulk near the surface or bottom. We have avoided theseproblems by setting the source at the sound channel axis and restricting our attention tomodes with turning points located well below the surface and above the bottom. The WKBapproximation and, hence, our approach can be applied for modes with turning points onthe waveguide boundaries. But in the present paper such modes were not considered.Numerical simulation has confirmed the prediction of Ref. [12] that the range-dependenciesof mode intensities and other statistical moments are not smooth. They manifest jump-likechanges at ranges corresponding to upper turning points of the mode rays.Important advantage of our approach is its applicability for evaluating the cross-correlationsof mode amplitudes at different frequencies needed for treating pulse propagation. It is clearthat analytical estimates for mode amplitudes at different frequencies are readily derivedalong the same lines as estimates for statistical moments given by Eqs. (17)–(20). But thisissue goes beyond the scope of the present work and it is not broached here.It should be emphasized that formula (14) is derived under assumption that in the pres-ence of perturbation the ray paths do not deviate from their unperturbed positions. Thisassumption is valid only at short enough ranges. In subsection V.B it is shown that the modecoupling causes the equalization of mean mode intensities. According to Eqs. (19) and (27)mean intensities (cid:10) | a m ( r ) | (cid:11) of modes with numbers close to m approach Q m . Numericalsimulation confirms this prediction (see Fig. 5). However, our approach cannot describesubsequent changes of mode intensities with distance. In particular, it cannot be used tostudy establishing the equipartition of energy among the modes in the limit r → ∞ predictedin Ref. [4].In Sec. IV, a numerical result of Refs. [8–10] that the master equations properly describesmoothed range-dependencies of mean mode intensities is explained from the view point ofour ray-based approach. It is shown that formula (27) for the mean mode intensity obtained25y smoothing the range-dependent parameters of Eq. (19) over the cycle of the mode raygives an approximate solution of the master equations valid in the high frequency limit.We assume that the expressions for h a m a ∗ n i given by Eq. (18) represent an approximatesolution of the complete system of M equations for these joint moments derived in theMarkov approximation in Refs. [6, 8]. However, for now, this assumption has not beenverified by direct substitution of Eqs. (18) in this system.In Refs. [9,10], it is shown for a shallow water waveguide, that the analytic expression forthe joint moment h b m b ∗ n i with m = n derived in the adiabatic approximation may be validif the sound speed fluctuations are non-adiabatic. In the present paper, this issue has notbeen studied in detail. However, we hope that our comment on applicability of the adiabaticresults in a non-adiabatic environment made at the end of subsection V.C may contributeto understanding this result of Refs. [9, 10].Finally, note that Eq. (14) for the mode amplitude and the expressions for statistical mo-ments following from this formula can be easily generalized to the case of a range-dependentunperturbed waveguide. This can be done using analytical relations expressing mode ampli-tudes in a range-dependent waveguide through parameters of ray paths [20–22]. Acknowledgment
The work was parially supported by the Program “Fundamentals of acoustic diagnostics ofartificial and natural media” of Physical Sciences Division of Russian Academy of Sciences,Grants No. 13-02-00932 and 13-02-97082 from the Russian Foundation for Basic Research,and Leading Scientific Schools grant N 339.2014.2.
References [1] S.M. Flatte, R. Dashen, W.M. Munk, K.M. Watson, and F. Zakhariasen,
Sound trans-mission through a fluctuating ocean (Cambringe U.P., London, 1979), Chaps. 7, 8, 11.[2] L.M. Brekhovskikh and Yu.P. Lysanov,
Fundamentals of Ocean Acoustics (Springer-Verlag, New York, 2003), Chaps. 6, 10. 263] W. Kohler and G.C. Papanicolaou, “Sound propagation in a randomly inhomogeneousocean,” in
Lecture Notes in Physics. V. 80. Wave propagation and underwater acoustics ,edited by J.B.Keller and J.S.Papadakis (Springer-Verlag, Berlin, 1977), pp. 153-223.[4] L.B. Dozier and F.D. Tappert, “Statistics of normal mode amplitudes in a random ocean.I. Theory ”, J. Acoust. Soc. Am., , 353–365 (1978).[5] L.B. Dozier and F.D. Tappert, “Statistics of normal mode amplitudes in a random ocean.II. Computations”, J. Acoust. Soc. Am., , 533–547 (1978).[6] D.B. Creamer, “Scintillating shallow-water waveguides,” J. Acoust. Soc. Am., , 2825–2838 (1996).[7] A.G. Sazontov, A.L. Matveyev, and N.K. Vdovicheva, “Acoustic coherence in shallowwater: Theory and observation,” IEEE J. Ocean. Eng., , 653–663 (2002).[8] J.A. Colosi and A.K. Morozov, “Statistics of normal mode amplitudes in an oceanwith random sound-speed perturbations: Cross-mode coherence and mean intensity,” J.Acoust. Soc. Am., , 1026–1035 (2009).[9] J.A. Colosi, T.F. Duda, and A.K. Morozov, “Statistics of low-frequency normal-modeamplitudes in an ocean with random sound-speed perturbations: Shallow-water envi-ronments,” J. Acoust. Soc. Am., , Pt. 2, 1026–1035 (2012).[10] K. Raghukumara and J.A. Colosi, “High frequency normal mode statistics in a shallowwater waveguide: The effect of random linear internal waves,” J. Acoust. Soc. Am., ,66-79 (2014).[11] A.L. Virovlyanskii and A.G. Kosterin, “Method of smooth perturbation for the descrip-tion of the fields in multimode waveguides” (in Russian), Akust. Zh. , 599-605 (1987);English transl.: Sov. Phys. Acoust., , 351–354 (1987).[12] A.L. Virovlyanskii, A.G. Kosterin, and A.N. Malakhov, “Mode fluctuations in a canon-ical underwater sound channel” (in Russian), Akust. Zh. , 229-235 (1989); Englishtransl.: Sov. Phys. Acoust., , 138–142 (1989).2713] A.L. Virovlyanskii, “Correlations of modes in a waveguide with large-scale random inho-mogeneities” (in Russian), Izv. Vuzov Radiofizika , 832-838 (1989); English. transl.:Radiophysics and Quantum electronics, , 619–624 (1989).[14] A.L. Virovlyansky, A.G. Kosterin, and A.N. Malakhov, “Fresnel zones for modes andanalysis of field fluctuations in random multimode waveguides,” Waves in Random Me-dia, , 409–481 (1991).[15] A.L. Virovlyansky, V.V. Kurin, N.V. Pronchatov-Rubtsov, and S.I. Simdyankin, “Fres-nel zones for modes,” J. Acoust. Soc. Am., , 163–173 (1997).[16] A.G. Voronovich and V.E. Ostashev, ”Low-frequency sound scattering by internal wavesin the ocean,” J. Acoust. Soc. Am., , 1406–1419 (2006).[17] A.G. Voronovich and V.E. Ostashev, “Coherence function of a sound field in an oceanicwaveguide with horizontally isotropic statistics,” J. Acoust. Soc. Am., textbf125, 99–110,(2009).[18] F.B. Jensen, W.A. Kuperman, M.B. Porter, and H. Schmidt, Computational OceanAcoustics (Springer, New York, 2011), Chaps. 5,6.[19] W. Munk and C. Wunsch, “Ocean acoustic tomography: A scheme for large scale mon-itoring”, Deep-Sea Res., , 123–161 (1979).[20] A.L. Virovlyansky, and G.M. Zaslavsky, “Wave chaos in terms of normal modes,” Phys.Rev. E, , 1656-1668 (1999).[21] A.L. Virovlyansky, A.Yu. Kazarova, and L.Ya. Lyubavin, ”Ray-based description ofnormal mode amplitudes in a range-dependent waveguide,” Wave motion, , 317–334(2005).[22] D. Makarov, S. Prants, A. Virovlyansky, and G. Zaslavsky. Ray and wave chaos in oceanacoustics (Word Scientific, New Jersey, 2010), Chap.3, pp. 150-168.[23] J.A. Colosi and M.G. Brown, “Efficient numerical simulation of stochastic internal-wave-induced sound-speed perturbation field,” J. Acoust. Soc. Am.,103