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