WaveQ3D: Fast and Accurate Acoustic Transmission Loss (TL) Eigenrays, in Littoral Environments

This study defines a new 3D Gaussian ray bundling acoustic transmission loss model in geodetic coordinates: latitude, longitude, and altitude. This approach is designed to lower the computation burden of computing accurate environmental effects in sonar training application by eliminating the need to transform the ocean environment into a collection of Nx2D Cartesian radials. This approach also improves model accuracy by incorporating real world 3D effects, like horizontal refraction, into the model. This study starts with derivations for a 3D variant of Gaussian ray bundles in this coordinate system. To verify the accuracy of this approach, acoustic propagation predictions of transmission loss, time of arrival, and propagation direction are compared to analytic solutions and other models. To validate the model’s ability to predict real world phenomena, predictions of transmission loss and propagation direction are compared to at-sea measurements, in an environment where strong horizontal refraction effect have been observed. This model has been integrated into U.S. Navy active sonar training system applications, where testing has demonstrated its ability to improve transmission loss calculation speed without sacrificing accuracy.


Introduction
Many underwater acoustic propagation models transform the three dimensional There are several ongoing efforts to deliver 3-D versions of Parabolic Equation 25;24 and Normal Mode 4 models. However, these models require large numbers of modes at frequencies above 1000 Hz, and this impedes their ability to provide real-time, active sonar results on low cost computer hardware. This paper develops a new acoustic transmission loss model using 3-D ray theory in geodetic coordinates. Unlike other efforts to model ray theory in geodetic coordinates, 23

Derivation
This section develops an implementation of 3-D refraction in spherical coordinates, interface reflection in this coordinate system, eigenray detection, and Gaussian ray bundle transmission loss.

3-D ray propagation in spherical coordinates
Ray theory is a high frequency approximation of the wave equation that decomposes the acoustic impulse response of the environment into surfaces of constant travel time (t) from the source ( Fig. 1.1). The rays are a vector field r that is normal to these surfaces at each point in space, and the path of these rays through the medium defines the direction of propagation. The fundamental equations of ray theory are derived by seeking power series solutions 21 to the Helmholtz equation.
where p( r) is the acoustic pressure as a function of location, r 0 is the source position, where α(t), β(t), and γ(t) are the r, θ, and φ components of ξ . The time derivatives ofr,θ, andφ are computed using a conversion into Cartesian coordinates.
Matching terms forr,θ, andφ yields a system of six scalar first-order differential equations.
where ∆t is the time step, and f is a vector of the positions, directions, and their derivatives.
When past values are cached instead of re-calculated, AB3 is much faster than other integrators with similar accuracy. However, because AB3 is not self-starting, a third order Runge-Kutta (RK3) algorithm 37 is used whenever the ray parameters must be initialized, or re-initialized as part of reflection. Interface reflection starts with a recognition that the next location in the iteration of Eq. (1.31) will put the wavefront location above the ocean surface or below the bottom. However, estimating the precise location of the point of impact requires an estimate of the point in time when the incident ray strikes the interface. Our derivation for estimating that time uses the geometry illustrated in Fig. 1.3 where I is the incident ray along directionÎ , R is the reflected ray along directionR, ζ is the incident grazing angle, andŝ is the surface normal. If the interface slope is nearly constant across the length of the incident ray, then the ratio of the time steps is equivalent to the ratio of the distances normal to the surface where h is the incident ray depth below the surface.
The direction of reflection uses the fact that incident and reflected rays share a common grazing angle. This means that the two vectors labeled B in Fig. 1  changes the depression/elevation angle of the ray path. The depression/elevation angle increases as the ray travels up slope, and decrease as it travels down slope. For each of these reflections, the ray path is also diverted away from the seamount peak in the horizontal. As shown in Fig. 1  targets. Fig. 1.6 is a side view of such a collision where r p is the position of a single acoustic contact, r njk is the position of a point on the wavefront, d njk is the distance from target to each point on wavefront, ∆t is the target offset along the direction of propagation, ∆µ is the target offset in the depression/elevation direction, and ∆ϕ is the target offset in the azimuthal direction.
A point on the wavefront r njk is the closest point of approach (CPA) for a specific target if it has the smallest distance to that target relative to the 26 wavefront points immediately surround it in the t, µ, and φ directions. The offset in each of these directions is computed by expressing d 2 p , the square of the distance to this target, as a second order Taylor series, relative to the CPA, in vector form.
where ρ is the target offset from CPA in vector form, g is the gradient of squared distance at CPA (3 elements), and H is the Hessian matrix of squared distance at CPA (3x3). One way to solve this equation would be to search for a value of ρ for which Eq. (1.48) was zero. However, since d 2 p is positive definite, searching for the minimum value of d 2 p , indicated by a zero in the first derivative, also solves this problem. Some eigenray products are computed directly from this offset vector.
where t p is the travel time from the source, µ p is the depression/elevation launch angle at source, and ϕ p is the azimuthal launch angle at source. The direction at the target is found by forward solving a 2nd order Taylor series for ξ in the neighborhood of the CPA. The eigenray data is completed by the calculation of transmission loss, which is discussed in the next section.
The computation of d 2 p in spherical coordinates is much less efficient than an equivalent calculation in Cartesian coordinates. However, the impact of this difference is minimized when the number of targets is small compared to the number ray tracing points. This is a good assumption for sonar training systems, but it makes our model much slower than existing methods for calculating transmission loss over millions of range, depth, and bearing combinations in tactical decision aids.

3-D Gaussian ray bundles
In conventional ray theory, the spreading of acoustic transmission loss is estimated by measuring the changes in ensonified area between ray paths. The intensity across the wavefront is inversely proportional to the change in a surface area segment compared to its area at the source. 21 The Gaussian beam approach uses a parabolic equation approximation normal to each ray to compute spreading in the form of a Gaussian profile for each ray. 9;35;34;5 Gaussian ray bundling models, such as CASS/GRAB, use the distance between rays (from conventional ray theory) to define the size of the Gaussian profile for each ray.
In 2-D Gaussian beam models, the intensity at the target location is a summation of contributions from rays above and below the eigenray target. To extend Gaussian ray bundles to three dimensions, we use an approximation that computes independent factors in the µ and ϕ directions ( Fig. 1.7).
where G( r p ) is the total Gaussian ray bundle intensity at the eigenray target, (j, k) are the index numbers of the cell containing the eigenray target, g j are the Gaussian ray bundle contributions along depression/elevation direction, g k are the Gaussian ray bundle contributions along the azimuthal direction, 2J + 1 is the number of beams used in the depression/elevation summation, and 2K + 1 is the number of beams used in the azimuthal summation. Computing a product of independent factors is equivalent to assuming that the width of the Gaussian ray bundles in the µ direction represents a local average in the ϕ direction, and vice versa.
The intensity of each Gaussian ray bundle contribution is a function of the width of each beam and the distance along the wavefront to the eigenray target, normalized to the average distance across the beam at the source where w j and w k are the half-widths of the Gaussian ray bundle in the µ and ϕ directions, and d j and d k are the distances in the µ and ϕ directions from the Gaussian ray bundle center to the target.
GRAB 45 models the frequency dependent component of the beamwidth by giving each beam a minimum width.
where λ is the wavelength of the signal being modeled, w j is the half cell width of beam j, and w j (f ) is the adjusted width of beam j. GRAB models beams centered on each ray and then between each ray to create a minimum overlap of 50% between Gaussian ray bundles.
A physical interpretation of Eq. (1.63) is that the λ term is the frequency dependent Gaussian spreading that GRAB expects for rays that are infinitely close together.
The w j terms can be interpreted as the Gaussian width created by discreetly sampling the launch angles. Instead of using the maximum of these two contributions, our model convolves these two sources of spreading and adds their Gaussian widths as the sum of squares.
(w j (f )) 2 = (2w j ) 2 + (2πλ) 2 (1.64) Eq. (1.64) produces results that are similar to (1.63), but there is a smooth transition between the domains dominated by the w j and λ terms. The factor of 2 in w j has been artificially introduced to produce the same 50% overlap as GRAB without doubling the number of beam calculations. Normalizing Eq. (1.59) and (1.60) by the combined effect of both spreading sources conserves energy across the wavefront.

Test Results
Although the primary goal of this research is to develop a transmission loss model that reduces computational burden, the training systems also require an accurate representation of real-world phenomena. This section presents the results of several key accuracy tests including for ray path refraction accuracy using a Munk profile (Section 3.1), Gaussian beam projection into the shadow zone for an n 2 linear profile (Section 3.2), and horizontal refraction from a 3-D analytic wedge (Section 3.3). An comparison to CASS/GRAB executions times is provided in Section 3.4.

Refraction accuracy benchmark
Because our model's calculation of transmission loss is closely tied to the location of ray paths, refraction accuracy is an important element of its overall accuracy.
Because the use of spherical coordinates incorporates the radius of the Earth (a large number) into the radial coordinate in Eqs. 1.31 through 1.36, we need to address concerns that our approach would suffer from numerical accuracy problems.
To evaluate the refraction accuracy, we compare an analytic solution for the Munk profile 28 representation of a deep sound channel to ray paths generated by our model.
The version of the Munk profile used in this test is defined by Eqs. (1.65) and Eq. (1.66) where z is the normalized depth (positive is down), z 1 is the depth of the deep sound channel axis (1300 m), B is a depth scaling factor (1300 m), c 1 is the sound speed on deep sound channel axis (1500 m/s), and is the profile scaling factor (7.37x10 −3 ).  z t is the target depth. Although these integrals only apply between the source and the first vertex or reflection, paths out to any range can be constructed by repeating this process after the vertex or reflection.
Pekeris' modified index of refraction, shown in Eq. (1.70), allows Cartesian models, like the Munk profile, to incorporate earth curvature effects into their calculations. 33 This process can also be inverted, as shown in Eq. (1.71), to allow spherical models, like ours, be compared to flat earth benchmarks.

2-D transmission loss benchmark
To compare the accuracy of our transmission loss estimates to existing 2-D Gaussian beam models, coherent transmission loss is calculated at the edge of a shadow zone     results. Prior to the shadow zone, all three models produce similar results. Although the loss predicted by our model is slightly higher than FFP in the shadow zone, its results are similar to the ones produced by GRAB. In both cases, the slight rise in the coherent transmission loss appears to be caused by phase inaccuracies in the multi-path travel times predicted in the shadow zone.

3-D transmission loss benchmark
To demonstrate 3-D effects in transmission loss, we examine the analytic solutions for a 3-D wedge. In this scenario, receivers are at the same distance from the wedge apex as the source, but offset in range across the slope. In an 2-D model, these receivers appear to exist in an environment of constant depth. Because the 3-D solution horizontally refracts acoustic energy down the slope, the 3-D solution has higher transmission loss, as a function of range across the slope, than the 2-D model. Using the method of images, we assume that each reflection gives rise to a source image, and that these images lie on a circle centered on the apex of the wedge. This derivation is similar to the Deane/Buckingham model 13 , but it simplifies that model by assuming that interface reflection coefficients are limited to 1. This simplification creates an analytic solution that is accurate at higher frequencies.  is the index number for each receiver; r q is the location of each receiver; R n,m,q is the slant range from each source image to each receiver; c is the speed of sound in water; f is the signal frequency; k is the acoustic wave number = 2πf /c; and p q is the total complex pressure for each receiver.
To compute R n,m,q , we define a cylindrical coordinate system whose axis travels along the wedge apex: R s is the slant range of original source from the wedge apex; α s is the angle of original source down from the ocean surface; α n,m is the angle of each source image, relative to the ocean surface, negative if above surface; R q is the slant range of each receiver from the wedge apex; α q is the angle of each receiver down from the ocean surface; and y q is the cross-slope distance of each receiver relative the vertical source/origin plane.
An inspection of the geometry in Fig. 1.12 allow us to compute α n,m and R n,m,q for the 3-D wedge.
Source images outside of the range α n,m ∈ [−π, π] result in imaginary images that contribute to the diffracted component of the acoustic field. For small wedge angles and locations far from the apex, the diffracted components are negligible and need not be considered. 13 An equivalent solution for an environment of constant depth uses an alternate version of R n,m,q .
where x s , z s is the range and depth of original source relative to ocean surface; x q , z q is the range and depth of each receiver relative to ocean surface; y q is the cross-slope distance of the receiver relative the vertical source/origin plane; z n,m is the depth of each source image; and z w is the water depth; Because our model uses geodetic coordinates, the simple wedge used in our analytic solution can only be approximated in this comparison. On a round earth, an interface with constant slope is a curved surface instead of a plane. To minimize the impact of this curvature, the wide wedge angle scenario is used to shorten the range over

Computational efficiency
A key premise of this paper is that, when target/sensor geometries are constantly evolving, it is more computationally efficient to perform acoustic transmission loss in the latitude, longitude, altitude coordinates of the underlying environmental databases, than it is to convert the 3-D environments into a series of Nx2-D radials. To evaluate Cross Slope Range (km)  targets. For 100 targets, our model is over 6 times faster than GRAB and 10 times faster than the speed of sound.

Acknowledgments
Testing and productization of this model were funded by the High Fidelity Active

Abstract
Acoustic transmission loss measurements from the Calibration Operations (CALOPS) experiment for the Shallow Water Array Performance (SWAP) program included horizontally refracted returns that were as much as 30 degrees away from the true bearing between source and receiver. In many cases, the in-shore refracted path was as much as 20 dB stronger than the true bearing path. In this study CALOPS transmission loss measurements at 415 Hz are compared to predictions from a 3D Gaussian Ray bundling model. The geoacoustic model that provides good modeldata comparison is consistent with the geologic and sediment core data collected at the location but differs slightly from the bottom model used at lower frequencies (206 Hz and 52.5 Hz) in a previous study.

Introduction
Several investigators have recently studied the presence of strong 3D propagation effects in experimental data on the continental shelf in the Florida Straits area. 30;22;41;19;20;4 Acoustic transmission loss measurements from the Calibration Operations (CALOPS) experiment for the Shallow Water Array Performance (SWAP) program included horizontally refracted returns that were as much as 30 degrees away from the true bearing between source and receiver. 19 In many cases, the in-shore refracted path was as much as 20 dB stronger than the true bearing path. CALOPS transmis-  The measured signal level is extracted from the frequency spectrum peak for each

Comparison of measurements with model predictions
The environmental conditions and geometry discussed in Section II were used for the 3D modeling of the acoustic propagation for comparison with the measurements.      It should be noted that the 415 Hz measurements used in this study may be more sensitive to surficial sediments than measurements at 52.5 Hz and 206 Hz. This sensitivity difference could explain our need to change the effective location of the boundary between sand and limestone bottom loss provinces.

Abstract
This paper describes the migration process undertaken by the US Navy to migrate from legacy, large footprint mainframe computer-based sonar simulation systems to next-generation sonar simulation systems with a smaller footprint, lower costs and better accuracy than the legacy models. The paper will describe the development efforts to create a faster and more accurate acoustic transmission loss (TL) and reverberation model for sonar simulation/stimulation systems in littoral environments based on ray theory for active sonar frequencies (above 1000 Hz). The paper also describes how the next-generation model augments ray theory with Gaussian beam techniques (based on the Gaussian Ray Bundling or GRAB), which enables simulation of frequencies as low as 150 Hz. The paper will detail the integration challenges faced by the US Navy to migrate from the legacy models to the next-generation sonar simulation model into the Navy's Live Virtual Constructive Modeling and Simulation (LVCMS) product line that includes PACT3, BATTT, and EFAAS simulators/simulations. The paper will also describe the results of these integration efforts, including the ability to provide trainees with improved training via more complex scenarios in the LVCMS training suite without increasing their hardware costs or footprint.

Introduction
The when they are in port, they train for the next operation. The Navy Continuous Training Environment (NCTE) is the hardware/software architecture that brings together the training systems for individual ships, submarines, aircraft, and commanders into a common virtual environment. Each combat system builds their own training system, and the NCTE integrates these systems across wide area networks.
Coordinated anti-submarine warfare (ASW) depends on a complex mosaic of diverse capabilities in a highly variable physical environment. 27 No single ASW plat-

Wavefront Queue 3-D Model
The lost. The amount of data to be analyzed is much smaller than in the TDA case, but predictions must be generated in real time for a dynamically evolving scenario.
The new model takes a new approach by leaving the data in latitude, longitude, and depth coordinates, and changing the transmission loss equations to operate in this 3-D coordinate system. The oceanographic data for large areas is cached in memory and re-used for multiple sensors. As illustrated in Figure 3.3, oceanographic

Speed and Accuracy Testing
To evaluate the speed potential of the new model, a scenario was developed to compare speeds to existing Navy models for high fidelity TDAs. One challenge in developing such a scenario was finding a case in which the speed of the TDA model This test was modified to create a variable number of targets, 100 km from the source, evenly spaced in azimuth. An equivalent scenario was then created using WaveQ3D.
The timing tests were run on a Dell Latitude Laptop with an E6520 Intel i5-2540M CPU running at 2.60GHz. As shown in Figure 3  • Verification relative to analytic solutions for simple conditions, • Verification relative to results of other standard models, • Validation relative to at-sea measurements, • Validation by sonar training experts as part of integration into actual training systems.
One of the key tests (Section A.5.5) is a classic benchmark that looks at the accuracy of transmission loss propagation into an acoustic shadow zone. 32 In this scenario, illustrate by the left side of Figure 3.5, a strong downward refracting environment limits the amount of energy reaching distant targets. The right side of  The resulting audio is then sent back to the Acoustic Server, which forwards it to the sonar for processing and presentation to the trainee. Integrating the new model into the AP module required the following changes: • A method was added to the Dynamic Ocean Generator (DOG) that allowed clients to request information for the whole gaming area in latitude, longitude, depth coordinates. This was already the native coordinate system of the data, but the DOG had no request mechanism for this format.
• A singleton Environment Manager was created to store and manage environmental data, in a WaveQ3D compatible format, for the whole gaming area.
Multiple processing threads, for multiple sensors, can share this single representation of the environmental data.
• The function that creates the requests for transmission loss was modified to create a single request for all targets, instead of a separate request for each target.
• WaveQ3D specific configuration options were added to the AP modules setup process.
Most of the engineer effort focused on fitting the new model into the existing architecture, and then testing the accuracy and robustness of the integration.  To quantify the differences between the two models, difference in the transmission loss, travel time, and source/receiver angles of each eigenray were measured and analyzed. 8 Table 3.1 presents the differences between the WaveQ3D and FeyRay models, for the scenario illustrated in Figure 3.8. Our analysis indicates that differences between the two models were primarily due to differences in the sound velocity profiles in the 2-D and 3-D environments. • Development of a reverberation model based on WaveQ3D, • Replacing the legacy model in all LVCMS/ASW VAST signal processing modules, • Continuing to support accuracy testing, and

Introduction
Acoustic tomography uses propagation models to predict the multi-path arrival structure of one-way transmissions received from a distant source. 40 Tomographic inversion repeatedly perturbs the ocean sound speed and re-computes the arrival structure until only trivial differences remain between the measured and modeled results. WaveQ3D is a new 3D Gaussian ray bundling model that is designed to reduce the computation burden of propagation loss calculations. This feature may have significant benefit in tomographic inversion, where thousands of modeling runs may be required to converge on a solution. However, speed is only a benefit if the model accurately predict the multi-path arrival structure for a given profile. This paper investigates WaveQ3D's ability to predict the multi-path arrival structure for the Munk profile, an idealized representation of deep sound channel conditions in the North Pacific. 28 where z is the normalized depth (positive is down), z 1 is the depth of the deep sound channel axis (1300 meters below ocean surface), B is a depth scaling factor (1300 meters), c 1 is the sound speed on the deep sound channel axis (1500 m/s), and is a profile scaling factor (7.37x10 -3 ).

Benchmark solution for Snell's Law in spherical media
WaveQ3D computes propagation on a spherical section of the Earth using geodetic coordinates: latitude, longitude, and altitude. This approach speeds up 3D applications by using environmental parameters and their derivatives directly in the latitude, longitude, and depth directions without reducing the problem to a series of Nx2D Cartesian projections. Snell's Law for spherical media (Eq. where η is the depression/elevation angle along the ray path; m is the ray parameter for spherical media (constant for each launch angle), τ is the travel time along the ray path; ρ s is the radial coordinate for the source depth; and ρ t is the radial coordinate The eigenrays are the discreet set of paths that connect the source to a target at the same depth. To find eigenrays for the benchmark solution, we break the each propagation path into a series of source depth crossings, and calculate the range at which each crossing occurs. After an even number of crossings, paths that left the source traveling up will arrive at their target from below. At the end of each cycle, the ranges for each path increase monotonically as a function of increasing launch angle.
An eigenray is present when these ranges span the target location. Interpolation allows us to invert range as a function of launch angle into launch angle at the target range. The same process is used for Eq. 4.6, the travel time integral . To find all eigenrays, the process is repeated for odd numbers of crossings, and paths that leave the source traveling down.

Equivalent benchmark solution in Cartesian coordinates
To evaluate the impact of numerical error in the discreet evaluation of Eqs.
where f is frequency of signal; c(ρ) is speed of sound as a function of depth in this coordinate system. We define an equivalent flat Earth using a change of variables (x, z).
where R is radius of the earth; x is range as change in distance along the surface of the earth; z is depth below the ocean surface (mean sea level); and c (z) is speed of sound as a function of depth in this coordinate system. To model acoustic propagation in the (x, z) system, we use another change of variables (x, u) to transform Eq. 4.7 into a Cartesian form.
where ρ is depth non-linear function of ρ; c (u) is speed of sound as a function of depth in this coordinate system. The relationship between the (x, z) and (x, u) coordinate systems is given by The relationship between z and u is often approximated.
where u 1 is the first order approximation; and u 2 is the second order approximation.
The analysis solution in (x, u) terms is implemented using Snell's Law in Cartesian coordinates.    Transform. Note that the results for "2nd order" are so close to the "no approx"

WaveQ3D accuracy
To     Step (

Conclusions
For a target at 200 km, the differences between WaveQ3D and the benchmark solu-

A.1 Abstract
The theory paper defined a new undersea acoustic propagation loss model that is specifically designed to support real-time, sonar simulation/stimulation systems, in littoral environments, at active sonar frequencies. This paper seeks to verify the model's implementation by comparing the modeled results to analytic solutions.

A.2 Introduction
The theory paper 38  This follow-up paper analyzes the capabilities and limitations of this implementation.
The Capability Maturity Model Integration (CMMI) 11 separates testing needed for such an effort into two phases: • Verification testing ensures that selected work products meet their specified requirements. For this analysis, the WaveQ3D model is decomposed into its component parts (such as ray tracing, reflection, eigenray finding, and propagation loss), and the results from each parts are compared to analytic solutions.
• Validation testing demonstrates that a product or product component fulfills its intended use when placed in its intended environment. For the WaveQ3D model, this will consist of comparisons to real-world results in a subsequent paper.
Decomposing the testing in this way is designed to ensure that any conclusions drawn from the modeled results rest on a firm foundation of understanding.

A.3 Ray Tracing Tests
The ray paths in this model use a third order Adams-Bashforth (AB3) marching solution 49 to create a circular queue of time domain wavefronts. The evolution of the wavefront shape, as it passed through the 3-D ocean environment, is defined by the following equations.
where c is the speed of sound as a function of location; t is the travel time; (r, θ, φ)

A.3.1 Comparisons to "flat earth" benchmarks
Because WaveQ3D is one of the few models to use non-Cartesian coordinate system, comparisons to other work often require translation before differences can be analyzed. This section analyzes the accuracy of the translation between spherical and Cartesian coordinate models. Cartesian coordinate propagation models frequently use a modified index of refraction 33 to incorporate earth curvature effects into their calculations.
where R is the radius of earth's curvature in this area of operations; r is the radial distance from the center of curvature (positive is up); z = R − r is the below the ocean surface (positive is down); n(z), n (z) are the original and modified index of refraction, and c(z), c (z) are the original and modified speed of sound. When a testing benchmark is specified in Cartesian coordinate, the inverse of this process must be implemented to create an equivalent environment in spherical earth coordinates.
where c(r) is the benchmark's original c(z) sound speed converted to a function of r.
The Munk profile 28 was used to evaluate the accuracy of testing Cartesian benchmarks in a model based on spherical earth coordinates. The Munk profile is an idealized representation of a typical deep sound channel, and it was chosen for its ability to support long range paths without interface reflection.
where z is the normalized depth (positive is down), z 1 is the depth of the deep sound channel axis (1300 meters), B is a depth scaling factor (1300 meters), c 1 is the sound speed on deep sound channel axis (1500 m/s), and is the profile scaling factor The second case is more interesting for this study, because it represents a widely accepted (but seldom mentioned) lower limit on ray path range accuracy, approximately 0.03% of the total range.

A.3.2 Ray path accuracy in a deep sound channel
In this test, the refraction accuracy of the WaveQ3D model was computed for a With a 100 ms step size, the WaveQ3D result deviates from the analytic solution by a maximum of -8.62 meters at cycle range of 129.95 km (0.007% error). However, 50 out of 58 samples (86%) exhibited errors less than ± 2 meters (0.002% error). Ray paths that were initially launched toward the surface, where the sound speed gradient is highest, had consistently larger errors than paths that were launched down. The fact that these errors were all significantly less than those of associated with the earth curvature correction (Fig. A.2) suggests that WaveQ3D model's cycle range estimate meets or exceeds the accuracy of Cartesian models used on a spherical earth.
To estimate the impact of WaveQ3D options on this result, the maximum error was also computed as a function of time step size. The circles in Step sizes as large as 150 yielded results that were at least as accurate as errors associated with the modified index of refraction. Given that some environments may have stronger gradients than the Munk profile, we conclude that a 100 ms step size should be adequate for most long range applications. there appeared to be a slight bias. We discovered that many of the launch angles had cycle range errors right around -0.5 m.
The sensitivity of these results to time step size is illustrated in Fig. A.7. The 150 ms step size had errors of about 2 meters. Errors for smaller time steps approached zero as expected. Larger time steps increase as a power law until the step size was larger than 350 ms, where the error quickly grew to hundreds of meters. This sudden change in error appears to be the result of an inability of model to properly sample the sound velocity profile field at the larger step sizes. From this, we conclude that our earlier 100 ms time step recommendation continues to be valid for this case.

A.3.4 Ray path accuracy along great circle routes
In this test, an ocean with a small amount of downward refraction was used to verify the WaveQ3D model's ability to propagate rays along great circle routes. Acoustic , where χ(t), φ(t) are latitude and longitude as a function of time; and ϕ a (t) is the analytic solution for the azimuthal launch angle for a target at χ(t), φ(t). The rays traveled at constant depth with a maximum deviation from ϕ a (t) of 2.81x10 -10 degrees. This level of accuracy should be more than adequate for most applications.

A.4 Interface Reflection Tests
The WaveQ3D model estimates the partial time step during a 3-D interface collision using the equation whereŝ is the surface normal;r is the unit vector in the radial direction; h is the incident ray height above bottom; δt is the time step needed to reach the interface; and dr dt is the radial ray tracing component defined by Eq. (A.1). The direction of the 3-D reflection is computed usingR whereÎ is the incident ray path direction; andR is the reflected ray path direction.
The WaveQ3D model also applies a second order Taylor where R is the radius for the ocean surface; D is the water depth; R b = R − D is the radius for the ocean bottom; ∆θ is the latitude change between surface bounces; L

A.4.2 Reflection accuracy with a sloped bottom
This test looked at the ability of WaveQ3D to predict the direction of reflection from a bottom that has a 1 degree up-slope in the latitude direction. At each bottom reflection, the depression/elevation angle of the WaveQ3D ray path should decrease by 2 o , just like the analytic result. Launching a ray at a the same depression/elevation angle as the last test (-5.183617057 o , down) produced the results shown in Fig. A.11.
Note that the time step for this test was set to 1 ms to make it easier to compute depression/elevation angle just before, and just after, each collision. The maximum deviation of any of the three bottom reflections from their analytic reflection direction result was 1x10 -5 degrees. From this, we conclude that WaveQ3D produces acceptable  (Fig. A.11) expected for sloped bottoms. What was new in this test was the face that ray paths were reflected into a new azimuthal direction each time that they interacted with the bottom. These out-of-plane reflections resulted in a down slope ray path that was offset by more than 21.9 km from the up slope path, after 14 bounces off of the bottom. From this, we conclude that WaveQ3D results will have a significant contributions from out-of-plane ray paths whenever there are multiple interactions with complex bottom bathymetry features.

A.5 Eigenray and Propagation Loss Tests
The tests discussed in this section compare WaveQ3D's eigenray and propagation loss calculation to analytic solutions.
The WaveQ3D eigenray estimation process establishes the relative geometry be-tween rays paths and targets. That geometry then allows the calculation of travel time (t), source angle (µ), and target angle (ϕ) for each multi-path arrival. WaveQ3D computes these eigenray products by searching for the offsets which minimize the squared distance from targets to points on the wavefront. Once the CPA has been determined, the 3-D offset to the target is computed using where ρ is the offset from CPA in vector form; d 2 is the squared distance from each point on the wavefront to the target; g is the gradient of squared distance, evaluated at CPA (3 elements), and H is the Hessian matrix of squared distance, evaluated at CPA (3x3).
The propagation loss at the target location is a summation of contributions from the rays that surround the eigenray target. To create a 3-D acoustic field across the wavefront, WaveQ3D uses independent Gaussian beams in the µ and ϕ directions and ignores the cross terms. WaveQ3D treats the beam width calculation as a convolution between Weinberg's frequency dependent "minimum width" term 45 and a second Gaussian that represents the spatial spreading created by the sampling of the wavefront.
(A. 35) where λ is the wavelength of the signal being modeled; w j,k is the cell width of beam j or k, and w j,k (f ) is the adjusted width of beam j or k. The factor of 2 in the w j,k term creates a minimum overlap of 50% between neighboring beams. Normalizing

A.5.1 Eigenray accuracy for a simple geometry
This test constructs a short range geometry in which travel time, source angle, and target angle can be computed analytically for direct path, surface reflected, and bottom reflected paths on a spherical earth. The geometry for this test is illustrated in

A.5.2 Eigenray accuracy for Lloyd's mirror on spherical earth
On a flat earth, the Lloyd's mirror geometry generates exactly two paths: a direct path and a surface reflection. However (as shown in Fig. A.14) isovelocity ray paths actually form unexpected caustics when the curvature of the earth is incorporated.
where t s is the surface-reflected travel time from source to target; µ s,source is the surface-reflected depression/elevation angle at source; and µ s,target is the surfacereflected depression/elevation angle at target.  was produced at each target location (see Fig. A.14). This choice allows the test to be run with a ray spacing that was much more typical than the fine scale D/E angles used in the previous section.
where r is the target range; z is the target depth (positive is down); z s is the source depth (positive is down); L 1 is the slant range to source; L 2 is the slant range to image source (above water); k is the acoustic wave number (2πf /c); and f is the signal frequency; c is the speed of sound; p(r, z) is the complex pressure.  To compare these results quantitatively, we will use a set of statistical measures In this test, the eigneray and propagation loss accuracy of the WaveQ3D model were analyzed for the extreme n 2 linear test case developed by Pedersen and Gordon. 32 Two geometries were supported in this test: • The shallow source geometry puts the source at a depth of 75 m and creates a series of targets at a depth of 75 m with ranges from 500-1000 m.  (Because this geometry does not support the formation of a caustic, neither WaveQ3D nor GRAB applies a −π/2 phase shift to this path.) Rays above the critical angle hit the surface before ensonifying the target. Both ray families included targets that were ensonified in their evanescent region, the region outside of the ray fan.   nique. 14;6 Note that, in all regions, this implementation of FFP was consistent with an ideal wave equation solution, except for the presence of some implementation jitter in the ranges above 880 m. In the region between 500-750 m, all three models produced almost identical results. As target passed into the shadow zone region, WaveQ3D and GRAB produced values that were similar to each other, but slightly higher than FFP. But we also noted that, if a coarser set of depression elevation launch angles was used for WaveQ3D, the surface reflected path would have disappeared prematurely, which would have manifested as shadow zone oscillations in the total propagation loss. Taken as a whole, we feel confident that Wave3D propagation loss errors are similar to those of GRAB for this scenario. Deep source geometry The big difference between the models appears to be the fact that the WaveQ3D transmission loss has a more gradual roll-off into the shadow zones than GRAB.
Quantitative comparisons for the deep source's individual eigenrays are summarized in Table A.7. The statistics for the WaveQ3D angles for the caustic path are skewed by  Table A.8. In the region from 3000 to 3040 meters, the FFP has stronger oscillations than either the WaveQ3D or GRAB models. This can be attributed to surface reflected contribution that is stronger in the FFP result than in either of the Gaussian beam models. At all other ranges, the WaveQ3D and GRAB results were both a very close fit to FFP's total propagation loss.
A WaveQ3D anomaly was discovered during the deep source geometry testing.
Errors in the total propagation loss within the shadow zone were frequently seen when the launch angles finer than 0.01 degrees were used. What we discovered was that at this spacing, the contribution in the shadow zone was the result of a summation of over 100 Gaussian beams for both the direct and caustic paths. Small errors in the   this anomaly results in a slight broadening of the transmission loss decay tails. But, since the shape of the outer edge of the total propagation loss depends on destructive interference between the paths, these errors often resulted in imperfect cancellation, on the order of -70 dB, in the region between 3090 and 3105 meters. Although these problems may have been aggravated by the extreme environment, developers should use extra caution when using WaveQ3D with super-fine ray spacings.

A.6 Derivations
A.6.1 Ray path derivation for concave ocean surface The analytic solution for the direct-path eigenrays was derived from the laws of sines and cosines. where t d is the direct-path travel time from source to target; µ d,source is the directpath depression/elevation angle at source; and µ d,target is the direct-path depression/elevation angle at target.
The analytic solution for the surface reflected solution also starts with the law of cosines.
This can be reduced to a simpler form using R 2 + b 2 n − r 2 n = R 2 + (R 2 + r 2 n − 2Rr n cos(∆θ n )) − r  Because sonar detection ranges are often much longer than the depths of interest, there is frequently a need to emphasis propagation paths near the horizontal. The GRAB model 17 manages this requirement by automatically adding rays in the de-pression/elevation directions needed to support caustics, surface ducts, and SOFAR channels. The WaveQ3D model currently does not support any such automatic ray adjustment; but it does support any ray spacing that is monotonically increasing.
Tangent spaced depression/elevation angles will be used frequently in the eigenray and propagation loss testing for WaveQ3D. As shown in the top panel of Fig

A.6.3 Propagation loss error statistics
We would like to define quantitative differences between modeled and analytic propagation losses in a way that illuminates the suitability of the model for real-time, sonar simulation/stimulation systems. To that end, we define the following statistical This testing also led to the development of a few general guidelines for the use of the WaveQ3D model.
• Using a time step size as coarse as 100 ms seems to produce accurate results for most applications. However, this step size should be decreased to 10 ms for applications that need high accuracy at ranges less than 2 km. This is particularly true for applications that are interested in effects near the ocean surface, or directly below the source.
• The selection of launch angles can have a significant effect on propagation loss results. Like all models that as based on ray theory, WaveQ3D can miss features in the environment because of spatial under-sampling. Tangent spaced depression/elevation angles appear to improve WaveQ3D model performance for scenarios dominated by horizontal paths. Uniform spacing is suggested for applications that need high accuracy at ranges less than 2 km.
• The WaveQ3D model is designed to be computationally efficient when computing transmission loss for up to 100 targets, at multiple frequencies, in a fully 3-D environment. But because the WaveQ3D eigenray detection process is less efficient than an equivalent calculation in Cartesian coordinates, WaveQ3D can actually be much slower than other models when thousands of range/depth combinations are required. Applying the WaveQ3D model in 2-D environments is also not very efficient.