Underwater propagation solvers

The closed-form propagation loss of Underwater sound propagation knows nothing of the sound-speed profile, the seabed or the surface. When refraction and boundaries decide the answer, the field has to be computed: this guide covers the three numerical solvers of the underwater module, the physics each one discretises, and how to choose between them and the closed forms. All three assume a range-independent (horizontally stratified) ocean with a pressure-release surface, take the same $c(z)$ profile as input, and follow Jensen, Kuperman, Porter & Schmidt, Computational Ocean Acoustics.

1. The three solvers at a glance

• normal_modes solves the depth-separated Sturm-Liouville eigenvalue problem by finite differences and sums the propagating modes into the propagation loss. Validated against the ideal (pressure-release) waveguide’s exact modes.
• ray_trace integrates the ray-trajectory equations (Runge-Kutta, vectorised over all rays at once) through a sound-speed profile, reflecting at the surface and bottom, and carries the travel time along each ray as a state of the same integration. Validated against the circular-arc paths of a linear gradient and the closed-form travel time along them.
• parabolic_equation marches the standard (Tappert) PE with the split-step Fourier algorithm. Validated against free-field spherical spreading; it agrees with the normal-mode propagation loss in trend.

Three panels, three different questions about the same water. Left, rays in a Munk profile whose minimum is 1500.0 m/s at 1271 m against 1548.5 m/s at the surface: a 48 m/s spread is enough to fold the ±12° fan back on itself into convergence zones, and the ±12° ray still sweeps from 143 m to 3997 m of depth. Centre, the same environment as a PE field, where the zones appear as interference rather than as lines. Right, the two solvers on a different, 200 m shallow waveguide: they track the same trend with a mean offset of 1.6 dB, but the per-range difference reaches 25 dB at the interference nulls, which is the honest statement of what “agree in trend” means.

Show the code for this figure

import matplotlib.pyplot as plt
import numpy as np
from phonometry import underwater

# A Munk deep-water sound-speed profile for the ray and PE panels.
z = np.linspace(0.0, 5000.0, 60)
eta = 2.0 * (z - 1300.0) / 1300.0
c = 1500.0 * (1.0 + 0.00737 * (eta - 1.0 + np.exp(-eta)))

fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# (a) rays through the Munk profile, out to 100 km.
underwater.ray_trace(z, c, source_depth=1000.0,
                     launch_angles_deg=np.linspace(-12.0, 12.0, 13),
                     max_range=100e3, n_steps=6000).plot(ax=axes[0])

# (b) the split-step Fourier PE field over the same environment.
underwater.parabolic_equation(50.0, z, c, source_depth=1000.0,
                              max_range=100e3, range_step=25.0,
                              n_depth_points=1024).plot(ax=axes[1])

# (c) modes against PE in a *different*, shallow 200 m waveguide.
r = np.linspace(100.0, 20_000.0, 400)
nm = underwater.normal_modes(50.0, [0.0, 200.0], [1500.0, 1530.0],
                             source_depth=30.0, receiver_depth=120.0,
                             ranges_m=r, n_depth_points=800)
pe = underwater.parabolic_equation(50.0, [0.0, 200.0], [1500.0, 1530.0],
                                   source_depth=30.0, max_range=20e3,
                                   range_step=20.0, n_depth_points=512)
zi = int(np.argmin(np.abs(pe.depths - 120.0)))
axes[2].plot(nm.ranges / 1000.0, nm.propagation_loss, label="Normal modes")
axes[2].plot(pe.ranges / 1000.0, pe.propagation_loss[zi], label="PE")
axes[2].invert_yaxis()
axes[2].set(xlabel="Range [km]", ylabel="Propagation loss [dB]")
axes[2].legend()
plt.show()

import numpy as np
from phonometry import underwater

# A Munk deep-water profile.
z = np.linspace(0.0, 5000.0, 60)
eta = 2.0 * (z - 1300.0) / 1300.0
c = 1500.0 * (1.0 + 0.00737 * (eta - 1.0 + np.exp(-eta)))

rays = underwater.ray_trace(z, c, source_depth=1000.0,
                    launch_angles_deg=np.linspace(-12.0, 12.0, 21), max_range=100e3)
rays.plot()   # ray paths / convergence zones (needs matplotlib)

# Shallow isovelocity waveguide: modes and PE.
modes = underwater.normal_modes(50.0, [0.0, 200.0], [1500.0, 1500.0],
                        source_depth=50.0, receiver_depth=100.0)
print(modes.wavenumbers.size, "propagating modes")
field = underwater.parabolic_equation(50.0, [0.0, 200.0], [1500.0, 1500.0],
                              source_depth=50.0, max_range=20e3)
field.plot()  # PL field over range x depth (needs matplotlib)

normal_modes returns a NormalModeResult (wavenumbers, mode_functions, propagation_loss); ray_trace a RayTraceResult (ranges, depths and travel_times per ray); parabolic_equation a ParabolicEquationResult (the propagation_loss field). The three differ in what they approximate, not in what they model, so choosing between them is a regime question — the table in section 5 is the short answer.

2. Normal modes: the waveguide as a sum of standing waves

In a horizontally stratified ocean the Helmholtz equation separates in cylindrical coordinates, $p(r,z)=\Phi (r)\Psi (z)$, and the depth factor obeys a Sturm-Liouville eigenvalue problem (Jensen Eq. 5.3):

$$\frac{d^2{\Psi }_m}{d{z}^2}+\left[\frac{{\omega }^2}{c^2(z)}-k_{rm}^2\right]{\Psi }_m=0,{\Psi }_m(0)=0,$$

with a pressure-release surface at $z=0$ and, at the bottom $z=D$, $\Psi (D)=0$ for a pressure-release bed or $d\Psi /dz{|}_D=0$ for a rigid one. Each eigenfunction ${\Psi }_m(z)$ is a standing wave in depth that travels in range as $e^{i{k}_{rm}r}$ with its own horizontal wavenumber $k_{rm}$; only the modes with real $k_{rm}$ propagate, the rest are evanescent and die within a few water depths. The field is the modal sum (Eq. 5.14),

$$p(r,z)\simeq \frac{i{e}^{-i\pi /4}}{\rho (z_s)\sqrt{8\pi r}}\sum _m{\Psi }_m(z_s){\Psi }_m(z)\frac{e^{i{k}_{rm}r}}{\sqrt{k_{rm}}},$$

each mode weighted by its excitation at the source depth ${\Psi }_m(z_s)$ and its amplitude at the receiver depth ${\Psi }_m(z)$, and the coherent propagation loss follows as $PL=-20{\log }_{10}|p(r,z)/p_0(1\mathrm{m})|$ (Eq. 5.15). normal_modes discretises the depth equation by finite differences (a symmetric tridiagonal eigenproblem) on a grid refined enough to keep the near-cutoff eigenvalues honest, and warns when a retained mode sits too close to its discretisation error band.

The mode count is the physics: an isovelocity channel of depth $D$ carries $M\approx kD/\pi$ propagating modes, so a 200 m channel at 50 Hz carries 13 and the same channel at 5 kHz would carry over 1300 — which is the frequency at which rays become the economical description. Low frequency and shallow water mean few modes and a compact, essentially exact description. The ideal pressure-release waveguide is the validation oracle, with $k_{rm}=\sqrt{k^2-(m\pi /D{)}^2}$ in closed form:

import numpy as np
from phonometry import underwater

# A 200 m isovelocity channel at 50 Hz: kD/pi = 2 f D / c = 13.3.
modes = underwater.normal_modes(50.0, [0.0, 200.0], [1500.0, 1500.0],
                                source_depth=50.0, receiver_depth=100.0)
print(modes.wavenumbers.size)                   # 13  propagating modes
k = 2 * np.pi * 50.0 / 1500.0
print(round(float(modes.wavenumbers[0]), 5))    # 0.20885  computed kr1
print(round(np.sqrt(k**2 - (np.pi / 200.0) ** 2), 5))   # 0.20885  exact

Mode $m$ has $m-1$ interior nulls, and the source depth decides which modes exist in the field at all: the fourth mode has a null at exactly 50 m in a 200 m channel, so a source there does not excite it. The staircase in the middle is the mode count cutting on one mode at a time as frequency rises, tracking $M\approx kD/\pi$; the right-hand panel is the interference of the 13 modes that survive at 50 Hz. A computed result draws all three with modes.plot().

Show the code for this figure

import matplotlib.pyplot as plt
import numpy as np

fig, axes = plt.subplots(1, 3, figsize=(16, 5.4))
r = np.linspace(100.0, 20_000.0, 500)
res = underwater.normal_modes(50.0, [0.0, 200.0], [1500.0, 1500.0],
                              source_depth=50.0, receiver_depth=100.0,
                              ranges_m=r, n_depth_points=800)
for m in range(4):
    axes[0].plot(res.mode_functions[m], res.mode_depths, label=f"m = {m + 1}")
axes[0].axhline(50.0, linestyle="--")
axes[0].invert_yaxis()
axes[0].set(xlabel="Mode function Psi_m(z)", ylabel="Depth [m]")
axes[0].legend()

freqs = np.arange(10.0, 201.0, 1.0)
counts = [underwater.normal_modes(f, [0.0, 200.0], [1500.0, 1500.0],
                                  source_depth=50.0, receiver_depth=100.0,
                                  n_depth_points=1200).wavenumbers.size
          for f in freqs]
axes[1].step(freqs, counts, where="post", label="modes returned")
axes[1].plot(freqs, 2.0 * freqs * 200.0 / 1500.0, "--", label="M = kD/pi")
axes[1].set(xlabel="Frequency [Hz]", ylabel="Propagating modes M")
axes[1].legend()

res.plot(ax=axes[2])
plt.show()

That ideal waveguide plays two roles, and they should not be confused. It is the exact oracle the solver is checked against — which is why it is used here — and it is also a physically extreme bottom: a pressure-release bed reflects with unit magnitude and a phase inversion, so it transmits nothing into the sediment and absorbs nothing. A modal or PE propagation loss computed over it therefore understates the loss of any real seabed, and the gap widens with range because it accumulates over bounces. The rigid option is the opposite extreme, offered for the same validation reason and not as a sediment model. For a real bottom in shallow water, price the seabed with the Weston regimes of Underwater sound propagation, whose critical angle and reflection-loss gradient encode exactly what these boundary conditions omit, and read the modal answer as an upper bound on the received level.

3. Ray tracing: turning points and travel times

In the high-frequency limit the Helmholtz equation collapses to the eikonal equation, and its characteristics are rays: trajectories integrated from the first-order system (Jensen Eqs. 3.23-3.24)

$$\frac{dr}{ds}=c\xi ,\frac{dz}{ds}=c\zeta ,\frac{d\zeta }{ds}=-\frac{1}{c^2}\frac{\partial c}{\partial z},$$

with $s$ the arc length and $(\xi ,\zeta )$ the ray slowness. In a range-independent profile the horizontal slowness is conserved along each ray, which is Snell’s law in continuous form, $\cos \theta (z)/c(z)=\cos {\theta }_0/c(z_0)$: a ray bends toward lower sound speed, flattens as $c$ grows, and turns where $c(z_t)=c(z_s)/\cos {\theta }_0$. In a linear gradient $c(z)=c_0+gz$ the arcs are exactly circular with radius $R=c_0/(g\cos {\theta }_0)$, the closed form the solver is validated against; in a deep-water profile the family of rays refocuses periodically into the convergence zones of the section 1 figure. ray_trace integrates all launch angles at once with a fixed-step fourth-order Runge-Kutta scheme, reflecting at the surface and the bottom, and carries the travel time along with them as a third state of the same step ($dt/dr=1/(\xi c^2)$, with $\xi =\cos {\theta }_0/c(z_s)$ the Snell invariant):

import numpy as np
from phonometry import underwater

# An isothermal deep layer: c rises 0.017 (m/s)/m with pressure, so a ray
# launched 6 degrees downward from 100 m turns back up where Snell gives
# c(z_t) = c(z_s)/cos(6 deg).
z = [0.0, 1000.0]
c = [1490.0, 1507.0]                        # linear gradient, g = 0.017 1/s
rays = underwater.ray_trace(z, c, source_depth=100.0,
                            launch_angles_deg=[6.0], max_range=40e3,
                            n_steps=20000)
z_turn = (1491.7 / np.cos(np.radians(6.0)) - 1490.0) / 0.017
print(round(z_turn, 1))                     # 583.3  analytic turning depth
print(round(float(rays.depths.max()), 1))   # 583.3  traced

# The arc also fixes the time to that turn: the sound speed cancels and only
# the launch angle and the gradient are left,
# t = (1/g) ln[(1 + sin theta_0) / cos theta_0].
th = np.radians(6.0)
r_turn = np.sin(th) / (np.cos(th) / 1491.7 * 0.017)   # 9222.6 m
t_turn = np.log((1.0 + np.sin(th)) / np.cos(th)) / 0.017
print(round(t_turn, 3))                                        # 6.171  analytic
traced = np.interp(r_turn, rays.ranges[0], rays.travel_times[0])
print(round(float(traced), 3))                                 # 6.171  traced

Steeper launch means deeper turn: in the same water the 2° ray turns at 154 m and the 10° ray at 1454 m. The 6° ray carries the page’s own validation, the 583.3 m analytic turning depth, which the trace reproduces to the tenth of a metre it prints, and the dashed overlay is the exact circular arc of radius $R=c_0/(g\cos {\theta }_0)$ the solver is checked against. A result draws its own fan with rays.plot().

Show the code for this figure

import matplotlib.pyplot as plt
import numpy as np

c_top, grad, z_source = 1490.0, 0.017, 100.0
c_source = c_top + grad * z_source
fig, (ax_c, ax_r) = plt.subplots(1, 2, figsize=(13, 6),
                                 gridspec_kw={"width_ratios": [1, 3.4]})
zz = np.linspace(0.0, 2000.0, 200)
ax_c.plot(c_top + grad * zz, zz)
ax_c.invert_yaxis()
ax_c.set(xlabel="c(z) [m/s]", ylabel="Depth [m]")

for angle in (2.0, 4.0, 6.0, 8.0, 10.0):
    rays = underwater.ray_trace([0.0, 2000.0], [c_top, c_top + grad * 2000.0],
                                source_depth=z_source,
                                launch_angles_deg=[angle], max_range=32e3,
                                n_steps=20000)
    ax_r.plot(rays.ranges[0] / 1000.0, rays.depths[0], label=f"{angle:.0f} deg")
    z_turn = (c_source / np.cos(np.radians(angle)) - c_top) / grad
    print(f"{angle:4.0f} deg -> analytic {z_turn:7.1f} m, "
          f"traced {float(rays.depths.max()):7.1f} m")
ax_r.set_ylim(1760.0, -70.0)
ax_r.set(xlabel="Range [km]", ylabel="Depth [m]")
ax_r.legend()
plt.show()

Rays buy geometry and timing: paths, turning depths, convergence-zone ranges and the travel time along every one of them, at a cost independent of frequency. Because the time rides the same four Runge-Kutta stages as the trajectory, and takes its sound speed from the interpolation those stages already do, it describes the path actually returned rather than a second reading of it, and it matches the closed form for a constant gradient (Medwin & Clay 1998, Eq. (3.3.20)) to about $10^{-14}$ s. What rays do not carry here is a full amplitude: the geometric ray-tube intensity diverges at caustics, so the result object leaves the level to the modal or PE field.

4. The parabolic equation: a one-way field, marched in range

The parabolic equation trades the boundary-value Helmholtz problem for an initial-value problem in range. Factor out the fast outgoing oscillation, $p(r,z)\simeq \psi (r,z)e^{i{k}_{0}r}/\sqrt{k_{0}r}$ with a reference wavenumber $k_0=\omega /c_0$, and for energy travelling within a small angle of the horizontal the envelope obeys the standard (Tappert) PE (Jensen §6.2):

$$2i{k}_0\frac{\partial \psi }{\partial r}+\frac{{\partial }^2\psi }{\partial z^2}+k_0^2\left[n^2(z)-1\right]\psi =0,n(z)=\frac{c_0}{c(z)}.$$

The split-step Fourier algorithm marches it by operator splitting, alternating two individually exact half-physics steps: diffraction is a multiplication by $e^{-i{k}_z^2\Delta r/2k_0}$ in the vertical-wavenumber domain, and refraction a phase screen $e^{i{k}_0(n^2-1)\Delta r/2}$ back in depth, with one transform pair per range step. parabolic_equation starts from a Gaussian field matched to the point source and uses a discrete sine transform in depth, which enforces the pressure-release surface and bottom by construction. The price is the paraxial approximation: the standard PE is accurate within roughly ±15-20° of the horizontal, and steeper energy carries a phase error that shows at short range in shallow waveguides (Jensen §6.2). The free-field calibration is the oracle: with no gradient at all, the marched field must reproduce spherical spreading,

import numpy as np
from phonometry import underwater

# Free field: in a 5000 m isovelocity column, before any boundary is felt,
# the PE must reproduce spherical spreading, PL = 20 lg R.
field = underwater.parabolic_equation(50.0, [0.0, 5000.0], [1500.0, 1500.0],
                                      source_depth=2500.0, max_range=2000.0,
                                      range_step=10.0)
iz = np.argmin(np.abs(field.depths - 2500.0))
ir = np.argmin(np.abs(field.ranges - 1000.0))
print(round(float(field.propagation_loss[iz, ir]), 2))   # 60.0 = 20 lg 1000

and it does so to about $10^{-4}$ dB at the default range step. That is the check the PE cannot fail: free field, no boundaries, no steep energy. The paraxial limit is only visible where those three are absent.

In a 100 m waveguide at 50 Hz, four of the six propagating modes leave the source at 27°, 37°, 49° and 64° — outside the band the standard PE is written for. The consequence is on the left: the range-averaged PE loss sits 1.6 dB above the modal reference, at 300 m and still at 5.7 km. The offset does not wash out with range, because an ideal waveguide is lossless and strips nothing away; in a real seabed the steep modes attenuate first and the gap does close. This is the failure mode to look for in your own output — a level that is uniformly a little too low, not a wrong shape.

Show the code for this figure

import matplotlib.pyplot as plt
import numpy as np

r = np.linspace(50.0, 6000.0, 700)
modes = underwater.normal_modes(50.0, [0.0, 100.0], [1500.0, 1500.0],
                                source_depth=30.0, receiver_depth=60.0,
                                ranges_m=r, n_depth_points=800)
pe_field = underwater.parabolic_equation(50.0, [0.0, 100.0], [1500.0, 1500.0],
                                         source_depth=30.0, max_range=6000.0,
                                         range_step=5.0, n_depth_points=1024)
iz2 = int(np.argmin(np.abs(pe_field.depths - 60.0)))
pe_pl = np.interp(r, pe_field.ranges, pe_field.propagation_loss[iz2])

def range_average(pl, width=61):        # incoherent, in the energy domain
    return -10.0 * np.log10(np.convolve(10.0 ** (-pl / 10.0),
                                        np.ones(width) / width, mode="same"))

inner = (r > 300.0) & (r < 5700.0)
sm_nm = range_average(np.asarray(modes.propagation_loss))
sm_pe = range_average(pe_pl)
print(round(float(np.mean(sm_pe[inner] - sm_nm[inner])), 2))   # 1.64 dB

fig, (ax_pl, ax_ang) = plt.subplots(1, 2, figsize=(13.5, 5.6),
                                    gridspec_kw={"width_ratios": [2, 1]})
ax_pl.plot(r / 1000.0, modes.propagation_loss, alpha=0.35)
ax_pl.plot(r / 1000.0, pe_pl, alpha=0.35)
ax_pl.plot(r[inner] / 1000.0, sm_nm[inner], label="Normal modes, averaged")
ax_pl.plot(r[inner] / 1000.0, sm_pe[inner], label="PE, averaged")
ax_pl.set_ylim(30.0, 92.0)
ax_pl.invert_yaxis()
ax_pl.set(xlabel="Range [km]", ylabel="Propagation loss [dB]")
ax_pl.legend()

k0 = 2.0 * np.pi * 50.0 / 1500.0
grazing = np.degrees(np.arccos(np.clip(modes.wavenumbers / k0, -1.0, 1.0)))
ax_ang.bar(np.arange(1, grazing.size + 1), grazing)
ax_ang.axhline(20.0, linestyle="--")
ax_ang.set(xlabel="Mode index m", ylabel="Grazing angle arccos(k_rm/k) [deg]")
plt.show()

Setting up a run

The three solvers are configured by the caller, and every snippet above sets a grid without saying why. The rules are short.

The depth grid. The sine transform of the PE, and the finite-difference eigenproblem of the modes, both alias a spectrum that is steep in depth, so the grid is sized in wavelengths. At 50 Hz, $\lambda =30$ m, and a grid near $\lambda /10$ over 5000 m of water means a few thousand interior points — which is why the section 1 snippet’s n_depth_points=512 is labelled the coarse, fast setting rather than the right one. Refine and re-run before trusting a number.

The split-step range step. The refraction half-step is a phase screen $e^{i{k}_0(n^2-1)\Delta r/2}$, so it should turn the phase by well under a radian per step, and be tightened wherever the profile has strong gradients. The honest test is the same one: halve $\Delta r$ and see whether the field moves.

When normal_modes warns. The warning says a retained mode sits inside the finite-difference error band. The remedy is to raise n_depth_points: the eigenvalue error is $O(\Delta z^2)$, and near-cutoff modes are its first casualties — which are exactly the shallow-angle modes that carry the field to long range, so ignoring the warning quietly costs the answer at the ranges you care about.

The ray fan and the ray step. The launch angles are not a free choice either: pick them from Snell’s turning condition so the fan spans the rays the channel traps, rather than the round $\pm 12°$ of the snippet above. The number of steps has to keep the arc length per step small against the vertical scale of the profile, or a turning point is stepped straight over.

The convergence discipline. Refine, re-run, compare; only then compare solvers. Agreement between the modes and the PE, which section 5 calls the practical convergence test, is worth nothing if neither of them has converged on its own grid.

5. Choosing a model

Every propagation function of the underwater module answers the same question, “how much level survives the path”, at a different price in physics. Terminology throughout follows ISO 18405:2017 (propagation loss, source level, levels re 1 µPa). The sound-speed and absorption models named below are implemented and referenced in Underwater sound propagation.

Sound speed. The four equations agree to within about 1 m/s inside the domain they share — 0.98 m/s at 25 °C, 35 ppt and 1000 m, and 0.52 m/s over the upper kilometre of a full profile, as the model comparison on the propagation page shows — so the choice is about validity range, not accuracy. The default UNESCO / Chen-Millero form (as recast by Wong & Zhu 1995) covers 0–40 °C, 0–40 ppt and 0–1000 bar, the widest envelope, and is the international standard. Del Grosso (1974) is restricted to 0–30 °C and 30–40 ppt but is preferred by some authors for deep-ocean work inside that domain (much of the SOFAR-channel literature uses it). Mackenzie (1981) trades pressure for depth directly (2–30 °C, 25–40 ppt, 0–8000 m), which makes it the convenient choice when you have an echo-sounder depth rather than a CTD pressure; the other two convert depth to pressure through Leroy & Parthiot (1998) internally. Medwin (1975) is a six-term simplification valid to about 1000 m in shallow, warm water; it is the coarsest of the four and the one behind the rules of thumb $\partial c/\partial T\approx 4.6-0.110T$ m/s per °C and $\partial c/\partial z\approx 0.016$ m/s per metre, so keep it for mental arithmetic and hand checks, not for a budget.

Absorption. Francois–Garrison (1982) is the reference and the default: it carries the boric-acid, magnesium-sulfate and pure-water relaxations with their full temperature, salinity, depth and pH-implicit dependences, and is trusted from about 100 Hz to 1 MHz. Ainslie–McColm (1998) is a deliberate simplification of the same physics that stays within about 10 % of it across that range; use it when a legible formula matters more than the last percent. Thorp (1967) depends on frequency only (it bakes in 4 °C water near 1000 m) and predates both; keep it for quick low-frequency estimates below a few tens of kHz and for comparison with older literature that used it.

Spreading law. Spherical spreading ($20{\log }_{10}R$) describes a wavefront that expands freely in three dimensions, before any boundary confines it; cylindrical spreading ($10{\log }_{10}R$) describes energy trapped between the surface and the bottom (or in the SOFAR channel) that can only expand in range. The "practical" law splices the two at a transition range $R_0$, which is physically of the order of the water (or channel) depth: spherical while the wavefront has not yet filled the duct, cylindrical once it has. In the 10 kHz example of the propagation-loss section the choice is not cosmetic: against the same figure of merit of 87 dB, spherical-only spreading predicts detection out to about 8.7 km while the practical law with $R_0=1000$ m stretches it to about 15.8 km. When the spreading law is the biggest uncertainty in the budget, that is the cue to stop using a closed form and compute the field.

Closed form or solver. The closed-form propagation loss knows nothing of the sound-speed profile, the seabed or the surface; it is honest for short, direct, boundary-free paths and for first-cut sonar budgets. When refraction and boundaries decide the answer, pick the solver by frequency and geometry (Jensen et al. 2011, Ch. 1):

Solver	Natural regime	What it buys you
ray_trace	High frequency (water depth ≫ λ), deep water	Ray-path geometry, turning depths, travel times, convergence zones; cost independent of frequency
normal_modes	Low frequency, shallow water, range-independent	Finite-difference modal sum with few propagating modes ($M\approx kD/\pi$); the reference solution for its regime, validated against the ideal waveguide’s exact modes
parabolic_equation	Low frequency, long one-way paths	Full-field PL($z$,$r$) with refraction, marched in range over the range-independent $c(z)$ all three solvers assume

The boundaries blur in practice: rays remain usable at surprisingly low frequencies for travel-time work, and the PE remains the workhorse well above its formal small-angle regime. When two of the three agree on a case, as the modes and the PE do in the section 1 figure, that agreement is the practical convergence test.

6. A worked sonar budget

Chain the pieces end to end: a 140 dB re 1 µPa²/Hz source at 10 kHz, a 60 dB ambient spectrum level, a 15 dB array gain and an 8 dB detection threshold give the figure of merit $FOM=140-(60-15)-8=87$ dB, computed by passive_sonar_equation in the sonar-equation section of the propagation guide. Every term then moves the crossing of the propagation-loss curve of that guide’s propagation-loss section (10 °C, 35 ppt, 100 m), and that crossing is the predicted detection range.

import numpy as np

env = {"temperature": 10.0, "salinity": 35.0, "depth": 100.0}
se = underwater.passive_sonar_equation(
    source_level=140.0, propagation_loss=np.linspace(40.0, 120.0, 400),
    noise_level=60.0, directivity_index=15.0, detection_threshold=8.0)
print(se.figure_of_merit)                                    # 87.0 dB

for law, r0 in (("spherical", None), ("practical", 1000.0)):
    res = underwater.detection_range(se.figure_of_merit, 10e3, law=law,
                                     transition_range=r0, **env)
    print(law, round(float(res.detection_range) / 1000.0, 1))  # 8.7 / 15.8 km

# Two ways to lose it: half the array gain, or double the frequency.
trim = underwater.passive_sonar_equation(
    source_level=140.0, propagation_loss=np.linspace(40.0, 120.0, 400),
    noise_level=60.0, directivity_index=7.5, detection_threshold=8.0)
print(trim.figure_of_merit)                                  # 79.5 dB
print(round(float(underwater.detection_range(
    trim.figure_of_merit, 10e3, law="practical",
    transition_range=1000.0, **env).detection_range) / 1000.0, 1))       # 10.0 km
print(round(float(underwater.detection_range(
    se.figure_of_merit, 20e3, law="practical",
    transition_range=1000.0, **env).detection_range) / 1000.0, 1))       # 5.9 km

The same source and the same sea, six answers between 4.1 and 15.8 km. The spreading law is the biggest lever: choosing spherical-only over the practical law with $R_0=1$ km costs 7.1 km at 10 kHz. Halving the array gain costs 5.8 km; doubling the frequency costs 9.9 km, because $\alpha$ goes from 0.95 to 3.30 dB/km. When the spreading law is the biggest uncertainty in a budget, that is the cue to stop using a closed form and compute the field.

Show the code for this figure

import matplotlib.pyplot as plt

ranges = np.linspace(50.0, 30_000.0, 800)
fig, ax = plt.subplots(figsize=(11, 6.4))
cases = (("10 kHz, spherical only", 10e3, "spherical", None),
         ("10 kHz, practical R0 = 1 km", 10e3, "practical", 1000.0),
         ("20 kHz, practical R0 = 1 km", 20e3, "practical", 1000.0))
for label, freq, law, r0 in cases:
    ax.plot(ranges / 1000.0,
            underwater.propagation_loss(ranges, freq, law=law,
                                         transition_range=r0, **env).pl,
            label=label)
for fom in (se.figure_of_merit, trim.figure_of_merit):
    ax.axhline(fom, linewidth=1.4, label=f"FOM = {fom:.1f} dB")
    for _label, freq, law, r0 in cases:
        r50 = float(underwater.detection_range(fom, freq, law=law,
                                               transition_range=r0,
                                               **env).detection_range)
        ax.plot([r50 / 1000.0], [fom], "o")
ax.set_ylim(40.0, 126.0)
ax.invert_yaxis()
ax.set(xlabel="Range [km]", ylabel="Propagation loss [dB]")
ax.legend()
plt.show()

One number in that budget deserves a second look. A 60 dB ambient spectrum level at 10 kHz is well above what the sea alone produces: ocean_ambient_noise gives 42 dB re 1 µPa²/Hz there in a 15 knot wind. Sixty decibels therefore describes a noisy site or a self-noise-limited receiver, not a quiet ocean — and noticing that is the difference between a budget and an arithmetic exercise. This coupling between the absorption model, the spreading law and the sonar equation is why they all live in one module.

What this guide covers

• Covered

  normal_modes, ray_trace and parabolic_equation implement the Jensen et al. (2011) numerical solvers for a range-independent stratified ocean: the finite-difference Sturm-Liouville modal sum (Ch. 5, Eqs. 5.3-5.17), the Runge-Kutta ray-trajectory integration (Ch. 3, Eqs. 3.23-3.24) and the split-step Fourier standard PE (Ch. 6), with ISO 18405:2017 terminology throughout. ray_trace also returns the travel time along each ray, integrated as a third state of the same Runge-Kutta step. Each is validated against an exact closed form: the ideal pressure-release waveguide’s modes, the circular-arc rays of a linear gradient together with the published travel time along them (Medwin & Clay 1998, Eq. (3.3.20)), and free-field spherical spreading, plus the mutual agreement of the modal and PE propagation loss.

• Not covered

  All three solvers assume a range-independent water column with a pressure-release (or, for the modes, optionally rigid) boundary: there is no absorbing or elastic bottom, no sediment attenuation and no real bathymetry, so range-dependent problems are out of scope. The ray solver returns paths and travel times, not amplitudes: ray-tube intensity with caustic corrections is not computed, and nothing searches for the eigenrays that join a given source and receiver. The PE is the standard small-angle (Tappert) form, not a wide-angle Padé variant. For the elastic seabed physics these fluid solvers leave out, see Elastic waves and fluid-solid coupling.

See also

• Underwater sound propagation: the closed forms these solvers replace when refraction and boundaries matter, and the sound-speed profiles they consume.
• Underwater acoustics: radiated noise and pile driving: the ISO 18405 reference levels in which every propagation loss here is expressed.
• Atmospheric refraction: rays and the GFPE: the airborne siblings of these solvers, with the same ray bending and a Green’s-function PE marched over ground impedance instead of a seabed.
• 2D FDTD wave simulation: the time-domain alternative behind the SOFAR ducting animation of the propagation guide.
• API reference: underwater.propagation.numerical.

Quick answers

Which underwater propagation solver should I use?

Pick by frequency and geometry (Jensen et al. 2011, Ch. 1): rays for high frequency and deep water (ray-path geometry, travel times and convergence zones at a cost independent of frequency), normal modes for low frequency in shallow water (few propagating modes, $M\approx kD/\pi$, the reference solution for its regime) and the parabolic equation for low-frequency, long one-way paths (the full PL($z$,$r$) field with refraction). When two solvers agree on a case, that agreement is the practical convergence test.

When is a closed-form propagation loss no longer enough?

When refraction or boundaries decide the answer: a sound-speed minimum that traps energy (the SOFAR channel), surface and bottom reflections in shallow water, or a detection range that swings with the choice of spreading law. The closed form $PL=\text{spreading}+\alpha R$ is honest for short, direct, boundary-free paths and first-cut sonar budgets; beyond that, compute the field.

References

• International Organization for Standardization. (2017). Underwater acoustics — Terminology (ISO 18405:2017). The standardized definitions (propagation loss, source level, sound pressure level re 1 µPa) behind the quantities of this page.
• Jensen, F. B., Kuperman, W. A., Porter, M. B., & Schmidt, H. (2011). Computational ocean acoustics (2nd ed.). Springer. https://doi.org/10.1007/978-1-4419-8678-8The reference monograph implemented here: the modal derivation of section 2 (Ch. 5, Eqs. 5.3-5.17), the ray equations of section 3 (Ch. 3, Eqs. 3.23-3.24), the split-step Fourier parabolic equation of section 4 (Ch. 6) and the model-selection guidance of section 5 (Ch. 1).
• Munk, W. H. (1974). Sound channel in an exponentially stratified ocean, with application to SOFAR. The Journal of the Acoustical Society of America, 55(2), 220-226. https://doi.org/10.1121/1.1914492The canonical deep-water sound-speed profile used by the section 1 figure and snippets.