Energy-conserving and reciprocal solutions for higher-order parabolic equations

Dmitry Yu. Mikhin

Atlantic Oceanographic and Meteorological Lab / NOAA
4301 Rickenbacker Causeway Miami, FL 33149, U.S.A.
Phone: (305) 361-4540, Fax: (305) 361-4449, E-mail:

Parabolic equations

Conservation of sound energy in non-absorbing media and acoustical reciprocity in stationary fluids at rest are the fundamental properties of acoustic field [1]. Preserving these properties in parabolic approximations is essential for accurate modeling of sound propagation in shallow water environments.

An energy-conserving and reciprocal One Way Wave Equation (OWWE) for motionless fluid was proposed in [2]. OWWE allows improved description of mode coupling without adversely affecting phase accuracy [2]. For moving media a closely related Generalized Claerbout PE (GCPE) was derived [3]. Its solutions preserve the acoustic energy in the absence of dissipation and satisfy the flow reversal theorem (FRT), which is an extension of reciprocity principle to moving media [1].

OWWE of [2] is a pseudo-differential equation for 2D problem

. (1)

Here is a complex envelope of sound pressure, is density, is a reference wavenumber, is the source amplitude, and is the horizontal coordinate along the propagation path. The operators are , and , where is a dimensionless vertical operator defined below.

The GCPE for motionless medium corresponds to [1, 1] Padé approximation of operator and two-term Taylor expansion of . In [3] it was obtained independently for moving medium by multi-scale asymptotic analysis:


where is a local Doppler factor, is an in-plane flow component. We write the equation in general form with arbitrary complex constants to reveal similarity between OWWE and GCPE solutions. In [3] and were the real coefficients of Padé and Taylor approximations and . Finally, the vertical operator is

. (3)

OWWE (1) is derived for motionless medium with . However, we can solve the equation with the vertical operator (3) with arbitrary and in (1) replaced by , in agreement with GCPE (2). Compared to GCPE, such solution does not improve asymptotic accuracy with respect to the Mach number. But it does not deteriorate the accuracy either, and has better wide-angle capabilities. Hence, this generalization is useful.

As with other PE's, there are several numerical approaches for solving (1,2) (e.g., [4, 5]). The existing algorithms can be applied to the new equations. However, the numerical scheme itself must be reciprocal and energy conserving, otherwise these advantages of the original differential PE's will be lost. This is especially important for reciprocity, because oceanic currents are slow and their effects are small. Even minor errors in numerical solution might become comparable with the current-induced effects. The existing IFD and split-step Fourier algorithms were found deficient for the new problem.

Numerical solution

An IFD algorithm preserving the conservation properties of the differential GCPE (2) with real-valued coefficients was proposed in [6]. The numerical solution satisfies discrete analogs of the acoustical reciprocity principle and FRT. Below it is generalized to handle OWWE, as well GCPE with source term and complex Padé coefficients.

Let us re-write (1,2) in terms of new dependent functions

and (4)

which are the local energy fluxes [2, 3] corresponding to (1,2). In terms of the equations become

and (1a)

. (2a)

To convert OWWE into a parabolic equation the operators in (1a) are approximated by Padé series

and . (5)

It is important to have the same set of coefficients in denominators of both approximations. Using decomposition (5), equation (1a) is solved by operator splitting method [7, Sect. 2.13]. The equation for each term of the Padé approximation is identical to (2a) with . Then solution of (1a) is , where each term in the sum is taken from the solution of (2a) with corresponding .

Consider equation (2a) on a uniform range-depth grid and apply Crank-Nicolson scheme to discretize the equation in range. Then introducing a new dependent function


the differential equation (2a) transforms into an equivalent system of two finite-difference equations [6]:


. Finally, the vertical operator is discretized as

where . Assume a point source of volume velocity is located at depth grid within range cell . Then the source term is presented in discreet form via Kroneker symbols as .

To solve the equation numerically one determines from using the first equation (7) and then substitutes to the second equation (7) to find . Each step in range requires one solution of a three-diagonal system of algebraic equations and one back-substitution.

Properties of the new IFD solution

Exact reciprocity of the numerical solution is an inherent property of the IFD scheme [6], which is valid for arbitrary coefficients . Compare propagation of sound between two point transceivers located in the grid cells and . It is assumed that direction of flow is changed to the opposite in the reverse propagation. The same derivation as in [6] proves that acoustical fields in direct and reverse transmissions are related by


where the value corresponds to in reverse transmission and are amplitudes of sources. In accordance with the FRT, the fields are symmetric with respect to interchange of source and receiver positions and simultaneous reversal of the flow direction. The same analysis gives the discrete boundary conditions on stair-step interface [6]. Reciprocity of the OWWE solution follows from reciprocity on each split-step. The reciprocal value for point sources of volume velocity, i.e., the discrete analog of acoustic pressure, is .

Energy conservation and reciprocity are equivalent in case of equations with real-valued coefficients [2, 3]. However, strict energy conservation is more a problem than an asset. Both the field components corresponding to discreet spectrum (normal modes) and the components corresponding to continuous spectrum are advanced in range without losses. In the solution of wave equation the modes should propagate (probably, with some attenuation), while the continuous spectrum should vanish. Cf. analysis in [8], propagation of continuous spectrum does not lead to infinite growth of the solution, because the energy is conserved, hence, always limited. However, the model does not reproduce an important feature of the wave equation and yields erroneous results.

To eliminate the evanescent spectrum in IFD-based PE's the operator is approximated with complex Padé coefficients . The energy of the discrete solution of proved to be non-increasing for coefficients with , . Attenuation of a field component with horizontal wavenumber is proportional to , where is a measure of how wide-angle is the modal spectrum. If is chosen close to horizontal wavenumbers of propagating modes, than additional attenuation of modes due to complex coefficients is small. Field components with horizontal wavenumbers located far from , including continuous spectrum, are effectively eliminated by IFD solution.

In terms of finite differences precise energy conservation or energy attenuation means the IFD scheme is conservative [9], hence, the discrete solution is stable. System (7) and corresponding discrete boundary conditions [6] are consistent with the original PE and its boundary conditions [3]. Together with stability this establishes convergence of the scheme to the true solution for small [9].

Analysis of energy conservation for OWWE also demonstrates non-increase of acoustic energy for properly chosen coefficients. However, derivation is beyond the volume of the present paper.


Solution with large range step

The discreet solution presented above uses Crank-Nicolson finite-difference scheme in range. The range step should be small for accurate solution. A numerical technique for solving a parabolic equation with very large steps in range was developed in [5]. Following [5], we write a formal solution of (2a) in operator form:

This operator equation is solved by the method of operator splitting [7]. The exponent is approximated by Padé series as . Then on each sub-step we define


Combining the terms these equations are transformed into the form of (2a) with coefficient replaced with . The solution is reciprocal for any set of and its energy does not increase if , .

Numerical examples

Due to the format limitations of the paper, all numerical examples are given in the accompanying paper [10].


Exact reciprocity, energy conservation, wide-angle capability and large steps in range are simultaneously achieved in IFD solution of one-way wave equation. The numerical scheme is robust, efficient, and allows generalization to 3D problems. The developed algorithm provides accurate description of acoustical effects of oceanic currents in complicated shallow-water environments. The software package is free for non-commercial and non-governmental use. It is available for download at This work is partially supported by INTAS and NRC.


[1] Godin, O. A. Reciprocity and energy theorems for waves in a compressible inhomogeneous moving fluid. Wave Motion, 25, p.143-167, 1997.

[2] Godin, O. A. Reciprocity and energy conservation within the parabolic approximation. Wave motion, 29 (2), p.175-194, 1999.

[3] Godin, O. A. A wide-angle, energy-conserving parabolic equation for sound in moving medium. Proceedings of the Third International Conference on Theoretical & Computational Acoustics, Newark, NJ, September 1997, World Scientific, Singapore, 1999.

[4] Hardin, R. H., F. D. Tappert. Applications of the split-step Fourier method to the numerical solution of nonlinear and variable coefficient wave equations. SIAM rev., 15, p.423, 1973.

[5] Collins, M. D. A split-step Padé solution for parabolic equation method. J. Acoust. Soc. Am., 93 (4, Pt. 1), p.1736-1742, 1993.

[6] Mikhin, D. Yu. Modeling Acoustic Tomography of Ocean Currents in Coastal Regions via Wide-Angle Parabolic Approximation. Fourth European Conference on Underwater Acoustics (ECUA), Rome, Italy, September 1998, p.581-586, 1998.

[7] Mitchell, A. R., D. F. Griffiths. The finite difference method in partial differential equation. Chichester; New York: Wiley, 1980.

[8] Wetton, B. T. R., G. H. Brooke. One-way wave equations for seismoacoustic propagation in elastic waveguides. J. Acoust. Soc. Am., 87 (2), p.624, 1990.

[9] Samarsky, A. A., A. V. Gulin. Stability of difference schemes. Nauka, Moskva, 1973.

[10] Mikhin, D. Yu. Simulations of full-field tomography of oceanic currents in shallow areas. In the same issue.