{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Tutorial 06 \u2014 Hydrometeor mixtures at C-band\n",
"\n",
"Real radar volumes often contain more than one species \u2014 rain with\n",
"melting aggregates, graupel in convective cores, pristine crystals\n",
"with aggregates in stratiform ice. The combined polarimetric\n",
"signature is the incoherent sum of the per-species amplitude (S) and\n",
"phase (Z) matrices. The *non*-linear observables (Z_dr, \u03c1_hv, \u03b4_hv)\n",
"cannot be averaged from per-species values; they must be recomputed\n",
"from the summed matrices.\n",
"\n",
"This tutorial demonstrates how the **mixture \u03c1_hv drops below\n",
"unity** whenever two populations contribute comparable power and\n",
"carry different polarimetric fingerprints \u2014 even when each\n",
"individual species has \u03c1_hv \u2248 1. Three pairings illustrate the\n",
"point:\n",
"\n",
"1. **rain + light snow** (baseline) \u2014 rain dominates Z_h, \u03c1_hv barely\n",
" moves.\n",
"2. **rain + heavy snow** \u2014 snow concentration bumped so its Z_h\n",
" contribution matches rain's. The Z_dr mismatch (rain is positive,\n",
" snow is near zero with wide canting) drags \u03c1_hv down.\n",
"3. **rain + graupel** \u2014 near-spherical graupel with \u03c1 = 0.4 g/cm\u00b3,\n",
" D_max = 8 mm and wide Gaussian canting (\u03c3 = 40\u00b0). C-band\n",
" non-Rayleigh resonance at the large-D tail adds a small \u03b4_hv\n",
" signature and further degrades \u03c1_hv.\n",
"\n",
"Graupel parameters follow Ryzhkov, Zrni\u0107, Burgess (2005, *JAMC*\n",
"44:557) and Kumjian (2013, *J. Operational Meteor.*): density\n",
"0.4 g/cm\u00b3, oblateness 0.8, tumbling canting \u03c3 \u2248 40\u00b0.\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",
" orientation, radar, psd as rs_psd)\n",
"from rustmatrix.tmatrix_aux import (dsr_thurai_2007, geom_horiz_back,\n",
" geom_horiz_forw, K_w_sqr, wl_C)\n",
"from rustmatrix.refractive import m_w_10C, mi\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Build one scatterer per species\n",
"\n",
"Each species gets its own `PSDIntegrator` and (for snow/graupel) an\n",
"orientation PDF that models tumbling. All integrators register the\n",
"same two geometries the mixture will query (backscatter for Z_h,\n",
"Z_dr, \u03c1_hv; forward for K_dp, A_i) and share the wavelength.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"def build_rain():\n",
" s = Scatterer(wavelength=wl_C, m=m_w_10C[wl_C],\n",
" Kw_sqr=K_w_sqr[wl_C], ddelt=1e-4, ndgs=2)\n",
" integ = rs_psd.PSDIntegrator()\n",
" integ.D_max = 6.0; integ.num_points = 64\n",
" integ.axis_ratio_func = lambda D: 1.0 / dsr_thurai_2007(D)\n",
" integ.geometries = (geom_horiz_back, geom_horiz_forw)\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" return s\n",
"\n",
"def build_tumbling(rho, axis_ratio, D_max):\n",
" s = Scatterer(wavelength=wl_C, m=mi(wl_C, rho),\n",
" Kw_sqr=K_w_sqr[wl_C], axis_ratio=axis_ratio,\n",
" ddelt=1e-4, ndgs=2)\n",
" s.orient = orientation.orient_averaged_fixed\n",
" s.or_pdf = orientation.gaussian_pdf(std=40.0, mean=90.0)\n",
" s.n_alpha = 6; s.n_beta = 12\n",
" integ = rs_psd.PSDIntegrator()\n",
" integ.D_max = D_max; integ.num_points = 64\n",
" integ.geometries = (geom_horiz_back, geom_horiz_forw)\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" return s\n",
"\n",
"rain = build_rain()\n",
"snow = build_tumbling(rho=0.2, axis_ratio=0.6, D_max=8.0)\n",
"graupel = build_tumbling(rho=0.4, axis_ratio=0.8, D_max=8.0)\n",
"\n",
"rain_psd = rs_psd.GammaPSD(D0=1.5, Nw=8e3, mu=4, D_max=6.0)\n",
"snow_psd_light = rs_psd.ExponentialPSD(N0=5e3, Lambda=2.0, D_max=8.0)\n",
"snow_psd_heavy = rs_psd.ExponentialPSD(N0=1.5e5, Lambda=2.0, D_max=8.0)\n",
"graupel_psd = rs_psd.ExponentialPSD(N0=4e3, Lambda=1.4, D_max=8.0)\n",
"\n",
"rain.psd = rain_psd\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Assemble the three mixtures and read the observables\n",
"\n",
"Each mixture is a list of `MixtureComponent(Scatterer, PSD)` pairs.\n",
"`HydroMix` exposes a Scatterer-shaped API so the usual `radar.*`\n",
"helpers work on it directly.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"def obs(x):\n",
" x.set_geometry(geom_horiz_back)\n",
" out = dict(Zh=10*np.log10(radar.refl(x)),\n",
" Zdr=10*np.log10(radar.Zdr(x)),\n",
" rho=radar.rho_hv(x),\n",
" delta=np.degrees(radar.delta_hv(x)))\n",
" x.set_geometry(geom_horiz_forw)\n",
" out['Kdp'] = radar.Kdp(x); out['Ai'] = radar.Ai(x)\n",
" return out\n",
"\n",
"mix_lightsnow = HydroMix([\n",
" MixtureComponent(rain, rain_psd, 'rain'),\n",
" MixtureComponent(snow, snow_psd_light, 'snow'),\n",
"])\n",
"mix_heavysnow = HydroMix([\n",
" MixtureComponent(rain, rain_psd, 'rain'),\n",
" MixtureComponent(snow, snow_psd_heavy, 'snow'),\n",
"])\n",
"mix_graupel = HydroMix([\n",
" MixtureComponent(rain, rain_psd, 'rain'),\n",
" MixtureComponent(graupel, graupel_psd, 'graupel'),\n",
"])\n",
"\n",
"# Evaluate each case eagerly \u2014 the scatterer objects are shared by\n",
"# mutation (snow used by both light- and heavy-snow cases), so we read\n",
"# out observables immediately after setting each psd.\n",
"results = {}\n",
"results['rain only'] = obs(rain)\n",
"snow.psd = snow_psd_light; results['light snow only'] = obs(snow)\n",
"snow.psd = snow_psd_heavy; results['heavy snow only'] = obs(snow)\n",
"graupel.psd = graupel_psd; results['graupel only'] = obs(graupel)\n",
"results['rain + light snow'] = obs(mix_lightsnow)\n",
"results['rain + heavy snow'] = obs(mix_heavysnow)\n",
"results['rain + graupel'] = obs(mix_graupel)\n",
"\n",
"hdr = f'{\"case\":<22} {\"Z_h\":>7} {\"Z_dr\":>7} {\"rho_hv\":>9} {\"delta_hv\":>9} {\"K_dp\":>9}'\n",
"print(hdr)\n",
"print(f'{\"\":<22} {\"[dBZ]\":>7} {\"[dB]\":>7} {\"\":>9} {\"[deg]\":>9} {\"[deg/km]\":>9}')\n",
"print('-' * len(hdr))\n",
"for name, r in results.items():\n",
" print(f'{name:<22} {r[\"Zh\"]:>7.2f} {r[\"Zdr\"]:>+7.3f} '\n",
" f'{r[\"rho\"]:>9.5f} {r[\"delta\"]:>+9.4f} {r[\"Kdp\"]:>+9.4f}')\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bar-chart comparison\n",
"\n",
"The heavy-snow and graupel mixtures both push \u03c1_hv below the\n",
"rain-only and light-snow baselines. Z_dr drops because the\n",
"horizontally-oriented rain signal is diluted by near-isotropic\n",
"tumbling-ice contributions. K_dp rises in the heavy-snow mix because\n",
"the snow itself contributes a small forward phase shift.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"mix_names = ['rain only', 'rain + light snow',\n",
" 'rain + heavy snow', 'rain + graupel']\n",
"colors = ['C0', 'C2', 'C4', 'C1']\n",
"\n",
"fig, axes = plt.subplots(1, 4, figsize=(12, 3.4))\n",
"for ax, key, ylab in zip(\n",
" axes,\n",
" ['Zh', 'Zdr', 'rho', 'Kdp'],\n",
" ['Z_h [dBZ]', 'Z_dr [dB]', '\u03c1_hv', 'K_dp [\u00b0/km]'],\n",
"):\n",
" vals = [results[n][key] for n in mix_names]\n",
" ax.bar(range(len(mix_names)), vals, color=colors)\n",
" ax.set_xticks(range(len(mix_names)))\n",
" ax.set_xticklabels(mix_names, rotation=20, ha='right', fontsize=8)\n",
" ax.set_ylabel(ylab)\n",
" ax.grid(True, axis='y', alpha=0.3)\n",
"axes[2].set_ylim(min(results[n]['rho'] for n in mix_names) - 0.003, 1.0005)\n",
"fig.suptitle('C-band mixtures: heavy-snow and graupel pull \u03c1_hv below rain-only baseline')\n",
"fig.tight_layout();\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sweep the snow fraction\n",
"\n",
"Scale the snow PSD's N0 from 0 up to well beyond the rain-matching\n",
"value and watch Z_dr and \u03c1_hv interpolate between the rain-only\n",
"limit (left edge) and the snow-dominated limit (right edge). The\n",
"minimum \u03c1_hv sits near the crossover where Z_h_snow \u2248 Z_h_rain \u2014\n",
"exactly where the Z_dr mismatch has maximum leverage.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"N0_snow_sweep = np.geomspace(1e3, 5e5, 16)\n",
"snow.psd = snow_psd_light\n",
"Zh_sw, Zdr_sw, rho_sw, Kdp_sw = [], [], [], []\n",
"for N0 in N0_snow_sweep:\n",
" psd_i = rs_psd.ExponentialPSD(N0=float(N0), Lambda=2.0, D_max=8.0)\n",
" m_sw = HydroMix([\n",
" MixtureComponent(rain, rain_psd, 'rain'),\n",
" MixtureComponent(snow, psd_i, 'snow'),\n",
" ])\n",
" o = obs(m_sw)\n",
" Zh_sw.append(o['Zh']); Zdr_sw.append(o['Zdr'])\n",
" rho_sw.append(o['rho']); Kdp_sw.append(o['Kdp'])\n",
"\n",
"fig, axes = plt.subplots(2, 2, figsize=(9, 6), sharex=True)\n",
"axes[0, 0].semilogx(N0_snow_sweep, Zh_sw, 'C0-o'); axes[0, 0].set_ylabel('Z_h [dBZ]')\n",
"axes[0, 1].semilogx(N0_snow_sweep, Zdr_sw, 'C1-o'); axes[0, 1].set_ylabel('Z_dr [dB]')\n",
"axes[1, 0].semilogx(N0_snow_sweep, rho_sw, 'C2-o'); axes[1, 0].set_ylabel('\u03c1_hv')\n",
"axes[1, 1].semilogx(N0_snow_sweep, Kdp_sw, 'C3-o'); axes[1, 1].set_ylabel('K_dp [\u00b0/km]')\n",
"for ax in axes.flat:\n",
" ax.set_xlabel('snow N\u2080 [m\u207b\u00b3 mm\u207b\u00b9]')\n",
" ax.grid(True, alpha=0.3, which='both')\n",
"fig.suptitle('C-band rain + tumbling snow \u2014 sweeping snow N\u2080 across five decades')\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
}