{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Tutorial 07 \u2014 W-band Mie Doppler spectrum in convective rain (Kollias 2002)\n",
"\n",
"Reference: Kollias, P., Albrecht, B. A., and Marks Jr., F. D., 2002.\n",
"*Cloud radar observations of vertical drafts and microphysics in\n",
"convective rain*, J. Geophys. Res., 107 (doi:10.1029/2001JD002033).\n",
"\n",
"At 94 GHz the raindrop backscattering cross-section \u03c3_b(D) is no\n",
"longer Rayleigh \u2014 it rings through a sequence of Mie maxima and\n",
"minima. Because each drop falls at a deterministic terminal velocity\n",
"v_t(D), these Mie features map directly onto the observed Doppler\n",
"spectrum: the first Mie minimum appears at v \u2248 5.9 m/s in still air,\n",
"and any displacement of that feature from its theoretical location\n",
"is a direct measurement of the mean vertical air motion \u2014 independent\n",
"of the drop-size distribution. This notebook reproduces the paper's\n",
"Figures 1\u20133 and demonstrates the air-motion retrieval.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from rustmatrix import Scatterer, SpectralIntegrator, radar, spectra\n",
"from rustmatrix.psd import ExponentialPSD, PSDIntegrator\n",
"from rustmatrix.refractive import m_w_20C\n",
"from rustmatrix.tmatrix_aux import (K_w_sqr, dsr_thurai_2007,\n",
" geom_vert_back, wl_W)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Figure 1 \u2014 single-drop \u03c3_b(D) at 94 GHz\n",
"\n",
"Each drop's backscatter cross-section at vertical incidence, as a\n",
"function of its equivalent diameter D. The Rayleigh D\u2076 regime holds\n",
"only below about 0.7 mm; beyond that the curve oscillates through\n",
"distinct Mie resonances.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"def build_drop(D_mm, oblate=True):\n",
" axis_ratio = 1.0 / dsr_thurai_2007(D_mm) if oblate else 1.0\n",
" s = Scatterer(radius=D_mm/2, wavelength=wl_W, m=m_w_20C[wl_W],\n",
" Kw_sqr=K_w_sqr[wl_W], axis_ratio=axis_ratio,\n",
" ddelt=1e-4, ndgs=2)\n",
" s.set_geometry(geom_vert_back)\n",
" return s\n",
"\n",
"D_grid = np.linspace(0.05, 5.0, 400) # mm\n",
"sigma_b_obl = np.array([radar.radar_xsect(build_drop(D, True))\n",
" for D in D_grid])\n",
"sigma_b_sph = np.array([radar.radar_xsect(build_drop(D, False))\n",
" for D in D_grid])\n",
"\n",
"fig, ax = plt.subplots(figsize=(7, 4))\n",
"ax.semilogy(D_grid, sigma_b_obl, 'C0-', label='oblate (Thurai 2007)')\n",
"ax.semilogy(D_grid, sigma_b_sph, 'k--', label='sphere')\n",
"ax.set_xlabel('diameter D [mm]')\n",
"ax.set_ylabel(r'$\\sigma_b$ [mm\u00b2]')\n",
"ax.set_title('Figure 1 \u2014 single-drop backscatter at 94 GHz, vertical incidence')\n",
"ax.legend(); ax.grid(True, which='both', alpha=0.3)\n",
"fig.tight_layout();\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Figure 3 \u2014 Mie oscillations mapped to terminal velocity\n",
"\n",
"Plotting \u03c3_b / (\u03c0 r\u00b2) versus the drop's terminal fall speed collapses\n",
"the Mie ripple structure into the velocity coordinate the radar\n",
"actually measures. Oblate drops fall slightly faster than the spheres\n",
"of the same equivalent diameter, so their first Mie minimum lands at\n",
"v \u2248 5.95 m/s instead of 5.88 m/s \u2014 a 7 cm/s shift that matters when\n",
"you're using it as a zero-air-motion fiducial.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"fall = spectra.fall_speed.atlas_srivastava_sekhon_1973\n",
"v_t = fall(D_grid)\n",
"r_mm = D_grid / 2.0\n",
"norm = np.pi * r_mm ** 2\n",
"\n",
"def first_min(D, sigma_b):\n",
" mask = (D > 1.3) & (D < 2.2)\n",
" i = np.argmin(sigma_b[mask])\n",
" return D[mask][i], v_t[mask][i]\n",
"\n",
"D_min_obl, v_min_obl = first_min(D_grid, sigma_b_obl)\n",
"D_min_sph, v_min_sph = first_min(D_grid, sigma_b_sph)\n",
"print(f'First Mie minimum, oblate: D = {D_min_obl:.3f} mm, '\n",
" f'v = {v_min_obl:.3f} m/s')\n",
"print(f'First Mie minimum, sphere: D = {D_min_sph:.3f} mm, '\n",
" f'v = {v_min_sph:.3f} m/s')\n",
"print(f'Shift: {100 * (v_min_obl - v_min_sph):+.1f} cm/s '\n",
" '(oblate vs sphere)')\n",
"\n",
"fig, ax = plt.subplots(figsize=(7, 4))\n",
"ax.semilogy(v_t, sigma_b_obl / norm, 'C0-', label='oblate spheroids')\n",
"ax.semilogy(v_t, sigma_b_sph / norm, 'k--', label='spheres')\n",
"ax.axvline(v_min_obl, color='C3', alpha=0.6,\n",
" label=f'oblate 1st min: {v_min_obl:.2f} m/s')\n",
"ax.axvline(v_min_sph, color='C2', alpha=0.6, ls=':',\n",
" label=f'sphere 1st min: {v_min_sph:.2f} m/s')\n",
"ax.set_xlim(0, 10); ax.set_ylim(1e-5, 1e1)\n",
"ax.set_xlabel('terminal velocity $v_t(D)$ [m/s]')\n",
"ax.set_ylabel(r'$\\sigma_b / (\\pi r^2)$')\n",
"ax.set_title('Figure 3 \u2014 Mie ripples in velocity space (94 GHz)')\n",
"ax.legend(); ax.grid(True, which='both', alpha=0.3)\n",
"fig.tight_layout();\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Full Doppler spectrum at 94 GHz (cf. Figure 2)\n",
"\n",
"Running an exponential rain DSD through `SpectralIntegrator` yields\n",
"the Mie-modulated Doppler spectrum the paper observes. With enough\n",
"PSD diameter samples, the Mie minima appear as *notches* in sZ_h(v)\n",
"at the velocities we just identified.\n",
"\n",
"### Note on diameter sampling\n",
"With zero turbulence and too few PSD diameters, the spectrum\n",
"*looks* spiky \u2014 every tabulated drop lands in a single velocity bin\n",
"and the bins between stay at zero. This is a sampling artifact,\n",
"**not** a physical feature: densify the diameter grid (bump\n",
"`num_points` from 64 to 256) or add a hair of turbulence and the\n",
"spikes collapse into the smooth Mie-modulated spectrum.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"def build_rain_W(num_points=256):\n",
" s = Scatterer(wavelength=wl_W, m=m_w_20C[wl_W],\n",
" Kw_sqr=K_w_sqr[wl_W], ddelt=1e-4, ndgs=2)\n",
" integ = PSDIntegrator()\n",
" integ.D_max = 5.0\n",
" integ.num_points = num_points\n",
" integ.axis_ratio_func = lambda D: 1.0 / dsr_thurai_2007(D)\n",
" integ.geometries = (geom_vert_back,)\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" s.psd = ExponentialPSD(N0=8e3, Lambda=2.2, D_max=5.0)\n",
" return s\n",
"\n",
"rain64 = build_rain_W(num_points=64)\n",
"rain256 = build_rain_W(num_points=256)\n",
"\n",
"def run(sc, turb, w=0.0):\n",
" return SpectralIntegrator(\n",
" sc, fall_speed=fall, turbulence=turb, w=w,\n",
" v_min=-1.0, v_max=12.0, n_bins=1024,\n",
" geometry_backscatter=geom_vert_back,\n",
" ).run()\n",
"\n",
"r_sparse = run(rain64, spectra.NoTurbulence())\n",
"r_dense = run(rain256, spectra.NoTurbulence())\n",
"r_turb = run(rain256, spectra.GaussianTurbulence(0.1))\n",
"\n",
"def dBZ(x):\n",
" return 10 * np.log10(np.maximum(x, 1e-10))\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 4))\n",
"ax.plot(r_sparse.v, dBZ(r_sparse.sZ_h), 'C1-', alpha=0.6,\n",
" label='num_points=64, no turb (spiky \u2014 sampling artifact)')\n",
"ax.plot(r_dense.v, dBZ(r_dense.sZ_h), 'C0-',\n",
" label='num_points=256, no turb (Mie notches resolved)')\n",
"ax.plot(r_turb.v, dBZ(r_turb.sZ_h), 'C3-', alpha=0.8,\n",
" label='num_points=256, \u03c3_t=0.1 m/s')\n",
"ax.axvline(v_min_obl, color='k', ls=':', alpha=0.6,\n",
" label=f'1st Mie min ({v_min_obl:.2f} m/s)')\n",
"ax.set_xlim(0, 10); ax.set_ylim(-50, 25)\n",
"ax.set_xlabel('Doppler velocity v [m/s]')\n",
"ax.set_ylabel('sZ_h [dBZ / (m/s)]')\n",
"ax.set_title('W-band Doppler spectrum \u2014 Mie notches and sampling artifact')\n",
"ax.legend(fontsize=9); ax.grid(True, alpha=0.3)\n",
"fig.tight_layout();\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Air-motion retrieval\n",
"\n",
"Simulate a 1 m/s downward air motion (positive *w* under our\n",
"convention): the whole spectrum shifts by *w*, so the first Mie\n",
"minimum moves from 5.95 m/s to ~6.95 m/s. Subtracting the theoretical\n",
"still-air Mie-minimum location from the observed one recovers the air\n",
"motion \u2014 this is the Kollias 2002 technique in a nutshell, and it\n",
"doesn't depend on knowing the DSD.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"r_w = run(rain256, spectra.GaussianTurbulence(0.1), w=1.0)\n",
"\n",
"band = (r_w.v > v_min_obl) & (r_w.v < v_min_obl + 2.5)\n",
"v_obs = r_w.v[band][int(np.argmin(r_w.sZ_h[band]))]\n",
"w_retrieved = v_obs - v_min_obl\n",
"print(f'Prescribed air motion w = +1.00 m/s (downward)')\n",
"print(f'Observed 1st Mie min v = {v_obs:.2f} m/s')\n",
"print(f'Still-air 1st Mie min v = {v_min_obl:.2f} m/s')\n",
"print(f'Retrieved w = {w_retrieved:+.2f} m/s')\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 4))\n",
"ax.plot(r_turb.v, dBZ(r_turb.sZ_h), 'C0-', label='w = 0')\n",
"ax.plot(r_w.v, dBZ(r_w.sZ_h), 'C3-', label='w = +1 m/s')\n",
"ax.axvline(v_min_obl, color='k', ls=':', alpha=0.6,\n",
" label=f'still-air min ({v_min_obl:.2f})')\n",
"ax.axvline(v_obs, color='C3', ls=':', alpha=0.6,\n",
" label=f'shifted min ({v_obs:.2f})')\n",
"ax.set_xlim(0, 10); ax.set_ylim(-50, 25)\n",
"ax.set_xlabel('Doppler velocity v [m/s]')\n",
"ax.set_ylabel('sZ_h [dBZ / (m/s)]')\n",
"ax.set_title('Kollias 2002 air-motion retrieval: Mie minimum as a fiducial')\n",
"ax.legend(fontsize=9); ax.grid(True, alpha=0.3)\n",
"fig.tight_layout();\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python",
"version": "3.11"
}
},
"nbformat": 4,
"nbformat_minor": 5
}