Aerosol activation in parcel

NEWS.md

@@ -5,6 +5,8 @@ main
 ------
 <!--- # Add changes since the most recent release here --->
 
+- Add aerosol activation to parcel. ([#429](https://github.com/CliMA/CloudMicrophysics.jl/pull/429))
+
 - Bugfix; fixed the evaporation scheme of SB2006 to prevent returning NaN values when limiters are not applied. ([#420](https://github.com/CliMA/CloudMicrophysics.jl/pull/415))
 
 v0.20.0

docs/src/IceNucleationParcel0D.md

@@ -164,6 +164,18 @@ A similar approach is used for the lognormal size distribution. We assume that t
 \end{equation}
 ```
 
+## Aerosol Activation
+Aerosol activation is described by ([see discussion](https://clima.github.io/CloudMicrophysics.jl/dev/AerosolActivation/#Number-and-mass-of-activated-particles)).
+  It is inherently assumed that the aerosols have a lognormal size distribution.
+  For simplicity, the parcel accepts one mode and one aerosol type at a time, therefore,
+  internal mixing is not needed. The maxiumum supersaturation as described in the above
+  mentioned documentation is replaced by the liquid supersaturation in the parcel as
+  it evolves over time.
+
+!!! note
+    Standard deviation and r_{mean} of the aerosol size distribution may change as aerosols activate.
+    For now, we will neglect these effects.
+
 ## Condensation growth
 
 The diffusional growth of individual cloud droplet is described by
@@ -288,7 +300,15 @@ Here we show various example simulation results from the adiabatic parcel
   liquid processes only, immersion freezing with condensation and deposition growth,
   and homogeneous freezing with deposition growth.
 
-We start with deposition freezing on dust.
+First, we check that aerosol activation works reasonably within the parcel.
+
+```@example
+include("../../parcel/Example_AerosolActivation.jl")
+```
+![](Parcel_Aerosol_Activation.svg)
+
+The following examples show ice nucleation, starting with deposition
+  freezing on dust.
 The model is run three times using the `"MohlerAF_Deposition"` approach
   for 30 minutes simulation time, (shown by three different colors on the plot).
 Between each run the water vapor specific humidity is changed,

parcel/Example_AerosolActivation.jl

@@ -0,0 +1,83 @@
+import OrdinaryDiffEq as ODE
+import CairoMakie as MK
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import ClimaParams as CP
+
+# definition of the ODE problem for parcel model
+include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
+FT = Float32
+
+# Get free parameters
+tps = TD.Parameters.ThermodynamicsParameters(FT)
+wps = CMP.WaterProperties(FT)
+
+# Initial conditions
+Nₐ = FT(5e8)
+Nₗ = FT(0)
+Nᵢ = FT(0)
+T₀ = FT(230)
+cᵥ₀ = FT(5 * 1e-5)
+ln_INPC = FT(0)
+
+# Constants
+ρₗ = wps.ρw
+R_v = TD.Parameters.R_v(tps)
+R_d = TD.Parameters.R_d(tps)
+ϵₘ = R_d / R_v
+eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
+qᵥ = ϵₘ / (ϵₘ - 1 + 1 / cᵥ₀)
+# Compute qₗ assuming that initially droplets are lognormally distributed
+# with N(r₀, σ). We are not keeping that size distribution assumption
+# in the simulation
+r₀ = FT(3e-7)
+σ = FT(2)
+qₗ = Nₗ * FT(4 / 3 * π) * exp((6 * log(r₀) + 9 * σ^2) / (2))
+qᵢ = FT(0)
+Sₗ = FT(0.99)
+e = Sₗ * eₛ
+p₀ = e / cᵥ₀
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, ln_INPC]
+
+# Simulation parameters passed into ODE solver
+w = FT(1.2)         # updraft speed
+const_dt = FT(1)    # model timestep
+t_max = FT(100)     # total time
+aerosol = CMP.Sulfate(FT)
+
+condensation_growth = "Condensation"
+aerosol_act = "AeroAct"     # turn on aerosol activation
+aero_σ_g = FT(2.3)
+r_nuc = r₀
+
+params = parcel_params{FT}(
+    w = w,
+    const_dt = const_dt,
+    aerosol_act = aerosol_act,
+    aerosol = aerosol,
+    aero_σ_g = aero_σ_g,
+    r_nuc = r_nuc,
+    condensation_growth = condensation_growth,
+)
+
+# solve ODE
+sol = run_parcel(IC, FT(0), t_max, params)
+
+# Plot results
+fig = MK.Figure(size = (1000, 800), fontsize = 20)
+ax1 = MK.Axis(fig[1, 1], ylabel = "Liquid Saturation [-]")
+ax2 = MK.Axis(fig[1, 2], xlabel = "Time [s]", ylabel = "Temperature [K]")
+ax3 = MK.Axis(fig[2, 1], ylabel = "N_aero [m^-3]", xlabel = "Time [s]")
+ax4 = MK.Axis(fig[2, 2], ylabel = "N_liq [m^-3]", xlabel = "Time [s]")
+ax5 = MK.Axis(fig[3, 1], ylabel = "q_vap [g/kg]", xlabel = "Time [s]")
+ax6 = MK.Axis(fig[3, 2], ylabel = "q_liq [g/kg]", xlabel = "Time [s]")
+
+MK.lines!(ax1, sol.t, (sol[1, :]))            # liq saturation
+MK.lines!(ax2, sol.t, sol[3, :])              # temperature
+MK.lines!(ax3, sol.t, sol[7, :])              # N_aero
+MK.lines!(ax4, sol.t, sol[8, :])              # N_liq
+MK.lines!(ax5, sol.t, sol[4, :] .* 1e3)       # q_vap
+MK.lines!(ax6, sol.t, sol[5, :] .* 1e3)       # q_liq
+
+
+MK.save("Parcel_Aerosol_Activation.svg", fig)

parcel/ParcelModel.jl

@@ -1,12 +1,15 @@
 import Thermodynamics as TD
 import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.AerosolActivation as AA
+import CloudMicrophysics.AerosolModel as AM
 
 include(joinpath(pkgdir(CM), "parcel", "ParcelParameters.jl"))
 
 """
     Parcel simulation parameters
 """
 Base.@kwdef struct parcel_params{FT} <: CMP.ParametersType{FT}
+    aerosol_act = "None"
     deposition = "None"
     heterogeneous = "None"
     homogeneous = "None"
@@ -15,9 +18,11 @@ Base.@kwdef struct parcel_params{FT} <: CMP.ParametersType{FT}
     liq_size_distribution = "Monodisperse"
     ice_size_distribution = "Monodisperse"
     aerosol = Empty{FT}()
+    aero_σ_g = FT(0)
     wps = CMP.WaterProperties(FT)
     aps = CMP.AirProperties(FT)
     tps = TD.Parameters.ThermodynamicsParameters(FT)
+    aap = CMP.AerosolActivationParameters(FT)
     ips = CMP.IceNucleationParameters(FT)
     H₂SO₄ps = CMP.H2SO4SolutionParameters(FT)
     const_dt = 1
@@ -44,7 +49,8 @@ function parcel_model(dY, Y, p, t)
     # Simulation parameters
     (; wps, tps, r_nuc, w) = p
     (; liq_distr, ice_distr) = p
-    (; dep_params, imm_params, hom_params, ce_params, ds_params) = p
+    (; aero_act_params, dep_params, imm_params, hom_params) = p
+    (; ce_params, ds_params) = p
     # Y values stored in a named tuple for ease of use
     state = (
         Sₗ = Y[1],
@@ -90,6 +96,11 @@ function parcel_model(dY, Y, p, t)
     # Mean radius, area and volume of liquid droplets and ice crystals
     PSD_liq = distribution_moments(liq_distr, qₗ, Nₗ, ρₗ, ρ_air)
     PSD_ice = distribution_moments(ice_distr, qᵢ, Nᵢ, ρᵢ, ρ_air)
+
+    # Aerosol activation
+    dNₗ_dt_act = aerosol_activation(aero_act_params, state)
+    dqₗ_dt_act = dNₗ_dt_act * 4 * FT(π) / 3 * r_nuc^3 * ρₗ / ρ_air
+
     # Deposition ice nucleation
     # (All deposition parameterizations assume monodisperse aerosol size distr)
     dNᵢ_dt_dep = deposition_nucleation(dep_params, state, dY)
@@ -112,17 +123,17 @@ function parcel_model(dY, Y, p, t)
 
     # number concentration and ...
     dNᵢ_dt = dNᵢ_dt_dep + dNᵢ_dt_imm + dNᵢ_dt_hom
-    dNₐ_dt = -dNᵢ_dt_dep
-    dNₗ_dt = -dNᵢ_dt_imm - dNᵢ_dt_hom
+    dNₐ_dt = -dNᵢ_dt_dep - dNₗ_dt_act
+    dNₗ_dt = dNₗ_dt_act - dNᵢ_dt_imm - dNᵢ_dt_hom
     # ... water mass budget
-    dqₗ_dt_v2l = dqₗ_dt_ce
+    dqₗ_dt_v2l = dqₗ_dt_ce + dqₗ_dt_act
     dqᵢ_dt_l2i = dqᵢ_dt_imm + dqᵢ_dt_hom
     dqᵢ_dt_v2i = dqᵢ_dt_dep + dqᵢ_dt_ds
 
     # Update the tendecies
     dqᵢ_dt = dqᵢ_dt_v2i + dqᵢ_dt_l2i
     dqₗ_dt = dqₗ_dt_v2l - dqᵢ_dt_l2i
-    dqᵥ_dt = -dqᵢ_dt - dqₗ_dt
+    dqᵥ_dt = -dqₗ_dt_v2l - dqᵢ_dt_v2i
 
     dSₗ_dt =
         a1 * w * Sₗ - (a2 + a3) * Sₗ * dqₗ_dt_v2l -
@@ -216,6 +227,16 @@ function run_parcel(IC, t_0, t_end, pp)
 
     info *= "Aerosol: $(chop( string(typeof(pp.aerosol)), head = 29, tail = 9))\n"
 
+    info *= "Aerosol Activation: $(pp.aerosol_act)\n"
+    if pp.aerosol_act == "None"
+        aero_act_params = Empty{FT}()
+    elseif pp.aerosol_act == "AeroAct"
+        aero_act_params =
+            AeroAct{FT}(pp.aap, pp.aerosol, pp.aero_σ_g, pp.r_nuc, pp.const_dt)
+    else
+        throw("Unrecognized aerosol activation mode")
+    end
+
     info *= "Deposition: $(pp.deposition)\n"
     if pp.deposition == "None"
         dep_params = Empty{FT}()
@@ -283,13 +304,16 @@ function run_parcel(IC, t_0, t_end, pp)
     p = (
         liq_distr = liq_distr,
         ice_distr = ice_distr,
+        aero_act_params,
         dep_params = dep_params,
         imm_params = imm_params,
         hom_params = hom_params,
         ce_params = ce_params,
         ds_params = ds_params,
         wps = pp.wps,
         tps = pp.tps,
+        aps = pp.aps,
+        aap = pp.aap,
         r_nuc = pp.r_nuc,
         w = pp.w,
     )

parcel/ParcelParameters.jl

@@ -3,6 +3,14 @@ import Thermodynamics.Parameters as TDP
 
 struct Empty{FT} <: CMP.ParametersType{FT} end
 
+struct AeroAct{FT} <: CMP.ParametersType{FT}
+    aap::CMP.ParametersType{FT}
+    aerosol::CMP.AerosolType{FT}
+    aero_σ_g::FT
+    r_nuc::FT
+    const_dt::FT
+end
+
 struct MohlerAF{FT} <: CMP.ParametersType{FT}
     ips::CMP.ParametersType{FT}
     aerosol::CMP.AerosolType{FT}

parcel/ParcelTendencies.jl

@@ -4,6 +4,35 @@ import CloudMicrophysics.HetIceNucleation as CMI_het
 import CloudMicrophysics.HomIceNucleation as CMI_hom
 import CloudMicrophysics.Parameters as CMP
 import Distributions as DS
+import SpecialFunctions as SF
+
+function aerosol_activation(::Empty, state)
+    FT = eltype(state)
+    return FT(0)
+end
+
+function aerosol_activation(params::AeroAct, state)
+    (; aap, aerosol, aero_σ_g, r_nuc, const_dt) = params
+    (; T, Sₗ, Nₐ) = state
+    FT = eltype(state)
+
+    ad = AM.Mode_κ(
+        r_nuc,
+        aero_σ_g,
+        Nₐ,
+        (FT(1.0),),
+        (FT(1.0),),
+        (aerosol.M,),
+        (aerosol.κ,),
+    )
+    all_ad = AM.AerosolDistribution((ad,))
+
+    smax = (Sₗ - 1) < 0 ? FT(0) : (Sₗ - 1)
+    sm = AA.critical_supersaturation(aap, all_ad, T)
+    u = 2 * log(sm[1] / smax) / 3 / sqrt(2) / log(ad.stdev)
+
+    return ad.N * (1 / 2) * (1 - SF.erf(u)) / const_dt
+end
 
 function deposition_nucleation(::Empty, state, dY)
     FT = eltype(state)

parcel/Project.toml

@@ -5,5 +5,6 @@ CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"