## Abstract

Simulations of laser-beam propagation through atmospheric turbulence and acoustic pulse propagation though ocean internal waves in the presence of a transverse waveguide are described. Both these problems are amenable to the parabolic wave equation for propagation in the forward direction. In the optical case, the question treated is the irradiance variance and spatial spectrum. In the ocean case, pulse propagation over long ranges is investigated. Determining the travel time of a pulse requires expanding the numerical simulation in a broadband, multifrequency calculation that takes even more time. Much effort has been expended in approximating the propagation by rays, so that the trajectory of the energy propagating from a given source to a given receiver can be tracked, and coherence calculations can be ray-based, using semi-classical formulas. This paper reviews the comparisons of several analytic ray-based approximations with numerical parabolic-equation simulations to determine their accuracies.

© 2002 Optical Society of America

## 1 Introduction

During the past ten years, significant progress in the field of wave propagation through random media has been obtained through the use of numerical simulation. The purpose of this article is to bring together for the first time these results from two quite different fields (atmospheric optics and ocean acoustics) to provide an opportunity for the optics community to have a review of their similarities and differences in this special issue of Optics Express.

#### Atmospheric Optics

In the case of optical propagation through atmospheric turbulence, it has been known for thirty years that the variance of irradiance as a function of turbulence strength has three regions: (1) weak fluctuations in which the variance rises linearly with turbulence strength; (2) strong focussing, in which the variance peaks at values well above unity; and (3) the asymptotic region, or strong fluctuations, in which the log of (variance – 1.0) decreases linearly as a function of the log of turbulence strength. In other words, the variance logarithmically approaches a limit of unity. (The explicit form of the approach to unity has only been found recently.)

In 1993, Consortini et al[1] measured irradiance variances for a spherical wave propagating over a Colorado mesa for a distance of about one kilometer, with the important addition that simultaneous measurements of both the strength and the inner scale of the turbulence in the propagation path were provided. The characteristics of the propagation put the propagation in the “strong-focussing” regime, near the peak of the variance curve as a function of turbulence strength.

A following paper[2] gave results for numerical calculations that simulated the geometry and wavelength of the Colorado experiment. The calculations provided irradiance variances under the assumption that the turbulence along the propagation path was homogeneous, isotropic, Kolmogorov turbulence, with the spectral cutoff at the Kolmogorov microscale provided by the fluid-dynamical form calculated by Hill and Clifford.[3] It was shown that the numerical simulation variances matched the experimental variances to an accuracy of about 5%. This success was in dramatic contrast to analytic theory for this region, that was shown to be inaccurate, with matches no better than 50–100%.

In 2000, simulations were carried further, to reach the asymptotic region in which the logarithmic behavior cited above is reached.[4] The slopes and intercepts of the behavior in the asymptotic region were evaluated, and compared with first-order asymptotic perturbation theories. The results were dramatically out of agreement; again the simulations had shown themselves to be dramatically superior to the analytic calculations.

The explanation for this disagreement is understood in terms of the failure of first-order asymptotic perturbation theory. Work was done to evaluate the second and third order contributions to a number of quantities, and the results were clear-cut; it is necessary to include third order terms in order to obtain the necessary accuracy. [5, 6, 7]

There is a question of the expected value of turbulence strength beyond which the analytic first-order theory should work. This is a matter of research.[8, 9]

Finally, the question of understanding these fluctuations, rather than just predicting them phenomenologically remains an issue.

#### Ocean Acoustics

In the case of ocean sound propagation, it has been realized for about 25 years that the medium fluctuations that cause most acoustic fluctuations are very different from turbulence; they are caused by internal gravity waves resulting from the continuous gradient in density with depth in the ocean.[10] Typical amplitudes of these waves at depths of several hundred meters are 5 to 10 meters. The density changes are primarily caused by temperature fluctuations, with a small effect from salinity. The internal wave field has a vertical correlation length of about one hundred meters, and a horizontal correlation length of several thousand meters; they are strongly anisotropic. They are also inhomogeneous, in that the variance of sound speed varies with depth, more or less exponentially with a scale length of 1000 meters.

In addition to the complications of the internal-wave field, the entire propagation occurs within a continuous-gradient waveguide. The deterministic sound speed in the ocean has a minimum at a depth of about 800 m, with higher values at the surface and bottom (whose depth is 4 to 5 km). The sound-speed values at top and bottom are about 3% greater than the minimum. Thus propagation of rays at angles within about 10 degrees from horizontal are trapped within the waveguide. In this way, energy can propagate thousands of kilometers and still have enough energy flux to be detected easily with present technology. In passing through these large distances, the internal waves induce the fluctuations that are observed in the receivers.[10, 11]

In the early 1990’s, the results and analysis of a detailed, 1000-km, 250-Hz experiment in the North Pacific were presented.[12, 13] Pulses were propagated every ten minutes over a 21 hour period, allowing the observation of the entire internal-wave spectrum. The wavefront resulting from each pulse was observed by fifty hydrophones vertically spanning 3000 m; because of the waveguide, the original pulse is observed as many individual wavefront segments, connected in a broken line.[12] The transmitted pulse width restricted the measurement accuracy for travel time to about 2 ms. The fluctuations due to internal waves were 10 ms. The asymmetric character of the waveguide resulted in a concentration of energy in the last second of the arrival pulse; a region in which separation into separate wavefronts was not possible. The use of parabolic-equation numerical simulation was quite successful for this experiment. The use of ray simulation was successful even in the last second of the arrived pulse.

To scale the 1000-km experiment relative to a similar optical experiment, note that the wavelength of a 250-Hz sound wave is 6 m. Therefore, the waveguide has a gradient scale length of 150λ, and the entire waveguide width is about 800λ. The propagation range is 150,000λ. For 0.5*μ* light, this would correspond to a fiber of diameter 400*μ* and length 70 mm.

An experiment from the year 2000 was similar except that the range was 3250 km; the observed arrival behavior was similar to the 1000-km experiment in its agreement with simulation.

Beginning in the mid-seventies, an attempt was made to analytically calculate the variance of travel time and the coherence functions as functions of vertical and horizontal separations, and of time and propagation frequency.[10, 11] These calculations were possible because in the strong fluctuation regime it is possible to average over the medium fluctuations first, before doing the path integral. The resulting formulas are in the form of integrals along a geometric-optics ray, with the Fresnel length providing a scale of transverse size beyond which any diffraction effects are ignored.

The path-integral formulas were quite successful in explaining detailed structure in experiments whose range was less than 100 km.[14, 15, 16, 17] For many years experiments with ranges in the thousands of kilometers were explained in terms of extensions of these formulas, for want of a better way. However, very recently it has been shown that while predictions of travel-time variances continue to be accurate at least to 1000 km, coherence functions and pulse spreads fail very strongly beyond a few hundred kilometers.[18]

#### Rays through random media

Another approach has been followed by a number of researchers that is different from trying to add diffraction effects onto deterministic rays. Instead, they follow rays through the random medium, allowing it to experience a random walk away from their deterministic trajectories, ignoring diffraction effects, and then use semi-classical interference to generate coherent effects at the receiver. This very different technique can be compared again with the (assumed accurate) parabolic equation simulation. The results of this work are very encouraging, as they not only can predict travel-time variance and bias, but also can predict the behavior in the last second of the pulse.[19] It has not been determined whether or not this technique can be used to predict coherence functions.

## 2 The Parabolic Equation

The standard parabolic equation is used:

where *ψ* is the wave function,
${q}_{0}=\frac{2\pi}{\lambda}$ and λ is the wavelength of the propagating wave. The quantity *U*
_{0} is the deterministic (average) deviation of the sound speed from a constant. (For optical propagation through turbulence *U*
_{0} = 0, while for acoustic propagation through the ocean, *U*
_{0} represents the waveguide.) The quantity *μ*(*x*,*y*,*z*) is the fractional sound-speed fluctuation from zero due to the statistically fluctuating medium. Any PE equation has the critical advantage that it can be advanced through a series of range steps (in *x*) without recourse to costly iterative or relaxation methods.[10, 20, 21, 11]

Note that in the optical propagation case, both transverse dimensions *z* and *y* are used. In the acoustic propagation in the ocean case, only the vertical coordinate *z* is used, because internal waves are so anisotropic. However, the ocean case is broadband, requiring of the order of 1024 frequencies, so the required computer times are similar.

## 3 Optical Propagation through Atmospheric Turbulence

Define ${\beta}_{0}^{2}$ as the weak-fluctuation irradiance variance. In the case of a plane-wave initial condition:

where ${C}_{n}^{2}$ is a measure of turbulence strength, *k* is the wavenumber of the propagating wave, and *L* is the range from the source to the final screen. In the case of a spherical wave:

Consider a simulation characterized by ${\beta}_{0}^{2}$, and imagine a number of screens *n*_{s}
of size *N*_{s}
. Also, the Fresnel length *R*_{f}
will be necessary. Atmospheric turbulence has two parameters: its strength ${C}_{n}^{2}$ and its inner scale *l*
_{0}. If the fluctuations are weak, the inner scale is negligible, so that it does not appear in the previous equations, but in the strong-fluctuation regime it is crucial. How to do an accurate simulation in the face of all these parameters is explained (with references) in Flatté and Gerber.[4]

#### 3.1 Spatial Spectra

Figure 1 shows variance preserving[22] spectra for both plane and spherical-wave initial conditions. This figure shows that the spectra have been well sampled so that the total variance (which is the integral under each curve) is well estimated.

Figures 2 and 3 show the same spectra with a logarithmic ordinate (and without the extra factor of *κ*). Each spectrum has been divided by its respective ${\beta}_{0}^{2}$ so that the asymptotic behavior at large *κ* for the cases with zero inner scale all fall on the same line.

#### 3.2 Irradiance Variances

Figures 4 and 5 show the variances minus unity for both plane and spherical wave cases.[4]

The quantity ${\sigma}_{I}^{2}$ – 1 is plotted on a log-log plot so that the data will follow a straight line if the expected asymptotic behavior is followed. Each of the data sets has been fit to functions of the following form:

where ${\beta}_{0}^{2}$ = 60.6 and the parameters *a* and *b* are different for each initial condition and inner scale. The results of these fits are shown in Tables 1 and 2. For finite inner scale, all the slopes for a given initial condition are close to being the same. Therefore lines with the average slope have been plotted over each data set corresponding to a given inner scale.

Figure 6 shows the value of ${\sigma}_{I}^{2}$ – 1 for each fitted line at a value of ${\beta}_{0}^{2}$ equal to ${\beta}_{c}^{2}$ as a function of inner scale.

From the figure it can be seen that simulation points do not lie on the straight line fits within uncertainties. However, these lines do represent the first-order behavior with inner scale reasonably well.

The predictions of analytic theory were first done for zero inner scale by Gochelashvili and Shishov,[23] and then for finite inner scale by Fante[24] and Frehlich.[25] Updates to some of the values of *a* have been made by Andrews et al.[8, 9] These predictions are shown as dashed lines on the figures. The simulations are dramatically different from the analytic theory.

Figures 2 and 3 show that the simulations are sampling large scales sufficiently. However, the very fact that finite-size screens are being used introduces a reduction in the spectrum at large scales. In the past, it was considered adequate to represent the Kolmogorov spectrum accurately down to wavenumbers of 0.01 the size of the Fresnel length. In the work of Flatté and Gerber, this was done.

In a recent publication, Andrews et al[26] have shown that the comparison of analytic theory with simulations in Figures 4 nd 5 needs modification when an outer-scale is considered. Their results were that near the peak of the scintillation index, the effect is relatively small, but at values near ${\beta}_{0}^{2}$ = 100 the effects of outer scales of about 1 m can be of order unity in the value of the scintillation index. Since the statistical errors in the simulations are of order 0.05 or less, this is crucial. The Fresnel zone for optical propagation over 1 km is 1 cm, so that the outer scale of importance is just at the ratio of 100. This implies that the simulations need to be compared not with the analytic theories assuming infinite outer scale, but rather to the analytic expressions of Andrews et al. This has not been done as yet. Alternatively, the screen sizes need to be increased by a factor of ten. This is not feasible in the near future.

Another related problem is whether the first-order analytic theory will be adequate at ${\beta}_{0}^{2}$ = 100 in any case. Previous work in this area by Dashen et al[27] has shown clearly that at ${\beta}_{0}^{2}$ values comparable to 100, the second and even third order terms are significant at this level. Thus we are still far from reaching the asymptotic limit in which analytic theory can be used to predict results to the accuracy available in the simulations.

#### 3.3 Empirical Power-Law Formulae

In order to bring the inner scale into the equation for variance, Flatté and Gerber generalized Equation 4 to:

Again the simulations dramatically differ from the analytic predictions. Table 3 lists the values of the variances at ${\beta}_{c}^{2}$ from the straight-line fits, as well as the parameter values for the straight-line fit to the behavior as a function of inner scale.

## 4 Acoustical Propagation through Ocean Internal Waves

#### 4.1 Introduction

Analytic, line-integral approximations derived from the acoustic path integral have been used to estimate the magnitude of the fluctuations in a signal travelling through the ocean.[28] These approximations for the bias and variance of travel time, the length scale of acoustic coherence in depth, and the spreading of the acoustic pulse have been compared to values from simulations that use the standard parabolic equation (PE).[18] Acoustic propagations with a frequency of 250 Hz were examined with a maximum range of 1000 km, and using two different temperate-latitude sound-speed profiles. The sound speed was perturbed by internal waves conforming to the Garrett-Munk (GM) spectral model with strengths of 0.5, 1 and 2 times the GM reference energy level. Estimates of the rms travel-time variability were on average within statistical uncertainty at 1000 km. Bias approximations were accurate to this level for the first few hundred kilometers of propagation, but deviated quickly beyond. Predicted depth-coherence lengths and pulse spreads at 1000 km differed significantly from simulation results beyond 100 km. Results imply that long-range acoustic tomography using identified rays is restricted to the use of travel-time variance.

Hardin and Tappert first applied the parabolic equation to the problems of ocean acoustics and presented a practical computer algorithm that can solve it to any desired accuracy.[29, 21] The use of the parabolic equation in numerical simulations of acoustic propagation through ocean internal waves was first described by Flatté and Tappert.[30] They investigated the fluctuations in a 100 Hz signal propagating to a range of 100 km.

At about the same time, it was shown that the path integral for propagation through a random media that is characterized by inhomogeneity, anisotropy, and a deterministic waveguide – as in the case of ocean acoustics – could be approximately solved in terms of line integrals along deterministic rays.[10, 15] The accuracy of the approximate integrals have been compared with computer-intensive parabolic-equation (PE) results.[18] The quantities compared include: the root-mean-square variation in travel time (*τ*), the internal-wave-induced travel-time bias (*τ*
_{1}), the induced spread of the acoustic pulse (*τ*
_{0}), and the length scale characterizing the coherence of the acoustic signal at different depths (*z*
_{0}). An important motivation for an investigation into the accuracy of these approximations is their computational efficiency and the relative simplicity of subsequent analysis. Calculations using both techniques were performed with two different sound-speed profiles – a Munk Canonical profile and a profile from the Slice89 experiment[12] -for two acoustic frequencies – 250 Hz and 450 Hz. Three different internal-wave strengths – 0.5, 1 and 2 times the reference Garrett-Munk (GM) level – were employed in both the integral approximations and parabolic-equation simulations for each situation. For the parabolic simulations, fifty different realizations of the stochastic internal-wave field were used for cases with 2 GM internal waves and a 250 Hz acoustic frequency; ten realizations were used in all other circumstances. Information was recorded from all calculations at every 100-km interval to a maximum propagation range of 1000 km.

#### 4.2 Ocean Model

The simulations described in the following treat the ocean environment as an average sound-speed profile in depth perturbed by internal gravity waves that conform to the Garrett-Munk (GM) spectral model. Two different profiles were used, each possessing features appropriate to the temperate latitudes. The Munk Canonical profile[31] is an analytic model, while the other sound-speed profile was measured during the Slice89 experiment that took place in the North Pacific.[12] The two sound-speed profiles are displayed in Figure 7.

Both aspects of the simulated ocean – the deterministic profile and the stochastic perturbations – are expressed in the sound-speed equation,[10]

where *U*
_{0}(*z*) represents the dependence of the background profile, *c*(*z*) = **c**
_{0}[1+*U*
_{0}(*z*)], on depth and *μ*(*x*, *t*) contains the effect of internal waves on the speed. A useful definition of *c*
_{0} is the average of *c*(*z*) over depth. Ocean internal waves are not affected by pressure. For adiabatic vertical displacements, the value of *μ* can be expressed,[10]

where *ζ*(*x*, *t*) is the magnitude of the displacement and [*∂*_{z}*U*
_{0}(*z*)]_{p} is the fractional potential gradient of the sound-speed profile. (“potential” means with the pressure effect removed)

In order to define the internal-wave spectrum of *ζ* several parameters are needed, as well as the profile of buoyancy frequency *N*(*z*). (See Figure 8.) (The buoyance frequency is a measure of the density gradient that allows internal waves.)

The Garrett-Munk (GM) spectrum of internal wave displacement is:

where *j* is the vertical mode number, *k* is the magnitude of the horizontal wave number (the spectrum is horizontally isotropic), *E* = 6.3 ∙ 10^{-5} (for the reference internal-wave energy referred to as the Garrett-Munk energy level), *N*
_{0} = 3 cycles per hour and *j*
_{*} = 3 are empirical constants,
${k}_{j}\equiv \frac{\pi {\omega}_{i}j}{{N}_{0}B}$ and M is the normalization constant for the sum over mode number, *M* = ${\sum}_{j=1}^{\infty}$(*j*
^{2} + ${j}_{*}^{2}$)^{-1}.[32, 33, 34]

The numerical internal-wave field for the parabolic simulations was constructed using a technique due to Colosi and Brown.[35]

For the analytic integrals, three derived moments of the spectrum are needed: 〈*μ*
^{2}〉; the correlation length parallel to a given ray direction, called *L*_{p}
; and a measure of the transverse (which means vertical) inverse correlation length squared, designated by {${k}_{\upsilon}^{2}$}.

An empirical expression for *L*_{p}
has been shown to be sufficiently accurate:[28, 36] Typical internal-wave correlation lengths are a few kilometers in the horizontal direction and *O*(100) m in the vertical.

The quantity {${k}_{\upsilon}^{2}$} is given by ${l}_{\upsilon}^{-2}$} where *l*_{υ}
depends on experimental constants like *j*
_{*} and ∫ *N*(*z*)*dz*.[34] Details are given in the reference; the result is:

and *l*_{o}
≈ 1000 m.

#### 4.3 Parabolic-Equation Simulations

Simulations with the Canonical profile had a source depth of -1000 m; the source in the Slice89 simulations was located at a depth of -800 m. A total of 1024 frequencies spanned the 100 Hz bandwidth. The step size in range was 25 m and the number of vertical grid points spanning the water column from surface to the bottom, *z*_{b}
= 5118.75 m, was 2048. In order to implement the reflecting surface boundary condition, an image ocean was constructed above the surface for a total of 4096 vertical points.

The simulated acoustic timefronts display the expected effects of internal-wave perturbations, including both small-scale acoustic fluctuations (Figures 9 and 10)[12] and the reshaping of the transmission finalé (Figures 11 and 12).[13]

#### 4.4 Integral-Approximation Techniques

The integral approximations presented here begin with the calculation of ray trajectories, subject only to the influence of the depth-dependent sound-speed profile. The results of the analytic approximations take the form of weighted integrals over the deterministic ray from the source to the receiver. That ray path is periodic in range, because the deterministic sound-speed profile is range independent. (Each ray path has an identifier, labelled *ID*, that counts the number of turning points from source to receiver.) The weighting functions that enter into the integrals are not range independent, nor even periodic in range. Nevertheless, the nature of the weighting functions is such that great gains in computational speed are achievable compared with straightforward integration.[28]

These integral expressions provide estimates of parameters in second moments of the acoustic field function. The second moment for changes in depth can be expressed as,

where *z*
_{0} is the characteristic depth scale that will be calculated using an integral along the ray trajectory. Similarly, for changes in frequency,

where *σ* is the acoustic frequency, *τ*
^{2} gives the variance of travel time, *τ*
_{1} is the internal-wave bias (the average difference in travel time caused by the presence of the speed disturbances), and *τ*
_{0} describes the spreading of an intensity peak in time due to internal waves.[37, 15]

The variance of travel time, *τ*
^{2}, is calculated from the properties of internal waves in the following way,[10]

As previously stated, the integral is performed along the ray trajectory. The sound speed variance, 〈*μ*
^{2}(*z*)〉, is a profile in depth. The quantity *L*_{p}
(*θ*, *z*) expresses the correlation length of the internal waves along the ray direction. *τ*
^{2} is the simplest of the integral approximations.

The integrals in the following expressions involve the depth separation of nearby unperturbed rays which differ by defined amounts at the source or receiver. For the ray, *z*_{ray}
(*x*), and a nearby ray, *z*(*x*), this separation is,

The behavior of the vertical separation in the sound channel, when assumed small, is governed by,[10, 28]

Of particular interest is the solution, *ξ*
_{1}, that yields the vertical separation as a function of range for two rays which begin at the same source location and arrive at the receiver with unit separation. This function will appear in the depth-coherence integral. The Green’s function generated by the related equation,* ∂*_{x′x′}*g*(*x*, *x*′)+*U*
_{0}″*g*(*x*, *x*′) = *δ*(*x*′-*x*), is another necessary component of the acoustic integrals to be performed.

The acoustic vertical coherence length can be expressed as an integral along the ray trajectory in the following way,

The internal-wave bias, *τ*
_{1}, is the average change in travel time caused by the presence of sound-speed fluctuations. It is calculated according to,

Analytic integral techniques also offer an estimate of the pulse spreading due to internal-wave fluctuations, *τ*
_{0}. This quantity characterizes the additional width in time of the intensity peak for a timefront ID segment, at the ray’s arrival depth, caused by the presence of internal-wave fluctuations. The pulse spreading is given by,

The inverse of *τ*
_{0} is also one contribution to the estimate of the coherent bandwidth.

These integral approximations were implemented using a computer code package referred to as “CAFI” (Computation of Acoustic Fluctuations from Internal waves) prepared by Flatté and Rovner.[28] Approximately 6500 ray paths were calculated in each case at launch-angle intervals of 0.04 degrees.

#### 4.5 Comparison of Results from Integral Expressions and the Parabolic Equation

The four quantities above (*τ*, *τ*
_{1}, *τ*
_{0}, and *z*
_{0}) that have been calculated from integral approximations have been compared to corresponding results from simulations using the parabolic equation: The results obtained at 450 Hz did not differ qualitatively from those obtained at 250 Hz.

### 4.5.1 Root-Mean-Square Travel-Time Variability: τ

Results for *τ* as a function of range, averaged over ID, are displayed in Figure 13. The values for individual timefront ID segments at a range of 1000 km are shown in Figure 14.

This quantity is among the most frequently used in the analysis of experimental data.[12] As can be seen in Equation 12, the integral approximation is not a function of the acoustic frequency. The predicted magnitude of *τ* increases with internal-wave strength and also varies somewhat with the choice of sound-speed profile and timefront ID number. The statistical uncertainty in the PE values was approximately 10% for the cases employing fifty realizations.

Slice89 *τ* analytic values were within statistical uncertainty of the 1000-km simulation results. Canonical *τ* analytic values were within the statistical error of the simulation results in a few cases, at low timefront ID numbers and internal-wave strengths. In other cases they overestimated by between 20% and 60%.

### 4.5.2 Internal-Wave-Induced Travel-Time Bias: τ_{1}

The calculated travel-time bias due to internal waves, *τ*
_{1}, is displayed in Figure 15.

The integral approximation is a function of acoustic frequency in this case, though only as an overall factor that is independent of range. The choice of sound-speed profile significantly affects the result, particularly for long ranges and strong internal waves. The approximations and the PE results disagree significantly by the end of the 1000-km propagation range. They are, however, in reasonable agreement for the first few hundred kilometers (approximately 400 – 500 km) in range. At 500 km, the Slice89 analytic approximations are within the statistical uncertainty of the simulations. Except for the case with 2 GM internal waves, the Canonical approximation is within the simulation error at 400 km. Beyond this range, the discrepancies grow quickly.

### 4.5.3 Depth-Coherence Length: z_{0}

Note that *z*
_{0}, given in Equation 15, describes the scale of a gaussian function.

Coherence-function results from both the simulations and integral techniques are displayed in Figure 16 for an acoustic frequency of 250 Hz and a range of 1000 km.

The end result uses Equation 10 with *z*
_{0} determined from Equation 15 averaged over ID. The parabolic-equation results are shown averaged over all ID segments in each case. The coherence functions predicted by integral methods decrease over substantially shorter length scales than the parabolic-simulation results in all cases

Analysis of the Slice89 experiment suggested an internal-wave strength of 0.5 GM.[38, 13] For those conditions, integral approximations of the coherence length are approximately 90 m and parabolic results are roughly 400 m. Neither is completely outside the experimental range, but the two methods differ strongly.

### 4.5.4 Internal-Wave-Induced Pulse Spreading: τ_{0}

The investigation into the accuracy of the integral approximation for the pulse spreading caused by the presence of internal waves, *τ*
_{0}, is based on a comparison of the mean pulse shapes, 〈*I*(*t* – *T*_{i,ID}
)〉, from parabolic simulations with internal waves to those without internal waves. The construction of these mean pulses depends on the use of multifrequency techniques. Pulses are pictured in Figures 17 and 18 for each speed profile, an acoustic frequency of 250 Hz, and an internal-wave strength of 1 GM at a range of 1000 km. The pulse-spread approximations, calculated according to Equation 17 and averaged over each timefront ID, are included as numbers within each box.

Internal-wave-influenced mean pulses from the parabolic-equation simulations are often nearly indistinguishable from those calculated without internal-wave disturbances. Even when the pulses at 1000 km display the sort of spreading that could potentially be explained in terms of a Gaussian function with a width increase given by a single τ_{0} parameter, the predicted value greatly exceeds the spreading seen in the PE results, often by an order of magnitude. An example of this situation occurs in the Slice89 profile for ID-37.

Similar findings for the integral approximation of pulse spreading were reported for the 3250 km, 75 Hz Acoustic Engineering Test (AET), part of the Acoustic Thermometry of Ocean Climate (ATOC) project.[39, 40] It was found that the *τ*
_{0} approximation was too large by more than an order of magnitude. It was also reported that the effects of internal waves on the mean pulse shape are more accurately characterized as increases in the mean intensity at times away from the maximum, rather than as spreading of the central peaks themselves.

## 5 Geometrical-Optics Rays: Energy Diffusion and Semi-classical Interference

An analysis of the effects of ocean internal waves on long-range acoustic pulse propagation from the geometrical-optics point of view has been done.[19] The chaotic behavior of rays and the micro-folding of timefronts were investigated. It was found that the extent of the region of the timefront in which strongly chaotic rays appear, and the strength of the rays’ sensitivity to initial conditions, depend on the average (range-independent) sound-speed profile, on the range from the source to the receiver, and on the internal-wave spectral model, but not on the specific realization of the internal waves. Figure 19 shows the timefront calculated with geometrical optics rays for a 1000-km propagation, first without internal waves, and then for three realizations with internal waves.

The diffusion of energy in the last second of the pulse is evident. This is observed in the experiment and in PE simulations. Thus, for this particular experiment (SLICE89), it was concluded that the observed depth diffusion of energy in the late-arriving portion of the timefront is a result of refraction (of geometrical-optics rays), not diffraction.

Given that following rays gives accurate results for the arrival time and depth at which energy arrives, tomography can be analyzed in detail. After following 65,000 rays from the source, the much smaller number of rays that arrived within one meter of a specific receiver depth were selected. The different paths of these rays defines the resolution of tomography. Figure 20 shows that the bundle of rays that start at a source point and end at a receiver point 1000 km away have extents that are about 100 m vertically and 10 km horizontally.

It is also of interest to attempt the synthesis of the coherent field at the receiver by adding these rays, giving each a semi-classical amplitude and phase determined from its travel time and divergence of another ray with a small difference in launch angle. An example of this synthesis is shown in Figure 21, where a multitude of rays arrive near enough in arrival time that they interfere strongly. The figure shows that the arrival time is well approximated by the semi-classical sum, and the amplitude is within 3 db of the PE solution. Another example is shown in Figure 22, where two peaks are shown in both PE and ray-wave solutions. It was found that the amplitudes agreed very well at 1000 Hz, but differed by up to 5 db at 250 Hz.[19]

The chaotic behavior of rays in certain regions of launch angles were investigated in some detail in Simmen et al. That discussion has been avoided due to lack of space here.

## 6 Summary and Conclusion

Optical propagation through atmospheric turbulence has been simulated for turbulence strength varying from (0 < ${\beta}_{0}^{2}$ < 100) and inner scale in the region (0 mm < *l*
_{0} < 8 mm). In the strong-focussing regime (near the variance peak) the simulations agree very well with experiment. In the strong fluctuation regime (where logarithmic scaling is observed) the simulations show scaling, but the scaling slopes do not agree at all with the analytic predictions from first order; it is known that third order is needed.

In the strong-fluctuation, asymptotic regime, analytic theory predicts that the irradiance variance difference from unity behaves like power laws as a function of turbulence strength (represented by 30 < ${\beta}_{0}^{2}$ < 100) and inner scale. It is found that the irradiance variance indeed behaves qualitatively in the expected manner, but that the power-law indices are dramatically different from the analytic predictions. The analytic theories do have a prediction as to the values of ${\beta}_{0}^{2}$ and *l*
_{0} at which they should be valid, and our simulations are in the region of expected validity.[8, 9] It then becomes an important avenue of future research to understand the cause of the discrepancies between analytic theory and simulations.

The ocean-acoustic comparisons that have been made have serious implications for the usefulness of the analytic integral expressions in providing accurate estimates of some fluctuation quantities. Acoustic tomography depends, to a large extent, on the use of such quantities, so these results might complicate the process of acoustic tomography. The fundamental quantity *τ* is found to be successful, to which we attribute the correctness of the basic assumptions that go into that analytical integral expression. The quantities *τ*
_{1}, *τ*
_{0}, and *z*
_{0} from the analytic integral expressions do not correspond to the PE results over the entire propagation range. Speculations have been made as to the reasons for this failure.[18]

The analytic formulas for the other quantities (*τ*
_{1} and *τ*
_{0}) make the assumption that internal-wave effects at separation distances of order a Fresnel length are accurate. But a Fresnel length at these small frequencies is of order several hundred meters. It seems that the expressions for the medium fluctuations fail at those large separations.

It is also quite important to note that the analytic integral expressions have been derived under the assumption of a single acoustic frequency at the source, while the PE simulations are done with a bandwidth of 100 Hz. Thus single-frequency expressions are being compared with broadband numerical simulations. Perhaps the PE simulations are providing an average over frequency that smooths out the effect. The same kind of ideas appear to apply to *τ*
_{0}.

It is useful to characterize the conditions in which the analytical integral expressions are useful. Situations that involve a single frequency and short range (such that a single ray can be identified by beamforming) should allow the expressions to be used. The successful comparisons with real data from the Azores Fixed Acoustic Range (AFAR) experiment are examples.[10, 14] In situations that involve pulse propagation, or that involve ranges beyond *O*(100) km, the analytic integrals appear not to be useful, with the single exception of the travel-time variance, which does not depend either on the Fresnel length or the internal-wave-induced vertical curvature of sound speed.

From the point of view of acoustic tomography, the use of travel-time fluctuations for determining the distribution of internal-wave spectral strength is well borne out. These results offer additional support for the utility of *τ* as a fundamental parameter in tomographic measurements of internal-wave strength.

Experiments have shown that there is another useful measurement that does not depend on analytic integral estimates of the type considered here. The behavior of the finale of the pulse, within which one cannot identify timefronts, has yielded clear measurements of the internal-wave strength.[13] This measurement depends on theoretical calculations which can be either of two types: PE simulations or ray trajectory following through the random medium. Thus we have two useful experimental measurements of internal-wave strength at a range of 1000 km (and presumably farther, for the reasons discussed above), and no other useful measurements that might provide other views of the internal-wave spectrum beyond a few hundred kilometers. We have not, in this paper, investigated the possibility of using intensity within the pulse.[40]

## 7 Acknowledgments

This research was supported in part by the U.S. Office of Naval Research, Code 321OA. The help of many people over a period of 25 years is gratefully acknowledged. The deaths of Roger Dashen, Fred Tappert, and Fred Zachariasen are considerable losses to the furtherance of this work. The help of recent graduate students James Gerber and Michael Vera in creating this body of work is gratefully acknowledged.

## References and links

**1. **A. Consortini, F. Cochetti, J. H. Churnside, and R. J. Hill, “Inner-scale effect on intensity variance measured for weak to strong atmospheric scintillation,” J. Opt. Soc. Am. A **10**, 2354–2362 (1993). [CrossRef]

**2. **S. M. Flatté, G. Y. Wang, and J. Martin, “Irradiance variance of optical waves through atmospheric turbulence by numerical simulation and comparison with experiment,” J. Opt. Soc. Am. A **10**, 2363–2370 (1993). [CrossRef]

**3. **R. Hill and S. Clifford, “Modified spectrum of atmospheric temperature fluctuations and its applications to optical propagation,” J. Opt. Soc. Am. **68**, 892–899 (1978). [CrossRef]

**4. **S. Flatté and J. Gerber, “Irradiance variance behavior for plane- and spherical-wave optical propagation through strong turbulence,” J. Opt. Soc. Am. A **17**, 1092–1097 (2000). [CrossRef]

**5. **R. Dashen and G.-Y. Wang, “Intensity fluctuation for waves behind a phase screen: a new asymptotic scheme,” J. Opt. Soc. Am. **A10**, 1219–1225 (1993). [CrossRef]

**6. **R. Dashen, G.-Y. Wang, S. M. Flatté, and C. Bracher, “Moments of intensity and log intensity: new asymptotic results for waves in power-law media,” J. Opt. Soc. Am. **A10**, 1233–1242 (1993). [CrossRef]

**7. **R. Dashen, G. Wang, S. M. Flatté, and C. Bracher, “Moments of intensity and log-intensity: New asymptotic results for waves in power-law media,” J. Opt. Soc. Am. A **10**, 1233–1242 (1993). [CrossRef]

**8. **L. Andrews and R. Phillips, *Laser beam propagation through random media* (SPIE Publishing Services, Bellingham, WA, 1999).

**9. **L. Andrews, R. Phillips, C. Hopen, and M. Al-Habash, “Theory of optical scintillation,” J. Opt. Soc. Am. A **16**, 1417–1429 (1999). [CrossRef]

**10. **S. Flatté, R. Dashen, W. Munk, K. Watson, and F. Zachariasen, *Sound Transmission Through a Fluctuating Ocean* (A 300 page monograph published by the Cambridge University Press in their series on Mechanics and Applied Mathematics, 1979).

**11. **S. Flatté, “Wave propagation through random media: Contributions from ocean acoustics,” Proc. of the IEEE **71**, 1267–1294 (1983). [CrossRef]

**12. **T. Duda, S. M. Flatté, J. Colosi, B. Cornuelle, J. Hildebrand, W. Hodgkiss Jr., P. Worcester, B. Howe, J. Mercer, and R. Spindel, “Measured wavefront fluctuations in 1000-km pulse propagation in the Pacific Ocean,” J. Acoust. Soc. Am. **92**, 939–955 (1992). [CrossRef]

**13. **J. Colosi, S. M. Flatté, and C. Bracher, “Internal-wave effects on 1000-km oceanic acoustic pulse propagation: Simulation and comparison with experiment,” J. Acoust. Soc. Am. **96**, 452–468 (1994). [CrossRef]

**14. **S. Reynolds, S. Flatté, R. Dashen, B. Buehler, and P. Maciejewski, “AFAR measurements of acoustic mutual coherence functions of time and frequency,” J. Acoust. Soc. Am. **77**, 1723–31 (1985). [CrossRef]

**15. **R. Dashen, S. Flatté, and S. Reynolds, “Path-integral treatment of acoustic mutual coherence functions for rays in a sound channel,” J. Acoust. Soc. Am. **77**, 1716–22 (1985). [CrossRef]

**16. **S. Flatté, S. Reynolds, R. Dashen, B. Buehler, and P. Maciejewski, “AFAR measurements of intensity and intensity moments,” J. Acoust. Soc. Am. **82**, 973–980 (1987). [CrossRef]

**17. **S. Flatté, S. Reynolds, and R. Dashen, “Path-integral treatment of intensity behavior for rays in a sound channel,” J. Acoust. Soc. Am. **82**, 967–972 (1987). [CrossRef]

**18. **S. Flatté and M. Vera, “Comparison between ocean-acoustic fluctuations in parabolic-equation simulations and estimates from integral approximations,” J. Acoust. Soc. Am. in review (2002).

**19. **J. Simmen, S. M. Flatté, and G.-Y. Wang, “Wavefront folding, chaos, and diffraction for sound propagation through ocean internal waves,” J. Acoust. Soc. Am. **102**, 239–255 (1997). [CrossRef]

**20. **F. Jensen, W. Kuperman, M. Porter, and H. Schmidt, *Computational Ocean Acoustics* (American Institute of Physics Press, 1994).

**21. **F. Tappert, “The parabolic approximation method,” In *Wave Propagation and Underwater Acoustics*, J. Keller and J. Papadakis, eds., pp. 224–287 (Springer-Verlag, 1977). [CrossRef]

**22. **S. M. Flatté, C. Bracher, and G.-Y. Wang, “Probability-density functions of irradiance for waves in atmospheric turbulence calculated by numerical simulation,” JOSA **A11**, 2080–2092 (1994).

**23. **K. Gochelashvili and V. Shishov, “Saturated fluctuations in the laser radiation intensity in a turbulent medium,” Sov. Phys. JETP **39**, 605–609 (1974).

**24. **R. L. Fante, “Inner-scale size effect on the scintillations of light in the turbulent atmosphere,” J. Opt. Soc. Am. **73**, 277–281 (1983). [CrossRef]

**25. **R. Frehlich, “Intensity covariance of a point source in a random medium with a Kolmogorov spectrum and an inner scale of turbulence,” J. Opt. Soc. Am. **A4**, 360–366 (1987). [CrossRef]

**26. **L. Andrews, M. Al-Habash, C. Hopen, and R. Phillips, “Theory of optical scintillation: Gaussian-beam wave model,” Waves Random Media **11**, 271–291 (2001). [CrossRef]

**27. **R. Dashen, G.Y. Wang, Stanley M. Flatté, and C. Bracher, “Moments of intensity and log-intensity: new asymptotic results from waves in power-law random media,” J. Opt. Soc. Am. A **10**, 1233–1242 (1993) [CrossRef]

**28. **S. Flatté and G. Rovner, “Calculations of internal-wave-induced fluctuations in ocean-acoustic propagation,” J. Acoust. Soc. Am. **108**, 526–534 (2000). [CrossRef]

**29. **R. Hardin and F. Tappert, “Applications of the Split-Step Fourier Method to the Numerical Solution of Nonlinear and Variable Coefficient Wave Equations,” Society for Industrial and Applied Mathematics Review **15**, 423 (1973).

**30. **S. Flatté and F. Tappert, “Calculation of the effect of internal waves on oceanic sound transmission,” J. Acoust. Soc. Am. **58**, 1151–1159 (1975). [CrossRef]

**31. **W. H. Munk, “Sound channel in an exponentially stratified ocean, with application to SOFAR,” J. Acoust. Soc. Am. **55**, 220–226 (1974). [CrossRef]

**32. **C. Garrett and W. Munk, “Space-time scales of ocean internal waves,” Geophys. Fluid Dyn. **2**, 225–264 (1972).

**33. **C. Garrett and W. Munk, “Space-time scales of internal waves: a progress report,” J. Geophys. Res. **80**, 291–297 (1975). [CrossRef]

**34. **S. Flatté and R. Esswein, “Calculation of the phase-structure function density from oceanic internal waves,” J. Acoust. Soc. Am. **70**, 1387–96 (1981). [CrossRef]

**35. **J. Colosi and M. Brown, “Efficient numerical simulation of stochastic internal-wave-induced sound-speed perturbation fields,” J. Acoust. Soc. Am. **103**, 2232–2235 (1998). [CrossRef]

**36. **S. Flatté and G. Rovner, “Path-integral expressions for fluctuations in acoustic transmission in the ocean waveguide,” In *Methods of Theoretical Physics Applied to Oceanography*, P. Müller, ed., pp. 167–174 (Proceedings of the Ninth ‘Aha Huliko’a Hawaiian Winter Workshop, 1997).

**37. **S. Flatté and R. Stoughton, “Predictions of internal-wave effects on ocean acoustic coherence, travel-time variance, and intensity moments for very long-range propagation,” J. Acoust. Soc. Am. **84**, 1414–1424 (1988). [CrossRef]

**38. **J. Colosi, Ph.D. thesis, University of California at Santa Cruz, 1993.

**39. **J. Colosi*et al*., “Comparisons of measured and predicted acoustic fluctuations for a 3250-km propagation experiment in the eastern North Pacific Ocean,” J. Acoust. Soc. Am. **105**, 3202–3218 (1999). [CrossRef]

**40. **J. Colosi, F. Tappert, and M. Dzieciuch, “Further analysis of intensity fluctuations from a 3252-km acoustic propagation experiment in the eastern North Pacific,” J. Acoust. Soc. Am. **110**, 163–169 (2001). [CrossRef]