H2O is zero in solubility of CO2

[Mohamadreza]

I used a multiphase reactive transport package at the pore scale, coupled with PHREEQC (phreeqcRM, IPhreeqc, IPhreeqcPhast, etc.), to simulate the solubility of supercritical CO₂ in a 2 molal CaCl₂ salt solution at a pressure of 10 MPa and a temperature of 60 ?C. The geometry of the model is at the microscale, with very small cells and a high concentration of salt ions. However, when running the simulation, I encountered the following error:
"A(H₂O): Activity of water has not converged."
Notably, when I ran the simulation in pure water or simulate with cacl2 solution with (Dilute calcium chloride solution)0.01 molal, no such error occurred.
Could you please advise on how to resolve this issue?

Code: [Select]

SOLUTION 0 formation water
pressure 100
temp 60
water 1 # kg
pe 4
redox pe
units   Mol/L
C 1e-13
Ca 2
Cl 4
KNOBS
    diagonal_scale true
    iterations         2000
    convergence_tolerance   1e-12

SOLUTION_MASTER_SPECIES
#
#element species alk gfw_formula element_gfw
#
H H+ -1.0 H 1.008
H(0) H2 0 H
H(1) H+ -1.0 0
E e- 0 0.0 0
O H2O 0 O 16.0
O(0) O2 0 O
O(-2) H2O 0 0
Ca Ca+2 0 Ca 40.08
Cl Cl- 0 Cl 35.453
C CO3-2 2.0 HCO3 12.0111
C(+4) CO3-2 2.0 HCO3
Alkalinity CO3-2 2.0 61.0173 61.0173

SOLUTION_SPECIES
H+ = H+
-gamma 9.0 0

-dw 9.31e-9
e- = e-
H2O = H2O
-dw 2.299e-9  -254  1e-10  1e-10
H2O + 0.01e- = H2O-0.01
        -log_k -9 # aids convergence
Ca+2 = Ca+2
-gamma 5.0 0.1650
-dw 0.793e-9
-Vm  -0.3456  -7.252  6.149  -2.479  1.239  5  1.60  -57.1  -6.12e-3  1 # ref. 1
Cl- = Cl-
-gamma 3.5   0.015
-gamma 3.63  0.017 # cf. pitzer.dat
-dw 2.03e-9
-Vm  4.465  4.801  4.325  -2.847  1.748  0  -0.331  20.16  0  1 # ref. 1
CO3-2 = CO3-2
-gamma 5.4 0
-dw 0.955e-9
-Vm  5.95  0  0  -5.67  6.85  0  1.37  106  -0.0343  1 # ref. 1
# aqueous species
H2O = OH- + H+
-analytic  293.29227  0.1360833  -10576.913  -123.73158  0  -6.996455e-5
-gamma 3.5 0
-dw 5.27e-9
-Vm  -9.66  28.5  80.0 -22.9 1.89 0 1.09 0 0 1 # ref. 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 + 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
-Vm  8.472  0  -11.5  0  1.56  0  0  146  3.16e-3  1 # ref. 1
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
-Vm   7.29  0.92  2.07  -1.23  -1.60 # ref. 1 + McBride et al. 2015, JCED 60, 171
2CO2 = (CO2)2 # activity correction for CO2 solubility at high P, T
-log_k -1.8
-analytical_expression  8.68  -0.0103  -2190
-Vm   14.58  1.84  4.14  -2.46  -3.20

-delta_h -0.396 kcal
-Vm  6.16  0  29.4  0  0.9 # ref. 2
Ca+2 + H2O = CaOH+ + H+
-log_k -12.78
Ca+2 + CO3-2 = CaCO3
-log_k 3.224
-delta_h 3.545 kcal
-analytic -1228.732 -0.299440 35512.75 485.818
-dw 4.46e-10 # complexes: calc'd with the Pikal formula
-Vm  -.2430  -8.3748  9.0417  -2.4328  -.0300 # supcrt
Ca+2 + CO3-2 + H+ = CaHCO3+
-log_k 11.435
-delta_h -0.871 kcal
-analytic 1317.0071 0.34546894 -39916.84 -517.70761 563713.9
-gamma 6.0 0
-dw 5.06e-10
-Vm  3.1911  .0104  5.7459  -2.7794  .3084 5.4 # supcrt
2Cl- + Ca+2  =  CaCl2   
        log_k           -0.6436
-delta_h -5.8325 kJ/mol # Calculated enthalpy of reaction CaCl2
# Enthalpy of formation: -211.06 kcal/mol
        -analytic 1.8178e+002 7.6910e-002 -3.1088e+003 -7.8760e+001 -4.8563e+001
#       -Range:  0-300

H+ + Cl- = HCl
-log_k  -0.5
-analytical_expression  0.334  -2.684e-3  1.015 # from Pitzer.dat, up to 15 M HCl, 0 - 50°C
-gamma  0  0.4256

END




Error information:

Code: [Select]

SOLUTION_RAW                 2304 Solution after simulation 1.
  -temp                      60
  -pressure                  100
  -total_h                   0.095811131926754
  -total_o                   0.047908152856899
  -cb                        -2.823231909479e-05
  -density                   5.3998185596486
  -totals
    C                        4.0364888717517e-06
    Ca                       0.043624548285507
    Cl                       0.087288301058543
    H2O                      0
  -pH                        2.4144543407868
  -pe                        1.6050376215329
  -mu                        14.291834424749
  -ah2o                      0.02504013734053
  -mass_water                0.00086320996138556
  -soln_vol                  0.0010861515257731
  -total_alkalinity          -4.9209075364978e-06
  -activities
    C(4)                     -14.557124780346
    Ca                       1.9142597622842
    Cl                       0.85511383656121
    E                        -1.6050376215329
    H(0)                     -11.369300642517
    O(0)                     -62.818347710313
  -gammas
  -species_map
    0 2.2368336837094e-06
    1 1.0276889800347e-12
    2 0.0037581979967257
    3 2.6578823195791e-14
    4 4.0759235929202
    5 3.4643761868298e-11
    6 36.855862360579
    7 0.00029962982747638
    8 9.9048503640988e-13
    9 8.1661248579579
    10 0.0048304248887352
    11 1.3617743923983e-13
    12 47.527553493865
    13 2.2064441141925e-11
    14 2.1056574872548e-07
    15 5.8398132073799e-09
    16 0
    17 1.2943077267615e-12
  -log_gamma_map
    0 1.4291834424749
    1 1.4291834424749
    2 0.94326107203345
    3 -1.0490566960506
    4 1.2366375606398
    5 1.4291834424749
    6 1.4291834424749
    7 -0.23909330015695
    8 1.8981769210125
    9 -0.12429834682639
    10 -0.16583584166974
    11 1.4291834424749
    12 0
    13 0.00018981769210125
    14 -0.26226417401264
    15 6.0826047311733
    16 1.4291834424749
    17 -0.37838544448794
SURFACE_RAW                  2304
  # SURFACE_MODIFY candidate identifiers #
  -type                      2
  -dl_type                   0
  -only_counter_ions         0
  -thickness                 1e-08
  -debye_lengths             0
  -DDL_viscosity             1
  -DDL_limit                 0.8
  # SURFACE_MODIFY candidates with new_def=true #
  -new_def                   0
  -sites_units               0
  -solution_equilibria       0
  -n_solution                -999
  # Surface workspace variables #
  -transport                 0
  -totals                   


[dlparkhurst]

PHREEQC works best with approximately 1 kg water, so it is best to scale your problem to larger water volumes.

However, I do not think the mass of water alone is the problem. The concentration in your solution is approximately 100 molal Cl. That is way past the  concentration range for any of the databases. So, somehow you have gone from 2 m CaCl2 to 50 m CaCl2. PHREEQC will fail with concentrations that large.

[Mohamadreza]

I did not set the concentrations of Cl and Ca to 100 and 50 m; could the non-convergence be due to this? I set them to 4 and 2 instead.
When the simulation does not converge, it reaches unrealistic values and fails. What could be the underlying issue causing this behavior? What's your opinion to solve this error? Thanks.

[dlparkhurst]

The values of SOLUTION_RAW should be those after transport, but before reaction. That is where I am getting the values to calculate concentration. Could your transport be incorrect?

Try conservative transport with no reactions to see if the CaCl2 concentrations are transported conservatively.

[Mohamadreza]

How do I change it to conservative transport with no reaction?

[dlparkhurst]

Only define define solutions--no surfaces, equilibrium_phases, etc. Make sure the Phreeqc output file for the solutions matches your expectations.


If the calculation is successful, then you can add back reactants to find the problem.

If the calculation fails, remove first RunCells from the time stepping, and if it still fails, SetConcentrations, GetConcentrations from the time stepping. That will leave you with only transport.

[Mohamadreza]

Thanks. How do I fixed mass or mole of water in geometry to maintain its amount in the geometry?
Do I need to explicitly inject water to maintain its amount in the geometry? Does PHREEQCRM (used for reactive transport) preserve the amount of water in the geometry?

[dlparkhurst]

PhreeqcRM uses relies completely on the concentrations given in SetConcentrations and the volume--as determined by RV, porosity, saturation, and sometimes density (when using mass fraction)--to arrive at the number of moles of each element in each cell. If there are no reactants (equilibirum_phases, surfaces, etc), then Phreeqc will preserve the number of moles of each element in solution, including H and O (and charge balance).

[Mohamadreza]

Thanks. Well, if we want to investigate the dissolution of carbon dioxide in a 2 molal CaCl₂ solution, the amount of water (in moles) will be consumed. And if a cell runs out of water, PHREEQCRM will encounter the same problem that I experienced.

[dlparkhurst]

I think it might depend on your database. Phreeqc.dat and pitzer.dat define CO2(aq), not H2CO3(aq); so, dissolution of CO2 alone would not consume water. Activity of water will decrease as CO2 dissolves, basically as a mole fraction of the system. Phreeqc.dat uses 1 - 0.017*sum(m(i)) for activity, so there is a specific sum of molalities that will result in activity of water going to zero. pitzer.dat uses a formulation for activity of water based on the osmotic coefficient.

Depending on your system, water could be lost to a gas phase, to minerals, or perhaps to other aqueous species

Solubility of CO2 has been fit fairly carefully, although perhaps not as much in CaCl2 solutions. However, I think you should be able to model CO2 solubility in 2 molal CaCl2 within limits of say a few hundred atmospheres and 200 to 300 C. I think you said PHREEQC calculations were successful.

[Mohamadreza]

I simulate CO₂ solubility in CaCl₂ or NaCl solutions with concentrations ranging from 2 to 4 molal, using a custom PHREEQC input and the standard PHREEQC database. The simulation runs without any errors when using the PHREEQC software directly. Additionally, I have tested the same setup in the reactive transport package with a dilute solution (0.01 molal), and it also runs successfully without any issues.

If the geometry is defined at the microscale and each cell has a very small volume due to low porosity and saturation, does PHREEQCRM still assume 1 kg of water per cell, or does it calculate the amount of water based on the actual cell size?

[dlparkhurst]

"PhreeqcRM relies completely on the concentrations given in SetConcentrations and the volume--as determined by RV, porosity, saturation, and sometimes density (when using mass fraction)--to arrive at the number of moles of each element in each cell. "

The representative volume is 1 L by default. Extremely low porosity could cause problems because the water volume is small. PHREEQC runs best when the mass of water is within a couple of orders of magnitude of 1. However, if the simulation runs in Phreeqc, it should run in PhreeqcRM if you replicate the conditions.

Does the transport work correctly when you skip the SetConcentrations, RunCells, and GetConcentrations in the time stepping loop?

[Mohamadreza]

Hi, I don't know how to turn off runcell and etc, so I must ask the developer of the package.

[dlparkhurst]

But it fails if you define just solutions with no other reactants?

[Mohamadreza]

Yes, the simulation failed and shows an error related to RunCell and the database. I think this solver does not include transport, as it investigates mass transfer at the interface. However, this should be confirmed with the developer to resolve the issue. Thank you and I appreciate your guidance.

[dlparkhurst]

The only other suggestion I have is (1) to do a MIX; n 1.0 in the input file for every SOLUTION n that you define. That will insure that there is no problem with the solution definitions when they react to redox equilibrium, and (2) to take a tiny or zero time step with your transport code. That should be almost equivalent to the MIX calculations, assuming there are no other reactants.

[Mohamadreza]

Have error again.

Code: [Select]

SOLUTION_RAW                 1614 Solution after simulation 1.
  -temp                      60
  -pressure                  100
  -total_h                   0.67130554146401
  -total_o                   0.85983634067137
  -cb                        0.073206320882858
  -density                   2.3142631433938
  -totals
    C                        0.24151535276651
    Ca                       0.13715378795726
    Cl                       0.11879552615305
    H2O                      0
  -pH                        13.909184647295
  -pe                        -13.519212114321
  -mu                        41.542710410569
  -ah2o                      4.5201078730186e-21
  -mass_water                0.015492546678764
  -soln_vol                  0.011645386226071
  -total_alkalinity          0.082006749529583
  -activities
    C(4)                     -5.5182198475323
    Ca                       6.8661503126831
    Cl                       1.603403724178
    E                        13.519212114321
    H(0)                     -4.0501771737544
    O(0)                     -114.9434380835
  -gammas
  -species_map
    0 8.4846372557809
    1 0.075446610889288
    2 2.0004782687933e-05
    3 8.2125646730294
    4 3.4945699251791
    5 0.052814963939728
    6 1.2627750463038e-19
    7 10.187970628791
    8 9.4895701644842e-15
    9 2.8204317339416e-09
    10 28.787544910685
    11 3.4140270775793e-30
    12 3.1779245935487e-10
    13 1.5961701272666e-31
    14 0
    15 4.2138593735943e-20
  -log_gamma_map
    0 4.1542710410569
    1 4.1542710410569
    2 -1.107580343223
    3 5.6658860609835
    4 4.1542710410569
    5 -0.25119351275409
    6 6.2954139521402
    7 0.30962053719855
    8 -0.17156818687514
    9 4.1542710410569
    10 0
    11 0.00062954139521402
    12 -0.27689508580575
    13 17.680577550738
    14 4.1542710410569
    15 -0.4096120145322
SURFACE_RAW                  1614
  # SURFACE_MODIFY candidate identifiers #
  -type                      2
  -dl_type                   0
  -only_counter_ions         0
  -thickness                 1e-08
  -debye_lengths             0
  -DDL_viscosity             1
  -DDL_limit                 0.8
  # SURFACE_MODIFY candidates with new_def=true #
  -new_def                   0
  -sites_units               0
  -solution_equilibria       0
  -n_solution                -999
  # Surface workspace variables #
  -transport                 0
  -totals                   

[dlparkhurst]

Sorry, I can't help you with your simulator. I can help debug if you reproduce a problem in Phreeqc, PHAST, or a test case in C++, Fortran, C, or Python.

You could give more context for your error. Did you run a short time step? Does it fail if you run no time steps? Are you using only SOLUTIONs?

[Mohamadreza]

No, I don't encounter any errors with short time steps. Yes, I am using only one solution. In the multiphase flow (VOF), the volume fraction of water in the porous medium is initially 1. When CO₂ is injected into the system, the volume fraction of CO₂ increases while the water fraction decreases. Because of this, I get an error related to the activity of water.

I also tested injecting CO₂ and water together into the system to fixed the fraction of water to 1, and in that case, I did not get any errors. However, this approach is not correct because I cannot observe the proper CO₂ phase distribution.

Why does PHREEQCRM require the moles of water to remain constant? This seems unrealistic.

[dlparkhurst]

If the simulation runs with small time steps, perhaps your transport solution is unstable. For rectangular grids and centered-in-time, centered-in-space, you can get oscillations if your time step is too large. The stability relations are deltaX/alpha < 2, and alpha*V*deltaT / deltaX^2 < 1. Upstream-in-space and backward-in-time are unconditionally stable, so if you can use fully implicit weighting you may avoid your problems. From what I remember, you generated extremely large concentrations, which sounds to me like oscillations in the transport solution. This is consistent with your simulator working with smaller time steps.

The moles of water are not required to be constant. The transport values are concentration. The conversion to PHREEQC units involves the volume of water in the cell, which is RV * porosity * saturation. If the saturation is less, then the volume of solution, and volume of water, will be less.

[Mohamadreza]

Thanks for your explanations.
I have a question: If a cell is completely occupied by CO2 gas and contains no water, what happens? Does phreeqcRM produce an error in this case?

[dlparkhurst]

If the saturation is zero, I think no calculation is made for the cell.

[Mohamadreza]

Thanks. As I mentioned, my geometry is at the pore scale (micrometer).
The water volume in each cell is approximately 1e-15 m3, based on porosity, saturation, and representative volume.
How can the high salt concentration  (1 or 2 mol/kgw) be scaled to this unit? Or scale phreeqc calculation to this unit? Maybe the error I get is due to the fact that when the number of moles of salt is high and the water volume is very small, it causes the error.

[dlparkhurst]

There is no reason the representative volume should be a tiny number. phreeqcRM expects concentrations, so you can use a larger representative  volume to get a larger volume of water. You just have to scale the other reactants appropriately.

Sorry, but I am done with this thread. Unless you can demonstrate your problem with PHREEQC or Phast (or possibly Python), I cannot help help you.

[Mohamadreza]

Thank you very much for your help.