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: [Select]

    Kinetics glass_dissolution
    equilibrium_phases basalt_glass_gel

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: [Select]

EQUILIBRIUM_PHASES 1
    CO2(g)    0 -3 dissolve_only

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

[ReaganK]

Hi, thanks for the response. I cleaned it up a bit but still no dissolution rxn. However,

DF (DF = ((ACT("H+") ^ 3) / (ACT("Al3+"))) ^ (1/3)) is actaully

1.77E-03 when considering speciation of solution with
aH+   3.21E-08
aAl3+   5.99E-15
aH4SiO4   2.54E-04

and 2.52E-03 after equilibriating with mentiond CO2 phase (before rxn step 1), with
aH+   1.11E-08
aAl3+   8.55E-17
aH4SiO4   2.54E-04

So in both cases DF is far from --> infinity. At least a resonable value is expected in the equation,  rate = rate_const * DF * (1 - SR("basalt_glass_solid"))

Infact when I calculate the rate on a spreadsheet I get a value so close to what  Gıslason and Oelkers, 2003 have (https://www.sciencedirect.com/science/article/pii/S0016703703001765#APP1), on page 3831. See calculation below, or attached spreadsheet, sheet3 (Attachment could not load)


Ageo   250   cm2/g
kgeo   2.51E-06   
EA   25500   kJ mol^-1
R   8.3145   
T   318.15   
R*T   2645.258175   
rate constant   1.63474E-10
   
RXN STEP 1 SOLID      
aH+   1.11E-08   
aAl3+   8.55E-17   
aH4SiO4   2.54E-04   
      
DF   2.52E-03   
Q(IAP of  SiAl0.35O2(OH)1.05)   4.670E-18   
Log K    25   
K   1E+25   
SR = Q/K   4.7E-43   
rate   4.12E-13        # log(rate = 4.12E-12) = -12.38, this is not far from what Gıslason and Oelkers, 2003 have on page 3831.

However, I am not sure about the ACT(Al3+). (PS: I also tried using DF = (ACT("H+")^1, and rate = rate_const * DF to try to minimise the influence of AL, but this changed not alot)

I have not calculated the log K of the glass as per now, but even when i use arbitrary numbers esp to favour forward rxn like log k = 25, there is still no dissolution.

Please look at the code below, what am I missing?

SOLUTION 1 Vellankatla
    temp      45
    pH        7.84 charge
    pe        4
    redox     pe
    units     umol/kgw
    density   1
    Al        10.9
    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    25


EQUILIBRIUM_PHASES 1
    CO2(g)    -3.5 0


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 151625
    -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

[MichaelZ20]

My fault

Code: [Select]

SOLUTION 1 Vellankatla
    temp      45
    pH        7.84 charge
    pe        4
    redox     pe
    units     umol/kgw
    density   1
    Al        10.9
    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    25


EQUILIBRIUM_PHASES 1
    CO2(g)    -3.5 0


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("Al+3")))^(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
220 PRINT "DF = ", ((ACT("H+")^3) / (ACT("Al+3")))^(1/3)
230 PRINT "Rate = ", (PARM(1) * EXP(-PARM(2) / (8.3144621 * (273.15 + TC))))*(((ACT("H+")^3) / (ACT("Al+3")))^(1/3)) * (1 - SR("basalt_glass_solid"))
240 PRINT "ACT H+ = ",ACT("H+")
250 PRINT "ACT Al+3 = ", ACT("Al+3")
260 PRINT "T = ", TC + 273.15
270
-end

KINETICS 1 glass_dissolution
glass_dissolution
    -formula  basalt_glass_solid  1
    -m        5
    -m0       5
    -parms    2.51189e-06 25500 151625
    -tol      1e-08
-time       315360000 in 60 steps # 10 years
-step_divide 5
-runge_kutta 3
-bad_step_max 500
END

[dlparkhurst]

Please use the # button above the input text box to set off your script.

The initial Basic statements to check for M=0 are incorrect. It is actually unnecessary to check for M=0 because the reaction will stop if it runs out of moles of reactant and moles is initialized to zero. The following is correct syntax:

Code: [Select]

RATES
    glass_dissolution
-start
 10 A0 = PARM(3)
 20 IF (M <= 0) THEN GOTO 210
 60 R = 8.3144621
...


Now you can work on your rate. All of the glass dissolves in 1 second.

[dlparkhurst]

The bare aluminum ion is Al+3, not Al3+, the plus or minus signs for ions go before the charge.