{ "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 }