Skip to content

Commit 2f94ce2

Browse files
frankfeifanjkgolla
frankfeifan
and
jkgolla
authored
fix: corrected rate parameters in carbonate system (#9)
* Update carbonate system reaction constants * Resolve redundant aqueous carbonate secondary species and update affected unit tests Co-authored-by: jkgolla <[email protected]>
1 parent f20e43a commit 2f94ce2

File tree

5 files changed

+135
-153
lines changed

5 files changed

+135
-153
lines changed

‎src/reactions/geochemistry/Carbonate.hpp‎

Lines changed: 62 additions & 72 deletions
Original file line numberDiff line numberDiff line change
@@ -24,92 +24,80 @@ namespace geochemistry
2424
namespace carbonate
2525
{
2626

27-
constexpr CArrayWrapper<double, 11, 18> stoichMatrix =
28-
{ // OH- CO2 CO3-2 H2CO3 CaHCO3+ CaSO4 CaCl+ CaCl2 MgSO4 NaSO4- CaCO3 H+ HCO3- Ca+2 SO4-2 Cl- Mg+2 Na+
29-
{ -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0 }, // OH- + H+ = H2O
30-
{ 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0 }, // CO2 + H2O = H+ + HCO3-
31-
{ 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 0, 0, 0 }, // CO3-2 + H+ = HCO3-
32-
{ 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0 }, // H2CO3 = H+ + HCO3-
33-
{ 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0 }, // CaHCO3+ = Ca+2 + HCO3-
34-
{ 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0 }, // CaSO4 = Ca+2 + SO4-2
35-
{ 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0 }, // CaCl+ = Ca+2 + Cl-
36-
{ 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 2, 0, 0 }, // CaCl2 = Ca+2 + 2Cl-
37-
{ 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 1, 0 }, // MgSO4 = Mg+2 + SO4-2
38-
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 1 }, // NaSO4- = Na+ + SO4-2
39-
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, 1, 1, 0, 0, 0, 0 } // CaCO3(s) + H+ = Ca+2 + HCO3- (kinetic)
40-
// { 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)
27+
constexpr CArrayWrapper<double, 10, 17> stoichMatrix =
28+
{ // OH- CO2 CO3-2 CaHCO3+ CaSO4 CaCl+ CaCl2 MgSO4 NaSO4- CaCO3 H+ HCO3- Ca+2 SO4-2 Cl- Mg+2 Na+
29+
{ -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0 }, // OH- + H+ = H2O
30+
{ 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0 }, // CO2 + H2O = H+ + HCO3-
31+
{ 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 0, 0, 0 }, // CO3-2 + H+ = HCO3-
32+
{ 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0 }, // CaHCO3+ = Ca+2 + HCO3-
33+
{ 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0 }, // CaSO4 = Ca+2 + SO4-2
34+
{ 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0 }, // CaCl+ = Ca+2 + Cl-
35+
{ 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 2, 0, 0 }, // CaCl2 = Ca+2 + 2Cl-
36+
{ 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 1, 0 }, // MgSO4 = Mg+2 + SO4-2
37+
{ 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 1 }, // NaSO4- = Na+ + SO4-2
38+
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, 1, 1, 0, 0, 0, 0 } // CaCO3(s) + H+ = Ca+2 + HCO3- (kinetic)
4139
};
4240

43-
constexpr CArrayWrapper<double, 11, 17> stoichMatrixNosolid =
44-
{ // OH- CO2 CO3-2 H2CO3 CaHCO3+ CaSO4 CaCl+ CaCl2 MgSO4 NaSO4- H+ HCO3- Ca+2 SO4-2 Cl- Mg+2 Na+
45-
{ -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0 }, // OH- + H+ = H2O
46-
{ 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0 }, // CO2 + H2O = H+ + HCO3-
47-
{ 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 0, 0, 0 }, // CO3-2 + H+ = HCO3-
48-
{ 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0 }, // H2CO3 = H+ + HCO3-
49-
{ 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0 }, // CaHCO3+ = Ca+2 + HCO3-
50-
{ 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0 }, // CaSO4 = Ca+2 + SO4-2
51-
{ 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0 }, // CaCl+ = Ca+2 + Cl-
52-
{ 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 2, 0, 0 }, // CaCl2 = Ca+2 + 2Cl-
53-
{ 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 1, 0 }, // MgSO4 = Mg+2 + SO4-2
54-
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 1 }, // NaSO4- = Na+ + SO4-2
55-
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 1, 0, 0, 0, 0 } // CaCO3(s) + H+ = Ca+2 + HCO3- (kinetic)
56-
// { 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)
41+
constexpr CArrayWrapper<double, 10, 16> stoichMatrixNosolid =
42+
{ // OH- CO2 CO3-2 CaHCO3+ CaSO4 CaCl+ CaCl2 MgSO4 NaSO4- H+ HCO3- Ca+2 SO4-2 Cl- Mg+2 Na+
43+
{ -1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0 }, // OH- + H+ = H2O
44+
{ 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0 }, // CO2 + H2O = H+ + HCO3-
45+
{ 0, 0, -1, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 0, 0, 0 }, // CO3-2 + H+ = HCO3-
46+
{ 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0 }, // CaHCO3+ = Ca+2 + HCO3-
47+
{ 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0 }, // CaSO4 = Ca+2 + SO4-2
48+
{ 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0 }, // CaCl+ = Ca+2 + Cl-
49+
{ 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 2, 0, 0 }, // CaCl2 = Ca+2 + 2Cl-
50+
{ 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 1, 0 }, // MgSO4 = Mg+2 + SO4-2
51+
{ 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 1 }, // NaSO4- = Na+ + SO4-2
52+
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 1, 0, 0, 0, 0 } // CaCO3(s) + H+ = Ca+2 + HCO3- (kinetic)
5753
};
5854

59-
// C^{n+1} - C^n - r( C^{n+1} ) * dt = 0
60-
constexpr CArrayWrapper<double, 11> equilibriumConstants =
55+
// thermodynamic constants derived from 'llnl.tdat' used by Geochemists' Workbench (originally from EQ36)
56+
constexpr CArrayWrapper<double, 10> equilibriumConstants =
6157
{
62-
9.77E+13, // OH- + H+ = H2O
63-
4.37E-07, // CO2 + H2O = H+ + HCO3-
64-
2.14E+10, // CO3-2 + H+ = HCO3-
65-
1.70E-04, // H2CO3 = H+ + HCO3-
66-
8.13E-02, // CaHCO3+ = Ca+2 + HCO3-
67-
6.92E-03, // CaSO4 = Ca+2 + SO4-2
68-
4.68E+00, // CaCl+ = Ca+2 + Cl-
69-
3.98E+00, // CaCl2 = Ca+2 + 2Cl-
70-
3.72E-03, // MgSO4 = Mg+2 + SO4-2
71-
1.51E-01, // NaSO4- = Na+ + SO4-2
72-
1.17E+07 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic)
73-
// 1
74-
}; // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic)
58+
9.89E+13, // OH- + H+ = H2O
59+
4.42E-07, // CO2 + H2O = H+ + HCO3-
60+
2.21E+10, // CO3-2 + H+ = HCO3-
61+
6.00E-02, // CaHCO3+ = Ca+2 + HCO3-
62+
4.79E-03, // CaSO4 = Ca+2 + SO4-2
63+
2.00E-01, // CaCl+ = Ca+2 + Cl-
64+
3.98E+00, // CaCl2 = Ca+2 + 2Cl-
65+
5.92E-03, // MgSO4 = Mg+2 + SO4-2
66+
2.02E-01, // NaSO4- = Na+ + SO4-2
67+
5.16E+01 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic)
68+
};
7569

76-
constexpr CArrayWrapper<double, 11> forwardRates =
70+
constexpr CArrayWrapper<double, 10> forwardRates =
7771
{
7872
1.4e11, // OH- + H+ = H2O
7973
0.039, // CO2 + H2O = H+ + HCO3-
80-
1.0e10, // CO3-2 + H+ = HCO3-
81-
0.57, // H2CO3 = H+ + HCO3-
74+
1.0e10, // CO3-2 + H+ = HCO3-
8275
1.5e6, // CaHCO3+ = Ca+2 + HCO3-
8376
1.0e5, // CaSO4 = Ca+2 + SO4-2
8477
1.0e8, // CaCl+ = Ca+2 + Cl-
8578
1.0e7, // CaCl2 = Ca+2 + 2Cl-
8679
1.0e5, // MgSO4 = Mg+2 + SO4-2
8780
1.0e7, // NaSO4- = Na+ + SO4-2
88-
1.0e5 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic)
89-
90-
// 1
91-
}; // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic)
81+
1.55E-06 // CaCO3 + H+ = Ca+2 + HCO3- (kinetic)
82+
};
9283

93-
constexpr CArrayWrapper<double, 11> reverseRates =
94-
{ 1.43E-03, // OH- + H+ = H2O
95-
8.92E+04, // CO2 + H2O = H+ + HCO3-
96-
4.67E-01, // CO3-2 + H+ = HCO3-
97-
3.35E+03, // H2CO3 = H+ + HCO3-
98-
1.85E+07, // CaHCO3+ = Ca+2 + HCO3-
99-
1.45E+07, // CaSO4 = Ca+2 + SO4-2
100-
2.14E+07, // CaCl+ = Ca+2 + Cl-
101-
2.51E+06, // CaCl2 = Ca+2 + 2Cl-
102-
2.69E+07, // MgSO4 = Mg+2 + SO4-2
103-
6.62E+07, // NaSO4- = Na+ + SO4-2
104-
8.55E-03 // CaCO3 + H+ = Ca+2 + HCO3-
105-
// 1 // Ca(OH)2​(s) + 2H+ = Ca2+ + 2H2​O (kinetic)
84+
constexpr CArrayWrapper<double, 10> reverseRates =
85+
{ 1.43E-03, // OH- + H+ = H2O
86+
8.92E+04, // CO2 + H2O = H+ + HCO3-
87+
4.67E-01, // CO3-2 + H+ = HCO3-
88+
1.85E+07, // CaHCO3+ = Ca+2 + HCO3-
89+
1.45E+07, // CaSO4 = Ca+2 + SO4-2
90+
2.14E+07, // CaCl+ = Ca+2 + Cl-
91+
2.51E+06, // CaCl2 = Ca+2 + 2Cl-
92+
2.69E+07, // MgSO4 = Mg+2 + SO4-2
93+
6.62E+07, // NaSO4- = Na+ + SO4-2
94+
3.00E-08 // CaCO3 + H+ = Ca+2 + HCO3-
10695
};
10796

108-
constexpr CArrayWrapper<int, 11> mobileSpeciesFlag =
97+
constexpr CArrayWrapper<int, 10> mobileSpeciesFlag =
10998
{ 1, // OH- + H+ = H2O
11099
1, // CO2 + H2O = H+ + HCO3-
111100
1, // CO3-2 + H+ = HCO3-
112-
1, // H2CO3 = H+ + HCO3-
113101
1, // CaHCO3+ = Ca+2 + HCO3-
114102
1, // CaSO4 = Ca+2 + SO4-2
115103
1, // CaCl+ = Ca+2 + Cl-
@@ -118,14 +106,16 @@ constexpr CArrayWrapper<int, 11> mobileSpeciesFlag =
118106
1, // NaSO4- = Na+ + SO4-2
119107
1 // CaCO3 + H+ = Ca+2 + HCO3-
120108
};
121-
}
122-
using carbonateSystemAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 18, 11, 0 >;
123-
using carbonateSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 18, 11, 11 >;
124-
using carbonateSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 17, 11, 10 >;
125109

126-
constexpr carbonateSystemAllKineticType carbonateSystemAllKinetic( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag );
127-
constexpr carbonateSystemAllEquilibriumType carbonateSystemAllEquilibrium( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag );
128-
constexpr carbonateSystemType carbonateSystem( carbonate::stoichMatrixNosolid, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag );
110+
}
111+
112+
using carbonateSystemAllKineticType = reactionsSystems::MixedReactionsParameters< double, int, int, 17, 10, 0 >;
113+
using carbonateSystemAllEquilibriumType = reactionsSystems::MixedReactionsParameters< double, int, int, 17, 10, 10 >;
114+
using carbonateSystemType = reactionsSystems::MixedReactionsParameters< double, int, int, 16, 10, 9 >;
115+
116+
constexpr carbonateSystemAllKineticType carbonateSystemAllKinetic( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag );
117+
constexpr carbonateSystemAllEquilibriumType carbonateSystemAllEquilibrium( carbonate::stoichMatrix, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag );
118+
constexpr carbonateSystemType carbonateSystem( carbonate::stoichMatrixNosolid, carbonate::equilibriumConstants, carbonate::forwardRates, carbonate::reverseRates, carbonate::mobileSpeciesFlag );
129119

130120
// *****UNCRUSTIFY-ON******
131121
} // namespace geochemistry

‎src/reactions/geochemistry/unitTests/testGeochemicalEquilibriumReactions.cpp‎

Lines changed: 27 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -38,12 +38,12 @@ TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium )
3838
{
3939
using namespace hpcReact::geochemistry;
4040

41-
double const initialSpeciesConcentration[18] =
41+
double const initialSpeciesConcentration[17] =
4242
{
4343
1.0e-16, // OH-
4444
1.0e-16, // CO2
4545
1.0e-16, // CO3-2
46-
1.0e-16, // H2CO3
46+
//1.0e-16, // H2CO3
4747
1.0e-16, // CaHCO3+
4848
1.0e-16, // CaSO4
4949
1.0e-16, // CaCl+
@@ -60,25 +60,24 @@ TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium )
6060
1.09 // Na+1
6161
};
6262

63-
double const expectedSpeciesConcentrations[18] =
64-
{ 2.327841695586879e-11, // OH-
65-
3.745973700632716e-01, // CO2
66-
3.956656978189456e-11, // CO3-2
67-
9.629355924567627e-04, // H2CO3
68-
6.739226982791492e-05, // CaHCO3+
69-
5.298329882666738e-03, // CaSO4
70-
5.844517547638333e-03, // CaCl+
71-
1.277319392670652e-02, // CaCl2
72-
6.618125707964991e-03, // MgSO4
73-
1.769217213462983e-02, // NaSO4-
74-
1.065032288527957e-09, // CaCO3
75-
4.396954721488358e-04, // H+
76-
3.723009698453808e-04, // HCO3-
77-
1.471656530812871e-02, // Ca+2
78-
2.491372274738741e-03, // SO4-2
79-
1.858609094598949e+00, // Cl-
80-
9.881874292035110e-03, // Mg+2
81-
1.072307827865370e+00 // Na+1
63+
double const expectedSpeciesConcentrations[17] =
64+
{ 2.1579694253441686e-11, // OH-
65+
0.3755789961165058, // CO2
66+
3.4214835005538611e-11, // CO3-2
67+
1.9160014879413049e-05, // CaHCO3+
68+
0.0025013721967110923, // CaSO4
69+
0.030083903919781853, // CaCl+
70+
0.0028032598559295028, // CaCl2
71+
0.0063383337566393907, // MgSO4
72+
0.019567697287351467, // NaSO4-
73+
4.754873524374959e-05, // CaCO3
74+
0.00046855267453226469, // H+
75+
0.00035429509915661095, // HCO3-
76+
0.0032447552774548935, // Ca+2
77+
0.0036925967592983458, // SO4-2
78+
1.8543095763683592, // Cl-
79+
0.01016166624336071, // Mg+2
80+
1.0704323027126488 // Na+1
8281
};
8382

8483
std::cout<<" RESIDUAL_FORM 0:"<<std::endl;
@@ -140,13 +139,13 @@ TEST( testEquilibriumReactions, testcarbonateSystemAllEquilibrium2 )
140139

141140
double const expectedPrimarySpeciesConcentrations[numPrimarySpecies] =
142141
{
143-
4.396954721488358e-04, // H+
144-
3.723009698453808e-04, // HCO3-
145-
1.471656530812871e-02, // Ca+2
146-
2.491372274738741e-03, // SO4-2
147-
1.858609094598949e+00, // Cl-
148-
9.881874292035110e-03, // Mg+2
149-
1.072307827865370e+00 // Na+1
142+
0.00046855267453254149, // H+
143+
0.00035429509915645743, // HCO3-
144+
0.0032447552774548518, // Ca+2
145+
0.0036925967592983211, // SO4-2
146+
1.8543095763683592, // Cl-
147+
0.010161666243360675, // Mg+2
148+
1.0704323027126488 // Na+1
150149
};
151150

152151
for( int r=0; r<numPrimarySpecies; ++r )

0 commit comments

Comments
 (0)