P3 shedding!

docs/src/P3Scheme.md

@@ -366,6 +366,26 @@ include("plots/P3Melting.jl")
 ```
 ![](P3_ice_melt.svg)
 
+## Shedding
+
+Shedding is a sink of ``L_{liq}`` that transfers liquid mass on large, mixed-phase particles to rain.
+  Following [Choletteetal2019](@cite), shedding occurs only for ``D > 9 mm``, and it is scaled
+  with ``F_{rim}`` such that more in-filling with rime translates to more liquid mass being
+  shed from the ice core. It is calculated as follows by numerically integrating over the
+  whole mixed-phase particle PSD:
+```math
+\begin{equation}
+ \left. \frac{dL}{dt} \right|_{shed} =  \int_{9 mm}^{\infty} F_{rim} F_{liq} m(D) N_{0} D^{\mu} \exp^{- \lambda D} dD
+\end{equation}
+```
+The source of rain number concentration from shedding is computed as described by [Choletteetal2019](@cite)
+  assuming mean raindrop diameter ``D = 1 mm``:
+```math
+\begin{equation}
+  \left. \frac{dN_{rai}}{dt} \right|_{shed} = \frac{1}{\frac{\pi}{6}\rho_{l}D^{3}} \left. \frac{dL}{dt} \right|_{shed}
+\end{equation}
+```
+
 ## Acknowledgments
 
 Click on the P3 mascot duck to be taken to the repository

src/P3_particle_properties.jl

@@ -102,10 +102,10 @@ function thresholds(p3::PSP3{FT}, ρ_r::FT, F_rim::FT) where {FT}
     @assert F_rim >= FT(0)   # rime mass fraction must be positive ...
     @assert F_rim < FT(1)    # ... and there must always be some unrimed part
 
-    if F_rim == FT(0)
+    if F_rim == FT(0) || ρ_r == FT(0)
         return (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0))
     else
-        @assert ρ_r > FT(0)   # rime density must be positive ...
+        @assert ρ_r >= FT(0)   # rime density must be positive ...
         @assert ρ_r <= p3.ρ_l # ... and as a bulk ice density can't exceed the density of water
 
         P3_problem(ρ_d) =

src/P3_processes.jl

@@ -1,3 +1,24 @@
+"""
+A structure containing the rates of change of the specific humidities and number
+densities of liquid and rain water.
+"""
+Base.@kwdef struct P3Rates{FT}
+    "Rate of change of total ice mass content"
+    dLdt_p3_tot::FT = FT(0)
+    "Rate of change of rime mass content"
+    dLdt_rim::FT = FT(0)
+    "Rate of change of mass content of liquid on ice"
+    dLdt_liq::FT = FT(0)
+    "Rate of change of rime volume"
+    ddtB_rim::FT = FT(0)
+    "Rate of change of ice number concentration"
+    dNdt_ice::FT = FT(0)
+    "Rate of change of rain mass content"
+    dLdt_rai::FT = FT(0)
+    "Rate of change of rain number concentration"
+    dNdt_rai::FT = FT(0)
+end
+
 """
     het_ice_nucleation(pdf_c, p3, tps, q, N, T, ρₐ, p, aerosol)
 
@@ -114,3 +135,78 @@ function ice_melt(
     dLdt = min(dLdt, L_ice / dt)
     return (; dNdt, dLdt)
 end
+
+"""
+   ice_shed(p3, L_ice, N_ice, F_rim, ρ_rim, F_liq, dt)
+
+- p3 - a struct containing p3 parameters
+- L_p3_tot - total ice mass content
+- N_ice - ice number concentration
+- F_rim - rime mass fraction (L_rim / L_ice)
+- ρ_rim - rime density (L_rim/B_rim)
+- F_liq - liquid fraction (L_liq / L_p3_tot)
+- dt - model time step (for limiting the tendnecy)
+
+Returns the sink of L_liq due to shedding—N_ice remains constant.
+"""
+function ice_shed(
+    p3::PSP3,
+    L_p3_tot::FT,
+    N_ice::FT,
+    F_rim::FT,
+    ρ_rim::FT,
+    F_liq::FT,
+    dt::FT,
+) where {FT}
+    dLdt = FT(0)
+    if L_p3_tot > eps(FT) && N_ice > eps(FT) && F_liq > eps(FT)
+        # process dependent on F_liq
+        # (we want whole particle shape params)
+        # Get the P3 diameter distribution...
+        (λ, N_0) = distribution_parameter_solver(
+            p3,
+            L_p3_tot,
+            N_ice,
+            ρ_rim,
+            F_rim,
+            F_liq,
+        )
+        N(D) = N′ice(p3, D, λ, N_0)
+
+        # ... and D_max for the integral
+        # TODO - once get_ice_bound() is changed,
+        # find a reasonable tolerance for which
+        # it's worth computing shedding over the PSD
+        # (i.e. if only 1% of the particles are bigger than 9 mm
+        # is it really worth doing a computation and returning a tiny dLdt?)
+        bound = get_ice_bound(p3, λ, FT(1e-6))
+
+        # critical size for shedding
+        shed_bound = FT(9e-3)
+
+        # check if there are particles big enough to shed
+        shedding = ifelse(bound > shed_bound, true, false)
+
+        # liquid mass
+        m_liq(D) = F_liq * mass_s(D, p3.ρ_l)
+        # integrand (mass shed is a function of F_rim)
+        f(D) = F_rim * m_liq(D) * N(D)
+
+        # Integrate
+        if shedding
+            (dLdt, error) =
+                QGK.quadgk(d -> f(d), shed_bound, bound, rtol = FT(1e-6))
+        end
+
+        # ... don't exceed the available liquid mass content
+        dLdt = min(dLdt, F_liq * L_p3_tot / dt)
+    end
+    # return rates struct, with dNdt
+    # assuming raindrops of D = 1 mm
+    return P3Rates{FT}(
+        dLdt_p3_tot = dLdt,
+        dLdt_liq = dLdt,
+        dLdt_rai = dLdt,
+        dNdt_rai = dLdt / mass_s(FT(1e-3), p3.ρ_l),
+    )
+end

test/p3_tests.jl

@@ -23,8 +23,8 @@ function test_p3_thresholds(FT)
         F_rim_good = (FT(0.5), FT(0.8), FT(0.95)) # representative F_rim values
 
         # test asserts
-        for _ρ_r in (FT(0), FT(-1))
-            TT.@test_throws AssertionError("ρ_r > FT(0)") P3.thresholds(
+        for _ρ_r in (FT(-1))
+            TT.@test_throws AssertionError("ρ_r >= FT(0)") P3.thresholds(
                 p3,
                 _ρ_r,
                 F_rim,
@@ -868,6 +868,53 @@ function test_p3_melting(FT)
     end
 end
 
+function test_p3_shedding(FT)
+
+    TT.@testset "Shedding Smoke Test" begin
+
+        p3 = CMP.ParametersP3(FT)
+
+        qᵢ = FT(5e-3)
+        ρₐ = FT(1.2)
+        Lᵢ = qᵢ * ρₐ
+        Nᵢ = FT(5e3) * ρₐ
+        F_rim = FT(0.8)
+        ρ_rim = FT(800)
+        dt = FT(1)
+        F_liq = FT(0)
+
+        rate = P3.ice_shed(p3, Lᵢ, Nᵢ, F_rim, ρ_rim, F_liq, dt)
+
+        TT.@test rate.dLdt_p3_tot == 0
+        TT.@test rate.dLdt_liq == 0
+        TT.@test rate.dLdt_rai == 0
+        TT.@test rate.dNdt_rai == 0
+
+        F_liq = FT(0.9)
+
+        rate = P3.ice_shed(p3, Lᵢ, Nᵢ, F_rim, ρ_rim, F_liq, dt)
+
+        TT.@test rate.dLdt_p3_tot >= 0
+        TT.@test rate.dLdt_liq >= 0
+        TT.@test rate.dLdt_rai >= 0
+        TT.@test rate.dNdt_rai >= 0
+
+        TT.@test rate.dLdt_p3_tot == 3.929198115623993e-6
+        TT.@test rate.dNdt_rai == 7.504215629867026
+
+        Lᵢ = FT(0)
+        Nᵢ = FT(0)
+
+        rate = P3.ice_shed(p3, Lᵢ, Nᵢ, F_rim, ρ_rim, F_liq, dt)
+
+        TT.@test rate.dLdt_p3_tot == 0
+        TT.@test rate.dLdt_liq == 0
+        TT.@test rate.dLdt_rai == 0
+        TT.@test rate.dNdt_rai == 0
+
+    end
+end
+
 
 println("Testing Float32")
 test_p3_thresholds(Float32)
@@ -886,4 +933,5 @@ test_bulk_terminal_velocities(Float64)
 test_integrals(Float64)
 test_p3_het_freezing(Float64)
 test_p3_melting(Float64)
+test_p3_shedding(Float64)
 #test_tendencies(Float64)

test/performance_tests.jl

@@ -199,6 +199,13 @@ function benchmark_test(FT)
             2e3,
             3,
         )
+        bench_press(
+            P3.ice_shed,
+            (p3, L, N, F_rim, ρ_r, F_liq, Δt),
+            1.5e5,
+            1.5e3,
+            2,
+        )
     end
 
     # aerosol activation