Normalizing chemical reactions

CADET Development
1

Background

In CADET, we can model the pH by numerically solving (de)protonation reactions:

\ce{HA <=>[k_{fwd}][k_{bwd}] A^- + H^+}\\

with

\frac{k_{fwd}}{k_{bwd}} = k_{eq} = k_{a} = 10^{-pK_a} = \frac{\ce{[A^-][H^+]}}{\ce{[HA]}}

The pH is then defined as:

pH = -log_{10}(\ce{[H^+]})

Typically, pK_a values range between 2 for strong acids and 12 for weak acids (assuming molar concentrations).

Issues

  • Proton concentrations at high pH are (by definition) very low.
    => Tight tolerances are required to solve the problem.
  • The ratio of forward and backward reaction is very high for high pK_a values.
    => Different time scales.
  • CADET cannot solve reactions in rapid equilibrium.
    => Forward and backward coefficient need to be scaled up which further increases the stiffness of the problem.

Idea: Normalization/Scaling of reaction rates and proton concentration

My question is now, is it possible to somehow scale/normalize the reacton rates and proton concentration to reduce stiffness or are the non-linearities of the system such that this is not easily possible?

For example, could we introduce some factor \alpha that scales the forward rate, as well as some coefficients that correspondingly scale the stoichiometry of the product side.

\ce{HA <=>[\alpha \cdot k_{fwd}][k_{bwd}] \beta A^- + \gamma H^+}\\

My initial thought was to set \beta = 1, \gamma = \alpha or \beta = \sqrt{\alpha}, \gamma = \sqrt{\alpha}, but so far I haven’t been able to yield the same results. Any help is highly appreciated!

Requirements:

  • Coupling between reactions via proton concentration must be conserved.
  • Adsorption processes must be conserved (e.g. by also scaling the equilibria).

Of course this would also imply some other aspects in pre- and postprocessing but for now I’m just interested whether this still computes the correct solution.

2

I still don’t get this normalization idea. Can you elaborate a bit more?

For the given example, the equations look like

\frac{\mathrm{d}}{\mathrm{d}t} \begin{pmatrix} c_{HA} \\ c_{A^-} \\ c_{H^+} \end{pmatrix} = \begin{pmatrix} -1 \\ 1 \\ 1 \end{pmatrix} r( c_{HA}, c_{A^-}, c_{H^+}),

where the rate function is given by

r( c_{HA}, c_{A^-}, c_{H^+}) = k_{\mathrm{fwd}}\, c_{HA} - k_{\mathrm{bwd}}\, c_{A^-}\, c_{H^+} = k_{\mathrm{bwd}} \left( \frac{k_{\mathrm{fwd}}}{k_{\mathrm{bwd}}}\, c_{HA} - c_{A^-}\, c_{H^+} \right)

I don’t see how the scaling factors \beta and \gamma could work.

3

I would not be surprised if it didn’t work. :sweat_smile:

The question is whether it is somehow possible to introduce some parameter transform by scaling the reaction rates and stoichiometric coefficients such that the problem becomes easier to solve.