Rain radar observables

This notebook turns single-particle T-matrices into simple rain radar and brightness-temperature observables. It is adapted from the radiative-transfer prototype contributed by @xiongyuup in PR #3.

The original contribution demonstrates a larger workflow with WRF input, hydrometeor microphysics, radar moments, and microwave radiative transfer. This example keeps only the self-contained core: water refractive index, rain-drop axis ratio, spheroidal T-matrices, a Marshall-Palmer-type drop-size distribution, and integrated dual-polarization moments.

begin
    import Pkg
    Pkg.activate(@__DIR__)
    Pkg.develop(; path = dirname(@__DIR__))
    Pkg.instantiate()
    using TransitionMatrices, Plots
end

Water, ice, and drop models

The dielectric fits below are compact microwave models for liquid water and ice. Rain drops are represented as oblate spheroids with diameter \(D\) and axis ratio \(c/a\).

begin
    c0_m_per_s() = 2.99792458e8
    c0_mm_per_s() = 1000 * c0_m_per_s()

    function water_refractive_index(freq_hz, T_celsius)
        λ_m = c0_m_per_s() / freq_hz
        ε0 = 8.854e-12
        εs = 78.54 * (1 - 4.579e-3 * (T_celsius - 25) +
                      1.19e-5 * (T_celsius - 25)^2 -
                      2.8e-8 * (T_celsius - 25)^3)
        ε∞ = 5.27137 + 2.16474e-2 * T_celsius - 1.31198e-3 * T_celsius^2
        α = -16.8129 / (T_celsius + 273) + 6.09265e-2
        λs = 3.3836e-6 * exp(2513.98 / (T_celsius + 273))
        σ = 1.1117e-4

        x = (λs / λ_m)^(1 - α)
        denominator = 1 + 2x * sinpi(α / 2) + (λs / λ_m)^(2 - 2α)
        εreal = ε∞ + (εs - ε∞) * (1 + x * sinpi(α / 2)) / denominator
        εimag = (εs - ε∞) * x * cospi(α / 2) / denominator +
                σ * λ_m / (2π * c0_m_per_s() * ε0)

        sqrt(complex(εreal, εimag))
    end

    function ice_refractive_index(freq_hz, T_celsius)
        λ_m = c0_m_per_s() / freq_hz
        ε0 = 8.854e-12
        εs = 203.168 + 2.5 * T_celsius + 0.15T_celsius^2
        ε∞ = 3.168
        α = 0.288 + 5.2e-3 * T_celsius + 2.3e-4 * T_celsius^2
        λs = 9.990288e-6 * exp(6643.5 / (T_celsius + 273))
        σ = 1.1146e-13 * exp(-6291.2 / (T_celsius + 273))

        x = (λs / λ_m)^(1 - α)
        denominator = 1 + 2x * sinpi(α / 2) + (λs / λ_m)^(2 - 2α)
        εreal = ε∞ + (εs - ε∞) * (1 + x * sinpi(α / 2)) / denominator
        εimag = (εs - ε∞) * x * cospi(α / 2) / denominator +
                σ * λ_m / (2π * c0_m_per_s() * ε0)

        sqrt(complex(εreal, εimag))
    end

    function rain_axis_ratio(D_mm)
        D_mm < 0.7 && return 1.0
        D_mm < 1.5 && return 1.173 - 0.5265D_mm + 0.4698D_mm^2 -
                            0.1317D_mm^3 - 8.5e-3D_mm^4
        1.065 - 6.25e-2D_mm - 3.99e-3D_mm^2 +
            7.66e-4D_mm^3 - 4.095e-5D_mm^4
    end
end

rain_axis_ratio (generic function with 1 method)

begin
    round_complex(z; digits = 4) =
        complex(round(real(z); digits), round(imag(z); digits))

    material_summary = [
        (; material = "water", freq_GHz = 9.375, T = 0.0,
            m = round_complex(water_refractive_index(9.375e9, 0.0))),
        (; material = "ice", freq_GHz = 9.375, T = -10.0,
            m = round_complex(ice_refractive_index(9.375e9, -10.0)))
    ]
end

2-element Vector{@NamedTuple{material::String, freq_GHz::Float64, T::Float64, m::ComplexF64}}:
 (material = "water", freq_GHz = 9.375, T = 0.0, m = 7.252 + 2.8558im)
 (material = "ice", freq_GHz = 9.375, T = -10.0, m = 1.7799 + 0.0001im)

Single-particle scattering table

Each diameter is solved as a spheroid with the IITM backend. The backward amplitudes feed reflectivity and differential reflectivity; the forward amplitudes feed differential phase; extinction feeds the brightness-temperature toy calculation.

begin
    function scattering_table(Ds, freq_hz; temperature_celsius = 0.0,
            solver = IITM(6, 10, 16))
        λ = c0_mm_per_s() / freq_hz
        m = water_refractive_index(freq_hz, temperature_celsius)
        map(Ds) do D
            a = D / 2
            axis_ratio = rain_axis_ratio(D)
            c = a * axis_ratio
            shape = Spheroid{Float64, ComplexF64}(a, c, ComplexF64(m))
            T = transition_matrix(shape, λ, solver)
            Sback = amplitude_matrix(T, π / 2, 0.0, π / 2, π; λ)
            Sfwd = amplitude_matrix(T, π / 2, 0.0, π / 2, 0.0; λ)
            (; D_mm = D,
                axis_ratio,
                Shh_back = Sback[2, 2],
                Svv_back = Sback[1, 1],
                Shh_forward = Sfwd[2, 2],
                Svv_forward = Sfwd[1, 1],
                Cext_mm2 = calc_Cext(T, λ))
        end
    end

    spacing(table) = length(table) > 1 ? table[2].D_mm - table[1].D_mm : 1.0
end

spacing (generic function with 1 method)

begin
    freq_hz = 9.375e9
    temperature_celsius = 0.0
    Ds = collect(0.75:0.75:3.75)
    solver = IITM(6, 10, 16)

    table = scattering_table(Ds, freq_hz; temperature_celsius, solver)
    table_summary = [
        (; D = row.D_mm,
            axis_ratio = round(row.axis_ratio; digits = 3),
            Cext = round(row.Cext_mm2; digits = 4))
        for row in table
    ]
end

5-element Vector{@NamedTuple{D::Float64, axis_ratio::Float64, Cext::Float64}}:
 (D = 0.75, axis_ratio = 0.984, Cext = 0.0052)
 (D = 1.5, axis_ratio = 0.965, Cext = 0.0705)
 (D = 2.25, axis_ratio = 0.912, Cext = 0.4337)
 (D = 3.0, axis_ratio = 0.859, Cext = 1.8037)
 (D = 3.75, axis_ratio = 0.807, Cext = 5.3644)

table_summary

5-element Vector{@NamedTuple{D::Float64, axis_ratio::Float64, Cext::Float64}}:
 (D = 0.75, axis_ratio = 0.984, Cext = 0.0052)
 (D = 1.5, axis_ratio = 0.965, Cext = 0.0705)
 (D = 2.25, axis_ratio = 0.912, Cext = 0.4337)
 (D = 3.0, axis_ratio = 0.859, Cext = 1.8037)
 (D = 3.75, axis_ratio = 0.807, Cext = 5.3644)

Bulk rain observables

For compactness this uses \(N(D, R) = N_0 \exp[-\Lambda(R)D]\) with \(D\) in millimeters. The constants are chosen to produce a plausible monotonic example, not a complete retrieval algorithm.

begin
    drop_distribution(D_mm, rain_rate_mm_h) =
        rain_rate_mm_h == 0 ? 0.0 : 8.0e3 * exp(-4.1 * rain_rate_mm_h^(-0.21) * D_mm)

    function dualpol_moments(table, freq_hz, rain_rate_mm_h;
            temperature_celsius = 0.0)
        λ = c0_mm_per_s() / freq_hz
        ΔD = spacing(table)
        m = water_refractive_index(freq_hz, temperature_celsius)
        Kw2 = abs2((m^2 - 1) / (m^2 + 2))
        weights = [
            drop_distribution(row.D_mm, rain_rate_mm_h) * ΔD
            for row in table
        ]

        hh = sum(abs2(row.Shh_back) * weights[i] for (i, row) in pairs(table))
        vv = sum(abs2(row.Svv_back) * weights[i] for (i, row) in pairs(table))
        Zhh = hh * λ^4 / (π^5 * Kw2)
        Zvv = vv * λ^4 / (π^5 * Kw2)
        cov = sum(conj(row.Shh_back) * row.Svv_back * weights[i]
            for (i, row) in pairs(table))
        ρhv = abs(cov) / sqrt(hh * vv)
        kdp = 180 / π * 1e-3 * λ *
              sum(real(row.Shh_forward - row.Svv_forward) * weights[i]
                  for (i, row) in pairs(table))

        (; R = rain_rate_mm_h,
            ZH = round(10log10(Zhh); digits = 2),
            ZDR = round(10log10(Zhh / Zvv); digits = 2),
            ρhv = round(ρhv; digits = 4),
            KDP = round(kdp; digits = 4))
    end

    function brightness_temperature(table, rain_rate_mm_h;
            path_km = 5.0, Tatm_K = 270.0)
        ΔD = spacing(table)
        α = sum(drop_distribution(row.D_mm, rain_rate_mm_h) *
                row.Cext_mm2 * ΔD for row in table) / 1e6
        τ = α * path_km
        Tatm_K * (1 - exp(-τ))
    end
end

brightness_temperature (generic function with 1 method)

begin
    rain_rates = [1.0, 10.0, 50.0]
    radar = [
        dualpol_moments(table, freq_hz, R; temperature_celsius)
        for R in rain_rates
    ]
    tb = [
        (; R, TB = round(brightness_temperature(table, R); digits = 4))
        for R in rain_rates
    ]
    (; radar, tb)
end

(radar = [(R = 1.0, ZH = 13.23, ZDR = 0.56, ρhv = 0.9987, KDP = 0.0412), (R = 10.0, ZH = 27.33, ZDR = 1.23, ρhv = 0.9963, KDP = 0.5779), (R = 50.0, ZH = 36.03, ZDR = 1.66, ρhv = 0.9962, KDP = 3.1766)], tb = [(R = 1.0, TB = 0.0036), (R = 10.0, TB = 0.0418), (R = 50.0, TB = 0.2254)])

Multi-frequency brightness-temperature lookup

The larger PR #3 prototype writes lookup tables to MAT files. Here the same idea stays in memory: build one scattering table per frequency, integrate over the rain distribution, and return a compact table that Pluto renders directly.

begin
    function brightness_temperature_lut(Ds, frequencies_hz, rain_rates;
            temperature_celsius = 0.0, solver = IITM(6, 10, 16),
            path_km = 5.0, Tatm_K = 270.0)
        rows = NamedTuple[]
        for f in frequencies_hz
            local_table = scattering_table(Ds, f; temperature_celsius, solver)
            for R in rain_rates
                push!(rows, (; freq_GHz = round(f / 1e9; digits = 3), R,
                    TB = round(brightness_temperature(local_table, R;
                        path_km, Tatm_K); digits = 4)))
            end
        end
        rows
    end

    lookup_frequencies = [10e9, 18e9]
    tb_lut = brightness_temperature_lut(Ds, lookup_frequencies, rain_rates;
        temperature_celsius, solver)
end

6-element Vector{NamedTuple}:
 (freq_GHz = 10.0, R = 1.0, TB = 0.0042)
 (freq_GHz = 10.0, R = 10.0, TB = 0.0496)
 (freq_GHz = 10.0, R = 50.0, TB = 0.2675)
 (freq_GHz = 18.0, R = 1.0, TB = 0.0167)
 (freq_GHz = 18.0, R = 10.0, TB = 0.1989)
 (freq_GHz = 18.0, R = 50.0, TB = 0.9967)

tb_lut

6-element Vector{NamedTuple}:
 (freq_GHz = 10.0, R = 1.0, TB = 0.0042)
 (freq_GHz = 10.0, R = 10.0, TB = 0.0496)
 (freq_GHz = 10.0, R = 50.0, TB = 0.2675)
 (freq_GHz = 18.0, R = 1.0, TB = 0.0167)
 (freq_GHz = 18.0, R = 10.0, TB = 0.1989)
 (freq_GHz = 18.0, R = 50.0, TB = 0.9967)

Plots

let
    p1 = plot(getfield.(table_summary, :D), getfield.(table_summary, :Cext);
        marker = :circle, lw = 2, label = "Cext",
        xlabel = "drop diameter D (mm)", ylabel = "extinction cross section (mm²)")
    p2 = plot(getfield.(radar, :R), getfield.(radar, :ZH);
        marker = :circle, lw = 2, xscale = :log10, label = "ZH",
        xlabel = "rain rate (mm h⁻¹)", ylabel = "reflectivity (dBZ)")
    p3 = plot(getfield.(radar, :R), getfield.(radar, :ZDR);
        marker = :circle, lw = 2, xscale = :log10, label = "ZDR",
        xlabel = "rain rate (mm h⁻¹)", ylabel = "differential reflectivity (dB)")
    p4 = plot(getfield.(tb, :R), getfield.(tb, :TB);
        marker = :circle, lw = 2, xscale = :log10, label = "TB",
        xlabel = "rain rate (mm h⁻¹)", ylabel = "brightness temperature (K)")
    plot(p1, p2, p3, p4; layout = (2, 2), size = (820, 620))
end

let
    p = plot(; xscale = :log10, xlabel = "rain rate (mm h⁻¹)",
        ylabel = "brightness temperature (K)", legend = :topleft)
    for f in getfield.(tb_lut, :freq_GHz) |> unique
        subset = filter(row -> row.freq_GHz == f, tb_lut)
        plot!(p, getfield.(subset, :R), getfield.(subset, :TB);
            marker = :circle, lw = 2, label = "$(f) GHz")
    end
    p
end

References

H. Li, Y. Xiong, and Y. Chen, "Simulation of Complex Meteorological Target Echoes for Airborne Dual-Polarization Weather Radar Based on Invariant Imbedding T-Matrix," IEEE Transactions on Geoscience and Remote Sensing, 62, 5105817, 2024.
G. Zhang, Weather Radar Polarimetry, China Meteorological Press, 2018, pp. 39-40.
B. R. Brown, M. M. Bell, and A. J. Frambach, "Validation of Simulated Hurricane Drop Size Distributions Using Polarimetric Radar," Geophysical Research Letters, 42, 2016.