isorropia progress

Project.toml

@@ -11,6 +11,7 @@ EarthSciMLBase = "e53f1632-a13c-4728-9402-0c66d48804b0"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"

docs/src/isorropia/implementation.md

@@ -38,6 +38,37 @@ to be effectively instantaneous.
 However, to make the equilibrium occur more quickly or slowly we could simply multiply both forward and reverse reaction 
 rate constants by a constant factor.
 
+## Solid phase activity
+
+The paper states that the "activity of each solid phase species is assumed to be unity." (Fountoukis and Nenes, 2007)
+However, this cannot be strictly true because if the activity were always equal to one the solid would continue to be consumed 
+even after its concentration had reached zero, resulting in negative concentrations which are physically impossible.
+Instead, the activity of a solid should be a function which asymptotes to one for normal concentrations and 
+asymptotes to infinity for concentrations approaching zero, with an inflection point at a concentration where the solid
+would not considered to be meaningfully present in the aerosol.
+Given that an atmospherically relevant aerosol concentration is typically ≥1 ``\textrm{\mu g m^{-3}}``, we choose an inflection
+concentration of 0.01 ``\mu g m^{-3}``, which is around 1e-10 mol ``\textrm{m^{-3}}`` for sulfate.
+Smooth function forms are preferable for numerical stability so we therefore choose the equation:
+
+```math
+a(c) = 10^{-(log10(c)+10)}+1
+```
+
+where ``a`` is activity, ``c`` is concentration, and ``+10`` represents the exponent of our threshold value. The resulting curve looks like this:
+
+```@example
+using Plots
+c = 10 .^(-13:0.1:-7)
+a(c) = 10^-(log10(c)+10)+1
+plot(c, a.(c), yscale=:log10, xscale=:log10, title="Solid Activity Coefficient",
+    xlabel="Concentration (μg/m³)", ylabel="Activity Coefficient", label=:none)
+
+
+logit(x) = 1 / (1 + exp(-(log10(x)+11)*5))
+x = 10 .^ (-13:0.1:-10)
+plot(x, logit.(x), xscale=:log10)
+```
+
 ## Deliquescence
 
 When interpolation between MRDH and DRH relative humidities as in Foungtoukis and Nenes (2007), we take ``\mathrm{RH_{wet}}``

src/isorropia/isorropia.jl

@@ -1,7 +1,7 @@
 #module ISORROPIA
 using EarthSciMLBase
 using ModelingToolkit, Catalyst, Unitful
-using IfElse
+using IfElse, NaNMath
 
 using Latexify
 using NonlinearSolve
@@ -107,16 +107,17 @@ struct Isorropia <: EarthSciMLODESystem
         active_salts = intersect(active_specs, values(salts))
         active_solids = intersect(active_specs, values(solids))
         active_gases = intersect(active_specs, values(gases))
+        active_ions_including_salts = unique(vcat(active_ions, [[s.cation, s.anion] for s ∈ active_salts]...))
 
         water = Water(t, active_salts)
-        ionic_strength = IonicStrength(t, active_ions, water.W)
+        ionic_strength = IonicStrength(t, active_ions_including_salts, water.W)
         salt_act = Activities(t, active_salts, (s) -> activity(s, active_salts, salts, ionic_strength.I, water.W))
         ion_act = Activities(t, active_ions, (i) -> activity(i, water.W))
         gas_act = Activities(t, active_gases, activity)
         solid_act = Activities(t, active_solids, activity)
         water_act = Activities(t, [H2Oaq], activity)
         activities = extend(water_act, extend(salt_act, extend(ion_act, extend(gas_act, solid_act))))
-        balance = Balance(t, T, active_specs, active_ions, active_salts, gases, solids, ions)
+        balance = Balance(t, T, active_specs, active_ions_including_salts, gases, solids, ions)
         #out = Output(active_vars)
 
         # Convert dictionaries of equations into an equation system by combining equations with the same left-hand side.
@@ -138,18 +139,25 @@ struct Isorropia <: EarthSciMLODESystem
         sys = ODESystem(eqs, t, active_vars, [T, RH]; systems=syss, name=name)
         sys = extend(sys, extend(ionic_strength, extend(activities, extend(balance, water))))
         #sys = sub_trans(sys, active_vars)
-        new(sys, mw, molecs, solids, gases, ions, salts)
+
+        filter_present(d, l) = filter(((k,v),) -> !isnothing(findfirst(isequal(v), l)), d)
+        new(sys, mw, molecs, 
+            filter_present(solids, active_solids), 
+            filter_present(gases, active_gases), 
+            filter_present(ions, active_ions_including_salts), 
+            filter_present(salts, active_salts),
+        )
     end
 end
 
 """
 Create a system of equations for mass and charge balances, 
 where slt, g, sld, and i are dictionaries of gases, solids, and ions, respectively.
 """
-function Balance(t, T, active_specs, active_ions, active_salts, g, sld, i)
+function Balance(t, T, active_specs, active_ions_including_salts, g, sld, i)
     # Return the given species, or zero if the species is not in the list.
-    all_active_ions = unique(vcat(active_ions, [[s.cation, s.anion] for s ∈ active_salts]...))
-    is(v) = ((v ∈ active_specs) || (v ∈ all_active_ions)) ? only(vars(v)) : 0
+
+    is(v) = ((v ∈ active_specs) || (v ∈ active_ions_including_salts)) ? only(vars(v)) : 0
 
     @constants R = 8.31446261815324 [unit = u"m_air^3*Pa/K/mol", description = "Universal gas constant"]
     @constants PaPerAtm = 101325 [unit = u"Pa/atm", description = "Number of pascals per atmosphere"]
@@ -184,12 +192,12 @@ function Balance(t, T, active_specs, active_ions, active_salts, g, sld, i)
         totalCl ~ is(i[:Cl]) + is(i[:HCl]) + atm2mol(is(g[:HCl])) +
                   is(sld[:NH4Cl]) + is(sld[:NaCl]) + 2is(sld[:CaCl2]) + is(sld[:KCl]) + 2is(sld[:MgCl2])
         # Charge balance
-        charge ~ sum([i.m * i.z for i in all_active_ions])
+        charge ~ sum([i.m * i.z for i in active_ions_including_salts])
 
         # Ratios from Section 3.1
-        R_1 ~ (totalNH + totalCa + totalK + totalMg + totalNa) / totalSO4
-        R_2 ~ (totalCa + totalK + totalMg + totalNa) / totalSO4
-        R_3 ~ (totalCa + totalK + totalMg) / totalSO4
+        # R_1 ~ (totalNH + totalCa + totalK + totalMg + totalNa) / totalSO4
+        # R_2 ~ (totalCa + totalK + totalMg + totalNa) / totalSO4
+        # R_3 ~ (totalCa + totalK + totalMg) / totalSO4
     ]
     nonzeros = [e.rhs !== 0 for e ∈ eqs]
     ODESystem(eqs[nonzeros], t, vs, [], name=:balance)

src/isorropia/reactions.jl

@@ -14,7 +14,8 @@ struct Rxn
     name
 end
 
-@constants tinyconc = 1e-12 [unit = u"mol/m_air^3", description = "Tiny concentration"]
+
+@constants tinyconc = 1e-11 [unit = u"mol/m_air^3", description = "Tiny concentration"]
 @constants tinypress = 1e-12 [unit = u"atm", description = "Tiny concentration"]
 
 """ Return whether the given concentration(s) are significantly greater than zero. """
@@ -24,6 +25,7 @@ present(s::Solid) = s.m > tinyconc
 present(s::SaltLike) = min(s.cation.m > tinyconc, s.anion.m > tinyconc) > 0
 present(v::AbstractVector) = min(present.(v)...) > 0
 present(_::H2O) = true
+bothabsent(r::Rxn) = max(present(r.reactant), present(r.product)) == 0
 
 """
 Create an equation system for this reaction based on the information in 
@@ -48,27 +50,33 @@ function rxn_sys(r::Rxn, t, activities::ModelingToolkit.AbstractSystem)
     @variables K_eq(t) [unit = r.K⁰units, description = "Equilibrium constant"]
     @variables rawrate(t) [unit = r.K⁰units, description = "Pseudo reaction rate"]
     @variables rate(t) [unit = r.K⁰units, description = "Normalized Pseudo reaction rate"]
-    @constants rateconst = 1 [unit = r.K⁰units, description = "Rate constant"]
-    @constants zerorate = 0 [unit = r.K⁰units, description = "Rate constant"]
-
-    rate_rhs = clamp(rawrate, -rateconst, rateconst) # Clip the mass transfer rate to avoid stiffness in the overall system.
-    # Products or reactants can only be consumed if their concentration is significantly greater than zero.
-    rate_rhs = ifelse(present(r.reactant), rate_rhs, max(rate_rhs, zerorate))
-    rate_rhs = ifelse(present(r.product), rate_rhs, min(rate_rhs, zerorate))
-
+    @constants rateconst =  1e-6 [unit = r.K⁰units, description = "Rate constant (chosen to manage stiffness)"]
+    @constants zerorate = 0 [unit = r.K⁰units, description = "Zero rate"]
+    @variables present(t) = 1 [description = "Whether the reactant is present (only used when reactant is a solid)"]
+    @constants unitconc = 1 [unit = u"mol/m_air^3", description = "Unit concentration"]
+    logit(x) = 1 / (1 + exp(-(NaNMath.log10(x/unitconc)+11)*5)) # Solid is not present if its concentration is less than ~1e-11
+
+#    rate_rhs = 
+ #   if (r.reactant isa Solid) 
+  #      rate_rhs = ifelse(present(r.reactant), rate_rhs, max(rate_rhs, zerorate))
+   # end
+
     sys = ODESystem([
-        # Equation 5: Equilibrium constant
-        K_eq ~ K⁰ * exp(-H_group * (T₀ / T - 1) - C_group * (1 + log(T₀ / T) - T₀ / T))
-        rawrate ~ (ap/ar - K_eq)
-        rate ~ rate_rhs
-    ], t, [K_eq, rawrate, rate], []; name=r.name)
+            # Equation 5: Equilibrium constant
+            K_eq ~ K⁰ * exp(-H_group * (T₀ / T - 1) - C_group * (1 + log(T₀ / T) - T₀ / T))
+            rawrate ~ (ap / ar - K_eq)
+            # Solids can only be consumed if their concentration is significantly greater than zero.
+            present ~ (r.reactant isa Solid) ? logit(r.reactant.m) : 1
+            # Clip the mass transfer rate to avoid stiffness in the overall system.
+            rate ~ tanh(rawrate / K⁰) * rateconst * present
+        ], t, [K_eq, rawrate, rate, present], []; name=r.name)
 
     # Equations to move toward equilibrium
     ode_eqs = Dict()
     for (v, s, sign) ∈ zip(vcat(pv, rv), vcat(ps, rs), vcat(fill(-1, length(pv)), fill(1, length(rv))))
-            x = Symbol(r.name, :conv_, Symbolics.tosymbol(v, escape=false))
-            conv = only(@constants $x=1 [unit = ModelingToolkit.get_unit(v/rate/t), description = "Unit onversion factor"])
-            ode_eqs[Dt(v)] = sign * sys.rate * s * conv
+        x = Symbol(r.name, :conv_, Symbolics.tosymbol(v, escape=false))
+        conv = only(@constants $x = 1 [unit = ModelingToolkit.get_unit(v / rate / t), description = "Unit onversion factor"])
+        ode_eqs[Dt(v)] = sign * sys.rate * s * conv
     end
     sys, ode_eqs
 end

src/isorropia/solid.jl

@@ -5,12 +5,12 @@ struct Solid <: Species
     m
 end
 
-#@constants unit_conc = 1.0 [unit = u"mol/m_air^3", description = "Unit concentration"]
-
 """
-From Section 2.2 in Fountoukis and Nenes (2007), the activity of a solid is 1.
+From Section 2.2 in Fountoukis and Nenes (2007), the activity of a solid is 1. 
+However, this cannot be strictly true because if the activity were always equal to one the solid would continue to be consumed 
+even after its concentration had reached zero, resulting in negative concentrations which are physically impossible.
 """
-activity(s::Solid) = 1.0
+activity(s::Solid) = 1
 
 vars(s::Solid) = [s.m]
 terms(s::Solid) = [s.m], [1]

test/isorropia/isorropia_test.jl

@@ -4,11 +4,17 @@ using Test
 using Plots
 
 include(joinpath(@__DIR__, "../../src/isorropia/isorropia.jl"))
-using .ISORROPIA
+#using .ISORROPIA
 
 @variables t [unit = u"s", description = "Time"]
 
-model = Isorropia(t, :all);
+#model = Isorropia(t, :all);
+#rxn_nums = [1]
+#rxn_nums = [1, 2, 3, 4, 6, 7, 8, 9, 10]
+#rxn_nums = [2, 7, 10]
+#rxn_nums = [7, 10]
+rxn_nums = 1:27
+model = Isorropia(t, rxn_nums);
 
 sys = structural_simplify(get_mtk(model))
 
@@ -22,32 +28,63 @@ u₀ = Dict{Any,Float64}([s => 1.0e-8 for s ∈ states(sys)])
 
 p = Dict{Any,Float64}([p => defaults[p] for p ∈ parameters(sys)])
 p[RH] = 0.30
-prob = ODEProblem(sys, u₀, (0.0, 3e-8), p)
+prob = ODEProblem(sys, u₀, (0.0, 1e-2), p)
+cb = PositiveDomain(zeros(length(prob.u0)))
 # Need low tolerance for mass balance checks to pass.
-@time sol = solve(prob, Rosenbrock23(), abstol=1e-12, reltol=1e-12)
+@time sol = solve(prob, Rosenbrock23(), abstol=1e-12, reltol=1e-12; callback = cb)
 
-begin
-    p1 = plot(title="Solids", xscale=:log10)
+let
+    xscale = :log10
+    yscale = :log10
+    p1 = plot(title="Solids", xscale=xscale, yscale=yscale, legend=:bottomleft)
     for (n, s) in model.solids
         plot!(sol.t[2:end], sol[s.m, 2:end], label=string(n))
     end
-    p2 = plot(title="Aqueous Ions", xscale=:log10)
+    p2 = plot(title="Aqueous Ions", xscale=xscale, yscale=yscale, legend=:bottomright)
     for (n, i) in model.ions
         plot!(sol.t[2:end], sol[i.m, 2:end], label=string(n))
     end
-    p3 = plot(title="Gases", xscale=:log10)
+    p3 = plot(title="Gases", xscale=xscale, yscale=yscale)
     for (n, g) in model.gases
         plot!(sol.t[2:end], sol[g.p, 2:end], label=string(n))
     end
-    p4 = plot(title="Reaction Rates", xscale=:log10)
-    for i ∈ 1:27
+    p4 = plot(title="Reaction Rates", xscale=xscale, legend=:outertopright)
+    for i ∈ rxn_nums
         r = Symbol(:rxn, i)
         y = eval(:(sol[sys.$r.rate]))
         plot!(sol.t[2:end], y[2:end], label="$r.rate")
-    end    
+    end
     plot(p1, p2, p3, p4, size=(1000, 800))
 end
 
+plot(
+    plot(sol.t, sol[sys.a_CaNO32_aqs]),
+    plot(sol.t, sol[sys.a_CaNO32_s]),
+    begin
+        plot(sol.t, sol[sys.a_CaNO32_aqs] ./ sol[sys.a_CaNO32_s], label="a ratio")
+        #plot!(sol.t, sol[sys.a_CaNO32_s] ./ sol[sys.a_CaNO32_aqs], yscale=:log10)
+        plot!(sol.t, sol[sys.rxn1.K_eq], label="K_eq")
+    end,
+)
+plot(
+    plot(sol.t, sol[sys.a_CaNO32_aqs] ./ sol[sys.a_CaNO32_s] - sol[sys.rxn1.K_eq], label="a ratio - k_eq"),
+    plot(sol.t, sol[sys.rxn1.rate], label="rate"),
+    plot(sol.t, sol[sys.rxn1.rawrate], label="rawrate"),
+)
+
+
+
+
+let
+    ps = []
+    for i ∈ rxn_nums
+        r = Symbol(:rxn, i)
+        y = eval(:(sol[sys.$r.rate]))
+        p = plot(sol.t[2:end], y[2:end], xscale=:log10, label="$r.rate")
+        push!(ps, p)
+    end
+    plot(ps..., size=(1000, 800))
+end
 
 plot(
     plot(sol[t], sol[f_CaNO32],