Adding radar reflectivity to 1moment scheme

Co-authored-by: Caterina Croci <catecroci.cc@gmail.com>

docs/src/API.md

@@ -25,6 +25,7 @@ Microphysics1M
 Microphysics1M.get_v0
 Microphysics1M.get_n0
 Microphysics1M.lambda
+Microphysics1M.radar_reflectivity
 Microphysics1M.terminal_velocity
 Microphysics1M.conv_q_liq_to_q_rai
 Microphysics1M.conv_q_ice_to_q_sno_no_supersat

docs/src/Microphysics1M.md

@@ -263,6 +263,13 @@ In other derivations cloud ice, similar to cloud liquid water,
 
      - Do we want to test different size distributions?
 
+Here we plot the Marshall-Palmer particle size distribution for 4 different values for the rain specific humidity (q_rai).
+
+```@example
+include("plots/MarshallPalmer_distribution.jl")
+```
+![](MarshallPalmer_distribution.svg)
+
 ## Parameterized processes
 
 Parameterized processes include:
@@ -749,6 +756,34 @@ If ``T > T_{freeze}``:
     \right)
 \end{equation}
 ```
+
+## Rain radar reflectivity
+
+The rain radar reflectivity factor (``Z``) is used to measure the power returned by a radar signal when it encounters rain particles, and it is defined as the sixth moment of the rain particles distribution:
+```math
+\begin{equation}
+Z = {\int_0^\infty r^{6} \, n(r) \, dr}.
+\label{eq:Z}
+\end{equation}
+```
+Integrating over the assumed Marshall-Palmer distribution (eq. 6) leads to
+```math
+\begin{equation}
+Z = {\frac{6! \, n_{0}^{rai}}{\lambda^7}},
+\end{equation}
+```
+where:
+ - ``n_{0}^{rai}`` - rain drop size distribution parameter,
+ - ``\lambda`` - as defined in eq. 7
+
+By dividing ``Z`` with the equivalent return of a ``1 mm`` drop in a volume of a meter cube (``Z_0``) and applying the decimal logarithm to the result, we obtains the logarithmic rain radar reflectivity ``L_Z``, which is the variable that is commonly used to refer to the radar reflectivity values: 
+```math
+\begin{equation}
+L_Z = {10 \, \log_{10}(\frac{Z}{Z_0})}.
+\end{equation}
+```
+The resulting logarithmic dimensionless unit is decibel relative to ``Z``, or ``dBZ``.
+
 ## Example figures
 
 ```@example

docs/src/plots/MarshallPalmer_distribution.jl

@@ -0,0 +1,74 @@
+import CairoMakie as MK
+
+import SpecialFunctions as SF
+
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import ClimaParams as CP
+
+import Thermodynamics.Parameters as TDP
+import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.Common as CO
+import CloudMicrophysics.Microphysics1M as CM1
+
+FT = Float64
+
+function Marshall_Palmer_distribution(
+    (; pdf, mass)::CMP.Rain{FT},
+    q::FT,
+    ρ::FT,
+    r::FT,
+) where {FT}
+
+    n₀::FT = CM1.get_n0(pdf)
+    λ = CM1.lambda(pdf, mass, q, ρ)
+    # distribution
+    n_r = n₀ * exp(-λ * r)
+
+    return n_r
+end
+
+ρ = FT(1.01)
+qᵣ_0 = FT(1.0e-5)
+qᵣ_1 = FT(2.0e-4)
+qᵣ_2 = FT(1.0e-4)
+qᵣ_3 = FT(1.0e-3)
+
+r_range = range(25e-6, stop = 1e-2, length = 1000)
+
+n_r_0 =
+    [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_0, ρ, r) for r in r_range]
+n_r_1 =
+    [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_1, ρ, r) for r in r_range]
+n_r_2 =
+    [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_2, ρ, r) for r in r_range]
+n_r_3 =
+    [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_3, ρ, r) for r in r_range]
+
+fig = MK.Figure(resolution = (1100, 600))
+ax1 = MK.Axis(
+    fig[1, 1],
+    title = "Marshall-Palmer distribution",
+    xlabel = "r [m]",
+    ylabel = "n(r) [mm/m]",
+    limits = ((0, 0.0015), nothing),
+)
+MK.lines!(ax1, r_range, n_r_0, label = "q⁰ᵣ = 10⁻⁵", color = :blue)
+MK.lines!(ax1, r_range, n_r_1, label = "q¹ᵣ = 2x10⁻⁴", color = :red)
+MK.lines!(ax1, r_range, n_r_2, label = "q²ᵣ = 10⁻⁴", color = :green)
+MK.lines!(ax1, r_range, n_r_3, label = "q³ᵣ = 10⁻³", color = :orange)
+ax2 = MK.Axis(
+    fig[1, 2],
+    title = "Marshall-Palmer distribution - semilog scale",
+    ylabel = "n(r) [mm/m]",
+    xlabel = "r [m]",
+    xscale = log10,
+)
+MK.lines!(ax2, r_range, n_r_0, label = "q⁰ᵣ = 10⁻⁵", color = :blue)
+MK.lines!(ax2, r_range, n_r_1, label = "q¹ᵣ = 2x10⁻⁴", color = :red)
+MK.lines!(ax2, r_range, n_r_2, label = "q²ᵣ = 10⁻⁴", color = :green)
+MK.lines!(ax2, r_range, n_r_3, label = "q³ᵣ = 10⁻³", color = :orange)
+MK.axislegend()
+
+MK.save("MarshallPalmer_distribution.svg", fig)
+#! format: on

src/Microphysics1M.jl

@@ -97,6 +97,29 @@ function lambda(
     )^FT(1 / (me + Δm + 1)) : FT(0)
 end
 
+""" 
+    radar_reflectivity(precip, q, ρ)
+
+    - `precip` - struct with rain free parameters
+    - `q` - specific humidity of rain
+    - `ρ` - air density
+
+Returns radar reflectivity from the assumed rain particle size distribuion 
+normalized by the reflectivty of 1 millimiter drop in a volume of one meter cube
+"""
+function radar_reflectivity(
+    (; pdf, mass)::CMP.Rain{FT},
+    q::FT,
+    ρ::FT,
+) where {FT}
+
+    n0 = get_n0(pdf)
+    λ = lambda(pdf, mass, q, ρ)
+    Z₀ = FT(1e-18)
+
+    return 10 * log10((720 * n0 / λ^7) / Z₀)
+end
+
 """
     terminal_velocity(precip, vel, ρ, q)
 

test/microphysics1M_tests.jl

@@ -52,6 +52,22 @@ function test_microphysics1M(FT)
         end
     end
 
+    TT.@testset "1M_microphysics RadarReflectivity" begin
+
+        # some example values
+        ρ_air = FT(1)
+        q_rai = FT(0.18e-3)
+
+        TT.@test CM1.radar_reflectivity(rain, q_rai, ρ_air) ≈ FT(12.17) atol =
+            0.2
+
+        q_rai = FT(0.89e-4)
+
+        TT.@test CM1.radar_reflectivity(rain, q_rai, ρ_air) ≈ FT(6.68) atol =
+            0.2
+
+    end
+
     TT.@testset "1M_microphysics - Chen 2022 rain terminal velocity" begin
         #setup
         ρ = FT(1.2)

test/performance_tests.jl

@@ -173,6 +173,7 @@ function benchmark_test(FT)
         (liquid, rain, blk1mvel.rain, ce, q_liq, q_rai, ρ_air),
         350,
     )
+    bench_press(CM1.radar_reflectivity, (rain, q_rai, ρ_air), 250)
 
     # 2-moment
     bench_press(