46 typedef ::Opm::IdealGas<Scalar>
IdealGas;
47 static const int liquidPhaseIdx = 0;
48 static const int gasPhaseIdx = 1;
58 template <
class Evaluation>
59 static Evaluation
gasDiffCoeff(
const Evaluation& temperature,
const Evaluation& pressure,
bool extrapolate =
false)
62 Scalar k = 1.3806504e-23;
64 Scalar R_h = 1.72e-10;
66 return k / (c * M_PI * R_h) * (temperature / mu);
75 template <
class Evaluation>
99 template <
class Evaluation>
101 const Evaluation& pg,
102 const Evaluation& salinity,
103 const int knownPhaseIdx,
106 const int& activityModel,
107 bool extrapolate =
false)
109 OPM_TIMEFUNCTION_LOCAL();
112 bool iterate =
false;
113 if ((activityModel == 1 && salinity > 0.0) || (activityModel == 2 && temperature > 372.15)) {
119 if (knownPhaseIdx < 0) {
120 Evaluation molalityNaCl = massFracToMolality_(salinity);
126 if (activityModel == 3) {
127 auto [xCO2, yH2O] = mutualSolubilitySpycherPruess2005_(temperature, pg, molalityNaCl, extrapolate);
135 auto [xCO2, yH2O] = fixPointIterSolubility_(temperature, pg, molalityNaCl, activityModel, extrapolate);
142 auto [xCO2, yH2O] = nonIterSolubility_(temperature, pg, molalityNaCl, activityModel, extrapolate);
152 else if (knownPhaseIdx == liquidPhaseIdx && activityModel == 3) {
153 Evaluation x_NaCl = salinityToMolFrac_(salinity);
154 const Evaluation& A = computeA_(temperature, pg, Evaluation(0.0), Evaluation(0.0),
false, extrapolate,
true);
155 ygH2O = A * (1 - xlCO2 - x_NaCl);
161 else if (knownPhaseIdx == gasPhaseIdx && activityModel == 3) {
163 Evaluation x_NaCl = salinityToMolFrac_(salinity);
164 const Evaluation& A = computeA_(temperature, pg, Evaluation(0.0), Evaluation(0.0),
false, extrapolate,
true);
165 xlCO2 = 1 - x_NaCl - ygH2O / A;
172 template <
class Evaluation>
173 static Evaluation
henry(
const Evaluation& temperature,
bool extrapolate =
false)
184 template <
class Evaluation>
186 const Evaluation& pg,
187 const Evaluation& yH2O,
189 bool extrapolate =
false,
190 bool spycherPruess2005 =
false)
192 OPM_TIMEFUNCTION_LOCAL();
193 Valgrind::CheckDefined(temperature);
194 Valgrind::CheckDefined(pg);
197 Evaluation pg_bar = pg / 1.e5;
201 Evaluation a_CO2 = aCO2_(temperature, highTemp);
202 Evaluation a_CO2_H2O = aCO2_H2O_(temperature, yH2O, highTemp);
203 Evaluation a_mix = aMix_(temperature, yH2O, highTemp);
204 Scalar b_CO2 = bCO2_(highTemp);
205 Evaluation b_mix = bMix_(yH2O, highTemp);
208 if (spycherPruess2005) {
209 lnPhiCO2 = log(V / (V - b_CO2));
210 lnPhiCO2 += b_CO2 / (V - b_CO2);
211 lnPhiCO2 -= 2 * a_CO2 / (R * pow(temperature, 1.5) * b_CO2) * log((V + b_CO2) / V);
215 * pow(temperature, 1.5)
218 * (log((V + b_CO2) / V)
219 - b_CO2 / (V + b_CO2));
220 lnPhiCO2 -= log(pg_bar * V / (R * temperature));
223 lnPhiCO2 = (b_CO2 / b_mix) * (pg_bar * V / (R * temperature) - 1);
224 lnPhiCO2 -= log(pg_bar * (V - b_mix) / (R * temperature));
225 lnPhiCO2 += (2 * (yH2O * a_CO2_H2O + (1 - yH2O) * a_CO2) / a_mix - (b_CO2 / b_mix)) *
226 a_mix / (b_mix * R * pow(temperature, 1.5)) * log(V / (V + b_mix));
228 return exp(lnPhiCO2);
239 template <
class Evaluation>
241 const Evaluation& pg,
242 const Evaluation& yH2O,
244 bool extrapolate =
false,
245 bool spycherPruess2005 =
false)
247 OPM_TIMEFUNCTION_LOCAL();
248 Valgrind::CheckDefined(temperature);
249 Valgrind::CheckDefined(pg);
252 const Evaluation& pg_bar = pg / 1.e5;
256 Evaluation a_H2O = aH2O_(temperature, highTemp);
257 Evaluation a_CO2_H2O = aCO2_H2O_(temperature, yH2O, highTemp);
258 Evaluation a_mix = aMix_(temperature, yH2O, highTemp);
259 Scalar b_H2O = bH2O_(highTemp);
260 Evaluation b_mix = bMix_(yH2O, highTemp);
263 if (spycherPruess2005) {
266 + b_H2O/(V - b_mix) - 2*a_CO2_H2O
267 / (R*pow(temperature, 1.5)*b_mix)*log((V + b_mix)/V)
268 + a_mix*b_H2O/(R*pow(temperature, 1.5)*b_mix*b_mix)
269 *(log((V + b_mix)/V) - b_mix/(V + b_mix))
270 - log(pg_bar*V/(R*temperature));
273 lnPhiH2O = (b_H2O / b_mix) * (pg_bar * V / (R * temperature) - 1);
274 lnPhiH2O -= log(pg_bar * (V - b_mix) / (R * temperature));
275 lnPhiH2O += (2 * (yH2O * a_H2O + (1 - yH2O) * a_CO2_H2O) / a_mix - (b_H2O / b_mix)) *
276 a_mix / (b_mix * R * pow(temperature, 1.5)) * log(V / (V + b_mix));
278 return exp(lnPhiH2O);
285 template <
class Evaluation>
286 static Evaluation aCO2_(
const Evaluation& temperature,
const bool& highTemp)
289 return 8.008e7 - 4.984e4 * temperature;
292 return 7.54e7 - 4.13e4 * temperature;
299 template <
class Evaluation>
300 static Evaluation aH2O_(
const Evaluation& temperature,
const bool& highTemp)
303 return 1.337e8 - 1.4e4 * temperature;
313 template <
class Evaluation>
314 static Evaluation aCO2_H2O_(
const Evaluation& temperature,
const Evaluation& yH2O,
const bool& highTemp)
318 Evaluation aCO2 = aCO2_(temperature, highTemp);
319 Evaluation aH2O = aH2O_(temperature, highTemp);
322 Evaluation K_CO2_H2O = 0.4228 - 7.422e-4 * temperature;
323 Evaluation K_H2O_CO2 = 1.427e-2 - 4.037e-4 * temperature;
324 Evaluation k_CO2_H2O = yH2O * K_H2O_CO2 + (1 - yH2O) * K_CO2_H2O;
327 return sqrt(aCO2 * aH2O) * (1 - k_CO2_H2O);
337 template <
class Evaluation>
338 static Evaluation aMix_(
const Evaluation& temperature,
const Evaluation& yH2O,
const bool& highTemp)
342 Evaluation aCO2 = aCO2_(temperature, highTemp);
343 Evaluation aH2O = aH2O_(temperature, highTemp);
344 Evaluation a_CO2_H2O = aCO2_H2O_(temperature, yH2O, highTemp);
346 return yH2O * yH2O * aH2O + 2 * yH2O * (1 - yH2O) * a_CO2_H2O + (1 - yH2O) * (1 - yH2O) * aCO2;
349 return aCO2_(temperature, highTemp);
356 static Scalar bCO2_(
const bool& highTemp)
369 static Scalar bH2O_(
const bool& highTemp)
382 template <
class Evaluation>
383 static Evaluation bMix_(
const Evaluation& yH2O,
const bool& highTemp)
387 Scalar bCO2 = bCO2_(highTemp);
388 Scalar bH2O = bH2O_(highTemp);
390 return yH2O * bH2O + (1 - yH2O) * bCO2;
393 return bCO2_(highTemp);
400 template <
class Evaluation>
401 static Evaluation V_avg_CO2_(
const Evaluation& temperature,
const bool& highTemp)
403 if (highTemp && (temperature > 373.15)) {
404 return 32.6 + 3.413e-2 * (temperature - 373.15);
414 template <
class Evaluation>
415 static Evaluation V_avg_H2O_(
const Evaluation& temperature,
const bool& highTemp)
417 if (highTemp && (temperature > 373.15)) {
418 return 18.1 + 3.137e-2 * (temperature - 373.15);
428 template <
class Evaluation>
429 static Evaluation AM_(
const Evaluation& temperature,
const bool& highTemp)
431 if (highTemp && temperature > 373.15) {
432 Evaluation deltaTk = temperature - 373.15;
433 return deltaTk * (-3.084e-2 + 1.927e-5 * deltaTk);
443 template <
class Evaluation>
444 static Evaluation Pref_(
const Evaluation& temperature,
const bool& highTemp)
446 if (highTemp && temperature > 373.15) {
447 const Evaluation& temperatureCelcius = temperature - 273.15;
448 static const Scalar c[5] = { -1.9906e-1, 2.0471e-3, 1.0152e-4, -1.4234e-6, 1.4168e-8 };
449 return c[0] + temperatureCelcius * (c[1] + temperatureCelcius * (c[2] +
450 temperatureCelcius * (c[3] + temperatureCelcius * c[4])));
460 template <
class Evaluation>
461 static Evaluation activityCoefficientCO2_(
const Evaluation& temperature,
462 const Evaluation& xCO2,
463 const bool& highTemp)
467 Evaluation AM = AM_(temperature, highTemp);
468 Evaluation lnGammaCO2 = 2 * AM * xCO2 * (1 - xCO2) * (1 - xCO2);
469 return exp(lnGammaCO2);
479 template <
class Evaluation>
480 static Evaluation activityCoefficientH2O_(
const Evaluation& temperature,
481 const Evaluation& xCO2,
482 const bool& highTemp)
486 Evaluation AM = AM_(temperature, highTemp);
487 Evaluation lnGammaH2O = (1 - 2 * (1 - xCO2)) * AM * xCO2 * xCO2;
488 return exp(lnGammaH2O);
500 template <
class Evaluation>
501 static Evaluation salinityToMolFrac_(
const Evaluation& salinity) {
502 OPM_TIMEFUNCTION_LOCAL();
504 const Scalar Ms = 58.44e-3;
506 const Evaluation X_NaCl = salinity;
508 const Evaluation x_NaCl = -Mw * X_NaCl / ((Ms - Mw) * X_NaCl - Ms);
517 template <
class Evaluation>
518 static Evaluation moleFracToMolality_(
const Evaluation& x_NaCl)
521 return 55.508 * x_NaCl / (1 - x_NaCl);
524 template <
class Evaluation>
525 static Evaluation massFracToMolality_(
const Evaluation& X_NaCl)
527 const Scalar MmNaCl = 58.44e-3;
528 return X_NaCl / (MmNaCl * (1 - X_NaCl));
536 template <
class Evaluation>
537 static Evaluation molalityToMoleFrac_(
const Evaluation& m_NaCl)
540 return m_NaCl / (55.508 + m_NaCl);
546 template <
class Evaluation>
547 static std::pair<Evaluation, Evaluation> fixPointIterSolubility_(
const Evaluation& temperature,
548 const Evaluation& pg,
549 const Evaluation& m_NaCl,
550 const int& activityModel,
551 bool extrapolate =
false)
553 OPM_TIMEFUNCTION_LOCAL();
556 Evaluation xCO2 = 0.009;
557 Evaluation gammaNaCl = 1.0;
560 if (m_NaCl > 0.0 && activityModel == 2) {
561 gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, Evaluation(0.0), activityModel);
567 bool highTemp =
true;
568 if (activityModel == 1) {
571 const bool iterate =
true;
574 for (
int i = 0; i < max_iter; ++i) {
576 if (m_NaCl > 0.0 && activityModel == 1) {
577 gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, xCO2, activityModel);
581 auto [xCO2_new, yH2O_new] = mutualSolubility_(temperature, pg, xCO2, yH2O, m_NaCl, gammaNaCl, highTemp,
582 iterate, extrapolate);
585 if (abs(xCO2_new - xCO2) < tol && abs(yH2O_new - yH2O) < tol) {
604 template <
class Evaluation>
605 static std::pair<Evaluation, Evaluation> nonIterSolubility_(
const Evaluation& temperature,
606 const Evaluation& pg,
607 const Evaluation& m_NaCl,
608 const int& activityModel,
609 bool extrapolate =
false)
612 Evaluation gammaNaCl = 1.0;
613 if (m_NaCl > 0.0 && activityModel > 0 && activityModel < 3) {
614 gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, Evaluation(0.0), activityModel);
619 const bool highTemp =
false;
620 const bool iterate =
false;
621 auto [xCO2, yH2O] = mutualSolubility_(temperature, pg, Evaluation(0.0), Evaluation(0.0), m_NaCl, gammaNaCl,
622 highTemp, iterate, extrapolate);
630 template <
class Evaluation>
631 static std::pair<Evaluation, Evaluation> mutualSolubility_(
const Evaluation& temperature,
632 const Evaluation& pg,
633 const Evaluation& xCO2,
634 const Evaluation& yH2O,
635 const Evaluation& m_NaCl,
636 const Evaluation& gammaNaCl,
637 const bool& highTemp,
639 bool extrapolate =
false)
642 const Evaluation& A = computeA_(temperature, pg, yH2O, xCO2, highTemp, extrapolate);
643 Evaluation B = computeB_(temperature, pg, yH2O, xCO2, highTemp, extrapolate);
649 Evaluation yH2O_new = (1. - B) * 55.508 / ((1. / A - B) * (2 * m_NaCl + 55.508) + 2 * m_NaCl * B);
652 xCO2_new = B * (1 - yH2O);
655 xCO2_new = B * (1 - yH2O_new);
658 return {xCO2_new, yH2O_new};
664 template <
class Evaluation>
665 static std::pair<Evaluation, Evaluation> mutualSolubilitySpycherPruess2005_(
const Evaluation& temperature,
666 const Evaluation& pg,
667 const Evaluation& m_NaCl,
668 bool extrapolate =
false)
671 const Evaluation& A = computeA_(temperature, pg, Evaluation(0.0), Evaluation(0.0),
false, extrapolate,
true);
672 const Evaluation& B = computeB_(temperature, pg, Evaluation(0.0), Evaluation(0.0),
false, extrapolate,
true);
675 Evaluation yH2O = (1 - B) / (1. / A - B);
676 Evaluation xCO2 = B * (1 - yH2O);
680 const Evaluation& gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, Evaluation(0.0), 3);
683 Evaluation mCO2 = (xCO2 * 55.508) / (1 - xCO2);
685 xCO2 = mCO2 / (m_NaCl + 55.508 + mCO2);
688 const Evaluation& xNaCl = molalityToMoleFrac_(m_NaCl);
689 yH2O = A * (1 - xCO2 - xNaCl);
703 template <
class Evaluation>
704 static Evaluation computeA_(
const Evaluation& temperature,
705 const Evaluation& pg,
706 const Evaluation& yH2O,
707 const Evaluation& xCO2,
708 const bool& highTemp,
709 bool extrapolate =
false,
710 bool spycherPruess2005 =
false)
712 OPM_TIMEFUNCTION_LOCAL();
714 const Evaluation& deltaP = pg / 1e5 - Pref_(temperature, highTemp);
715 Evaluation v_av_H2O = V_avg_H2O_(temperature, highTemp);
716 Evaluation k0_H2O = equilibriumConstantH2O_(temperature, highTemp);
717 Evaluation phi_H2O =
fugacityCoefficientH2O(temperature, pg, yH2O, highTemp, extrapolate, spycherPruess2005);
718 Evaluation gammaH2O = activityCoefficientH2O_(temperature, xCO2, highTemp);
722 if ( temperature > 372.15 && temperature < 382.15 && !spycherPruess2005) {
723 const Evaluation weight = (382.15 - temperature) / 10.;
724 const Evaluation& k0_H2O_low = equilibriumConstantH2O_(temperature,
false);
725 const Evaluation& phi_H2O_low =
fugacityCoefficientH2O(temperature, pg, Evaluation(0.0),
false, extrapolate);
726 k0_H2O = k0_H2O * (1 - weight) + k0_H2O_low * weight;
727 phi_H2O = phi_H2O * (1 - weight) + phi_H2O_low * weight;
731 const Evaluation& pg_bar = pg / 1.e5;
733 return k0_H2O * gammaH2O / (phi_H2O * pg_bar) * exp(deltaP * v_av_H2O / (R * temperature));
744 template <
class Evaluation>
745 static Evaluation computeB_(
const Evaluation& temperature,
746 const Evaluation& pg,
747 const Evaluation& yH2O,
748 const Evaluation& xCO2,
749 const bool& highTemp,
750 bool extrapolate =
false,
751 bool spycherPruess2005 =
false)
753 OPM_TIMEFUNCTION_LOCAL();
755 const Evaluation& deltaP = pg / 1e5 - Pref_(temperature, highTemp);
756 Evaluation v_av_CO2 = V_avg_CO2_(temperature, highTemp);
757 Evaluation k0_CO2 = equilibriumConstantCO2_(temperature, pg, highTemp, spycherPruess2005);
758 Evaluation phi_CO2 =
fugacityCoefficientCO2(temperature, pg, yH2O, highTemp, extrapolate, spycherPruess2005);
759 Evaluation gammaCO2 = activityCoefficientCO2_(temperature, xCO2, highTemp);
763 if ( temperature > 372.15 && temperature < 382.15 && !spycherPruess2005) {
764 const Evaluation weight = (382.15 - temperature) / 10.;
765 const Evaluation& k0_CO2_low = equilibriumConstantCO2_(temperature, pg,
false, spycherPruess2005);
766 const Evaluation& phi_CO2_low =
fugacityCoefficientCO2(temperature, pg, Evaluation(0.0),
false, extrapolate);
767 k0_CO2 = k0_CO2 * (1 - weight) + k0_CO2_low * weight;
768 phi_CO2 = phi_CO2 * (1 - weight) + phi_CO2_low * weight;
772 const Evaluation& pg_bar = pg / 1.e5;
774 return phi_CO2 * pg_bar / (55.508 * k0_CO2 * gammaCO2) * exp(-deltaP * v_av_CO2 / (R * temperature));
780 template <
class Evaluation>
781 static Evaluation activityCoefficientSalt_(
const Evaluation& temperature,
782 const Evaluation& pg,
783 const Evaluation& m_NaCl,
784 const Evaluation& xCO2,
785 const int& activityModel)
787 OPM_TIMEFUNCTION_LOCAL();
794 if (activityModel == 1) {
795 lambda = computeLambdaRumpfetal_(temperature);
797 Evaluation m_CO2 = xCO2 * (2 * m_NaCl + 55.508) / (1 - xCO2);
798 convTerm = (1 + (m_CO2 + 2 * m_NaCl) / 55.508) / (1 + m_CO2 / 55.508);
800 else if (activityModel == 2) {
801 lambda = computeLambdaSpycherPruess2009_(temperature);
802 xi = computeXiSpycherPruess2009_(temperature);
803 convTerm = 1 + 2 * m_NaCl / 55.508;
805 else if (activityModel == 3) {
806 lambda = computeLambdaDuanSun_(temperature, pg);
807 xi = computeXiDuanSun_(temperature, pg);
811 throw std::runtime_error(
"Activity model for salt-out effect has not been implemented!");
815 const Evaluation& lnGamma = 2 * lambda * m_NaCl + xi * m_NaCl * m_NaCl;
818 return convTerm * exp(lnGamma);
824 template <
class Evaluation>
825 static Evaluation computeLambdaSpycherPruess2009_(
const Evaluation& temperature)
828 static const Scalar c[3] = { 2.217e-4, 1.074, 2648. };
831 return c[0] * temperature + c[1] / temperature + c[2] / (temperature * temperature);
837 template <
class Evaluation>
838 static Evaluation computeXiSpycherPruess2009_(
const Evaluation& temperature)
841 static const Scalar c[3] = { 1.3e-5, -20.12, 5259. };
844 return c[0] * temperature + c[1] / temperature + c[2] / (temperature * temperature);
850 template <
class Evaluation>
851 static Evaluation computeLambdaRumpfetal_(
const Evaluation& temperature)
854 static const Scalar c[4] = { 0.254, -76.82, -10656, 6312e3 };
856 return c[0] + c[1] / temperature + c[2] / (temperature * temperature) +
857 c[3] / (temperature * temperature * temperature);
867 template <
class Evaluation>
868 static Evaluation computeLambdaDuanSun_(
const Evaluation& temperature,
const Evaluation& pg)
870 static const Scalar c[6] =
871 { -0.411370585, 6.07632013E-4, 97.5347708, -0.0237622469, 0.0170656236, 1.41335834E-5 };
873 Evaluation pg_bar = pg / 1.0E5;
874 return c[0] + c[1]*temperature + c[2]/temperature + c[3]*pg_bar/temperature + c[4]*pg_bar/(630.0 - temperature)
875 + c[5]*temperature*log(pg_bar);
885 template <
class Evaluation>
886 static Evaluation computeXiDuanSun_(
const Evaluation& temperature,
const Evaluation& pg)
888 static const Scalar c[4] =
889 { 3.36389723E-4, -1.98298980E-5, 2.12220830E-3, -5.24873303E-3 };
891 Evaluation pg_bar = pg / 1.0E5;
892 return c[0] + c[1]*temperature + c[2]*pg_bar/temperature + c[3]*pg_bar/(630.0 - temperature);
901 template <
class Evaluation>
902 static Evaluation equilibriumConstantCO2_(
const Evaluation& temperature,
903 const Evaluation& pg,
904 const bool& highTemp,
905 bool spycherPruess2005 =
false)
907 OPM_TIMEFUNCTION_LOCAL();
908 Evaluation temperatureCelcius = temperature - 273.15;
909 std::array<Scalar, 4> c;
911 c = { 1.668, 3.992e-3, -1.156e-5, 1.593e-9 };
916 if (temperatureCelcius < 31 && pg > psat && !spycherPruess2005) {
917 c = { 1.169, 1.368e-2, -5.38e-5, 0.0 };
920 c = { 1.189, 1.304e-2, -5.446e-5, 0.0 };
923 Evaluation logk0_CO2 = c[0] + temperatureCelcius * (c[1] + temperatureCelcius *
924 (c[2] + temperatureCelcius * c[3]));
925 Evaluation k0_CO2 = pow(10.0, logk0_CO2);
935 template <
class Evaluation>
936 static Evaluation equilibriumConstantH2O_(
const Evaluation& temperature,
const bool& highTemp)
938 Evaluation temperatureCelcius = temperature - 273.15;
939 std::array<Scalar, 5> c;
941 c = { -2.1077, 2.8127e-2, -8.4298e-5, 1.4969e-7, -1.1812e-10 };
944 c = { -2.209, 3.097e-2, -1.098e-4, 2.048e-7, 0.0 };
946 Evaluation logk0_H2O = c[0] + temperatureCelcius * (c[1] + temperatureCelcius * (c[2] +
947 temperatureCelcius * (c[3] + temperatureCelcius * c[4])));
948 return pow(10.0, logk0_H2O);