{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Tutorial 10 \u2014 Supercooled liquid water vs. snow at cloud-radar frequencies\n",
"\n",
"Reference: Billault-Roux, A.-C., Georgakaki, P., Grazioli, J.,\n",
"Romanens, G., Sotiropoulou, G., Phillips, V., Nenes, A., and Berne,\n",
"A., 2023. *Distinct secondary ice production processes observed in\n",
"radar Doppler spectra: insights from a case study*, Atmos. Chem.\n",
"Phys., 23, 10207\u201310234 (doi:10.5194/acp-23-10207-2023).\n",
"\n",
"Mixed-phase cloud volumes contain supercooled liquid droplets (SLW)\n",
"and ice particles side by side. The two populations have completely\n",
"different scattering and fall-speed signatures:\n",
"\n",
"* **SLW**: small (\u2272200 \u00b5m), spherical, cold refractive index of\n",
" water, falls at a few cm/s. Rayleigh at both X- and W-band \u21d2 very\n",
" low Z_h, Z_dr \u2248 0 dB, no polarimetric signature.\n",
"* **Snow**: millimetre-scale oblate low-density aggregates, falls at\n",
" ~0.5\u20131 m/s. Moderately non-Rayleigh at W-band \u21d2 higher Z_h, small\n",
" positive Z_dr from the habit anisotropy.\n",
"\n",
"This notebook builds dedicated scatterers for each population, shows\n",
"their \u03c3_b(D) and bulk radar observables at X and W band, and then\n",
"combines them in a `HydroMix` to produce the bimodal Doppler spectrum\n",
"that motivates the Billault-Roux et al. case-study analysis.\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 (HydroMix, MixtureComponent, Scatterer,\n",
" SpectralIntegrator, radar, spectra)\n",
"from rustmatrix.psd import ExponentialPSD, GammaPSD, PSDIntegrator\n",
"from rustmatrix.refractive import m_w_0C, mi\n",
"from rustmatrix.tmatrix_aux import (K_w_sqr, geom_vert_back,\n",
" wl_X, wl_W)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Build a supercooled-liquid-water scatterer\n",
"\n",
"Small spherical drops (D in 10\u2013200 \u00b5m), refractive index at 0\u00b0C,\n",
"log-normal-ish gamma PSD centred around 30 \u00b5m. We use mm units\n",
"throughout to match `rustmatrix` conventions.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"SLW_DMAX = 0.2 # mm\n",
"\n",
"def build_slw(wl):\n",
" s = Scatterer(wavelength=wl, m=m_w_0C[wl],\n",
" Kw_sqr=K_w_sqr[wl], axis_ratio=1.0,\n",
" ddelt=1e-4, ndgs=2)\n",
" integ = PSDIntegrator()\n",
" integ.D_max = SLW_DMAX\n",
" integ.num_points = 128\n",
" integ.geometries = (geom_vert_back,)\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" # Gamma PSD: D0 = 0.03 mm (30 \u00b5m), Nw = 1e11 m\u207b\u00b3 mm\u207b\u00b9, \u03bc = 4\n",
" s.psd = GammaPSD(D0=0.03, Nw=1e11, mu=4, D_max=SLW_DMAX)\n",
" return s\n",
"\n",
"slw_X = build_slw(wl_X)\n",
"slw_W = build_slw(wl_W)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Build a snow scatterer (low-density oblate aggregates)\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"SNOW_DMAX = 8.0\n",
"RHO_SNOW = 0.2\n",
"AXIS_SNOW = 0.6\n",
"\n",
"def build_snow(wl):\n",
" s = Scatterer(wavelength=wl, m=mi(wl, RHO_SNOW),\n",
" Kw_sqr=K_w_sqr[wl], axis_ratio=AXIS_SNOW,\n",
" ddelt=1e-4, ndgs=2)\n",
" integ = PSDIntegrator()\n",
" integ.D_max = SNOW_DMAX\n",
" integ.num_points = 256\n",
" integ.geometries = (geom_vert_back,)\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" s.psd = ExponentialPSD(N0=5e3, Lambda=1.0, D_max=SNOW_DMAX)\n",
" return s\n",
"\n",
"snow_X = build_snow(wl_X)\n",
"snow_W = build_snow(wl_W)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## \u03c3_b(D) \u2014 Rayleigh vs. Mie across the two populations\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"def sigma_b(wl, D_mm, m, axis_ratio):\n",
" s = Scatterer(radius=D_mm/2, wavelength=wl, m=m,\n",
" Kw_sqr=K_w_sqr[wl], axis_ratio=axis_ratio,\n",
" ddelt=1e-4, ndgs=2)\n",
" s.set_geometry(geom_vert_back)\n",
" return radar.radar_xsect(s)\n",
"\n",
"D_slw = np.linspace(0.005, SLW_DMAX, 60)\n",
"D_snow = np.linspace(0.1, SNOW_DMAX, 200)\n",
"\n",
"sb_slw_X = np.array([sigma_b(wl_X, d, m_w_0C[wl_X], 1.0) for d in D_slw])\n",
"sb_slw_W = np.array([sigma_b(wl_W, d, m_w_0C[wl_W], 1.0) for d in D_slw])\n",
"sb_snow_X = np.array([sigma_b(wl_X, d, mi(wl_X, RHO_SNOW), AXIS_SNOW)\n",
" for d in D_snow])\n",
"sb_snow_W = np.array([sigma_b(wl_W, d, mi(wl_W, RHO_SNOW), AXIS_SNOW)\n",
" for d in D_snow])\n",
"\n",
"fig, ax = plt.subplots(figsize=(7.5, 4.5))\n",
"ax.loglog(D_slw, sb_slw_X, 'C0-', label='SLW, X-band')\n",
"ax.loglog(D_slw, sb_slw_W, 'C0--', label='SLW, W-band')\n",
"ax.loglog(D_snow, sb_snow_X, 'C3-', label='snow, X-band')\n",
"ax.loglog(D_snow, sb_snow_W, 'C3--', label='snow, W-band')\n",
"ax.set_xlabel('diameter D [mm]')\n",
"ax.set_ylabel(r'$\\sigma_b$ [mm\u00b2]')\n",
"ax.set_title('Single-particle backscatter: SLW droplets vs. snow aggregates')\n",
"ax.legend(fontsize=9); ax.grid(True, which='both', alpha=0.3)\n",
"fig.tight_layout();\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bulk radar observables at X and W band\n",
"\n",
"Z_h(SLW) is ~40+ dB below Z_h(snow) even though the droplet\n",
"concentration is high \u2014 the D\u2076 Rayleigh penalty at 30-\u00b5m scale kills\n",
"the signal. Z_dr is zero for SLW (spheres) but positive for snow\n",
"(oblate). Dual-wavelength ratio DWR = Z_X \u2212 Z_W is near zero for SLW\n",
"and positive for snow (non-Rayleigh at W-band).\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"def bulk(s):\n",
" s.set_geometry(geom_vert_back)\n",
" return 10 * np.log10(radar.refl(s))\n",
"\n",
"print(f'{\"\":<10} {\"Z_h X [dBZ]\":>12} {\"Z_h W [dBZ]\":>12} '\n",
" f'{\"DWR [dB]\":>10}')\n",
"for name, sX, sW in [('SLW', slw_X, slw_W),\n",
" ('snow', snow_X, snow_W)]:\n",
" zX = bulk(sX); zW = bulk(sW)\n",
" print(f'{name:<10} {zX:>12.2f} {zW:>12.2f} {zX - zW:>10.2f}')\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bimodal Doppler spectrum of the SLW + snow mixture\n",
"\n",
"A `HydroMix` with per-component fall-speed and turbulence captures\n",
"the two populations cleanly: SLW pins to v \u2248 0 with a narrow spread,\n",
"snow sits near 0.5\u20131.0 m/s with a broader tail. Plotted in dBZ space\n",
"the two modes are clearly separable even with realistic noise.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"# Use W-band for the spectrum \u2014 best velocity-wise resolution of the\n",
"# SLW mode, and non-Rayleigh snow signature is most evident here.\n",
"mix = HydroMix([\n",
" MixtureComponent(slw_W, slw_W.psd, 'slw'),\n",
" MixtureComponent(snow_W, snow_W.psd, 'snow'),\n",
"])\n",
"\n",
"integ = SpectralIntegrator(\n",
" mix,\n",
" component_kinematics={\n",
" 'slw': (lambda D: 0.003 * (D / 0.01) ** 2,\n",
" spectra.GaussianTurbulence(0.1)),\n",
" 'snow': (spectra.fall_speed.locatelli_hobbs_1974_aggregates,\n",
" spectra.GaussianTurbulence(0.2)),\n",
" },\n",
" v_min=-1.0, v_max=3.0, n_bins=1024,\n",
" geometry_backscatter=geom_vert_back,\n",
" noise='realistic',\n",
").run()\n",
"\n",
"def dBZ(x):\n",
" return 10 * np.log10(np.maximum(x, 1e-12))\n",
"\n",
"# Run each population on its own (no noise) so each mode's true peak\n",
"# is unambiguous. The HydroMix spectrum is dominated by snow (~50 dB\n",
"# louder than SLW), so looking for the SLW peak on the combined,\n",
"# noise-added spectrum would just pick up the snow tail crossing\n",
"# through the near-zero velocity range.\n",
"slw_only = SpectralIntegrator(\n",
" slw_W,\n",
" fall_speed=lambda D: 0.003 * (D / 0.01) ** 2,\n",
" turbulence=spectra.GaussianTurbulence(0.1),\n",
" v_min=-1.0, v_max=3.0, n_bins=1024,\n",
" geometry_backscatter=geom_vert_back,\n",
").run()\n",
"snow_only = SpectralIntegrator(\n",
" snow_W,\n",
" fall_speed=spectra.fall_speed.locatelli_hobbs_1974_aggregates,\n",
" turbulence=spectra.GaussianTurbulence(0.2),\n",
" v_min=-1.0, v_max=3.0, n_bins=1024,\n",
" geometry_backscatter=geom_vert_back,\n",
").run()\n",
"\n",
"v_slw = slw_only.v[int(np.argmax(slw_only.sZ_h))]\n",
"v_snow = snow_only.v[int(np.argmax(snow_only.sZ_h))]\n",
"print(f'SLW mode peak v = {v_slw:+.3f} m/s')\n",
"print(f'snow mode peak v = {v_snow:+.3f} m/s')\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 4))\n",
"ax.plot(slw_only.v, dBZ(slw_only.sZ_h), 'C2-', lw=1.3, label='SLW only')\n",
"ax.plot(snow_only.v, dBZ(snow_only.sZ_h), 'C3-', lw=1.3, label='snow only')\n",
"ax.plot(integ.v, dBZ(integ.sZ_h), 'C0-', lw=2.0,\n",
" label='SLW + snow (HydroMix, with noise)')\n",
"ax.axvline(v_slw, color='C2', alpha=0.4, ls=':',\n",
" label=f'SLW peak ({v_slw:+.2f})')\n",
"ax.axvline(v_snow, color='C3', alpha=0.4, ls=':',\n",
" label=f'snow peak ({v_snow:+.2f})')\n",
"ax.axhline(dBZ(integ.noise_h / (integ.v[-1] - integ.v[0])),\n",
" color='k', ls='--', alpha=0.5, label='noise floor')\n",
"ax.set_xlabel('Doppler velocity v [m/s]')\n",
"ax.set_ylabel('sZ_h [dBZ / (m/s)]')\n",
"ax.set_title('W-band SLW vs. snow vs. combined Doppler spectrum')\n",
"ax.legend(fontsize=9, loc='upper right'); ax.grid(True, alpha=0.3)\n",
"fig.tight_layout();\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Takeaway\n",
"\n",
"The Billault-Roux et al. 2023 case study relies on the *spectral*\n",
"separability demonstrated here: a narrow near-zero Doppler peak is\n",
"SLW, a broader fall-velocity peak is ice; the polarimetric spectral\n",
"variables further distinguish columnar (high SLDR) from planar and\n",
"aggregated ice. With the two populations built as independent\n",
"rustmatrix scatterers you can plug them into `HydroMix` and recover\n",
"bulk radar observables that match what a vertically-pointing radar\n",
"would measure, or analyse the spectral signatures by looking at the\n",
"individual modes.\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python",
"version": "3.11"
}
},
"nbformat": 4,
"nbformat_minor": 5
}