PQI script for a specific chemical system

[Sara Tabrizi]

Hello Dear Dr. Parkhurst,
I am sorry for the long question. But it has been a while since I am stuck with it, so any help is deeply appreciated. I would like to simulate a reactive transport problem which is published in the paper by Ahusborde et al 2018 (https://doi.org/10.2516/ogst/2018033) in the application of CO2 sequestration. The transport code is coupled with phreeqcrm. The reaction system is as follows:

H2O(g) = H2O(l)
CO2(g) = CO2(l)
HCO3- = CO2(l) + H2O(l) - H+
CO32- = CO2(l) + H2O(l) – 2H+
OH- = H2O(l) – H+
Anorthite = 4H2O(l)  + Ca2+  + 2Al3+ + 2SiO2(l) – 8H+
Calcite = Ca2+ + HCO3- - H+
Kaolinite = 5H2O(l) + 2Al3+ + 2SiO2(l) - 6H+

The first five reactions are considered in equilibrium phase and three mineral dissolutions as kinetic reactions and the chemical species defined in the system consist of six primary species (H2O(l), CO2(l), H+, Ca2+ , Al3+  SiO2(l) ) and eight secondary species (H2O(g), CO2(g), OH-, CO32-, HCO3-, Anorthite, Calcite, Kaolinite). Therefore, there are 14 species in this system. The following is my phreeqc script in the pqi file to model the exact chemical system that has been mentioned above. However, with this script, I received 21 species and I cannot find the reactions corresponding to all these species. When I try to print the concentration of these extra species, the values are not negligible. Therefore, it is affecting the results and the transport code stops converging in the second time step after the first reactive transport step. Because in the transport code, only the transport of species in the above-mentioned chemical system is taken into account.  I would be really grateful if you could let me know what to change in this script to receive exactly the same chemical system as what I mentioned above.


SOLUTION_MASTER_SPECIES
    H           H+              -1.0     1                1
    O           H2O             0        18             16
    C           CO3-2           2.0   HCO3        60                12
    Ca          Ca+2            0        40.08             40.08
    Al          Al+3             0        27               27
    Si    H4SiO4             0      SiO2      28.0843

SOLUTION_SPECIES
H+ = H+
  log_k     0
  -gamma   9.0   0

H4SiO4=H4SiO4
  log_k     0
  -gamma   9.0   0

H2O = H2O
  log_k     0
  -gamma    1000000000 0

Al+3 = Al+3
  log_k     0
  -gamma    1000000000 0

CO3-2 = CO3-2
  log_k     0
 -gamma   5.4   0

Ca+2 = Ca+2
  log_k     0
  -gamma   5.0   0.1650

H2O = H+ + OH-
  log_k     -13.26
  -gamma   3.5   0

H4SiO4=SiO2 + 2 H2O
   -log_k   2.71

 CO2 + H2O = HCO3- + H+
  -log_k   -10.23


 H+ + CO3-2 =HCO3-
  log_k     -6.32

 SOLUTION 0
units  mol/kgw
temp             25.0
pH 6.24 charge
Al  1.28E-9
C   0.00569
Ca  0.00252
Si  4.73E-4

Phases
Calcite
   CaCO3 + H+ = Ca+2 + HCO3-
   log_k   1.6

Anorthite
   CaAl2Si2O8 + 8H+ =4H2O + Ca+2 + 2Al+3 + 2SiO2
   log_k   25.82

 Kaolinite
   Al2Si2O5(OH)4+ 6H+ = 5H2O + 2Al+3 + 2SiO2
   log_k   6.82


RATES
    Calcite
-start

  1   REM   PARM(1) = specific surface area of calcite, 1/dm calcite
  2   REM   PARM(2) = exponent for M/M0
 10  si_cc = SR("Calcite")
 20  IF (M <= 0  and si_cc < 0) THEN GOTO 200

 80  IF M0 > 0 THEN area = PARM(1)*(M/M0)^PARM(2) ELSE area = PARM(1)*M

130 rate = area * 1.58e-9

140 rate = rate * (1 - SR("Calcite"))
150 moles = rate  * TIME
200 SAVE moles
-end

    Anorthite
-start

  1   REM   PARM(1) = specific surface area of calcite, 1/dm calcite
  2   REM   PARM(2) = exponent for M/M0
 10  si_cc = SR("Anorthite")
 20  IF (M <= 0  and si_cc < 0) THEN GOTO 200

 80  IF M0 > 0 THEN area = PARM(1)*(M/M0)^PARM(2) ELSE area = PARM(1)*M

130 rate = area * 1.0e-12

140 rate = rate * (1 - SR("Anorthite"))
150 moles = rate  * TIME
200 SAVE moles
-end


   Kaolinite
-start
  1   REM   PARM(1) = specific surface area of calcite, 1/dm calcite
  2   REM   PARM(2) = exponent for M/M0
 10  si_cc = SR("Kaolinite")
 20  IF (M <= 0  and si_cc < 0) THEN GOTO 200

 80  IF M0 > 0 THEN area = PARM(1)*(M/M0)^PARM(2) ELSE area = PARM(1)*M

130 rate = area * 1.0e-13

140 rate = rate * (1 - SR("Kaolinite"))
150 moles = rate  * TIME
200 SAVE moles
-end
END

USE SOLUTION 0
KINETICS 1
Calcite
    -m0       238
    -parms    8.8 1
    -tol      1e-8

KINETICS 2
Anorthite
    -m0       87
    -parms    8.8 1
    -tol      1e-8


KINETICS 3
Kaolinite
    -m0       88
    -parms    1760 1
    -tol      1e-8

[dlparkhurst]

I don't know about that paper. The reaction system seems a little simplistic.

As far as using PHREEQC, there is a minimum set of species that must be included: H2O, H+, OH-, O2, H2, and e-. I think you will have to make a new database that contains just the species and phases that you want. If you use one of the distributed databases, you will get a lot more species than are specified in the paper. Here is a first start on a minimum database that copies data from phreeqc.dat. For simplicity, you can probably remove the -Dw, -Vm and analytical expressions lines.

Code: [Select]

SOLUTION_MASTER_SPECIES
H H+ -1.0 H 1.008
H(0) H2 0 H
H(1) H+ -1.0 0           1
E e- 0 0.0 0
O H2O 0 O 16.0
O(0) O2 0 O
O(-2) H2O 0 0
C CO3-2 2.0 HCO3 12.0111       60                12
Ca          Ca+2            0        40.08             40.08
Al          Al+3             0        27               27
Si          H4SiO4             0      SiO2      28.0843
SOLUTION_SPECIES
H+ = H+
-gamma 9.0 0
-dw 9.31e-9  1000  0.46  1e-10 # The dw parameters are defined in ref. 3.
# Dw(TK) = 9.31e-9 * exp(1000 / TK - 1000 / 298.15) * viscos_0_25 / viscos_0_tc
# Dw(I) = Dw(TK) * exp(-0.46 * DH_A * |z_H+| * I^0.5 / (1 + DH_B * I^0.5 * 1e-10 / (1 + I^0.75)))
e- = e-
H2O = H2O
# H2O + 0.01e- = H2O-0.01; -log_k -9 # aids convergence
Ca+2 = Ca+2
-gamma 5.0 0.1650
-dw 0.793e-9  97  3.4  24.6
-Vm  -0.3456  -7.252  6.149  -2.479  1.239  5  1.60  -57.1  -6.12e-3  1 # The apparent volume parameters are defined in ref. 1 & 2
Al+3 = Al+3
-gamma 9.0 0
-dw 0.559e-9
-Vm   -2.28  -17.1  10.9  -2.07  2.87  9  0  0  5.5e-3  1 # ref. 2 and Barta and Hepler, 1986, Can. J.C. 64, 353.
H4SiO4 = H4SiO4
-dw 1.10e-9
-Vm  10.5  1.7  20  -2.7  0.1291 # supcrt + 2*H2O in a1
H2O = OH- + H+
-analytic  293.29227  0.1360833  -10576.913  -123.73158  0  -6.996455e-5
-gamma 3.5 0
-dw 5.27e-9  548  0.52  1e-10
-Vm  -9.66  28.5  80.0 -22.9 1.89 0 1.09 0 0 1
2 H2O = O2 + 4 H+ + 4 e-
-log_k -86.08
-delta_h 134.79 kcal
-dw 2.35e-9
-Vm  5.7889  6.3536  3.2528  -3.0417  -0.3943 # supcrt
2 H+ + 2 e- = H2
-log_k -3.15
-delta_h -1.759 kcal
-dw 5.13e-9
-Vm 6.52  0.78  0.12 # supcrt
CO3-2 = CO3-2
-gamma 5.4 0
        -dw 0.955e-9  28.9  14.3  98.1
-Vm  8.69  -10.2  -20.31  -0.131  4.65  0  3.75  0  -4.04e-2  0.678
CO3-2 + H+ = HCO3-
-log_k 10.329
-delta_h -3.561 kcal
-analytic  107.8871  0.03252849  -5151.79  -38.92561  563713.9
-gamma 5.4      0
-dw 1.18e-9  -182  0.351  -4.94
-Vm  9.03  -7.03e-2  -13.38  0  2.05  0  0  128  0  0.8242
CO3-2 + 2 H+ = CO2 + H2O
-log_k 16.681
-delta_h -5.738 kcal
-analytic 464.1965 0.09344813 -26986.16 -165.75951 2248628.9
-dw 1.92e-9  -120 # TK dependence from Cadogan et al. 2014, , JCED 59, 519
-Vm   7.29  0.92  2.07  -1.23  -1.60 # McBride et al. 2015, JCED 60, 171
-gamma 0 0.066 # Rumpf et al. 1994, J. Sol. Chem. 23, 431
END



[Sara Tabrizi]

Hello again Dear Dr. Parkhurst,
Thank you for your reply. I took your advice and modified the database, so the number of species was reduced. However, as you mentioned, the existence of H2 and O2 is inevitable and since they do not exist in the reactive system that I introduced in the previous question, I think the obtained concentration for them should be zero. With the reactive input file that was introduced in the previous question, I have received zero concentration for O2, but for H2 the concentration is not zero and the values are not negligible. I think this is affecting the concentration of other species such as H+ which is a primary species considered in the transport system and therefore, the transport system stops converging. I would really appreciate it if you could help me understand how to obtain zero concentration for H2 and hopefully the right values for primary and secondary species considered in the system.
Regrds,
Sara Tabrizinejadas

[dlparkhurst]

H2 and O2 are the only redox species in your system, so either your initial conditions produce H2, or some of your reactions are incorrect.  H2 should not be produced by any of the mineral or solution reactions.

First check whether there is H2 in your initial conditions. If it is present, adjust the pe  to a point where there is neither O2 or H2.

Then equilibrate each of your minerals with pure water. If one of them produces H2, there is an error in your reactions.

I'm a little suspicious of your handling of H4SiO4. SOLUTION_MASTER_SPECIES had H4SiO4 as the master species, but the minerals had SiO2 in the reactions. You should only use SiO2 or H4SiO4 as an aqueous species, not both. Use the reaction

[Sara Tabrizi]

Hello Dear Dr. Parkhurst,
Thank you very much for replying to my questions. I have looked at the initial situations and there is no problem of generating redox species wrongly there. I think that the problem with that is that I did not apply the injection of CO2 properly to the system in the phreeqc input file that I included in the first question. When I include CO2(g) with the following lines the values of species concentration are getting closer to the right results.
EQUILIBRIUM_PHASES 0-400
CO2(g) -2 10

However, it is not still the correct result. I think the problem is that I should introduce CO2 directly in the aqueous phase. That is because I have decoupled the two-phase and one-phase flow problems and therefore, the multicomponent reactive transport problem in which the chemistry part is solved with Phreeqc is a single-phase problem and I do not think that I should equilibrate the gaseous CO2. But I should include the dissolved aqueous CO2 which is introduced to the one-phase reactive problem from the two-phase problem. I would appreciate it if you could help me with the expression that exists in the phreeqc batch file for including the aqueous CO2 instead of the gaseous CO2.

[dlparkhurst]

When you use EQUILIBRIUM_PHASES, there is no explicit gas phase. You are specifying a fixed partial pressure for CO2(g) in the solution. In response to any other reactions, CO2 will be added or removed from solution to maintain the fixed partial pressure of the solution.  You will have to consider whether this is the correct approach for your calculations. The approach would not be appropriate if the solution is (1) not in contact with a gas phase, or (2) if the partial pressure of CO2 is variable. In case (1), you should not include CO2 in EQUILIBRIUM_PHASES for cells not in contact with a gas. In case (2), you would need to either adjust the target SI for CO2(g) at each iteration with EQUILIBRIUM_PHASES_MODIFY, or use a GAS_PHASE definition where the moles of gas is adjusted at each iteration.

I don't think that the handling of CO2 should have been the cause of the accumulation of H2 in solution. Unfortunately, it is difficult to debug in your reactive transport code, which has many sources of complication. If you have problems, I suggest you develop a simple PHREEQC transport problem to investigate issues. Once you are satisfied with the PHREEQC calculation, you can reproduce it in your reactive transport model before expanding to your complete simulation.

[Sara Tabrizi]

Thank you very much for your fast reply. I think that my model is considered as the second case that you mentioned with the variable partial gas pressure. That is why I have tried to use the EQUILIBRIUM_PHASES_MODIFY keyword. For the chemical system, I have tried to simplify the chemical system and I am using the following system in which the only mineral is CaCO3 and equilibrium dissolution is considered for that with the injection of CO2. Here is the simplified chemical system:
HCO3- = CO2(l) + H2O(l) - H+
OH- = H2O(l) – H+
Calcite = Ca2+ + HCO3- - H+

Initial liquid CO2 is For the pqi batch file, I am using the following:
SOLUTION 0-400
    pH        7.38
    units     mol/L
    C         0.0001223
    Ca        1.708E-3

EQUILIBRIUM_PHASES 0-400
Calcite 0.0 24.3
EQUILIBRIUM_PHASES_MODIFY 0-400
-component      CO2(g)
-si         0
-moles      0.0001223

However, I am not sure if this input file is precisely describing an equilibrium calcite dissolution due to injection of CO2 to the domain. Is there a way that I can precisely define the injection rate for CO2? I have difficulty understanding how to convert the inputs about initial value and injection rate of CO2 into EQUILIBRIUM_PHASES_MODIFY keyword correctly.  I would really appreciate if you could help me about it.
Regards,
Sara Tabrizinejadas

[dlparkhurst]

I don't quite understand how you want to do this calculation.

You can add CO2 kinetically at a rate by using KINETICS and RATES.

You can add a specified amount of CO2 with REACTION.

You can use EQUILIBRIUM_PHASES to fix a constant partial pressure subject to an amount of CO2 that is available to dissolve.

You can equilibrate between a solution and a GAS_PHASE at a fixed total pressure (or total gas volume).

I still don't trust the paper that you cited. Tables 6 and 8 do not seem to be consistent and the reaction network is not realistic (Al(OH)4- tends to be a dominant Al species). Still, if you can figure out what they did, I guess you can use it as a test case.

[Sara Tabrizi]

Hello again Dr.Parkhurst,
Thank you very much for your response. About the inconsistency between the tables, I agree with you. I have contacted the authors of that paper and Al(OH)2+  in this network of reaction does not exist and it was written in table 8 by mistake. The reactive system in this paper is also derived from a known reactive network in the paper: https://doi.org/10.1016/j.advwatres.2012.03.012
The simple reactive network in this paper is what I am trying to model. My question is regarding the input of CO2 injection during the time in Phreeqc. It is assumed that a total mass of 18.9E9 kg of CO2 is injected during a 20-year period and an initial molality of CO2 in the liquid phase is 3.55E-3 mol/kg. I was wondering, among the propositions you made for inputting the injection of CO2 in Phreeqc, which one would be the most accurate?
Thank you in advance,
Sara

[dlparkhurst]

You could add CO2 at a constant rate for 20 years with KINETICS and RATES. You would probably define this at a cell or cells that represent the injection location.

As for a gas phase, I don't know if you are accounting for gas transport. For starters, without gas transport, you could simply use EQUILIBRIUM_PHASES including CO2(g) at a fixed pressure with zero or some non-zero initial moles.  The CO2 that you add with kinetics will distribute between the solution and EQUILIBRIUM_PHASES. If you consume enough CO2 with mineral reactions, the EQUILIBRIUM_PHASES CO2(g) could go to zero.

It sounds like the reports you cite have fully coupled gas transport. If you are trying to develop a multiphase flow model, you probably need to use GAS_PHASE definitions. In this case, assuming you are using a sequential approach of liquid transport,  gas transport, and reactions, you would update the moles in the GAS_PHASE after gas transport. I have no experience with this approach. It seems difficult because the gas transport could be quite fast, requiring small time steps.

I have no doctorate, so no Dr. please.

[Sara Tabrizi]

Thank you, sir. Sorry for my mistake.

[Sara Tabrizi]

Hello again dear sir Parkhurst,
I have simplified the chemical system and I am only considering the dissolution of Calcite. I have succeded in getting good results in Phreeqcrm coupled with my transport code for both kinetic and equilibrium dissolution of Calcite. However, for kinetic dissolution, the simulation is very slow. I am not providing any time step for the integration of kinetics in the pqi file. So the time step for transport is given to phreeqcrm with the function SetTimeStep(timestep). I think that Phreeqc is integrating the kinetic reactions over this timestep. My question is if this time step is divided into smaller time steps in Phreeqc by default? because I was not expecting to have such a long simulation time, so I am suspicious that maybe Phreeqc is dividing the time step for kinetic calculation by default. Thank you very much again for your help.
Regards,
Sara

[Sara Tabrizi]

To better explain my problem, Following is my Phreeqc script for this simple chemical system:

(1) H2O=OH- + H+
(2) CO2 + H2O = HCO3- + H+
(3) HCO3- = H+ + CO3--
(4) CaCO3 + H+ = HCO3- + Ca++

The fourth reaction is considered kinetically in the following phreeqc script. 
However, as I mentioned in the previous question, the simulation is very slow and I am receiving this kind of warnings:
WARNING: Negative moles in solution xx for C, -1.026552e-03. Recovering...
WARNING: Negative moles in solution xx for Ca, -2.703841e-05. Recovering...


[
SOLUTION_MASTER_SPECIES
    H           H+              -1.0     1                1
    O           H2O             0        18             16
    C           CO3-2           2.0   HCO3        60                12
    Ca          Ca+2            0        40.08             40.08

SOLUTION_SPECIES
H+ = H+
  log_k     0
  -gamma   9.0   0

H2O = H2O
  log_k     0
  -gamma    1000000000 0

CO3-2 = CO3-2
  log_k     0
  -gamma   5.4   0

Ca+2 = Ca+2
  log_k     0
  -gamma   5.0   0.1650

 H2O = H+ + OH-
  log_k     -13.95
  -gamma   3.5   0

CO2 + H2O = HCO3- + H+
  -log_k   -6.293

HCO3- = H+ + CO3-2
  log_k     -10.279

 Phases
    Calcite
      CaCO3  + H+ = Ca+2 + HCO3-
      log_k   -1.899

SOLUTION 0-400
units  mol/kgw
temp             25.0
pH 9.8844
C   0.00012051043469
Ca  0.0001204690124
END

RATES
    Calcite
-start
  1   REM   PARM(1) = specific surface area of calcite, 1/dm calcite
  2   REM   PARM(2) = exponent for M/M0
 10  si_cc = SR("Calcite")
 20  IF (M <= 0  and si_cc < 0) THEN GOTO 200
 80  IF M0 > 0 THEN area = PARM(1)*(M/M0)^PARM(2) ELSE area = PARM(1)*M
 130 rate = area * 1.6e-09
 140 rate = rate * (1 - SR("Calcite"))
 150 moles = rate  * TIME
 200 SAVE moles
-end

USE SOLUTION 0-400
KINETICS 1-400
Calcite
-m       10
-m0       0.1
-parms    8.8 1
-tol      1e-08


SELECTED_OUTPUT 1
-reset false
-ph                   true
-kin Calcite
]

[dlparkhurst]

Normally, M0 is greater than or equal to M, and the M/M0 factor is used to account for changing surface area as the amount of reactant decreases. You have M0 100 times smaller than M, which both increases the rate and apparently makes the calculation more unstable because small changes in reaction cause exaggerated changes in rate. If you really love that value for M0, use the -cvode option. The implicit method is better at anticipating how the rate will change due to reaction.

With your rate, calcite equilibrates in about 5 seconds. In the laboratory, calcite equilibrates in hours to days. If your simulation time is long relative to kinetic equilibration, you are better to include Calcite in EQUILIBRIUM_PHASES. You calculation will be more robust and will run faster.

Note that the Warnings are just warnings. Unless the calculation stops with an Error, the results should be correct even if the kinetic calculation issues a lot of warnings.

[Sara Tabrizi]

Hi again,
Thanks for correcting my mistake about m0. I have tried with m0=m=10 and still have the same warning with the same speed of simulation. Since I would like include other minerals with slower dissolution in the next codes, I would like to use Kinetics and my question is if the kinetics are exactly integrated over the input of "SetTimeStep" function in Phreeqcrm or this value is divided by a default number to increase the precision?
Thank you in advance,
Sara

[dlparkhurst]

I can only help you with a PHREEQC input file; things get too complicated with PhreeqcRM and another transport code. I get no warnings when I run the file you posted using M0=10 and -cvode. -cvode can be slower or faster than the default Runge-Kutta method, you will just have to experiment.

Kinetic reactions are integrated for the entire time step defined in SetTimeStep for each cell that has KINETICS. The integration algorithm uses at least 6 sub time steps in the integration process and perhaps more if a convergence criterion is not met. You have a definition of 400 solutions in your file (I would think only one copy of the solution is needed for PhreeqcRM). I don't think that kinetics on 400 cells would be too slow, especially if you are using the parallelization of PhreeqcRM. If you make a PHREEQC input file, I will look at it, although 400 may a lot for PHREEQC's transport method depending on the cell size and dispersion. You should probably compare your PhreeqcRM run to equivalent PHREEQC runs in any event (starting with a simple case).

[Sara Tabrizi]

Thanks you very much for the advice. For Phreeqcrm also -cvode option helps a lot, it becomes much faster.