Connecting element-based flow and transport solver to PHREEQC

Beginners > PHREEQC basics

(1/1)

anovikov:
Hello everyone,

I am modeling the injection of carbonated water into 100% limestone (calcite) by connecting my solver for reactive flow and transport in porous media to PHREEQC. PHREEQC is supposed to do geochemical modeling within pore space only, while I take care myself of the dissolution/precipitation kinetics based on PHREEQC output. To reduce the number of mass balance equations my solver is element-based, i.e. solving mass balance equations written for elements. So, my unknowns are pressure and molar fractions of elements. Therefore, I specify moles of elements as input for PHREEQC along with given pressure and water mass.

My input logic is the following. Given molar fractions of elements, i mix them in 1000 moles of elements in the following order. First, I calculate the water mass consuming hydrogen and oxygen until one runs out. Second, I provide the remaining elements as input to the REACTION keyword in moles.

--- Code: --- USER_PUNCH
-headings H(mol) O(mol) C(mol) Ca(mol) Vol_aq SI SR ACT("H+") ACT("CO2") ACT("H2O") {species_headings}
10 PUNCH TOTMOLE("H") TOTMOLE("O") TOTMOLE("C") TOTMOLE("Ca") SOLN_VOL SI("Calcite") SR("Calcite") ACT("H+") ACT("CO2") ACT("H2O") {species_punch}

SELECTED_OUTPUT
-selected_out true
-user_punch true
-reset false
-high_precision true
-gases CO2(g) H2O(g)

SOLUTION 1
temp {{temperature:.2f}}
pressure {{pressure:.4f}}
pH 7 charge
-water {{water_mass:.10f}} # kg
REACTION 1
H {{hydrogen:.10f}}
O {{oxygen:.10f}}
C {{carbon:.10f}}
Ca {{calcium:.10f}}
1
KNOBS
-convergence_tolerance 1e-10
END

--- End code ---

The problem is that the following input calculated for about H2O: 98.5%, CO2: 1.5% molar fractions and 100 bar pressure, i.e. CO2 must be completely dissolved in water

--- Code: --- USER_PUNCH
headings H(mol) O(mol) C(mol) Ca(mol) Vol_aq SI SR ACT("H+") ACT("CO2") ACT("H2O") MOL("OH-") MOL("H+") MOL("H2O") MOL("C(-4)") MOL("CH4") MOL("C(4)") MOL("HCO3-") MOL("CO2") MOL("CO3-2") MOL("CaHCO3+") MOL("CaCO3") MOL("(CO2)2") MOL("Ca+2") MOL("CaOH+") MOL("H(0)") MOL("H2") MOL("O(0)") MOL("O2")
10 PUNCH TOTMOLE("H") TOTMOLE("O") TOTMOLE("C") TOTMOLE("Ca") SOLN_VOL SI("Calcite") SR("Calcite") ACT("H+") ACT("CO2") ACT("H2O") MOL("OH-") MOL("H+") MOL("H2O") MOL("C(-4)") MOL("CH4") MOL("C(4)") MOL("HCO3-") MOL("CO2") MOL("CO3-2") MOL("CaHCO3+") MOL("CaCO3") MOL("(CO2)2") MOL("Ca+2") MOL("CaOH+") MOL("H(0)") MOL("H2") MOL("O(0)") MOL("O2")
SELECTED_OUTPUT
selected_out true
user_punch true
reset false
high_precision true
gases CO2(g) H2O(g)
SOLUTION 1
temp 323.15
pressure 98.6923
pH 7 charge
water 5.9200054446 # kg
REACTION 1
H 0.0000000000
O 9.4724565506
C 4.7362282753
Ca 0.0000000100
1
KNOBS
convergence_tolerance 1e-10
PRINT
species true
totals true
END

--- End code ---

out of sudden results in high O2 molality, which is unphysical:

--- Code: -------------------------------Distribution of species----------------------------

Log Log Log mole V
Species Molality Activity Molality Activity Gamma cm?/mol

H+ 2.940e-005 2.887e-005 -4.532 -4.539 -0.008 0.00
OH- 3.288e-008 3.228e-008 -7.483 -7.491 -0.008 -349.34
H2O 5.551e+001 1.523e-001 1.744 -0.817 0.000 26.33
C(-4) 0.000e+000
CH4 0.000e+000 0.000e+000 -80.891 -80.891 0.000 315.32
C(4) 8.000e-001
CO2 7.233e-001 7.233e-001 -0.141 -0.141 0.000 330.42
(CO2)2 3.838e-002 3.838e-002 -1.416 -1.416 0.000 660.84
HCO3- 2.937e-005 2.884e-005 -4.532 -4.540 -0.008 -256.58
CaHCO3+ 1.689e-009 1.658e-009 -8.772 -8.780 -0.008 -44.15
CO3-2 6.856e-012 6.372e-012 -11.164 -11.196 -0.032 -1254.26
CaCO3 3.852e-025 3.852e-025 -24.414 -24.414 0.000 -8.31
Ca 1.689e-009
CaHCO3+ 1.689e-009 1.658e-009 -8.772 -8.780 -0.008 -44.15
Ca+2 4.606e-013 4.281e-013 -12.337 -12.368 -0.032 -241.99
CaOH+ 2.323e-022 2.281e-022 -21.634 -21.642 -0.008 (0)
CaCO3 3.852e-025 3.852e-025 -24.414 -24.414 0.000 -8.31
H(0) 2.804e-024
H2 1.402e-024 1.402e-024 -23.853 -23.853 0.000 28.50
O(0) 9.820e+001
O2 4.910e+001 4.910e+001 1.691 1.691 0.000 106.74

--- End code ---

I appreciate if someone could give me a hint on how to fix either the input logic or the PHREEQC request to get more meaningful result.

Best,
Aleks

dlparkhurst:
Pe of 4 at 300 C puts you outside the stability field of water for SOLUTION 1. It produces 98 moles of O2(aq). Use a lower pe.

Better yet, use PhreeqcRM to do your calculations, instead of PHREEQC (or IPhreeqc). PhreeqcRM is designed to be coupled with flow and transport models. One of the transport units that it accepts is mole fraction. So, you can use a method SetConcentrations with mole fractions, and PhreeqcRM will convert to PHREEQC concentrations, run reactions, and return mole fraction concentrations with the method GetConcentrations.

dlparkhurst:
Oops, I meant mass fraction instead of mole fraction.

anovikov:
I set the wrong temperature in Kelvins instead of Celcius. Stupid mistake.
Regarding coupling, phreeqpy from Python is very convenient, but it would indeed be interesting to use PhreeqcRM.

Thanks a lot, David!

Navigation

[0] Message Index

Go to full version