47 typedef ::Opm::IdealGas<Scalar>
IdealGas;
48 static const int liquidPhaseIdx = 0;
49 static const int gasPhaseIdx = 1;
59 template <
class Evaluation,
class CO2Params>
60 OPM_HOST_DEVICE
static Evaluation
gasDiffCoeff(
const CO2Params& params,
const Evaluation& temperature,
const Evaluation& pressure,
bool extrapolate =
false)
63 Scalar k = 1.3806504e-23;
65 Scalar R_h = 1.72e-10;
66 const Evaluation& mu =
CO2::gasViscosity(params, temperature, pressure, extrapolate);
67 return k / (c * M_PI * R_h) * (temperature / mu);
76 template <
class Evaluation>
77 OPM_HOST_DEVICE
static Evaluation
liquidDiffCoeff(
const Evaluation& ,
const Evaluation& )
100 template <
class Evaluation,
class CO2Params>
102 const Evaluation& temperature,
103 const Evaluation& pg,
104 const Evaluation& salinity,
105 const int knownPhaseIdx,
108 const int& activityModel,
109 bool extrapolate =
false)
111 OPM_TIMEFUNCTION_LOCAL();
114 bool iterate =
false;
115 if ((activityModel == 1 && salinity > 0.0) || (activityModel == 2 && temperature > 372.15)) {
121 if (knownPhaseIdx < 0) {
122 Evaluation molalityNaCl = massFracToMolality_(salinity);
128 if (activityModel == 3) {
129 auto [xCO2, yH2O] = mutualSolubilitySpycherPruess2005_(params, temperature, pg, molalityNaCl, extrapolate);
137 auto [xCO2, yH2O] = fixPointIterSolubility_(params, temperature, pg, molalityNaCl, activityModel, extrapolate);
144 auto [xCO2, yH2O] = nonIterSolubility_(params, temperature, pg, molalityNaCl, activityModel, extrapolate);
154 else if (knownPhaseIdx == liquidPhaseIdx && activityModel == 3) {
155 Evaluation x_NaCl = salinityToMolFrac_(salinity);
156 const Evaluation& A = computeA_(params, temperature, pg, Evaluation(0.0), Evaluation(0.0),
false, extrapolate,
true);
157 ygH2O = A * (1 - xlCO2 - x_NaCl);
163 else if (knownPhaseIdx == gasPhaseIdx && activityModel == 3) {
165 Evaluation x_NaCl = salinityToMolFrac_(salinity);
166 const Evaluation& A = computeA_(params, temperature, pg, Evaluation(0.0), Evaluation(0.0),
false, extrapolate,
true);
167 xlCO2 = 1 - x_NaCl - ygH2O / A;
174 template <
class Evaluation>
175 static Evaluation
henry(
const Evaluation& temperature,
bool extrapolate =
false)
186 template <
class Evaluation,
class CO2Params>
188 const Evaluation& temperature,
189 const Evaluation& pg,
190 const Evaluation& yH2O,
192 bool extrapolate =
false,
193 bool spycherPruess2005 =
false)
195 OPM_TIMEFUNCTION_LOCAL();
196 Valgrind::CheckDefined(temperature);
197 Valgrind::CheckDefined(pg);
200 const Evaluation pg_bar = pg / 1.e5;
204 const Evaluation a_CO2 = aCO2_(temperature, highTemp);
205 const Evaluation a_CO2_H2O = aCO2_H2O_(temperature, yH2O, highTemp);
206 const Evaluation a_mix = aMix_(temperature, yH2O, highTemp);
207 const Scalar b_CO2 = bCO2_(highTemp);
208 const Evaluation b_mix = bMix_(yH2O, highTemp);
209 const Evaluation Rt15 = R * pow(temperature, 1.5);
212 if (spycherPruess2005) {
213 const Evaluation logVpb_V = log((V + b_CO2) / V);
214 lnPhiCO2 = log(V / (V - b_CO2));
215 lnPhiCO2 += b_CO2 / (V - b_CO2);
216 lnPhiCO2 -= 2 * a_CO2 / (Rt15 * b_CO2) * logVpb_V;
223 - b_CO2 / (V + b_CO2));
224 lnPhiCO2 -= log(pg_bar * V / (R * temperature));
227 lnPhiCO2 = (b_CO2 / b_mix) * (pg_bar * V / (R * temperature) - 1);
228 lnPhiCO2 -= log(pg_bar * (V - b_mix) / (R * temperature));
229 lnPhiCO2 += (2 * (yH2O * a_CO2_H2O + (1 - yH2O) * a_CO2) / a_mix - (b_CO2 / b_mix)) *
230 a_mix / (b_mix * Rt15) * log(V / (V + b_mix));
232 return exp(lnPhiCO2);
243 template <
class Evaluation,
class CO2Params>
245 const Evaluation& temperature,
246 const Evaluation& pg,
247 const Evaluation& yH2O,
249 bool extrapolate =
false,
250 bool spycherPruess2005 =
false)
252 OPM_TIMEFUNCTION_LOCAL();
253 Valgrind::CheckDefined(temperature);
254 Valgrind::CheckDefined(pg);
257 const Evaluation& pg_bar = pg / 1.e5;
261 const Evaluation a_H2O = aH2O_(temperature, highTemp);
262 const Evaluation a_CO2_H2O = aCO2_H2O_(temperature, yH2O, highTemp);
263 const Evaluation a_mix = aMix_(temperature, yH2O, highTemp);
264 const Scalar b_H2O = bH2O_(highTemp);
265 const Evaluation b_mix = bMix_(yH2O, highTemp);
266 const Evaluation Rt15 = R * pow(temperature, 1.5);
269 if (spycherPruess2005) {
270 const Evaluation logVpb_V = log((V + b_mix) / V);
273 + b_H2O/(V - b_mix) - 2*a_CO2_H2O
274 / (Rt15*b_mix)*logVpb_V
275 + a_mix*b_H2O/(Rt15*b_mix*b_mix)
276 *(logVpb_V - b_mix/(V + b_mix))
277 - log(pg_bar*V/(R*temperature));
280 lnPhiH2O = (b_H2O / b_mix) * (pg_bar * V / (R * temperature) - 1);
281 lnPhiH2O -= log(pg_bar * (V - b_mix) / (R * temperature));
282 lnPhiH2O += (2 * (yH2O * a_H2O + (1 - yH2O) * a_CO2_H2O) / a_mix - (b_H2O / b_mix)) *
283 a_mix / (b_mix * Rt15) * log(V / (V + b_mix));
285 return exp(lnPhiH2O);
292 template <
class Evaluation>
293 OPM_HOST_DEVICE
static Evaluation aCO2_(
const Evaluation& temperature,
const bool& highTemp)
296 return 8.008e7 - 4.984e4 * temperature;
299 return 7.54e7 - 4.13e4 * temperature;
306 template <
class Evaluation>
307 OPM_HOST_DEVICE
static Evaluation aH2O_(
const Evaluation& temperature,
const bool& highTemp)
310 return 1.337e8 - 1.4e4 * temperature;
320 template <
class Evaluation>
321 OPM_HOST_DEVICE
static Evaluation aCO2_H2O_(
const Evaluation& temperature,
const Evaluation& yH2O,
const bool& highTemp)
325 Evaluation aCO2 = aCO2_(temperature, highTemp);
326 Evaluation aH2O = aH2O_(temperature, highTemp);
329 Evaluation K_CO2_H2O = 0.4228 - 7.422e-4 * temperature;
330 Evaluation K_H2O_CO2 = 1.427e-2 - 4.037e-4 * temperature;
331 Evaluation k_CO2_H2O = yH2O * K_H2O_CO2 + (1 - yH2O) * K_CO2_H2O;
334 return sqrt(aCO2 * aH2O) * (1 - k_CO2_H2O);
344 template <
class Evaluation>
345 OPM_HOST_DEVICE
static Evaluation aMix_(
const Evaluation& temperature,
const Evaluation& yH2O,
const bool& highTemp)
349 Evaluation aCO2 = aCO2_(temperature, highTemp);
350 Evaluation aH2O = aH2O_(temperature, highTemp);
351 Evaluation a_CO2_H2O = aCO2_H2O_(temperature, yH2O, highTemp);
353 return yH2O * yH2O * aH2O + 2 * yH2O * (1 - yH2O) * a_CO2_H2O + (1 - yH2O) * (1 - yH2O) * aCO2;
356 return aCO2_(temperature, highTemp);
363 OPM_HOST_DEVICE
static Scalar bCO2_(
const bool& highTemp)
376 OPM_HOST_DEVICE
static Scalar bH2O_(
const bool& highTemp)
389 template <
class Evaluation>
390 OPM_HOST_DEVICE
static Evaluation bMix_(
const Evaluation& yH2O,
const bool& highTemp)
394 Scalar bCO2 = bCO2_(highTemp);
395 Scalar bH2O = bH2O_(highTemp);
397 return yH2O * bH2O + (1 - yH2O) * bCO2;
400 return bCO2_(highTemp);
407 template <
class Evaluation>
408 OPM_HOST_DEVICE
static Evaluation V_avg_CO2_(
const Evaluation& temperature,
const bool& highTemp)
410 if (highTemp && (temperature > 373.15)) {
411 return 32.6 + 3.413e-2 * (temperature - 373.15);
421 template <
class Evaluation>
422 OPM_HOST_DEVICE
static Evaluation V_avg_H2O_(
const Evaluation& temperature,
const bool& highTemp)
424 if (highTemp && (temperature > 373.15)) {
425 return 18.1 + 3.137e-2 * (temperature - 373.15);
435 template <
class Evaluation>
436 OPM_HOST_DEVICE
static Evaluation AM_(
const Evaluation& temperature,
const bool& highTemp)
438 if (highTemp && temperature > 373.15) {
439 Evaluation deltaTk = temperature - 373.15;
440 return deltaTk * (-3.084e-2 + 1.927e-5 * deltaTk);
450 template <
class Evaluation>
451 OPM_HOST_DEVICE
static Evaluation Pref_(
const Evaluation& temperature,
const bool& highTemp)
453 if (highTemp && temperature > 373.15) {
454 const Evaluation& temperatureCelcius = temperature - 273.15;
455 static const Scalar c[5] = { -1.9906e-1, 2.0471e-3, 1.0152e-4, -1.4234e-6, 1.4168e-8 };
456 return c[0] + temperatureCelcius * (c[1] + temperatureCelcius * (c[2] +
457 temperatureCelcius * (c[3] + temperatureCelcius * c[4])));
467 template <
class Evaluation>
468 OPM_HOST_DEVICE
static Evaluation activityCoefficientCO2_(
const Evaluation& temperature,
469 const Evaluation& xCO2,
470 const bool& highTemp)
474 Evaluation AM = AM_(temperature, highTemp);
475 Evaluation lnGammaCO2 = 2 * AM * xCO2 * (1 - xCO2) * (1 - xCO2);
476 return exp(lnGammaCO2);
486 template <
class Evaluation>
487 OPM_HOST_DEVICE
static Evaluation activityCoefficientH2O_(
const Evaluation& temperature,
488 const Evaluation& xCO2,
489 const bool& highTemp)
493 Evaluation AM = AM_(temperature, highTemp);
494 Evaluation lnGammaH2O = (1 - 2 * (1 - xCO2)) * AM * xCO2 * xCO2;
495 return exp(lnGammaH2O);
507 template <
class Evaluation>
508 OPM_HOST_DEVICE
static Evaluation salinityToMolFrac_(
const Evaluation& salinity) {
509 OPM_TIMEFUNCTION_LOCAL();
511 const Scalar Ms = 58.44e-3;
513 const Evaluation X_NaCl = salinity;
515 const Evaluation x_NaCl = -Mw * X_NaCl / ((Ms - Mw) * X_NaCl - Ms);
525 template <
class Evaluation>
526 OPM_HOST_DEVICE
static Evaluation moleFracToMolality_(
const Evaluation& x_NaCl)
529 return 55.508 * x_NaCl / (1 - x_NaCl);
533 template <
class Evaluation>
534 OPM_HOST_DEVICE
static Evaluation massFracToMolality_(
const Evaluation& X_NaCl)
536 const Scalar MmNaCl = 58.44e-3;
537 return X_NaCl / (MmNaCl * (1 - X_NaCl));
545 template <
class Evaluation>
546 OPM_HOST_DEVICE
static Evaluation molalityToMoleFrac_(
const Evaluation& m_NaCl)
549 return m_NaCl / (55.508 + m_NaCl);
555 template <
class Evaluation,
class CO2Parameters>
556 OPM_HOST_DEVICE
static std::pair<Evaluation, Evaluation> fixPointIterSolubility_(
const CO2Parameters& params,
557 const Evaluation& temperature,
558 const Evaluation& pg,
559 const Evaluation& m_NaCl,
560 const int& activityModel,
561 bool extrapolate =
false)
563 OPM_TIMEFUNCTION_LOCAL();
566 Evaluation xCO2 = 0.009;
567 Evaluation gammaNaCl = 1.0;
570 if (m_NaCl > 0.0 && activityModel == 2) {
571 gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, Evaluation(0.0), activityModel);
577 bool highTemp =
true;
578 if (activityModel == 1) {
581 const bool iterate =
true;
584 for (
int i = 0; i < max_iter; ++i) {
586 if (m_NaCl > 0.0 && activityModel == 1) {
587 gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, xCO2, activityModel);
591 auto [xCO2_new, yH2O_new] = mutualSolubility_(params, temperature, pg, xCO2, yH2O, m_NaCl, gammaNaCl, highTemp,
592 iterate, extrapolate);
595 if (abs(xCO2_new - xCO2) < tol && abs(yH2O_new - yH2O) < tol) {
614 template <
class Evaluation,
class CO2Parameters>
615 OPM_HOST_DEVICE
static std::pair<Evaluation, Evaluation> nonIterSolubility_(
const CO2Parameters& params,
616 const Evaluation& temperature,
617 const Evaluation& pg,
618 const Evaluation& m_NaCl,
619 const int& activityModel,
620 bool extrapolate =
false)
623 Evaluation gammaNaCl = 1.0;
624 if (m_NaCl > 0.0 && activityModel > 0 && activityModel < 3) {
625 gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, Evaluation(0.0), activityModel);
630 const bool highTemp =
false;
631 const bool iterate =
false;
632 auto [xCO2, yH2O] = mutualSolubility_(params, temperature, pg, Evaluation(0.0), Evaluation(0.0), m_NaCl, gammaNaCl,
633 highTemp, iterate, extrapolate);
641 template <
class Evaluation,
class CO2Parameters>
642 OPM_HOST_DEVICE
static std::pair<Evaluation, Evaluation> mutualSolubility_(
const CO2Parameters& params,
643 const Evaluation& temperature,
644 const Evaluation& pg,
645 const Evaluation& xCO2,
646 const Evaluation& yH2O,
647 const Evaluation& m_NaCl,
648 const Evaluation& gammaNaCl,
649 const bool& highTemp,
651 bool extrapolate =
false)
654 const Evaluation& A = computeA_(params, temperature, pg, yH2O, xCO2, highTemp, extrapolate);
655 Evaluation B = computeB_(params, temperature, pg, yH2O, xCO2, highTemp, extrapolate);
661 Evaluation yH2O_new = (1. - B) * 55.508 / ((1. / A - B) * (2 * m_NaCl + 55.508) + 2 * m_NaCl * B);
664 xCO2_new = B * (1 - yH2O);
667 xCO2_new = B * (1 - yH2O_new);
670 return {xCO2_new, yH2O_new};
676 template <
class Evaluation,
class CO2Parameters>
677 OPM_HOST_DEVICE
static std::pair<Evaluation, Evaluation> mutualSolubilitySpycherPruess2005_(
const CO2Parameters& params,
678 const Evaluation& temperature,
679 const Evaluation& pg,
680 const Evaluation& m_NaCl,
681 bool extrapolate =
false)
684 const Evaluation& A = computeA_(params, temperature, pg, Evaluation(0.0), Evaluation(0.0),
false, extrapolate,
true);
685 const Evaluation& B = computeB_(params, temperature, pg, Evaluation(0.0), Evaluation(0.0),
false, extrapolate,
true);
688 Evaluation yH2O = (1 - B) / (1. / A - B);
689 Evaluation xCO2 = B * (1 - yH2O);
693 const Evaluation& gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, Evaluation(0.0), 3);
696 Evaluation mCO2 = (xCO2 * 55.508) / (1 - xCO2);
698 xCO2 = mCO2 / (m_NaCl + 55.508 + mCO2);
701 const Evaluation& xNaCl = molalityToMoleFrac_(m_NaCl);
702 yH2O = A * (1 - xCO2 - xNaCl);
716 template <
class Evaluation,
class CO2Params>
717 OPM_HOST_DEVICE
static Evaluation computeA_(
const CO2Params& params,
718 const Evaluation& temperature,
719 const Evaluation& pg,
720 const Evaluation& yH2O,
721 const Evaluation& xCO2,
722 const bool& highTemp,
723 bool extrapolate =
false,
724 bool spycherPruess2005 =
false)
726 OPM_TIMEFUNCTION_LOCAL();
728 const Evaluation& deltaP = pg / 1e5 - Pref_(temperature, highTemp);
729 Evaluation v_av_H2O = V_avg_H2O_(temperature, highTemp);
730 Evaluation k0_H2O = equilibriumConstantH2O_(temperature, highTemp);
731 Evaluation phi_H2O =
fugacityCoefficientH2O(params, temperature, pg, yH2O, highTemp, extrapolate, spycherPruess2005);
732 Evaluation gammaH2O = activityCoefficientH2O_(temperature, xCO2, highTemp);
736 if ( temperature > 372.15 && temperature < 382.15 && !spycherPruess2005) {
737 const Evaluation weight = (382.15 - temperature) / 10.;
738 const Evaluation& k0_H2O_low = equilibriumConstantH2O_(temperature,
false);
739 const Evaluation& phi_H2O_low =
fugacityCoefficientH2O(params, temperature, pg, Evaluation(0.0),
false, extrapolate);
740 k0_H2O = k0_H2O * (1 - weight) + k0_H2O_low * weight;
741 phi_H2O = phi_H2O * (1 - weight) + phi_H2O_low * weight;
745 const Evaluation& pg_bar = pg / 1.e5;
747 return k0_H2O * gammaH2O / (phi_H2O * pg_bar) * exp(deltaP * v_av_H2O / (R * temperature));
758 template <
class Evaluation,
class CO2Parameters>
759 static Evaluation computeB_(
const CO2Parameters& params,
760 const Evaluation& temperature,
761 const Evaluation& pg,
762 const Evaluation& yH2O,
763 const Evaluation& xCO2,
764 const bool& highTemp,
765 bool extrapolate =
false,
766 bool spycherPruess2005 =
false)
768 OPM_TIMEFUNCTION_LOCAL();
770 const Evaluation& deltaP = pg / 1e5 - Pref_(temperature, highTemp);
771 Evaluation v_av_CO2 = V_avg_CO2_(temperature, highTemp);
772 Evaluation k0_CO2 = equilibriumConstantCO2_(temperature, pg, highTemp, spycherPruess2005);
773 Evaluation phi_CO2 =
fugacityCoefficientCO2(params, temperature, pg, yH2O, highTemp, extrapolate, spycherPruess2005);
774 Evaluation gammaCO2 = activityCoefficientCO2_(temperature, xCO2, highTemp);
778 if ( temperature > 372.15 && temperature < 382.15 && !spycherPruess2005) {
779 const Evaluation weight = (382.15 - temperature) / 10.;
780 const Evaluation& k0_CO2_low = equilibriumConstantCO2_(temperature, pg,
false, spycherPruess2005);
781 const Evaluation& phi_CO2_low =
fugacityCoefficientCO2(params, temperature, pg, Evaluation(0.0),
false, extrapolate);
782 k0_CO2 = k0_CO2 * (1 - weight) + k0_CO2_low * weight;
783 phi_CO2 = phi_CO2 * (1 - weight) + phi_CO2_low * weight;
787 const Evaluation& pg_bar = pg / 1.e5;
789 return phi_CO2 * pg_bar / (55.508 * k0_CO2 * gammaCO2) * exp(-deltaP * v_av_CO2 / (R * temperature));
795 template <
class Evaluation>
796 OPM_HOST_DEVICE
static Evaluation activityCoefficientSalt_(
const Evaluation& temperature,
797 const Evaluation& pg,
798 const Evaluation& m_NaCl,
799 const Evaluation& xCO2,
800 const int& activityModel)
802 OPM_TIMEFUNCTION_LOCAL();
809 if (activityModel == 1) {
810 lambda = computeLambdaRumpfetal_(temperature);
812 Evaluation m_CO2 = xCO2 * (2 * m_NaCl + 55.508) / (1 - xCO2);
813 convTerm = (1 + (m_CO2 + 2 * m_NaCl) / 55.508) / (1 + m_CO2 / 55.508);
815 else if (activityModel == 2) {
816 lambda = computeLambdaSpycherPruess2009_(temperature);
817 xi = computeXiSpycherPruess2009_(temperature);
818 convTerm = 1 + 2 * m_NaCl / 55.508;
820 else if (activityModel == 3) {
821 lambda = computeLambdaDuanSun_(temperature, pg);
822 xi = computeXiDuanSun_(temperature, pg);
826 throw std::runtime_error(
"Activity model for salt-out effect has not been implemented!");
830 const Evaluation& lnGamma = 2 * lambda * m_NaCl + xi * m_NaCl * m_NaCl;
833 return convTerm * exp(lnGamma);
839 template <
class Evaluation>
840 OPM_HOST_DEVICE
static Evaluation computeLambdaSpycherPruess2009_(
const Evaluation& temperature)
843 static const Scalar c[3] = { 2.217e-4, 1.074, 2648. };
846 return c[0] * temperature + c[1] / temperature + c[2] / (temperature * temperature);
852 template <
class Evaluation>
853 OPM_HOST_DEVICE
static Evaluation computeXiSpycherPruess2009_(
const Evaluation& temperature)
856 static const Scalar c[3] = { 1.3e-5, -20.12, 5259. };
859 return c[0] * temperature + c[1] / temperature + c[2] / (temperature * temperature);
865 template <
class Evaluation>
866 OPM_HOST_DEVICE
static Evaluation computeLambdaRumpfetal_(
const Evaluation& temperature)
869 static const Scalar c[4] = { 0.254, -76.82, -10656, 6312e3 };
871 return c[0] + c[1] / temperature + c[2] / (temperature * temperature) +
872 c[3] / (temperature * temperature * temperature);
882 template <
class Evaluation>
883 OPM_HOST_DEVICE
static Evaluation computeLambdaDuanSun_(
const Evaluation& temperature,
const Evaluation& pg)
885 static const Scalar c[6] =
886 { -0.411370585, 6.07632013E-4, 97.5347708, -0.0237622469, 0.0170656236, 1.41335834E-5 };
888 Evaluation pg_bar = pg / 1.0E5;
889 return c[0] + c[1]*temperature + c[2]/temperature + c[3]*pg_bar/temperature + c[4]*pg_bar/(630.0 - temperature)
890 + c[5]*temperature*log(pg_bar);
900 template <
class Evaluation>
901 OPM_HOST_DEVICE
static Evaluation computeXiDuanSun_(
const Evaluation& temperature,
const Evaluation& pg)
903 static const Scalar c[4] =
904 { 3.36389723E-4, -1.98298980E-5, 2.12220830E-3, -5.24873303E-3 };
906 Evaluation pg_bar = pg / 1.0E5;
907 return c[0] + c[1]*temperature + c[2]*pg_bar/temperature + c[3]*pg_bar/(630.0 - temperature);
916 template <
class Evaluation>
917 static Evaluation equilibriumConstantCO2_(
const Evaluation& temperature,
918 const Evaluation& pg,
919 const bool& highTemp,
920 bool spycherPruess2005 =
false)
922 OPM_TIMEFUNCTION_LOCAL();
923 Evaluation temperatureCelcius = temperature - 273.15;
924 std::array<Scalar, 4> c;
926 c = { 1.668, 3.992e-3, -1.156e-5, 1.593e-9 };
937 c = { 1.169, 1.368e-2, -5.38e-5, 0.0 };
940 c = { 1.189, 1.304e-2, -5.446e-5, 0.0 };
943 Evaluation logk0_CO2 = c[0] + temperatureCelcius * (c[1] + temperatureCelcius *
944 (c[2] + temperatureCelcius * c[3]));
945 Evaluation k0_CO2 = pow(10.0, logk0_CO2);
955 template <
class Evaluation>
956 OPM_HOST_DEVICE
static Evaluation equilibriumConstantH2O_(
const Evaluation& temperature,
const bool& highTemp)
958 Evaluation temperatureCelcius = temperature - 273.15;
959 std::array<Scalar, 5> c;
961 c = { -2.1077, 2.8127e-2, -8.4298e-5, 1.4969e-7, -1.1812e-10 };
964 c = { -2.209, 3.097e-2, -1.098e-4, 2.048e-7, 0.0 };
966 Evaluation logk0_H2O = c[0] + temperatureCelcius * (c[1] + temperatureCelcius * (c[2] +
967 temperatureCelcius * (c[3] + temperatureCelcius * c[4])));
968 return pow(10.0, logk0_H2O);