Dmitry Yu. Mikhin

P. P. Shirshov Oceanology Institute of the Russian Academy of Sciences
36, Nakhimovskii prospect, Moscow 117218, Russia

Generalized Claerbout Parabolic Equation for Moving media (GCPE-M) ensures reciprocity and energy conservation of its solutions. The developed numerical scheme preserves these important properties in the finite-difference solution of GCPE-M. The presented model is applied to simulate full-field inversion of ocean currents in shallow sea under the presence of uncertainties in the sound speed field and bottom topography.

Current implementations of one-way methods proved to be erroneous in many range-dependent problems [1]. Errors exist even if the back-scattered field is negligible. The problem results from improper treatment of interface conditions which insure continuity of acoustic pressure and not of acoustic energy flux [1]. Available PE models provide at best only approximate conservation of acoustic energy [2, 3].

Exact solution of wave equation should also conform the principle of acoustic reciprocity or, in moving media, the flow reversal theorem (FRT) [4]. Preserving this property is extremely important when approximate one-way solutions are applied to simulate subtle non-reciprocal effects of ocean currents on sound propagation. A simple way to check accuracy of PE models in motionless media is to verify that the transmission loss is not affected by interchange in source and receiver locations. However, usually no such tests are reported. As for moving media, all available PE models are based on narrow-angle PE's [e.g. 5, 6, 7] with limited applicability.

Recently Godin [8] proposed a Generalized Claerbout PE for Moving medium (GCPE-M). It has higher asymptotic accuracy than existing PE's for moving media, possesses energy conservation and comply with FRT. Our goal is to develop an IFD scheme for solving GCPE-M, which preserves the conservation properties of the original equation.

For monochromatic 2D acoustic field of frequency propagating in the positive direction of x axis we define a complex envelope according to , where is a representative wavenumber, is a reference sound speed, x and z are the propagation range and depth. The time dependence is assumed. In the regions free from sound sources the GCPE-M is written as:

 , (1)

where is the medium density, is the local Doppler factor, is the horizontal flow component in the source-receiver direction. The dimensionless second-order differential operator is given by

 , (2)

is local wavenumber. Together with corresponding boundary conditions (BC's) at sloping




and horizontal


interfaces GCPE-M ensure exact conservation of acoustic energy and symmetry of the solution with respect to interchange of source and receiver positions. For motionless layered medium with constant density GCPE-M coincides with the standard Claerbout PE. In a general case, the most important difference is that (1) includes range derivatives of the medium parameters.

We start derivation of our IFD scheme with introducing a new dependent function

 . (4)

In terms of U GCPE-M becomes


We look for the solution on a uniform grid with the steps . First, the equation (5) is discretized in range by the Crank-Nicolson scheme. Let us define one more dependent function

 , (6)

where are the values of U at range steps n, n+1. Note that for layered medium and in general . Applying (6), the equation (5) is transformed into an equivalent system of two differential equations

 . (7a)

Alternatively, we can re-write (7a) as

 , (7b)

which is the form convenient for computer implementation: using the first equation one determines from , and then substitutes to the second equation to find . Below the indexes in and the operators in r.h.s. are omitted. Finally, the vertical operator is discretized as

 . (8)

Now we will demonstrate that the derived finite-difference scheme (7, 8) preserves the reciprocal and energy-conserving properties of the original differential equation (1). At the same time we will also establish the discrete analogs of the BC's (3). Let us consider sound propagation in the medium with reversed flow, i.e., the medium with the same fields and the opposite direction of flow velocity . The medium is taken non-absorbing with real-valued c, u, . The equations for the wave propagating in such a medium towards negative x are obtained from (7, 8) by substituting for . Hence, and are unchanged. The functions for the reversed problem are denoted as . The equations for V, differ from (7) by i replaced with -i. We multiply by , add the product and use (7a) to obtain:

 . (9)

Assume the waveguide is bounded by ideal (pressure-release or rigid) horizontal surfaces and consists of two liquid media with continuous variation of parameters separated with a stair-case boundary. Across this interface c, u, and/or can be discontinuous. The locations of the grid nodes with respect to the boundaries are shown in Fig. 1. We added artificial points below and above the physical boundaries. The medium parameters at these horizons are found by extrapolation and the functions , are determined from the discrete BC's to be derived below. The interface between the upper and lower media is assumed to be horizontal at the step and vertical at . Summing up (9) over yield:

 , (10)

where indexes (<) and (>) refer to the lower and upper media and


So, if the BC at the pressure-release surface (bottom) is chosen as

 , (11a)

the BC at the rigid surface (bottom) as

 , (11b)

and the BC and the interface of two liquids as

 , (11c)

the first two terms in the r.h.s. of (10) disappear and the last two terms cancel each other. (The same BC's are valid for .) Obviously, the discrete BC's (11a-c) approximate the continuous BC's at soft () and hard () boundaries and at a horizontal interface of two liquids (3c). The discrete BC at stairs (e.g., at ) is simply:

 . (11d)

This discrete counterpart of (3b) guarantees conservation of across a vertical interface at the step n+1.

Fig. 1. Finite-difference grid and placement of boundaries and interfaces

We proved that l.h.s. of (10) is zero. This is a discrete formulation of the FRT for moving medium. It holds exactly for the developed IFD scheme (7, 8) under BC's (11). The same derivation is applicable for media with multiple horizontal interfaces and/or higher stairs providing the field at each interface satisfies the discrete BC's. The energy conservation in the absence of dissipation also follows from (10). In this case the complex-conjugate field obeys the same equations and BC's as the field propagating in the opposite direction in the medium with reversed flow. Hence, in (10) and the value is exactly conserved at runs of the interface. From (11d) it is also conserved across the stairs. This establishes conservation of acoustic energy in the proposed IFD.

For an absorbing medium E might attenuate with range. Assume the density is real-valued, while the sound speed is complex with . In this case and r.h.s. of (9) has an extra term (now ). Summing up over j and applying BC's at the interfaces gives

 . (12)

Non-zero attenuation at any point might cause only decrease of sound energy as it should be.

In terms of finite differences precise energy conservation (energy attenuation for absorbing media) means the IFD scheme is conservative. This results in stability of the discrete solution. Equations (7, 8) and BC's (11) are consistent with the original PE (1) and BC's (3). Together with stability this establishes convergence of the scheme to the true solution for small .

Fig. 2. Transmission loss vs. propagation range for reciprocal propagations upslope (solid) and downslope (dashed). The bottom topography is shown in the lower plot.

At each fluid-fluid interface we have two non-physical values and two equations (11c) to exclude them from (7, 8). At the lower and upper boundaries there are two more non-physical values , which are eliminated using BC's (11a-b). Given N values , the total number of unknowns in (7, 8) equals the total number of equations. Using (7b, 8) one step in range requires solving two three-diagonal systems of linear equations cf. solving one such system for the standard numerical implementation of Claerbout equation. It seems to be a moderate price for conservative scheme, exact reciprocity, and energy conservation.

Usually the initial sound field and the sought-for PE solution are expressed in terms of acoustic pressure p. Hence, to initiate our marching scheme the initial field should be converted into U using (4). To obtain the acoustic pressure at some step of PE one would have to solve (4) for . This additional effort can be avoided by using the relationship . It can be done since (i) the introduced error is local and does not accumulate with range, and (ii) the approximation has high asymptotic accuracy with respect to both the numerical step and the problem parameters controlling the GCPE-M accuracy (i.e., the ratio of representative vertical and horizontal space scales of the medium inhomogeneities and the typical relative deviation of normal mode phase velocities from ).

Finally, we verify reciprocal properties of the developed model in a wedge-shaped example (Fig. 2). The 25 Hz source is located at a range of 0 km for upslope propagation and at a range of 6 km for downslope case. Source and receiver depth is 50 m. The water sound speed is 1.5 km/s while the bottom speed is 2 km/s. The density ratio between bottom and water is 2. The medium is non-absorbing (except for 500 m strongly attenuating layer near the deep bottom boundary at 2 km) and motionless. The starting fields are simulated using a normal mode algorithm. Coincidence of the TL plots at a range of 6 km proves reciprocity of the PE solution.

GCPE-M-based models can be used to simulate input data for modeling full-wave acoustic monitoring of currents in shallow water regions with strong bottom and volume inhomogeneities. Specific examples of flow inversion via Matched Non-reciprocity Tomography [9] will be presented in the report.


This work was supported, in part, by INTAS and RFBR, project 95-1002, and CRDF, project RG1-190.


[1] M. B. Porter, F. B. Jensen and C. M. Ferla, "The problem of energy conservation in one-way models", J. Acoust. Soc. Am., 89, 1058-1067 (1991)

[2] M. D. Collins, E. K. Westwood, "A higher order energy conserving parabolic equation for range-dependent ocean depth, sound speed and density", J. Acoust. Soc. Am., 89, 1068-1075 (1991)

[3] M. D. Collins, "An energy-conserving parabolic equation for elastic media", J. Acoust. Soc. Am., 94, 975-982 (1993)

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

[5] J. S. Robertson, W. L. Siegmann and M. J. Jacobson, "Current and current shear effects in the parabolic approximation for underwater sound channels", J. Acoust. Soc. Am., 77, 1768-1780 (1985)

[6] Nghiem-Phu Lan and F. D. Tappert, "Parabolic equation modeling of the effects of ocean currents on sound transmission and reciprocity in time domain", J. Acoust. Soc. Am., 78, 642-648 (1985)

[7] O. A. Godin and A. V. Mokhov, "Small-angle parabolic approximation for sound fields in a moving fluid", Sov. Phys. Acoust., 38, 243-246 (1992)

[8] O. A. Godin, "A wide-angle, energy-conserving parabolic equation for sound in moving medium", submitted for publication in Proceedings of the Third International Conference on Theoretical & Computational Acoustics, Newark, N. J. (1997)

[9] O. A. Godin and D. Yu. Mikhin, "An opportunity for improved observation of ocean currents in the coastal zone", Proc. OCEANS '96 MTS/IEEE Conference, Fort Lauderdale, Florida, pp. 345-350 (1996)