Conceptual Models > Kinetics and rate controlling factors
Mineral dissolution in carbonated water reaction not proceeding
ReaganK:
Hello,
I am trying to dissolve basaltic glass in carbonated water. I am representing the glass as an already leached amorphous gel of aluminium hydroxide and silica. I have set up the rates and kinectics as seen below. But I cannot seem to get the gel to actually dissolve when I look at the output. Here is the code, what is potentially the problem, please
SOLUTION 1 Vellankatla
temp 45
pH 6.5 charge
pe 4
redox pe
units mmol/kgw
density 1
Al 0.001 as Al3+
C 0.5 as CO2
Ca 0.75
Cl 0.3
Fe 0.002 as Fe2+
K 0.05
Mg 1.1
Na 0.45
S 0.4 as SO4
Si 0.07 as SiO2
-water 1 # kg
EQUILIBRIUM_PHASES 1
CO2(g) 0 -3 dissolve_only
PHASES
basalt_glass_gel
SiAl0.35O2(OH)1.05 + 1.05H+ + 0.95H2O = 0.35Al+3 + H4SiO4
log_k -7.47
Magnesite
MgCO3 + H+ = HCO3- + Mg+2
log_k -4.48
Chalcedony
SiO2 + 2H2O = H4SiO4
log_k -3.93
Allophane
Al2O3SiO2(H2O)2.5 + 6H+ = 2Al+3 + 3.5H2O + H4SiO4
log_k 0
Nontronite
Fe2.5Mg0.5Si4O10(OH)2 + 6H+ + 4H2O = 2.5Fe+2 + 4H4SiO4 + 0.5Mg+2
log_k 0
Saponite
Mg3Si4O10(OH)2 + 6H+ + 4H2O = 4H4SiO4 + 3Mg+2
log_k 0
RATES
glass_dissolution
-start
10 A0 = PARM(3)
20 if (M <= 0) then goto 130
30 R = 8.3144621
40 T = 273.15 + TC
50 area = A0 * (M/M0)**(2/3)
60 rate_const = PARM(1)*EXP(-PARM(2)/(R*T))
70 DF = ((ACT("H+")^3)/(ACT("Al3+")))**(1/3)
80 rate = rate_const * DF * (1 - SUPSAT("basalt_glass_gel"))
90 rate = rate * area
100 moles = rate * time_step
110 PUT(area, 1)
120 PUT(rate, 2)
130 SAVE moles
-end
KINETICS 1 glass_dissolution
glass_dissolution
-formula glass_dissolution 1
-m 1
-m0 1
-parms 2.5e-10 25500 250
-tol 1e-08
-time_step 3600
-step_divide 5
-runge_kutta 3
-bad_step_max 500
REACTION 1 basalt_path
basalt_glass_gel 1
0.5 moles in 100 steps
Kinetics glass_dissolution
equilibrium_phases basalt_glass_gel
SELECTED_OUTPUT 1
-file basalt_glass_gel.sel
-time true
-step true
-pH true
-alkalinity true
-molalities SO4-2 S-2 Al+3
-equilibrium_phases basalt_glass_gel CO2(g)
-saturation_indices Magnesite Calcite Chalcedony Allophane
Saponite K-feldspar
-kinetic_reactants glass_dissolution
END
Plot
import pandas as pd
import matplotlib.pyplot as plt
df = pd.read_csv(
'basalt_glass_gel.sel',
delim_whitespace=True,
comment='#'
)
TIME_STEP = 3600 # seconds per step as defined in your KINETICS block
df['time_s'] = df['step'] * TIME_STEP
df['time_h'] = df['time_s'] / 3600 # optional hours axis
df_rxn = df[df['step'] > 0].copy()
plt.figure()
plt.plot(
df_rxn['time_s'],
df_rxn['k_glass_dissolution'],
marker='o'
)
plt.xlabel('Time (s)')
plt.ylabel('Dissolution rate (mol/s)')
plt.title('Glass Dissolution Rate')
plt.grid(True)
plt.show()
df_rxn['cumulative_moles'] = df_rxn['dk_glass_dissolution'].cumsum()
plt.figure()
plt.plot(
df_rxn['time_s'],
df_rxn['cumulative_moles'],
marker='o'
)
plt.xlabel('Time (s)')
plt.ylabel('Cumulative moles dissolved')
plt.title('Total Glass Dissolved Over Time')
plt.grid(True)
plt.show()
plt.figure()
plt.plot(df_rxn['time_s'], df_rxn['pH'], marker='o')
plt.xlabel('Time (s)')
plt.ylabel('pH')
plt.title('pH Evolution Over Time')
plt.grid(True)
plt.show()
si_cols = [col for col in df_rxn.columns if col.startswith('si_')]
plt.figure()
for col in si_cols:
plt.plot(df_rxn['time_s'], df_rxn[col], label=col[3:]) # strip the 'si_' prefix
plt.xlabel('Time (s)')
plt.ylabel('Saturation Index')
plt.title('Mineral Saturation Indices Over Time')
plt.legend()
plt.grid(True)
plt.show()
dlparkhurst:
I'm always surprised when you are having trouble that you set steps to 100. Use one short step until you have things working.
These two statements redefine KINETICS 1 to have no kinetic reactions and EQUILIBRIUM_PHASES 1 to have no equilibrium phases. I think you want to eliminate these lines.
--- Code: --- Kinetics glass_dissolution
equilibrium_phases basalt_glass_gel
--- End code ---
Here you have defined a partial pressure of CO2 of 1 atm (10^0) with -3 moles available to react. With no moles available, there is no way for the phase to dissolve.
--- Code: ---EQUILIBRIUM_PHASES 1
CO2(g) 0 -3 dissolve_only
--- End code ---
The Basic operator for exponentiation is ^, not **.
There is not Basic function SUPSAT, perhaps you meant SR. See The Basic Interpreter in the documentation.
If you fix these things, you may still have errors. Use PRINT statements in RATES to debug any errors in the calculation of the rate expression.
ReaganK:
Hello, thanks a lot fr the previous response. I have tried amending the code accordingly. I have increased the total Al, changed dissolving phase to basalt_glass_solid, handling the dissolution in the RATES and KINETIC BLOCKS only (not in the EQUILIBRIUM_PHASE or REACTION blocks) debuging the RATES BLOCK with PRINT statements, etc. However, I still have no dissolution and it looks like phreeqc is not even reaching the RATES block before the simulation ends. Kindly have a look and and advise. Thanks
SOLUTION 1 Vellankatla
temp 45
pH 7.8 charge
pe 4
redox pe
units umol/kgw
density 1
Al 10.9 as Al3+
C 354 as CO2
Ca 71
Cl 120
Fe 0.16
K 11.9
Mg 38
Na 269
S 15 as SO4
Si 256 as H4SiO4
-water 1 # kg
PHASES
basalt_glass_gel
SiAl0.35O2(OH)1.05 + 1.05H+ + 0.95H2O = 0.35Al+3 + H4SiO4
log_k -7.47
Magnesite
MgCO3 + H+ = HCO3- + Mg+2
log_k -4.48
Allophane
Al2O3SiO2(H2O)2.5 + 6H+ = 2Al+3 + 3.5H2O + H4SiO4
log_k 0
Nontronite
Na.66Fe4Al.66Si7.34O24H4 + 14.64H+ + 5.36H2O = 0.66Al+3 + 4Fe+3 + 7.34H4SiO4 + 0.66Na+
log_k -21.188
Saponite-Ca
Ca.165Mg3Al.33Si3.67O10(OH)2 + 7.32H+ + 2.68H2O = 0.33Al+3 + 0.165Ca+2 + 3.67H4SiO4 + 3Mg+2
log_k 26.29
delta_h -207.971 kJ
-analytical_expression -46.904 0.0062555 22572 5.3198 -1572500 0
Saponite-H
H.33Mg3Al.33Si3.67O10(OH)2 + 6.99H+ + 2.68H2O = 0.33Al+3 + 3.67H4SiO4 + 3Mg+2
log_k 25.3321
delta_h -200.235 kJ
-analytical_expression -39.828 0.0089566 22165 2.3941 -1593300 0
Saponite-K
K.33Mg3Al.33Si3.67O10(OH)2 + 7.32H+ + 2.68H2O = 0.33Al+3 + 3.67H4SiO4 + 0.33K+ + 3Mg+2
log_k 26.0075
delta_h -196.402 kJ
-analytical_expression 32.113 0.018392 17918 -22.874 -1354200 0
Saponite-Mg
Mg3.165Al.33Si3.67O10(OH)2 + 7.32H+ + 2.68H2O = 0.33Al+3 + 3.67H4SiO4 + 3.165Mg+2
log_k 26.2523
delta_h -210.822 kJ
-analytical_expression 9.8888 0.01432 19418 -15.259 -1371600 0
Saponite-Na
Na.33Mg3Al.33Si3.67O10(OH)2 + 7.32H+ + 2.68H2O = 0.33Al+3 + 3.67H4SiO4 + 3Mg+2 + 0.33Na+
log_k 26.3459
delta_h -201.401 kJ
-analytical_expression -67.611 0.0047327 23586 12.868 -1649300 0
basalt_glass_solid
K0.008Na0.08Ca0.27Mg0.26Mn0.003Fe0.18Al0.35SiO3.282 + 2.564H+ + 0.718H2O = 0.35Al+3 + 0.27Ca+2 + 0.18Fe+2 + H4SiO4 + 0.008K+ + 0.26Mg+2 + 0.003Mn+2 + 0.08Na+
log_k -15
END
EQUILIBRIUM_PHASES 1
CO2(g) -3.5 0
END
RATES
glass_dissolution
-start
10 A0 = PARM(3)
20 IF (M <= 0) THEN
30 moles = 0
40 GOTO 210
50 ENDIF
60 R = 8.3144621
70 T = 273.15 + TC
80 area = A0 * (M / M0) ^ (2/3)
90 rate_const = PARM(1) * EXP(-PARM(2) / (R * T))
100 DF = ((ACT("H+") ^ 3) / (ACT("Al3+"))) ^ (1/3)
110 rate = rate_const * DF * (1 - SR("basalt_glass_solid"))
120 rate = rate * area
130 PRINT "DEBUG_TIME=", TIME, "M=", M, "A0=", A0
140 PRINT "R=", R, " T=", T, " area=", area, "ACT H+=",ACT("H+"), "ACT Al3+=", ACT("Al3+")
150 PRINT "rate_const=", rate_const, " DF=", DF
160 PRINT "SR(basalt_glass_solid)=", SR("basalt_glass_solid")
170 PRINT "rate (before TIME)=", rate
180 moles = rate * TIME
190 PUT(area, 1)
200 PUT(rate, 2)
210 SAVE moles
-end
KINETICS 1 glass_dissolution
glass_dissolution
-formula basalt_glass_solid 1
-m 5
-m0 5
-parms 2.51189e-06 25500 15.1625
-tol 1e-08
-steps 3600
-step_divide 5
-runge_kutta 3
-bad_step_max 500
SELECTED_OUTPUT 1
-file basalt_glass_solid.sel
-time true
-step true
-pH true
-alkalinity true
-molalities SO4-2 S-2 Al+3
-equilibrium_phases CO2(g)
-saturation_indices Magnesite Calcite Chalcedony Allophane
K-feldspar Nontronite Ca-Montmorillonite Saponite-Mg
basalt_glass_gel Illite
-kinetic_reactants glass_dissolution
END
dlparkhurst:
Calculations are done when an END is encountered.
Calculations only occur when there is (1) a SOLUTION definition, (2) a USE solution n definition, (3) a MIX definition, or (4) a TRANSPORT definition before the next END statement.
MichaelZ20:
Your rate is very high because of the DF that equal 1.1e25
Navigation
[0] Message Index
[#] Next page
Go to full version