## S V nL expt kor

which is the asymptotic expansion for large arguments. The equation for V(r,z) [Equation (4.25)] can then be simplified to:

d r2

Further assume that: d2V dV

d r2 d r which is the paraxial approximation. Then, Equation (4.28) reduces to:

d2V dV

d z2

which is the parabolic wave equation. In this equation, n depends on depth (z), range (r) and azimuth (0). This equation can be numerically solved by "marching solutions" when the initial field is known (e.g. Tappert, 1977). The computational advantage of the parabolic approximation lies in the fact that a parabolic differential equation can be marched in the range dimension whereas the elliptic reduced wave equation must be numerically solved in the entire range-depth region simultaneously. Typically, a Gaussian field or a normal-mode solution is used to generate the initial solution. Additional methods for generating accurate starting fields have been described by Greene (1984) and by Collins (1992). Collins (1999) improved the "self-starter" (a PE technique for generating initial conditions) by removing a stability problem associated with evanescent modes.

### 4.7.2 Numerical techniques

Existing PE models employ one of four basic numerical techniques: (1) split-step Fourier algorithm; (2) implicit finite-difference (IFD); (3) ordinary differential equation (ODE); or (4) finite element (FE). The original split-step algorithm developed by Tappert (1977) solves the parabolic wave equation by imposing an artificial zero bottom-boundary condition and pressure-release surface condition.

At the time Tappert introduced the PE method to the underwater acoustics community, there was a critical need for a capability to predict long-range, low-frequency sound propagation, as would occur in the vicinity of the sound channel axis. Since this type of propagation is characterized by low-angle, non-boundary interacting energy, the PE method was ideally suited to this purpose. Thus, the first introduction of the PE method to naval sonar applications was a celebrated event. Subsequent work focused on making the PE method more robust so that it could be applied to a wider range of problems in underwater acoustics.

While the split-step algorithm is an efficient method for solving a pure initial-value problem (Perkins et al., 1982), several difficulties arise when there is significant interaction with the sea floor (Bucker, 1983). These difficulties are due to density and compressional sound-speed discontinuities at the water-sediment interface, strong gradients of the compressional sound speed in the sediment layers and rigidity in the sediment layers that produces shear waves. Thus, when bottom interaction is strong, a more generalpurpose solution is desirable. For this reason, implicit finite-difference (IFD) schemes (Lee et al., 1981; Robertson et al., 1989) and ordinary differential equation (ODE) methods (Lee and Papadakis, 1980) have also been developed to solve the parabolic equation. Lee and McDaniel (1987) provided an in-depth development of the IFD technique. Their development was made more useful by the inclusion of benchmark test examples together with a listing of the computer program (comprising approximately 1,700 lines of FORTRAN code). A microcomputer implementation of the IFD model was reported by Robertson et al. (1991), who included a listing of their computer program. Collins (1988a) described an FE solution for the PE.

Recent modeling developments utilizing parabolic approximations have been directed at refining and expanding the capabilities of existing techniques. Useful summaries of this progress are available in the literature (Lee, 1983; Scully-Power and Lee, 1984; Lee and Pierce, 1995; Lee et al., 2000). Actual comparisons of computer model results for a set of four ocean acoustic test cases was documented by Davis et al. (1982). A more recent comparison of PE models involving seven test cases was described by Chin-Bing et al. (1993b). One of the most interesting test cases described by Chin-Bing et al. (1993b) involved a leaky surface duct. So-called "ducted precursors" are generated by modal energy in the surface duct that leaks out and travels as CZ or bottom-bounce paths or both, before coupling back into the surface duct down range. The importance of this problem is that a small phase error in the refracted (leaky) path can produce large changes in the predicted sound level in the duct beyond the first CZ (Porter and Jensen, 1993). This pathology was particularly evident in some wide-angle, split-step PE models. Yevick and Thomson (1994) explored this problem further and formulated a new propagation operator that retained the computational efficiency of the split-step algorithm but which was more accurate.

Traditionally, PE models have been applied to anelastic ocean-bottom regions by treating the shear waves as an additional loss mechanism. This approximation breaks down when anelastic propagation becomes significant, as in shear conversion due to backscattered acoustic fields. In general, PE models propagate the acoustic field only in the forward direction, thus excluding backscatter. A two-way PE model (Collins and Evans, 1992) was developed to improve the computation of backscattered energy for use in reverberation simulations (e.g. Schneider, 1993). A modified version of this two-way PE model (called spectral PE) treats backscattering problems in 3D geometries (Orris and Collins, 1994). Lingevitch and Collins (1998) argued that a poro-acoustic medium is, in fact, the limiting case of a poro-elastic medium in which the shear wave speed vanishes. Collins and Siegmann (1999) extended energy-conservation corrections from the acoustic case to the elastic case.

By using operator formalism, the parabolic equation can be derived in another fashion (Lee and McDaniel, 1987: 314-15). By assuming commutation of operators:

d r2

128 Propagation II: mathematical models (Part One) can be factored as:

By considering only the outgoing wave, Equation (4.33) reduces to:

which can be rewritten as d

where

The square-root operator (Q) can be approximated using a rational functional representation:

where

Using Equation (4.37), Equation (4.35) can be rewritten as:

For [AB C D] = [111 0], Equation (4.39) reverts back to that of Tappert [Equation (4.30)]. For [A B C D] = [1 4 1 11, the form attributed to Claerbout (1976) is obtained.

### 4.7.3 Wide-angle and 3D adaptations

Tappert's (1977) original split-step PE method handled small-angle (< 15°) propagation paths, a restriction imposed by the basic paraxial assumption. Unless the split-step algorithm is modified for wide angles (e.g. Thomson and Chapman, 1983), large phase errors may be introduced into the solution. In order to extend the maximum angle of half-beamwidth propagation in the PE model, alternate forms of the square-root operator have been explored. For split-step PE formulations, Thomson and Chapman's (1983) wide-angle approximation has been used. The Claerbout (1976: 2o6) wide-angle approximation is used in some finite-difference PE formulations with values of [A B C D] = [1414] in Equation (4.39). Other approximations of the square-root operator have been explored including higher-order Pade forms, which have been incorporated in finite-difference and finite-element PE formulations to achieve near 9o° half-beamwidth propagation (Chin-Bing etal., 1993b). Thomson and Wood (1987) investigated a post-processing method for correcting phase errors in the parabolic approximation approach. This method was later extended by Thomson and Mayfield (1994) to include an exact, non-local boundary condition at the sea floor. Since strong boundary (surface and bottom) interactions are associated with wide-angle propagation, the accurate treatment of irregular interfaces and rough boundaries assumed greater importance in PE models (Bucker, 1983; Lee and McDaniel, 1983; Dozier, 1984; Collins and Chin-Bing, 199o; Brooke and Thomson, 2ooo). Adaptations to under-ice environments have also been of interest to sonar modelers.

Extensions to 3Ds are generally implemented in an approximate manner (e.g. N x 2D), although the primary parabolic wave equation is fully 3D (Perkins et al., 1983; Siegmann et al., 1985; Lee and Siegmann, 1986). Such 3D extensions are often coupled with wide-angle modifications to achieve maximum utility from the models (Botseas et al., 1983; Thomson and Chapman, 1983).

### 4.7.4 Range-refraction corrections

Range-refraction corrections are introduced to accommodate propagation through strong oceanic fronts. The standard parabolic approximation is known to have intrinsic phase errors that will degrade the accuracy of any PE solution for long-range propagation in the ocean (Tappert and Lee, 1984; Jensen and Martinelli, 1985). The accuracy can be improved using updated mean phase speeds as a function of range. Tolstoy et al. (1985) suggested a simple transformation of the sound-speed profile to compensate for the inability of the PE approach to correctly locate the range of the signal turning points (e.g. Brock etal., 1977). Schurman etal. (1991) developed an energy-conserving PE model that incorporated range refraction. A model called LOGPE (Berman et al., 1989) used a logarithmic expression for the index of refraction in the standard PE in order to closely approximate solutions of the more exact Helmholtz equation in weakly range-dependent ocean environments.

### 4.7.5 High-frequency adaptations

At frequencies higher than about 500 Hz, PE models become impractical due to excessive execution times. Computational intensity is proportional to the number of range-interval steps. As the frequency increases, the step size decreases and more steps are thus required to achieve the desired prediction range. High-frequency approximations can be obtained by introducing a hybrid approach that combines aspects of parabolic approximations and ray theory (Tappert et al., 1984). A special-purpose microcomputer model (PESOGEN) was developed to compute the acoustic field in a fully range-dependent environment over a wide range of frequencies for both deep-water and shallow-water regions with greatly reduced execution times (Nghiem-Phu et al., 1984).

The so-called calculation frequency method (CFM) can be used to extend PE solutions to high frequencies with execution speeds typically associated with those at low frequencies (Moore-Head etal., 1989). The CFM produces range-averaged TL information by substituting a (low) calculation frequency for a (high) prediction frequency. This substitution is accomplished by sacrificing volume attenuation and boundary-loss phase information under the assumption that diffraction and interference effects are not important at the desired (high) prediction frequency.

### 4.7.6 Time-domain applications

Only frequency-domain applications of the PE have been explored thus far. However, the parabolic approximation has been adapted to consider timedomain applications. Specifically, the acoustic pressure is advanced in time by the so-called progressive wave equation (PWE). However, the PWE is, for practical purposes, limited to narrow-angle propagation paths. Collins (1988b) derived an inverse Fourier transform of the wide-angle parabolic equation referred to as the time-domain PE (TDPE). The TDPE advances the acoustic pressure field in range. Collins (1988b) used the TDPE to investigate the effects of sediment dispersion on pulse (or broadband) propagation in the ocean. Orchard et al. (1992) described the development of a TDPE model appropriate for use in 3D ocean acoustic propagation problems.

The progressive wave equation has also been adapted to model the propagation of a non-linear acoustic pulse that is subject to refraction and diffraction. McDonald and Kuperman (1987) examined the problem of acoustic propagation of finite-amplitude pulses and weak shocks in the ocean where refraction can lead to caustic formation. Their result was a first-order non-linear PWE (NPE), which is the non-linear time-domain counterpart of the linear frequency-domain PE. When the non-linear term is omitted, the NPE reduces to the linear frequency-domain PE. The simplicity of their formulation suggests additional applications to broadband linear acoustic problems in the ocean.

### 4.8 The RAYMODE model - a specific example

The RAYMODE model, originally developed by Leibiger (1968), numerically solves the wave equation for underwater acoustic propagation using the multipath expansion approach (Section 4.5). RAYMODE is actually a hybrid approach involving elements of both ray theory and wave theory. RAYMODE is a stand-alone propagation model intended for application to passive sonar performance prediction problems in range-independent ocean environments.

Since its original development, the RAYMODE passive acoustic propagation model was widely used in the naval sonar modeling community. The US Navy later brought RAYMODE under configuration management where it was maintained as an interim standard for many years. While use of this model has subsequently been eclipsed by the development of more capable models, RAYMODE continues to provide an instructive example of acoustic model construction. Additional sources of information for RAYMODE are provided in Section 4.9.

The model accommodates beam patterns at both the source and the receiver. Boundary interaction losses are also included. Sea-surface losses are computed by a user-supplied table or by internal algorithms. Surface-scattering losses assume two mechanisms: (1) a high-frequency, large-roughness loss (Beckmann and Spizzichino, 1963: 89; Clay and Medwin, 1964) and (2) a low-frequency loss (Marsh and Schulkin, 1962a). Sea-floor losses are computed either through identification of bottom province types in order to access stored values or by a user-supplied table. The volume attenuation expression of Thorp (1967) is used. Transmission loss is calculated both coherently and incoherently. The basic assumptions incorporated in RAYMODE are:

1 range-independent sound speed and bathymetry;

2 sound-speed profile fit with segments such that the index-of-refraction squared is a linear function of depth;

3 multiple reflections due to sound-speed discontinuities are ignored except in surface duct situations;

4 plane wave reflection coefficients are assumed;

5 only a finite number of ray cycles is considered in the multipath evaluation of pressure integrals;

6 a harmonic (single-frequency) source is assumed; and

### 7 a water density of unity is assumed.

The specific implementation of the multipath expansion technique is explained as follows. The infinite integral in Equation (4.13) is expressed as a sum of finite integrals, each of which is associated with a particular ray path. In Equation (4.13), 0 is taken to represent the acoustic field pressure and the integration is performed piecewise over particular ray-path regions defined by the limiting rays (R.L. Deavenport, 1978: unpublished manuscript):

where £i = to/ci, ci is the sound speed at vertexes, SD the surface duct region, CZ the convergence zone region and BB the bottom-bounce region.

Figure 4.12 illustrates a typical geometry. Here, zs is the source depth, C0 is the sound speed at the sea surface, C1 is the sound speed at the sonic layer depth (SLD) and c2 is the sound speed at the sea floor (at the base of the water column). The particular sound-speed profile and source depth geometry illustrated here would support all three ray-path families described (i.e. surface duct, CZ and bottom bounce). Note that:

cos 0s c s

Figure 4.12 Typical environment for the RAYMODE propagation model showing a sound-speed profile with a source located at depth zs.

where 0s is the maximum ray angle leaving the source and cs is the sound speed at the source depth. The integral for each region is expanded into four parts:

0 0