Simulating CO₂-Dissolved Solution Interaction with Rock and Cement Min

[ahmalisha]

Hello everyone,

I am working on a PHREEQC simulation and need help fine-tuning the process. My goal is to:

Simulate the dissolution of CO₂ in water.
Model the interaction of the resulting CO₂-dissolved solution with rock minerals (e.g., calcite, quartz, dolomite).
Use the resulting solution to interact with cement minerals (e.g., Portlandite, C-S-H, Ettringite) to study chemical changes and porosity evolution.
Below is the code I have developed so far: I have got this ERROR: Numerical method failed on all combinations of convergence parameters, cell/soln/mix 1. Could you please help me   

[dlparkhurst]

I don't have your database. Here is a modified version of your script with phases that I found in other databases. Runs for me.

If you have problems, please simplify your script as much as possible.

Code: [Select]

RATES
    Calcite
-start
 10 REM PARM(1) = SA (BET surface area/kg water)
 20 REM PARM(2) = k25a (rate constant at 25?C in mol/m?s acid mechanism)
 30 REM PARM(3) = k25n (rate constant at 25?C in mol/m?s neutral mechanism)
 40 REM PARM(4) = k25b (rate constant at 25?C in mol/m?s base mechanism)
 50 REM PARM(5) = Eaa (activation energy in J/mol acid mechanism)
 60 REM PARM(6) = Ean (activation energy in J/mol neutral mechanism)
 70 REM PARM(7) = Eab (activation energy in J/mol base mechanism)
 80 REM PARM(8) = na (order of H+ catalysis acid mechanism)
 90 REM PARM(9) = nb (order of H+ catalysis base mechanism)
100 REM PARM(10) = upper pH limit acid mechanism
110 REM PARM(11) = upper pH limit neutral mechanism
120 IF (SI("Calcite") < 0) THEN GOTO 130 ELSE GOTO 310
130 REM Dissolution Mechanism
140 IF (M <= 0) THEN GOTO 310
150 t = 1
160 IF (M0 > 0) THEN t = M / M0
170 IF ((-LA("H+")) < PARM(10)) THEN GOTO 230
180 IF ((-LA("H+")) > PARM(11)) THEN GOTO 270
190 REM neutral mechanism
200 k1 = PARM(3)*EXP((PARM(6) / -8.314) * ((1 / TK) - (1 / 298.15)))
210 rate = k1 * PARM(1) * t * (10^(SI("Calcite")) - 1)
220 GOTO 300
230 REM acid mechanism
240 k1 = PARM(2)*EXP((PARM(5) / -8.314) * ((1 / TK) - (1 / 298.15)))
250 rate = k1 * PARM(1) * t * ((ACT("H+"))^(PARM(8))) * (10^(SI("Calcite")) - 1)
260 GOTO 300
270 REM base mechanism
280 k1 = PARM(4)*EXP((PARM(7) / -8.314) * ((1 / TK) - (1 / 298.15)))
290 rate = k1 * PARM(1) * t * ((ACT("H+"))^(PARM(9))) * (10^(SI("Calcite")) - 1)
300 moles = -rate * TIME
310 SAVE moles
-end

    Dolomite
-start
 10 REM Similar setup to Calcite but with different parameters
120 IF (SI("Dolomite") < 0) THEN GOTO 130 ELSE GOTO 310
130 REM Dissolution Mechanism
140 IF (M <= 0) THEN GOTO 310
150 t = 1
160 IF (M0 > 0) THEN t = M / M0
170 IF ((-LA("H+")) < PARM(10)) THEN GOTO 230
180 IF ((-LA("H+")) > PARM(11)) THEN GOTO 270
190 REM neutral mechanism
200 k1 = PARM(3)*EXP((PARM(6) / -8.314) * ((1 / TK) - (1 / 298.15)))
210 rate = k1 * PARM(1) * t * (10^(SI("Dolomite")) - 1)
220 GOTO 300
230 REM acid mechanism
240 k1 = PARM(2)*EXP((PARM(5) / -8.314) * ((1 / TK) - (1 / 298.15)))
250 rate = k1 * PARM(1) * t * ((ACT("H+"))^(PARM(8))) * (10^(SI("Dolomite")) - 1)
260 GOTO 300
270 REM base mechanism
280 k1 = PARM(4)*EXP((PARM(7) / -8.314) * ((1 / TK) - (1 / 298.15)))
290 rate = k1 * PARM(1) * t * ((ACT("H+"))^(PARM(9))) * (10^(SI("Dolomite")) - 1)
300 moles = -rate * TIME
310 SAVE moles
-end

    Quartz
-start
 10 REM Similar setup to Calcite but with different parameters
120 IF (SI("Quartz") < 0) THEN GOTO 130 ELSE GOTO 310
130 REM Dissolution Mechanism
140 IF (M <= 0) THEN GOTO 310
150 t = 1
160 IF (M0 > 0) THEN t = M / M0
170 IF ((-LA("H+")) < PARM(10)) THEN GOTO 230
180 IF ((-LA("H+")) > PARM(11)) THEN GOTO 270
190 REM neutral mechanism
200 k1 = PARM(3)*EXP((PARM(6) / -8.314) * ((1 / TK) - (1 / 298.15)))
210 rate = k1 * PARM(1) * t * (10^(SI("Quartz")) - 1)
220 GOTO 300
230 REM acid mechanism
240 k1 = PARM(2)*EXP((PARM(5) / -8.314) * ((1 / TK) - (1 / 298.15)))
250 rate = k1 * PARM(1) * t * ((ACT("H+"))^(PARM(8))) * (10^(SI("Quartz")) - 1)
260 GOTO 300
270 REM base mechanism
280 k1 = PARM(4)*EXP((PARM(7) / -8.314) * ((1 / TK) - (1 / 298.15)))
290 rate = k1 * PARM(1) * t * ((ACT("H+"))^(PARM(9))) * (10^(SI("Quartz")) - 1)
300 moles = -rate * TIME
310 SAVE moles
-end

    Kaolinite
-start
 10 REM Similar setup to Calcite but with different parameters
120 IF (SI("Kaolinite") < 0) THEN GOTO 130 ELSE GOTO 310
130 REM Dissolution Mechanism
140 IF (M <= 0) THEN GOTO 310
150 t = 1
160 IF (M0 > 0) THEN t = M / M0
170 IF ((-LA("H+")) < PARM(10)) THEN GOTO 230
180 IF ((-LA("H+")) > PARM(11)) THEN GOTO 270
190 REM neutral mechanism
200 k1 = PARM(3)*EXP((PARM(6) / -8.314) * ((1 / TK) - (1 / 298.15)))
210 rate = k1 * PARM(1) * t * (10^(SI("Kaolinite")) - 1)
220 GOTO 300
230 REM acid mechanism
240 k1 = PARM(2)*EXP((PARM(5) / -8.314) * ((1 / TK) - (1 / 298.15)))
250 rate = k1 * PARM(1) * t * ((ACT("H+"))^(PARM(8))) * (10^(SI("Kaolinite")) - 1)
260 GOTO 300
270 REM base mechanism
280 k1 = PARM(4)*EXP((PARM(7) / -8.314) * ((1 / TK) - (1 / 298.15)))
290 rate = k1 * PARM(1) * t * ((ACT("H+"))^(PARM(9))) * (10^(SI("Kaolinite")) - 1)
300 moles = -rate * TIME
310 SAVE moles
-end

Aragonite
-start
 40 SR_mineral=SR("Aragonite")
 41 if (M<0) then goto 200
 42 if (M=0 and SR_mineral<1) then goto 200
 43 SA=PARM(1)
 50 if (SA<=0) then SA=1
 60 R=8.31451
 70 E=31379.792
 80 Rate=5.91e-4*EXP(-E/R/TK)
 90 Rate=Rate*(1-SR_mineral)*SA
100 moles= rate*Time
200 save moles
-end

    Vaterite
-start
 40 SR_mineral=SR("Vaterite")
 41 if (M<0) then goto 200
 42 if (M=0 and SR_mineral<1) then goto 200
 43 SA=PARM(1)
 50 if (SA<=0) then SA=1
 90 Rate=1.491e-5*(1-SR_mineral)*SA
100 moles= rate*Time
200 save moles
-end

    C2S
-start
 40 SR_mineral=SR("C2S")
 41 if (M<0) then goto 200
 42 if (M=0 and SR_mineral<1) then goto 200
 43 SA=PARM(1)
 50 if (SA<=0) then SA=1
 90 Rate=1.98e-2*(1-SR_mineral)*SA
100 moles= rate*Time
200 save moles
-end

    C3S
-start
 40 SR_mineral=SR("C3S")
 41 if (M<0) then goto 200
 42 if (M=0 and SR_mineral<1) then goto 200
 43 SA=PARM(1)
 50 if (SA<=0) then SA=1
 90 Rate=1.98e-5*(1-SR_mineral)*SA
100 moles= rate*Time
200 save moles
-end

    ECSH1-TobCa
-start
 10 rem acid solution parameters
 11 a1=6.32E-04
 12 E1=23000
 13 n1=0.275
 20 rem neutral solution parameters
 21 a2=1.71E-14
 22 E2=23000
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("CSH(0.8)")
 40 SR_mineral=SR("ECSH1-TobCa")
 41 if (M<0) then goto 200
 42 if (M=0 and SR_mineral<1) then goto 200
 50 if (SA<=0) then SA=1
 60 R=8.31451
 75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1
 80 Rate2=a2*EXP(-E2/R/TK)               
 90 Rate=(Rate1+Rate2)*(1-SR_mineral)*SA
100 moles= rate*Time
200 save moles
-end
END
PHASES
C2S
(CaO)2SiO2 + 4H+ + 2H2O = 2Ca+2 + 2H2O + H4SiO4
-Vm 51.79
-analytical_expression -4.75838 0 12467.437992 0.569296 0 0 0
-log_K 38.567691

C3S
(CaO)3SiO2 + 6H+ + 2H2O = 3Ca+2 + 3H2O + H4SiO4
-Vm 73.18
-analytical_expression -6.720801 0 23294.175088 0.748984 0 0 0
-log_K 73.405906

ECSH1-TobCa
((Ca(OH)2)0.8333SiO2H2O)1 + 1.6666H+ +2H2O = 0.8333Ca+2 + 2.6666H2O + H4SiO4
-Vm 68
-analytical_expression -13.776918 0 3023.19863 5.923868 0 0 0
-log_K 11.019995

Portlandite
Ca(OH)2 + 2H+ = Ca+2 + 2H2O
-Vm 33.06
-analytical_expression -11.299363 0 7301.394065 3.883957 0 0 0
-log_K 22.799937

Vaterite
CaCO3 = Ca+2 + CO3-2
log_k -7.913 # PB82
-analytic -172.1295 -0.077993 3074.688 71.595 # PB82
-Vm 37.628 # Webmineral
END
SOLUTION 1
    temp      140
    pH        5.5  # Slightly acidic to drive carbonate dissolution
    pe        4
    redox     pe
    units     mg/l
    C(4)      10000  # Higher initial CO2 to prompt dissolution
    Ca        54.6
    Cl        18927
    Fe        0.02
    K         111
    Mg        12.1
    Na        8504
    S(6)      1300
    Sr        0.05
    -water    1000 # kg

GAS_PHASE 1
    -fixed_volume               # Switch to a fixed volume to allow CO2 to dissolve
    -volume 1                   # Volume in liters
    CO2(g)    81.5292           # CO2 partial pressure in atm
    H2S(g)    0.126             # H2S partial pressure in atm
    -temperature 140            # Temperature in ?C

EQUILIBRIUM_PHASES 1
    Calcite                  0 78200
    Dolomite                 0 60900
    Quartz                   0 648400
    Kaolinite                0 84100

SAVE SOLUTION 2
END

USE SOLUTION 2

GAS_PHASE 1
    -fixed_pressure
    -pressure 81.6552
    -volume 1
    -temperature 140
    CO2(g)    810.5292
    H2S(g)    0.126

KINETICS 1
Calcite
    -formula  CaCO3  1
    -m        78.2
    -m0       78.2
    -parms    67.1 0.501 1.55e-06 0.0003311 14400 23500 35400 1 1 6 8
    -tol      1e-08
Dolomite
    -formula  CaMg(CO3)2  1
    -m        60.9
    -m0       60.9
    -parms    63.6 6.457e-06 2.9512e-08 7.76e-06 36100 52200 34800 1 1 6 8
    -tol      1e-08
Quartz
    -formula  SiO2  1
    -m        648.4
    -m0       648.4
    -parms    68.6 0 1.023e-14 0 0 87700 0 0 0 6 8
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        84.1
    -m0       84.1
    -parms    1160 4.9e-12 6.61e-14 8.9e-18 65900 22200 17900 0.78 -0.47 6 8
    -tol      1e-08
SAVE SOLUTION 2
END
USE SOLUTION 2
EQUILIBRIUM_PHASES 2
    CO2(g)    2.45 10
    C2S       0 0.2
    C3S       0 0.3
    ECSH1-TobCa 0 0.2
    Calcite   0 0.1
    Quartz    0 0.2
    Aragonite 0 0.1
    Vaterite  0 0.1
    Portlandite 0 1.15

KINETICS 2
Calcite
    -formula  CaCO3  1
    -m        78.2
    -m0       78.2
    -parms    67.1 0.501 1.55e-06 0.0003311 14400 23500 35400 1 1 6 8
    -tol      1e-08

Quartz
    -formula  SiO2  1
    -m        648.4
    -m0       648.4
    -parms    68.6 0 1.023e-14 0 0 87700 0 0 0 6 8
    -tol      1e-08

C2S
    -formula  (CaO)2SiO2  1
    -m        1
    -m0       20
    -parms    10
    -tol      1e-08
C3S
    -formula  (CaO)3SiO2  1
    -m        1
    -m0       20
    -parms    10
    -tol      1e-08
ECSH1-TobCa
    #-formula  ((Ca(OH)2)0.8333SiO2H2O)1  1
    -m        1
    -m0       20
    -parms    10
    -tol      1e-08
Aragonite
    -formula  Aragonite  1
    -m        0
    -m0       0
    -parms    91.5 1
    -tol      1e-08
Vaterite
    -formula  CaCO3  1
    -m        0
    -m0       0
    -parms    91.5 1
    -tol      1e-08

-steps       1e0 1e1 1e2 1e3 1e4 1e5 1e6 1e7 6.3072e7 # Seconds up to exactly 2 years
-step_divide 1
-runge_kutta 3
-bad_step_max 5000
-cvode true
-cvode_steps 1000
-cvode_order 5
#INCREMENTAL_REACTIONS True

USER_GRAPH 3
    -headings               time Volume_Change_Calcite Volume_change_Dolomite Volume_change_Quartz Volume_change_Kaolinite
    -axis_titles            "Time (Days)" "Phase Volume Changes (cm?/100g)"
    -chart_title            "Phase Volume Changes for Each Mineral Over Time"
    -axis_scale x_axis      auto auto auto auto log
-start
10 GRAPH_X TOTAL_TIME / 86400
20 v_calcite = KIN_DELTA("Calcite") * 36.89
30 v_dolomite = KIN_DELTA("Dolomite") * 64.5
40 v_quartz = KIN_DELTA("Quartz") * 22.69
50 v_kaolinite = KIN_DELTA("Kaolinite") * 99.35
60 GRAPH_Y v_calcite, v_dolomite, v_quartz, v_kaolinite
-end

USER_GRAPH 4
    -headings               "Time (Days)" "Porosity (%)"
    -axis_titles            "Time (Days)" "Porosity (%)"
    -chart_title            "Porosity Change Over Time"
-start
10 GRAPH_X TOTAL_TIME / 86400  # Time in days
20 v_calcite = KIN_DELTA("Calcite") * 36.89
30 v_dolomite = KIN_DELTA("Dolomite") * 64.5
40 v_quartz = KIN_DELTA("Quartz") * 22.69
50 v_kaolinite = KIN_DELTA("Kaolinite") * 99.35

60 v_total = v_calcite + v_dolomite + v_quartz + v_kaolinite
70 initial_volume = 1000  # Assumed initial rock volume in cm?
80 initial_porosity = 10  # Initial porosity in percentage
90 porosity = initial_porosity + (v_total / initial_volume) * 100  # Adjusted porosity
100 GRAPH_Y porosity
-end

-active true
USER_GRAPH 5
    -headings               "Time (Days)" "SI_Calcite" "SI_Dolomite" "SI_Quartz" "SI_Kaolinite" "SI_Pyrite" "SI_Illite" "SI_Gypsum"
    -axis_titles            "Time (Days)" "Saturation Index"
    -chart_title            "Saturation Index of Minerals Over Time"
    -axis_scale x_axis      auto auto auto auto log
-start
10 GRAPH_X TOTAL_TIME / 86400  # Time in days
20 SI_Calcite = SI("Calcite")
30 SI_Dolomite = SI("Dolomite")
40 SI_Quartz = SI("Quartz")
50 SI_Kaolinite = SI("Kaolinite")
60 GRAPH_Y SI_Calcite, SI_Dolomite, SI_Quartz, SI_Kaolinite
-end
-active true

[ahmalisha]

Thank you so much for your help. Can you please share the database that you modify?

[dlparkhurst]

I used phreeqc.dat (version 3.8.5, but any version should work). All the additions are in the input file.

[dlparkhurst]

You should not put minerals in both EQUILIBRIUM_PHASES and KINETICS. EQUILIBRIUM_PHASES will cause the minerals to immediately react to equilibrium, so what is the point of kinetics? Choose one or the other.

Calcite is stable relative to vaterite and aragonite in most cases. Putting all three in EQUILIBRIUM_PHASES will cause vaterite and aragonite to dissolve completely to form calcite. Please look at your output file.