From 94792c7c5f36bd746af4936acb8a6a52c79055de Mon Sep 17 00:00:00 2001 From: Christopher Tessum Date: Tue, 31 Oct 2023 16:58:06 -0500 Subject: [PATCH] isorropia progress --- Project.toml | 1 + docs/src/isorropia/implementation.md | 31 +++++++++++++++ src/isorropia/isorropia.jl | 30 ++++++++------ src/isorropia/reactions.jl | 42 ++++++++++++-------- src/isorropia/solid.jl | 8 ++-- test/isorropia/isorropia_test.jl | 59 ++++++++++++++++++++++------ 6 files changed, 128 insertions(+), 43 deletions(-) diff --git a/Project.toml b/Project.toml index f5418fb2..772a7459 100644 --- a/Project.toml +++ b/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" diff --git a/docs/src/isorropia/implementation.md b/docs/src/isorropia/implementation.md index 8b372166..eb1c9186 100644 --- a/docs/src/isorropia/implementation.md +++ b/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}}`` diff --git a/src/isorropia/isorropia.jl b/src/isorropia/isorropia.jl index fe727836..e540d15b 100644 --- a/src/isorropia/isorropia.jl +++ b/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,7 +139,14 @@ 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 @@ -146,10 +154,10 @@ 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) diff --git a/src/isorropia/reactions.jl b/src/isorropia/reactions.jl index 18c8f6eb..8eaad983 100644 --- a/src/isorropia/reactions.jl +++ b/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 diff --git a/src/isorropia/solid.jl b/src/isorropia/solid.jl index 71a0488d..e2c0ef13 100644 --- a/src/isorropia/solid.jl +++ b/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] diff --git a/test/isorropia/isorropia_test.jl b/test/isorropia/isorropia_test.jl index 2147582d..11f33356 100644 --- a/test/isorropia/isorropia_test.jl +++ b/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],