diff --git a/src/isorropia/aqueous.jl b/src/isorropia/aqueous.jl index 8a41b6ae..fc3b575a 100644 --- a/src/isorropia/aqueous.jl +++ b/src/isorropia/aqueous.jl @@ -11,32 +11,30 @@ function Base.nameof(i::Ion) end vars(i::Ion) = [i.m] +terms(i::Ion) = [i.m], [1] #== The activity coefficient of an ion is assumed to be one (Fountoukis and Nenes (2007), Section 3.3). ==# activity(i::Ion, W) = i.m / W -# Generate the aqueous ions. -# Each ion has an associated MTK variable named -# _aq, where is the name of the compound, and -# a Ion struct named _ion. -ion_names = [:NH4, :Na, :H, :Cl, :NO3, :SO4, :HNO3, :NH3, :HCl, :HSO4, :Ca, :K, :Mg, :OH] -ion_valence = [1, 1, 1, -1, -1, -2, 0, 0, 0, -1, 2, 1, 2, -1] -ion_charge = [1, 1, 1, -1, -1, -2, 0,] -all_ions = [] -all_Ions = [] -for i ∈ eachindex(ion_names) - n = Symbol(ion_names[i], "_ion") - varname = Symbol(ion_names[i], "_aq") - s = "Aqueous $(ion_names[i])" - eval(quote - @species $varname = 1e-8 [unit = u"mol/m_air^3", description = $s] - $varname = ParentScope($varname) - push!(all_ions, $varname) - $n = Ion($varname, $(ion_valence[i])) - push!(all_Ions, $n) - end) +""" +Generate the aqueous ions. +Each ion has an associated MTK variable named +_aq, where is the name of the compound. +""" +function generate_ions(t) + ion_names = [:NH4, :Na, :H, :Cl, :NO3, :SO4, :HNO3, :NH3, :HCl, :HSO4, :Ca, :K, :Mg, :OH] + ion_valence = [1, 1, 1, -1, -1, -2, 0, 0, 0, -1, 2, 1, 2, -1] + ions = Dict() + for i ∈ eachindex(ion_names) + varname = Symbol(ion_names[i], "_aq") + s = "Aqueous $(ion_names[i])" + var = only(@variables $varname(t) = 1e-8 [unit = u"mol/m_air^3", description = s]) + var = ParentScope(var) + ions[ion_names[i]] = Ion(var, ion_valence[i]) + end + ions end abstract type SaltLike <: Species end @@ -76,42 +74,47 @@ function Base.nameof(s::SaltLike) end vars(s::SaltLike) = [s.cation.m, s.anion.m] - -# Salts from Table 4. -CaNO32_aqs = Salt(Ca_ion, 1, NO3_ion, 2, 0.4906, 509.4, 0.93) -CaCl2_aqs = Salt(Ca_ion, 1, Cl_ion, 2, 0.2830, 551.1, 2.4) -# CaSO4 and KHSO4 are below as SpecialSalts. -K2SO4_aqs = Salt(K_ion, 2, SO4_ion, 1, 0.9751, 35.6, -0.25) -KNO3_aqs = Salt(K_ion, 1, NO3_ion, 1, 0.9248, missing, -2.33) -KCl_aqs = Salt(K_ion, 1, Cl_ion, 1, 0.8426, 158.9, 0.92) -MgSO4_aqs = Salt(Mg_ion, 1, SO4_ion, 1, 0.8613, -714.5, 0.15) -MgNO32_aqs = Salt(Mg_ion, 1, NO3_ion, 2, 0.5400, 230.2, 2.32) -MgCl2_aqs = Salt(Mg_ion, 1, Cl_ion, 2, 0.3284, 42.23, 2.90) -NaCl_aqs = Salt(Na_ion, 1, Cl_ion, 1, 0.7528, 25.0, 2.23) -Na2SO4_aqs = Salt(Na_ion, 2, SO4_ion, 1, 0.9300, 80.0, -0.19) -NaNO3_aqs = Salt(Na_ion, 1, NO3_ion, 1, 0.7379, 304.0, -0.39) -NH42SO4_aqs = Salt(NH4_ion, 2, SO4_ion, 1, 0.7997, 80.0, -0.25) -NH4NO3_aqs = Salt(NH4_ion, 1, NO3_ion, 1, 0.6183, 852.0, -1.15) -NH4Cl_aqs = Salt(NH4_ion, 1, Cl_ion, 1, 0.7710, 239.0, 0.82) -# NH4HSO4, NaHSO4, and NH43HSO42 are below as SpecialSalts. -H2SO4_aqs = Salt(H_ion, 2, SO4_ion, 1, 0.000, missing, -0.1) -HHSO4_aqs = Salt(H_ion, 1, HSO4_ion, 1, 0.000, missing, 8.00) -HNO3_aqs = Salt(H_ion, 1, NO3_ion, 1, NaN, missing, 2.60) # There is no aqueous to solid conversion for HNO3. -HCl_aqs = Salt(H_ion, 1, Cl_ion, 1, NaN, missing, 6.00) # There is no aqueous to solid conversion for HCl. -all_salts = SaltLike[CaNO32_aqs, CaCl2_aqs, K2SO4_aqs, KNO3_aqs, KCl_aqs, MgSO4_aqs, - MgNO32_aqs, MgCl2_aqs, NaCl_aqs, Na2SO4_aqs, NaNO3_aqs, NH42SO4_aqs, NH4NO3_aqs, - NH4Cl_aqs, H2SO4_aqs, HHSO4_aqs, HNO3_aqs, HCl_aqs] +terms(s::SaltLike) = [s.cation.m, s.anion.m], [s.ν_cation, s.ν_anion] + +"Generate Salts from Table 4, where `i` is a dictionary of ions." +function generate_salts(i) + Dict( + :CaNO32 => Salt(i[:Ca], 1, i[:NO3], 2, 0.4906, 509.4, 0.93), + :CaCl2 => Salt(i[:Ca], 1, i[:Cl], 2, 0.2830, 551.1, 2.4), + :CaSO4 => CaSO4aqs(i[:Ca], 1, i[:SO4], 1, 0.9700, missing, NaN), + :KHSO4 => KHSO4aqs(i[:K], 1, i[:HSO4], 1, 0.8600, missing, NaN), + :K2SO4 => Salt(i[:K], 2, i[:SO4], 1, 0.9751, 35.6, -0.25), + :KNO3 => Salt(i[:K], 1, i[:NO3], 1, 0.9248, missing, -2.33), + :KCl => Salt(i[:K], 1, i[:Cl], 1, 0.8426, 158.9, 0.92), + :MgSO4 => Salt(i[:Mg], 1, i[:SO4], 1, 0.8613, -714.5, 0.15), + :MgNO32 => Salt(i[:Mg], 1, i[:NO3], 2, 0.5400, 230.2, 2.32), + :MgCl2 => Salt(i[:Mg], 1, i[:Cl], 2, 0.3284, 42.23, 2.90), + :NaCl => Salt(i[:Na], 1, i[:Cl], 1, 0.7528, 25.0, 2.23), + :Na2SO4 => Salt(i[:Na], 2, i[:SO4], 1, 0.9300, 80.0, -0.19), + :NaNO3 => Salt(i[:Na], 1, i[:NO3], 1, 0.7379, 304.0, -0.39), + :NH42SO4 => Salt(i[:NH4], 2, i[:SO4], 1, 0.7997, 80.0, -0.25), + :NH4NO3 => Salt(i[:NH4], 1, i[:NO3], 1, 0.6183, 852.0, -1.15), + :NH4Cl => Salt(i[:NH4], 1, i[:Cl], 1, 0.7710, 239.0, 0.82), + :NH4HSO4 => NH4HSO4aqs(i[:NH4], 1, i[:HSO4], 1, 0.4000, 384.0, NaN), + :NaHSO4 => NaHSO4aqs(i[:Na], 1, i[:HSO4], 1, 0.5200, -45.0, NaN), + :NH43HSO42 => NH43HSO42aqs(i[:NH4], 3, i[:HSO4], 2, 0.6900, 186.0, Inf), + :H2SO4 => Salt(i[:H], 2, i[:SO4], 1, 0.000, missing, -0.1), + :HHSO4 => Salt(i[:H], 1, i[:HSO4], 1, 0.000, missing, 8.00), + :HNO3 => Salt(i[:H], 1, i[:NO3], 1, NaN, missing, 2.60), # There is no aqueous to solid conversion for HNO3. + :HCl => Salt(i[:H], 1, i[:Cl], 1, NaN, missing, 6.00), # There is no aqueous to solid conversion for HCl. + ) +end """ Add additional salts as necessary to all the calculation of the activity -of the SpecialSalts. +of the SpecialSalts, where s is a dictionary of *all* salts. """ -function augmented_salts(active_salts) +function augmented_salts(active_salts, s) extra_salts = Dict( # Additional species needed to calculate activity of the SpecialSalts - KHSO4_aqs => [HHSO4_aqs, KCl_aqs, HCl_aqs], - NH4HSO4_aqs => [HHSO4_aqs, NH4Cl_aqs, HCl_aqs], - NaHSO4_aqs => [HHSO4_aqs, NaCl_aqs, HCl_aqs], - NH43HSO42_aqs => [NH42SO4_aqs, NH4HSO4_aqs], + s[:KHSO4] => [s[:HHSO4], s[:KCl], s[:HCl]], + s[:NH4HSO4] => [s[:HHSO4], s[:NH4Cl], s[:HCl]], + s[:NaHSO4] => [s[:HHSO4], s[:NaCl], s[:HCl]], + s[:NH43HSO42] => [s[:NH42SO4], s[:NH4HSO4]], ) for es in keys(extra_salts) if es in active_salts @@ -124,16 +127,16 @@ end """ Find all salts that have the same cation as the given salt. """ -function same_cation(s::SaltLike, active_salts) - aug_salts = augmented_salts(active_salts) +function same_cation(s::SaltLike, active_salts, all_salt_dict) + aug_salts = augmented_salts(active_salts, all_salt_dict) aug_salts[[s.cation == ss.cation for ss in aug_salts]] end """ Find all salts that have the same anion as the given salt. """ -function same_anion(s::SaltLike, active_salts) - aug_salts = augmented_salts(active_salts) +function same_anion(s::SaltLike, active_salts, all_salt_dict) + aug_salts = augmented_salts(active_salts, all_salt_dict) aug_salts[[s.anion == ss.anion for ss in aug_salts]] end @@ -145,23 +148,22 @@ end @constants I_one = 1 [unit = u"mol/kg_water", description = "An ionic strength of 1"] # Equation 6 -logγ₁₂T⁰(s::Salt, active_salts, I, W) = -Aᵧ * (abs(s.cation.z) * abs(s.anion.z) * √I) / (√I_one + √I) + - (abs(s.cation.z) * abs(s.anion.z)) / (abs(s.cation.z) + abs(s.anion.z)) * - (F₁(s, active_salts, I, W) / abs(s.cation.z) + F₂(s, active_salts, I, W) / abs(s.anion.z)) +logγ₁₂T⁰(s::Salt, active_salts, sss, I, W) = -Aᵧ * (abs(s.cation.z) * abs(s.anion.z) * √I) / (√I_one + √I) + + (abs(s.cation.z) * abs(s.anion.z)) / (abs(s.cation.z) + abs(s.anion.z)) * + (F₁(s, active_salts, sss, I, W) / abs(s.cation.z) + F₂(s, active_salts, sss, I, W) / abs(s.anion.z)) # Equation 7 -F₁(s::Salt, active_salts, I, W) = sum([ +F₁(s::Salt, active_salts, sss, I, W) = sum([ Y(ss, I, W) * logγ⁰₁₂(ss, I) * √I_one + Aᵧ * √I / (√I_one + √I) * abs(ss.cation.z) * abs(ss.anion.z) * Y(ss, I, W) - for ss ∈ same_cation(s, active_salts) + for ss ∈ same_cation(s, active_salts, sss) ]) - # Equation 8 -F₂(s::Salt, active_salts, I, W) = sum([ +F₂(s::Salt, active_salts, sss, I, W) = sum([ X(ss, I, W) * logγ⁰₁₂(ss, I) * √I_one + Aᵧ * √I / (√I_one + √I) * abs(ss.cation.z) * abs(ss.anion.z) * X(ss, I, W) - for ss ∈ same_anion(s, active_salts) + for ss ∈ same_anion(s, active_salts, sss) ]) # Supplemental equations after 7 and 8 @@ -184,15 +186,15 @@ C(q, I) = 1 + 0.055q * exp(-0.023I^3 / I_one^3) # Equation 14 A(I) = -((0.41√I / (√I_one + √I)) + 0.039(I / I_one)^0.92) -logγ₁₂(s::Salt, active_salts, I, W) = (1.125 - c_1 * (T - T₀₂)) * logγ₁₂T⁰(s, active_salts, I, W) / - √I_one - (0.125 - c_1 * (T - T₀₂)) * A(I) +logγ₁₂(s::Salt, active_salts, sss, I, W) = (1.125 - c_1 * (T - T₀₂)) * logγ₁₂T⁰(s, active_salts, sss, I, W) / + √I_one - (0.125 - c_1 * (T - T₀₂)) * A(I) """ Calculate the activity coefficient of a salt based on Section 2.2 in Fountoukis and Nenes (2007). """ -activity(s::Salt, active_salts, I, W) = (s.cation.m/W)^s.ν_cation * (s.anion.m/W)^s.ν_anion * - min(exp(logγ₁₂(s, active_salts, I, W)),1)^(s.ν_cation + s.ν_anion) +activity(s::Salt, active_salts, sss, I, W) = (s.cation.m / W)^s.ν_cation * (s.anion.m / W)^s.ν_anion * + min(exp(logγ₁₂(s, active_salts, sss, I, W)), 1)^(s.ν_cation + s.ν_anion) # TODO(CT): Clipping γ to 1. This may cause problems. ### Special cases @@ -203,23 +205,14 @@ abstract type SpecialSalt <: SaltLike end The activity of a SpecialSalt is the same as for a salt except that it has a special activity coefficient as defined in the footnotes to Table 4. """ -activity(s::SpecialSalt, active_salts, I, W) = (s.cation.m/W)^s.ν_cation * (s.anion.m/W)^s.ν_anion * - min(γ₁₂(s, active_salts, I, W),1)^(s.ν_cation + s.ν_anion) +activity(s::SpecialSalt, active_salts, sss, I, W) = (s.cation.m / W)^s.ν_cation * (s.anion.m / W)^s.ν_anion * + min(γ₁₂(s, active_salts, sss, I, W), 1)^(s.ν_cation + s.ν_anion) # TODO(CT): Clipping γ to 1. This may cause problems. logγ⁰₁₂(s::SpecialSalt, I) = abs(s.cation.z) * abs(s.anion.z) * log(Γ⁰(0, I) / I_one) # Assuming q = 0 for SpecialSalts +# Generate special salt structs for multiple dispatch. specialsaltnames = [:CaSO4, :KHSO4, :NH4HSO4, :NaHSO4, :NH43HSO42] -# Data from Fountoukis and Nenes (2007) Table 4 -specialsalts = [ - Salt(Ca_ion, 1, SO4_ion, 1, 0.9700, missing, NaN), - Salt(K_ion, 1, HSO4_ion, 1, 0.8600, missing, NaN), - Salt(NH4_ion, 1, HSO4_ion, 1, 0.4000, 384.0, NaN), - Salt(Na_ion, 1, HSO4_ion, 1, 0.5200, -45.0, NaN), - Salt(NH4_ion, 3, HSO4_ion, 2, 0.6900, 186.0, Inf)] - for (i, name) ∈ enumerate(specialsaltnames) - s = specialsalts[i] - varname = Symbol(name, "_aqs") structname = Symbol(name, "aqs") eval(quote """ @@ -236,71 +229,71 @@ for (i, name) ∈ enumerate(specialsaltnames) drh "Enthalpy term (-18/1000R L_s m_s)" l_term + "Kusik-Meissner Binary activity parameter" + q::Number end - $varname = $structname($(s.cation), $(s.ν_cation), $(s.anion), $(s.ν_anion), $(s.drh), $(s.l_term)) - push!(all_salts, $varname) end) end """ From Table 4 footnote a, CaSO4 has an activity coefficient of zero. """ -γ₁₂(s::CaSO4aqs, active_salts, I, W) = 1.0e-20 +γ₁₂(s::CaSO4aqs, active_salts, sss, I, W) = 1.0e-20 """ From Table 4 footnote b, KHSO4 has a unique activity coefficient """ -γ₁₂(s::KHSO4aqs, active_salts, I, W) = (exp(logγ₁₂(HHSO4_aqs, active_salts, I, W)) * - exp(logγ₁₂(KCl_aqs, active_salts, I, W)) / - exp(logγ₁₂(HCl_aqs, active_salts, I, W)))^(1 / 2) +γ₁₂(s::KHSO4aqs, active_salts, sss, I, W) = (exp(logγ₁₂(sss[:HHSO4], active_salts, sss, I, W)) * + exp(logγ₁₂(sss[:KCl], active_salts, sss, I, W)) / + exp(logγ₁₂(sss[:HCl], active_salts, sss, I, W)))^(1 / 2) """ From Table 4 footnote c, NH4HSO4 has a unique activity coefficient """ -γ₁₂(s::NH4HSO4aqs, active_salts, I, W) = (exp(logγ₁₂(HHSO4_aqs, active_salts, I, W)) * - exp(logγ₁₂(NH4Cl_aqs, active_salts, I, W)) / - exp(logγ₁₂(HCl_aqs, active_salts, I, W)))^(1 / 2) +γ₁₂(s::NH4HSO4aqs, active_salts, sss, I, W) = (exp(logγ₁₂(sss[:HHSO4], active_salts, sss, I, W)) * + exp(logγ₁₂(sss[:NH4Cl], active_salts, sss, I, W)) / + exp(logγ₁₂(sss[:HCl], active_salts, sss, I, W)))^(1 / 2) """ From Table 4 footnote d, NaHSO4 has a unique activity coefficient """ -γ₁₂(s::NaHSO4aqs, active_salts, I, W) = (exp(logγ₁₂(HHSO4_aqs, active_salts, I, W)) * - exp(logγ₁₂(NaCl_aqs, active_salts, I, W)) / - exp(logγ₁₂(HCl_aqs, active_salts, I, W)))^(1 / 2) +γ₁₂(s::NaHSO4aqs, active_salts, sss, I, W) = (exp(logγ₁₂(sss[:HHSO4], active_salts, sss, I, W)) * + exp(logγ₁₂(sss[:NaCl], active_salts, sss, I, W)) / + exp(logγ₁₂(sss[:HCl], active_salts, sss, I, W)))^(1 / 2) """ From Table 4 footnote e, NH43HSO42 has a unique activity coefficient """ -γ₁₂(s::NH43HSO42aqs, active_salts, I, W) = (exp(logγ₁₂(NH42SO4_aqs, active_salts, I, W))^3 * - γ₁₂(NH4HSO4_aqs, active_salts, I, W))^(1 / 5) +γ₁₂(s::NH43HSO42aqs, active_salts, sss, I, W) = (exp(logγ₁₂(sss[:NH42SO4], active_salts, sss, I, W))^3 * + γ₁₂(sss[:NH4HSO4], active_salts, sss, I, W))^(1 / 5) """ Create a system of equations for the ionic strength of the multicomponent solution, as specified by Fountoukis and Nenes (2007), between equations 8 and 9: ``I = \\frac{1}{2} \\sum_i m_i z_i^2`` """ -function IonicStrength(all_Ions, W) - @variables I = 1.0e-4 [unit = u"mol/kg_water", description = "Ionic strength"] +function IonicStrength(t, all_Ions, W) + @variables I(t) = 1.0e-4 [unit = u"mol/kg_water", description = "Ionic strength"] I = ParentScope(I) @constants I_one = 1 [unit = u"mol/kg_water", description = "An ionic strength of 1"] - @named ionic_strength = NonlinearSystem([ - I ~ 1 / 2 * sum([ion.m/W * ion.z^2 for ion in all_Ions]) - ], [I], []) + @named ionic_strength = ODESystem([ + I ~ 1 / 2 * sum([ion.m / W * ion.z^2 for ion in all_Ions]) + ], t, [I], []) end """ Create an equation system for the activities of all the species listed in all_species, where `f` is a function that takes a species as its argument and returns its activity. """ -function Activities(all_species, f) +function Activities(t, all_species, f) eqs = Equation[] vars = [] for s ∈ all_species rhs = f(s) unit = ModelingToolkit.get_unit(rhs) x = Symbol("a_", nameof(s)) - a = only(@variables $x=1 [unit = unit, description = "Activity of $(nameof(s))"]) + a = only(@variables $x(t) = 1 [unit = unit, description = "Activity of $(nameof(s))"]) push!(vars, a) push!(eqs, a ~ rhs) end - NonlinearSystem(eqs, vars, []; name=:activities) + ODESystem(eqs, t, vars, []; name=:activities) end \ No newline at end of file diff --git a/src/isorropia/gas.jl b/src/isorropia/gas.jl index efa0baf5..178b3292 100644 --- a/src/isorropia/gas.jl +++ b/src/isorropia/gas.jl @@ -14,23 +14,21 @@ function Base.nameof(g::Gas) end vars(g::Gas) = [g.p] +terms(g::Gas) = [g.p], [1] -# Generate the gases. -# Each gas has an associated MTK variable named -# _g, where is the name of the compound, and -# a Gas struct named g. -# All SO4 immediately goes to aerosol phase as per Section 3.3 (item 1) of Fountoukis and Nenes (2007), so we don't include it here. -all_gases = [] -all_Gases = [] -for (s, v) ∈ [:HNO3 => 1e-8, :HCl => 1e-8, :NH3 => 1e-8] - varname = Symbol(s, "_g") - gasname = Symbol(s, "g") - description = "Gasous $s" - eval(quote - @species $varname = $v [unit = u"atm", description = $description] - $varname = ParentScope($varname) - push!($all_gases, $varname) - $gasname = $Gas($varname) - push!(all_Gases, $gasname) - end) +""" +Generate the gases. +Each gas has an associated MTK variable named _g, where is the name of the compound. +All SO4 immediately goes to aerosol phase as per Section 3.3 (item 1) of Fountoukis and Nenes (2007), so we don't include it here. +""" +function generate_gases(t) + gases = Dict() + for (s, v) ∈ [:HNO3 => 1e-8, :HCl => 1e-8, :NH3 => 1e-8] + varname = Symbol(s, "_g") + description = "Gasous $s" + var = only(@species $varname(t) = v [unit = u"atm", description = description]) + var = ParentScope(var) + gases[s] = Gas(var) + end + gases end \ No newline at end of file diff --git a/src/isorropia/isorropia.jl b/src/isorropia/isorropia.jl index c939ac65..fe727836 100644 --- a/src/isorropia/isorropia.jl +++ b/src/isorropia/isorropia.jl @@ -12,7 +12,7 @@ export Isorropia include("units.jl") include("species.jl") include("aqueous.jl") -include("deliquescence.jl") +#include("deliquescence.jl") include("solid.jl") include("gas.jl") include("water.jl") @@ -26,14 +26,41 @@ Substitute in transformed concentrations to avoid negative values, where `sys` is an equation system and `active_vars` is a vector of variables to be substituted. """ -function sub_trans(sys::ModelingToolkit.NonlinearSystem, active_vars::AbstractVector) +function sub_trans(sys::ModelingToolkit.ODESystem, active_vars::AbstractVector) trans_subs = [v => trans(v) for v ∈ active_vars] - subbed_eqs = substitute(equations(sys), Dict(trans_subs)) - NonlinearSystem(subbed_eqs, states(sys), parameters(sys); name=nameof(sys)) + subbed_eqs = substitute(equations(sys), Dict(trans_subs)) # Type 1 reactions (R_1 < 1.0) + ODESystem(subbed_eqs, ModelingToolkit.get_iv(sys), states(sys), parameters(sys); name=nameof(sys)) end +# Molar mass in g/mol +mw = Dict(:Na_aq => 22.989769, :SO4_aq => 96.0636, :NH3_aq => 17.03052, :NH3_g => 17.03052, + :NO3_aq => 62.0049, :Cl_aq => 35.453, :NaCl_s => 58.44, + :Ca_aq => 40.078, :K_aq => 39.0983, :Mg_aq => 24.305, :H_aq => 1.00784, :NH4_aq => 18.04, :HCl_g => 36.46, + :K2SO4_s => :174.259, :KNO3_s => 101.1032, :CaNO32_s => 164.1, :HNO3_g => 63.01, :HNO3_aq => 63.01, + :KHSO4_s => 136.169, :KCl_s => 74.5513, :NH4NO3_s => 80.043) + +# Species containing each molecule for mass balancing +molecs = Dict( + :K => [(1, :K_aq), (1, :KHSO4_s), (2, :K2SO4_s), (1, :KNO3_s), (1, :KCl_s)], + :Ca => [(1, :Ca_aq), (1, :CaSO4_s), (1, :CaNO32_s), (1, :CaCl2_s)], + :Mg => [(1, :Mg_aq), (1, :MgSO4_s), (1, :MgNO32_s), (1, :MgCl2_s)], + :NH => [(1, :NH4_aq), (1, :NH3_aq), (1, :NH3_g), (1, :NH4HSO4_s), + (2, :NH42SO4_s), (3, :NH43HSO42_s), (1, :NH4Cl_s), (1, :NH4NO3_s)], + :Na => [(1, :Na_aq), (1, :NaHSO4_s), (2, :Na2SO4_s), (1, :NaCl_s), (1, :NaNO3_s)], + :SO4 => [(1, :SO4_aq), (1, :HSO4_aq), + (1, :KHSO4_s), (1, :NaHSO4_s), (1, :NH4HSO4_s), (1, :CaSO4_s), (1, :Na2SO4_s), (1, :NH42SO4_s), + (2, :NH43HSO42_s), (1, :K2SO4_s), (1, :MgSO4_s)], + :NO3 => [(1, :NO3_aq), (1, :HNO3_aq), (1, :HNO3_g), (1, :NH4NO3_s), (1, :NaNO3_s)], + :Cl => [(1, :Cl_aq), (1, :HCl_aq), (1, :HCl_g), (1, :NH4Cl_s), + (1, :NaCl_s), (2, :CaCl2_s), (1, :KCl_s), (2, :MgCl2_s)], + :H => [(1, :H_aq), (1, :HNO3_g), (1, :HCl_g)], +) + """ - Isorropia(rxn_nums) + Isorropia(t) + Isorropia(t, 1:3) # Only include the first three of the 27 reactions. + Isorropia(t, 1) # Only include the first reaction. + Isorropia(t, :all) # Choose a pre-specified set of reactions. Current options are `:all` and `:type1` (referring to first aerosol type in Fountoukis and Nenes (2007) Table 3). An implementation of ISORROPIA II, a model for the thermodynamic equilibrium of gas-aerosol interactions, as described in: @@ -41,82 +68,123 @@ An implementation of ISORROPIA II, a model for the thermodynamic equilibrium of """ struct Isorropia <: EarthSciMLODESystem - sys::NonlinearSystem - Isorropia(sys::ModelingToolkit.NonlinearSystem) = new(sys) - function Isorropia(rxn_nums::AbstractVector, name=:isorropia) - rxns = reactions[rxn_nums] + "ModelingToolkit ODESystem" + sys + "Molar mass for each species in g/mol" + mw + "Species containing each molecule for mass balancing" + molecs + "Dictionary of solids" + solids + "Dictionary of gases" + gases + "Dictionary of aqueous ions" + ions + "Dictionary of salts" + salts + + Isorropia(t, name=:isorropia) = Isorropia(t, :all; name=name) + function Isorropia(t, type::Symbol; name=:isorropia) + sys_types = Dict( + :all => 1:27, # Use all 27 reactions + :type1 => [3, 5, 11, 12, 13, 14, 15, 16, 17, 18, 24, 26], # Type 1 reactions (R_1 < 1.0) + ) + @assert (type ∈ keys(sys_types)) "invalid system type $type, valid options are $(keys(sys_types))" + Isorropia(t, sys_types[type]; name=name) + end + Isorropia(t, rxn_num::Int; name=:isorropia) = Isorropia(t, [rxn_num]; name=name) + function Isorropia(t, rxn_nums::AbstractVector; name=:isorropia) + ions = generate_ions(t) + salts = generate_salts(ions) + gases = generate_gases(t) + solids = generate_solids(t) + H2Oaq = H2O(RH) + + rxns = generate_reactions(salts, gases, solids, ions, H2Oaq)[rxn_nums] active_specs = active_species(rxns) active_vars = unique(vcat(vars.(active_specs)...)) - active_ions = intersect(active_specs, all_Ions) - active_salts = intersect(active_specs, all_salts) - active_solids = intersect(active_specs, all_Solids) - active_gases = intersect(active_specs, all_Gases) - - water = Water(active_salts) - ionic_strength = IonicStrength(active_ions, water.W) - salt_act = Activities(active_salts, (s) -> activity(s, active_salts, ionic_strength.I, water.W)) - ion_act = Activities(active_ions, (i) -> activity(i, water.W)) - gas_act = Activities(active_gases, activity) - solid_act = Activities(active_solids, activity) - water_act = Activities([H2Oaq], activity) + active_ions = intersect(active_specs, values(ions)) + active_salts = intersect(active_specs, values(salts)) + active_solids = intersect(active_specs, values(solids)) + active_gases = intersect(active_specs, values(gases)) + + water = Water(t, active_salts) + ionic_strength = IonicStrength(t, active_ions, 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, active_vars, active_ions, active_salts) + balance = Balance(t, T, active_specs, active_ions, active_salts, gases, solids, ions) #out = Output(active_vars) - sys = NonlinearSystem([], active_vars, [T, RH], - systems=[rxn_sys(x, activities) for x ∈ rxns]; name=name, - ) + # Convert dictionaries of equations into an equation system by combining equations with the same left-hand side. + eq_dict = Dict() + syss = [] + for rx ∈ rxns + sys, ode_eqs_dict = rxn_sys(rx, t, activities) + push!(syss, sys) + for lhs ∈ keys(ode_eqs_dict) + if lhs in keys(eq_dict) + eq_dict[lhs] += ode_eqs_dict[lhs] + else + eq_dict[lhs] = ode_eqs_dict[lhs] + end + end + end + eqs = [lhs ~ rhs for (lhs, rhs) ∈ pairs(eq_dict)] + + 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) + #sys = sub_trans(sys, active_vars) + new(sys, mw, molecs, solids, gases, ions, salts) end end """ -Create a system of equations for mass and charge balances. +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, active_vars, active_ions, active_salts) +function Balance(t, T, active_specs, active_ions, active_salts, g, sld, i) # Return the given species, or zero if the species is not in the list. - g(v) = !isnothing(findfirst(isequal(v), active_vars)) ? v : 0 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 @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"] atm2mol(p) = p * PaPerAtm / R / T - totals = @parameters begin - totalK = 1.0e-7, [unit = u"mol/m_air^3", description = "Total concentration of K"] - totalCa = 1.0e-7, [unit = u"mol/m_air^3", description = "Total concentration of Ca"] - totalMg = 1.0e-7, [unit = u"mol/m_air^3", description = "Total concentration of Mg"] - totalNH = 1.0e-7, [unit = u"mol/m_air^3", description = "Total concentration of NH"] - totalNa = 1.0e-7, [unit = u"mol/m_air^3", description = "Total concentration of Na"] - totalSO4 = 1.0e-7, [unit = u"mol/m_air^3", description = "Total concentration of SO4"] - totalNO3 = 1.0e-7, [unit = u"mol/m_air^3", description = "Total concentration of NO3"] - totalCl = 1.0e-7, [unit = u"mol/m_air^3", description = "Total concentration of Cl"] - end - - ratios = @variables begin - R_1 = 1.0, [description = "Total sulfate ratio from Section 3.1 of Fountoukis and Nenes (2007)"] - R_2 = 1.0, [description = "Crustal species and sodium ratio from Section 3.1 of Fountoukis and Nenes (2007)"] - R_3 = 1.0, [description = "Crustal species ratio from Section 3.1 of Fountoukis and Nenes (2007)"] - end + @variables totalK(t) = 1.0e-7 [unit = u"mol/m_air^3", description = "Total concentration of K"] + @variables totalCa(t) = 1.0e-7 [unit = u"mol/m_air^3", description = "Total concentration of Ca"] + @variables totalMg(t) = 1.0e-7 [unit = u"mol/m_air^3", description = "Total concentration of Mg"] + @variables totalNH(t) = 1.0e-7 [unit = u"mol/m_air^3", description = "Total concentration of NH"] + @variables totalNa(t) = 1.0e-7 [unit = u"mol/m_air^3", description = "Total concentration of Na"] + @variables totalSO4(t) = 1.0e-7 [unit = u"mol/m_air^3", description = "Total concentration of SO4"] + @variables totalNO3(t) = 1.0e-7 [unit = u"mol/m_air^3", description = "Total concentration of NO3"] + @variables totalCl(t) = 1.0e-7 [unit = u"mol/m_air^3", description = "Total concentration of Cl"] + @variables R_1(t) = 1.0 [description = "Total sulfate ratio from Section 3.1 of Fountoukis and Nenes (2007)"] + @variables R_2(t) = 1.0 [description = "Crustal species and sodium ratio from Section 3.1 of Fountoukis and Nenes (2007)"] + @variables R_3(t) = 1.0 [description = "Crustal species ratio from Section 3.1 of Fountoukis and Nenes (2007)"] + @variables charge(t) = 0 [unit = u"mol/m_air^3", description = "Total ion charge"] + vs = [totalK, totalCa, totalMg, totalNH, totalNa, totalSO4, totalNO3, totalCl, R_1, R_2, R_3, charge] eqs = [ # Mass balances - totalK ~ g(K_aq) + g(KHSO4_s) + 2g(K2SO4_s) + g(KNO3_s) + g(KCl_s) - totalCa ~ g(Ca_aq) + g(CaSO4_s) + g(CaNO32_s) + g(CaCl2_s) - totalMg ~ g(Mg_aq) + g(MgSO4_s) + g(MgNO32_s) + g(MgCl2_s) - totalNH ~ g(NH4_aq) + g(NH3_aq) + atm2mol(g(NH3_g)) + g(NH4HSO4_s) + - 2g(NH42SO4_s) + 3g(NH43HSO42_s) + g(NH4Cl_s) + g(NH4NO3_s) - totalNa ~ g(Na_aq) + g(NaHSO4_s) + 2g(Na2SO4_s) + g(NaCl_s) + g(NaNO3_s) - totalSO4 ~ g(SO4_aq) + g(HSO4_aq) + - g(KHSO4_s) + g(NaHSO4_s) + g(NH4HSO4_s) + g(CaSO4_s) + g(Na2SO4_s) + g(NH42SO4_s) + - 2g(NH43HSO42_s) + g(K2SO4_s) + g(MgSO4_s) - totalNO3 ~ g(NO3_aq) + g(HNO3_aq) + atm2mol(g(HNO3_g)) + g(NH4NO3_s) + g(NaNO3_s) - totalCl ~ g(Cl_aq) + g(HCl_aq) + atm2mol(g(HCl_g)) + - g(NH4Cl_s) + g(NaCl_s) + 2g(CaCl2_s) + g(KCl_s) + 2g(MgCl2_s) + totalK ~ is(i[:K]) + is(sld[:KHSO4]) + 2is(sld[:K2SO4]) + is(sld[:KNO3]) + is(sld[:KCl]) + totalCa ~ is(i[:Ca]) + is(sld[:CaSO4]) + is(sld[:CaNO32]) + is(sld[:CaCl2]) + totalMg ~ is(i[:Mg]) + is(sld[:MgSO4]) + is(sld[:MgNO32]) + is(sld[:MgCl2]) + totalNH ~ is(i[:NH4]) + is(i[:NH3]) + atm2mol(is(g[:NH3])) + is(sld[:NH4HSO4]) + + 2is(sld[:NH42SO4]) + 3is(sld[:NH43HSO42]) + is(sld[:NH4Cl]) + is(sld[:NH4NO3]) + totalNa ~ is(i[:Na]) + is(sld[:NaHSO4]) + 2is(sld[:Na2SO4]) + is(sld[:NaCl]) + is(sld[:NaNO3]) + totalSO4 ~ is(i[:SO4]) + is(i[:HSO4]) + + is(sld[:KHSO4]) + is(sld[:NaHSO4]) + is(sld[:NH4HSO4]) + is(sld[:CaSO4]) + is(sld[:Na2SO4]) + is(sld[:NH42SO4]) + + 2is(sld[:NH43HSO42]) + is(sld[:K2SO4]) + is(sld[:MgSO4]) + totalNO3 ~ is(i[:NO3]) + is(i[:HNO3]) + atm2mol(is(g[:HNO3])) + is(sld[:NH4NO3]) + is(sld[:NaNO3]) + 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 - 0 ~ sum([i.m * i.z for i in all_active_ions]) + charge ~ sum([i.m * i.z for i in all_active_ions]) # Ratios from Section 3.1 R_1 ~ (totalNH + totalCa + totalK + totalMg + totalNa) / totalSO4 @@ -124,7 +192,7 @@ function Balance(T, active_vars, active_ions, active_salts) R_3 ~ (totalCa + totalK + totalMg) / totalSO4 ] nonzeros = [e.rhs !== 0 for e ∈ eqs] - NonlinearSystem(eqs[nonzeros], ratios, totals, name=:balance) + ODESystem(eqs[nonzeros], t, vs, [], name=:balance) end """ @@ -132,41 +200,9 @@ Create a system of equations for output variables. """ function Output(active_vars) names = Symbolics.tosymbol.(active_vars, escape=false) - ovars = [only(@variables $(n) [unit=ModelingToolkit.get_unit(v), description="Output for $n"]) for - (n, v) ∈ zip(names, active_vars)] + ovars = [only(@variables $(n) [unit = ModelingToolkit.get_unit(v), description = "Output for $n"]) for + (n, v) ∈ zip(names, active_vars)] # This ends up with a double transform (e.g. ||v||), but is necessary to avoid cancelling. eqs = [ov ~ trans(v) for (ov, v) ∈ zip(ovars, active_vars)] NonlinearSystem(eqs, ovars, [], name=:out) -end - -# Type 1 reactions (R_1 < 1.0) -rxn_nums = [3, 5, 11, 12, 13, 14, 15, 16, 17, 18, 24, 26] -sys = get_mtk(Isorropia(rxn_nums)) - -render(latexify(sys)) -simple_sys = structural_simplify(sys) - -render(latexify(simple_sys)) -render(latexify(observed(simple_sys))) - -ps = [simple_sys.totalSO4 => 1e-7, simple_sys.totalNa => 1e-20, simple_sys.totalNH => 1e-8, - simple_sys.totalNO3 => 1e-7, simple_sys.totalCl => 1e-20, simple_sys.totalCa => 5e-9, - simple_sys.totalK => 5e-9, simple_sys.totalMg => 1e-20] -prob = NonlinearProblem(simple_sys, [], ps) -sol = solve(prob, TrustRegion()) -sol[sys.R_1] # Confirm that R_1 < 1.0 -show(sol.resid) # Residuals should be small, but they are not. This means that the solver hasn't converged -sol[sys.W] -sol[sys.I] - -sts = [sys.CaSO4_s, sys.Ca_aq, sys.SO4_aq, sys.KHSO4_s, sys.K_aq, - sys.HSO4_aq, sys.H_aq, sys.NH3_g, sys.NH3_aq, sys.NH4_aq, - sys.OH_aq, sys.HNO3_g, sys.NO3_aq, sys.HNO3_aq, sys.HCl_g, - sys.Cl_aq, sys.HCl_aq, sys.NaHSO4_s, sys.Na_aq, sys.NH4HSO4_s] -show([v => trans(sol[v]) for v in sts]) # Negative numbers are fine because we take the absolute value. - -sts = [sys.a_CaSO4_s, sys.a_KHSO4_s, sys.a_NaHSO4_s, sys.a_NH4HSO4_s, - sys.a_NH3_g, sys.a_HNO3_g, sys.a_HCl_g, sys.a_HSO4_aq, sys.a_H_aq, sys.a_SO4_aq, - sys.a_NH3_aq, sys.a_NH4_aq, sys.a_OH_aq, sys.a_HNO3_aq, sys.a_Cl_aq, sys.a_HCl_aq, - sys.a_CaSO4_aqs, sys.a_KHSO4_aqs, sys.a_HNO3_aqs, sys.a_NaHSO4_aqs, sys.a_NH4HSO4_aqs, sys.a_H2O] -show([v => sol[v] for v in sts]) \ No newline at end of file +end \ No newline at end of file diff --git a/src/isorropia/reactions.jl b/src/isorropia/reactions.jl index 3c001a25..18c8f6eb 100644 --- a/src/isorropia/reactions.jl +++ b/src/isorropia/reactions.jl @@ -14,22 +14,63 @@ struct Rxn name end +@constants tinyconc = 1e-12 [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. """ +present(s::Gas) = s.p > tinypress +present(s::Ion) = s.m > tinyconc +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 + """ Create an equation system for this reaction based on the information in Equation 5 and Table 2 of Fountoukis and Nenes (2007). + +Rather than setting up for nonlinear rootfinding as is done in the original paper, +we create equations for a mass-action based solution where mass is moved between +the product and reactant species based on different between the current ratio and +the equilibrium ratio. """ -function rxn_sys(r::Rxn, activities::ModelingToolkit.AbstractSystem) +function rxn_sys(r::Rxn, t, activities::ModelingToolkit.AbstractSystem) + ap = speciesactivity(activities, r.product) + ar = speciesactivity(activities, r.reactant) + pv, ps = terms(r.product) # Product variables and stoichiometry coefficients + rv, rs = terms(r.reactant) # Reactant variables and stoichiometry coefficients + Dt = Differential(t) @constants T₀ = 293.15 [unit = u"K", description = "Standard temperature"] # These are the variables from Fountoukis and Nenes (2007) Table 2 @constants K⁰ = r.K⁰ [unit = r.K⁰units, description = "Equilibrium constant at 298.15 K"] @constants H_group = r.hgroup [description = "ΔH⁰ / (R * T₀) (unitless)"] @constants C_group = r.cgroup [description = "ΔC⁰ₚ / R (unitless)"] - @variables K_eq [unit = r.K⁰units, description = "Equilibrium constant"] - NonlinearSystem([ + @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)) + + sys = ODESystem([ # Equation 5: Equilibrium constant K_eq ~ K⁰ * exp(-H_group * (T₀ / T - 1) - C_group * (1 + log(T₀ / T) - T₀ / T)) - K_eq ~ speciesactivity(activities, r.product) / speciesactivity(activities, r.reactant) - ], [K_eq], [K⁰, H_group, C_group, T₀, T₀₂, c_1, I_one]; name=r.name) + rawrate ~ (ap/ar - K_eq) + rate ~ rate_rhs + ], t, [K_eq, rawrate, rate], []; 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 + end + sys, ode_eqs end @@ -53,37 +94,43 @@ as shown in Table 2 of Fountoukis and Nenes (2007). """ speciesactivity(activities::ModelingToolkit.AbstractSystem, s::AbstractVector) = reduce(*, speciesactivity.((activities,), s)) -# NOTE: Assuming that H_group and C_group are zero when they are left out of Table 2. -reactions = [ - Rxn(CaNO32s, CaNO32_aqs, 6.067e5, u"mol^3/kg_water^3", -11.299, 0.0, :rxn1) - Rxn(CaCl2s, CaCl2_aqs, 7.974e11, u"mol^3/kg_water^3", -14.087, 0.0, :rxn2) - Rxn(CaSO4s, CaSO4_aqs, 4.319e-5, u"mol^2/kg_water^2", 0.0, 0.0, :rxn3) - Rxn(K2SO4s, K2SO4_aqs, 1.569e-2, u"mol^3/kg_water^3", -9.589, 45.807, :rxn4) - Rxn(KHSO4s, KHSO4_aqs, 24.016, u"mol^2/kg_water^2", -8.423, 17.964, :rxn5) - Rxn(KNO3s, KNO3_aqs, 0.872, u"mol^2/kg_water^2", 14.075, 19.388, :rxn6) - Rxn(KCls, KCl_aqs, 8.680, u"mol^2/kg_water^2", -6.167, 19.953, :rxn7) - Rxn(MgSO4s, MgSO4_aqs, 1.079e5, u"mol^2/kg_water^2", 36.798, 0.0, :rxn8) - Rxn(MgNO32s, MgNO32_aqs, 2.507e15, u"mol^3/kg_water^3", -8.754, 0.0, :rxn9) - Rxn(MgCl2s, MgCl2_aqs, 9.557e21, u"mol^3/kg_water^3", -1.347, 0.0, :rxn10) - Rxn(HSO4_ion, [H_ion, SO4_ion], 1.015e-2, u"mol/kg_water", 8.85, 25.14, :rxn11) - Rxn(NH3g, NH3_ion, 5.764e1, u"mol/kg_water/atm", 13.79, -5.39, :rxn12) - Rxn([NH3_ion, H2Oaq], [NH4_ion, OH_ion], 1.805e-5, u"mol/kg_water", -1.50, 26.92, :rxn13) - Rxn(HNO3g, HNO3_aqs, 2.511e6, u"mol^2/kg_water^2/atm", 29.17, 16.83, :rxn14) - Rxn(HNO3g, HNO3_ion, 2.1e5, u"mol/kg_water/atm", 29.17, 16.83, :rxn15) - Rxn(HClg, [H_ion, Cl_ion], 1.971e6, u"mol^2/kg_water^2/atm", 30.20, 19.91, :rxn16) - Rxn(HClg, HCl_ion, 2.5e3, u"mol/kg_water/atm", 30.20, 19.91, :rxn17) - Rxn(H2Oaq, [H_ion, OH_ion], 1.010e-14, u"mol^2/kg_water^2", -22.52, 26.92, :rxn18) - Rxn(Na2SO4s, Na2SO4_aqs, 4.799e-1, u"mol^3/kg_water^3", 0.98, 39.75, :rxn19) - Rxn(NH42SO4s, NH42SO4_aqs, 1.87e0, u"mol^3/kg_water^3", -2.65, 38.57, :rxn20) - Rxn(NH4Cls, [NH3g, HClg], 1.086e-16, u"atm^2", -71.00, 2.40, :rxn21) - Rxn(NaNO3s, NaNO3_aqs, 1.197e1, u"mol^2/kg_water^2", -8.22, 16.01, :rxn22) - Rxn(NaCls, NaCl_aqs, 3.766e1, u"mol^2/kg_water^2", -1.56, 16.90, :rxn23) - Rxn(NaHSO4s, NaHSO4_aqs, 2.413e4, u"mol^2/kg_water^2", 0.79, 14.75, :rxn24) - Rxn(NH4NO3s, [NH3g, HNO3g], 4.199e-17, u"atm^2", -74.375, 6.025, :rxn25) - Rxn(NH4HSO4s, NH4HSO4_aqs, 1.383e0, u"mol^2/kg_water^2", -2.87, 15.83, :rxn26) - Rxn(NH43HSO42s, NH43HSO42_aqs, 2.972e1, u"mol^5/kg_water^5", -5.19, 54.40, :rxn27) + +""" +generate reactions, where slt, g, sld, and i are dictionaries of aqueous salts, gases, solids, and ions, respectively, +and H2Oaq is water. +""" +generate_reactions(slt, g, sld, i, H2Oaq) = [ + # NOTE: Assuming that H_group and C_group are zero when they are left out of Table 2. + Rxn(sld[:CaNO32], slt[:CaNO32], 6.067e5, u"mol^3/kg_water^3", -11.299, 0.0, :rxn1) + Rxn(sld[:CaCl2], slt[:CaCl2], 7.974e11, u"mol^3/kg_water^3", -14.087, 0.0, :rxn2) + Rxn(sld[:CaSO4], slt[:CaSO4], 4.319e-5, u"mol^2/kg_water^2", 0.0, 0.0, :rxn3) + Rxn(sld[:K2SO4], slt[:K2SO4], 1.569e-2, u"mol^3/kg_water^3", -9.589, 45.807, :rxn4) + Rxn(sld[:KHSO4], slt[:KHSO4], 24.016, u"mol^2/kg_water^2", -8.423, 17.964, :rxn5) + Rxn(sld[:KNO3], slt[:KNO3], 0.872, u"mol^2/kg_water^2", 14.075, 19.388, :rxn6) + Rxn(sld[:KCl], slt[:KCl], 8.680, u"mol^2/kg_water^2", -6.167, 19.953, :rxn7) + Rxn(sld[:MgSO4], slt[:MgSO4], 1.079e5, u"mol^2/kg_water^2", 36.798, 0.0, :rxn8) + Rxn(sld[:MgNO32], slt[:MgNO32], 2.507e15, u"mol^3/kg_water^3", -8.754, 0.0, :rxn9) + Rxn(sld[:MgCl2], slt[:MgCl2], 9.557e21, u"mol^3/kg_water^3", -1.347, 0.0, :rxn10) + Rxn(i[:HSO4], [i[:H], i[:SO4]], 1.015e-2, u"mol/kg_water", 8.85, 25.14, :rxn11) + Rxn(g[:NH3], i[:NH3], 5.764e1, u"mol/kg_water/atm", 13.79, -5.39, :rxn12) + Rxn([i[:NH3], H2Oaq], [i[:NH4], i[:OH]], 1.805e-5, u"mol/kg_water", -1.50, 26.92, :rxn13) + Rxn(g[:HNO3], slt[:HNO3], 2.511e6, u"mol^2/kg_water^2/atm", 29.17, 16.83, :rxn14) + Rxn(g[:HNO3], i[:HNO3], 2.1e5, u"mol/kg_water/atm", 29.17, 16.83, :rxn15) + Rxn(g[:HCl], [i[:H], i[:Cl]], 1.971e6, u"mol^2/kg_water^2/atm", 30.20, 19.91, :rxn16) + Rxn(g[:HCl], i[:HCl], 2.5e3, u"mol/kg_water/atm", 30.20, 19.91, :rxn17) + Rxn(H2Oaq, [i[:H], i[:OH]], 1.010e-14, u"mol^2/kg_water^2", -22.52, 26.92, :rxn18) + Rxn(sld[:Na2SO4], slt[:Na2SO4], 4.799e-1, u"mol^3/kg_water^3", 0.98, 39.75, :rxn19) + Rxn(sld[:NH42SO4], slt[:NH42SO4], 1.87e0, u"mol^3/kg_water^3", -2.65, 38.57, :rxn20) + Rxn(sld[:NH4Cl], [g[:NH3], g[:HCl]], 1.086e-16, u"atm^2", -71.00, 2.40, :rxn21) + Rxn(sld[:NaNO3], slt[:NaNO3], 1.197e1, u"mol^2/kg_water^2", -8.22, 16.01, :rxn22) + Rxn(sld[:NaCl], slt[:NaCl], 3.766e1, u"mol^2/kg_water^2", -1.56, 16.90, :rxn23) + Rxn(sld[:NaHSO4], slt[:NaHSO4], 2.413e4, u"mol^2/kg_water^2", 0.79, 14.75, :rxn24) + Rxn(sld[:NH4NO3], [g[:NH3], g[:HNO3]], 4.199e-17, u"atm^2", -74.375, 6.025, :rxn25) + Rxn(sld[:NH4HSO4], slt[:NH4HSO4], 1.383e0, u"mol^2/kg_water^2", -2.87, 15.83, :rxn26) + Rxn(sld[:NH43HSO42], slt[:NH43HSO42], 2.972e1, u"mol^5/kg_water^5", -5.19, 54.40, :rxn27) ] + """ Return a vector of the Species that are products or reactants in the given reactions. """ diff --git a/src/isorropia/solid.jl b/src/isorropia/solid.jl index e6d7e6b8..71a0488d 100644 --- a/src/isorropia/solid.jl +++ b/src/isorropia/solid.jl @@ -13,28 +13,28 @@ From Section 2.2 in Fountoukis and Nenes (2007), the activity of a solid is 1. activity(s::Solid) = 1.0 vars(s::Solid) = [s.m] +terms(s::Solid) = [s.m], [1] function Base.nameof(s::Solid) string(Symbolics.tosymbol(s.m, escape=false)) end -# Generate the solid compounds. -# Each solid compound has an associated MTK variable named -# _s, where is the name of the compound, and -# a Solid struct named s. -all_solids = [] -all_Solids = [] -for s = (:CaNO32, :CaCl2, :CaSO4, :KHSO4, :K2SO4, :KNO3, :KCl, - :MgSO4, :MgNO32, :MgCl2, :NaCl, :NaNO3, :Na2SO4, :NaHSO4, :NH4Cl, - :NH4NO3, :NH42SO4, :NH4HSO4, :NH43HSO42) - varname = Symbol(s, "_s") - solidname = Symbol(s, "s") - description = "Solid $s" - eval(quote - @species $varname = 1e-8 [unit = u"mol/m_air^3", description = $description] - $varname = ParentScope($varname) - push!($all_solids, $varname) - $solidname = $Solid($varname) - push!(all_Solids, $solidname) - end) +""" +Generate the solid compounds. +Each solid compound has an associated MTK variable named +_s, where is the name of the compound, and +a Solid struct named s. +""" +function generate_solids(t) + solids = Dict() + for s = (:CaNO32, :CaCl2, :CaSO4, :KHSO4, :K2SO4, :KNO3, :KCl, + :MgSO4, :MgNO32, :MgCl2, :NaCl, :NaNO3, :Na2SO4, :NaHSO4, :NH4Cl, + :NH4NO3, :NH42SO4, :NH4HSO4, :NH43HSO42) + varname = Symbol(s, "_s") + description = "Solid $s" + v = only(@variables $varname(t) = 1e-8 [unit = u"mol/m_air^3", description = description]) + v = ParentScope(v) + solids[s] = Solid(v) + end + solids end diff --git a/src/isorropia/species.jl b/src/isorropia/species.jl index 353a4b95..93951056 100644 --- a/src/isorropia/species.jl +++ b/src/isorropia/species.jl @@ -1,8 +1,18 @@ """ A species represents a chemical species in the system. + +It should have a `terms` method which returns the variable(s) and +stoichiometry coefficient(s) associated with the species, and +a `vars` method which returns only the variable(s) associated with the species. """ abstract type Species end +"Return the combined terms for a vector of species." +function terms(s::AbstractVector) + tt = terms.(s) + vcat([t[1] for t ∈ tt]...), vcat([t[2] for t ∈ tt]...) +end + # Miscellaneous variables and parameters @parameters T = 293.15 [unit = u"K", description = "Temperature"] diff --git a/src/isorropia/water.jl b/src/isorropia/water.jl index 8e44b4d4..851c1e21 100644 --- a/src/isorropia/water.jl +++ b/src/isorropia/water.jl @@ -4,9 +4,8 @@ end Base.nameof(w::H2O) = "H2O" -H2Oaq = H2O(RH) - vars(w::H2O) = [] +terms(w::H2O) = [], [] # From Equation 15 in Fountoukis and Nenes (2007), the activity of water is # equal to the relative humidity @@ -14,33 +13,35 @@ activity(w::H2O) = w.RH @constants unit_molality=1.0 [unit = u"mol/kg_water", description = "Unit molality"] -# Polynomial fit for molality at given water activity a_w from Table 7 of Fountoukis and Nenes (2007). -m_aw_coeffs = Dict( - CaNO32_aqs => [36.356, -165.66, 447.46, -673.55, 510.91, -155.56, 0] * unit_molality, - CaCl2_aqs => [20.847, -97.599, 273.220, -422.120, 331.160, -105.450, 0] * unit_molality, - KHSO4_aqs => [1.061, -0.101, 1.579 * 10^-2, -1.950 * 10^-3, 9.515 * 10^-5, -1.547 * 10^-6, 0] * unit_molality, - K2SO4_aqs => [1061.51, -4748.97, 8096.16, -6166.16, 1757.47, 0, 0] * unit_molality, - KNO3_aqs => [1.2141 * 10^4, -5.1173 * 10^4, 8.1252 * 10^4, -5.7527 * 10^4, 1.5305 * 10^4, 0, 0] * unit_molality, - KCl_aqs => [179.721, -721.266, 1161.03, -841.479, 221 / 943, 0, 0] * unit_molality, - MgSO4_aqs => [-0.778, 177.74, -719.79, 1174.6, -863.44, 232.31, 0] * unit_molality, - MgNO32_aqs => [12.166, -16.154, 0, 10.886, 0, -6.815, 0] * unit_molality, - MgCl2_aqs => [11.505, -26.518, 34.937, -19.829, 0, 0, 0] * unit_molality, - # TODO(CT): We are missing some salts here. -) + """ Create an equation system for the water content of the aerosol based on Fountoukis and Nenes (2007) Eq. 16. """ -function Water(active_salts) - @variables W = 1.0e-8 [unit = u"kg_water/m_air^3", description = "Aerosol water content"] +function Water(t, active_salts) + # Polynomial fit for molality at given water activity a_w from Table 7 of Fountoukis and Nenes (2007). + m_aw_coeffs = Dict( + :CaNO32_aqs => [36.356, -165.66, 447.46, -673.55, 510.91, -155.56, 0] * unit_molality, + :CaCl2_aqs => [20.847, -97.599, 273.220, -422.120, 331.160, -105.450, 0] * unit_molality, + :KHSO4_aqs => [1.061, -0.101, 1.579 * 10^-2, -1.950 * 10^-3, 9.515 * 10^-5, -1.547 * 10^-6, 0] * unit_molality, + :K2SO4_aqs => [1061.51, -4748.97, 8096.16, -6166.16, 1757.47, 0, 0] * unit_molality, + :KNO3_aqs => [1.2141 * 10^4, -5.1173 * 10^4, 8.1252 * 10^4, -5.7527 * 10^4, 1.5305 * 10^4, 0, 0] * unit_molality, + :KCl_aqs => [179.721, -721.266, 1161.03, -841.479, 221 / 943, 0, 0] * unit_molality, + :MgSO4_aqs => [-0.778, 177.74, -719.79, 1174.6, -863.44, 232.31, 0] * unit_molality, + :MgNO32_aqs => [12.166, -16.154, 0, 10.886, 0, -6.815, 0] * unit_molality, + :MgCl2_aqs => [11.505, -26.518, 34.937, -19.829, 0, 0, 0] * unit_molality, + # TODO(CT): We are missing some salts here. + ) + @variables W(t) = 1.0e-8 [unit = u"kg_water/m_air^3", description = "Aerosol water content"] W = ParentScope(W) w = 0 for s ∈ active_salts - if s ∈ keys(m_aw_coeffs) - m_aw_coeff = m_aw_coeffs[s] - # Water content (Fountoukis and Nenes Eq. 16). - w += (s.cation.m + s.anion.m ) / sum(m_aw_coeff .* RH.^collect(0:6)) + n = Symbol(nameof(s)) + if n ∈ keys(m_aw_coeffs) + m_aw_coeff = m_aw_coeffs[n] + # Water content (Fountoukis and Nenes Eq. 16). + w += (s.cation.m + s.anion.m) / sum(m_aw_coeff .* RH.^collect(0:6)) end end - @named water = NonlinearSystem([W ~ w], [W], []) + @named water = ODESystem([W ~ w], t, [W], []) end \ No newline at end of file diff --git a/test/isorropia/isorropia_test.jl b/test/isorropia/isorropia_test.jl index 5545baba..2147582d 100644 --- a/test/isorropia/isorropia_test.jl +++ b/test/isorropia/isorropia_test.jl @@ -1,5 +1,5 @@ using EarthSciMLBase -using ModelingToolkit, Catalyst, DifferentialEquations, Unitful +using ModelingToolkit, DifferentialEquations, Unitful using Test using Plots @@ -8,64 +8,46 @@ using .ISORROPIA @variables t [unit = u"s", description = "Time"] -model = Isorropia(t); +model = Isorropia(t, :all); sys = structural_simplify(get_mtk(model)) @test all([ModelingToolkit.check_units(eq) for eq in equations(get_mtk(model))]) -@unpack Na_aq, SO4_aq, SO4_g, NH3_aq, NH3_g, NO3_aq, Cl_aq, NaCl_s, MgNO32_s, HNO3_g, - Ca_aq, K_aq, Mg_aq, H_aq, NH4_aq, HCl_g, K2SO4_s, KNO3_s, CaNO32_s, HNO3_g, HNO3_aq, - KHSO4_s, KCl_s, NH4NO3_s, CaSO4_s, CaCl2_s, MgSO4_s, MgCl2_s, NH4HSO4_s, - NH42SO4_s, NH43HSO42_s, NH4Cl_s, NaHSO4_s, Na2SO4_s, NaNO3_s, HSO4_aq, HCl_aq, - RH, metastable, W, I, T, f_CaNO32, rxn1₊k_rev, rxn1₊T₀, - rxn1₊K⁰, rxn1₊H_group, rxn1₊C_group, T₀₂, c_1, I_one = sys - -mw = Dict(Na_aq => 22.989769, SO4_aq => 96.0636, SO4_g => 96.0636, NH3_aq => 17.03052, NH3_g => 17.03052, - NO3_aq => 62.0049, Cl_aq => 35.453, NaCl_s => 58.44, - Ca_aq => 40.078, K_aq => 39.0983, Mg_aq => 24.305, H_aq => 1.00784, NH4_aq => 18.04, HCl_g => 36.46, - K2SO4_s => 174.259, KNO3_s => 101.1032, CaNO32_s => 164.1, HNO3_g => 63.01, HNO3_aq => 63.01, - KHSO4_s => 136.169, KCl_s => 74.5513, NH4NO3_s => 80.043) # g/mol - -# Molecules for checking mass balance -K_molecs = [(1, K_aq), (1, KHSO4_s), (2, K2SO4_s), (1, KNO3_s), (1, KCl_s)] -Ca_molecs = [(1, Ca_aq), (1, CaSO4_s), (1, CaNO32_s), (1, CaCl2_s)] -Mg_molecs = [(1, Mg_aq), (1, MgSO4_s), (1, MgNO32_s), (1, MgCl2_s)] -NH_molecs = [(1, NH4_aq), (1, NH3_aq), (1, NH3_g), (1, NH4HSO4_s), - (2, NH42SO4_s), (3, NH43HSO42_s), (1, NH4Cl_s), (1, NH4NO3_s)] -Na_molecs = [(1, Na_aq), (1, NaHSO4_s), (2, Na2SO4_s), (1, NaCl_s), (1, NaNO3_s)] -SO4_molecs = [(1, SO4_aq), (1, HSO4_aq), (1, SO4_g), - (1, KHSO4_s), (1, NaHSO4_s), (1, NH4HSO4_s), (1, CaSO4_s), (1, Na2SO4_s), (1, NH42SO4_s), - (2, NH43HSO42_s), (1, K2SO4_s), (1, MgSO4_s)] -NO3_molecs = [(1, NO3_aq), (1, HNO3_aq), (1, HNO3_g), (1, NH4NO3_s), (1, NaNO3_s)] -Cl_molecs = [(1, Cl_aq), (1, HCl_aq), (1, HCl_g), (1, NH4Cl_s), - (1, NaCl_s), (2, CaCl2_s), (1, KCl_s), (2, MgCl2_s)] -H_molecs = [(1, H_aq), (1, HNO3_g), (1, HCl_g)] - - defaults = ModelingToolkit.get_defaults(sys) -u₀ = Dict{Any,Float64}([s => 1.0e-15 for s ∈ states(sys)]) -u₀[CaNO32_s] = 0.8 / 1e6 / mw[CaNO32_s] -u₀[Ca_aq] = 0.8 / 1e6 / mw[Ca_aq] -u₀[NO3_aq] = 10.0 / 1e6 / mw[NO3_aq] +u₀ = Dict{Any,Float64}([s => 1.0e-8 for s ∈ states(sys)]) +#u₀[CaNO32_s] = 0.8 / 1e6 / mw[CaNO32_s] +#u₀[Ca_aq] = 0.8 / 1e6 / mw[Ca_aq] +#u₀[NO3_aq] = 10.0 / 1e6 / mw[NO3_aq] p = Dict{Any,Float64}([p => defaults[p] for p ∈ parameters(sys)]) p[RH] = 0.30 -prob = ODEProblem(sys, u₀, (0.0, 100.0), p) +prob = ODEProblem(sys, u₀, (0.0, 3e-8), p) # Need low tolerance for mass balance checks to pass. @time sol = solve(prob, Rosenbrock23(), abstol=1e-12, reltol=1e-12) -# eq = [] -# for i in eachindex(sol.u) -# u = [x => sol[x][i] for x in [I, W, states(sys)...]] -# push!(u, T => sol[T]) -# γ_p = substitute(ModelingToolkit.subs_constants(γ(CaNO32_aqs)), u) -# a_r = substitute(ModelingToolkit.subs_constants(activity(CaNO32s)), u) -# a_p = substitute(ModelingToolkit.subs_constants(activity(CaNO32_aqs)), u) -# keq = sol[rxn1.sys.K_eq][i] -# @info :a_p => a_p, :γ_p => γ_p, :Ca_aq => sol[Ca_aq][i], :a_r => a_r, :keq => keq, :ratio => a_p / a_r / keq -# push!(eq, a_p / a_r / keq) -# end +begin + p1 = plot(title="Solids", xscale=:log10) + 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) + 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) + 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 + r = Symbol(:rxn, i) + y = eval(:(sol[sys.$r.rate])) + plot!(sol.t[2:end], y[2:end], label="$r.rate") + end + plot(p1, p2, p3, p4, size=(1000, 800)) +end + plot( plot(sol[t], sol[f_CaNO32],