Connecting element-based flow and transport solver to PHREEQC

[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: [Select]

            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

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: [Select]

    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

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

Code: [Select]

----------------------------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

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!