Table Of Content

Ray dynamics in a long-range acoustic propagation experiment Francisco J. Beron-Vera,1 Michael G. Brown,1 John A. Colosi,2 Steven Tomsovic,3 Anatoly L. Virovlyansky, 4 Michael A. Wolfson,5 George M. Zaslavsky 6 1 Rosenstiel School of Marine and Atmospheric Science, University of Miami, Miami, Florida. 2 Woods Hole Oceanographic Institution, Woods Hole, Massachusetts. 3 3 Department of Physics, Washington State University, Pullman, Washington. 0 4 Institute of Applied Physics, Russian Academy of Science, 6003600 Nizhny Novgorod, Russia. 0 5 Applied Physics Laboratory, University of Washington, Seattle, Washington. 2 6 Courant Institute of Mathematical Sciences, New York University, New York, New York. n (Dated: August 2002) a J Aray-basedwavefielddescriptionisemployedintheinterpretationofmeasurementsmadeduring 2 the November 1994 Acoustic Engineering Test (the AET experiment). In this experiment phase- 2 coded pulse-likesignals with 75 Hzcenterfrequency and37.5 Hzbandwidth weretransmitted near the sound channel axis in the eastern North Pacific Ocean. The resulting acoustic signals were ] recordedonamooredverticalreceivingarrayatarangeof3252km. Inouranalysisbothmesoscale D andinternal-wave-inducedsoundspeedperturbationsaretakenintoaccount. Muchofthisanalysis C exploits results that relate to the subject of ray chaos; these results follow from the Hamiltonian . structureoftherayequations. Wepresentevidencethatalloftheimportantfeaturesofthemeasured n AET wavefields, including their stability, are consistent with a ray-based wavefield description in i l which ray trajectories are predominantly chaotic. n [ 1 I. INTRODUCTION are consistent with measurements and simulations (both v ray- and parabolic-equation-based) at 250 Hz and 1000 6 kmrange[27]-[29]andat133Hzand3700kmrange[30]. The chaotic dynamics of ray trajectories in ocean 2 0 acoustics have been exploredin a number of recent pub- The analysis of Colosi et al. [24] focused on statistical properties of the early resolved AET arrivals. Measure- 1 lications ([1]-[22]), including two recent reviews [21, 22]. 0 In these publications a variety of theoretical results are ments of time spread, travel time variance and proba- 3 presented and illustrated, mostly using idealized mod- bility density functions (PDFs) of peak intensity were 0 els of the ocean sound channel. In contrast, in the presentedand comparedto theoreticalpredictions based n/ present study a ray-based wavefield description is em- on a path integral formulation as described in Ref. [31]. i ployed, in conjunction with a realistic environmental Travel time variance was well predicted by the theory, l n model, to interpret a set of underwater acoustic mea- but time spreads were two to three orders of magnitude v: surements. For this purpose we use both oceanographic smaller than theoretical predictions, and peak intensity i and acoustic measurements collected during the Novem- PDFswerecloseto lognormal,inmarkedcontrastto the X ber 1994 Acoustic Engineering Test (the AET experi- predictedexponentialPDF.The surprisingconclusionof r ment). In this experiment phase-coded pulse-like signals this analysis was that the measured AET pulse statis- a with75Hzcenterfrequencyand37.5Hzbandwidthwere tics suggested propagation in or near the unsaturated transmitted near the sound channel axis in the eastern regime (characterized physically by the absence of mi- NorthPacificOcean. The resultingacousticsignalswere cromultipaths and mathematically by use of a perturba- recordedon a moored vertical receiving array at a range tion analysis based on the Rytov approximation), while of 3252 km. We present evidence that all of the impor- thetheorypredictedpropagationinthesaturatedregime tantfeaturesofthemeasuredAETwavefieldsareconsis- (characterized by the presence of a large number of mi- tent with a ray-based wavefield description in which ray cromultipaths whose phases are random). trajectories are predominantly chaotic. In Ref. [25] it was shown that in the finale region of the measured AET wavefields, where no timefronts are The AET measurements have previously been ana- resolved,theintensityPDFisclosetothefullysaturated lyzed in Refs. [23]-[25]. In the analysis of Worcester exponential distribution. This result is not unexpected et al. [23] it was shown that the early AET arrivals because the finale regionis characterizedby strong scat- could be temporally resolved, were stable over the du- tering of both rays [14] and modes [32]. ration of the experiment and could be identified with specific ray paths. These features of the AET observa- The work reported here seeks to elucidate both the tions are consistent with those of other long-range data physics underlying the AET measurements, and the sets [26,27]. TheAET arrivalsin the last1 to 2sec (the causes of the successes and failures of the analyses re- arrival finale) could not be temporally resolved or iden- portedinRefs. [23]-[25]. Weshowthataray-basedanal- tified with specific ray paths, however. Also, significant ysis, in which ideas associated with ray chaos play an verticalscatteringofacousticenergywasobservedinthe important role, can account for the stability of the early arrival finale. These features of the AET arrival finale arrivals, their small time spreads, the associated near- 2 lognormal PDF for intensity peaks, the vertical scatter- assume propagation in the vertical plane with cartesian ing of acoustic energy in the reception finale, and the coordinates z and r representing depth and range, re- near-exponential intensity PDF in the reception finale. spectively. Ray methods can be introduced when the Part of the success of this interpretation comes about wavelength 2πc/σ of all waves of interest is small com- due to the presence of a mixture of stable and unstable pared to the length scales that characterize variations in ray trajectories. Also, an explanation is given for differ- c. Substituting the geometric ansatz, ing intensity statistics in the early and late portions of the measured wavefields, and the fairly rapid transition u(z,r;σ)= A (z,r)eiσTj(z,r), (2) j between these regimes. j X The remainder of this paper is organized as follows. In the next section we describe the AET environment representing a sum of locally plane waves, into the and the most important features of the measured and Helmholtz equation gives, after collecting terms in de- simulated wavefields. In the simulated wavefields both scending powers of σ, the eikonal equation, measuredmesoscaleandsimulatedinternal-wave-induced ( T)2 =c−2, (3) sound speed perturbations are taken into account. It is ∇ shownthatevenintheabsenceofinternalwaves,thelate and the transport equation, arrivalsarenotresolvedandtheassociatedraypathsare chaotic. In section II we present results that relate to (A2 T)=0. (4) the structure of the early portion of the measured time- ∇ ∇ front and its stability. The micromultipathing process is Fornotationalsimplicity,thesubscriptsonAandT have discussed,andaquantitativeexplanationfortheremark- been dropped in (3) and (4). ably small time spreads ofthe early ray arrivalsis given. Thesolutiontotheeikonalequation(3)canbereduced Section III is concerned with intensity statistics. An ex- tosolvingasetofrayequations. Ifaone-wayformulation planation is givenfor the cause of the different wavefield (see, e.g.,Ref. [22])is invoked,with r playingthe roleof intensity statistics in the early and late portions of the the independent variable, the ray equations are arrival pattern. In the final section we summarize and discuss our results. dz ∂H dp ∂H = , = , (5) dr ∂p dr − ∂z II. MEASURED AND SIMULATED and WAVEFIELDS dT dz =L=p H, (6) dr dr − Figure1showsameasuredAETwavefieldinthetime- depth plane and ray-based simulations of such a wave- where fieldwithandwithoutinternalwaves. Figure2isacom- plementary display that shows measured and simulated H(p,z,r)= c−2(z,r) p2. (7) − − wavefieldsafter plane wavebeamforming. In this section p These equations constitute a nonautonomous Hamilto- we describe the most important qualitative features of nian system with one degree of freedom; z and p are measured and simulated AET wavefields. We begin by canonicallyconjugatepositionandmomentumvariables, reviewingfundamentalraytheoreticalresultssothatour r is the time-like variable, H is the Hamiltonian and the raysimulationscanbefullyunderstood. Also,thismate- travel time T corresponds to the classical action. Note rialprovidesafoundationforsomeofthelaterdiscussion. also that p p = ∂T/∂z and p = ∂T/∂r = H are We then briefly describe the AET environment and the z r ≡ − the z- and r-components, respectively, of the ray slow- signal processing steps that were performed to produce ness vector. The ray angle relative to the horizontal ϕ themeasuredwavefieldshowninFig. 1. Next,ourtreat- satisfiesdz/dr=tanϕ,or,equivalently,cp=sinϕ. Fora mentofinternal-wave-inducedsoundspeedperturbations point source rayscan be labeled by their initial slowness is briefly described. Finally, we return to describing and p . The solution to the ray equations is then z(r;p ), comparing qualitative features of the measured and sim- 0 0 p(r;p ), T(r;p ). The terms in (2) correspond to eigen- ulated wavefields. 0 0 rays; these satisfy z(r ;p ) = z where the depth and Fixed-frequency (cw) acoustic wavefields satisfy the R 0 R range of the receiver are z and r , respectively. The Helmholtz equation, R R transport equation (4) can be reduced to a statement of 2u+σ2c−2(z,r)u=0, (1) constancy of energy flux in ray tubes. The solution, for ∇ the j-th eigenray, can be written where σ = 2πf is the angular frequency of the wave- tfirealdnsiaenndtwc(azv,erfi)elidsstihsestsroauinghdtfsopreweadr.d(uTsihnegeFxotuerniseironsytno- Aj(z,r)=A0|H/q21|1j/2e−iµjπ2 (8) thesis. It is not necessary to consider this complication where H, the matrix element q and the Maslov index 21 topresentthemostimportantresultsneededbelow.) We µ are evaluated at the endpoint of the ray, and A is a 0 3 constant. The stability matrix this coherent processing yields improved signal-to-noise is limited by ocean fluctuations, mostly due to internal ∂p ∂p waves. The statistics of the acoustic fluctuations were Q= q11 q12 = ∂p0 z0 ∂z0 p0 (9) showninRef. [25]nottobeadverselyaffectedbythe12.7 (cid:18)q21 q22 (cid:19)  ∂∂pz0(cid:12)(cid:12)(cid:12)z0 ∂∂zz0(cid:12)(cid:12)(cid:12)p0  minute coherentintegration. Nearly concurrentwith the  (cid:12) (cid:12)  week-long experiment, temperature profiles in the upper (cid:12) (cid:12) quantifies ray spreading, (cid:12) (cid:12) 700 m were measured with XBT’s every 30 km along the transmission path. An objective mapping technique δp δp =Q 0 . (10) was applied [23] to these measurements to produce a δz δz 0 mapofthe backgroundsoundspeedstructure (including (cid:18) (cid:19) (cid:18) (cid:19) mesoscale structure) along the transmission path. The Elements of this matrix evolve according to sampling interval in range of the resulting sound speed d field is 30 km so that no structure with horizontalwave- Q=KQ (11) lengths less than 60 km is present in this field. More dr details on the AET experiment, environment and pro- where Q at r =0 is the identity matrix, and cessing of the acoustic signals can be found in Ref. [23]. ∂2H ∂2H Internal-wave-induced sound speed perturbations are K = −∂z∂p −∂z2 . (12) taken into account in most of our simulations. These ∂2H ∂2H ∂p2 ∂z∂p ! are assumed to satisfy the relationship δc = (∂c/∂z)pζ and the statistics of ζ, the internal-wave-induced ver- At caustics q vanishes and, for waves propagating in 21 tical displacement of a water parcel, were assumed to twospacedimensions,the Maslovindex advancesbyone be described by the empirical Garrett-Munk spectrum unit. (see, e.g., Ref. [33]). The potential sound speed gradi- In ocean environments with realistic range- ent (∂c/∂z) , the buoyancy frequency N, and the sound p dependence, including the AET environment, ray speed c were estimated directly from hydrographic mea- trajectories are predominantly chaotic. Chaotic rays surements using a procedure similar to that described in diverge fromneighboring rays at anexponential rate, on the appendix of Ref. [30]; it was found that a good ap- average, and are characterized by a positive Lyapunov proximation for the AET environment is (∂c/∂z) /c p exponent, (1.25 s2/m)N2. The vertical displacement ζ(r,z) wa≈s computed using equation 19 of Ref. [34] with the vari- lim 1 lim (r) νL = r r (0) 0 lnD(0) . (13) able x replaced by r, and y = t = 0. Physically this →∞ (cid:18) D → D (cid:19) correspondstoafrozenverticalsliceofthe internalwave field that includes the influence of transversely propa- Here (r) isameasureofthe separationbetweenraysat D gating internal wave modes. A Fourier method, with ranger; suitablechoicesfor areanyofthe fourmatrix D ∆k = 2π/3276.8 km and 0 < k 2π/1 km, was used elementsofQorthetraceofQ. Thechaoticmotionofray x x ≤ to numerically generate sound speed perturbation fields. trajectories leads to some limitations on predictability. (Ray calculations were also performed using an internal This will be discussed in more detail below. Additional wave field in which the hard upper bound on k was background material relating to this topic can be found x replaced by an exponential damping of the large k en- in [21, 22] and references therein. x ergy. The qualitative ray behavior described below did We digress now to briefly describe the AET experi- not depend on the manner in which the large k cut-off ment and relevant signal processing. In this experiment x was treated.) A mode number cut-off of j = 30 was a source, suspended at 652 m depth from the floating max usedinallthesimulationsshown. Contributionsfromap- instrument platform R/P FLIP in deep water west of proximately 230 (215/5 215/5 30 where the factors of SanDiego,transmittedsequencesofaphase-codedsignal × × 1/5arepresentbecausetheassumedwavenumbercut-off whose autocorrelation is pulse-like (with a pulse width wasone-fifthoftheassumedNyquistwavenumber)inter- of approximately 27 ms). The source center frequency nal wave modes were included in our simulated internal- was 75 Hz and the bandwidth was approximately 37.5 wave-inducedsound speed perturbation fields. Measure- Hz. The resulting acoustic signals were recorded east of ments of sound speed variance made during the experi- Hawaiiatarangeof3252.38kmonamoored20element ment (see Fig. 3) indicate that the internal wave energy receivingarraybetweenthedepthsof900mand1600m. strength parameter E was close to E , the nominal After correcting for the motion of the source and receiv- GM Garrett-Munk value. Simulations were performed using ingarray,including removalofassociatedDopplershifts, both E = 0.5E and E = 1.0E . All of the simu- the receivedacousticsignalswerecross-correlatedwitha GM GM lations shown in this paper use E = 0.5E . The dif- replica of the transmitted signal (to achieve the desired GM ference between these simulations and those performed pulse compression) and complex demodulated. To im- using E =1.0E will be discussed where appropriate. prove the signal-to-noise ratio, 28 consecutive processed GM receptions, extending over a total duration of 12.7 min- With the foregoing discussion as background, we are utes, were coherently added. The duration over which prepared to discuss our ray simulations in the AET en- 4 vironment. These are shown in Figs. 1, 2, 4 through field. Our simulations show that the stability exponents 11, and 14. The points plotted in Figs. 1, 2 and 4 were of rays with axial angles of approximately 5o or less ap- computedbyintegratingthecoupledrayequations(5,6) proximately double when j increases from 30 to 100, max from the source to the range of the receiving array. Sev- while steeper rays show very little sensitivity to j . max eral fixed and variable step-size integration algorithms In contrast, stability exponents of rays with axial angles were tested. In the presence of internal waves, conver- fromapproximately8oto15oapproximatelydoublewhen gence using a fixed step fourth order Runge Kutta in- k increases from 2π/1000 m to 2π/250 m, while the max tegrator required that the step size not exceed approxi- flatter rays show little sensitivity to k . With these max mately one meter. All ray tracing calculations presented comments in mind, one should not attach too much sig- inthispaperwereperformedusingdoubleprecisionarith- nificancetothe detailsofFig. 5. We note,however,that metic (64 bit floating point wordsize). The wavefield in- rayswith axialangles of approximately5o or less consis- tensity is approximately proportional to the density of tentlyhavehigherstabilityexponentsthanraysinthe5o rays (dots) that are plotted in Figs. 1 and 2. A some- to 10o band. whatmoredifficultcalculationinwhichthecontributions Note that what appears to be temporally resolved ar- tothewavefieldfrommanyraysarecoherentlyaddedwill rivals in the measured wavefieldcorrespond in our simu- be discussed below. lationstocontributionsfromanexponentiallylargenum- Inthe presenceofinternalwaves,rayswithlaunchan- berofraypaths. (Onlyasmallfractionofthetotalnum- gles between approximately 5o form the diffuse finale ber are included in our simulations, however, because of ± regionofthearrivalpatternarrivingafterapproximately the relativesparsenessofourinitialsetofrays.) The ob- 2196 sec. Steeper rays contribute to the earlier, mostly servation that the travel times of chaotic ray paths may resolved arrivals. It is seen in Fig. 1 that our ray sim- clusterandberelativelystablewasfirstmadeinRef. [8]. ulation in the presence of internal waves underestimates Simmenetal. [14]havepreviouslyproducedplotssim- the verticalspreadof the finale region. Simulations with ilar to our Figs. 1 and 4 for ray motion in a deep ocean a stronger internal wave field, E =1.0EGM, give a bet- model consisting of a background sound speed struc- ter fit to the vertical spread of energy (but a poorer fit ture very similarto oursonwhich internal-wave-induced to other wavefield features, as described below). Note, soundspeedperturbationsweresuperimposed. Although however,thatdiffractiveeffects, missinginourraysimu- ourresults are similarin many respects,it is noteworthy lations, will contribute to enhanced broadening in depth thatthe rightpanelofourFig. 4showsmorechaoticbe- of the arrival pattern. A striking feature of Figs. 1, 2 havior than is present in the corresponding plot in Ref. and 4 is the contrast between the seemingly structure- [14]. Thisdifferencepersistsifour3252.38kmrangesim- less distribution of steep (large p) rays in (p,z) seen in ulations are replaced by simulations at the same range | | Fig. 4wheninternalwavesarepresentandtherelatively (1000 km) that was used in Ref. [14]. Our rays, espe- structureddistributionofthesameraysinthetime-depth ciallythesteeperrays,aremorechaoticthanthoseinRef. (Fig. 1) and time-angle (Fig. 2) plots. [14]. We believe that this difference is primarily due to The cause of the complexity seen in Fig. 4 is ray the increased complexity of our internal wave field over chaos [21, 22]; most ray trajectories diverge from ini- that of Simmen et al. [14] who included contributions tially infinitesimally perturbed rays at an exponential fromten internal wavemodes in their simulated internal rate. 48,000 rays with launch angles between 12o were wave fields. ± traced to produce Figs. 1, 2 and 4 in the presence of in- Figures 6 and 7 show plots of ray depth at the AET ternalwaves. This fan of raysis far too sparse to resolve range vs. launch angle. In Fig. 6 the same rays (exclud- whatshouldbeanunbrokensmoothcurve–aLagrangian ing those with negative launch angles)that were used to manifold – which does not intersect itself in phase space produceFigs. 1,2and4areplotted. InFig. 7twosmall (Fig. 4). Under chaotic conditions the separation be- angularbandsofthe upperpanelofFig. 6areblownup; tween neighboring rays grows exponentially, on average, theraydensityisfourtimesgreaterthanthatusedinFig. in range. The complexity of the Lagrangian manifold 6. Fig. 6 strongly suggests that, even in the absence of growsatthesameexponentialrate. TheLyapunovexpo- internal waves, the near-axial rays are chaotic. Not sur- nentis the recriprocalofthe e-foldingdistance (see, e.g., prisingly, in the presence of internal waves, steeper rays [22]). Finite range numerical estimates of Lyapunov ex- arealsopredominantlychaotic. Note,however,thatFigs. ponents (hereafter referredto as stability exponents)are 5, 6 and 7 show that the near-axial rays are evidently shown as a function of launch angle in Fig. 5. It is seen more chaotic than the steeper rays. An explanation for thatinthisenvironmenttheflatrays(ϕ0 <∼5o)havesta- this difference in behavior will be given in section III. | | bility exponents of approximately 1/(100 km), while the Inthissectionwehavebrieflydescribedandcompared steeper rays(6o <∼ ϕ0 <∼11o) havestability exponentsof thegrossfeaturesofmeasuredandsimulatedAETwave- | | approximately 1/(300 km). It follows, for example, that fields. Raytrajectoriesinoursimulatedwavefieldsinthe the complexity of the flat ray portion of the Lagrangian presence ofinternalwavesare predominantly chaotic. In manifold grows approximately like exp(r/100km). spite of this, Figs. 1 and 2 show that many features of Not surprisingly, stability exponents show some sen- our simulated wavefields are stable and in good quali- sitivity to parameters that describe the internal wave tative agreement with the observations. Predicted and 5 measured spreads of acoustic energy in time, depth and seen in Fig. 1. angle are generally in good agreement, both in the early Animportant(andsurprising)featureofsetsofmicro- and late portions of the arrival pattern. We shall not multipaths is that they are highly nonlocal in the sense discuss further the spreadsofenergyin depth andangle, that interspersed (i.e., having intermediate launch an- excepttonotethatFigs. 1and2showthatthesespreads gles)amongagroupofmicromultipathsarerayswithM can be accounted for using ray methods in the presence values that differ by severalunits. This behavior is illus- ofrealistic(including internal-wave-inducedsoundspeed trated in Fig. 10. In view of the observationthat, in the perturbations) ocean structure. In the two sections that presence of internal waves in the AET environment, ray follow we shall consider in more detail the time spreads motion is strongly chaotic and the function M(ϕ ) has o of the early ray arrivals, and intensity statistics in both local oscillations of several units, it is remarkable that the early and late portions of the arrival pattern. the early portion of the timefront (see Fig. 1) is not de- stroyed. The nonlocality of a set of micromultipaths in the range-depth plane is seen in Fig. 8. Note that there III. CHAOS, MICROMULTIPATHS AND are significant differences in the upper and lower turn- TIMEFRONT STABILITY ingdepthsoftheplottedmicromultipathsandthatthese raysarespreadin rangeby a significantfractionof a ray cycle distance. In this section we consider in some detail eigenrays, timefront stability and time spreads. We focus our at- Theresultofperformingaray-basedsynthesisofwhat tention on the early branches of the timefront, where appears in the measurements to be an isolated arrival measured time spreads were, surprisingly, only approxi- is shown in Fig. 11. This involves finding a complete mately 2 ms, in sharp contrast to the theoretical predic- set of micromultipaths, including their travel times and tion of approximately 1 sec [24]. We focus our attention Maslov indices, and coherently summing their contribu- ontheearlyarrivaltimespreads,primarilybecausequan- tions. Unfortunately,owingtothepredominantlychaotic titative estimates of the unresolved late arrivals are not motion of ray trajectories in the AET environment in available. the presence of internal waves, it is extremely difficult to find a complete set of micromultipaths. Indeed, the Eigenrays at a fixed depth z correspond to the inter- constituent ray arrivals used in the Fig. 11 synthesis do sections of a horizontal line at depth z with the curve notconstituteacompletesetofmicromultipaths. Thisis z(ϕ ); these intersections are the roots of the equation o lessimportantthanmightbe expectedbecausestandard z z(ϕ ) = 0. Although the discrete samples of z(ϕ ) o o − eigenrayfindingtechniqueseasilylocatethestrongestmi- plotted in Figs. 6 and 7 are too sparse to reveal what cromultipaths; the missing micromultipaths in the Fig. should be a smooth curve, it is evident the number of 11 synthesis are highly chaotic raysthat have very small eigenraysat almostalldepths is verylarge;because rays amplitudes. Aninterestingresultofthissynthesisisthat are predominantly chaotic, this number grows exponen- the micromultipath-induced time spread in the synthe- tially in range. In principle, eigenrays can be found us- sized pulse is only about 1 ms, which is in very good ingacombinationofinterpolationanditeration,starting agreement with the AET measurements [24, 25]. Note with a set of discrete samples of z(ϕ ). In practice, this o that the total time spread among the micromultipaths procedure reliably finds only those eigenrays in regions found – about 11 ms – is much greater than the spread where z(ϕ ) has relatively little structure. These rays o of the synthesized pulse. This difference arises because have the highestintensity and arethe leastchaotic. It is thetimespreadofthedominantmicromultipathsismuch seen in Fig. 1 that the early (steeper) ray arrivals come smaller than the total micromultipath time spread. in clusters with small travel time spreads. The set of eigenrays within such a cluster is often referred to as a Simulations (not shown) using E = 1.0 EGM result setofmicromultipaths. Anincompletesetofmicromulti- in synthesized steep ray pulses that are spread in time paths,foundusingtheprocedurejustdescribed,isshown much more than is shown in Fig. 11. With the stronger inFig. 8. Aninterestingandimportantfeatureofmicro- internalwavefield,the totalmicromultipathtimespread multipaths in our simulations is that all of the rays that increases by approximately a factor of two – to about make up a set of micromultipaths have the same iden- 25 ms. More importantly, however, the dominant mi- tifier M. Here is the sign of the launch angle and cromultipaths have time spreads comparable to the to- M is t±he number o±f ray turning points (where p changes tal micromultipath time spread, leading to synthesized sign) following a ray. In Fig. 9 rays with two fixed val- pulses that are spread by 5 to 10 ms, rather than 1 to ues of ray identifier (+137 and +151) are plotted in the 2 ms. Thus, measured early arrival time spreads are in time-depth plane. In both cases the points plotted are better agreement with E = 0.5EGM simulations than a subset of those plotted in Fig. 1. It is seen that the with E = 1.0EGM simulations. The question of which relatively steep +137 rays have a small time spread and value of the internal wave strength E in our simulations form one of the branches of the timefront seen in Fig. 1 gives the best fit to the observations will be revisited in (similar behavior was noted in Ref. [14]), while the rela- the next section. tivelyflat+151rayshaveamuchlargertime spreadand An interesting and unexpected feature of our simula- fall within the diffuse finale region of the arrival pattern tionsisthattheMaslovindexwasconsistentlylowerthan 6 the number of turns M made by the same ray. The dif- apex scattering events; the other ray is an unscattered ference was typically three units in the E = 0.5 E ray with the same launch angle as the scattered ray. GM simulations and four units in the E = 1.0E simula- Our strategy to building the eigenray constraint into GM tions. We have no explanation for this behavior. anestimateoftraveltimespreadscannowbestated. We Sofarinthissectionwehaveseenthatassociatedwith consider a scattered and an unperturbed (that is, in the the chaotic motionof ray trajectoriesis extensive micro- absence of internal-wave-induced scattering events) ray multipathing, and that the micromultipathing process is withthesamerayidentifier;theraysshownschematically highly nonlocal. The many micromultipaths that add in Fig. 13 have identifier +3. In addition to having the at the receiver to produce what appears to be a single same ray identifier, the scattered and unperturbed rays arrival may sample the ocean very differently. Surpris- are assumed to start from the same point and end at ingly, on the early arrival branches the highly nonlocal the same depth, but they will generally end at different micromultipathing process causes only very small time ranges. Weassumethatthescatteredrayhastraveltime spreads and does not lead to a mixing of ray identifiers. Tr and is one of many scatteredeigenrayswith the same To complete this picture it is necessary to explain why rayidentifieratranger. WeestimatethetraveltimeT(r) the time spreads of the early (steep) ray arrivals are so ofanunperturbedeigenrayatrangerwhoserayidentifier small. We shall now provide such an explanation. Re- is the same as that of the scattered eigenray. It will be call first, however, that the measured time spread of the shown that δT(r) = Tr T(r) vanishes to first order if − early arrivalsis only approximately 2 ms [24], in marked the apex approximation is exploited. Because the same contrastto the theoretical prediction[24, 31] of approxi- result applies to all of the scattered eigenrays, it follows mately 1 sec. Thus our quantitative explanation for the that the time spreadvanishes to first order,independent smallearlyarrivaltimespreadsservesbothtogiveinsight of r, if the apex approximation is exploited. Note that into the underlyingphysics,andprovidesanexplanation this result is not expected to apply, even approximately, for a critically importnat feature of the measured AET tothelatenear-axialrayswheretheapexapproximation wavefields. is known to fail. The travel time T(r) of the eigenray in the unper- Computing time spreads is conceptually straightfor- turbed ocean can, of course, be computed numerically ward using ray methods. In a known environment one but this provides little help in deriving an analytical es- finds a large ensemble of rays with the same identifier timate ofδT. Instead,considera rayinthe unperturbed that solve the same two point boundary value problem; oceanwiththe samelaunchangleasoneofthe scattered the spread in travel times over many ensembles of such eigenrays, as shown in Fig. 13. In general, the range rays, each ensemble computed using in a different real- of this ray at the receiver depth, after making the ap- ization of the environment, is the desired quantity. The propriate number of turns, will be r = r. But T(r) constraint that all of the computed rays have the same o 6 can be estimated from T(r ). This follows from the ob- fixed endpoints complicates this calculation. In the fol- o servation that ray travel time (assuming a point source lowing we exploit an approximate form of the eigenray initial condition, for instance) is a continuously differen- constraint that simplifies the calculation; the approxi- tible functionofz andr with T =pwherepis the ray mate form of the eigenray constraint is first order accu- ∇ slowness vector. Because of this property T(z,r) can be rate in a sense to be described below. expandedin a Taylorseries. Ifwe considerrayswith the Before describing the manner in which the eigenray samerayidentifierandfixthe finalraydepthtocoincide constraint is imposed, it is useful to describe the physi- with the receiver depth, then calsettinginwhichtheconstraintwillbeapplied. Figure 12shows H vsr followingamoderatelysteepray(axial T(r) T(r )+p (r )(r r ) (14) angleapp|rox|imately11o)inacanonicalenvironment[35] ≈ o r o − o where p = H is the r-component of the ray slow- on which the previously described (based on the mea- r − ness vector. (More generally, T(z,r) consists of multiple sured AET N(z) profile) internal-wave-induced sound smooth sheets that are joined at cusped ridges; in spite speed perturbation field was superimposed. It is seen ofthis complication,the propertyofcontinuousdifferen- thatH(r)consistsofasequenceofessentiallyconstantH tiabilityismaintained,eveninahighlystructuredocean. segments, separatedby near-step-likejumps. The jumps Forourpurposes,weneedonlyconsideronesheetofthis occur at the ray’s upper turning points. A simple model multi-valued function. Note also that our application of that captures this behavior – the so-called apex approx- (14) in the backgroundocean makesthis expressionpar- imation [31] – assumes that the transition regions can ticularly simple to use.) It follows from (14) that be approximated as step functions. (The validity of the apexapproximationislinkedtotheanisotropyandinho- δT(r)=T T(r) T T(r ) p (r )(r r ) (15) r r o r o o mogeneity of oceanic internal waves; details of the back- − ≈ − − − ground sound speed profile are not important. For this is, correct to first order, the travel time difference at reason we have chosen to illustrate this effect using a range r between eigenrays with and without internal background canonical profile rather than the AET envi- waves. ronment.) Figure 13 shows a schematic diagram of two ToevaluateδT(r)weshallassumethatthebackgound rays. One ray has undergone two internal-wave-induced environmentis range-independent. Also, consistent with 7 the use of the apex approximation, we may treat the Several comments concerning the preceeding calcula- perturbed environment as piecewise range-independent. tion are noteworthy. First, we note that although it was Within each range-independent segment it is advanta- assumed that the background sound speed structure is geousto make use ofthe action-anglevariableformalism range-independent, the preceeding argument also holds (see. e.g., [22] or [36]). In terms of the action-anglevari- in the presence of slow background range-dependence, ables (I,θ), H˜ =H˜(I) and the ray equations are i.e.,withstructurewhosehorizontalscalesarelargerela- tivetoatypicalraydoublelooplength. Adiabaticinvari- dθ ∂H˜ anceinsuchenvironmentsguaranteesthatwhileH isnot = =ω(I), (16) dr ∂I constant following rays between apex scattering events, I isnearlyconstant. Second,afternrandomkicksI I 1 0 − in(19)and(20)isreplacedbyI I √n(∆I) ,and dI ∂H˜ n− 0 ≈ rms = =0, (17) themagnitudeof(19)underAETconditionsisO(1sec). dr − ∂θ This is an example of a traveltime spread estimate that and fails to enforce the eigenray constraint. This calculation shows that the difference between travel time estimates dT dθ =I H˜(I)=Iω(I) H˜(I), (18) that do and do not enforce the eigenray constraint can dr dr − − be quite significant. (The constrainedestimate is refined where the physical interpretation of H˜ as p is main- below.) Third, it should be emphasized that Eqs. (19) r − and(20)apply(approximately)onlytothesteepraysbe- tained. The angle variable can be defined to be zero cause the apex approximation applies only to the steep at the upper turning depth of a ray and increases by rays. We have not addressed the time spreads of near- 2π over a ray cycle (double loop). It follows that the axial rays. Note, however, that Fig. 9 shows that time frequency ω(I) = 2π/R(I) where R(I) is a ray cycle spreadsare largerfor flatter rays. Fourth, the abovecal- distance. Integrating the ray equations over a com- culation shows that, to lowest order in ∆I, there is no plete ray cycle then gives R(I) = 2π/ω(I) and T(I) = 2π(I H˜(I)/ω(I)) with I constant following each ray. internal-wave-scattering-induced travel time bias if the − apexapproximationisstrictlyapplied. Andfifth,thear- In the apex approximation I jumps discontinuously at gumentsleadingtoEqs. (19)and(20)applywhetherthe each ray’s upper turning depth; for the perturbed ray scattered rays are chaotic or not. R(I +∆I) 2π/ω(I) 2π∆Iω′(I)/(ω(I))2 and T(I + ≈ − A relaxed form of the apex approximation in which ∆I) 2π(I H˜(I)/ω(I)) + 2π∆IH˜(I)ω′(I)/(ω(I))2 ≈ − the action jump transition region has width ∆θ gives a where ω′(I) = dω(I)/dI. Note that, like these expres- nonzero travel time spread estimate. In the transition sions, Eqs. 14 and 15 are first order accurate in ∆I. region, taken for convenience to be 0 θ ∆θ, one Consider again Eq. 15 and Fig. 13 with I equal to may choose a Hamiltonian of the form≤H˜ =≤H˜(I sθ) I , I and I in the left center and right ray segments, − 0 1 2 where s = ∆I/∆θ, and H˜ = H˜(I) elsewhere. It follows respectively. The left segment gives no contribution to that in the transition region I(θ) = I + sθ and I = δT as the ray has not yet been perturbed. In the center 0 constantelsewhere. Asimplegeneralizationofthe above ray segment calculation then gives for a complete ray cycle H˜(I )ω′(I ) 0 0 T T(r ) 2π(I I ) (19) T T(r ) r− o ≈ 1− 0 (ω(I0))2 r− o ≈ H˜(I )ω′(I ) ∆θ 0 0 and (I1−I0)"(2π−∆θ) (ω(I0))2 + 2 # (21) ω′(I ) p (r )(r r ) H˜(I ) 2π(I I ) 0 . (20) − r o − o ≈ 0 − 1− 0 (ω(I ))2 and (cid:20) 0 (cid:21) These terms are seen to cancel. Note that if additional p (r )(r r ) r o o − − ≈ complete cycle raysegmentsare addedto the center sec- ω′(I ) tionoftheray,Eqs. 19and20areunchangedexceptthat H˜(I ) (2π ∆θ)(I I ) 0 ,(22) I1−I0 = ∆I1 is replaced by ni=1∆Ii = In−I0; again 0 (cid:20)− − 1− 0 ω(I0)2(cid:21) the two terms cancel. In the final (incomplete cycle) ray so the sum (recall [15]) is P segmentthe differencebetweenthe termsT T(r )and r o pr(ro)(r ro)canbeshowntobeO((∆I)2);t−hetermsdo ∆θ not exac−tly cancel because the θ-values of the perturbed δT (I1 I0). (23) ≈ 2 − and unperturbed rays are generally not identical at the receiver depth. Thus, to first order in ∆I, δT = 0, in- Itshouldbeemphasizedthatthisexpressionisfirstorder dependent of range, if the apex approximation is valid. accurate in ∆I =I I , but that no assumption about 1 0 − This simple calculation provides an explanation of why, the smallness of ∆θ has been made. As noted above, in spite of extensive ray chaos, the time spreads of the incomplete cycle end segment pieces give O((∆I)2) cor- early AET ray arrivals are quite small. rections to δT. If ∆θ has approximately the same value 8 at each upper turn, then one has after n upper turns, The latter regime is characterized by a lognormal inten- correct to O(∆I), sity PDF. This fix is conceptually problematic inasmuch as, in this theory, the unsaturated regime is character- n ∆θ ∆θ ized by the absence of micromultipaths, which seriously δT (I I ) = (I I ) ≈ 2 i− i−1 2 n− 0 conflicts with the numericalsimulationspresentedin the Xi=1 previous section where the number of micromultipaths ∆θ √n(∆I) . (24) is very large. In this section we provide self-consistent rms ≈ 2 explanations for both the late arrival exponential distribution and the early arrivallognormaldistribution. The Consistent with the numerical simulations shown in Fig. challenge is to reconcile the early arrival near-lognormal 12, ∆θ 0.8 radians and ∆I 4 ms. (We, and in- ≈ ≈ intensity PDF with the presence of a large number of dependently F. Henyey [personal communication], have micromultipaths. Some of the arguments presented are confirmed that these estimates also apply under AET- heuristic,andsomebuildontheresultsofnumericalsim- like conditons.) With these numbers and √n = 8, ap- ulations. A complete theoreticalunderstanding of inten- propriate for AET, (24) gives a time spread estimate sity statistics has proven difficult. of approximately 13 ms, in approximate agreement with the numerical simulations shown in Fig. 11. Note that Our approach to describing wavefield intensity statis- (24) does not preclude a travel time bias. A cautionary ticsbuildsonthesemiclassicalconstructiondescribedby remark concerning the use of (24) is that Virovlyansky Eq. (2). At a fixed location it is seen that the wave- [37] has pointed out that, owing to secular growth, the field amplitude distribution is determined by the distri- O((∆I)2) contribution to δT may dominate the O(∆I) bution of ray amplitudes and their relative phases. Note contribution at long range. that both travel times and Maslov indices influence the Figure9showsthatintheAETenvironmentsimulated phasesofrayarrivals. Complexitiesassociatedwithtran- near-axial ray time spreads are greater than simulated sient wavefields and caustic corrections will be discussed steep ray time spreads. We do not fully understand this below. behavior. Apossibleexplanationforthisbehavioristhat An important observation is that in the AET envi- timespreadsincreaseasraysbecomeincreasinglyflatow- ronment, including internal-wave-induced sound speed ingtothebreakdownoftheapexapproximation. Butone perturbations, simulated geometric amplitudes of both would expect that this trend should be offset, in part or steep and flat ray arrivals approximately fit lognormal whole,bytherelativesmallnessofinternal-wave-induced PDFs. This is shown in Fig. 14. Previously it has sound speed perturbations near the sound channel axis. been shown [20] that ray intensities in a very different In addition, we have seen some evidence that an addi- chaotic system also fit a lognormal PDF; in that system tionalfactormaybeimportantintheAETenvironment. single scale isotropic fluctuations are superimposed on Namely,we haveobservedapositive correlationbetween a homogeneous background. (Note that all powers of a traveltime spreads and stability exponents; stability ex- lognormallydistributedvariablearealsolognormallydis- ponents in the AET environment are shown in Fig. 5. tributed,sorayamplitudeshavethispropertyifandonly We have chosen not to dwell on flat ray time spreads in if ray intensities have this property.) The apparent gen- this paper because there are no AET measurements of erality of the near-lognormalray intensity PDF suggests these spreads to which simulations can be compared. It that it applies generallyto raysystems that are far from is clear, however, that the issues just raised need to be integrable; the arguments presented in Ref. [20] suggest better understood. that this should be the case. The intensity distributions oftheearlyandlateAETarrivals,correspondingtosteep and flat rays, respectively, will be discussed separately. IV. WAVEFIELD INTENSITY STATISTICS Weconsiderfirsttheearlyarrivals. Wesawinthepre- vioussectionthatthemicromultipathsthatmakeupone In this section we consider the statistical distribution of these arrivals have the same ray identifier and have oftheintensitiesofboththeearlyandlateAETarrivals. a very small spread in travel time. The dominant mi- Recall that experimentally the early arrival intensities cromultipaths were seen (see Fig. 11) to have a time have been shown [24, 25] to approximately fit a lognor- spread of approximately 1 ms which is a small fraction mal probability density function (PDF) and the late ar- of the approximately 13 ms period of the 75 Hz carrier rival intensities have been shown [25] to fit an exponen- wave. Also, our simulations show that the Maslov in- tial distribution. The late arrival exponential distribu- dices of the dominantmicromultipaths differ by no more tionis notsurprisingasthis distributionischaracteristic thanoneunit. Theseconditionsdictatethatinterference of saturated statistics. The early arrival near-lognormal among the dominant micromultipaths is predominantly distribution is surprising, however, inasmuch as theory constructive. Becausetraveltimedifferencesaresosmall, [24, 31] predicts saturated statistics, i.e., an exponential the pulse shape should have negligible influence on the intensity PDF. It has been argued [38] that the theory distribution of peak intensities. In a model of the early can be modified in such a way as to move the early ar- AET arrivals consisting of a superposition of interfer- rival predictionfrom saturatedto unsaturated statistics. ingmicromultipaths,the micromultipathpropertiesthat 9 play a critical role in controlling peak wavefield intensi- value ofq , andν isthe true Lyapunovexponent. The 21 L ties are thus: 1) their amplitudes have a near-lognormal new(unbiasedinthesensedescribedabove)PDF,ρ′,for distribution; 2)the dominantmicromultipathshaveneg- uniformsamplinginδz isrelatedtothe previousoneby r ligible travel time differences; and 3) the dominant mi- x cromultipathshaveMaslovindicesthatdifferbynomore ρ′ (x) = ρ (x) |q21| x |q21| than one unit. It should be noted that very different h i behavior would have been observed if: the source band- 1 1 width were significantly more narrow as this would have = caused micromultipaths with different ray identifiers to s2πr(ν¯−νL)x · interferewithoneanother;orthesourcecenterfrequency (ln(x) ν¯r)2 weresignificantlyhigherasphasedifferences betweenin- exp − − , x 0 (26) " 2r(ν¯ νL) # ≥ terfering micromultipaths would then have been signifi- − cant. wherethefactorxaccountsfortheextracountingweight An additional subtlety must be introduced now: the of q , and x just preservesnormalizationandis calcu- PDF of the intensities of the constituent micromulti- 21 h i lated using Eq. (25). This calculation shows that the paths that make up a single arrival is not identical to unbiased micromultipath ray intensity PDF also has a the PDF described in Ref. [20] and shown in the upper lognormal distribution; the only change relative to the and middle panels of Fig. 14. The latter PDF describes biased PDF is an increase in the mean from ν r to ν¯r. the distribution of intensities of randomly (with uni- L Because lognormality is maintained, this correction rep- form probability) selected rays leaving the source within resents only a trivial change to the problem. some small angular band and whose range is fixed. This Approximate lognormality of the constrained (eigen- PDF is biased in the sense that it overcounts the micro- ray) PDF of ray intensity is shown in the lower panel multipaths with large intensities and undercounts those of Fig. 14. This PDF was constructed using the same with small intensities. Unbiased micromultipath inten- eigenraysthatwereusedtoproduce Figs. 8 and11. The sity PDFs can be constructed from the biased PDFs correspondingunconstrainedrayintensityPDFisshown shown. To doso,considera manifold(a smoothcurvein in the middle panel of Fig. 14. (Two constraints – fixed phase space corresponding to a fan of initial rays)which receiver depth and fixed ray identifier – are built into begins near (z ,p ) and arrives in the neighborhood of 0 0 the lower panel PDF. A very similar constrained PDF (z ,p ) at range r. Sampling in fixed steps of the differ- r r results if only the receiver depth constraint is applied, ential δp (as was done to produce the upper and mid- 0 provided ray launch angles are limited to the ‘steep’ ray dle panels of Fig. 14; this is equivalent to uniform ran- band used to construct the middle panel.) It should be dom sampling) leads to a highly nonuniform density of noted that the constrained(eigenray)PDF shown in the points on the final manifold since (z ,p + δp ) prop- 0 0 0 lower panel of Fig. 14 was constructed from numerically agates to (z + q δp ,p + q δp ). The greater q , r 21 0 r 11 0 21 found eigenrays; because weak eigenrays are difficult to the lower the density of points locally at final range. find numerically they are undercounted and the PDF is A uniform sampling in final position is acheived instead biased. In a practical sense this bias is of little conse- by considering the initial conditions that would lead to quence because the weak eigenrays that are difficult to (z +δz ,p +q /q δz ). Itsdensityofpointsontheini- r r r 11 21 r find contribute negligibly to the wavefield. tial manifold can be deduced from its initial condition, We return now to the problem of simulating the early (z ,p + δz /q ). It is necessary to sample q times 0 0 r 21 21 AET arrivals. It is tempting to think that because the more densely on the initial manifold in order to acheive constituent micromultipaths have a near-lognormal dis- uniform sampling in δz at range r. To account for this r tribution, the sum of their contributions should also be effect, we need to know the PDF for q with uniform 21 near-lognormallydistributed. Unfortunately, this is gen- initial sampling. Roughly speaking, the PDF of the ab- erally not the case. Consider, for example, the special solute values of the individual matrix elements of q have case in which phase, including Maslov index, differences the same form as for Tr(Q), apart from a shift of the | | are negligible. Then all micromultipaths interfere con- centroid that is lower order in range than the leading structively and peak wavefield amplitudes can be mod- term. From the results of Ref. [20], the (biased in the elled as the sum of many lognormally distributed vari- sense described above) probability that q falls in the 21 ables. Because all moments of the lognormal distribu- interval between x and x+dx is tion are finite, the central limit theorem applies. Under 1 1 these conditions, if sufficiently many contributions are ρ|q21|(x)=s2πr(ν¯ νL)x · summed, the distribution of the sums – the wavefield − amplitudes – would be a gaussian. (ln(x) ν r)2 TosimulatethestatisticsoftheearlyAETarrivals(re- L exp − − , x 0. (25) call Fig. 14 and the accompanying discussion) we have " 2r(ν¯ νL) # ≥ − used several variations of a simple model. An arrival Here ν¯ is a finite-range estimate of the Lyapunov ex- wasmodelledasasumofn interferingmicromultipaths m ponent based on an average (over an ensemble of rays) whose: 1) amplitudes are lognormally distributed; and 10 2) phases, ωT µ π/2mod2π, have a clearly identifi- arrivals,correspondingtothenear-axialrays. Here,time i i − ablepeak. Micromultipathcontributionswerecoherently spreads are sufficiently large that micromultipaths with added. Thepeakintensityofthesum–whosetraveltime different ray identifiers are not temporally resolved, i.e., is not known a priori – was then recorded. Using an en- arenotseparatedintime bymorethan(∆f)−1. Ateach semble of104 peak intensities,a peakintensity PDFwas (z,T) in the tail of the arrival pattern the wavefield can then constructed. Peak intensity PDFs constructed in be modelled as a superposition of micromultipath con- this fashion were found to be very close to lognormal; a tributions with random phases. The quadrature compo- typical example is shown in Fig. 15. In this example nents of the wavefield have the form of sums of terms of f = ω/2π = 75 Hz, the T ’s were identical, µ j,j+1 the form i i ∈ with equal probability (note that choice of the integer j x = a cos(φ ), and is unimportant) andn =5. Other combinationsof dis- i i i m tributions for Ti (either a Gaussian or the limiting case yi = aisin(φi) (27) of a delta distribution), µ (taken either from j,j+1 or i where φ is a random variable uniformly distributed on j 1,j,j+1withequalprobability),andthe parameter i − [0,2π). Thedistributionofa isclosetolognormal,buta n (between 2 and 100) were tested. These simulations i m correctionmust be applied to accountfor pulse shape as showed that provided the phase constraint noted above many of the interfering micromultipaths partially over- was satisfied, a near log-normal peak intensity PDF re- lap in time. This correction is unimportant inasmuch as sulted. the central limit theorem guarantees that, provided the Two points regarding this simple model are notewor- distributions of x and y have finite moments, the dis- i i thy. First, this model does not constitute a theory of tributions of sums of x and y converge to zero mean i i wavefield peak intensity statistics, but it does serve to gaussians. Thus, wavefield intensity is expected to have demonstrate that our ray-based simulations of the early anexponential distribution, consistentwith the observa- AET arrivals are consistent with the measured distributions. The comments made earlier about caustics apply tionofpeakintensities. Second,simulations(notshown), here as well. performed with E =1.0E yield sets of dominant mi- GM The question of what causes the transition from the cromultipathsthatviolateassumption2);phasesareuni- structured early portion of the AET arrival pattern to formly distributed and summing micromultipath contri- theunstructuredfinaleregiondeservesfurtherdiscussion. butions yields a distribution of peak intensities that is Recall that the early resolved arrivals have small time not close to lognormal. Thus, our simulations suggest spreads and peak intensities that are near-lognormally that a near-lognormal distribution of early arrival peak distributed, while the finale region is characterized by intensity requires a relatively weak internal wave field. unresolvedarrivalsandnear-exponentiallydistributedin- Acomplicationnotaccountedforinthepreceedingdis- tensities. In both regions the time spreads and intensity cussionis the presenceofcaustics. At causticsgeometric statistics are consistent with each other inasmuch as in amplitudes (8) diverge and diffractive corrections must our simulations a near-lognormal intensity distribution beapplied. Atshortrange(ontheorderofthefirstfocal is obtained only when there is a preferred phase, while distance – a few tens of kmin deep oceanenvironments) theexponentialdistributionisgeneratedwhenphasesare we expect that intensity fluctuations will be dominated random, i.e., when phases are uniformly distributed on by diffractive effects. The entire wavefield should be or- [0,2π). ganized by certain high-order caustics which leads to a With these comments in mind, it is evident that the PDFofwavefieldintensitywithlongtails[39,40,41]. In most important factor in causing the transition to the spiteoftheimportanceofdiffractiveeffectsatshortrange finale region is the increase in internal-wave-scattering- –andprobablyalsoatverylongrange–webelievethat, induced time spreads as rays become less steep; as inthetransitionalregimedescribedabove,intensityfluc- time spreads increase, neighboring timefront branches tuations are not dominated by diffractive effects. This blend together and the phases of interfering micromulti- somewhat counterintuitive behavior can be understood pathsgetrandomized. Thetrendtowardincreasingtime by noting that in the vicinity of caustics the importance spreads as rays become less steep is evident in Fig. 9. of diffractive corrections to (8) decreases as the curva- The surprising result is that the scattering-induced time ture of the caustic increases. Under chaotic conditions spreads of the steep arrivals are so small; we have seen the curvature of caustics increases, on average, with in- thatthis canbe explainedbymakinguseoftheapexap- creasing range, so the fraction of the total number of proximation. Asnotedattheendoftheprevioussection, multipaths that require caustic corrections decreases, on wedonotfullyunderstandthecauseofthetrendtoward average, with increasing range. This is true even as the larger time spreads as rays flatten. In the finale region number of caustics grows exponentially, on average, in internal-wave-scattering-inducedtimespreadsexceedthe range. This argument leads to the somewhat paradoxi- time difference between neighboring timefront branches calconclusionthat,priortosaturationatleast,weexpect thatwouldhavebeenobservedinthe absenceofinternal thattheimportanceofcausticcorrectionsdecreaseswith waves. Fig. 1 shows that these time gaps decrease as increasing range. rays become increasingly flat. Indeed, this figure shows WeturnourattentionnowtothelateunresolvedAET that even in the absence of internal waves there would