From d4cbd4812009a54d17b0172ffb8e4a78d38c0c2b Mon Sep 17 00:00:00 2001 From: jkgolla Date: Wed, 2 Jul 2025 11:25:09 -0700 Subject: [PATCH 1/2] Update carbonate system reaction constants --- src/reactions/geochemistry/Carbonate.hpp | 53 +++++++++++------------- 1 file changed, 25 insertions(+), 28 deletions(-) diff --git a/src/reactions/geochemistry/Carbonate.hpp b/src/reactions/geochemistry/Carbonate.hpp index 137903c..c7aea78 100644 --- a/src/reactions/geochemistry/Carbonate.hpp +++ b/src/reactions/geochemistry/Carbonate.hpp @@ -41,35 +41,34 @@ constexpr CArrayWrapper stoichMatrix = }; constexpr CArrayWrapper stoichMatrixNosolid = - { // OH- CO2 CO3-2 H2CO3 CaHCO3+ CaSO4 CaCl+ CaCl2 MgSO4 NaSO4- H+ HCO3- Ca+2 SO4-2 Cl- Mg+2 Na+ - { -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0 }, // OH- + H+ = H2O - { 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0 }, // CO2 + H2O = H+ + HCO3- - { 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 0, 0, 0 }, // CO3-2 + H+ = HCO3- - { 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0 }, // H2CO3 = H+ + HCO3- - { 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0 }, // CaHCO3+ = Ca+2 + HCO3- - { 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0 }, // CaSO4 = Ca+2 + SO4-2 - { 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0 }, // CaCl+ = Ca+2 + Cl- - { 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 2, 0, 0 }, // CaCl2 = Ca+2 + 2Cl- - { 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 1, 0 }, // MgSO4 = Mg+2 + SO4-2 - { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 1 }, // NaSO4- = Na+ + SO4-2 - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 1, 0, 0, 0, 0 } // CaCO3(s) + H+ = Ca+2 + HCO3- (kinetic) - // { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 1, 0, 0, 0, 0 } // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic) + { // OH- CO2 CO3-2 CaHCO3+ CaSO4 CaCl+ CaCl2 MgSO4 NaSO4- H+ HCO3- Ca+2 SO4-2 Cl- Mg+2 Na+ + { -1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0 }, // OH- + H+ = H2O + { 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0 }, // CO2 + H2O = H+ + HCO3- + { 0, 0, -1, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 0, 0, 0 }, // CO3-2 + H+ = HCO3- + { 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0 }, // CaHCO3+ = Ca+2 + HCO3- + { 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0 }, // CaSO4 = Ca+2 + SO4-2 + { 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0 }, // CaCl+ = Ca+2 + Cl- + { 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 2, 0, 0 }, // CaCl2 = Ca+2 + 2Cl- + { 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 1, 0 }, // MgSO4 = Mg+2 + SO4-2 + { 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 1 }, // NaSO4- = Na+ + SO4-2 + { 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 1, 0, 0, 0, 0 } // CaCO3(s) + H+ = Ca+2 + HCO3- (kinetic) + // { 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 1, 0, 0, 0, 0 } // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic) }; // C^{n+1} - C^n - r( C^{n+1} ) * dt = 0 +// thermodynamic constants derived from 'llnl.tdat' used by Geochemists' Workbench (originally from EQ36) constexpr CArrayWrapper equilibriumConstants = { - 9.77E+13, // OH- + H+ = H2O - 4.37E-07, // CO2 + H2O = H+ + HCO3- - 2.14E+10, // CO3-2 + H+ = HCO3- - 1.70E-04, // H2CO3 = H+ + HCO3- - 8.13E-02, // CaHCO3+ = Ca+2 + HCO3- - 6.92E-03, // CaSO4 = Ca+2 + SO4-2 - 4.68E+00, // CaCl+ = Ca+2 + Cl- + 9.89E+13, // OH- + H+ = H2O + 4.42E-07, // CO2 + H2O = H+ + HCO3- + 2.21E+10, // CO3-2 + H+ = HCO3- + 6.00E-02, // CaHCO3+ = Ca+2 + HCO3- + 4.79E-03, // CaSO4 = Ca+2 + SO4-2 + 2.00E-01, // CaCl+ = Ca+2 + Cl- 3.98E+00, // CaCl2 = Ca+2 + 2Cl- - 3.72E-03, // MgSO4 = Mg+2 + SO4-2 - 1.51E-01, // NaSO4- = Na+ + SO4-2 - 1.17E+07 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) + 5.92E-03, // MgSO4 = Mg+2 + SO4-2 + 2.02E-01, // NaSO4- = Na+ + SO4-2 + 5.16E+01 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) // 1 }; // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic) @@ -77,15 +76,14 @@ constexpr CArrayWrapper forwardRates = { 1.4e11, // OH- + H+ = H2O 0.039, // CO2 + H2O = H+ + HCO3- - 1.0e10, // CO3-2 + H+ = HCO3- - 0.57, // H2CO3 = H+ + HCO3- + 1.0e10, // CO3-2 + H+ = HCO3- 1.5e6, // CaHCO3+ = Ca+2 + HCO3- 1.0e5, // CaSO4 = Ca+2 + SO4-2 1.0e8, // CaCl+ = Ca+2 + Cl- 1.0e7, // CaCl2 = Ca+2 + 2Cl- 1.0e5, // MgSO4 = Mg+2 + SO4-2 1.0e7, // NaSO4- = Na+ + SO4-2 - 1.0e5 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) + 1.55E-06 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic) // 1 }; // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic) @@ -94,14 +92,13 @@ constexpr CArrayWrapper reverseRates = { 1.43E-03, // OH- + H+ = H2O 8.92E+04, // CO2 + H2O = H+ + HCO3- 4.67E-01, // CO3-2 + H+ = HCO3- - 3.35E+03, // H2CO3 = H+ + HCO3- 1.85E+07, // CaHCO3+ = Ca+2 + HCO3- 1.45E+07, // CaSO4 = Ca+2 + SO4-2 2.14E+07, // CaCl+ = Ca+2 + Cl- 2.51E+06, // CaCl2 = Ca+2 + 2Cl- 2.69E+07, // MgSO4 = Mg+2 + SO4-2 6.62E+07, // NaSO4- = Na+ + SO4-2 - 8.55E-03 // CaCO3 + H+ = Ca+2 + HCO3- + 3.00E-08 // CaCO3 + H+ = Ca+2 + HCO3- // 1 // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic) }; } From 453dbc683ae60c95a2eb5a12be29dafe8ec0828a Mon Sep 17 00:00:00 2001 From: jkgolla Date: Wed, 2 Jul 2025 11:29:02 -0700 Subject: [PATCH 2/2] Update carbonate reaction system --- src/reactions/geochemistry/Carbonate.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/reactions/geochemistry/Carbonate.hpp b/src/reactions/geochemistry/Carbonate.hpp index c7aea78..6a86d9c 100644 --- a/src/reactions/geochemistry/Carbonate.hpp +++ b/src/reactions/geochemistry/Carbonate.hpp @@ -40,7 +40,7 @@ constexpr CArrayWrapper stoichMatrix = // { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 1, 0, 0, 0, 0 } // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic) }; -constexpr CArrayWrapper stoichMatrixNosolid = +constexpr CArrayWrapper stoichMatrixNosolid = { // OH- CO2 CO3-2 CaHCO3+ CaSO4 CaCl+ CaCl2 MgSO4 NaSO4- H+ HCO3- Ca+2 SO4-2 Cl- Mg+2 Na+ { -1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0 }, // OH- + H+ = H2O { 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0 }, // CO2 + H2O = H+ + HCO3- @@ -104,7 +104,7 @@ constexpr CArrayWrapper reverseRates = } using carbonateSystemAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 18, 11, 0 >; using carbonateSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 18, 11, 11 >; - using carbonateSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 17, 11, 10 >; + using carbonateSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 16, 11, 10 >; constexpr carbonateSystemAllKineticType carbonateSystemAllKinetic( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates ); constexpr carbonateSystemAllEquilibriumType carbonateSystemAllEquilibrium( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates );