Conceptual Models > Design of conceptual models
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: ---[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]
--- End code ---
dlparkhurst:
--- Code: ---pH 7 charge #PHREEQC will adjust pH to maintain charge balance
--- End code ---
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: ---[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
--- End code ---
kennedyantwi1:
--- Code: ---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]
--- End code ---
Navigation
[0] Message Index
[#] Next page
Go to full version