No element or master species given for concentration input

[kennedyantwi1]

Dear Dr. Parkhurst,
I've been trying to run a CO2-rock experiment in PHREEQC and there's an Error message as follows:
                                          No element or master species given for concentration input
Please note, there is no pH value given for the experiment and I may want to get outputs for the pH with Time as well. The code is included for you suggestions and guidance.

Code: [Select]

[SOLUTION_SPECIES
CO3-2 + 2 H+ = CO2 + H2O
     log_k 16.681
     delta_h -5.738 kcal
    -analytical_expression 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
END

PHASES
CO2(g)
CO2 = CO2
     log_k -1.468
     delta_h -4.776 kcal
     -analytical_expression 10.5624 -2.3547e-2 -3972.8 0 5.8746e5 1.9194e-5
     -T_c 304.2 # critical T, K
     -P_c 72.86 # critical P, atm
     -Omega 0.225 # acentric factor
Ankerite
CaMg0.3Fe0.7(CO3)2 + 2H+ = Ca++ + 0.3Mg++ + 0.7Fe++ + 2HCO3-
     log_k 1.54
     delta_h 0
     -analytical_expression -1.8649e+03 -2.9583e-01 1.0468e+05 6.7554e+02 -6.0514e+06
END


SOLUTION 1 #Brine with auto-calculated pH
temp      44
units     mol/kgw
    K         0.08655
    Ca        0.1166     
    Mg        0.01778
    Na        1.305
    Cl        1.2397
    C(4)      0.002
    S(6)      0.003
    water     1.0 #1.0 kg
    charge   #PHREEQC will adjust pH to maintain charge balance
    pressure 148.038 atm #Experimental pressure


EQUILIBRIUM_PHASES
    CO2(g)        1.398  10.0  # 25 bar PCO2
    Calcite       0.0  0.017742 
    Dolomite      0.0  0.0.0321001
    K-feldspar    0.0  0.017014
    Albite        0.0  0.040634
    Ankerite      0.0  0.022945
    Quartz        0.0  0.472921
    Illite        0.0  0.0082214
END


EQUILIBRIUM_PHASES 1
CO2(g)  1.398  10.0  # 25 bar PCO2 Partial pressure of carbon dioxide at 50 atm, and injected mol

RATES
Calcite
   -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_calc = SI("Calcite")
30 if (M <= 0 and si_calc < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_calc > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-0.30)*EXP(-14.40e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(1.000)
70 k_neut = 10^(-5.81)*EXP(-23.50e+03/8.314*(1.0/TK-1.0/298.15))
80 k_carb = 10^(-3.48)*EXP(-35.40e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(1.000)
90 k_rateconst = k_acid + k_neut + k_carb
100 r = k_rateconst * SA * (1-(10^si_calc))
190 moles = r * TIME
200 SAVE moles
   -end

Dolomite
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_dolo = SI("Dolomite")
30 if (M <= 0 and si_dolo < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_dolo > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-3.76)*EXP(-56.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)
70 k_neut = 10^(-8.60)*EXP(-95.30e+03/8.314*(1.0/TK-1.0/298.15))
80 k_carb = 10^(-5.37)*EXP(-45.70e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(0.500)
90 k_rateconst = k_acid + k_neut + k_carb
100 r = k_rateconst * SA * (1-(10^si_dolo))
190 moles = r * TIME
200 SAVE moles
    -end

K-Feldspar
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_kfeld = SI("K-Feldspar")
30 if (M <= 0 and si_kfeld < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_kfeld > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-10.06)*EXP(-51.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)
70 k_neut = 10^(-12.41)*EXP(-38.00e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-21.20)*EXP(-94.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.823)
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_kfeld))
190 moles = r * TIME
200 SAVE moles
    -end

Albite
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_alb = SI("Albite")
30 if (M <= 0 and si_alb < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_alb > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-10.16)*EXP(-65.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("H+")^(0.457))
70 k_neut = 10^(-12.56)*EXP(-69.80e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-15.60)*EXP(-71.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("OH-")^(-0.572))
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA *(1-(10^si_alb))
190 moles = r * TIME
200 SAVE moles
    -end
Quartz
  -start
1  REM  Specific rate k from Rimstidt and Barnes, 1980, GCA 44,1683
2  REM  k = 10^-13.7 mol/m2/s (25 C), Ea = 90 kJ/mol
3  REM  sp. rate * parm(2) due to salts (Dove and Rimstidt, MSA Rev. 29, 259)
4  REM  PARM(1) = Specific area of Quartz, m^2/mol Quartz
5  REM  PARM(2) = salt correction: (1 + 1.5 * c_Na (mM)), < 35
10 dif_temp = 1/TK - 1/298
20 pk_w = 13.7 + 4700.4 * dif_temp
40 moles = PARM(1) * M0 * PARM(2) * (M/M0)^0.67 * 10^-pk_w * (1 -  SR("Quartz"))
# Integrate...
50 SAVE moles * TIME
  -end

KINETICS 1
Calcite
    -formula CaCO3 1
    -m 10
    -m0 10
    -parms 2.2661
    -tol 1e-08
Dolomite
     -formula CaMg(CO3)2 1
     -m 10
     -m0 10
     -parms 3.8685
     -tol 1e-08
K-Feldspar
    -formula KAlSi3O8 1
    -m 17.11
    -m0 17.11
    -parms 6.523
    -tol 1e-08
Albite
    -formula NaAlSi3O8 1
    -m 0.5650
    -m0 0.5650
    -parms 6.02
    -tol 1e-08

Quartz
    -formula SiO2 1
    -m 207.05
    -m0 207.05
    -parms 1.36
    -tol 1e-08


    -steps    0 86400 172800 345600 691200 2073600 4147200 6220800 8208000 10368000 13392000  # Converted days to seconds


USER_PUNCH
    -headings Time pH SI_Calcite SI_Dolomite SI_K_feldspar SI_Albite SI_Ankerite SI_Quartz SI_Illite new_phi
    -start
    10  REM Initial conditions
    20  phi_0 = 0.1244   # Initial porosity
    40  rho_CO2 = 65.7        # CO2 density at 50?C, 25 bar (kg/m^3)
   
    50  REM Calculate mineral volumes directly
    60  vol_calc = TOT("Calcite")*36.9 # molar volume of calcite = molar mass/ density = 100.09g/mol/ 2.71g/cm3 = 36.9cm3/mol
    80  vol_dol = TOT("Dolomite")*64.365
    90  vol_K_fel = KIN("K-feldspar")*108.87
    100 vol_alb = KIN("Albite")*100.389
    110 vol_ank = KIN("Ankerite")*65.56
    113 vol_qtz = KIN("Quartz")*22.67
    114 vol_ill = KIN("Illite")*141.47

    120 REM Calculate total mineral volume in cm?
    130 vol_total = vol_calc + vol_dol + vol_fel + vol_alb + vol_ank + vol_qtz + vol_ill
   
    140 REM Set representative volume to 1000 cm? #1.0 L
    150 RV = 1000
   
    160 REM Update porosity (ensure it doesn't go negative)
    170 new_phi = phi_0 - vol_total/RV
   
    240 REM Print debugging info to verify calculations
    250 PUT(vol_total, 1)
   
    260 REM Output results with formatting to ensure visibility
    270 punch TOTAL_TIME
    280 punch -LA("H+")
    290 punch SI("Calcite")
    300 punch SI("Dolomite")
    310 punch SI("K-feldspar")
    320 punch SI("Albite")
    330 punch SI("Ankerite")
    340 punch SI("Quartz")
    350 punch SI("Illite")
    360 punch new_phi
    -end

SELECTED_OUTPUT
    -file               EXP_CO2_correction 12.04.2025_003.xls
    -reset              false
    -time               true
    -pH                 true
    -alkalinity         true
    -ionic_strength     true
    -totals             K Na Ca Mg Cl S(6) C(4)
    -molalities  Ca+2 Cl- K+ Mg+2 Na+ SO4-2 CO3-2 HCO3-
    -saturation_indices Calcite Dolomite K-feldspar Albite Ankerite Quartz Illite
    -equilibrium_phases Calcite Dolomite K-feldspar Albite Ankerite Quartz Illite

    -kinetic_reactants Calcite Dolomite K-feldspar Albite Quartz
    -kinetics
    -gases CO2(g) H2O(g)
    -water
    -charge_balance true

PRINT
    -saturation_indices Calcite Dolomite K-feldspar Albite Ankerite Quartz Illite

REACTION 1
CO2(g)
END]

[dlparkhurst]

Code: [Select]

pH 7 charge   #PHREEQC will adjust pH to maintain charge balance


[kennedyantwi1]

Thanks Dr. Parkhurst,
I've been able to run the model after your corrections.
My results are however not showing in the Excel output.
I'd appreciate your suggestions once more, please.
 

[dlparkhurst]

The following runs. I added "USE solution 1", and parm(2) for quartz.

Code: [Select]

[SOLUTION_SPECIES
CO3-2 + 2 H+ = CO2 + H2O
     log_k 16.681
     delta_h -5.738 kcal
    -analytical_expression 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
END

PHASES
CO2(g)
CO2 = CO2
     log_k -1.468
     delta_h -4.776 kcal
     -analytical_expression 10.5624 -2.3547e-2 -3972.8 0 5.8746e5 1.9194e-5
     -T_c 304.2 # critical T, K
     -P_c 72.86 # critical P, atm
     -Omega 0.225 # acentric factor
Ankerite
CaMg0.3Fe0.7(CO3)2 + 2H+ = Ca++ + 0.3Mg++ + 0.7Fe++ + 2HCO3-
     log_k 1.54
     delta_h 0
     -analytical_expression -1.8649e+03 -2.9583e-01 1.0468e+05 6.7554e+02 -6.0514e+06
END


SOLUTION 1 #Brine with auto-calculated pH
temp      44
units     mol/kgw
    K         0.08655
    Ca        0.1166     
    Mg        0.01778
    Na        1.305
    Cl        1.2397
    C(4)      0.002
    S(6)      0.003
    water     1.0 #1.0 kg
    pH 7 charge   #PHREEQC will adjust pH to maintain charge balance
    pressure 148.038 atm #Experimental pressure


EQUILIBRIUM_PHASES
    CO2(g)        1.398  10.0  # 25 bar PCO2
    Calcite       0.0  0.017742
    Dolomite      0.0  0.0.0321001
    K-feldspar    0.0  0.017014
    Albite        0.0  0.040634
    Ankerite      0.0  0.022945
    Quartz        0.0  0.472921
    Illite        0.0  0.0082214
END
USE solution 1

EQUILIBRIUM_PHASES 1
CO2(g)  1.398  10.0  # 25 bar PCO2 Partial pressure of carbon dioxide at 50 atm, and injected mol

RATES
Calcite
   -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_calc = SI("Calcite")
30 if (M <= 0 and si_calc < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_calc > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-0.30)*EXP(-14.40e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(1.000)
70 k_neut = 10^(-5.81)*EXP(-23.50e+03/8.314*(1.0/TK-1.0/298.15))
80 k_carb = 10^(-3.48)*EXP(-35.40e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(1.000)
90 k_rateconst = k_acid + k_neut + k_carb
100 r = k_rateconst * SA * (1-(10^si_calc))
190 moles = r * TIME
200 SAVE moles
   -end

Dolomite
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_dolo = SI("Dolomite")
30 if (M <= 0 and si_dolo < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_dolo > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-3.76)*EXP(-56.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)
70 k_neut = 10^(-8.60)*EXP(-95.30e+03/8.314*(1.0/TK-1.0/298.15))
80 k_carb = 10^(-5.37)*EXP(-45.70e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(0.500)
90 k_rateconst = k_acid + k_neut + k_carb
100 r = k_rateconst * SA * (1-(10^si_dolo))
190 moles = r * TIME
200 SAVE moles
    -end

K-Feldspar
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_kfeld = SI("K-Feldspar")
30 if (M <= 0 and si_kfeld < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_kfeld > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-10.06)*EXP(-51.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)
70 k_neut = 10^(-12.41)*EXP(-38.00e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-21.20)*EXP(-94.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.823)
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_kfeld))
190 moles = r * TIME
200 SAVE moles
    -end

Albite
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_alb = SI("Albite")
30 if (M <= 0 and si_alb < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_alb > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-10.16)*EXP(-65.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("H+")^(0.457))
70 k_neut = 10^(-12.56)*EXP(-69.80e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-15.60)*EXP(-71.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("OH-")^(-0.572))
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA *(1-(10^si_alb))
190 moles = r * TIME
200 SAVE moles
    -end
Quartz
  -start
1  REM  Specific rate k from Rimstidt and Barnes, 1980, GCA 44,1683
2  REM  k = 10^-13.7 mol/m2/s (25 C), Ea = 90 kJ/mol
3  REM  sp. rate * parm(2) due to salts (Dove and Rimstidt, MSA Rev. 29, 259)
4  REM  PARM(1) = Specific area of Quartz, m^2/mol Quartz
5  REM  PARM(2) = salt correction: (1 + 1.5 * c_Na (mM)), < 35
10 dif_temp = 1/TK - 1/298
20 pk_w = 13.7 + 4700.4 * dif_temp
40 moles = PARM(1) * M0 * PARM(2) * (M/M0)^0.67 * 10^-pk_w * (1 -  SR("Quartz"))
# Integrate...
50 SAVE moles * TIME
  -end

KINETICS 1
Calcite
    -formula CaCO3 1
    -m 10
    -m0 10
    -parms 2.2661
    -tol 1e-08
Dolomite
     -formula CaMg(CO3)2 1
     -m 10
     -m0 10
     -parms 3.8685
     -tol 1e-08
K-Feldspar
    -formula KAlSi3O8 1
    -m 17.11
    -m0 17.11
    -parms 6.523
    -tol 1e-08
Albite
    -formula NaAlSi3O8 1
    -m 0.5650
    -m0 0.5650
    -parms 6.02
    -tol 1e-08

Quartz
    -formula SiO2 1
    -m 207.05
    -m0 207.05
    -parms 1.36 1
    -tol 1e-08


    -steps    0 1 #86400 172800 345600 691200 2073600 4147200 6220800 8208000 10368000 13392000  # Converted days to seconds


USER_PUNCH
    -headings Time pH SI_Calcite SI_Dolomite SI_K_feldspar SI_Albite SI_Ankerite SI_Quartz SI_Illite new_phi
    -start
    10  REM Initial conditions
    20  phi_0 = 0.1244   # Initial porosity
    40  rho_CO2 = 65.7        # CO2 density at 50?C, 25 bar (kg/m^3)
   
    50  REM Calculate mineral volumes directly
    60  vol_calc = TOT("Calcite")*36.9 # molar volume of calcite = molar mass/ density = 100.09g/mol/ 2.71g/cm3 = 36.9cm3/mol
    80  vol_dol = TOT("Dolomite")*64.365
    90  vol_K_fel = KIN("K-feldspar")*108.87
    100 vol_alb = KIN("Albite")*100.389
    110 vol_ank = KIN("Ankerite")*65.56
    113 vol_qtz = KIN("Quartz")*22.67
    114 vol_ill = KIN("Illite")*141.47

    120 REM Calculate total mineral volume in cm?
    130 vol_total = vol_calc + vol_dol + vol_fel + vol_alb + vol_ank + vol_qtz + vol_ill
   
    140 REM Set representative volume to 1000 cm? #1.0 L
    150 RV = 1000
   
    160 REM Update porosity (ensure it doesn't go negative)
    170 new_phi = phi_0 - vol_total/RV
   
    240 REM Print debugging info to verify calculations
    250 PUT(vol_total, 1)
   
    260 REM Output results with formatting to ensure visibility
    270 punch TOTAL_TIME
    280 punch -LA("H+")
    290 punch SI("Calcite")
    300 punch SI("Dolomite")
    310 punch SI("K-feldspar")
    320 punch SI("Albite")
    330 punch SI("Ankerite")
    340 punch SI("Quartz")
    350 punch SI("Illite")
    360 punch new_phi
    -end

SELECTED_OUTPUT
    -file               EXP_CO2_correction 12.04.2025_003.xls
    -reset              false
    -time               true
    -pH                 true
    -alkalinity         true
    -ionic_strength     true
    -totals             K Na Ca Mg Cl S(6) C(4)
    -molalities  Ca+2 Cl- K+ Mg+2 Na+ SO4-2 CO3-2 HCO3-
    -saturation_indices Calcite Dolomite K-feldspar Albite Ankerite Quartz Illite
    -equilibrium_phases Calcite Dolomite K-feldspar Albite Ankerite Quartz Illite

    -kinetic_reactants Calcite Dolomite K-feldspar Albite Quartz
    -kinetics
    -gases CO2(g) H2O(g)
    -water
    -charge_balance true

PRINT
    -saturation_indices Calcite Dolomite K-feldspar Albite Ankerite Quartz Illite

REACTION 1
CO2(g)
END


[kennedyantwi1]

Code: [Select]

Thanks for getting the model run, Dr. Parkhurst. Really appreciated.
The next issue is, I'm getting the same values in the Excel output and I'm not very sure what input I've missed or wrongly entered.
The current code is attached for your attention and guidance.
Thanks in advence.[SOLUTION_SPECIES
CO3-2 + 2 H+ = CO2 + H2O
     log_k 16.681
     delta_h -5.738 kcal
    -analytical_expression 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
END

PHASES
CO2(g)
CO2 = CO2
     log_k -1.468
     delta_h -4.776 kcal
     -analytical_expression 10.5624 -2.3547e-2 -3972.8 0 5.8746e5 1.9194e-5
     -T_c 304.2 # critical T, K
     -P_c 72.86 # critical P, atm
     -Omega 0.225 # acentric factor
Ankerite
CaMg0.3Fe0.7(CO3)2 + 2H+ = Ca++ + 0.3Mg++ + 0.7Fe++ + 2HCO3-
     log_k 1.54
     delta_h 0
     -analytical_expression -1.8649e+03 -2.9583e-01 1.0468e+05 6.7554e+02 -6.0514e+06
END


SOLUTION 1 #Brine with auto-calculated pH
temp      44
units     mol/kgw
    K         0.08655
    Ca        0.1166     
    Mg        0.01778
    Na        1.305
    Cl        1.2397
    C(4)      0.002
    S(6)      0.003
    water     1.0 #1.0 kg
    pH 7 charge   #PHREEQC will adjust pH to maintain charge balance
    pressure 148.038 atm #Experimental pressure


EQUILIBRIUM_PHASES 1
    CO2(g)        1.398  10.0  # 25 bar PCO2 Partial pressure of carbon dioxide at 50 atm, and injected mol
    Calcite       0.0  0.017742
    Dolomite      0.0  0.0.0321001
    K-feldspar    0.0  0.017014
    Albite        0.0  0.040634
    Ankerite      0.0  0.022945
    Quartz        0.0  0.472921
    Illite        0.0  0.0082214

Save SOLUTION 1
Save EQUILIBRIUM_PHASES 1
END

USE SOLUTION 1
USE EQUILIBRIUM_PHASES 1


RATES
Calcite
   -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_calc = SI("Calcite")
30 if (M <= 0 and si_calc < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_calc > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-0.30)*EXP(-14.40e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(1.000)
70 k_neut = 10^(-5.81)*EXP(-23.50e+03/8.314*(1.0/TK-1.0/298.15))
80 k_carb = 10^(-3.48)*EXP(-35.40e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(1.000)
90 k_rateconst = k_acid + k_neut + k_carb
100 r = k_rateconst * SA * (1-(10^si_calc))
190 moles = r * TIME
200 SAVE moles
   -end

Dolomite
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_dolo = SI("Dolomite")
30 if (M <= 0 and si_dolo < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_dolo > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-3.76)*EXP(-56.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)
70 k_neut = 10^(-8.60)*EXP(-95.30e+03/8.314*(1.0/TK-1.0/298.15))
80 k_carb = 10^(-5.37)*EXP(-45.70e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(0.500)
90 k_rateconst = k_acid + k_neut + k_carb
100 r = k_rateconst * SA * (1-(10^si_dolo))
190 moles = r * TIME
200 SAVE moles
    -end

K-Feldspar
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_kfeld = SI("K-Feldspar")
30 if (M <= 0 and si_kfeld < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_kfeld > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-10.06)*EXP(-51.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)
70 k_neut = 10^(-12.41)*EXP(-38.00e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-21.20)*EXP(-94.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.823)
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_kfeld))
190 moles = r * TIME
200 SAVE moles
    -end

Albite
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_alb = SI("Albite")
30 if (M <= 0 and si_alb < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_alb > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-10.16)*EXP(-65.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("H+")^(0.457))
70 k_neut = 10^(-12.56)*EXP(-69.80e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-15.60)*EXP(-71.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("OH-")^(-0.572))
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA *(1-(10^si_alb))
190 moles = r * TIME
200 SAVE moles
    -end
Quartz
  -start
1  REM  Specific rate k from Rimstidt and Barnes, 1980, GCA 44,1683
2  REM  k = 10^-13.7 mol/m2/s (25 C), Ea = 90 kJ/mol
3  REM  sp. rate * parm(2) due to salts (Dove and Rimstidt, MSA Rev. 29, 259)
4  REM  PARM(1) = Specific area of Quartz, m^2/mol Quartz
5  REM  PARM(2) = salt correction: (1 + 1.5 * c_Na (mM)), < 35
10 dif_temp = 1/TK - 1/298
20 pk_w = 13.7 + 4700.4 * dif_temp
40 moles = PARM(1) * M0 * PARM(2) * (M/M0)^0.67 * 10^-pk_w * (1 -  SR("Quartz"))
# Integrate...
50 SAVE moles * TIME
  -end

KINETICS 1
Calcite
    -formula CaCO3 1
    -m 10
    -m0 10
    -parms 2.2661
    -tol 1e-08
Dolomite
     -formula CaMg(CO3)2 1
     -m 10
     -m0 10
     -parms 3.8685
     -tol 1e-08
K-Feldspar
    -formula KAlSi3O8 1
    -m 17.11
    -m0 17.11
    -parms 6.523
    -tol 1e-08
Albite
    -formula NaAlSi3O8 1
    -m 0.5650
    -m0 0.5650
    -parms 6.02
    -tol 1e-08

Quartz
    -formula SiO2 1
    -m 207.05
    -m0 207.05
    -parms 1.36 1
    -tol 1e-08
KNOBS
-convergence_tolerance 1e-10


USER_PUNCH
    -headings Time pH SI_Calcite SI_Dolomite SI_K_feldspar SI_Albite SI_Ankerite SI_Quartz SI_Illite new_phi
    -start
    10  REM Initial conditions
    20  phi_0 = 0.1244   # Initial porosity
    40  rho_CO2 = 65.7        # CO2 density at 50?C, 25 bar (kg/m^3)
   
    50  REM Calculate mineral volumes directly
    60  vol_calc = TOT("Calcite")*36.9 # molar volume of calcite = molar mass/ density = 100.09g/mol/ 2.71g/cm3 = 36.9cm3/mol
    80  vol_dol = TOT("Dolomite")*64.365
    90  vol_K_fel = KIN("K-feldspar")*108.87
    100 vol_alb = KIN("Albite")*100.389
    110 vol_ank = KIN("Ankerite")*65.56
    113 vol_qtz = KIN("Quartz")*22.67
    114 vol_ill = KIN("Illite")*141.47

    120 REM Calculate total mineral volume in cm?
    130 vol_total = vol_calc + vol_dol + vol_fel + vol_alb + vol_ank + vol_qtz + vol_ill
   
    140 REM Set representative volume to 1000 cm? #1.0 L
    150 RV = 1000
   
    160 REM Update porosity (ensure it doesn't go negative)
    170 new_phi = phi_0 - vol_total/RV
   
    240 REM Print debugging info to verify calculations
    250 PUT(vol_total, 1)
   
    260 REM Output results with formatting to ensure visibility
    270 punch TOTAL_TIME
    280 punch -LA("H+")
    290 punch SI("Calcite")
    300 punch SI("Dolomite")
    310 punch SI("K-feldspar")
    320 punch SI("Albite")
    330 punch SI("Ankerite")
    340 punch SI("Quartz")
    350 punch SI("Illite")
    360 punch new_phi
    -end

SELECTED_OUTPUT
    -file               EXP_CO2_correction 12.04.2025_003_PARKHURST_01.xls
    -reset              false
    -time               true
    -pH                 true
    -alkalinity         true
    -ionic_strength     true
    -totals             K Na Ca Mg Cl S(6) C(4)
    -molalities  Ca+2 Cl- K+ Mg+2 Na+ SO4-2 CO3-2 HCO3-
    -saturation_indices Calcite Dolomite K-feldspar Albite Ankerite Quartz Illite
    -equilibrium_phases Calcite Dolomite K-feldspar Albite Ankerite Quartz Illite

    -kinetic_reactants Calcite Dolomite K-feldspar Albite Quartz
    -kinetics
    -gases CO2(g) H2O(g)
    -water
    -charge_balance true

PRINT
    -saturation_indices Calcite Dolomite K-feldspar Albite Ankerite Quartz Illite


GAS_PHASE 1
-fixed_volume
-volume 0.6
-pressure 148.038 atm
-temperature 44.00


REACTION 1
CO2(g)
20 in 20 steps
ENDcode]

[dlparkhurst]

No time step is defined for KINETICS, so it defaults to 1 second.

However, you have defined both EQUILIBRIUM_PHASES and KINETICS for the same minerals. Think about it. EQUILIBRIUM_PHASES ensures equilibrium with the minerals provided the moles are enough to produce equilibrium. What are the rates for a kinetic reaction when a mineral is at equilibrium?

Similarly, you have CO2 in EQUILIBRIUM_PHASES, which will set the CO2 concentration, regardless of the amount of CO2 added by REACTION.