# Broadband properties of potential and kinetic energies in an oceanic waveguide

Julien Flamant, Julien Bonnel

► To cite this version:

Julien Flamant, Julien Bonnel. Broadband properties of potential and kinetic energies in an oceanic waveguide. Journal of the Acoustical Society of America, 2023, 153 (5), pp.3012. 10.1121/10.0019545 . hal-04104816

HAL Id: hal-04104816

https://hal.science/hal-04104816

Submitted on 24 May 2023

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
# Broadband properties of potential and kinetic energies in an oceanic waveguide

Julien Flamant¹,ª ID and Julien Bonnel² ID

¹Université de Lorraine, Centre National de la Recherche Scientifique, Centre de Recherche en Automatique de Nancy, F-54000 Nancy, France

²Applied Ocean Physics & Engineering, Woods Hole Oceanographic Institution, Woods Hole, Massachusetts 02540, USA

## ABSTRACT:

The energetic properties of an acoustic field can be quantified through the potential ($E_p$) and kinetic ($E_k$) energies. This article derives broadband properties of $E_p$ and $E_k$ in an oceanic waveguide, with restriction to a far-field context under which the acoustic field can be described by a set of propagating trapped modes. Using a set of reasonable assumptions, it is analytically demonstrated that, when integrated over a wide enough frequency-band, $E_p = E_k$ everywhere in the waveguide, except at four specific depths: $z=0$ (sea surface), $z=D$ (seafloor), $z=z_s$ (source depth), and $z = D - z_s$ (mirrored source depth). Several realistic simulations are also presented to show the relevance of the analytical derivation. It is notably illustrated that, when integrated over third-octave bands, $E_p \approx E_k$ within 1 dB everywhere in the far-field waveguide, except in the first few meters of the water column (on a dB scale, no significant difference is found between $E_p$ and $E_k$ for $z=D$, $z=z_s$, and $z = D - z_s$).

© 2023 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1121/10.0019545

(Received 12 January 2023; revised 29 March 2023; accepted 5 May 2023; published online 23 May 2023)

[Editor: Oleg A. Godin]

Pages: 3012–3024

## I. INTRODUCTION

Historically, ocean acoustics relies on the measure and use of underwater sound pressure $p$, with a wide variety of applications, from marine ecology to underwater warfare. However, it is well known that the acoustic field can also be described using other quantities, such as the particle velocity $v$. The pressure $p$ and velocity $v$ are usually derived in the linear acoustic regime; as such, they represent small oscillations around their background values. In a static fluid (an accepted approximation in ocean acoustics), they are related through

$$ \frac{\partial v}{\partial t} = -\frac{1}{\rho} \nabla p, \qquad (1) $$

with $\rho$ the ocean density and $t$ the time. Further noting $c$ as the sound speed in the oceanic waveguide, the energetic properties of the associated acoustic field are usually described through the potential ($E_p$) and kinetic ($E_k$) energy densities, with

$$ E_p = \frac{1}{2} \frac{1}{\rho c^2} |p|^2 \qquad (2) $$

interpreted as the sound pressure level and

$$ E_k = \frac{1}{2} \rho ||v||^2 \qquad (3) $$

as the particle motion level. The main objective of this article is to derive and illustrate broadband properties of $E_p$ and $E_k$ in an oceanic waveguide. The scope is restricted to the context of a single sound source and propagation over ranges that exceed twice the water depth. In such a context, the acoustic field can be described as a sum of propagating trapped modes, and a fundamental result is that, when integrated over a wide enough frequency band, $E_p \approx E_k$. This result opens the door to simple broadband underwater measurements of $E_k$ through $E_p$. Note that a similar topic, restricted to the narrowband case, is discussed in Nedelec et al. (2021). In short, if the field behaves like a plane wave, then $E_k$ can be obtained through $E_p$. The present article extends the comparison to broadband measure in complex oceanic waveguides.

Further, the article provides a theoretical justification of previous empirical observations, e.g., $E_p = E_k$ (within calibration uncertainty) for ship noise measured in third octave bands (Dahl and Bonnel, 2022), for broadband explosive sounds (Dahl and Dall'Osto, 2020), or for pile driving sounds (Dahl et al., 2023). In other words, if one is interested into broadband particle motion level (PML) from an acoustic source in an ocean waveguide, one can make a pressure measurement, compute the broadband sound pressure level (SPL), and infer PML from SPL. This is important because making a pressure measurement with a hydrophone is significantly simpler than making a particle motion measurement with a vector sensor. This result is useful for practitioners that are exclusively interested in the energetic properties of the particle motion field, such as researchers

<footnote>ªElectronic mail: julien.flamant@cnrs.fr</footnote>
studying the impact of noise pollution on fish and crustaceans (Nedelec et al., 2016; Popper and Hawkins, 2018). However, it is important to remember that this result is obtained by describing the acoustic field using a finite set of propagating modes (i.e., the source/receiver range exceeds ≈ 2 water depth). An interesting counterexample is provided in Dahl and Dall'Osto (2022): if a ship transits directly over a vector sensor, then the modal assumption does not hold, and $E_p$ can be drastically different from $E_k$.

The remainder of the article is organized as follows. Section II formalizes the problem, presents the physical background and basic equations required for the article, and introduces the simulated scenarios that will be used through the article. Section III introduces a non-stationary phase approximation that drastically simplifies the formulas for $E_k$ and $E_p$. The final relationship between $E_p$ and $E_k$ obtained after frequency integration [i.e., the demonstration of Eq. (7)] is derived in Sec. IV. Numerical simulations that illustrate the relevance and generality of the proposed derivation are given in Sec. V. Last, Sec. VI concludes the article.

## II. PRELIMINARIES

### A. Problem definition

The scope of this paper is restricted to azimuthally symmetric waveguides, i.e., propagation occurs in a 2D $r-z$ plan, where $r$ is the range and $z$ the depth. The particle velocity is thus a 2D vector: $v = [v_r, v_z]$, with $v_r$ and $v_z$ the horizontal and vertical components. The pressure and particle velocity field depend on range $r$, depth $z$ and frequency $f$, and so do the potential and kinetic energy densities. Further assuming a range-independent waveguide, these are properly defined as

$$E_p(r, z, f) \triangleq \frac{1}{4\rho(z)c^2(z)} |p(r, z, f)|^2, \quad (4)$$

$$E_k(r, z, f) \triangleq \frac{1}{4} \rho(z) \left( |v_r(r, z, f)|^2 + |v_z(r, z, f)|^2 \right), \quad (5)$$

with $\rho(z)$ and $c(z)$ the density and sound speed profiles. Compared to Eqs. (2) and (3), the definitions (4) and (5) include an extra 1/2 prefactor to account that $E_p$ and $E_k$ will be computed using complex harmonic quantities in the subsequent sections.

In such a context, it is known that for any range $r$ and frequency $f$,

$$\int_{0}^{\infty} E_p(r, z, f) dz = \int_{0}^{\infty} E_k(r, z, f) dz. \quad (6)$$

A succinct demonstration of this result can be found in Pierce (1985). In short, the key idea to derive this property is to expand the acoustic field into normal modes (see Sec. II B), and to use the mode orthonormality to simplify the expression. A full derivation of Eq. (6) is provided in Appendix A for completeness. Although Eq. (6) is very elegant, it has no simple practical application. Since integration

is over depth, using this property to infer $E_k$ from $E_p$ requires the use of a dense vertical line array (VLA) spanning the whole water column (and, in theory, the upper part of the seafloor). VLAs are traditional assets for acoustical oceanographers, but their deployment is cumbersome, and their price prohibitive for many.

In this study, we instead investigate the relationship between $E_p$ and $E_k$ when integrated along the frequency axis. As stated in the introduction, empirical observations by Dahl and Dall'Osto (2020) and Dahl and Bonnel (2022) suggest that, at a given $(r, z)$ location,

$$\int_{f_1}^{f_2} E_p(r, z, f) df \approx \int_{f_1}^{f_2} E_k(r, z, f) df, \quad (7)$$

where $f_1$ and $f_2$ define the frequency band on which the averaging is performed. This approximation is easily verified using numerical simulations (a few will be presented in Sec. V). This article studies the relation (7) from a theoretical viewpoint. It notably provides a good understanding of its validity conditions (without resorting to the numerical evaluation of $E_k$ and $E_p$).

### B. Normal modes in shallow water acoustic propagation

The theory of normal modes enables a convenient modelling of low-frequency acoustic propagation in shallow water. Consider a broadband point source signal $\Omega(f)$ emitted at depth $z_v$, the pressure field $p$ received at distance $r \gg 1$, depth $z$ and frequency $f$ reads [Jensen et al. (2011), Chap. 5]

$$p(r, z, f) = \Omega(f) \frac{e^{i\pi/4}}{\rho(z_s)\sqrt{8\pi r}} \times \sum_{m=1}^{M} \Psi_m(z_s, f) \Psi_m(z, f) \frac{e^{ik_{rm}(f)r}}{\sqrt{k_{rm}(f)}}, \quad (8)$$

where $M$ is the number of propagating modes, $\Psi_m$ and $k_{rm}$ are, respectively, the modal depth function and horizontal wavenumber associated with mode $m$. For theoretical derivations, it is also assumed that there is no attenuation in the waveguide, i.e., the wavenumbers are real numbers. However, numerical simulations that account for seafloor attenuation will be presented later in the article (see Sec. V). To shorten notations, the modal amplitude

$$A_m(z_s, f) \triangleq \Omega(f) \frac{e^{i\pi/4}}{\rho(z_s)\sqrt{8\pi}} \Psi_m(z_s, f) \quad (9)$$

is introduced. As a result, Eq. (8) becomes

$$p(r, z, f) = r^{-1/2} \sum_{m=1}^{M} A_m(z_s, f) \Psi_m(z, f) \frac{e^{ik_{rm}(f)r}}{\sqrt{k_{rm}(f)}}. \quad (10)$$

Plugging Eq. (10) into Eq. (1) gives the radial and vertical components of the particle velocity,
$$v_r(r,z,f) \approx r^{-1/2} \sum_{m=1}^{M} A_m(z_s,f) \Psi_m(z,f) \sqrt{k_{rm}(f)} \frac{e^{ik_{rm}(f)r}}{2\pi f \rho(z)}, \quad (11)$$

$$v_z(r,z,f) = r^{-1/2} \sum_{m=1}^{M} A_m(z_s,f) \Psi'_m(z,f) \frac{1}{\sqrt{k_{rm}(f)2\pi f \rho(z)}}, \quad (12)$$

where for simplicity $\Psi'_m = \partial\Psi_m/\partial z$. Note that the approximation in Eq. (11) arises from neglecting higher order terms that decrease faster than $r^{-1}$, which is justified in our normal mode context.

### C. Potential and kinetic energies

We now fix $r$ and drop the range dependence in the sequel. Combining Eqs. (10) and (4), gives the expression of the potential energy in terms of normal modes,

$$E_p(z,f) = \frac{1}{4r\rho(z)c^2(z)} \left( \sum_{m=1}^{M} \frac{|A_m(z_s,f)|^2}{k_{rm}(f)} \Psi_m^2(z,f) \right. \\ \left. + \sum_{\substack{m,n=1 \\ m \neq n}}^{M} \frac{A_m(z_s,f)A_n^*(z_s,f)}{\sqrt{k_{rm}(f)k_{rn}(f)}} \times \Psi_m(z,f)\Psi_n(z,f) e^{i\Delta k_{mn}(f)r} \right), \quad (13)$$

where $\Delta k_{mn}(f) \triangleq k_{rm}(f) - k_{rn}(f)$. Similarly, the kinetic energy is obtained by combining the radial (11) and vertical (12) particle velocity components with Eq. (5), leading to

$$E_k(z,f) = \frac{1}{16r\pi^2 f^2 \rho(z)} \left( \sum_{m=1}^{M} |A_m(z_s,f)|^2 \right. \\ \left. \times \left[ \Psi_m^2(z,f) k_{rm}(f) + \frac{\Psi_m'^2(z,f)}{k_{rm}(f)} \right] \right. \\ \left. + \sum_{\substack{m,n=1 \\ m \neq n}}^{M} A_m(z_s,f) A_n^*(z_s,f) e^{i\Delta k_{mn}(f)r} \right. \\ \left. \times \sqrt{k_{rm}(f)k_{rn}(f)} \Psi_m(z,f) \Psi_n(z,f) \right. \\ \left. + \sum_{\substack{m,n=1 \\ m \neq n}}^{M} A_m(z_s,f) A_n^*(z_s,f) e^{i\Delta k_{mn}(f)r} \right. \\ \left. \times \frac{\Psi_m'(z,f) \Psi_n'(z,f)}{\sqrt{k_{rm}(f)k_{rn}(f)}} \right). \quad (14)$$

Equations (13) and (14) are valid for any shallow water acoustic propagation scenario conveniently described by Eqs. (10)–(12). Moreover, Eqs. (13) and (14) show that energies are comprised of two terms: one gathering the sum of energies of each mode, and a second term that encodes the sum of energies due to interferences between a mode pair $(m, n)$. This second term complicates the analysis of the properties of $E_p(z,f)$ and $E_k(z,f)$. Fortunately, we will show

in Sec. III that integrating over a sufficiently large frequency band permits to get rid of this cumbersome term, therefore enabling the analysis of the properties of frequency-integrated potential and kinetic energies.

### D. Benchmark scenarios

The aim of this study is to provide closed-form expressions to relate $E_p$ and $E_k$. Still, the article will rely on numerical simulations to illustrate the mathematical derivations. In practice, two scenarios will be considered: an ideal isovelocity waveguide, and a realistic range-independent waveguide with a non-constant sound-speed profile in the water column and a layered seafloor.

The first ideal waveguide consists of an isovelocity (sound speed $c(z) = 1500 \text{ m s}^{-1}$) isodensity ($\rho(z) = 1000 \text{ kg m}^{-3}$) water column (depth $D = 100 \text{ m}$) with an upper pressure release sea-surface and a bottom rigid seafloor. In this simple waveguide, the modal quantities ($\Psi_m$ and $k_{rm}$) have closed-form expression. This environment will be used as a proxy to establish some theoretical properties and to obtain simple physical insights.

The second layered waveguide also has water depth $D = 100 \text{ m}$ with an upper pressure release sea-surface. In the water column, the sound speed profile has respective values of 1525, 1525, 1485 and 1495 m/s for depths 0, 10, 30 and 100 m (the square of the index of refraction is linearly interpolated between those points). The seafloor consists of a sediment layer (thickness $h = 10 \text{ m}$, constant sound speed $c_1 = 1650 \text{ m/s}$, density $\rho_1 = 1600 \text{ kg m}^{-3}$, attenuation $\alpha_1 = 0.1 \text{ dB/}\lambda$) over a semi-infinite basement (constant sound speed $c_b = 2000 \text{ m/s}$, density $\rho_b = 2000 \text{ kg m}^{-3}$, attenuation $\alpha_b = 0.25 \text{ dB/}\lambda$). This realistic waveguide will be used to confirm the insights obtained with the isovelocity approximation. Note that the waveguide parameters have been arbitrarily chosen. They can be seen as a nominal model for the New England shelf, as characterized during the Shallow Water 2006 (SW06) experiment (Ballard et al., 2010; Bonnel and Chapman, 2011; Tang et al., 2007). Numerical simulations in this SW06 waveguide are run using the normal mode code ORCA (Westwood et al., 1996).

Most numerical simulations will be run for frequencies below 500 Hz. Various frequency bands will be considered to assess the influence of bandwidth on the results. In Sec. V, third octave bands centered on 63 and 125 Hz will be considered. Those have been chosen because they are the ones monitored in the European Union to assess noise pollution (Erbe, 2013). All the simulations assume flat source spectrum over the frequency bands of interest.

## III. RELATION BETWEEN POTENTIAL AND KINETIC ENERGIES: A NONSTATIONARY PHASE APPROXIMATION

### A. Frequency integration of potential and kinetic energies

Let us introduce the potential and kinetic energies integrated over the frequency band $[f_1, f_2]$ with $f_1 < f_2$,
$$E_p^{[f_1 f_2]}(z) \triangleq \int_{f_1}^{f_2} E_p(z, f) df, \qquad (15)$$

$$E_k^{[f_1 f_2]}(z) \triangleq \int_{f_1}^{f_2} E_k(z, f) df, \qquad (16)$$

where $E_p(z, f)$ and $E_k(z, f)$ are given by Eqs. (13) and (14), respectively. Note that the integrals are defined on $[f_1, f_2]$ even if the integrand only exists on a smaller interval, e.g., if the source spectrum is contained in a smaller bandwidth. This mathematical trick will be useful later on.

Separating the contributions of the sum of energies of each mode and the sum of energies of interferences, the frequency-integrated pressure can be decomposed as

$$E_p^{[f_1 f_2]}(z) = \frac{1}{4r\rho(z)c^2(z)} \left( \sum_{m=1}^{M} \mathcal{I}_{p,m}^{[f_1 f_2]}(z) + \sum_{\substack{m, n = 1 \\ m \neq n}}^{M} \mathcal{J}_{p,mn}^{[f_1 f_2]}(z) \right), \qquad (17)$$

where

$$\mathcal{I}_{p,m}^{[f_1 f_2]}(z) = \int_{f_1}^{f_2} \frac{|A_m(z_s, f)|^2}{k_{rm}(f)} \Psi_m^2(z, f) df, \qquad (18)$$

$$\mathcal{J}_{p,mn}^{[f_1 f_2]}(z) = \int_{f_1}^{f_2} \frac{A_m(z_s, f) A_n^*(z_s, f)}{\sqrt{k_{rm}(f) k_{rn}(f)}} \times \Psi_m(z, f) \Psi_n(z, f) e^{i \Delta k_{mn}(f) r} df. \qquad (19)$$

Remark that the cross-terms satisfy the symmetry $\mathcal{J}_{p,mn}^{[f_1 f_2]}(z) = (\mathcal{J}_{p,mn}^{[f_1 f_2]}(z))^*$, which guarantees that the integrated potential energy is real-valued. A similar decomposition can be obtained for the kinetic energy integrated over the frequency band $[f_1, f_2]$, i.e.,

$$E_k^{[f_1 f_2]}(z) = \frac{1}{16r\pi^2 f^2 \rho(z)} \left( \sum_{m=1}^{M} \mathcal{I}_{k,m}^{[f_1 f_2]}(z) + \sum_{\substack{m, n = 1 \\ m \neq n}}^{M} \mathcal{J}_{k,mn}^{[f_1 f_2]}(z) \right), \qquad (20)$$

with definitions and properties of $\mathcal{I}_{k,m}^{[f_1 f_2]}(z)$ and $\mathcal{J}_{k,mn}^{[f_1 f_2]}(z)$ following directly from Eq. (14). They are omitted here for convenience.

## B. Nonstationary phase approximation

Explicit expressions of integrated energies (17) and (20) are somewhat cumbersome. Therefore, it is crucial to simplify those expressions to make amenable further analysis and understanding of the properties of $E_p^{[f_1 f_2]}(z)$ and $E_k^{[f_1 f_2]}(z)$. Fortunately, the use of a nonstationary phase argument enables to neglect cross-terms $\mathcal{J}_{p,mn}^{[f_1 f_2]}(z)$ and $\mathcal{J}_{k,mn}^{[f_1 f_2]}(z)$, so that frequency integrated potential and kinetic energies are very well approximated by the sum of (integrated) energies of their individual modes. It is essential to remark

that this approximation is made possible by integrating out over a given frequency range: hence, it is not valid for $E_p(z, f)$ and $E_k(z, f)$ defined in Eqs. (13) and (14), respectively.

The approach relies on the far-field assumption ($r \gg 1$) that makes integrals $\mathcal{J}_{p,mn}^{[f_1 f_2]}(z)$ and $\mathcal{J}_{k,mn}^{[f_1 f_2]}(z)$ oscillatory integrals of the form, for $z$ and $z_s$ being fixed,

$$J(r) = \int_{f_1}^{f_2} g_{mn}(f) e^{i \Delta k_{mn}(f) r} df, \qquad (21)$$

where $g_{mn}$ can be expressed in terms of modal amplitudes, modal functions, radial wavenumbers, and source spectrum. The term oscillatory integral refers to the fact that Eq. (21) consists in the integral of a smooth, slowly varying function $g_{mn}(f)$ multiplied by a rapidly oscillating complex exponential $\exp(i\Delta k_{mn}(f)r)$. Intuitively, as $r$ increases, the complex exponential oscillates faster. As a result, the integrand $g_{mn}(f) e^{i\Delta k_{mn}(f)r}$ successively takes positive and negative values, which tend to cancel each other after integration. As a result, the overall behavior of $J(r)$ is to decrease with $r$.

Formally, the principle of non-stationary phase can be stated as follows [Stein and Murphy (1993), Chap. 8, Prop. 1]:¹ if $g_{mn}(f)$ is a smooth function compactly supported on $[f_1, f_2]$ [i.e., $g_{mn}(f_1) = g_{mn}(f_2) = 0$] and that $(\Delta k_{mn}(f))' \neq 0$ for every $f \in [f_1, f_2]$, then the integral $J(r)$ decreases faster than $r^{-\gamma}$ for any $\gamma > 1$ as $r \to \infty$. In particular, $J(r) \to 0$ as $r \to \infty$.

In practice, if one makes the following assumptions:

* the source signal $\Omega(f)$ vanishes outside the frequency band $[f_1, f_2]$, such that $\Omega(f_1) = \Omega(f_2) = 0$ (alternatively, $f_1$ may be low enough to be below the waveguide cutoff frequency, such that no modal energy is propagating);

* the group speeds of two distinct modes do not cross on the frequency band, i.e., for all $m \neq n$, $dk_m/df \neq dk_n/df$. This assumption is usually verified in relatively simple shallow water waveguides. Even if crossings do sometimes appear, e.g., around Airy phase, their associated acoustic energy will be very low²—meaning that such cases can be neglected in a first approximation;

then one can directly apply the nonstationary phase principle to integrals $\mathcal{J}_{p,mn}^{[f_1 f_2]}(z)$ and $\mathcal{J}_{k,mn}^{[f_1 f_2]}(z)$, leading to $\mathcal{J}_{p,mn}^{[f_1 f_2]}(z), \mathcal{J}_{k,mn}^{[f_1 f_2]}(z) \to 0$ as $r \to \infty$.

Therefore, frequency integrated potential and kinetic energies can be approximated for $r \gg 1$ as

$$E_p^{[f_1 f_2]}(z) \approx \frac{1}{4r\rho(z)c^2(z)} \sum_{m=1}^{M} \mathcal{I}_{p,m}^{[f_1 f_2]}(z), \qquad (22)$$

$$E_k^{[f_1 f_2]}(z) \approx \frac{1}{16r\pi^2 \rho(z)} \sum_{m=1}^{M} \mathcal{I}_{k,m}^{[f_1 f_2]}(z), \qquad (23)$$

i.e., they can be approximated by the sum of energies of individual modes. This property is essential to obtain a
closed form expression of the difference between integrated potential and kinetic energies, as we show in the Sec. III C.

### C. Numerical validation of the non-stationary phase approximation

Starting from Eqs. (22)–(23) standard algebra manipulations yield

$$E_k^{[f_1, f_2]}(z) - E_p^{[f_1, f_2]}(z) \approx \varepsilon(z) \triangleq \sum_{m=1}^{M} \varepsilon_m(z), \quad (24)$$

where $\varepsilon_m(z)$ is the energy difference per mode given by

$$\varepsilon_m(z) = \frac{1}{4r\rho(z)} \int_{f_1}^{f_2} \frac{|A_m(z_s, f)|^2 [\Psi_m^2(z, f) - \Psi_m^2(z, f)k_{zm}^2(z, f)]}{4\pi^2 f^2 k_{rm}(f)} df \quad (25)$$

with $k_{zm}(f)$ the vertical wavenumber, such that $(2\pi f)^2/c^2(z) = k_{rm}^2(f) + k_{zm}^2(z, f)$ for $m = 1, 2, \dots, M$. The quantity $\varepsilon(z)$

precisely controls the difference between integrated potential and kinetic energies. It is a function of the depth $z$, but also of the source depth $z_s$ (although not explicitly mentioned). Before going further, we illustrate the generality of the nonstationary phase approximation on the two benchmark scenarios (see Sec. II D) by comparing direct computations of $E_k^{[f_1, f_2]}(z) - E_p^{[f_1, f_2]}(z)$ against numerical evaluations of $\varepsilon(z)$.

Figure 1 depicts the energy difference $E_k^{[f_1, f_2]}(z) - E_p^{[f_1, f_2]}(z)$ and its nonstationary phase approximation $\varepsilon(z)$ for the two benchmark scenarios (see Sec. II D) and three source depths $z_s = 10, 25, 40$ m. Note that the energy difference is given using an arbitrary linear unit (a.u.), proportional to J/m³, as obtained for a unit source. We considered the frequency band [1, 250] Hz. Integrals over this frequency band were computed using a standard trapezoidal approximation with frequency step of 0.25 Hz. Moreover, the number of contributing modes was determined as $M = 33$ for the ideal isovelocity waveguide and $M = 24$ for the realistic SW06

<figure><img src="image_2.png" /><figcaption>FIG. 1. (Color online) Comparison of the energy difference $E_k^{[f_1, f_2]}(z) - E_p^{[f_1, f_2]}(z)$ (red lines) and its (nonstationary) phase approximation $\varepsilon(z)$ (black lines) for three source depths $z_s = 10, 25, 40$ m in the two benchmark scenarios described in Sec. II D, with $f_1 = 1$ Hz and $f_2 = 250$ Hz and $r = 7.5$ km. Top row: ideal isovelocity environment. Bottom row: realistic SW06 environment for the same three source depths. In each subplot, black horizontal lines give source depth $z = z_s$ and its water column mirrored version $z = D - z_s$.</figcaption></figure>
scenario. For both benchmark scenarios, one observes a very good match between direct evaluations of the energy difference $E_k^{[f_1, f_2]}(z) - E_p^{[f_1, f_2]}(z)$ and its approximation $\varepsilon(z)$, especially around peaks located at interfaces ($z=0$ and $z=D$) and mirrored source depths ($z=z_s$ and $z=D-z_s$). The specific shape of those plots (and notably the peak positions) will be discussed later. The important point here is that the numerical simulation validates the relevance of the nonstationary phase approximation approach described in Sec. III B. Notably, it makes it possible a detailed analysis of the properties of the energy difference $E_k^{[f_1, f_2]}(z) - E_p^{[f_1, f_2]}(z)$ by studying the function $\varepsilon(z)$.

# IV. RELATION BETWEEN POTENTIAL AND KINETIC ENERGIES: THE ISOVELOCITY CASE

## A. Analytic derivation of $\varepsilon(z)$

The examination of Fig. 1 shows that, for both scenarios, $E_k^{[f_1, f_2]}(z) \approx E_p^{[f_1, f_2]}(z)$ for almost all depths, excepted around interfaces ($z=0$ and $z=D$) and, more interestingly, around source depth $z=z_s$ and its water column mirrored version $z=D-z_s$. To understand such behavior, this section provides a theoretical analysis of the properties of the function $\varepsilon(z)$ for the ideal isovelocity waveguide. Remarkably, these insights translate well to the more realistic SW06 case: see Sec. V for further details and numerical experiments.

In the isovelocity case, the sound speed is constant $c(z) = c_0$ as well as the water density $\rho(z) = \rho_0$. Moreover, modes functions are frequency independent $\Psi_m(z, f) = \Psi_m(z)$ and vertical wavenumbers solely depend on the mode index $m$, i.e., $k_{zm}(z, f) = k_{zm}$. Further, the simple boundary conditions yield explicit expressions [Jensen et al. (2011), Sec. 5.4],

$$ \Psi_m(z) = \sqrt{\frac{2\rho_0}{D}} \sin k_{zm}z, \quad k_{zm} = \frac{\pi}{D} \left(m - \frac{1}{2}\right). \quad (26) $$

Plugging Eq. (26) into the energy difference per mode (25) leads to a simple closed form formula which reads

$$ \varepsilon_m(z) = C_m^{[f_1, f_2]} k_{zm}^2 \sin^2(k_{zm}z_s) \cos(2k_{zm}z), \quad (27) $$

where $C_m^{[f_1, f_2]}$ carries integration of frequency dependent terms over the band $[f_1, f_2]$ such that

$$ C_m^{[f_1, f_2]} \triangleq \frac{1}{8\pi\rho_0 r D^2} \int_{f_1}^{f_2} \frac{|\Omega(f)|^2}{4\pi^2 f^2 k_{rm}(f)} df. \quad (28) $$

As already mentioned, the energy difference $\varepsilon(z)$ directly depends on the receiver depth $z$ and source depth $z_s$. Observe that $\varepsilon_m(z)$ is symmetric around $z = D/2$ such that $\varepsilon_m(D/2 - z) = -\varepsilon_m(D/2 + z)$. As a consequence, the total energy difference is also symmetric around $z = D/2$, i.e., $\varepsilon(D/2 - z) = -\varepsilon(D/2 + z)$. This explains why $\varepsilon(D/2) = \varepsilon_m(D/2) = 0$ for every $m$, meaning that the total energy difference is always equal to 0 in the center of the water

column. Moreover, remark that $\cos(2k_{zm}z)$ has local extrema at $z=0$ and $z=D$ for any $m$. This implies that $\varepsilon(z)$ admits local extrema at those locations—a first hint explaining the large peaks (in absolute value) observed at interfaces. Further properties of $\varepsilon(z)$, in particular, the expected peaks at $z=z_s$ and $z=D-z_s$ are not directly apparent from looking at Eq. (27). Instead, it requires to understand the role played by the frequency integration over the frequency band $[f_1, f_2]$. This is discussed in Sec. IV B below.

## B. Broadband modeling: Understanding the properties of $\varepsilon(z)$

First, let us rewrite Eq. (27) using standard trigonometric formulas,

$$ \varepsilon_m(z) = \frac{1}{4} C_m^{[f_1, f_2]} k_{zm}^2 [2 \cos(2k_{zm}z) - \cos(2k_{zm}(z+z_s)) - \cos(2k_{zm}(z-z_s))]. \quad (29) $$

Summing over all $M$ modes, we observe that the total energy difference can be written as

$$ \varepsilon(z) = 2\varepsilon_0(z) + \varepsilon_0(z - z_s) + \varepsilon_0(z + z_s), \quad (30) $$

i.e., $\varepsilon(z)$ is the sum of depth-translated versions of the same elementary function $\varepsilon_0(z)$ given by

$$ \varepsilon_0(z) = \frac{1}{4} \sum_{m=1}^{M} C_m^{[f_1, f_2]} k_{zm}^2 \cos(2k_{zm}z). \quad (31) $$

Finding an explicit expression of $\varepsilon_0(z)$ is cumbersome since computation of coefficients $C_m^{[f_1, f_2]}$ given by Eq. (28) cannot be carried without further assumptions. Recall that we have only assumed so far that the source $\Omega(f)$ was band limited on $[f_1, f_2]$, i.e., such that $\Omega(f) = 0$ outside this frequency interval. Indeed, the nonstationary phase argument of Sec. III B is valid for any choice of $f_1, f_2$ as long as the range $r \gg 1$. To properly understand the properties of $\varepsilon_0(z)$ [hence that of $\varepsilon(z)$] one shall make the following (reasonable) assumptions:

* $f_1 \le f_{c1}$, where $f_{c1}$ is the cut-off frequency of the first mode [in an ideal isovelocity waveguide, $f_{c1} = c_0 k_{z1} / (2\pi)$];

* $f_2 \gg f_1$, i.e., integration is performed over a sufficiently large frequency interval;

* the source has a flat spectrum $|\Omega(f)| = 1$ for $f \in [f_1, f_2]$.

As explained in Appendix B, these assumptions allow to represent $\varepsilon_0(z)$ as the convolution of a Dirichlet-like kernel $D_M(z)$ with a smoothing operator $K(z)$ such that

$$ \varepsilon_0(z) = (K * D_M)(z), \quad (32) $$

where $K(z)$ can be computed explicitly [it involves Bessel functions; see Eq. (B19)] and

$$ D_M(z) = \cos \left[ \frac{\pi}{D} Mz \right] \frac{\sin \left( \frac{\pi}{D} Mz \right)}{\sin \left( \frac{\pi}{D} z \right)}. \quad (33) $$
The function $\mathcal{K}(z)$ plays the same role as data tapers in spectral analysis (Percival and Walden, 1993). It does not change the main properties (e.g., mainlobe location) of $\mathcal{D}_M(z)$ but tends to smooth out sidelobes due to the finite frequency integration interval $[f_1, f_2]$. Its spread diminishes as $f_2$ increases, and asymptotically $\mathcal{K}(z) \to \delta(z)$ when $f_2 \to \infty$, i.e., $\mathcal{K}(z)$ converges to Dirac's delta function if integration is done over a wide enough bandwidth. In essence, this means that the general properties of $\varepsilon_0(z)$ can be inferred from that of $\mathcal{D}_M(z)$.

The function $\mathcal{D}_M(z)$ is similar to kernels usually encountered in spectral analysis (Percival and Walden, 1993) when considering truncated discrete time Fourier transforms. It has a global maximum in $z=0$ and a global minimum in $z=D$ with opposite values $\pm M$. Away from these (global) extrema, $\mathcal{D}_M(z)$ shows oscillations of rapidly decreasing amplitudes. Concentration around global extrema increases with $M$, such that $\mathcal{D}_M(z)$ tends to two Dirac delta functions located at $z=0$ and $z=D$ for $M \to \infty$. Since $\mathcal{D}_M(z)$ is 2D-periodic, translation by $z \to z+z_s$ reveals only one global extrema (a minima) on $[0, D]$, located at $z=D-z_s$, with the same decaying properties as $|z+z_s|$ increases. Similarly, by translating $z \to z-z_s$, we obtain a single global extrema (a maxima) located at $z=z_s$.

Having identified the key properties of $\mathcal{D}_M(z)$, one can now analyze the behavior of the integrated energy difference $\varepsilon(z)$, which is given by

$$ \varepsilon(z) = 2(\mathcal{K} * \mathcal{D}_M)(z) - (\mathcal{K} * \mathcal{D}_M)(z+z_s) - (\mathcal{K} * \mathcal{D}_M)(z-z_s). \quad (34) $$

We have the following properties:

* $\varepsilon(z)$ shows large deviations (i.e., peaks) at interfaces ($z=0$ and $z=D$) and at sources mirrored locations ($z=z_s$ and $z=D-z_s$);

* besides these peaks, $\varepsilon(z)$ tends to cancel out, with a residual ripple due to oscillations in the Dirichlet-like kernels;

* the finesse of peaks (at $z=0, z_s, D-z_s$ and $D$) increases with the number of propagating modes $M$;

* for $M$ sufficiently large, $|\varepsilon(0)| = |\varepsilon(D)| \approx 2\mathcal{K}(0)M$ and $|\varepsilon(z_s)| = |\varepsilon(D-z_s)| \approx \mathcal{K}(0)M$ so that the energy difference at interfaces is twice that at mirrored source depths.

Remark that peaks properties of $\varepsilon(z)$ are closely related to constructive interferences between mode energy differences, which become more and more significant as the number of propagating modes $M$ increases. The latter is directly linked to the choice of the upper band frequency $f_2$: the larger $f_2$ is, the larger $M$ is. In particular, in the limit $f_2 \to \infty$, $E_p^{[f_1, \infty]}(z) \approx E_k^{[f_1, \infty]}(z)$ for almost every $z$, excepted for $z=0, z_s, D-z_s, D$ where the energy difference tends to $\pm \infty$! On the other hand, if one considers a small integration bandwidth, it is important for $f_2$ to be large enough for at least 2 modes to be propagating, so that interference between modes can occur. If less than two modes are propagating, none of the discussion above applies.

### C. Numerical illustration

Figure 2 displays integrated energy difference $\varepsilon(z)$ for frequency bands of increasing size, i.e., $[0, 125]$ Hz, $[0, 250]$ Hz, and $[0, 500]$ Hz. They confirm the predictions from the theoretical analysis presented in Sec. IV B. Namely, $\varepsilon(z)$ is symmetric with respect the center of the water column $z=D/2$. It exhibits peaks at interfaces and mirrored sources locations. When the integration bandwidth increases, the peak amplitude increases and the mainlobe width decreases, and vice versa. Indeed, when the upper band limit $f_2$ increases, the number of propagating modes $M$ also increases. This reinforces the destructive and constructive interference patterns of $\varepsilon(z)$, leading to the behavior explained previously.

The peaks in Fig. 2 (as well as in Figs. 1 and 3) suggest that, on a linear scale, $E_p^{[f_1, f_2]}$ and $E_k^{[f_1, f_2]}$ are vastly different when $z=0, z_s, D-z_s, D$. However, this difference stays relatively small when frequency integration is done over a reasonable band, particularly if results are reported in terms of an energy ratio in dB. This will be illustrated in Secs. V B and V C. Also, note that the peaks in Figs. 1 and 2 (and later Fig. 3) imply that the energy difference is high at those specific depths; it does not inform about individual $E_p$ or $E_k$ values. As an example, while the peak at $z=0$ is particularly pronounced, this should not suggest that $E_k$ increases around the sea surface. Rather, the opposite happens: $E_p$ drastically drops while $E_k$ is relatively constant.

### V. SIMULATED EXAMPLE

The theoretical results presented above are now evaluated on realistic simulations. All simulations are run in the SW06 waveguide presented in Sec. II D. The objective here is to go beyond the isovelocity case that has been used to theoretically derive the results. This section shows the relevance and illustrates the generality of the approximations made in Sec. IV.

#### A. Comparison with the isovelocity waveguide

Section IV C presented a study of the energy difference $\varepsilon(z)$ when integrated in three different frequency bands: $[0, 125]$ Hz, $[0, 250]$ Hz, and $[0, 500]$ Hz. A similar study is performed for the SW06 simulated data, with results presented in Fig. 3. Overall, the behavior of $\varepsilon(z)$ is similar in the SW06 waveguide (Fig. 3) than in the ideal waveguide (Fig. 2): it tends to oscillate with small values over most of the water column, but has four strong peaks at $z=0, z=D, z=z_s$, and $z=D-z_s$. Interestingly, in the SW06 waveguide, $\varepsilon(z)$ loses the symmetric property around the center of the water column $z=D/2$. This is because the SW06 seafloor boundary condition is more complex than the rigid case: sound can now penetrate in the seafloor. This tends to complicate the behavior of $\varepsilon(z)$. However, it also decreases the amplitude of the two peaks that are in the bottom part of the water column. Other simulations, not shown here, illustrate that the amplitude of
<figure><img src="image_3.png" /><figcaption>FIG. 2. Evolution of $\epsilon(z)$ in the isovelocity waveguide case with $r = 7.5$ km for integration frequency bands of increasing size. Note that the range of values in the right panel is twice that of left and center panels for better readability. The horizontal gray lines give the source depth $z = z_s$ and its mirror location $z = D - z_s$.</figcaption></figure>

those peaks decrease when the seafloor boundary condition becomes less reflective (i.e., slower sound speed and/or smaller density).

An empirical observation of this property is shown in Dahl and Bonnel (2022): ship noise is measured one meter above the seafloor ($z_r \approx D$) by a vector sensor and a hydrophone. Potential and kinetic energies are evaluated and integrated into third octave bands. Figure 2 in Dahl and Bonnel (2022) shows that, for this dataset, $E_p^{[f_1 f_2]} \approx E_k^{[f_1 f_2]}$, and thus $\tilde{\epsilon}(z) \approx 0$. This may seem counter-intuitive: since measurement is done with $z_r \approx D$, one may expect a strong peak in $\tilde{\epsilon}(z)$ and thus $E_p^{[f_1 f_2]} \neq E_k^{[f_1 f_2]}$. However, measurements were done in the New England Mud Patch, an area characterized by a thick mud layer, i.e., a very soft sediment. As stated above, this smooths the peaks of $\epsilon(z)$ in the bottom part of the water column, and leads to the empirical observation $E_p^{[f_1 f_2]} \approx E_k^{[f_1 f_2]}$ even with $z_r \approx D$.

## B. Energy difference: Importance of frequency integration

The broadband behavior of $E_p$ and $E_k$ is now studied through the energy ratio $E_p/E_k$, a quantity that is used in practice for geocoustic inversion (Dahl and Dall'Osto, 2022; Ren and Hermand, 2012) and/or when doing bioacoustic experiments in tanks (Olivier et al., 2022; Popper and Fay, 2011).

Figure 4 shows the ratio $E_p/E_k$ as simulated in the SW06 waveguide for range from 5 to 10 km, with $z_s = 20$ m and $z_r = 50$ m. Results computed at a single frequency are shown as blue lines, for $f = 63$ Hz (left panel) and $f = 125$ Hz (right panel). In such a monochromatic case, $E_p/E_k$ drastically oscillates with range. This illustrates the importance of frequency integration to (approximately) equate $E_p$ and $E_k$. Results obtained with integration over third-octave bands are shown as orange lines, for central frequencies 63 and

<figure><img src="image_4.png" /><figcaption>FIG. 3. Evolution of $\epsilon(z)$ in the SW06 waveguide case with $r = 7.5$ km for integration frequency bands of increasing size. As in Fig. 2, the range of values in the right panel is twice that of left and center panels, and the horizontal gray lines give the source depth $z = z_s$ and its mirror location $z = D - z_s$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 4. (Color online) Energy ratio $E_p/E_k$ simulated in the SW06 waveguide for $z_s = 20$ m and $z_r = 50$ m, as computed with a single frequency (blue lines), in a third octave band (orange line), or from 0 to 250 Hz (black line). The left panel shows simulations centered around 63 Hz while the right panel shows simulations centered around 125 Hz.</figcaption></figure>

125 Hz. Although the frequency bands are relatively narrow (respectively, [56.1, 70.7] and [111.4, 140.3] Hz), the ratio $E_p^{[f_1 f_2]}/E_k^{[f_1 f_2]}$ is smooth, with values contained between 0.94 and 1.09. Last, results computed for the band [0, 250] Hz (black line, right panel) shows a ratio $E_p^{[f_1 f_2]}/E_k^{[f_1 f_2]}$ that is virtually constant and equal to 1.

This illustrates that, although frequency integration is required to obtain $E_p/E_k \simeq 1$, one can integrate into relatively small frequency bands, like the third octave bands that are used to characterize noise pollution in the European Union.

### C. Energy difference: Influence of source position and range

The study of the energy ratio $E_p^{[f_1 f_2]}/E_k^{[f_1 f_2]}$ is extended to various source/receiver configurations. As in Sec. V B, the focus is on three frequency bands: the third octave bands centered on 63 and 125 Hz, as well as the 0–250 Hz band. Results will be given in dB, with 0 dB meaning that $E_p^{[f_1 f_2]} = E_k^{[f_1 f_2]}$.

Figure 5 shows $E_p^{[f_1 f_2]}/E_k^{[f_1 f_2]}$ for various source and receiver depths covering the whole water column, and a constant range $r = 7.5$ km. Figure 6 further shows $E_p^{[f_1 f_2]}/E_k^{[f_1 f_2]}$ for receiver depths covering the whole water column, range from 5 to 10 km, and a constant source depth $z_s = 20$ m. Both Figs. 5 and 6 have the same layout: the first line shows the ratio $E_p/E_k$ in dB with a color scale saturated between –2 and 1 dB. The second line shows the same plot but with a threshold of 1 dB: yellow color means that $|E_p^{[f_1 f_2]}/E_k^{[f_1 f_2]}| > 1$ dB, while blue color means that $|E_p^{[f_1 f_2]}/E_k^{[f_1 f_2]}| \le 1$ dB. Note that the 1 dB threshold has been chosen as representative of calibration uncertainty for vector

sensors, e.g., see Dahl and Bonnel (2022). Overall, the two figures illustrate that $E_p^{[f_1 f_2]} \simeq E_k^{[f_1 f_2]}$ within 1 dB (i.e., calibration uncertainty) nearly everywhere in the water column, except for really shallow depths ($z < 5$ m).

### VI. CONCLUSION

The paper presents a thorough theoretical derivation of the broadband properties of $E_k$ and $E_p$ in a shallow-water waveguide. In such context, the acoustic field is expended on a finite set of propagating normal modes. Using several assumptions, the normal mode model enables an analytical derivation of important properties for $E_k$ and $E_p$. The main result of the paper is that, when integrated over a wide enough frequency-band $[f_1, f_2]$, $E_k^{[f_1 f_2]} = E_p^{[f_1 f_2]}$ nearly everywhere in the waveguide.

The assumptions required to derive the theoretical results are as follows. First, a non-stationary phase approximation is done to cancel the cross-terms in $E_k^{[f_1 f_2]}$ and $E_p^{[f_1 f_2]}$. This requires the following hypotheses:

(1) far-field ($r \gg 1$) (see Sec. III B),

(2) the source signal $\Omega(f)$ is flat (see Sec. III B) and vanishes outside of the frequency band of interest (see Sec. IV B),

(3) modal group speeds do not cross in the frequency band of interest (see Sec. III B).

Second, a closed-form expression for $E_k^{[f_1 f_2]} - E_p^{[f_1 f_2]}$ is obtained (see Sec. IV) by assuming that the environment is an ideal isovelocity waveguide (i.e., perfectly reflecting surface and perfectly rigid seabed). This concludes the demonstration and shows that, if the considered frequency band is wide enough, then $E_k^{[f_1 f_2]} = E_p^{[f_1 f_2]}$ everywhere in the
<figure><img src="image_6.png" /><figcaption>FIG. 5. (Color online) First line: energy ratio $E_p/E_k$ simulated in the SW06 waveguide with $r=7.5$ km. Second line: similar plots but with a threshold of 1 dB (yellow: $|E_p/E_k| > 1$ dB; blue: $|E_p/E_k| \le 1$ dB).</figcaption></figure>

waveguide, except at four specific receiver depths: $z=0$, $z=D$, $z=z_s$, and $z=D-z_s$.

The theoretical derivations are backed with numerical simulations showing the relevance of the approach. Those

simulations also illustrate the generality of the results, demonstrating that the approximations made in the theoretical derivation do not need to be fully respected. It is notably observed that, in a realistic waveguide with a non-constant

<figure><img src="image_7.png" /><figcaption>FIG. 6. (Color online) First line: energy ratio $E_p/E_k$ simulated in the SW06 waveguide with $z_s=20$ m. Second line: similar plots but with a threshold of 1 dB (yellow: $|E_p/E_k| > 1$ dB; blue: $|E_p/E_k| \le 1$ dB).</figcaption></figure>
sound speed profile in the water column and a layered sea-floor, the property $E_k^{[f_1, f_2]} \simeq E_p^{[f_1, f_2]}$ holds true within 1 dB for the entire far-field waveguide when integration is done in third-octave bands, except for depths around the air/water interface.

The derivations and simulations presented in this article demonstrate that, if one is exclusively interested in broadband energetic properties of the acoustic field in an oceanic waveguide, it is possible to infer $E_k^{[f_1, f_2]}$ from $E_p^{[f_1, f_2]}$, and thus to use a hydrophone instead of a vector sensor for at-sea measurements. However, let us remember that particle velocity is a 3D vector field. It contains much more information than just energetic levels, e.g., directionality (Rogers et al., 2021; Thode et al., 2019; Zeddies et al., 2010) or polarization (Bonnel et al., 2021; Dahl and Bonnel, 2022; Dall'Osto et al., 2012). Measuring, understanding, and using the vector properties of the particle motion field requires the use of dedicated vector sensors.

## ACKNOWLEDGMENTS

The authors are grateful to David Brie for fruitful discussions that yield to the explicit convolution representation of the integrated energy difference. The authors thank the reviewers for their comments which helped to better describe the scope of the paper. This work was supported by the Office of Naval Research (USA), the Délégation Générale de l'Armement (France), and the Agence Nationale de la Recherche (France).

## APPENDIX A: VERTICAL INTEGRATION PROPERTY OF $E_p$ AND $E_k$

Equation (6) states that potential and kinetic energies are equal when vertically integrated (along the water column and within the seafloor). To the best of our knowledge, this elegant result first appeared in Pierce (1985). However, the main topic in Pierce (1985) is not on $E_p$ nor $E_k$, and thus the derivation of the result is scarcely outlined. For completeness, we give below a complete detailed proof of this result. To simplify the presentation, the explicit frequency dependence of quantities is removed throughout.

Recall that for a shallow-water waveguide, the potential $E_p$ and kinetic energies $E_k$ are given by Eqs. (13) and (14), respectively. These are obtained by plugging the modal expressions (10)–(12) of the pressure, radial and vertical components of the particle velocity into their respective definitions (2) and (3). To show the equality (6), we will exploit the orthonormality of the modal functions along $z$, i.e.,

$$ \int_{0}^{\infty} \rho^{-1}(z) \Psi_{m}(z) \Psi_{n}(z) dz = \delta_{mn}, \quad (A1) $$

with $\delta_{mn} = 1$ if $m=n$ and $\delta_{mn} = 0$ if $m \neq n$. This property holds since $\Psi_m$ are eigenfunctions (associated with distinct eigenvalues $k_{rm}$) of the Sturm-Liouville problem

$$ \rho(z) \left[ \frac{1}{\rho(z)} \Psi'_{m}(z) \right]' + \left[ \frac{(2\pi f)^2}{c^2(z)} - k_{rm}^2 \right] \Psi_{m}(z) = 0 \quad (A2) $$

with boundary conditions $\Psi(z=0) = 0$ and $\Psi'(z \to \infty) = 0$.

Using the equations above, one can demonstrate the following below (whose relevance will come soon):

$$ \begin{aligned} & \int_{0}^{\infty} \rho^{-1}(z) \Psi'_{m}(z) \Psi'_{n}(z) dz \\ &= \underbrace{[\rho^{-1}(z) \Psi'_{m}(z) \Psi'_{n}(z)]_{0}^{\infty}}_{=0 \text{ (boundary conditions)}} \\ & \quad - \int_{0}^{\infty} \Psi_{n}(z) [\rho^{-1}(z) \Psi'_{m}(z)]' dz \end{aligned} \quad (A3) $$

$$ = \int_{0}^{\infty} \rho^{-1}(z) \left[ \frac{(2\pi f)^2}{c^2(z)} - k_{rm}^2 \right] \Psi_{m}(z) \Psi_{n}(z) dz \quad (A4) $$

$$ = (2\pi f)^2 \int_{0}^{\infty} \frac{\Psi_{m}(z) \Psi_{n}(z)}{\rho(z) c^2(z)} dz - k_{rm}^2 \delta_{mn}, \quad (A5) $$

where the second line was obtained using the expression of $[\rho^{-1}\Psi'_{m}]'$ given by Eq. (A2) and the last line was obtained using the modal orthonormality (A1).

One is now ready to integrate the potential and kinetic energies over $z$. Grouping auto- and cross-terms in a single double sum, one gets for the potential energy

$$ \int_{0}^{\infty} E_{p}(z) dz = \frac{1}{4r} \sum_{m,n=1}^{M} \frac{A_{m}(z_{s}) A_{n}^{*}(z_{s}) e^{i \Delta k_{mn} r}}{\sqrt{k_{rm} k_{rn}}} \times \int_{0}^{\infty} \frac{\Psi_{m}(z) \Psi_{n}(z)}{\rho(z) c^2(z)} dz. \quad (A6) $$

Similarly, for the kinetic energy,

$$ \begin{aligned} & \int_{0}^{\infty} E_{k}(z) dz \\ &= \frac{1}{16r\pi^2 f^2} \sum_{m,n=1}^{M} A_{m}(z_{s}) A_{n}^{*}(z_{s}) e^{i \Delta k_{mn} r} \\ & \quad \times \underbrace{\int_{0}^{\infty} \rho^{-1}(z) \left[ \sqrt{k_{rm} k_{rn}} \Psi_{m}(z) \Psi_{n}(z) + \frac{\Psi'_{m}(z) \Psi'_{n}(z)}{\sqrt{k_{rm} k_{rn}}} \right] dz}_{=A_{mn}} \end{aligned} \quad (A7) $$

The integral $A_{mn}$ can be further simplified by using Eqs. (A1) and (A5) such that

$$ \begin{aligned} A_{mn} = \sqrt{k_{rm} k_{rn}} & \underbrace{\int_{0}^{\infty} \rho^{-1}(z) \Psi_{m}(z) \Psi_{n}(z) dz}_{=\delta_{mn}} \\ & + \frac{1}{\sqrt{k_{rm} k_{rn}}} \int_{0}^{\infty} \rho^{-1}(z) \Psi'_{m}(z) \Psi'_{n}(z) dz \end{aligned} \quad (A8) $$
$$ = \sqrt{k_{rm}k_{rn}}\delta_{mn} + \frac{(2\pi f)^2}{\sqrt{k_{rm}k_{rn}}} \int_0^\infty \frac{\Psi_m(z)\Psi_n(z)}{\rho(z)c^2(z)} dz - \frac{k_{rm}^2\delta_{mn}}{\sqrt{k_{rm}k_{rn}}} \quad (A9) $$

$$ = \frac{(2\pi f)^2}{\sqrt{k_{rm}k_{rn}}} \int_0^\infty \frac{\Psi_m(z)\Psi_n(z)}{\rho(z)c^2(z)} dz. \quad (A10) $$

As a result, one gets

$$ \int_0^\infty E_k(z) dz = \frac{1}{4r} \sum_{m,n=1}^M \frac{A_m(z_s)A_n^*(z_s)e^{i\Delta k_{mn}r}}{\sqrt{k_{rm}k_{rn}}} \times \int_0^\infty \frac{\Psi_m(z)\Psi_n(z)}{\rho(z)c^2(z)} dz = \int_0^\infty E_p(z) dz, \quad (A11) $$

which concludes the proof of Eq. (6).

## APPENDIX B: INTEGRATED ENERGY DIFFERENCE COMPUTATIONS

First, we remark that the frequency domain assumptions made in Sec. IV B allow for an explicit computation of the coefficients $C_m^{[f_1, f_2]}$. Letting $Q = 1/(8\pi\rho_0 r D^2)$ and with our assumption that the source has flat spectrum $|\Omega(f)| = 1$, one has

$$ C_m^{[f_1, f_2]} = Q \int_{f_1}^{f_2} \frac{1}{4\pi^2 f^2 k_{rm}(f)} df. \quad (B1) $$

By assumption, $f_1 \le f_{c1}$, with $f_{c1}$ the cutoff frequency of mode 1 (more generally, the cutoff frequency of mode $m$ is defined as $f_{cm}$). Moreover, for $m > 1$, $f_{c1} < f_{cm}$, and mode $m$ has no energy propagating for $f < f_{cm}$. Therefore, Eq. (B1) can be recast as

$$ C_m^{[f_1, f_2]} = Q \int_{f_{cm}}^{f_2} \frac{1}{4\pi^2 f^2 k_{rm}(f)} df \quad (B2) $$

$$ = Q \int_{f_{cm}}^{f_2} \frac{1}{4\pi^2 f^2 \sqrt{\frac{4\pi^2 f^2}{c^2} - k_{zm}^2}} df \quad (B3) $$

$$ = Q \left[ \frac{\sqrt{\frac{4\pi^2 f^2}{c^2} - k_{zm}^2}}{4\pi^2 f k_{zm}^2} \right]_{f_{cm}}^{f_2} \quad (B4) $$

$$ = Q \left[ \frac{\sqrt{\frac{4\pi^2 f_2^2}{c^2} - k_{zm}^2}}{4\pi^2 f_2 k_{zm}^2} - \underbrace{\sqrt{\frac{4\pi^2 f_{cm}^2}{c^2} - k_{zm}^2}}_{=0 \text{ by definition of } f_{cm}} \right] \quad (B5) $$

$$ = \frac{Q}{2\pi c k_{zm}^2} \sqrt{1 - \frac{c^2 k_{zm}^2}{4\pi^2 f_2^2}} \quad (B6) $$

Back to the expression (31) of $\varepsilon_0(z)$, we observe that the coefficients $\kappa_m \triangleq C_m^{[f_1, f_2]} k_{zm}^2$ essentially play the role of Fourier coefficients of a finite half-integer cosine Fourier series on $[0, D]$, i.e.,

$$ \varepsilon_0(z) = \sum_{m=1}^{M} \kappa_m \cos \left[ \frac{2\pi}{D} \left( m - \frac{1}{2} \right) z \right], \quad (B7) $$

where we have used that $k_{zm} = (m - 1/2)\pi/D$ for the isovelocity waveguide. The somewhat complex explicit expression of Fourier coefficients $\kappa_m$ makes it cumbersome to compute the sum (B7) over the $M$ modes in a simple manner. However, by computing the Fourier transform of Eq. (B7) with respect to the variable $z$ [exploiting that $\varepsilon_0(z)$ can be defined for any value of $z$], we will show that $\varepsilon_0(z)$ can be expressed as a smoothed version of a Dirichlet-like kernel.

Let us compute the Fourier transform of Eq. (B7) and denote by $\hat{z}$ the dual variable of $z$,

$$ \hat{\varepsilon}_0(\hat{z}) \triangleq \int_{\mathbb{R}} \varepsilon_0(z) e^{-i2\pi z \hat{z}} dz = \int_{\mathbb{R}} \sum_{m=1}^{M} \kappa_m \cos \left[ \frac{2\pi}{D} \left( m - \frac{1}{2} \right) z \right] e^{-i2\pi z \hat{z}} dz \quad (B8) $$

$$ = \sum_{m=1}^{M} \frac{\kappa_m}{2} \left( \delta\left(\hat{z} - \frac{2m-1}{2D}\right) + \delta\left(\hat{z} + \frac{2m-1}{2D}\right) \right), \quad (B9) $$

where $\delta(\hat{z})$ are Dirac's delta distributions. Properties of these distributions makes it possible to rewrite Eq. (B9) as the product of two functions such that

$$ \hat{\varepsilon}_0(\hat{z}) = \hat{\mathcal{K}}(\hat{z}) \hat{\mathcal{D}}_M(\hat{z}), \quad (B10) $$

where

$$ \hat{\mathcal{K}}(\hat{z}) = \text{rect}_a(\hat{z}) \frac{Q}{2\pi c} \sqrt{1 - \frac{\hat{z}^2}{a^2}} \quad \text{with } a = \frac{2\pi f_2}{c}, \quad (B11) $$

$$ \hat{\mathcal{D}}_M(\hat{z}) = \sum_{m=1}^{M} \frac{1}{2} \delta\left(\hat{z} - \frac{2m-1}{2D}\right) + \frac{1}{2} \delta\left(\hat{z} + \frac{2m-1}{2D}\right), \quad (B12) $$

and $\text{rect}_a(\hat{z})$ is the rectangular function which is equal to 1 for $\hat{z} \in [-a, a]$ and zero otherwise. Computing the inverse Fourier transform of Eq. (B10), one obtains the following convolution product representation for $\varepsilon_0(z)$:

$$ \varepsilon_0(z) = (\mathcal{K} * \mathcal{D}_M)(z), \quad (B13) $$

where both functions $\mathcal{K}(z)$ and $\mathcal{D}_M(z)$ can be computed explicitly.

Start with $\mathcal{K}(z)$. Carrying out explicit Fourier transform computations one recognizes the integral representation of Bessel functions of the first kind [Abramowitz and Stegun (1972), Sec. 9.1.20]
$$K(z) = \frac{Q}{2\pi c} \int_{\mathbb{R}} \mathrm{rect}_a(\hat{z}) \sqrt{1 - \frac{\hat{z}^2}{a^2} e^{i2\pi z \hat{z}}} d\hat{z} \quad (B14)$$

$$= \frac{Q}{2\pi c} \int_{-a}^{a} \sqrt{1 - \frac{\hat{z}^2}{a^2} e^{i2\pi z \hat{z}}} d\hat{z} \quad (B15)$$

$$= \frac{Q}{2\pi c} \times 2 \int_{0}^{a} \sqrt{1 - \frac{\hat{z}^2}{a^2} \cos(2\pi z \hat{z})} d\hat{z} \quad (B16)$$

$$= \frac{Q}{2\pi c} \times 2a \int_{0}^{1} \sqrt{1 - \hat{z}^2 \cos(2\pi a z \hat{z})} d\hat{z} \quad (B17)$$

$$= \frac{Q}{2\pi c} \times 2a \frac{J_1(2\pi a z) \pi}{2\pi a z \frac{1}{2}} \quad (B18)$$

$$= \frac{Q}{2\pi c} \frac{J_1(2\pi a z)}{2z} \quad (B19)$$

with $J_1$ the first-order Bessel function of the first kind. Note that $K(z)$ is well defined for $z=0$, since $J_1(z) \sim_0 z/2$, leading to $K(0) = Qa/(4c) = Q\pi f_2/(2c^2)$.

The computation of $D_M(z)$ follows a different line. First, we remark that the inverse Fourier transform of $\hat{D}_M(\hat{z})$ given by Eq. (B12) can be computed directly as

$$D_M(z) = \sum_{m=1}^{M} \cos \left[ \frac{2\pi}{D} \left( m - \frac{1}{2} \right) z \right], \quad (B20)$$

which is in essence the expression of $\epsilon_0(z)$ where Fourier series coefficients $\kappa_m$ have been replaced by a vector of ones. The sum (B20) can be computed by using the following identity:

$$\sum_{m=1}^{M} e^{i(2\pi/D)(m-\frac{1}{2})z} = e^{i(\pi/D)Mz} \frac{\sin\left(\frac{\pi}{D}Mz\right)}{\sin\left(\frac{\pi}{D}z\right)}. \quad (B21)$$

This property is easily recovered by observing that the left hand side of Eq. (B21) is simply, up to a unimodular multiplicative constant, the sum of the first $M$ terms of a geometric sequence. Thus, taking the real part of Eq. (B21) yields

$$D_M(z) = \cos \left[ \frac{\pi}{D} Mz \right] \frac{\sin \left( \frac{\pi}{D} Mz \right)}{\sin \left( \frac{\pi}{D} z \right)}. \quad (B22)$$

<footnote n="1">¹It is perhaps more common to encounter the principle of stationary phase, which allows to approximate the integral (21), for $r \to \infty$, as the (discrete) sum of weighted evaluations of the integrand at stationary points of the phase [i.e., such that $(\Delta k_{mn}(f)' = 0)$]. For more details, see, e.g., Bleistein and Handelsman (1975).</footnote>

<footnote n="2">²This is not far from saying that the waveguide invariant $\beta \approx 1$.</footnote>

<footnote>Abramowitz, M., and Stegun, I. (1972). *Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables* (Dover, New York).</footnote>

<footnote>Ballard, M., Becker, K., and Goff, J. (2010). “Geoacoustic inversion for the New Jersey shelf: 3-D sediment model,” *IEEE J. Oceanic Eng.* **35**(1), 28–42.</footnote>

<footnote>Bleistein, N., and Handelsman, R. A. (1975). *Asymptotic Expansions of Integrals* (Ardent Media, New York).</footnote>

<footnote>Bonnel, J., and Chapman, N. (2011). “Geoacoustic inversion in a dispersive waveguide using warping operators,” *J. Acoust. Soc. Am.* **130**, EL101–EL107.</footnote>

<footnote>Bonnel, J., Flamant, J., Dall’Osto, D. R., Le Bihan, N., and Dahl, P. H. (2021). “Polarization of ocean acoustic normal modes,” *J. Acoust. Soc. Am.* **150**(3), 1897–1911.</footnote>

<footnote>Dahl, P., and Dall’Osto, D. (2020). “Vector acoustic analysis of time-separated modal arrivals from explosive sound sources during the 2017 seabed characterization experiment,” *IEEE J. Oceanic Eng.* **45**, 131–143.</footnote>

<footnote>Dahl, P. H., and Bonnel, J. (2022). “Vector acoustic and polarization properties of underwater ship noise,” *J. Acoust. Soc. Am.* **151**(6), 3818–3827.</footnote>

<footnote>Dahl, P. H., and Dall’Osto, D. R. (2022). “Potential and kinetic energy of underwater noise measured below a passing ship and response to sub-bottom layering,” *J. Acoust. Soc. Am.* **152**(6), 3648–3658.</footnote>

<footnote>Dahl, P. H., MacGillivray, A., and Racca, R. (2023). “Vector acoustic properties of underwater noise from impact pile driving measured within the water column,” *Front. Mar. Sci.* **10**, 1146095.</footnote>

<footnote>Dall’Osto, D. R., Dahl, P. H., and Woong Choi, J. (2012). “Properties of the acoustic intensity vector field in a shallow water waveguide,” *J. Acoust. Soc. Am.* **131**(3), 2023–2035.</footnote>

<footnote>Erbe, C. (2013). “International regulation of underwater noise,” *Acoust. Australia* **41**(1), 12–19.</footnote>

<footnote>Jensen, F., Kuperman, W., Porter, M., and Schmidt, H. (2011). *Computational Ocean Acoustics*, 2nd ed. (AIP, New York).</footnote>

<footnote>Nedelec, S. L., Ainslie, M. A., Andersson, M., Cheong, S., Halvorsen, M., Linné, M., Martin, B., Nöjd, A., Robinson, S., Simpson, S. D., Wang, L., and Ward, J. (2021). “Best practice guide for underwater particle motion measurement for biological applications,” technical report by the University of Exeter for the IOGP Marine Sound and Life Joint Industry Programme.</footnote>

<footnote>Nedelec, S. L., Campbell, J., Radford, A. N., Simpson, S. D., and Merchant, N. D. (2016). “Particle motion: The missing link in underwater acoustic ecology,” *Methods Ecol. Evol.* **7**(7), 836–842.</footnote>

<footnote>Olivier, F., Gigot, M., Mathias, D., Jezequel, Y., Meziane, T., L’Her, C., Chauvaud, L., and Bonnel, J. (2022). “Assessing the impacts of anthropogenic sounds on early stages of benthic invertebrates: The ‘Larvosonic system,’” *Limnol. Oceanogr.* **21**, 53–68.</footnote>

<footnote>Percival, D. B., and Walden, A. T. (1993). *Spectral Analysis for Physical Applications* (Cambridge University Press, Cambridge).</footnote>

<footnote>Pierce, A. D. (1985). “The natural reference wavenumber for parabolic approximations in ocean acoustics,” *Comput. Math. Appl.* **11**(7-8), 831–841.</footnote>

<footnote>Popper, A. N., and Fay, R. R. (2011). “Rethinking sound detection by fishes,” *Hear. Res.* **273**(1-2), 25–36.</footnote>

<footnote>Popper, A. N., and Hawkins, A. D. (2018). “The importance of particle motion to fishes and invertebrates,” *J. Acoust. Soc. Am.* **143**(1), 470–488.</footnote>

<footnote>Ren, Q., and Hermand, J.-P. (2012). “Passive geoacoustic inversion using waveguide characteristic impedance: A sensitivity study,” *Proc. Mtgs. Acoust.* **17**, 070095.</footnote>

<footnote>Rogers, P., Debusschere, E., Haan, D. D., Martin, B., and Slabbekoorn, H. (2021). “North Sea soundscapes from a fish perspective: Directional patterns in particle motion and masking potential from anthropogenic noise,” *J. Acoust. Soc. Am.* **150**(3), 2174–2188.</footnote>

<footnote>Stein, E. M., and Murphy, T. S. (1993). *Harmonic Analysis: Real-Variable Methods, Orthogonality, and Oscillatory Integrals* (Princeton University Press, Princeton), Vol. 3.</footnote>

<footnote>Tang, D., Moum, J. N., Lynch, J. F., Abbot, P., Chapman, R., Dahl, P. H., Duda, T. F., Gawarkiewicz, G., Glenn, S., Goff, J. A., Graber, H., Kemp, J., Maffei, A., Nash, J. D., and Newhall, A. (2007). “Shallow Water ’06: A joint acoustic propagation/nonlinear internal wave physics experiment,” *Oceanography* **20**(4), 156–167.</footnote>

<footnote>Thode, A. M., Sakai, T., Michalec, J., Rankin, S., Soldevilla, M. S., Martin, B., and Kim, K. H. (2019). “Displaying bioacoustic directional information from sonobuoys using ‘azigrams,’” *J. Acoust. Soc. Am.* **146**(1), 95–102.</footnote>

<footnote>Westwood, E. K., Tindle, C., and Chapman, N. (1996). “A normal mode model for acousto-elastic ocean environments,” *J. Acoust. Soc. Am.* **100**(6), 3631–3645.</footnote>

<footnote>Zeddies, D. G., Fay, R. R., Alderks, P. W., Shaub, K. S., and Sisneros, J. A. (2010). “Sound source localization by the plainfin midshipman fish, *Porichthys notatus*,” *J. Acoust. Soc. Am.* **127**(5), 3104–3113.</footnote>