PekerisModeSolver

Model	UnderwaterAcoustics.PekerisModeSolver
Description	Normal mode model for constant depth iso-velocity environments
Language	Julia
Advantages	Fast, differentiable (forward mode), multi-threaded
Limitations	Iso-velocity, range-independent, pressure-release surface, fluid half-space seabed, no seabed absorption (no leaky modes)
Differentiability	ForwardDiff

PekerisModeSolver(env; ngrid=0)

A fast differentiable mode propagation model that only supports iso-velocity constant depth environments.

ngrid is the number of grid points to use for modal root finding for fluid bottom environments. If ngrid is too small, the mode solver may miss some modes. If ngrid is too large, the mode solver may take a long time to converge. The default value of ngrid of 0 will use a heuristic to automatically determine the number of grid points to use.

Experimental

Work in progress. Not benchmarked. Use at your own risk!

Notes:

• Differentiability currently works only with ForwardDiff, and excludes impulse_response().

Implementation based on mathematical description in Theory below.

Example

using UnderwaterAcoustics
using Plots

env = UnderwaterEnvironment(
  bathymetry = 5000,
  soundspeed = 1500,
  density = 1000,
  seabed = FluidBoundary(2000, 2000)
)
pm = PekerisModeSolver(env)

tx = AcousticSource(0, -500, 10)
rx = AcousticReceiver(200000, -2500)
modes = arrivals(pm, tx, rx)[1:7]      # first 7 modes

7-element Vector{UnderwaterAcoustics.ModeArrival{ComplexF64, UnderwaterAcoustics.Mode{Float64}, Float64}}:
   mode 1: kᵣ = 0.041883 + 0.0im rad/m (1499.84 m/s)
   mode 2: kᵣ = 0.04187 + 0.0im rad/m (1499.36 m/s)
   mode 3: kᵣ = 0.041847 + 0.0im rad/m (1498.56 m/s)
   mode 4: kᵣ = 0.041815 + 0.0im rad/m (1497.45 m/s)
   mode 5: kᵣ = 0.041773 + 0.0im rad/m (1496.01 m/s)
   mode 6: kᵣ = 0.041723 + 0.0im rad/m (1494.24 m/s)
   mode 7: kᵣ = 0.041663 + 0.0im rad/m (1492.15 m/s)

# plot the modes to a depth of 5 km
plot(modes[1:7], 5000)

rxs = AcousticReceiverGrid2D(200000:10:220000, -2500)
x = transmission_loss(pm, tx, rxs)

plot(200:0.01:220, x; ylims=(70,110), yflip=true, legend=false,
  xlabel="Range (km)", ylabel="Transmission loss (dB)")

Theory

The mathematical derivation of modal propagation in this section is based on:

• M. B. Porter, “The KRAKEN normal mode program”, Technical report NRL/MR/5120-92-6920, Naval Research Laboratory, 1992 (pdf).

Modal propagation

The acoustic wave equation is: \[ \nabla^2p = \frac{1}{c^2}\frac{\partial^2p}{\partial t^2}. \] In a Pekeris waveguide with a pressure-release surface and a pressure-release bottom (good approximation for a long range propagation problem, as it turns out), the boundary conditions in cylindrical coordinates are: \[ \begin{align*} p(r,0) &= 0 \quad \forall \; r, \\ p(r,-D) &= 0. \end{align*} \] For a monochromatic source of frequency \(f\) Hz located at \((0,-z_s)\) boundary condition is: \[ p(0,-z_s) = -A\exp(-i\omega t), \] where \(A\) is the source amplitude, and \(\omega = 2\pi f\). Due to symmetry of the Pekeris waveguide, we have dropped the azimuthal coordinate.

The wave equation with these boundary conditions admits an approximate modal solution of the form: \[ p(r,z,t) = \frac{i}{4}\sum_m\psi_m(z_s)\psi_m(z)H_0^{(1)}(k_mr)\exp(-i\omega t), \] where \(\psi_m(z)\) is the mode of order \(m\) and \(H_0^{(1)}\) is the Hankel function of the first kind and order 0, \(k_m\) is the horizontal wavenumber associated with mode. The mode \(\psi_m\) is given by: \[ \psi_m(z) = \sqrt{\frac{2}{D}}\sin(\gamma_m z), \] where \(\gamma_m\) is the vertical wavenumber associated with mode \(m\): \[ \gamma_m = \sqrt{k_0^2 - k_m^2}, \] for wavenumber \(k_0 = \omega/c\).

\(\gamma_m\) is chosen to ensure boundary conditions are satisfied: \[ \psi_m(0) = \psi_m(D) = 0, \] i.e., \[ \gamma_m = m\frac{\pi}{D}. \] If the bottom boundary condition is chosen for a rigid bottom: \[ \begin{align*} & \left.\frac{d}{dz}\psi_m(z)\right|_{z=D} = 0, \\ & \therefore \; \gamma_m = \left(m-\frac{1}{2}\right)\frac{\pi}{D}. \end{align*} \]

Substituting \(\psi_m\) back, we have: \[ p(r,z,t) = \frac{i}{2D}\sum_m\sin(\gamma_m z_s)\sin(\gamma_m z)H_0^{(1)}(k_mr)\exp(-i\omega t), \] The transmission loss is then given by: \[ \text{TL}(r,z) = -20\log\left| \frac{p(r,z)}{p^0_{r=1}} \right|, \] where \(p^0_{r=1}\) is the pressure due to the source at 1 m from the acoustic center of the source: \[ p^0_{r=1} = \frac{\exp(ik_0r)}{4\pi r} = \frac{\exp(ik_0)}{4\pi}. \]

Acousto-elastic boundary condition

If we use an acousto-elastic boundary condition for the seabed with sound speed \(c_b\) and density \(\rho_b\), we have a boundary condition: \[ \left.\rho_b\frac{d\psi_m}{dz}\right|_{z=D} = -\rho\sqrt{k_m^2 - k_b^2} \; \psi_m(r,D). \] for \(k_b = \omega/c_b\). Substituting \(\psi_m(r,z)\): \[ \rho_b\gamma_m\cos(\gamma_m D) + \rho\sqrt{k_m^2 - k_b^2} \; \sin(\gamma_m D) = 0. \]

Group velocity

The group velocity of mode \(m\) is \(d\omega/dk_m\): \[ \begin{align*} k_m^2 &= k_0^2 - \gamma_m^2 \\ \therefore \; 2k_m\frac{dk_m}{d\omega} &= \frac{2k_0}{c} - 2\gamma_m\frac{d\gamma_m}{d\omega} \\ \therefore \; \frac{dk_m}{d\omega} &= \frac{k_0}{c k_m} - \frac{\gamma_m}{k_m}\frac{d\gamma_m}{d\omega} \end{align*} \] We can compute \(\frac{d\gamma_m}{d\omega}\): \[ \begin{align*} &\rho_b\gamma_m\cos(\gamma_m D) + \rho\sqrt{k_m^2 - k_b^2} \; \sin(\gamma_m D) = 0 \\ \therefore \quad &\rho_b\gamma_m\cos(\gamma_m D) + \rho\sqrt{k_0^2 - \gamma_m^2 - k_b^2} \; \sin(\gamma_m D) = 0 \\ \therefore \quad &\rho_b\frac{d\gamma_m}{d\omega}\cos(\gamma_m D) - \rho_b\gamma_m\sin(\gamma_m D)D\frac{d\gamma_m}{d\omega} + \rho\zeta\cos(\gamma_m D)D\frac{d\gamma_m}{d\omega} \\ &- \frac{\rho}{2\zeta}\sin(\gamma_m D)\left(\frac{2k_0}{c} - 2\gamma_m\frac{d\gamma_m}{d\omega} - \frac{2k_b}{c}\right) = 0 \\ \therefore \quad &\left[\rho_b\cos(\gamma_m D) - \rho_bD\gamma_m\sin(\gamma_m D) + \rho D\zeta\cos(\gamma_m D) - \frac{\rho\gamma_m}{\zeta}\sin(\gamma_m D)\right]\frac{d\gamma_m}{d\omega} \\ &= -\frac{\rho k_0}{c\zeta}\sin(\gamma_m D) - \frac{\rho k_b}{c\zeta} \; \sin(\gamma_m D) \\ \therefore \quad &\frac{d\gamma_m}{d\omega} = -\frac{\rho\sin(\gamma_m D)(k_0 + k_b)}{c\zeta[\rho_b\cos(\gamma_m D) - \rho_bD\gamma_m\sin(\gamma_m D) + \rho D\zeta\cos(\gamma_m D) - \frac{\rho\gamma_m}{\zeta}\sin(\gamma_m D)]} \\ &= -\frac{\rho\sin (\gamma_m D)\omega\left(c^2-c_b^2\right)}{c^2 c_b^2 \gamma_m \sin (\gamma_m D) \left(D\text{$\rho $b} \zeta+\rho \right)+\cos(\gamma_m D) \left(c^2 \left(D\rho\omega ^2-c_b^2 \text{$\rho $b} \zeta\right)-c_b^2 D\rho\omega^2\right)+c^2 c_b^2 D\rho\gamma_m^2\cos(\gamma_m D)}. \end{align*} \] where \(\zeta = \sqrt{k_0^2 - \gamma_m^2 - k_b^2}\).