{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Tutorial 11 \u2014 Dual-frequency radar signatures across hydrometeor classes\n",
"\n",
"Reference: Honeyager, R., 2013. *Investigating the use of the T-matrix\n",
"method as a proxy for the discrete dipole approximation*, M.S. thesis,\n",
"Florida State University.\n",
"\n",
"Honeyager's thesis argues that the T-matrix applied to a size-matched\n",
"spheroid, with a dielectric built from an ice / air mixing formula at\n",
"the right *volume fraction*, reproduces the single-scattering\n",
"properties of a geometrically-complex hydrometeor (bullet rosettes,\n",
"plates, dendrites, aggregates) computed with DDA \u2014 provided the size\n",
"parameter \u03c7 = 2\u03c0 r / \u03bb and the effective density are right.\n",
"\n",
"That framing is directly useful for dual-frequency radar work: the\n",
"whole zoo of ice habits collapses, to first order, onto a\n",
"two-parameter family (*effective density \u03c1_eff*, *axis ratio*) that\n",
"`rustmatrix` already supports via `refractive.mi(wl, rho)` and the\n",
"`axis_ratio` keyword. This notebook walks through four representative\n",
"hydrometeor classes:\n",
"\n",
"| class | \u03c1_eff [g/cm\u00b3] | axis ratio |\n",
"|----------------------|---------------|------------|\n",
"| rain | 1.00 (water) | Thurai 2007|\n",
"| low-density aggregate| 0.10 | 0.70 |\n",
"| graupel | 0.50 | 0.90 |\n",
"| high-density ice | 0.90 | 1.00 |\n",
"\n",
"and shows how each class signs itself into single-particle \u03c3_b(D),\n",
"bulk DWR vs. D\u2080, and the spectral DWR profile sDWR(v).\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_10C, mi\n",
"from rustmatrix.tmatrix_aux import (K_w_sqr, dsr_thurai_2007,\n",
" geom_vert_back, wl_X, wl_W)\n",
"\n",
"CLASSES = {\n",
" 'rain': dict(rho=1.00, axis_ratio=None), # Thurai below\n",
" 'low-\u03c1 agg':dict(rho=0.10, axis_ratio=0.70),\n",
" 'graupel': dict(rho=0.50, axis_ratio=0.90),\n",
" 'high-\u03c1 ice':dict(rho=0.90, axis_ratio=1.00),\n",
"}\n",
"\n",
"def refractive_index(wl, cls):\n",
" return m_w_10C[wl] if cls == 'rain' else mi(wl, CLASSES[cls]['rho'])\n",
"\n",
"def class_axis_ratio(cls, D_mm):\n",
" ar = CLASSES[cls]['axis_ratio']\n",
" return (1.0 / dsr_thurai_2007(D_mm)) if ar is None else ar\n",
"\n",
"COLORS = {'rain': 'C0', 'low-\u03c1 agg': 'C2',\n",
" 'graupel': 'C1', 'high-\u03c1 ice': 'C3'}\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Single-particle \u03c3_b(D) at X and W\n",
"\n",
"For each class, tabulate the single-particle backscatter cross-section\n",
"at vertical incidence across 0.2\u20138 mm equivalent diameter. The key\n",
"effects to watch:\n",
"\n",
"* **Rayleigh regime** \u2014 all classes follow \u03c3_b \u221d D\u2076 at small D, but\n",
" offset vertically by |K(\u03c1)|\u00b2 differences. Low-\u03c1 ice sits ~30 dB below\n",
" water at equal D.\n",
"* **Non-Rayleigh onset at W-band** \u2014 the first Mie minimum appears once\n",
" \u03c7 = \u03c0D/\u03bb \u2248 1, i.e. D \u2248 1 mm at W-band. That happens for all four\n",
" classes because \u03c7 is set by size, not \u03c1.\n",
"* **Mie oscillation amplitude** \u2014 *this* is where \u03c1 matters. Denser\n",
" classes (rain, high-\u03c1 ice) show sharp Mie minima/maxima; low-density\n",
" aggregates barely oscillate because their refractive index contrast\n",
" with air is small.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"def sigma_b(wl, cls, D_mm):\n",
" s = Scatterer(radius=D_mm/2, wavelength=wl,\n",
" m=refractive_index(wl, cls),\n",
" Kw_sqr=K_w_sqr[wl],\n",
" axis_ratio=class_axis_ratio(cls, D_mm),\n",
" ddelt=1e-4, ndgs=2)\n",
" s.set_geometry(geom_vert_back)\n",
" return radar.radar_xsect(s)\n",
"\n",
"D_grid = np.linspace(0.2, 8.0, 100)\n",
"sigma = {cls: {b: np.array([sigma_b(wl, cls, D) for D in D_grid])\n",
" for b, wl in [('X', wl_X), ('W', wl_W)]}\n",
" for cls in CLASSES}\n",
"\n",
"fig, (axX, axW) = plt.subplots(1, 2, figsize=(11, 4), sharey=True)\n",
"for cls in CLASSES:\n",
" axX.loglog(D_grid, sigma[cls]['X'], color=COLORS[cls], label=cls)\n",
" axW.loglog(D_grid, sigma[cls]['W'], color=COLORS[cls], label=cls)\n",
"axX.set_title('\u03c3_b(D) at X-band (9 GHz)')\n",
"axW.set_title('\u03c3_b(D) at W-band (94 GHz)')\n",
"for ax in (axX, axW):\n",
" ax.set_xlabel('D [mm]')\n",
" ax.grid(True, which='both', alpha=0.3)\n",
" ax.legend(fontsize=9)\n",
"axX.set_ylabel(r'$\\sigma_b$ [mm\u00b2]')\n",
"fig.tight_layout();\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## DWR(D) \u2014 single-particle dual-wavelength ratio\n",
"\n",
"Re-plotting as DWR(D) = 10\u00b7log\u2081\u2080(\u03c3_b^X / \u03c3_b^W) collapses the key\n",
"information: the *size where each class first walks out of Rayleigh*.\n",
"This is Honeyager's thesis Table 2.2 idea generalised to arbitrary\n",
"habit \u2014 the \u03c7 parameter that sets Rayleigh validity is hidden inside\n",
"\u03c1_eff via the refractive index contrast.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"# Equivalent reflectivity normalisation: Ze = \u03bb\u2074 / (\u03c0\u2075 |K_w|\u00b2) \u00b7 \u03c3_b.\n",
"# In the Rayleigh regime Ze_X = Ze_W, so DWR_{X-W} = 0 dB; any positive\n",
"# excursion is a direct signature of non-Rayleigh scattering at W-band.\n",
"def single_Ze(sig, wl, Kw2):\n",
" return wl**4 / (np.pi**5 * Kw2) * sig\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 4.5))\n",
"dwr_curves_1p = {}\n",
"for cls in CLASSES:\n",
" zeX = single_Ze(sigma[cls]['X'], wl_X, K_w_sqr[wl_X])\n",
" zeW = single_Ze(sigma[cls]['W'], wl_W, K_w_sqr[wl_W])\n",
" dwr = 10 * np.log10(zeX / zeW)\n",
" dwr_curves_1p[cls] = dwr\n",
" ax.plot(D_grid, dwr, color=COLORS[cls], lw=1.8, label=cls)\n",
"ax.axhline(0.0, color='k', ls=':', alpha=0.5, label='Rayleigh (DWR=0)')\n",
"ax.set_xlabel('diameter D [mm]')\n",
"ax.set_ylabel('single-particle DWR$_{X-W}$ [dB]')\n",
"ax.set_title('Where does each habit walk out of Rayleigh?')\n",
"ax.set_xlim(0, 8); ax.grid(True, alpha=0.3); ax.legend(fontsize=9)\n",
"fig.tight_layout();\n",
"\n",
"print('Approximate D at which single-particle DWR crosses +3 dB:')\n",
"for cls, dwr in dwr_curves_1p.items():\n",
" above = np.where(dwr > 3.0)[0]\n",
" label = f'{D_grid[above[0]]:.2f} mm' if len(above) else '> 8 mm (stays Rayleigh)'\n",
" print(f' {cls:<12} {label}')\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bulk DWR vs. median-volume diameter D\u2080\n",
"\n",
"The single-particle curves fold through a PSD in the bulk radar\n",
"measurement. Hold the habit fixed, sweep the exponential PSD slope \u039b\n",
"(so D\u2080 = 3.67/\u039b varies), and plot the resulting bulk DWR. The curve\n",
"tells you how sensitive each habit's dual-frequency signature is to\n",
"the *population-average* particle size.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"def bulk_Zh(wl, cls, N0, lam, Dmax):\n",
" s = Scatterer(wavelength=wl, m=refractive_index(wl, cls),\n",
" Kw_sqr=K_w_sqr[wl], ddelt=1e-4, ndgs=2)\n",
" integ = PSDIntegrator()\n",
" integ.D_max = Dmax; integ.num_points = 128\n",
" integ.geometries = (geom_vert_back,)\n",
" if CLASSES[cls]['axis_ratio'] is None:\n",
" integ.axis_ratio_func = lambda D: 1.0 / dsr_thurai_2007(D)\n",
" else:\n",
" s.axis_ratio = CLASSES[cls]['axis_ratio']\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" s.psd = ExponentialPSD(N0=N0, Lambda=lam, D_max=Dmax)\n",
" s.set_geometry(geom_vert_back)\n",
" return radar.refl(s)\n",
"\n",
"# Hold N0 fixed, sweep \u039b so D0 = 3.67/\u039b spans 0.4\u20133 mm\n",
"Lambdas = np.linspace(1.2, 9.0, 10) # mm\u207b\u00b9\n",
"D0s = 3.67 / Lambdas\n",
"DMAX = 8.0; N0 = 8e3\n",
"\n",
"dwr_curves = {}\n",
"for cls in CLASSES:\n",
" dwr = []\n",
" for lam in Lambdas:\n",
" zX = bulk_Zh(wl_X, cls, N0, lam, DMAX)\n",
" zW = bulk_Zh(wl_W, cls, N0, lam, DMAX)\n",
" dwr.append(10 * np.log10(zX / zW))\n",
" dwr_curves[cls] = np.asarray(dwr)\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 4.5))\n",
"for cls, d in dwr_curves.items():\n",
" ax.plot(D0s, d, color=COLORS[cls], lw=1.8, marker='o',\n",
" markersize=4, label=cls)\n",
"ax.axhline(0.0, color='k', ls=':', alpha=0.5)\n",
"ax.set_xlabel('median-volume diameter $D_0$ [mm]')\n",
"ax.set_ylabel('bulk DWR$_{X-W}$ [dB]')\n",
"ax.set_title('Bulk DWR vs. PSD size \u2014 one line per hydrometeor class')\n",
"ax.grid(True, alpha=0.3); ax.legend(fontsize=9)\n",
"fig.tight_layout();\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Dual-frequency Doppler spectra \u2014 aggregate vs. graupel\n",
"\n",
"Now pull the spectral version together. Two bulk populations that\n",
"*can't* be told apart by reflectivity alone \u2014 a low-density aggregate\n",
"PSD and a graupel PSD tuned to give the same X-band Z_h \u2014 show very\n",
"different sDWR(v) profiles. Aggregates have a mostly-monotonic rise\n",
"because their \u03c3_b(D) curve stays smooth; graupel's spectrum carries\n",
"the Mie-resonance structure of its single-particle \u03c3_b(D) through into\n",
"velocity space, since every drop size maps deterministically to its\n",
"own fall speed.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"def build(cls, N0, lam, Dmax, wl):\n",
" s = Scatterer(wavelength=wl, m=refractive_index(wl, cls),\n",
" Kw_sqr=K_w_sqr[wl], ddelt=1e-4, ndgs=2)\n",
" integ = PSDIntegrator()\n",
" integ.D_max = Dmax; integ.num_points = 256\n",
" integ.geometries = (geom_vert_back,)\n",
" if CLASSES[cls]['axis_ratio'] is None:\n",
" integ.axis_ratio_func = lambda D: 1.0 / dsr_thurai_2007(D)\n",
" else:\n",
" s.axis_ratio = CLASSES[cls]['axis_ratio']\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" s.psd = ExponentialPSD(N0=N0, Lambda=lam, D_max=Dmax)\n",
" return s\n",
"\n",
"# PSDs chosen so both give ~moderate snowfall X-band reflectivity.\n",
"# Aggregate: low density, broader PSD to Dmax=6 mm.\n",
"# Graupel: medium density, narrower PSD (smaller max D).\n",
"agg_X = build('low-\u03c1 agg', N0=2e4, lam=2.0, Dmax=6.0, wl=wl_X)\n",
"agg_W = build('low-\u03c1 agg', N0=2e4, lam=2.0, Dmax=6.0, wl=wl_W)\n",
"gra_X = build('graupel', N0=5e3, lam=3.0, Dmax=4.0, wl=wl_X)\n",
"gra_W = build('graupel', N0=5e3, lam=3.0, Dmax=4.0, wl=wl_W)\n",
"\n",
"fall_agg = spectra.fall_speed.locatelli_hobbs_1974_aggregates\n",
"fall_gra = spectra.fall_speed.locatelli_hobbs_1974_graupel_hex\n",
"turb = spectra.GaussianTurbulence(0.2)\n",
"V_MIN, V_MAX = -0.5, 4.0\n",
"\n",
"def run(sc, fall):\n",
" return SpectralIntegrator(\n",
" sc, fall_speed=fall, turbulence=turb,\n",
" v_min=V_MIN, v_max=V_MAX, n_bins=1024,\n",
" geometry_backscatter=geom_vert_back,\n",
" ).run()\n",
"\n",
"r_agg_X = run(agg_X, fall_agg); r_agg_W = run(agg_W, fall_agg)\n",
"r_gra_X = run(gra_X, fall_gra); r_gra_W = run(gra_W, fall_gra)\n",
"\n",
"def dBZ(x):\n",
" return 10 * np.log10(np.maximum(x, 1e-12))\n",
"\n",
"print(f'{\"population\":<15} {\"Z_h X [dBZ]\":>12} {\"Z_h W [dBZ]\":>12} '\n",
" f'{\"DWR [dB]\":>10}')\n",
"for name, rX, rW in [('low-\u03c1 agg', r_agg_X, r_agg_W),\n",
" ('graupel', r_gra_X, r_gra_W)]:\n",
" zX = 10 * np.log10(np.trapezoid(rX.sZ_h, rX.v))\n",
" zW = 10 * np.log10(np.trapezoid(rW.sZ_h, rW.v))\n",
" print(f'{name:<15} {zX:>12.2f} {zW:>12.2f} {zX - zW:>10.2f}')\n",
"\n",
"fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6.5), sharex=True)\n",
"ax1.plot(r_agg_X.v, dBZ(r_agg_X.sZ_h), 'C2-', label='agg, X-band')\n",
"ax1.plot(r_agg_W.v, dBZ(r_agg_W.sZ_h), 'C2--', label='agg, W-band')\n",
"ax1.plot(r_gra_X.v, dBZ(r_gra_X.sZ_h), 'C1-', label='graupel, X-band')\n",
"ax1.plot(r_gra_W.v, dBZ(r_gra_W.sZ_h), 'C1--', label='graupel, W-band')\n",
"ax1.set_ylabel('sZ_h [dBZ / (m/s)]')\n",
"ax1.set_title('Dual-frequency Doppler spectra')\n",
"ax1.legend(fontsize=9); ax1.grid(True, alpha=0.3)\n",
"\n",
"sdwr_agg = dBZ(r_agg_X.sZ_h) - dBZ(r_agg_W.sZ_h)\n",
"sdwr_gra = dBZ(r_gra_X.sZ_h) - dBZ(r_gra_W.sZ_h)\n",
"mask_agg = (r_agg_X.sZ_h > 1e-8) & (r_agg_W.sZ_h > 1e-10)\n",
"mask_gra = (r_gra_X.sZ_h > 1e-8) & (r_gra_W.sZ_h > 1e-10)\n",
"ax2.plot(r_agg_X.v[mask_agg], sdwr_agg[mask_agg], 'C2-',\n",
" lw=1.6, label='low-\u03c1 agg')\n",
"ax2.plot(r_gra_X.v[mask_gra], sdwr_gra[mask_gra], 'C1-',\n",
" lw=1.6, label='graupel')\n",
"ax2.axhline(0.0, color='k', ls=':', alpha=0.5)\n",
"ax2.set_xlabel('Doppler velocity v [m/s]')\n",
"ax2.set_ylabel('sDWR$_{X-W}$ [dB]')\n",
"ax2.set_title('Spectral DWR \u2014 aggregates smooth, graupel resonant')\n",
"ax2.set_xlim(V_MIN, V_MAX)\n",
"ax2.legend(fontsize=9); ax2.grid(True, alpha=0.3)\n",
"fig.tight_layout();\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Takeaway\n",
"\n",
"Following Honeyager 2013's framing, we've used a single T-matrix\n",
"spheroid code path with four choices of (\u03c1_eff, axis ratio) to span\n",
"four physically distinct hydrometeor classes. Two observations pop\n",
"out that matter for interpreting dual-frequency radar data:\n",
"\n",
"* **The *size* at which DWR departs from zero is set by \u03c7 = \u03c0D/\u03bb**\n",
" (\u2248 1 mm at W-band), not by \u03c1_eff. But **the amplitude and shape**\n",
" of the DWR(D) curve beyond that threshold is a direct fingerprint\n",
" of \u03c1_eff: dense classes oscillate strongly, low-\u03c1 aggregates rise\n",
" smoothly and slowly.\n",
"* **Bulk DWR sweeps out very different trajectories with D\u2080**, so a\n",
" single (Z_h, DWR) measurement can be mapped to a habit-dependent\n",
" D\u2080 estimate provided you commit to a class \u2014 this is exactly the\n",
" lookup-table strategy that papers like Mason et al. 2019 and\n",
" Tridon et al. 2019 build their retrievals on.\n",
"* **Spectral DWR resolves habit ambiguity** that bulk DWR can't:\n",
" graupel's Mie resonances appear as wiggles in sDWR(v), while\n",
" low-\u03c1 aggregates produce a smooth monotonic rise.\n",
"\n",
"Swapping (\u03c1_eff, axis_ratio) for any other habit class (dendrites,\n",
"rimed aggregates, hail) is one `CLASSES` dict entry away.\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python",
"version": "3.11"
}
},
"nbformat": 4,
"nbformat_minor": 5
}