Tutorial 13 — Wind + turbulence sensitivity of Doppler moments#
A vertically pointing radar sees a Doppler spectrum broadened by two independent mechanisms:
Turbulence — eddies smaller than the resolution volume rearrange drop velocities, convolving the fall-speed spectrum with a Gaussian of width \(\sigma_t\).
Horizontal wind through a finite beam — any pixel off the boresight contributes a radial component \(u_h\sin\theta\cos(\phi-\phi_w)\) to the line-of-sight. Integrated over a Gaussian beam of one-way HPBW \(\theta_b\) this produces a deterministic Gaussian broadening of width
\[\sigma_\mathrm{beam} = \frac{|u_h|\,\theta_b}{2\sqrt{2\ln 2}}\](Doviak & Zrnić 1993, §5.3). It is the Doppler equivalent of a standard beam-filling correction.
The two widths add in quadrature: \(\sigma^2 = \sigma_t^2 + \sigma_\mathrm{beam}^2\). Neither biases the first moment (both are symmetric about the boresight); both inflate the second moment.
This notebook sweeps \(u_h \in \{0, 5, 10, 20\}\) m/s and \(\theta_b \in \{1°, 3°, 5°\}\) at fixed \(\sigma_t^2 = 0.5\) m²/s² (\(\sigma_t \approx 0.707\) m/s) for a W-band vertical-pointing radar, and tabulates the observed moments against the closed-form \(\sigma_\mathrm{beam}\) prediction.
Scene assumption — horizontal wind is constant in magnitude over the beam (no shear, no spatial gradient). Tutorial 14 takes up the spatially heterogeneous case.
import numpy as np
import matplotlib.pyplot as plt
from rustmatrix import Scatterer, SpectralIntegrator, spectra
from rustmatrix.psd import ExponentialPSD, PSDIntegrator
from rustmatrix.refractive import m_w_10C
from rustmatrix.tmatrix_aux import (K_w_sqr, dsr_thurai_2007,
geom_vert_back, wl_W)
BEAMWIDTHS_DEG = (1.0, 3.0, 5.0)
SIGMA_T_SQ = 0.5
SIGMA_T = float(np.sqrt(SIGMA_T_SQ))
U_H_LIST = (0.0, 5.0, 10.0, 20.0)
V_MIN, V_MAX, N_BINS = -1.0, 14.0, 1024
FWHM = 2.0 * np.sqrt(2.0 * np.log(2.0))
W-band rain scatterer (Marshall-Palmer DSD)#
rain = Scatterer(wavelength=wl_W, m=m_w_10C[wl_W],
Kw_sqr=K_w_sqr[wl_W], ddelt=1e-4, ndgs=2)
integ = PSDIntegrator()
integ.D_max = 5.0
integ.num_points = 64
integ.axis_ratio_func = lambda D: 1.0 / dsr_thurai_2007(D)
integ.geometries = (geom_vert_back,)
rain.psd_integrator = integ
rain.psd_integrator.init_scatter_table(rain)
rain.psd = ExponentialPSD(N0=8000.0, Lambda=2.2, D_max=5.0)
Sweep \(u_h\) and \(\theta_b\)#
We first compute a reference spectrum with no wind and a pencil beam (so only the turbulence Gaussian broadens the native DSD spectrum), then quadrature-add \(\sigma_\mathrm{beam}\) for every \((u_h, \theta_b)\) combination.
def moments(v, sZh):
sZh = np.clip(sZh, 0.0, None)
P = sZh.sum()
if P == 0:
return np.nan, np.nan
mu = float((v * sZh).sum() / P)
var = float(((v - mu) ** 2 * sZh).sum() / P)
return mu, float(np.sqrt(max(var, 0.0)))
def run_case(u_h, hpbw_rad):
return SpectralIntegrator(
rain,
fall_speed=spectra.brandes_et_al_2002,
turbulence=spectra.GaussianTurbulence(sigma_t=SIGMA_T),
v_min=V_MIN, v_max=V_MAX, n_bins=N_BINS,
u_h=u_h, beamwidth=hpbw_rad,
).run()
r_ref = run_case(0.0, 0.0)
mu_ref, sig_ref = moments(r_ref.v, r_ref.sZ_h)
print(f'reference: mu = {mu_ref:.3f} m/s, sigma = {sig_ref:.3f} m/s')
reference: mu = 4.636 m/s, sigma = 1.655 m/s
/tmp/ipykernel_5346/3705155406.py:17: UserWarning: Spectral power extends beyond v_bins: expected range [-2.12, 11.2] m/s, v_bins covers [-1, 14] m/s. Bulk-sum identity will be degraded by the leakage.
).run()
results = {}
for theta_deg in BEAMWIDTHS_DEG:
theta = np.deg2rad(theta_deg)
for u_h in U_H_LIST:
r = run_case(u_h, theta)
mu, sig = moments(r.v, r.sZ_h)
sig_beam = u_h * theta / FWHM
sig_pred = float(np.sqrt(sig_ref ** 2 + sig_beam ** 2))
results[(theta_deg, u_h)] = dict(r=r, mu=mu, sig=sig,
sig_beam=sig_beam, sig_pred=sig_pred)
header = f'{"theta_b":>7} {"u_h":>5} {"sigma_beam":>10} {"sigma_pred":>10} {"mu_obs":>7} {"sigma_obs":>9}'
print(header)
print('-' * len(header))
for (theta_deg, u_h), r in results.items():
print(f'{theta_deg:6.1f}° {u_h:5.1f} {r["sig_beam"]:10.3f} '
f'{r["sig_pred"]:10.3f} {r["mu"]:7.3f} {r["sig"]:9.3f}')
theta_b u_h sigma_beam sigma_pred mu_obs sigma_obs
----------------------------------------------------------
1.0° 0.0 0.000 1.655 4.636 1.655
1.0° 5.0 0.037 1.655 4.636 1.655
1.0° 10.0 0.074 1.656 4.636 1.656
1.0° 20.0 0.148 1.661 4.636 1.661
3.0° 0.0 0.000 1.655 4.636 1.655
3.0° 5.0 0.111 1.658 4.636 1.658
3.0° 10.0 0.222 1.670 4.636 1.670
3.0° 20.0 0.445 1.713 4.636 1.713
5.0° 0.0 0.000 1.655 4.636 1.655
5.0° 5.0 0.185 1.665 4.636 1.665
5.0° 10.0 0.371 1.696 4.636 1.696
5.0° 20.0 0.741 1.813 4.636 1.813
/tmp/ipykernel_5346/3705155406.py:17: UserWarning: Spectral power extends beyond v_bins: expected range [-2.12, 11.2] m/s, v_bins covers [-1, 14] m/s. Bulk-sum identity will be degraded by the leakage.
).run()
/tmp/ipykernel_5346/3705155406.py:17: UserWarning: Spectral power extends beyond v_bins: expected range [-2.13, 11.3] m/s, v_bins covers [-1, 14] m/s. Bulk-sum identity will be degraded by the leakage.
).run()
/tmp/ipykernel_5346/3705155406.py:17: UserWarning: Spectral power extends beyond v_bins: expected range [-2.17, 11.3] m/s, v_bins covers [-1, 14] m/s. Bulk-sum identity will be degraded by the leakage.
).run()
/tmp/ipykernel_5346/3705155406.py:17: UserWarning: Spectral power extends beyond v_bins: expected range [-2.15, 11.3] m/s, v_bins covers [-1, 14] m/s. Bulk-sum identity will be degraded by the leakage.
).run()
/tmp/ipykernel_5346/3705155406.py:17: UserWarning: Spectral power extends beyond v_bins: expected range [-2.22, 11.3] m/s, v_bins covers [-1, 14] m/s. Bulk-sum identity will be degraded by the leakage.
).run()
/tmp/ipykernel_5346/3705155406.py:17: UserWarning: Spectral power extends beyond v_bins: expected range [-2.51, 11.6] m/s, v_bins covers [-1, 14] m/s. Bulk-sum identity will be degraded by the leakage.
).run()
/tmp/ipykernel_5346/3705155406.py:17: UserWarning: Spectral power extends beyond v_bins: expected range [-2.19, 11.3] m/s, v_bins covers [-1, 14] m/s. Bulk-sum identity will be degraded by the leakage.
).run()
/tmp/ipykernel_5346/3705155406.py:17: UserWarning: Spectral power extends beyond v_bins: expected range [-2.39, 11.5] m/s, v_bins covers [-1, 14] m/s. Bulk-sum identity will be degraded by the leakage.
).run()
/tmp/ipykernel_5346/3705155406.py:17: UserWarning: Spectral power extends beyond v_bins: expected range [-3.07, 12.2] m/s, v_bins covers [-1, 14] m/s. Bulk-sum identity will be degraded by the leakage.
).run()
Plot the spectra and the width-scaling curve#
Left — \(sZ_h(v)\) for \(u_h = 0\) vs \(20\) m/s at each beamwidth: the wind broadens the fall-speed peak without shifting it. Right — the observed spectral width traced against the analytic \(\sqrt{\sigma_t^2 + \sigma_\mathrm{beam}^2}\) prediction.
colors = {1.0: 'tab:blue', 3.0: 'tab:orange', 5.0: 'tab:red'}
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4.5))
for theta_deg in BEAMWIDTHS_DEG:
r_quiet = results[(theta_deg, 0.0)]['r']
r_wind = results[(theta_deg, 20.0)]['r']
c = colors[theta_deg]
ax1.semilogy(r_quiet.v, np.clip(r_quiet.sZ_h, 1e-6, None),
c=c, ls='--', alpha=0.7,
label=f'θ_b = {theta_deg:.0f}°, u_h = 0')
ax1.semilogy(r_wind.v, np.clip(r_wind.sZ_h, 1e-6, None),
c=c, lw=1.8,
label=f'θ_b = {theta_deg:.0f}°, u_h = 20 m/s')
ax1.set_xlim(0, 10)
ax1.set_xlabel('Doppler velocity [m/s]')
ax1.set_ylabel('$sZ_h$ [mm$^6$/m$^3$ / (m/s)]')
ax1.set_title('Spectra: wind broadens, does not shift')
ax1.legend(fontsize=8)
ax1.grid(alpha=0.3)
for theta_deg in BEAMWIDTHS_DEG:
uhs = np.array(U_H_LIST)
obs = np.array([results[(theta_deg, u)]['sig'] for u in uhs])
pred = np.array([results[(theta_deg, u)]['sig_pred'] for u in uhs])
c = colors[theta_deg]
ax2.plot(uhs, pred, c=c, lw=1.0, alpha=0.6,
label=f'θ_b = {theta_deg:.0f}° (predicted)')
ax2.plot(uhs, obs, 'o', c=c, ms=8,
label=f'θ_b = {theta_deg:.0f}° (observed)')
ax2.axhline(SIGMA_T, c='k', ls=':', alpha=0.4, label='σ_t only')
ax2.set_xlabel('$u_h$ [m/s]')
ax2.set_ylabel('Spectral width σ [m/s]')
ax2.set_title('σ scales linearly with $u_h\\cdot\\theta_b$')
ax2.legend(fontsize=8)
ax2.grid(alpha=0.3)
fig.tight_layout();
Doppler velocity is unaffected by wind#
The first moment of every spectrum is within a fraction of a velocity bin of the reference value — horizontal wind in a symmetric beam cannot bias \(\mu\). The second moment, by contrast, scales exactly as the Doviak–Zrnić formula predicts.
fig, ax = plt.subplots(figsize=(7, 4))
for theta_deg in BEAMWIDTHS_DEG:
uhs = np.array(U_H_LIST)
mus = np.array([results[(theta_deg, u)]['mu'] for u in uhs])
ax.plot(uhs, mus - mu_ref, 'o-', c=colors[theta_deg],
label=f'θ_b = {theta_deg:.0f}°')
ax.axhline(0, c='k', ls=':', alpha=0.5)
ax.set_xlabel('$u_h$ [m/s]')
ax.set_ylabel(r'$\mu - \mu_\mathrm{ref}$ [m/s]')
ax.set_title('First-moment bias is zero within grid discretisation')
ax.legend()
ax.grid(alpha=0.3);
Take-aways#
Mean Doppler velocity is insensitive to horizontal wind when the beam is circularly symmetric and the scene is uniform — the contributions from the \(+\phi\) and \(-\phi\) sides of the beam cancel exactly.
Spectral width increases in quadrature with \(\sigma_\mathrm{beam} = u_h \theta_b / (2\sqrt{2\ln 2})\). The observed widths match the analytic prediction to the velocity-grid resolution.
Narrow beams are wind-immune. A 1° cloud radar adds only 0.15 m/s to \(\sigma\) even at \(u_h = 20\) m/s — well below the intrinsic DSD width. A 5° beam doubles the apparent width at the same wind, enough to mis-attribute a retrieval.
Scene structure changes the story — when the beam straddles a reflectivity gradient or a sheared updraft, the closed form breaks down and the moments develop real biases. That is Tutorial 14.