Correctness of henry's law, phase transfer algorithms

[K20shores]

#887 brought in support for two reactions, but we have questions on both

1. Are reversible reactions assumed to have a constant rate, and can therefore be specified entirely in the configuration file, or should they allowed to be updated throughout a simulation?
2. Is the henry's law equatio we used correct? It was for diprotic acids and, at least to me, seems to be some kind of taylor series expansion. It does not match what camp does, and varies from the equations from seinfeld and pandis. What equation should be used

For reference, this is the seinfeld and pandis equation, and what camp uses

$$ H_A(T_2)=H_A(T_1)~\exp\left(\frac{\Delta H_A}{R}\left(\frac{1}{T_1}-\frac{1}{T_2}\right)\right) $$

[K20shores]

@mattldawson would you happen to know any of these answers?

[mattldawson]

For the first question, the answer depends on whether the equilibrium constant is temperature and/or pressure dependent. In general:

K_eq = k_f / k_r

where K_eq is the equilibrium rate constant (in this case the Henry's Law constant), k_f is the forward (condensation) rate constant, and k_r is the reverse (evaporation) rate constant. CAMP uses a temperature-dependent Henry's Law Constant with this equation:

K_eq(T) = HLC(298K) * e^(C * ( 1/T - 1/298K ))

where HLC(298K) is the Henry's Law constant at 298K, C is a constant that captures the temperature dependence of a particular species' HLC, and T is temperature (K). (This is essentially the Seinfeld and Pandis equation with T1 = 298K and C = -dH_A/R)

Also, the forward (condensation) rate constant k_f can be calculated based on the current conditions (we can discuss this more if needed). So, the only unknown for a given set of environmental/aerosol conditions is k_r. The first equation can be rearranged to get:

k_r = k_f(T, r_eff, D_g) / K_eq(T)

where r_eff is the effective radius of the particle and D_g is the diffusion coefficient of the gas-phase species. (Note that these equations are for condensation/evaporation from a single particle. For collections of particles (e.g., modes) the number concentration will also be involved, but how depends on how the aerosol state is defined.)

So, long story short, if you want to do what CAMP does, you can include HLC(298K) and C to define a species' Henry's Law constant, and D_g for its diffusion coefficient, and everything else could be calculated from the current model state at each time step.

For the second question, CAMP uses actual Henry's Law Constants, whereas Henry's Law Constants that are labelled as being for mono- or di-protic acids, or for bases are actually "Effective" Henry's Law Constants (Mary can tell you more about these). In essence, the Effective Henry's Law Constant, say for a di-protic acid, combines the actual Henry's Law Constant with the first and second acid dissociation reactions. This treats the (typically very fast) acid-dissociation steps as being at quasi-steady-state. The advantage is that it makes the system less stiff. The disadvantage is that it requires a different treatment for each type of acid/base, and it requires that the proton concentration [H+] be known (or estimated). (This is what is done in the iterative loop in the cloud chemistry code in CAM.)

If David and Jiwon are still planning on implementing the DAE solver, my suggestion would be to use actual HLC equations (like in CAMP) with separate acid-base dissociation steps. You can include any/all of the acid-base dissociation reactions as algebraic constraints in the DAE solver using their equilibrium constants. (David can describe what this would look like in more detail). This way, you don't need multiple HLC equations, you don't need to estimate the [H+], and you can solve everything together with one solver (no process-splitting; no nested iterative loops).

Hope this helps! Let me know if anything needs more clarification or seems wrong (I'm mostly going from memory).

[K20shores]

@mattldawson ah, thank you!

@boulderdaze given this information, it appears we will need to modify your original PR

[boulderdaze]

Thank you for the detailed explanations, @mattldawson. I've created two issues based on his suggestions:

• Temperature-Dependent Henry’s Law and Condensation, Evaporation Rates #892,
• Using HLC Equations with separate acid-base dissociation steps #893

Tagging @dwfncar for visibility and input

[boulderdaze]

Hi @mattldawson, can I ask a question?

This is the current implementation for carbonic acid formation based on the MIAM README example. The species listed in CondensedSpecies() are CO32- + 2H+, which correspond to the second dissociation step. This suggests that the first dissociation process is not being accounted for, since its species are not listed.

  // Build an Effective Henry's Law PhaseTransferProcess
  // (will be for a diprotic acid: CO2(g) <-> H2CO3(aq) <-> HCO3- + H+ <-> CO32- + 2H+)
  // K_a1 = first acid dissociation constant
  // K_a2 = second acid dissociation constant
  // K_H = A * exp(C / T) = Henry's Law constant
  auto hlc_co2_parms = HenrysLawCoefficientParameters{ .A_ = 32.4, .C_ = -3.2e-3, .K_a1_ = 1.0e-5, .K_a2_ = 2.0e-5 };

  Process co2_phase_transfer = PhaseTransferProcessBuilder()
                               .SetGasSpecies(gas_phase, co2 )
                               .SetCondensedSpecies(aqueous_phase, { Yield(hp, 2.0), Yield(co32m) })
                               .SetSolvent(aqueous_phase, h2o )
                               .SetTransferCoefficient(HenrysLawCoefficient(hlc_co2_parms))
                               .Build();

When you suggest "to use actual HLC equations (like in CAMP) with separate acid-base dissociation steps", this is what I had in mind. Does that look like what you were suggesting?

 Process co2_phase_transfer_first_dissociation = 
   PhaseTransferProcessBuilder()
   ...
   .SetCondensedSpecies(aqueous_phase, { Yield(H2CO3), Yield(HCO3-), Yield(H+) });

 Process co2_phase_transfer_second_dissociation = 
   PhaseTransferProcessBuilder()
   ...
   .SetCondensedSpecies(aqueous_phase, { Yield(HCO3-), Yield(H+), Yield(CO32-), Yield(2H+) });

[mattldawson]

Similar, I think, but maybe a couple minor modifications:

1. There should be one phase transfer process (CO(g) <-> H2CO3(aq)) that uses the HLC to calculate the rates of condensation/evaporation of the undissociated species. (In this particular case it's a little confusing because there is an additional assumption about some fast oxidation reactions that convert CO2(aq) into H2CO3(aq), but generally these reactions should be something like foo(g) <-> foo(aq)). Also, there aren't yields for these processes - they're 1:1 gas-species:condensed-species.
2. The acid dissociation reactions are not actually related to phase transfer. They are just normal reversible reactions that are entirely within an aqueous phase.

So, maybe something like this:

  // Henry's Law partitioning of CO2, assuming rapid oxidation in the aqueous phase to form H2CO3(aq)
  auto hlc_co2_parms = HenrysLawCoefficientParameters{ .A_ = 32.4, .C_ = -3.2e-3, .K_a1_ = 1.0e-5, .K_a2_ = 2.0e-5 };
  Process co2_phase_transfer = PhaseTransferProcessBuilder()
                               .SetGasSpecies(gas_phase, co2 )
                               .SetCondensedPhase(aqueous_phase)
                               .SetCondensedSpecies(h2co3)
                               .SetSolvent(h2o)
                               .SetTransferCoefficient(HenrysLawCoefficient(hlc_co2_parms))
                               .Build();
                               
  // First reversible acid-dissociation reaction `H2CO3(aq) <-> HCO3-(aq) + H+(aq)`
  auto h2co3_hco3m_hp = DissolvedReversibleProcessBuilder()
                                .SetPhase(aqueous_phase)
                                .SetReactants({ h2co3 })
                                .SetProducts({ Yield(hco3m, 1.0), Yield( hp, 1.0) })
                                .SetSolvent(h2o)
                                .SetForwardRateConstant(ArrheniusRateConstant({ .A_ = 123.4, .B_ = 123.4 , <etc>})
                                .SetEquilibriumConstant(<might need a new algorithm, or it could be handled by Arrhenius equation>)
                                .Build()

  // Second reversible acid-dissociation reaction `HCO3-(aq) <-> CO3--(aq) + H+(aq)`
  auto hco3m_co32m_hp = DissolvedReversibleProcessBuilder()
                                .SetPhase(aqueous_phase)
                                .SetReactants({ hco3m })
                                .SetProducts({ Yield(co32m, 1.0), Yield( hp, 1.0) })
                                .SetSolvent(h2o)
                                .SetForwardRateConstant(ArrheniusRateConstant({ .A_ = 123.4, .B_ = 123.4 , <etc>})
                                .SetEquilibriumConstant(<might need a new algorithm, or it could be handled by Arrhenius equation>)
                                .Build()

Let me know if that looks reasonable. The forward rate constant and the equilibrium constant (sometimes called Ka1, Ka2) should be available for the two dissociation steps. (Let me know if you want me to track them down)

[boulderdaze]

Yes, this makes sense to me. So the plan is to create a new builder that handles the reversible dissolution process. Thank you for sharing the concrete APIs.

[mattldawson]

The DissolvedProcessBuilder() could accept any two of the following:
SetForwardRateConstant()
SetReverseRateConstant()
SetEquilibriumConstant()

Two of the three are sufficient to calculate the rates

[dwfncar]

See our open PR (1/3) for the DAE constraints. PR #900 (Add constraint base classes for DAE support) can achieve mathematically equivalent results to Effective Henry's Law Constants, but through separate equilibrium constraints solved together by the DAE solver.

Concrete example using EquilibriumConstraint:

For the CO2/carbonic acid system, instead of using an Effective HLC formula that requires knowing [H+] in advance:

H* = H × (1 + K_a1/[H+] + K_a1×K_a2/[H+]²)

Express the acid dissociation steps as separate constraints:

// First dissociation: H2CO3 <-> HCO3- + H+
// Equilibrium: K_a1 = [HCO3-][H+] / [H2CO3]
EquilibriumConstraint("Ka1",
    {Yield(H2CO3, 1.0)},                    // reactants
    {Yield(HCO3m, 1.0), Yield(Hp, 1.0)},    // products
    K_a1)

// Second dissociation: HCO3- <-> CO32- + H+
// Equilibrium: K_a2 = [CO32-][H+] / [HCO3-]
EquilibriumConstraint("Ka2",
    {Yield(HCO3m, 1.0)},                    // reactants
    {Yield(CO32m, 1.0), Yield(Hp, 1.0)},    // products
    K_a2)

When both constraints are satisfied (residual = 0), total dissolved carbon is:

[H2CO3] + [HCO3-] + [CO32-] = H × P_CO2 × (1 + K_a1/[H+] + K_a1×K_a2/[H+]²)

Equivalent to an Effective HLC:

Effective HLC	Separate Constraints
Different formulas for mono/di/triprotic acids	Same EquilibriumConstraint class for all
Must estimate [H+] beforehand	[H+] solved as part of the system
Process splitting / nested loops	Single unified DAE solve
Hard to extend to bases or mixed systems	Just add more constraints

I think what we are doing is consistent with @mattldawson's recommendation to use actual HLC equations with separate acid-base dissociation steps as algebraic constraints. The EquilibriumConstraint class in PR #900 handles the dissociation steps; still need a phase transfer process for CO2(g) <-> H2CO3(aq) using the actual (non-effective) HLC, which could be kinetic or treated as an additional equilibrium constraint if instantaneous.

[boulderdaze]

As Matt pointed out, reactants probably shouldn't be defined asYield type. Would it make sense to rename Yield to something like StoichSpecies (short for Stoichiometric Species)?

// First dissociation: H2CO3 <-> HCO3- + H+
// Equilibrium: K_a1 = [HCO3-][H+] / [H2CO3]
EquilibriumConstraint("Ka1",
    {Yield(H2CO3, 1.0)},                    // reactants    < ------ here
    {Yield(HCO3m, 1.0), Yield(Hp, 1.0)},    // products
    K_a1)

struct StoichSpecies 
{
  Species species_;
  double coefficient_;
};

[K20shores]

I'm fine with that. We could also consider Stoichiometry. I'm indifferent on the name

[boulderdaze]

Hi @mattldawson, looking at your implementation of the stub aerosol model, the model would know the constraint conditions, since it determines the species within a given mode. For example, in the case of an equilibrium constraint, it could be expressed as

$$K_eq = [H+] [model.mode.aqueous.baz] / ([gas.foo] )$$

Would it make sense for ExternalModel to own Constraint object? How do you envision the constraint being defined and configured?

This is the current version of the solver build (still in progress), where the constraints are separate from the aerosol models. I'm wondering whether the model itself should define and own the constraint. If that is the case, the API will likely need to be revised.

  // Current
  auto solver = builder.SetSystem(system)
                       .AddExternalModelProcesses(aerosol_1)
                       .AddExternalModelProcesses(aerosol_2)
                       .SetConstraints(constraints)
                       .SetIgnoreUnusedSpecies(true)
                       .Build();

I'd also like to hear what you @K20shores, @dwfncar think.

[mattldawson]

I haven't looked in detail at the DAE PRs, but in general I agree with the approach to structure the ConstraintSet similar to the ProcessSet. Constraints aren't unique to external models, there could also be constraints just for gas-phase systems. But, just like for processes, external models may introduce their own constraints. For the initial implementation of the DAE solver, I wouldn't worry too much about the external models - just organize constraints as consistently with processes as possible. After the DAE solver is in place, we can update the external_model.hpp classes to also allow constraints to be provided from external models. (This is when we would add the acid-base/Henry's Law type constraints.)

Hope this helps!