# An immersed interface method for the solution of the standard parabolic equation in range-dependent ocean environments

Roberto Sabatini, Paul Cristini

► To cite this version:

Roberto Sabatini, Paul Cristini. An immersed interface method for the solution of the standard parabolic equation in range-dependent ocean environments. Journal of the Acoustical Society of America, 2018, 143 (4), pp.EL243. 10.1121/1.5029394 . hal-01624959

HAL Id: hal-01624959

https://hal.science/hal-01624959

Submitted on 26 Oct 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An immersed interface method for the solution of the standard parabolic equation in range-dependent ocean environments

**Roberto Sabatini**ª) and **Paul Cristini**

Aix Marseille Univ, CNRS, Centrale Marseille, LMA, Marseille, France
sabatini@lma.cnrs-mrs.fr, cristini@lma.cnrs-mrs.fr

(Dated: 26 October 2017)

**Abstract:** A novel approach for the treatment of irregular ocean bottoms within the framework of the standard parabolic equation is proposed. The present technique is based on the immersed interface method originally developed by LeVeque and Li [SIAM J. Numer. Anal. **31**(4), 1019–1044, (1994)]. It is intrinsically energy-conserving and allows to naturally handle generic range-dependent bathymetries, without requiring any additional specific numerical procedure. An illustration of its capabilities is provided by solving the well-known wedge problem.

© 2017 Acoustical Society of America.

<footnote>ª) Author to whom correspondence should be addressed.</footnote>
## 1. Introduction

Propagation models based on different types of parabolic equations have been extensively used during the last four decades in underwater acoustics. Interested readers can find an exhaustive review in the book of Jensen et al. (2011) as well as in the paper of Xu et al. (2016). Since the earliest investigations of ocean sound propagation, one of the main issues has been related to the correct treatment of the interface between the water column and the seabed. Indeed, when solving a parabolic equation by approximating the seafloor as a series of stair-steps, a fundamental problem of energy conservation arises (see Jensen et al. (2011) and references therein). Current methods which allow to properly handle range-dependent bottoms essentially include: stair-step approximations within energy-conserving parabolic models (Collins and Westwood, 1991); domain rotations (Collins, 1990); and mapping techniques (Metzler et al., 2014). In this letter, a novel approach, based on the immersed interface method (IIM) originally developed by LeVeque and Li (1994), is proposed. The parabolic equation is solved on a regular Cartesian grid, in which the bathymetry is “immersed”. Away from the bottom interface, standard centered finite difference schemes are employed to compute the derivatives along the vertical axis and the Crank-Nicolson method is used for the integration in range. Conversely, at grid points lying across the bottom, the aforementioned numerical schemes are slightly modified in order for the acoustic field to satisfy not only the governing equation but also the physical interface conditions. As a result, the proposed approach is intrinsically energy-conserving and allows to handle generic range-dependent ocean
bottoms. Therefore, it can be viewed as a generalization of the IFD method developed by Lee and McDaniel (1988).

In this letter, the new technique is developed for the standard parabolic equation and for fluid-fluid interfaces. Extensions to more wide-angle parabolic equations and to fluid-solid interfaces will be considered in future works. The letter is organized as follows: after a brief review of the standard parabolic model (Section 2), the new method is presented (Section 3); an example of application is then described (Section 4); concluding remarks are finally drawn.

## 2. The standard parabolic equation and interface conditions

In the context of this work, wave propagation is assumed to be azimuthally symmetric. A cylindrical coordinate system $Orz$, with the origin $O$ on the sea surface, is thus considered. The bottom interface $z = \xi(r)$ is supposed to be irregular and the seabed is modeled as an equivalent fluid medium. A point source is placed on the $z$-axis at depth $z_s$. A sketch of the problem is illustrated in Fig. 1. In both the water column (medium 1) and the seabed (medium 2), the acoustic field in the far field can be described, in the frequency domain, by the standard parabolic equation¹ (Jensen et al., 2011)

$$ \psi_r = \mathcal{F}, \qquad \mathcal{F} = \frac{ik_0}{2} (\varepsilon^2 - 1) \psi + \frac{i}{2k_0} \psi_{zz}, \qquad (1) $$

with $\hat{p}'(r, z) = \psi(r, z)\mathcal{H}_0^{(1)}(k_0 r)$, where $\hat{p}'$ is the temporal Fourier transform of the perturbation of pressure, $\psi$ an envelope function, $\mathcal{H}_0^1$ the zeroth-order Hankel function of the first kind, $k_0 = \omega/c_0$ a reference wavenumber computed with respect to a reference speed of sound
<figure><img src="image_2.png" /><figcaption>Fig. 1. Sketch of the problem</figcaption></figure>

c₀ and ω the angular frequency. Finally, the term ε = c₀/c represents the index of refraction, where c is the speed of sound.

At the interface z = ξ(r) between the water column and the seabed, two conditions must be satisfied: the continuity of pressure and the continuity of the normal component of the particle velocity, which can be expressed in terms of the envelope function ψ as (see also Lee and McDaniel (1988))

$$ \psi^{-} = \psi^{+}, \qquad (2a) $$

$$ \psi_z^{-} - \psi_r^{-} \xi_r + k_0 \psi^{-} \frac{\mathcal{H}_1^{(1)}(k_0 r)}{\mathcal{H}_0^{(1)}(k_0 r)} \xi_r = \frac{\rho^{-}}{\rho^{+}} \left[ \psi_z^{+} - \psi_r^{+} \xi_r + k_0 \psi^{+} \frac{\mathcal{H}_1^{(1)}(k_0 r)}{\mathcal{H}_0^{(1)}(k_0 r)} \xi_r \right], \quad (2b) $$

where ρ is the density of the medium, $\mathcal{H}_1^{(1)}$ is the first-order Hankel function of the first kind and the superscripts ± indicate the limits $\lim_{z \to \xi(r) \pm}$ for a given range r.

### 3. An IIM method for the standard parabolic equation
<figure><img src="image_3.png" /><figcaption>Fig. 2. Seabed interface immersed in the computational grid. Intersection between the interface and (a) the radial direction or (b) the vertical axis.</figcaption></figure>

A uniform mesh $r^n = n\Delta r$, $z_j = j\Delta z$, with $n = 0, 1, \dots, N_r$, $j = 0, 1, \dots, N_z$, is employed, where $\Delta r$ and $\Delta z$ are the step sizes in the radial and vertical directions respectively. In what follows, the subscript $j$ and the superscript $n$ will be used to refer to point $(r^n, z_j)$. In the present approach, the bathymetry is “immersed” in the computational domain and, as schematically illustrated in Fig. 2(a) and Fig. 2(b), might cross the grid both in the radial direction and on the vertical axis. To introduce the new technique, a node $(r^n, z_j)$ away from the interface is first considered. At this regular mesh point, the second derivative $\psi_{zz,j}^n$ is approximated through the standard second order finite difference scheme

$$ \psi_{zz,j}^n = \sum_{m=-1}^{1} b_m \psi_{j+m}^n, \quad \text{with} \quad b_{-1} = b_{+1} = \frac{1}{\Delta z^2}, \quad b_0 = -\frac{2}{\Delta z^2}, $$

and the solution $\psi_j^{n+1}$ at range $r^{n+1}$ is integrated using the Crank-Nicolson method

$$ \frac{\psi_j^{n+1} - \psi_j^n}{\Delta r} = \frac{1}{2} [\mathcal{F}_j^n + \mathcal{F}_j^{n+1}] . $$
The resulting algorithm is second order accurate both in depth and in range. Nevertheless, at nodes $(r^n, z_j)$ close to the seafloor, the aforementioned schemes cannot be employed. As described in the following two paragraphs, their coefficients are then modified in such a way that the unknown function satisfies not only the governing equation but also the jump conditions.

## 3.1 Range-marching

Let $\varrho \in [r^n, r^{n+1}]$ be the interface position on the line $z = z_j$. To integrate the solution $\psi_j$ between ranges $r^n$ and $r^{n+1}$, Li (1997) elaborated the following first-order accurate scheme

$$ \frac{\psi_j^{n+1} - \psi_j^n}{\Delta r} - Q_k^{n+1/2} = \frac{1}{2} (\mathcal{F}_j^{n+1} + \mathcal{F}_j^n), $$

where the correction term $Q_j^{n+1/2}$ is given by

$$ Q_j^{n+1/2} = -\frac{r^n + \Delta r/2 - \varrho}{\Delta r} \xi_r(\varrho) \times \begin{cases} +\psi_{z,j}^{n+1} - \psi_{z,j}^n & \xi^n \ge \xi^{n+1} \\ -\psi_{z,j}^{n+1} + \psi_{z,j}^n & \xi^n < \xi^{n+1} \end{cases}. \quad (3) $$

Depending on the interface location along the vertical axis, the first derivatives appearing in Eq. (3) are computed using a standard or a modified finite difference scheme.

## 3.2 Depth derivative

Let $\xi^n \in [z_j, z_{j+1}]$ be the interface position at range $r^n$ (cf. Fig. 2(b)). As previously mentioned, modified standard finite difference methods, which take into account the jump conditions, are needed at the irregular nodes $z_j$ and $z_{j+1}$. In what follows, the derivation of such methods shall be treated in details only for the grid point $z_j$. To begin with, schemes for the derivatives $\psi_{z,j}^n$ and $\psi_{zz,j}^n$ are sought in the forms $\psi_{z,j}^n = \sum_{m=-1}^1 a_m^{(n,j)} \psi_{j+m}^n$ and
$\psi_{zz,j}^n = \sum_{m=-1}^{1} b_m^{(n,j)} \psi_{j+m}^n$. Second, up to first order accuracy, the terms $\psi_{z,j}^n$ and $\psi_{zz,j}^n$ can be written as $\psi_{z,j}^n = \psi_z^{n-} + \mathcal{O}(\Delta z)$ and $\psi_{zz,j}^n = \psi_{zz}^{n-} + \mathcal{O}(\Delta z)$. As a consequence, determining the $a_m^{(n,j)}$ and the $b_m^{(n,j)}$ amounts to expressing $\psi_z^{n-}$ and $\psi_{zz}^{n-}$ as functions of the grid values $\psi_{j-1}, \psi_j, \psi_{j+1}$. Since $\psi_z^{n-}$ and $\psi_{zz}^{n-}$ are linked to the eight jump values $\psi^{n\pm}, \psi_z^{n\pm}, \psi_r^{n\pm}, \psi_{zz}^{n\pm}$, eight equations are required to compute the unknowns $\psi^{n\pm}, \psi_z^{n\pm}, \psi_r^{n\pm}, \psi_{zz}^{n\pm}$. Two relations are provided by the jump conditions (2a) and (2b),

$$ \psi^{n-} = \psi^{n+}, \qquad (4a) $$

$$ \psi_z^{n-} - \xi_r^n \psi_r^{n-} + k_0 \frac{\mathcal{H}_1^{(1)}(k_0 r^n)}{\mathcal{H}_0^{(1)}(k_0 r^n)} \xi_r^n \psi^{n-} = \frac{\rho^{n-}}{\rho^{n+}} \left[ \psi_z^{n+} + \xi_r^n \psi_r^{n+} + k_0 \frac{\mathcal{H}_1^{(1)}(k_0 r^n)}{\mathcal{H}_0^{(1)}(k_0 r^n)} \xi_r^n \psi^{n+} \right]. \quad (4b) $$

According to Li (1997), a supplementary expression can be obtained by deriving Eq. (2a)

with respect to $r$. Using the chain rule, it follows that

$$ \psi_r^{n-} + \xi_r^n \psi_z^{n-} = \psi_r^{n+} + \xi_r^n \psi_z^{n+}. \qquad (5) $$

Furthermore, Eq. (1) must be satisfied on both sides of the interface, i.e.

$$ \psi_r^{n-} = \mathcal{F}^{-}, \quad (6a) \qquad \psi_r^{n+} = \mathcal{F}^{+}, \quad (6b) $$

The last three equations are given by the following truncated Taylor expansions

$$ \psi_{j-1}^n = \psi^{n-} + (z_{j-1} - \xi^n) \psi_z^{n-} + \frac{1}{2} (z_{j-1} - \xi^n)^2 \psi_{zz}^{n-}, \quad (7a) $$

$$ \psi_j^n = \psi^{n-} + (z_j - \xi^n) \psi_z^{n-} + \frac{1}{2} (z_j - \xi^n)^2 \psi_{zz}^{n-}, \quad (7b) $$

$$ \psi_{j+1}^n = \psi^{n+} + (z_{j+1} - \xi^n) \psi_z^{n+} + \frac{1}{2} (z_{j+1} - \xi^n)^2 \psi_{zz}^{n+}. \quad (7c) $$

Finally, solving the system (4-5-6-7) allows to express the terms $\psi_z^{-}$ and $\psi_{zz}^{-}$ as functions of

the grid values $\psi_{j-1}, \psi_j, \psi_{j+1}$ and thus to identify the coefficients $a_m^{(n,j)}, b_m^{(n,j)}, m = -1, \dots, 1$.
88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

In a similar manner, the derivatives $\psi_{z,j+1}^n$ and $\psi_{zz,j+1}^n$ at node $z_{j+1}$ are computed as

$\psi_{z,j+1}^n = \psi_z^{n+} = \sum_{m=0}^2 a_m^{(n,j+1)} \psi_{j+m}^n$ and $\psi_{zz,j+1}^n = \psi_{zz}^{n+} = \sum_{m=0}^2 b_m^{(n,j+1)} \psi_{j+m}^n$. The coefficients $a_m^{(n,j+1)}, b_m^{(n,j+1)}, m = 0, \dots, 2$ are determined from a linear system analogous to the previous one, where Eq. (7a) is replaced by a Taylor expansion for the term $\psi_{j+2}^n$.

It is worth emphasizing that, since the interface position $\xi(r)$ depends on the range $r$,

the terms $a_m^{(n,j)}, b_m^{(n,j)}, m = -1, \dots, 1$, and $a_m^{(n,j+1)}, b_m^{(n,j+1)}, m = 0, \dots, 2$, must be computed at each step $n$.

It is also worth noting that, although the local truncation error near the bottom becomes one order lower than at regular points, the global second order accuracy of the solution remains unaffected (Li, 1997).

To conclude, as in the IFD method, the implicit finite-difference equations which are obtained at grid points away from and close to the seafloor can be recast into a tridiagonal form, allowing standard fast linear solver to be employed. In addition, it is straightforward (although tedious) to show that the present numerical algorithm reduces to the IFD method in the case of an horizontal interface located on the line $z = z_j$.

## 4. Example of application

In order to show the capabilities of the new approach, the second wedge problem proposed by Jensen and Ferla (1990) and graphically illustrated in Fig. 3 is solved. The environment consists of a homogeneous water column ($c_1 = 1500 \text{ m.s}^{-1}$, $\rho_1 = 1000 \text{ kg.m}^{-3}$), limited above by a pressure-release flat sea surface ($\psi(r, 0) = 0$) and below by a sloping seafloor. The water depth is equal to 200 m at the source position and decreases to zero at 4 km range. The
<figure><img src="image_4.png" /><figcaption>Fig. 3. Test environment</figcaption></figure>

bottom is modeled as a homogeneous fluid half-space with a sound speed of $c_2 = 1700 \text{ m.s}^{-1}$ and a density of $\rho_2 = 1500 \text{ kg.m}^{-3}$. A source of frequency equal to $f = 25 \text{ Hz}$ is placed at $z_s = 100 \text{ m}$ depth. Finally, the Gaussian starter $\psi(0, z) = \sqrt{k_0}e^{-k_0^2(z-z_s)^2/2}$ is used (Jensen et al., 2011), where the reference wavenumber $k_0$ is defined with respect to the speed of sound in the water column, $k_0 = 2\pi f/c_1$. The computational domain is truncated at $D = 350 \text{ m}$ depth by a pressure-release false bottom ($\psi(r, D) = 0$). In order to avoid spurious reflections, the PML technique developed by Lu and Zhu (2007) is employed. The absorbing layer is located below 300 m depth. Finally, for the present calculations, the grid steps $\Delta r$ and $\Delta z$ are both set equal to 1 m.

As an illustration, the envelope function $\psi(r, z)$ is displayed in Fig. 4. Acoustic energy penetrates into the bottom at short ranges, where the incidence angle of the beam on the interface is close to $\pi/2$, and around 3.5 km. Besides, the PML technique clearly proves to be effective: in the PML layer, outgoing waves are absorbed without generating spurious reflections toward the water column. The transmission losses $\text{TL}(r, z) = -20 \log_{10}(|\psi(r, z)|/\sqrt{r})$
<figure><img src="image_5.png" /><figcaption>Fig. 4. Envelope function $\psi(r, z)$.</figcaption></figure>

computed at 30 m and 150 m depth are plotted in Fig. 5(a) and Fig. 5(b) respectively, along with the curves obtained by a coordinate rotation (Collins, 1990) and using a standard stair-step approximation of the bottom. At both depth, a very good agreement with the reference solution determined with the rotated PE equation is observed, which means that, as expected, the present results are not affected by energy losses.

## 5. Conclusion

A novel approach for the correct treatment of irregular fluid-fluid interfaces in parabolic wave equation models has been presented. The proposed technique is intrinsically energy-conserving and allows to consider generic range-dependent ocean floors. It is based on the immersed interface method, which consists in modifying the numerical algorithm in such a way that the acoustic field near an interface satisfies not only the governing equation but also the jump conditions. The present approach has been derived for the standard parabolic equation and has been validated with a well-known test case. This work represents a first
<figure><img src="image_6.png" /><figcaption>Fig. 5. Transmission loss at (a) $z = 30$ m and (b) $z = 150$ m: present results (solid lines), results obtained with a stair-step approximation (dashed lines), solution computed from a rotated equation (red crosses).</figcaption></figure>

step toward the development of a new methodology for the proper handling of irregular bottoms in the context of generic wide-angle parabolic equations.

## Acknowledgments

The authors would like to thank Dr. Bruno Lombard for fruitful discussions.

## References and links

¹The partial derivative of a function $\psi$ with respect to a variable $r$ is denoted by $\psi_r$.

Collins, M. D. (1990). “The rotated parabolic equation and sloping ocean bottoms,” J. Acoust. Soc. Am. **87**(3), 1035–1037.
147 Collins, M. D., and Westwood, E. K. (1991). “A higher-order energy-conserving parabolic equation for range-dependent ocean depth, sound speed, and density,” J. Acoust. Soc. Am. **89**(3), 1068–1075.

150 Jensen, F. B., and Ferla, C. M. (1990). “Numerical solutions of range-dependent benchmark problems in ocean acoustics,” J. Acoust. Soc. Am. **87**(4), 1499–1510.

152 Jensen, F. B., Kuperman, W. A., Porter, M. B., and Schmidt, H. (2011). *Computational Ocean Acoustics* (Springer, New York).

154 Lee, D., and McDaniel, S. T. (1988). *Ocean acoustic propagation by finite difference method* (Pergamon Press).

156 LeVeque, R. J., and Li, Z. (1994). “The immersed interface method for elliptic equations with discontinuous coefficients and singular sources,” SIAM J. Numer. Anal. **31**(4), 1019–1044.

159 Li, Z. (1997). “Immersed interface methods for moving interface problems,” Numerical Algorithms **14**(4), 269–293.

161 Lu, Y. Y., and Zhu, J. (2007). “Perfectly matched layer for acoustic waveguide modeling — benchmark calculations and perturbation analysis,” CMES **22**(3), 235–247.

163 Metzler, A. M., Moran, D., Collis, J. M., Martin, P. A., and Siegmann, W. L. (2014). “A scaled mapping parabolic equation for sloping range-dependent environments,” J. Acoust. Soc. Am. **135**(3), EL172–EL178.

166 Xu, C.-X., Tang, J., Piao, S.-C., Liu, J.-Q., and Zhang, S.-Z. (2016). “Developments of parabolic equation method in the period of 2000-2016,” Chin. Phys. B **11**(2), 1–12.