PROCEEDINGS of the
23rd International Congress on Acoustics
9 to 13 September 2019 in Aachen, Germany

The Analysis of Parabolic Equation Model Solutions using SplitStep/Finite Difference Method
Mustafa ASLAN1
1

Independent Scholar

ABSTRACTParabolic equation models solved using the split-step Fourier (SSF) algorithm is commonly studied for themodeling of wave propagation in the underwater environment. The smoothing function, in the solution ofrealistic problems, is introduced in the SSF algorithm to treat the density change between the water column
and the sediment layer. The algorithm with the smoothing approach accumulates the phase error with therange. As an alternative method, by Yevick and Thomson, the hybrid split-step/finite difference method wasdeveloped to accomplish the density discontinuity issue numerically, and it significantly improved the phaseerror. The additional operator containing density dependent terms are introduced to implement the hybrid
approach in this method. The results of the SSF parabolic equation model with the smoothing approach andwith the hybrid approach are compared and discussed to see the effect of numerical parameters on the model
implementation in this research.
Keywords: Parabolic equation, Split-Step Fourier algorithm, Smoothing function, Density discontinuity,
Split-step/finite difference method.

1. INTRODUCTION
Parabolic equation (PE) approximation to the modeling of radio wave propagation in theatmosphere was introduced by Leontovich and Fock (1). It applies to various fields related to wave
propagation such as optics, laser-beam propagation, microwave propagation and seismic propagation(2). It was used for the modelling of acoustic wave propagation in underwater acoustic, and PE (3,4)
was solved by the SSF algorithm as a numerical technique. Alternatively, the techniques, such as
finite-difference (FD) and finite-element (FE), are also applied to PE to model the ocean acoustic
propagation (5,6).
PE method, introduced by Tappert, is known small-angle PE (SPE) because the model gives theincrease of error with increasing of propagation angle. For more accurate solutions at higher angles of
propagation, the wide-angle parabolic equation (WAPE) method was introduced. The solutions of
WAPE with SSF model gives more stable and efficient solutions (7,8).
The general form of PE, the first order differential equation, is obtained by introducing the acoustic
field function in Helmholtz equation. It was a simplified differential form that expressed in terms of
the operator called the square root operator. The reduced wave equation contains the index ofrefraction which has the density and its derivative terms. In numerical solutions, the terms related tothe density are approached to describe the density discontinuity caused by the vertical density
variations at interface between two half space. Various numerical models treat density discontinuity indifferential form of wave equation by different approach. The boundary discontinuity leads to phaseerror with range in SSF numerical solutions and the appropriate approaches can be treated to decrease
the numerical errors (9). In PE/SSF methods, the smoothing functions such as the hyperbolic functions
are implemented to treat the boundary discontinuity density, but solutions accumulates the phase error
with the range, so by Yevick and Thomson (10), the hybrid split-step/finite difference method wasdeveloped to treat the density discontinuity issue, and it significantly improved the phase error.
In this paper, we will examine the smoothing approach and split-step Fourier / finite-difference(SSF/FD) approach. The results of PE/SSF model with smoothing functions and with hybrid approach
will be compared to benchmark solutions for shallow water environments.

1

aslanm1976@yahoo.com

1913

2. THEORY
2.1 The Parabolic Equation (PE) Approximation for Acoustic Wave Propagation
For the acoustic wave propagation, the complex pressure is defined as a solution with timeharmonic at frequency f in following equation.
(1)
P(r,M , z, Zt ) p(r ,M , z )eiZt
This equation is substituted in the Helmholtz equation, and it leads to another equation for underwater
acoustic propagation that is given by free-field Helmholtz equation in the following form with 3-D
cylindrical coordinates,

1 w § wp · 1 w 2 p w 2 p
 2  k02 n2 (r ,M , z ) p 0,(2)
¨r ¸  2
2
r wr © wr ¹ r wM
wz
where n(r ,M , z ) c0 / c(r , M , z ) is the acoustic index of refraction, k0 2S f / c0 is reference wave
number with frequency f, c0 is the reference sound speed.
In far field propagation, the media is assumed as a cylindrical waveguide and the p ressure field is
proportional to reversed by r because of energy conserving. To obtain more general PE form, thefield function \ is defined by the acoustic pressure field along a single radial as,\ ( r, z )k0r p(r, z )eik0r(3)
The PE is obtained by the substituting the field function into the Helmholtz equation in cylindricalcoordinates,w\
ik0 Q  1 \wr
where Q is generally called as the “square-root operator,”

Q

n2 

(4)

1 w §1 w ·U ¨¸
k02 wz © U wz ¹

(5)

where n r , z c0 c r , z is the acoustic index of refraction and c0 is typically assumed as 1500 m/s.
In the studies of modeling of underwater acoustic wave propagation based on PE, the first orderdifferential equation given above is numerically solved using various methods.
As an alternative approach taking account of the influence of density in the Helmholtz equation, thefield function is defined in the following form (8),

\

U\

(6)

which is substituted in the Helmholtz equation, and a PE form is also satisfied by this approach,w\
ik0 Q  1 \(7)
wr
The second order differential equation, the Helmholtz equation, is reduced to the first orderdifferential equation that can be solved one-way marching algorithm which is given by\ (r  'r , z ) expp ª¬ik0 'r (Q  1)) º¼\ ((r,
r z)

(8)

where the operator Q contains an effective index of refraction,

1 ª 1 w 2 U 3 § 1 wU · º ¨n  2«¸ »2k0 «¬ U wz 2 2 © U wz ¹ »¼
2

n

2

(9)

In this approach, the numeric solution treats the density contrast into the effective index ofrefraction.

1914

2.2 Split-Step Fourier (SSF) Method
The SSF method based on a marching algorithm which can be represented by the following form,\ (r  'r, z) )(r )\ (r, z)(10)where )(r ) is propagator function, PE can be written in terms of Hamiltonian -like operator, H op .

w\
ik0 H op\(11)
wr
for the numerical solution, the Dyson evolution operator (11) is introduced by
r 'r
ª
º
\ (r  'r , z )  exp « ik0 ³ dr cH op (r c) »\ (r )
(12)
r
¬
¼
where the expansion is analyzed by lower-order corrections and higher-order corrections is neglected.
The propagator operator advanced the solution with small range step is)(r ) | exp ª¬ ik0 'rH op (r ) º¼

(13)

where range-averaged operator, H op ( r ) ,r 'r

H op (r )1
dr cH op (r c)
'r ³r

(14)

and H op consists of scalar and differential operators, and these operators could be implemented bymultiplication, so the basis of SSF algorithm is based on separated terms of H op and approximations
to the square-root operator. In the model (8), the azimuthal coupling term is assumed small, and that isincorporated by binomial expansion. For the terms within the square root, the WAPE approximation ofThomson and Chapman (7) is utilized.

Q

n2 

1 w2
o QWWAPE
k02 wz 2

1 P  1 H 1

(15)

1 w2
k02 wz 2(16)

H op | Top  U op  Vop ,

(17)

H

n2  1 , P

The complete approximation for H op is
where Top , U op and Vop ,ª
1 w2 º1w2
1  «1  2 2 » , U op (n  1) , Vop  2 2(18)
2k0 r wM 2
¬ k0 wz ¼
For range-independent environment, operator H op ( r ) can be chosen as H op ( r ) | H op ( r ) which canTopbe simply expressed in terms of operators Top and U op (due to azimuthal symmetry, Vop o 0 )) (r ) e§ w2 ·
 ik0 'rTop ¨ 2 ¸¨ wz ¸  ik 'rU ( r )
0
op
©
¹

e

or ) ( r )

e

 ik0

§ w2 ·'r'r
U op ( r 'r )  ik0 'rTop ¨¨ wz 2 ¸¸  ik0 U op ( r )
©
¹
2
2

e

e

(19)

In SSF algorithm, U op is calculated by a multiplication operation in z-space, but the calculation ofoperator Top is not a simple operation that the field is transformed to vertical wave number space kz.
SSF method is essentially based on the separation of influence of the operators U op and Top . In thefundamental scheme of this algorithm, the field \ (r , z ) is firstly defined over depth at some range (r),then the operator e

 ik0 'rU op ( r )is applied on the initial field, and the result field is transformed by Fouriertransformation. The transformed field is multiplied by k-space operator, e ik0 ' rTˆop k2z. The k-space field1915

is transformed back to z-space, and the solution field, \ (r  'r , z ) is obtained at the end of process.

^\ (r  'r , z ) IFFT e

 ik0 'rTˆop ( k z )x FFT ª¬e

 ik0 'rU op ( r , z )`x \ (r , z ) º¼

(20)

Boundary Conditions and Smoothing Approach in SSF Method
In PE/SSF model, the surface boundary is assumed as a perfect reflector to provide the pressure
release boundary condition, \ ( z 0) 0 , and the image ocean method is used to realize the boundary
condition. In image ocean method, the real ocean is in the positive z-direction, and the image ocean
lays over the real ocean in the negative z-direction.
\ ( z ) \ ( z )(21)
Eq. (21) shows that the acoustic field of image ocean is equal to opposite sign of the acoustic fieldof the real ocean. In bottom boundary, the discontinuity occurs between two different media due to thecontrasting sound speed and density. These discontinuities create the numerical noise in the SSF
algorithm, so a function is needed to treat the discontinuity. The operator Uop depends on the effective
index of refraction term which includes the sound speed discontinuity and the density discontinuity inbottom interface. The discontinuities may describe as follows
c z cwater ( z )  cbottom ( z )  cwater ( z ) H z  zbottom
(22)
2.3

U z

U water  Ubottom  U water H z  zbottom(23)

The terms in the effective index of refraction Eq. (9) depend on the density change, U z , it is also

wU
w2 U
and
in numericalwz
wz 2
calculation. This mixing function is introduced analytically (8) by defining the density change at thewater-bottom interface according to the density profile Eq. (23) .

needed to a generalized function to represent the differential terms

2.4 The Hybrid Split-Step Fourier / Finite Difference Approach
For density discontinuity in bottom layer interface, the solutions of SSF algorithm based on densitysmoothing approach have phase error issue to be produced with range. Instead of smoothing approach,
Yevick and Thomson (10) were developed a hybrid technique which is the relying on the solving of
density dependent terms by finite difference calculations. The implementation of this approach canalso be adapted to existing SSF codes. They introduced another operator to split the operators as
density dependent and density independent that the density dependent operator is effectively appliedby using finite difference technique with Pade approximation. In this approach, as the density
dependent operator, J is defined and the modified square root operator is considered to be of theform

Q1 P  H  J(24)

1 wU wand now density discontinuity is treated by finite difference technique applied
k02 U wz wz
on the operator J , but the sound speed discontinuity in the layer will be still used the same treatment
without any modification as before. WAPE approximation is applied on new modified
square-root-operator that gives

where J

QWAPE

1 J  1 P  1 H  2(25)Thus, for the solution of one-way marching algorithm given in Eq. (8), an additional operation isrequired due to a new operator J , which is\ (r  'r , z ) exp ª¬ik0 'r ( 1  J  1) º¼\ (r  'r , z )

(26)

The
additional
operationisrealizedbyafirst-order
Padé
expansion.
In the final step, finite-difference approximation is applied to the expression defined by Eq. (20).

1916

ª 1
º
«1  4 1  iG J » < (r  ' r, z)
¬
¼

ª 1
º
z).
«1  4 1  iG J » < (r  ' r, z)
¬
¼

(27)

If J operator is directly implemented, it causes problematic calculation process. Yevick and Thomson
introduced an alternative operator to accomplish this issue and to solve the finite difference equations
effectively.

3. RESULTS
For numerical calculations, the shallow water environment is assumed to have flat surface andbottom boundary at a depth of 300 m from the water surface. The calculation environment consist ofwater column with sound speed of c w=1500 m/sec, density of ρ w=1.0 g/cm3 and fluid-like sedimentwith sound speed of cs=1700 m/s, density of ρs=1.5 g/cm3 . The benchmark solution is required theresults obtained from normal mode model (12) for our solutions which have the environment with flat,pressure-release surfaces.
Transmission loss is modeled for signals transmitted from a point source at zero range with a single
frequency of 100 Hz at a depth of 180 m. Transmission loss traces are computed for a receiver locatedat a distance of 20 km from source and at a depth of 100 m. The discretization values of environment
for numeric calculation are described by the range step size (dr) and depth mesh size (dz) and also
defined fast Fourier transform (FFT) size (NZ). The scaling factor (rz) is a multiplier between the dr
and dz. The depth mesh size was defined to be associated with acoustic wave length which is the rangeof λ/6–λ/30. Subsequently, the range step size was defined as some integer multiplier of dz.
Figure 1,2 and 3 display the solutions of basic SSF algorithm that uses the smoothing approach.
Different values of reference sound speed c0 , scaling factor rz and depth mesh size dz were used withthe goal of producing good match with the benchmark solution. The model results match relatively
well with the benchmark solution for c 0=1490, but the solutions in still accumulate phase error with therange depending on the different values of c 0, rz and dz.

Figure 1 – Smoothing approach with rz=5, dz=1,025 m for different values of sound speed c0

Figure 2 – Smoothing approach with c0=1500 m/s, dz=1,025 m for different values of scaling factor rz

1917

Figure 3 – Smoothing approach with c0=1500 m/s, rz=5 for different values of depth mesh size dz
In hybrid SSF/FD algorithm, Joperator is implemented before the operator term containing H andafter operator containing term P and the result is given by red trace of transmission loss in Figure 4.
J  is also applied after or before from both Hand P to test the effect of operator order and the result isgiven by blue trace of transmission loss in in Figure 4.

Figure 4 – Hybrid SSF/FD approach related to operator order
The solutions treated the hybrid SSF/FD approach in the presence of the bottom interface densitychange is displayed in Figures 5, 6 and 7. The model results match well with the benchmark solution
for c 0=1500 in Figure 5. Figure 6 shows that the solutions don’t importantly change with scaling factorrz.

Figure 5 – Hybrid approach with rz=5, dz=1,025 m for different values of sound speed c0

1918

Figure 6 – Hybrid approach with c0=1500 m/s, dz=1,025 m for different values of scaling factor rz
Figure 7 shows that the transmission loss trace at receiver depth of 100m for different depth meshvalues ( dz = 0.5m, 1.025m 1.5m and 2.5m) to be clear the effect of depth mesh size on the solution andthe best solution is at depth mesh size value of dz=1.025m.

Figure 7 – Hybrid approach with c0=1500 m/s, rz=5 for different values of depth mesh size dz
It can be observed from Figure 8 that the results from SSF/FD approach gives well agreement withthe benchmark solution, but SSF approach shows some phase and amplitude errors at certain ranges.
SSF/FD approach to boundary density discontinuity treats phase error for a depth mesh size of dz =1.025 m.

Figure 8 – Density smoothing approach and Hybrid approach compared with benchmark solution

4. CONCLUSIONS
The figures have two panels that the second panel shows clearly phase errors. The panel in left1919

shows full range, and the panel in right shows expanded view of final 3 km. It was shown that SSF
algorithm gives well solutions for a few kilometers, but then it accumulates phase error with the range.
Model solutions are sensitive to reference sound speed c 0. The solutions converge to benchmark
solution for c0=1490m/s, but the solutions significantly give the phase error.
Alternatively, the bottom density discontinuity was treated hybrid approach that was introduced by
Yevick-Thomson. Thus, standard SSF algorithm was updated with a hybrid SSF/FD algorithm that
improves the solutions. We take into account the operator order for operator applying in Eq. (20).
SSF/FD algorithm was compared to the benchmark solution and it was shown that produces highquality solutions for depth mesh size value of dz=1.025 m, sound speed value of c 0 =1500 m/s and first
choice of operator order. It is also shown that both the solutions of density smoothing approach and thesolutions of hybrid SSF/FD approach are not affected by the choice of scaling factor.

ACKNOWLEDGEMENTS
The author wish to acknowledge to DEGA for the financial support of attending conference.

REFERENCES
1. Leontovich M. A., Fock V. A. Solution of the problem of propagation of electromagnetic waves along
the earth's surface by the method of parabolic equations. J Phys of the USSR 1946; 10:13-24.
2. Jensen F. B., Kuperman W. A., Porter M. B., Schmidt H. Computational Ocean Acoustics, Woodbury,
NY: AIP Press; 1993. p. 343–405.
3. Hardin, R. H., Tappert, F. D. Applications of the split-step Fourier method to the numerical solution of
nonlinear and variable coefficient wave equations. SIAM Rev. 1973. p. 423.
4. Tappert, F. Parabolic equation method in underwater acoustics. J Acoust Soc Am. 1974; 63.
5. Lee, D. Ocean Acoustic Propagation by Finite Difference Methods, New York Pergamon; 1988.
6. Collins, M.D. Applications and time-domain solutions of higher-order parabolic equations inunderwater acoustics, J Acoust Soc Am. 1989; 1101-1102.
7. Thomson, D. J., Chapman, N. R. A wide-angle split-step algorithm for the parabolic equation.
J.Acoust.Soc.Am. 1983; 1848-1854.
8. Smith, K. B. Convergence, stability, and vaiability of shallow water acoustic predictions using a
split-step fourier parabolic equation model. J. Comput. Acous. 2001; 243-285.
9. Smith, K. B., Thomson, D. J., Hursky, P. An investigation into the bottom interface treatment inparabolic equation models utilizing split-step Fourier and finite-difference algorithms. J Acoust Soc Am.
2014; 2430.
10.Yevick, D., Thomson, D. J. A hybrid split-step/finite-difference PE algorithm for variable-density media.
J Acoust Soc Am. 1997; 1328-1335.
11. Sakurai J. J. Modern Quantum Mechanics. Benjamin/Cummings, Menlo Park, CA, 1985; 72-73.
12. R. B. Evans, “A coupled mode solution for acoustic propagation in a waveguide with stepwise depth
variations of a penetrable bottom,” J Acoust Soc Am. 1983; (74) 1: 188–195.

1920