Métodos numéricos de propagación submarina

La pérdida de propagación en forma cerrada de la propagación submarina del sonido no sabe nada del perfil de velocidad del sonido, del fondo ni de la superficie. Cuando la refracción y los contornos deciden la respuesta, el campo hay que calcularlo: esta guía cubre los tres métodos numéricos del módulo underwater, la física que cada uno discretiza, y cómo elegir entre ellos y las formas cerradas. Los tres suponen un océano independiente de la distancia (estratificado horizontalmente) con superficie de presión nula, toman como entrada el mismo perfil $c(z)$ y siguen a Jensen, Kuperman, Porter y Schmidt, Computational Ocean Acoustics.

1. Los tres métodos de un vistazo

• normal_modes resuelve el problema de autovalores de Sturm-Liouville separado en profundidad por diferencias finitas y suma los modos propagantes en la pérdida de propagación. Validado contra los modos exactos de la guía de ondas ideal (de presión nula).
• ray_trace integra las ecuaciones de trayectoria de rayos (Runge-Kutta, vectorizadas sobre todos los rayos a la vez) a través de un perfil de velocidad del sonido, reflejando en la superficie y el fondo, y arrastra el tiempo de propagación de cada rayo como un estado más de la misma integración. Validado contra los arcos de círculo de un gradiente lineal y contra la forma cerrada del tiempo de propagación a lo largo de ellos.
• parabolic_equation hace avanzar en distancia la PE estándar (Tappert) con el algoritmo split-step de Fourier. Validado contra la divergencia esférica en campo libre; coincide en tendencia con la pérdida de propagación de los modos normales.

Tres paneles, tres preguntas distintas sobre la misma agua. A la izquierda, rayos en un perfil de Munk cuyo mínimo es de 1500,0 m/s a 1271 m frente a los 1548,5 m/s de la superficie: 48 m/s de diferencia bastan para plegar sobre sí mismo el abanico de ±12° y formar zonas de convergencia, y el rayo de ±12° todavía barre de 143 m a 3997 m de profundidad. En el centro, el mismo entorno como campo de la PE, donde las zonas aparecen como interferencia y no como líneas. A la derecha, los dos métodos sobre una guía de ondas distinta, somera, de 200 m: siguen la misma tendencia con un desvío medio de 1,6 dB, pero la diferencia distancia a distancia llega a 25 dB en los mínimos de interferencia, que es la manera honesta de decir qué significa «coincidir en tendencia».

Mostrar el código de esta figura

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

# Un perfil de Munk de velocidad del sonido en aguas profundas, para los
# paneles de rayos y de PE.
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) rayos a través del perfil de Munk, hasta 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], language="es")

# (b) el campo de la PE split-step de Fourier sobre el mismo entorno.
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], language="es")

# (c) modos frente a PE en una guía de ondas *distinta*, somera, de 200 m.
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="Modos normales")
axes[2].plot(pe.ranges / 1000.0, pe.propagation_loss[zi], label="PE")
axes[2].invert_yaxis()
axes[2].set(xlabel="Distancia [km]", ylabel="Pérdida de propagación [dB]")
axes[2].legend()
plt.show()

import numpy as np
from phonometry import underwater

# Un perfil de Munk de aguas profundas.
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()   # trayectorias de rayos / zonas de convergencia

# Guía de ondas somera isovelocidad: modos y 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, "modos propagantes")
field = underwater.parabolic_equation(50.0, [0.0, 200.0], [1500.0, 1500.0],
                              source_depth=50.0, max_range=20e3)
field.plot()  # campo PL sobre distancia x profundidad

normal_modes devuelve un NormalModeResult (wavenumbers, mode_functions, propagation_loss); ray_trace un RayTraceResult (ranges, depths y travel_times por rayo); parabolic_equation un ParabolicEquationResult (el campo propagation_loss). Los tres difieren en lo que aproximan, no en lo que modelan, así que elegir entre ellos es una cuestión de régimen: la tabla de la sección 5 es la respuesta corta.

2. Modos normales: la guía de ondas como suma de ondas estacionarias

En un océano estratificado horizontalmente la ecuación de Helmholtz se separa en coordenadas cilíndricas, $p(r,z)=\Phi (r)\Psi (z)$, y el factor de profundidad obedece un problema de autovalores de Sturm-Liouville (Jensen, Ec. 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,$$

con superficie de presión nula en $z=0$ y, en el fondo $z=D$, $\Psi (D)=0$ para un lecho de presión nula o $d\Psi /dz{|}_D=0$ para uno rígido. Cada autofunción ${\Psi }_m(z)$ es una onda estacionaria en profundidad que viaja en distancia como $e^{i{k}_{rm}r}$ con su propio número de onda horizontal $k_{rm}$; solo propagan los modos con $k_{rm}$ real, el resto son evanescentes y se extinguen al cabo de unas pocas profundidades de agua. El campo es la suma modal (Ec. 5.14),

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

cada modo ponderado por su excitación en la profundidad de la fuente ${\Psi }_m(z_{\mathrm{s}})$ y su amplitud en la del receptor ${\Psi }_m(z)$, y la pérdida de propagación coherente sale como $PL=-20{\log }_{10}|p(r,z)/p_0(1\mathrm{m})|$ (Ec. 5.15). normal_modes discretiza la ecuación en profundidad por diferencias finitas (un problema de autovalores tridiagonal simétrico) sobre una malla lo bastante fina para mantener honestos los autovalores cercanos al corte, y avisa cuando un modo retenido queda demasiado cerca de su banda de error de discretización.

El número de modos es la física: un canal isovelocidad de profundidad $D$ transporta $M\approx kD/\pi$ modos propagantes, así que un canal de 200 m a 50 Hz transporta 13 y el mismo canal a 5 kHz transportaría más de 1300, que es la frecuencia a la que los rayos pasan a ser la descripción económica. Baja frecuencia y aguas someras significan pocos modos y una descripción compacta y esencialmente exacta. La guía de ondas ideal de presión nula es el oráculo de validación, con $k_{rm}=\sqrt{k^2-(m\pi /D{)}^2}$ en forma cerrada:

import numpy as np
from phonometry import underwater

# Un canal isovelocidad de 200 m a 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  modos propagantes
k = 2 * np.pi * 50.0 / 1500.0
print(round(float(modes.wavenumbers[0]), 5))    # 0.20885  kr1 calculado
print(round(np.sqrt(k**2 - (np.pi / 200.0) ** 2), 5))   # 0.20885  exacto

El modo $m$ tiene $m-1$ ceros interiores, y la profundidad de la fuente decide qué modos existen siquiera en el campo: el cuarto modo tiene un cero justo a 50 m en un canal de 200 m, así que una fuente ahí no lo excita. La escalera del centro es el número de modos entrando en propagación de uno en uno conforme sube la frecuencia, siguiendo $M\approx kD/\pi$; el panel de la derecha es la interferencia de los 13 modos que sobreviven a 50 Hz. Un resultado calculado dibuja los tres con modes.plot().

Mostrar el código de esta figura

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="Función modal Psi_m(z)", ylabel="Profundidad [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="modos devueltos")
axes[1].plot(freqs, 2.0 * freqs * 200.0 / 1500.0, "--", label="M = kD/pi")
axes[1].set(xlabel="Frecuencia [Hz]", ylabel="Modos propagantes M")
axes[1].legend()

res.plot(ax=axes[2], language="es")
plt.show()

Esa guía de ondas ideal cumple dos papeles, y no conviene confundirlos. Es el oráculo exacto contra el que se comprueba el cálculo, que es para lo que se usa aquí, y es además un fondo físicamente extremo: un lecho de presión nula refleja con módulo unidad e inversión de fase, así que no transmite nada al sedimento y no absorbe nada. Una pérdida de propagación modal o PE calculada sobre él subestima, por tanto, la pérdida de cualquier fondo marino real, y la diferencia se ensancha con el alcance porque se acumula rebote a rebote. La opción rígida es el extremo contrario, ofrecida por la misma razón de validación y no como modelo de sedimento. Para un fondo real en aguas someras, tasa el lecho con los regímenes de Weston de propagación submarina del sonido, cuyo ángulo crítico y gradiente de pérdida por reflexión codifican exactamente lo que estas condiciones de contorno omiten, y lee la respuesta modal como una cota superior del nivel recibido.

3. Trazado de rayos: puntos de giro y tiempos de propagación

En el límite de alta frecuencia la ecuación de Helmholtz colapsa en la ecuación eikonal, y sus características son los rayos: trayectorias integradas del sistema de primer orden (Jensen, Ecs. 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},$$

con $s$ la longitud de arco y $(\xi ,\zeta )$ la lentitud del rayo. En un perfil independiente de la distancia la lentitud horizontal se conserva a lo largo de cada rayo, que es la ley de Snell en forma continua, $\cos \theta (z)/c(z)=\cos {\theta }_0/c(z_0)$: un rayo se curva hacia la velocidad del sonido más baja, se aplana según crece $c$ y gira donde $c(z_{\mathrm{t}})=c(z_{\mathrm{s}})/\cos {\theta }_0$. En un gradiente lineal $c(z)=c_0+gz$ los arcos son exactamente circulares con radio $R=c_0/(g\cos {\theta }_0)$, la forma cerrada contra la que se valida el trazado de rayos; en un perfil de aguas profundas la familia de rayos se reenfoca periódicamente en las zonas de convergencia de la figura de la sección 1. ray_trace integra todos los ángulos de lanzamiento a la vez con un esquema Runge-Kutta de cuarto orden y paso fijo, reflejando en la superficie y en el fondo, y arrastra con ellos el tiempo de propagación como un tercer estado del mismo paso ($dt/dr=1/(\xi c^2)$, con $\xi =\cos {\theta }_0/c(z_{\mathrm{s}})$ el invariante de Snell):

import numpy as np
from phonometry import underwater

# Una capa profunda isoterma: c sube 0,017 (m/s)/m con la presión, así que
# un rayo lanzado 6 grados hacia abajo desde 100 m gira de vuelta donde
# Snell da c(z_t) = c(z_s)/cos(6 grados).
z = [0.0, 1000.0]
c = [1490.0, 1507.0]                        # gradiente lineal, 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  profundidad de giro analítica
print(round(float(rays.depths.max()), 1))   # 583.3  trazada

# El arco fija también el tiempo hasta ese giro: la velocidad del sonido se
# cancela y solo quedan el ángulo de lanzamiento y el gradiente,
# 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  analítico
traced = np.interp(r_turn, rays.ranges[0], rays.travel_times[0])
print(round(float(traced), 3))                                 # 6.171  trazado

Cuanto más empinado es el lanzamiento, más hondo gira el rayo: en la misma agua el rayo de 2° gira a 154 m y el de 10° a 1454 m. El rayo de 6° lleva la validación de la propia página, los 583,3 m de profundidad de giro analítica, que el trazado reproduce hasta la décima de metro que imprime, y la superposición discontinua es el arco de círculo exacto de radio $R=c_0/(g\cos {\theta }_0)$ contra el que se comprueba el cálculo. Un resultado dibuja su propio abanico con rays.plot().

Mostrar el código de esta figura

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="Profundidad [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} grados")
    z_turn = (c_source / np.cos(np.radians(angle)) - c_top) / grad
    print(f"{angle:4.0f} grados -> analitica {z_turn:7.1f} m, "
          f"trazada {float(rays.depths.max()):7.1f} m")
ax_r.set_ylim(1760.0, -70.0)
ax_r.set(xlabel="Distancia [km]", ylabel="Profundidad [m]")
ax_r.legend()
plt.show()

Los rayos compran geometría y tiempos: caminos, profundidades de giro, distancias de las zonas de convergencia y el tiempo de propagación a lo largo de cada uno, a un coste independiente de la frecuencia. Como el tiempo viaja en las mismas cuatro etapas Runge-Kutta que la trayectoria, y toma su velocidad del sonido de la interpolación que esas etapas ya hacen, describe el camino que de verdad se devuelve y no una segunda lectura de él, y coincide con la forma cerrada de un gradiente constante (Medwin y Clay 1998, Ec. (3.3.20)) hasta unos $10^{-14}$ s. Lo que aquí no llevan es una amplitud completa: la intensidad geométrica del tubo de rayos diverge en las cáusticas, así que el objeto de resultado deja el nivel al campo modal o de la PE.

4. La ecuación parabólica: un campo de un solo sentido que avanza en distancia

La ecuación parabólica cambia el problema de contorno de Helmholtz por un problema de valores iniciales en distancia. Factoriza la oscilación saliente rápida, $p(r,z)\simeq \psi (r,z)e^{i{k}_{0}r}/\sqrt{k_{0}r}$ con un número de onda de referencia $k_0=\omega /c_0$, y para la energía que viaja dentro de un ángulo pequeño respecto a la horizontal la envolvente obedece la PE estándar (Tappert) (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)}.$$

El algoritmo split-step de Fourier la hace avanzar por descomposición de operadores, alternando dos semipasos físicos, exacto cada uno por separado: la difracción es una multiplicación por $e^{-i{k}_z^2\Delta r/2k_0}$ en el dominio del número de onda vertical, y la refracción una pantalla de fase $e^{i{k}_0(n^2-1)\Delta r/2}$ de vuelta en profundidad, con un par de transformadas por paso de distancia. parabolic_equation arranca de un campo gaussiano adaptado a la fuente puntual y usa una transformada discreta de senos en profundidad, que impone por construcción la superficie y el fondo de presión nula. El precio es la aproximación paraxial: la PE estándar es precisa dentro de ±15-20° de la horizontal, aproximadamente, y la energía más empinada arrastra un error de fase que se nota a corta distancia en guías de ondas someras (Jensen §6.2). La calibración en campo libre es el oráculo: sin gradiente alguno, el campo así avanzado debe reproducir la divergencia esférica,

import numpy as np
from phonometry import underwater

# Campo libre: en una columna isovelocidad de 5000 m, antes de sentir
# ningún contorno, la PE debe reproducir la divergencia esférica,
# 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

y lo hace hasta unos $10^{-4}$ dB con el paso de distancia por defecto. Esa es la comprobación que la PE no puede fallar: campo libre, sin contornos, sin energía empinada. El límite paraxial solo se ve donde faltan esas tres cosas.

En una guía de ondas de 100 m a 50 Hz, cuatro de los seis modos propagantes salen de la fuente a 27°, 37°, 49° y 64°: fuera de la banda para la que está escrita la PE estándar. La consecuencia está a la izquierda: la pérdida de la PE promediada en distancia queda 1,6 dB por encima de la referencia modal, a 300 m y todavía a 5,7 km. El sesgo no se diluye con la distancia, porque una guía de ondas ideal no tiene pérdidas y no descama ningún modo; en un fondo marino real los modos empinados se atenúan antes y la diferencia sí se cierra. Este es el modo de fallo que hay que buscar en la propia salida: un nivel uniformemente un poco bajo, no una forma equivocada.

Mostrar el código de esta figura

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):        # incoherente, en el dominio de la energía
    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="Modos normales, promediado")
ax_pl.plot(r[inner] / 1000.0, sm_pe[inner], label="PE, promediada")
ax_pl.set_ylim(30.0, 92.0)
ax_pl.invert_yaxis()
ax_pl.set(xlabel="Distancia [km]", ylabel="Pérdida de propagación [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="Índice de modo m",
           ylabel="Ángulo rasante arccos(k_rm/k) [grados]")
plt.show()

Preparar un cálculo

Los tres métodos los configura quien llama, y todos los fragmentos anteriores fijan una malla sin decir por qué. Las reglas son cortas.

La malla en profundidad. La transformada de senos de la PE, y el autoproblema en diferencias finitas de los modos, producen ambos aliasing de un espectro que es empinado en profundidad, así que la malla se dimensiona en longitudes de onda. A 50 Hz, $\lambda =30$ m, y una malla cercana a $\lambda /10$ sobre 5000 m de agua son unos cuantos miles de puntos interiores, que es por lo que el n_depth_points=512 del fragmento de la sección 1 está etiquetado como el ajuste grueso y rápido y no como el correcto. Refina y vuelve a calcular antes de fiarte de un número.

El paso en distancia del split-step. El semipaso de refracción es una pantalla de fase $e^{i{k}_0(n^2-1)\Delta r/2}$, así que debe girar la fase bastante menos de un radián por paso, y hay que apretarlo allí donde el perfil tenga gradientes fuertes. La prueba honesta es la misma: divide $\Delta r$ entre dos y mira si el campo se mueve.

Cuando normal_modes avisa. El aviso dice que un modo retenido queda dentro de la banda de error de las diferencias finitas. El remedio es subir n_depth_points: el error de autovalor es $O(\Delta z^2)$, y los modos cercanos al corte son sus primeras víctimas, que son justo los modos de ángulo pequeño que llevan el campo a larga distancia, así que ignorar el aviso cuesta en silencio la respuesta a las distancias que interesan.

El abanico de rayos y el paso de rayo. Los ángulos de lanzamiento tampoco son una elección libre: se toman de la condición de giro de Snell para que el abanico abarque los rayos que el canal atrapa, en vez de los ±12° redondos del fragmento de arriba. El número de pasos tiene que mantener la longitud de arco por paso pequeña frente a la escala vertical del perfil, o un punto de giro se salta de una zancada.

La disciplina de convergencia. Refinar, recalcular, comparar; y solo entonces comparar métodos. Que los modos y la PE coincidan, que es lo que la sección 5 llama la prueba práctica de convergencia, no vale nada si ninguno de los dos ha convergido en su propia malla.

5. Elegir modelo

Todas las funciones de propagación del módulo underwater responden a la misma pregunta, «cuánto nivel sobrevive al camino», a distinto precio en física. La terminología sigue en toda la página la ISO 18405:2017 (pérdida por propagación, nivel de fuente, niveles re 1 µPa). Los modelos de velocidad del sonido y de absorción nombrados abajo están implementados y referenciados en la propagación submarina del sonido.

Velocidad del sonido. Las cuatro ecuaciones coinciden a menos de 1 m/s dentro del dominio que comparten: 0,98 m/s a 25 °C, 35 ppt y 1000 m, y 0,52 m/s en el kilómetro superior de un perfil completo, como muestra la comparación de modelos de la página de propagación, así que la elección va de rango de validez, no de exactitud. La forma UNESCO / Chen-Millero por defecto (en la reformulación de Wong y Zhu de 1995) cubre de 0 a 40 °C, de 0 a 40 ppt y de 0 a 1000 bar, la envolvente más amplia, y es el estándar internacional. Del Grosso (1974) queda restringida a 0-30 °C y 30-40 ppt, pero algunos autores la prefieren para trabajo de océano profundo dentro de ese dominio (buena parte de la literatura del canal SOFAR la usa). Mackenzie (1981) cambia presión por profundidad directamente (2-30 °C, 25-40 ppt, 0-8000 m), lo que la convierte en la opción cómoda cuando se dispone de una profundidad de ecosonda y no de una presión de CTD; las otras dos convierten profundidad a presión internamente mediante Leroy y Parthiot (1998). Medwin (1975) es una simplificación de seis términos válida hasta unos 1000 m en aguas someras y cálidas; es la más basta de las cuatro y la que hay detrás de las reglas prácticas $\partial c/\partial T\approx 4,6-0,110T$ m/s por °C y $\partial c/\partial z\approx 0,016$ m/s por metro, así que resérvala para el cálculo mental y las comprobaciones a mano, no para un presupuesto.

Absorción. Francois-Garrison (1982) es la referencia y el valor por defecto: incorpora las relajaciones del ácido bórico, el sulfato de magnesio y el agua pura con sus dependencias completas de temperatura, salinidad, profundidad y pH implícito, y es fiable desde unos 100 Hz hasta 1 MHz. Ainslie-McColm (1998) es una simplificación deliberada de la misma física que se mantiene a menos de un 10 % de ella en ese rango; úsese cuando una fórmula legible importe más que el último punto porcentual. Thorp (1967) depende solo de la frecuencia (lleva incorporada agua a 4 °C cerca de 1000 m) y es anterior a ambas; resérvese para estimaciones rápidas de baja frecuencia por debajo de unas decenas de kHz y para comparar con literatura antigua que la usó.

Ley de divergencia. La divergencia esférica ($20{\log }_{10}R$) describe un frente de onda que se expande libremente en tres dimensiones, antes de que ningún contorno lo confine; el cilíndrico ($10{\log }_{10}R$) describe energía atrapada entre la superficie y el fondo (o en el canal SOFAR) que solo puede expandirse en distancia. La ley "practical" empalma ambos en una distancia de transición $R_0$, físicamente del orden de la profundidad del agua (o del canal): esférico mientras el frente no ha llenado el conducto, cilíndrico cuando ya lo ha hecho. En el ejemplo de 10 kHz de la sección de pérdida de propagación la elección no es cosmética: frente a la misma figura de mérito de 87 dB, la divergencia solo esférica predice detección hasta unos 8,7 km, mientras que la ley práctica con $R_0=1000$ m la estira hasta unos 15,8 km. Cuando la ley de divergencia es la mayor incertidumbre del presupuesto, esa es la señal de dejar la forma cerrada y calcular el campo.

¿Forma cerrada o cálculo numérico? La pérdida de propagación en forma cerrada no sabe nada del perfil de velocidad del sonido, del fondo ni de la superficie; es honesta para caminos cortos, directos y sin contornos, y para presupuestos de sonar de primera pasada. Cuando la refracción y los contornos deciden la respuesta, elíjase el método por frecuencia y geometría (Jensen et al. 2011, cap. 1):

Método	Régimen natural	Qué aporta
ray_trace	Alta frecuencia (profundidad ≫ λ), océano profundo	Geometría de los caminos de rayos, profundidades de giro, tiempos de propagación, zonas de convergencia; coste independiente de la frecuencia
normal_modes	Baja frecuencia, aguas someras, independiente de la distancia	Suma modal por diferencias finitas con pocos modos propagantes ($M\approx kD/\pi$); la solución de referencia en su régimen, validada frente a los modos exactos de la guía de ondas ideal
parabolic_equation	Baja frecuencia, caminos largos de un solo sentido	PL($z$,$r$) de campo completo con refracción, avanzada en distancia sobre el $c(z)$ independiente de la distancia que suponen los tres métodos

En la práctica las fronteras se difuminan: los rayos siguen siendo útiles a frecuencias sorprendentemente bajas para trabajo de tiempos de propagación, y la PE sigue siendo el caballo de batalla bastante por encima de su régimen formal de ángulos pequeños. Cuando dos de los tres coinciden en un caso, como hacen los modos y la PE en la figura de la sección 1, esa coincidencia es la prueba práctica de convergencia.

6. Un presupuesto de sonar resuelto

Encadenando las piezas de punta a punta: una fuente de 140 dB re 1 µPa²/Hz a 10 kHz, un nivel espectral de ruido ambiental de 60 dB, una ganancia de array de 15 dB y un umbral de detección de 8 dB dan la figura de mérito $FOM=140-(60-15)-8=87$ dB, que calcula passive_sonar_equation en la sección de la ecuación del sonar de la guía de propagación. Cada término mueve después el cruce de la curva de pérdida de propagación de la sección de pérdida de propagación de esa guía (10 °C, 35 ppt, 100 m), y ese cruce es la distancia de detección predicha.

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

# Dos maneras de perderla: la mitad de ganancia de array, o el doble de frecuencia.
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

La misma fuente y el mismo mar, seis respuestas entre 4,1 y 15,8 km. La ley de divergencia es la palanca mayor: elegir la esférica sola en lugar de la práctica con $R_0=1$ km cuesta 7,1 km a 10 kHz. Reducir a la mitad la ganancia de array cuesta 5,8 km; doblar la frecuencia cuesta 9,9 km, porque $\alpha$ pasa de 0,95 a 3,30 dB/km. Cuando la ley de divergencia es la mayor incertidumbre de un presupuesto, esa es la señal de dejar la forma cerrada y calcular el campo.

Mostrar el código de esta figura

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, solo esférica", 10e3, "spherical", None),
         ("10 kHz, práctica R0 = 1 km", 10e3, "practical", 1000.0),
         ("20 kHz, práctica 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="Distancia [km]", ylabel="Pérdida de propagación [dB]")
ax.legend()
plt.show()

Un número de ese presupuesto merece una segunda mirada. Un nivel espectral ambiental de 60 dB a 10 kHz está bastante por encima de lo que produce el mar por sí solo: ocean_ambient_noise da ahí 42 dB re 1 µPa²/Hz con un viento de 15 nudos. Sesenta decibelios describen por tanto un emplazamiento ruidoso o un receptor limitado por su propio ruido, no un océano tranquilo, y darse cuenta de eso es la diferencia entre un presupuesto y un ejercicio de aritmética. Este acoplamiento entre el modelo de absorción, la ley de divergencia y la ecuación del sonar es la razón de que todos vivan en un mismo módulo.

Qué cubre esta guía

• Cubierto

  normal_modes, ray_trace y parabolic_equation implementan los métodos numéricos de Jensen et al. (2011) para un océano estratificado independiente de la distancia: la suma modal de Sturm-Liouville por diferencias finitas (cap. 5, Ecs. 5.3-5.17), la integración Runge-Kutta de las trayectorias de rayos (cap. 3, Ecs. 3.23-3.24) y la PE estándar split-step de Fourier (cap. 6), con la terminología de la ISO 18405:2017 en toda la página. ray_trace devuelve además el tiempo de propagación de cada rayo, integrado como un tercer estado del mismo paso Runge-Kutta. Cada uno está validado contra una forma cerrada exacta: los modos de la guía de ondas ideal de presión nula, los rayos en arco de círculo de un gradiente lineal junto con el tiempo de propagación publicado a lo largo de ellos (Medwin y Clay 1998, Ec. (3.3.20)) y la divergencia esférica en campo libre, más la coincidencia mutua de la pérdida de propagación modal y de la PE.

• No cubierto

  Los tres métodos suponen una columna de agua independiente de la distancia con un contorno de presión nula (o, para los modos, opcionalmente rígido): no hay fondo absorbente ni elástico, ni atenuación del sedimento, ni batimetría real, así que los problemas dependientes de la distancia quedan fuera de alcance. El trazado de rayos devuelve caminos y tiempos de propagación, no amplitudes: la intensidad de tubo de rayos con correcciones de cáusticas no se calcula, y nada busca los eigenrayos que unen una fuente y un receptor dados. La PE es la forma estándar de ángulos pequeños (Tappert), no una variante de gran angular de Padé. Para la física del lecho elástico que estos métodos fluidos dejan fuera, véase Ondas elásticas y acoplamiento fluido-sólido.

Véase también

• Propagación submarina del sonido: las formas cerradas que estos métodos sustituyen cuando la refracción y los contornos importan, y los perfiles de velocidad del sonido que consumen.
• Acústica submarina: ruido radiado e hincado de pilotes: los niveles de referencia de ISO 18405 en los que se expresa toda pérdida de propagación de esta página.
• Refracción atmosférica: rayos y la GFPE: los hermanos aéreos de estos métodos, con la misma curvatura de rayos y una PE de función de Green que avanza sobre impedancia de suelo en lugar de un lecho marino.
• Simulación de ondas FDTD 2D: la alternativa en el dominio del tiempo tras la animación del canal SOFAR de la guía de propagación.
• Referencia de la API: underwater.propagation.numerical.

Respuestas rápidas

¿Qué método de propagación submarina debo usar?

Elige por frecuencia y geometría (Jensen et al. 2011, cap. 1): rayos para alta frecuencia y océano profundo (geometría de los caminos de rayos, tiempos de propagación y zonas de convergencia a un coste independiente de la frecuencia), modos normales para baja frecuencia en aguas someras (pocos modos propagantes, $M\approx kD/\pi$, la solución de referencia en su régimen) y la ecuación parabólica para caminos largos de un solo sentido a baja frecuencia (el campo completo PL($z$,$r$) con refracción). Cuando dos métodos coinciden en un caso, esa coincidencia es la prueba práctica de convergencia.

¿Cuándo deja de bastar la pérdida de propagación en forma cerrada?

Cuando la refracción o los contornos deciden la respuesta: un mínimo de velocidad del sonido que atrapa energía (el canal SOFAR), reflexiones en la superficie y el fondo en aguas someras, o una distancia de detección que oscila con la elección de la ley de divergencia. La forma cerrada $PL=\text{divergencia}+\alpha R$ es honesta para caminos cortos, directos y sin contornos, y para presupuestos de sonar de primera pasada; más allá, calcula el campo.

Referencias

• International Organization for Standardization. (2017). Underwater acoustics — Terminology (ISO 18405:2017). Las definiciones normalizadas (pérdida por propagación, nivel de fuente, nivel de presión acústica re 1 µPa) tras las magnitudes de esta página.
• Jensen, F. B., Kuperman, W. A., Porter, M. B. y Schmidt, H. (2011). Computational ocean acoustics (2.ª ed.). Springer. https://doi.org/10.1007/978-1-4419-8678-8La monografía de referencia implementada aquí: la derivación modal de la sección 2 (cap. 5, Ecs. 5.3-5.17), las ecuaciones de rayos de la sección 3 (cap. 3, Ecs. 3.23-3.24), la ecuación parabólica split-step de Fourier de la sección 4 (cap. 6) y la guía de selección de modelo de la sección 5 (cap. 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.1914492El perfil canónico de velocidad del sonido en aguas profundas usado por la figura y los ejemplos de la sección 1.