Graph of pH/ time when cement leachate interact with soil minerals

[helenhuntpvt@gmail.com]

Hello

I feel so good coming across this platform on the web as I have been able to learn a lot. I have been trying to model the interaction of hyper-alkaline leachate comprising of young cement leachate with a pH of 13.1, intermediate leachate with a pH of 12.3 and old leachate with a pH of 10.4 interacting with different soil minerals as seen in the code. All effort has not been productive as I am not able to produce a graph of pH/time for this mineral to aid my proposed groundwater remediation strategy.  Please help my check my user graph prompt to know where the error lies.

ps I learnt using phreeqc code using this platform. Thanks  [/quote]

Code: [Select]

SOLUTION 2 Young Cement Leachate (YCL)
    temp      25
    pH        13.1
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    K         0.092 mol/kgw
    Na        0.095 mol/kgw
    Ca        0.000135 mol/kgw

    -water    1 # kg

End
SOLUTION 3 Intermediate Leachate (IL)
    temp      25
    pH        12.3
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Ca        0.0162 mol/kgw
    Cl        0.004141 mol/kgw
    K         0.00397 mol/kgw
    Na        0.000171 mol/kgw
    -water    1 # kg

End
SOLUTION 4 Old Leachate (OL)
    temp      25
    pH        10.4
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Ca        0.000202 mol/kgw
    -water    1 # kg

End

Mix 1
    2  1
    3  1
    4  1


SOLUTION 1

RATES
 Quartz
 # d qu / dt = (A0 / V) * (m / m0)^0.67 * k * (1 - SI_qu)
 # Specific rate k from Rimstidt and Barnes, 1980, GCA 44, 1683
 #  k = 10^-13.7 mol/m2/s (25 C)
 #  A0, initial surface of quartz (m2) recalc's to mol/s
 #  V, solution volume in contact with A0
   -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Quartz") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Quartz") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0                 
    70 k_neut = 1.02e-14*EXP((-87700/8.3145)*(1/TK-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-SR("Quartz"))
    190 moles = r * TIME
    200 SAVE moles
     -end

K-Feldspar
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("K-Feldspar") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("K-Feldspar") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 8.71e-11*EXP((-51700/8.3145)*(1/Tk-1/298.15))                 
    70 k_neut = 3.89e-13*EXP((-38000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 6.31e-22*EXP((-94100/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.5 + k_neut + k_base*act("H+")^(-0.823)
    100 r = k_rateconst * SA * (1-SR("K-Feldspar"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Muscovite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Muscovite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Muscovite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 1.41e-12*EXP((-22000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-14*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.82e-15*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.37 + k_neut + k_base*act("H+")^(-0.22)
    100 r = k_rateconst * SA * (1-SR("Muscovite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Kaolinite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Kaolinite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Kaolinite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 7.7e-14*EXP((-65900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 7.7e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
    80 k_base = 4.70e-18*EXP((-17900/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.777 + k_neut + k_base*act("H+")^(-0.472)
    100 r = k_rateconst * SA * (1-SR("Kaolinite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Albite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Albite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Albite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 6.92e-11*EXP((-65000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.75e-13*EXP((-69800/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.51e-16*EXP((-71000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.457 + k_neut + k_base*act("H+")^(-0.572)
    100 r = k_rateconst * SA * (1-SR("Albite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Pyrite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Pyrite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Pyrite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 3.02e-8*EXP((-56900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-5*EXP((-56900/8.3145)*(1/Tk-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid*(act("H+")^(-0.5))*(act("Fe+2")^(0.5))  + k_neut*act("O2")^0.5  + k_base*act("H+")^0
    100 r = (k_rateconst) * SA * (1-SR("Pyrite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Calcite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Calcite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Calcite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0.5012*EXP(-14.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")                 
    70 k_neut = 1.5488e-6*EXP(-23.5e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 3.31e-4*EXP(-35.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("CO2")
    90 k_rateconst = k_acid*act("H+")^1  + k_neut + k_base*act("H+")^1
    100 r = k_rateconst * SA * (1-SR("Calcite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
KINETICS 1
Quartz
 -formula SiO2
 -m0  158.8                        # initial moles of quartz
 -parms 23.13  0.16
 -time_step 5 year in 15           # 15 time steps
INCREMENTAL_REACTIONS true
KINETICS 1
Quartz
    -formula  SiO2  1
    -m        158.3
    -m0       158.3
    -parms    4.122
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        21.52
    -m0       21.52
    -parms    299.027
    -tol      1e-08
K-Feldspar
    -formula  KAlSi3O8  1
    -m        20.32
    -m0       20.32
    -parms    19.768
    -tol      1e-08
Muscovite     
    -formula  KAl3Si3O10(OH)2  1
    -m        11.4
    -m0       11.4
    -parms    25.607
    -tol      1e-08
Calcite       
    -formula  CaCO3  1
    -m        27.248
    -m0       27.248
    -parms    6.715
    -tol      1e-08

-steps   100*30 10*300  #31.536 3153.6 #3153.6 31536 31536 31536 31536 31536 315360 315360 3153600 3153600 3153600 3153600 3153600 3153600 3153600 31536000 31536000 31536000 31536000 31536000 31536000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 103280400
-step_divide 1
-runge_kutta 3
-bad_step_max 4000
-cvode true
-cvode_steps 500
-cvode_order 5

USER_GRAPH 1
    -headings               time Quartz Kaolinite K-Feldspar Muscovite Calcite Pyrite
    -axis_titles            "pH" "Time" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 graph_x total_pH
20 GRAPH_Y SI("Quartz"), SI("Kaolinite"), SI("K-Feldspar"), SI("Muscovite"), SI("Calcite")
  -end
    -active                 true

END

#KINETICS 1
#K-feldspar
# k0 * A/V = 1e-16 mol/cm2/s * (10% fsp, 0.1mm cubes) 136/cm = 136.e-13 mol/dm3/s
        -parms 1.36e-11
        -m0 2.16
        -m  1.94
        -time_step 5 year in 15           # 15 time steps
#        -step_divide 1e-6
#        -steps    1e2 1e3 1e4 1e5 1e6 1e7 1e8
#        -steps    1e2 1e3 1e4 1e5 63240.0 64950.0 1347610.0 1010300.0 45242800.0

#INCREMENTAL_REACTIONS true
#RATES
#K-feldspar
#-start
  10  REM store the initial amount of K-feldspar
  20  IF EXISTS(1) = 0 THEN PUT(M, 1)
  30  REM calculate moles of reaction
  40  SR_kfld = SR("K-feldspar")
  50  moles = PARM(1) * (M/M0)^0.67 * (1 - SR_kfld) * TIME
  60  REM The following is for printout of phase transitions
  80  REM      Start Gibbsite
  90  if ABS(SI("Gibbsite")) > 1e-3 THEN GOTO 150
  100   i = 2
  110   GOSUB 1500
  150 REM      Start Gibbsite -> Kaolinite
  160 if ABS(SI("Kaolinite")) > 1e-3 THEN GOTO 200
  170   i = 3
  180   GOSUB 1500
  200 REM      End Gibbsite -> Kaolinite
  210 if ABS(SI("Kaolinite")) > 1e-3 OR EQUI("Gibbsite") > 0 THEN GOTO 250
  220   i = 4
  230   GOSUB 1500
  250 REM      Start Kaolinite -> K-mica
  260 if ABS(SI("K-mica")) > 1e-3 THEN GOTO 300
  270   i = 5
  280   GOSUB 1500
  300 REM      End Kaolinite -> K-mica
  310 if ABS(SI("K-mica")) > 1e-3 OR EQUI("Kaolinite") > 0 THEN GOTO 350
  320   i = 6
  330   GOSUB 1500
  350 REM      Start K-mica -> K-feldspar
  360 if ABS(SI("K-feldspar")) > 1e-3 THEN GOTO 1000
  370   i = 7
  380   GOSUB 1500
  1000 SAVE moles
  1010 END
  1500 REM subroutine to store data
  1510 if GET(i) >= M THEN RETURN
  1520     PUT(M, i)
  1530     PUT(TOTAL_TIME, i, 1)
  1540     PUT(LA("K+")-LA("H+"), i, 2)
  1550     PUT(LA("H4SiO4"), i, 3)
  1560 RETURN
-end







[dlparkhurst]

pH is -LA("H+"), negative log activity of H+.

You only get one KINETICS definition per calculation.

Code: [Select]

SOLUTION 2 Young Cement Leachate (YCL)
    temp      25
    pH        13.1
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    K         0.092 mol/kgw
    Na        0.095 mol/kgw
    Ca        0.000135 mol/kgw

    -water    1 # kg

End
SOLUTION 3 Intermediate Leachate (IL)
    temp      25
    pH        12.3
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Ca        0.0162 mol/kgw
    Cl        0.004141 mol/kgw
    K         0.00397 mol/kgw
    Na        0.000171 mol/kgw
    -water    1 # kg

End
SOLUTION 4 Old Leachate (OL)
    temp      25
    pH        10.4
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Ca        0.000202 mol/kgw
    -water    1 # kg

End

Mix 1
    2  1
    3  1
    4  1


SOLUTION 1

RATES
 Quartz
 # d qu / dt = (A0 / V) * (m / m0)^0.67 * k * (1 - SI_qu)
 # Specific rate k from Rimstidt and Barnes, 1980, GCA 44, 1683
 #  k = 10^-13.7 mol/m2/s (25 C)
 #  A0, initial surface of quartz (m2) recalc's to mol/s
 #  V, solution volume in contact with A0
   -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Quartz") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Quartz") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0                 
    70 k_neut = 1.02e-14*EXP((-87700/8.3145)*(1/TK-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-SR("Quartz"))
    190 moles = r * TIME
    200 SAVE moles
     -end

K-Feldspar
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("K-Feldspar") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("K-Feldspar") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 8.71e-11*EXP((-51700/8.3145)*(1/Tk-1/298.15))                 
    70 k_neut = 3.89e-13*EXP((-38000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 6.31e-22*EXP((-94100/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.5 + k_neut + k_base*act("H+")^(-0.823)
    100 r = k_rateconst * SA * (1-SR("K-Feldspar"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Muscovite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Muscovite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Muscovite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 1.41e-12*EXP((-22000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-14*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.82e-15*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.37 + k_neut + k_base*act("H+")^(-0.22)
    100 r = k_rateconst * SA * (1-SR("Muscovite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Kaolinite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Kaolinite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Kaolinite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 7.7e-14*EXP((-65900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 7.7e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
    80 k_base = 4.70e-18*EXP((-17900/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.777 + k_neut + k_base*act("H+")^(-0.472)
    100 r = k_rateconst * SA * (1-SR("Kaolinite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Albite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Albite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Albite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 6.92e-11*EXP((-65000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.75e-13*EXP((-69800/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.51e-16*EXP((-71000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.457 + k_neut + k_base*act("H+")^(-0.572)
    100 r = k_rateconst * SA * (1-SR("Albite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Pyrite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Pyrite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Pyrite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 3.02e-8*EXP((-56900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-5*EXP((-56900/8.3145)*(1/Tk-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid*(act("H+")^(-0.5))*(act("Fe+2")^(0.5))  + k_neut*act("O2")^0.5  + k_base*act("H+")^0
    100 r = (k_rateconst) * SA * (1-SR("Pyrite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Calcite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Calcite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Calcite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0.5012*EXP(-14.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")                 
    70 k_neut = 1.5488e-6*EXP(-23.5e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 3.31e-4*EXP(-35.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("CO2")
    90 k_rateconst = k_acid*act("H+")^1  + k_neut + k_base*act("H+")^1
    100 r = k_rateconst * SA * (1-SR("Calcite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
#KINETICS 1
#Quartz
# -formula SiO2
# -m0  158.8                        # initial moles of quartz
# -parms 23.13  0.16
# -time_step 5 year in 15           # 15 time steps
INCREMENTAL_REACTIONS true
KINETICS 1
Quartz
    -formula  SiO2  1
    -m        158.3
    -m0       158.3
    -parms    4.122
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        21.52
    -m0       21.52
    -parms    299.027
    -tol      1e-08
K-Feldspar
    -formula  KAlSi3O8  1
    -m        20.32
    -m0       20.32
    -parms    19.768
    -tol      1e-08
Muscovite     
    -formula  KAl3Si3O10(OH)2  1
    -m        11.4
    -m0       11.4
    -parms    25.607
    -tol      1e-08
Calcite       
    -formula  CaCO3  1
    -m        27.248
    -m0       27.248
    -parms    6.715
    -tol      1e-08

-steps   100*30 10*300  #31.536 3153.6 #3153.6 31536 31536 31536 31536 31536 315360 315360 3153600 3153600 3153600 3153600 3153600 3153600 3153600 31536000 31536000 31536000 31536000 31536000 31536000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 103280400
-step_divide 1
-runge_kutta 3
-bad_step_max 4000
-cvode true
-cvode_steps 500
-cvode_order 5

USER_GRAPH 1
    -headings               time Quartz Kaolinite K-Feldspar Muscovite Calcite Pyrite
    -axis_titles            "pH" "Time" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 graph_x -LA("H+")
20 GRAPH_Y SI("Quartz"), SI("Kaolinite"), SI("K-Feldspar"), SI("Muscovite"), SI("Calcite")
  -end
    -active                 true

END

#KINETICS 1
#K-feldspar
# k0 * A/V = 1e-16 mol/cm2/s * (10% fsp, 0.1mm cubes) 136/cm = 136.e-13 mol/dm3/s
        -parms 1.36e-11
        -m0 2.16
        -m  1.94
        -time_step 5 year in 15           # 15 time steps
#        -step_divide 1e-6
#        -steps    1e2 1e3 1e4 1e5 1e6 1e7 1e8
#        -steps    1e2 1e3 1e4 1e5 63240.0 64950.0 1347610.0 1010300.0 45242800.0

#INCREMENTAL_REACTIONS true
#RATES
#K-feldspar
#-start
  10  REM store the initial amount of K-feldspar
  20  IF EXISTS(1) = 0 THEN PUT(M, 1)
  30  REM calculate moles of reaction
  40  SR_kfld = SR("K-feldspar")
  50  moles = PARM(1) * (M/M0)^0.67 * (1 - SR_kfld) * TIME
  60  REM The following is for printout of phase transitions
  80  REM      Start Gibbsite
  90  if ABS(SI("Gibbsite")) > 1e-3 THEN GOTO 150
  100   i = 2
  110   GOSUB 1500
  150 REM      Start Gibbsite -> Kaolinite
  160 if ABS(SI("Kaolinite")) > 1e-3 THEN GOTO 200
  170   i = 3
  180   GOSUB 1500
  200 REM      End Gibbsite -> Kaolinite
  210 if ABS(SI("Kaolinite")) > 1e-3 OR EQUI("Gibbsite") > 0 THEN GOTO 250
  220   i = 4
  230   GOSUB 1500
  250 REM      Start Kaolinite -> K-mica
  260 if ABS(SI("K-mica")) > 1e-3 THEN GOTO 300
  270   i = 5
  280   GOSUB 1500
  300 REM      End Kaolinite -> K-mica
  310 if ABS(SI("K-mica")) > 1e-3 OR EQUI("Kaolinite") > 0 THEN GOTO 350
  320   i = 6
  330   GOSUB 1500
  350 REM      Start K-mica -> K-feldspar
  360 if ABS(SI("K-feldspar")) > 1e-3 THEN GOTO 1000
  370   i = 7
  380   GOSUB 1500
  1000 SAVE moles
  1010 END
  1500 REM subroutine to store data
  1510 if GET(i) >= M THEN RETURN
  1520     PUT(M, i)
  1530     PUT(TOTAL_TIME, i, 1)
  1540     PUT(LA("K+")-LA("H+"), i, 2)
  1550     PUT(LA("H4SiO4"), i, 3)
  1560 RETURN
-end


[helenhuntpvt@gmail.com]

 Thank you for your response. However, I was expecting to see how the different pH of the leachate interact with the soil minerals  over time after mixing them with solution 1 containing the kinetics of the different minerals. "You only get one KINETICS definition per calculation" does this means I have to produce different user graph for the kinetic of each mineral.

Additionally, is there any user graph prompt to produce key chemical element from these minerals?

Thank you   

pH is -LA("H+"), negative log activity of H+.

You only get one KINETICS definition per calculation.

Code: [Select]

SOLUTION 2 Young Cement Leachate (YCL)
    temp      25
    pH        13.1
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    K         0.092 mol/kgw
    Na        0.095 mol/kgw
    Ca        0.000135 mol/kgw

    -water    1 # kg

End
SOLUTION 3 Intermediate Leachate (IL)
    temp      25
    pH        12.3
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Ca        0.0162 mol/kgw
    Cl        0.004141 mol/kgw
    K         0.00397 mol/kgw
    Na        0.000171 mol/kgw
    -water    1 # kg

End
SOLUTION 4 Old Leachate (OL)
    temp      25
    pH        10.4
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Ca        0.000202 mol/kgw
    -water    1 # kg

End

Mix 1
    2  1
    3  1
    4  1


SOLUTION 1

RATES
 Quartz
 # d qu / dt = (A0 / V) * (m / m0)^0.67 * k * (1 - SI_qu)
 # Specific rate k from Rimstidt and Barnes, 1980, GCA 44, 1683
 #  k = 10^-13.7 mol/m2/s (25 C)
 #  A0, initial surface of quartz (m2) recalc's to mol/s
 #  V, solution volume in contact with A0
   -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Quartz") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Quartz") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0                 
    70 k_neut = 1.02e-14*EXP((-87700/8.3145)*(1/TK-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-SR("Quartz"))
    190 moles = r * TIME
    200 SAVE moles
     -end

K-Feldspar
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("K-Feldspar") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("K-Feldspar") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 8.71e-11*EXP((-51700/8.3145)*(1/Tk-1/298.15))                 
    70 k_neut = 3.89e-13*EXP((-38000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 6.31e-22*EXP((-94100/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.5 + k_neut + k_base*act("H+")^(-0.823)
    100 r = k_rateconst * SA * (1-SR("K-Feldspar"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Muscovite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Muscovite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Muscovite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 1.41e-12*EXP((-22000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-14*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.82e-15*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.37 + k_neut + k_base*act("H+")^(-0.22)
    100 r = k_rateconst * SA * (1-SR("Muscovite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Kaolinite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Kaolinite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Kaolinite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 7.7e-14*EXP((-65900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 7.7e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
    80 k_base = 4.70e-18*EXP((-17900/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.777 + k_neut + k_base*act("H+")^(-0.472)
    100 r = k_rateconst * SA * (1-SR("Kaolinite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Albite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Albite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Albite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 6.92e-11*EXP((-65000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.75e-13*EXP((-69800/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.51e-16*EXP((-71000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.457 + k_neut + k_base*act("H+")^(-0.572)
    100 r = k_rateconst * SA * (1-SR("Albite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Pyrite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Pyrite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Pyrite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 3.02e-8*EXP((-56900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-5*EXP((-56900/8.3145)*(1/Tk-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid*(act("H+")^(-0.5))*(act("Fe+2")^(0.5))  + k_neut*act("O2")^0.5  + k_base*act("H+")^0
    100 r = (k_rateconst) * SA * (1-SR("Pyrite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Calcite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Calcite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Calcite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0.5012*EXP(-14.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")                 
    70 k_neut = 1.5488e-6*EXP(-23.5e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 3.31e-4*EXP(-35.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("CO2")
    90 k_rateconst = k_acid*act("H+")^1  + k_neut + k_base*act("H+")^1
    100 r = k_rateconst * SA * (1-SR("Calcite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
#KINETICS 1
#Quartz
# -formula SiO2
# -m0  158.8                        # initial moles of quartz
# -parms 23.13  0.16
# -time_step 5 year in 15           # 15 time steps
INCREMENTAL_REACTIONS true
KINETICS 1
Quartz
    -formula  SiO2  1
    -m        158.3
    -m0       158.3
    -parms    4.122
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        21.52
    -m0       21.52
    -parms    299.027
    -tol      1e-08
K-Feldspar
    -formula  KAlSi3O8  1
    -m        20.32
    -m0       20.32
    -parms    19.768
    -tol      1e-08
Muscovite     
    -formula  KAl3Si3O10(OH)2  1
    -m        11.4
    -m0       11.4
    -parms    25.607
    -tol      1e-08
Calcite       
    -formula  CaCO3  1
    -m        27.248
    -m0       27.248
    -parms    6.715
    -tol      1e-08

-steps   100*30 10*300  #31.536 3153.6 #3153.6 31536 31536 31536 31536 31536 315360 315360 3153600 3153600 3153600 3153600 3153600 3153600 3153600 31536000 31536000 31536000 31536000 31536000 31536000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 103280400
-step_divide 1
-runge_kutta 3
-bad_step_max 4000
-cvode true
-cvode_steps 500
-cvode_order 5

USER_GRAPH 1
    -headings               time Quartz Kaolinite K-Feldspar Muscovite Calcite Pyrite
    -axis_titles            "pH" "Time" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 graph_x -LA("H+")
20 GRAPH_Y SI("Quartz"), SI("Kaolinite"), SI("K-Feldspar"), SI("Muscovite"), SI("Calcite")
  -end
    -active                 true

END

#KINETICS 1
#K-feldspar
# k0 * A/V = 1e-16 mol/cm2/s * (10% fsp, 0.1mm cubes) 136/cm = 136.e-13 mol/dm3/s
        -parms 1.36e-11
        -m0 2.16
        -m  1.94
        -time_step 5 year in 15           # 15 time steps
#        -step_divide 1e-6
#        -steps    1e2 1e3 1e4 1e5 1e6 1e7 1e8
#        -steps    1e2 1e3 1e4 1e5 63240.0 64950.0 1347610.0 1010300.0 45242800.0

#INCREMENTAL_REACTIONS true
#RATES
#K-feldspar
#-start
  10  REM store the initial amount of K-feldspar
  20  IF EXISTS(1) = 0 THEN PUT(M, 1)
  30  REM calculate moles of reaction
  40  SR_kfld = SR("K-feldspar")
  50  moles = PARM(1) * (M/M0)^0.67 * (1 - SR_kfld) * TIME
  60  REM The following is for printout of phase transitions
  80  REM      Start Gibbsite
  90  if ABS(SI("Gibbsite")) > 1e-3 THEN GOTO 150
  100   i = 2
  110   GOSUB 1500
  150 REM      Start Gibbsite -> Kaolinite
  160 if ABS(SI("Kaolinite")) > 1e-3 THEN GOTO 200
  170   i = 3
  180   GOSUB 1500
  200 REM      End Gibbsite -> Kaolinite
  210 if ABS(SI("Kaolinite")) > 1e-3 OR EQUI("Gibbsite") > 0 THEN GOTO 250
  220   i = 4
  230   GOSUB 1500
  250 REM      Start Kaolinite -> K-mica
  260 if ABS(SI("K-mica")) > 1e-3 THEN GOTO 300
  270   i = 5
  280   GOSUB 1500
  300 REM      End Kaolinite -> K-mica
  310 if ABS(SI("K-mica")) > 1e-3 OR EQUI("Kaolinite") > 0 THEN GOTO 350
  320   i = 6
  330   GOSUB 1500
  350 REM      Start K-mica -> K-feldspar
  360 if ABS(SI("K-feldspar")) > 1e-3 THEN GOTO 1000
  370   i = 7
  380   GOSUB 1500
  1000 SAVE moles
  1010 END
  1500 REM subroutine to store data
  1510 if GET(i) >= M THEN RETURN
  1520     PUT(M, i)
  1530     PUT(TOTAL_TIME, i, 1)
  1540     PUT(LA("K+")-LA("H+"), i, 2)
  1550     PUT(LA("H4SiO4"), i, 3)
  1560 RETURN
-end


[dlparkhurst]

Sorry, I don't understand what you are trying to do. In particular, I don't understand this statement " solution 1 containing the kinetics of the different minerals". A solution can be the result of of a kinetic reaction, but the solution itself does not contain any of the kinetic reactants. SOLUTIONs and KINETICS definitions are separate entities that can be reacted. 

You can react one SOLUTION (or one MIX) with one KINETICS data block (and other reactants like EQUILIBRIUM_PHASES). You can SAVE the resulting solution to react in some other way. The kinetic reaction moles will also change when reacted with a SOLUTION and are automatically saved.

So you need to develop a sequence of reaction calculations. Each calculation must have a SOLUTION, a MIX, or a USE solution definition. You may have at most one of each reactant data blocks (one KINETICS, and or one EQUILIBRIUM_PHASES, and or one EXCHANGE, etc), and the calculation definition ends with an END statement. If you want to use the solution composition produced by a calculation in a subsequent calculation, you can SAVE solution n, and then USE solution n in the subsequent calculation.

By default, a USER_GRAPH definition, once defined, will graph results for the current (the one it is defined in) and each subsequent calculation.

Perhaps you can give me a narrative of what you want to do. I want solution 1 to react with these minerals kinetically, then I want to take that solution and mix it with solution 2, ...

[helenhuntpvt@gmail.com]

Thank you for your response. I want to mix three solutions one containing K, Na, Ca with a pH of 13.1 second solution containing Ca, Cl, K, Na with a pH of 12.3 and another solution containing Ca with a pH of 10.4 as I have done in the phreeqc code presented earlier. Now, I want these three mixed solution to react/interact with quartz, kaolinite, k-feldspar, muscovite and calcite as I have provided their kinetics and rate to aid this interaction. My end goal is to generate a graph of pH/time of these interactions and a graph of possible element that will be released from this mineral when they react with the mixed solution but the output I am having after including the pH prompt isn't showing how this minerals interact with different pH.

[dlparkhurst]

Maybe this is what you had in mind?

Code: [Select]

RATES
 Quartz
 # d qu / dt = (A0 / V) * (m / m0)^0.67 * k * (1 - SI_qu)
 # Specific rate k from Rimstidt and Barnes, 1980, GCA 44, 1683
 #  k = 10^-13.7 mol/m2/s (25 C)
 #  A0, initial surface of quartz (m2) recalc's to mol/s
 #  V, solution volume in contact with A0
   -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Quartz") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Quartz") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0                 
    70 k_neut = 1.02e-14*EXP((-87700/8.3145)*(1/TK-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-SR("Quartz"))
    190 moles = r * TIME
    200 SAVE moles
     -end

K-Feldspar
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("K-Feldspar") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("K-Feldspar") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 8.71e-11*EXP((-51700/8.3145)*(1/Tk-1/298.15))                 
    70 k_neut = 3.89e-13*EXP((-38000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 6.31e-22*EXP((-94100/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.5 + k_neut + k_base*act("H+")^(-0.823)
    100 r = k_rateconst * SA * (1-SR("K-Feldspar"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Muscovite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Muscovite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Muscovite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 1.41e-12*EXP((-22000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-14*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.82e-15*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.37 + k_neut + k_base*act("H+")^(-0.22)
    100 r = k_rateconst * SA * (1-SR("Muscovite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Kaolinite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Kaolinite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Kaolinite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 7.7e-14*EXP((-65900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 7.7e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
    80 k_base = 4.70e-18*EXP((-17900/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.777 + k_neut + k_base*act("H+")^(-0.472)
    100 r = k_rateconst * SA * (1-SR("Kaolinite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Albite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Albite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Albite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 6.92e-11*EXP((-65000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.75e-13*EXP((-69800/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.51e-16*EXP((-71000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.457 + k_neut + k_base*act("H+")^(-0.572)
    100 r = k_rateconst * SA * (1-SR("Albite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Pyrite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Pyrite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Pyrite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 3.02e-8*EXP((-56900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-5*EXP((-56900/8.3145)*(1/Tk-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid*(act("H+")^(-0.5))*(act("Fe+2")^(0.5))  + k_neut*act("O2")^0.5  + k_base*act("H+")^0
    100 r = (k_rateconst) * SA * (1-SR("Pyrite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Calcite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Calcite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Calcite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0.5012*EXP(-14.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")                 
    70 k_neut = 1.5488e-6*EXP(-23.5e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 3.31e-4*EXP(-35.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("CO2")
    90 k_rateconst = k_acid*act("H+")^1  + k_neut + k_base*act("H+")^1
    100 r = k_rateconst * SA * (1-SR("Calcite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
END

SOLUTION 2 Young Cement Leachate (YCL)
    temp      25
    pH        13.1
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    K         0.092 mol/kgw
    Na        0.095 mol/kgw
    Ca        0.000135 mol/kgw

    -water    1 # kg

End
SOLUTION 3 Intermediate Leachate (IL)
    temp      25
    pH        12.3
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Ca        0.0162 mol/kgw
    Cl        0.004141 mol/kgw
    K         0.00397 mol/kgw
    Na        0.000171 mol/kgw
    -water    1 # kg

End
SOLUTION 4 Old Leachate (OL)
    temp      25
    pH        10.4
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Ca        0.000202 mol/kgw
    -water    1 # kg

End

Mix 1
    2  1
    3  1
    4  1

SAVE solution 5
END
USE solution 5

INCREMENTAL_REACTIONS true
KINETICS 1
Quartz
    -formula  SiO2  1
    -m        158.3
    -m0       158.3
    -parms    4.122
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        21.52
    -m0       21.52
    -parms    299.027
    -tol      1e-08
K-Feldspar
    -formula  KAlSi3O8  1
    -m        20.32
    -m0       20.32
    -parms    19.768
    -tol      1e-08
Muscovite     
    -formula  KAl3Si3O10(OH)2  1
    -m        11.4
    -m0       11.4
    -parms    25.607
    -tol      1e-08
Calcite       
    -formula  CaCO3  1
    -m        27.248
    -m0       27.248
    -parms    6.715
    -tol      1e-08

-steps   100*30 10*300  #31.536 3153.6 #3153.6 31536 31536 31536 31536 31536 315360 315360 3153600 3153600 3153600 3153600 3153600 3153600 3153600 31536000 31536000 31536000 31536000 31536000 31536000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 103280400
-step_divide 1
-runge_kutta 3
-bad_step_max 4000
-cvode true
-cvode_steps 500
-cvode_order 5

USER_GRAPH 1
    -headings               time Quartz Kaolinite K-Feldspar Muscovite Calcite pH
    -axis_titles            "Time, seconds" "SI" "pH"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 graph_x TOTAL_TIME
20 GRAPH_Y SI("Quartz"), SI("Kaolinite"), SI("K-Feldspar"), SI("Muscovite"), SI("Calcite")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true

END


[helenhuntpvt@gmail.com]

Thank you so much exactly what I wanted. However after running the block of code I get this error message  "Muscovite can't be found" the rates and kinetics I used for muscovite was found from examples here. Also could you help check why new chemicals are not being generated from the minerals using the user graph 2. Thanks   

Code: [Select]

RATES
 Quartz
 # d qu / dt = (A0 / V) * (m / m0)^0.67 * k * (1 - SI_qu)
 # Specific rate k from Rimstidt and Barnes, 1980, GCA 44, 1683
 #  k = 10^-13.7 mol/m2/s (25 C)
 #  A0, initial surface of quartz (m2) recalc's to mol/s
 #  V, solution volume in contact with A0
   -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Quartz") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Quartz") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0                 
    70 k_neut = 1.02e-14*EXP((-87700/8.3145)*(1/TK-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-SR("Quartz"))
    190 moles = r * TIME
    200 SAVE moles
     -end

K-Feldspar
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("K-Feldspar") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("K-Feldspar") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 8.71e-11*EXP((-51700/8.3145)*(1/Tk-1/298.15))                 
    70 k_neut = 3.89e-13*EXP((-38000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 6.31e-22*EXP((-94100/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.5 + k_neut + k_base*act("H+")^(-0.823)
    100 r = k_rateconst * SA * (1-SR("K-Feldspar"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Muscovite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Muscovite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Muscovite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 1.41e-12*EXP((-22000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-14*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.82e-15*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.37 + k_neut + k_base*act("H+")^(-0.22)
    100 r = k_rateconst * SA * (1-SR("Muscovite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Kaolinite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Kaolinite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Kaolinite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 7.7e-14*EXP((-65900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 7.7e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
    80 k_base = 4.70e-18*EXP((-17900/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.777 + k_neut + k_base*act("H+")^(-0.472)
    100 r = k_rateconst * SA * (1-SR("Kaolinite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Albite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Albite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Albite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 6.92e-11*EXP((-65000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.75e-13*EXP((-69800/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.51e-16*EXP((-71000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.457 + k_neut + k_base*act("H+")^(-0.572)
    100 r = k_rateconst * SA * (1-SR("Albite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Pyrite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Pyrite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Pyrite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 3.02e-8*EXP((-56900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-5*EXP((-56900/8.3145)*(1/Tk-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid*(act("H+")^(-0.5))*(act("Fe+2")^(0.5))  + k_neut*act("O2")^0.5  + k_base*act("H+")^0
    100 r = (k_rateconst) * SA * (1-SR("Pyrite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Calcite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Calcite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Calcite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0.5012*EXP(-14.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")                 
    70 k_neut = 1.5488e-6*EXP(-23.5e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 3.31e-4*EXP(-35.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("CO2")
    90 k_rateconst = k_acid*act("H+")^1  + k_neut + k_base*act("H+")^1
    100 r = k_rateconst * SA * (1-SR("Calcite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
END

SOLUTION 2 Young Cement Leachate (YCL)
    temp      25
    pH        13.1
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    K         0.092 mol/kgw
    Na        0.095 mol/kgw
    Ca        0.000135 mol/kgw

    -water    1 # kg

End
SOLUTION 3 Intermediate Leachate (IL)
    temp      25
    pH        12.3
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Ca        0.0162 mol/kgw
    Cl        0.004141 mol/kgw
    K         0.00397 mol/kgw
    Na        0.000171 mol/kgw
    -water    1 # kg

End
SOLUTION 4 Old Leachate (OL)
    temp      25
    pH        10.4
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Ca        0.000202 mol/kgw
    -water    1 # kg

End

Mix 1
    2  1
    3  1
    4  1

SAVE solution 5
END
USE solution 5

INCREMENTAL_REACTIONS true
KINETICS 1
Quartz
    -formula  SiO2  1
    -m        158.3
    -m0       158.3
    -parms    4.122
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        21.52
    -m0       21.52
    -parms    299.027
    -tol      1e-08
K-Feldspar
    -formula  KAlSi3O8  1
    -m        20.32
    -m0       20.32
    -parms    19.768
    -tol      1e-08
Muscovite     
    -formula  KAl3Si3O10(OH)2  1
    -m        11.4
    -m0       11.4
    -parms    25.607
    -tol      1e-08
Calcite       
    -formula  CaCO3  1
    -m        27.248
    -m0       27.248
    -parms    6.715
    -tol      1e-08

-steps   15*1 5*1  #31.536 3153.6 #3153.6 31536 31536 31536 31536 31536 315360 315360 3153600 3153600 3153600 3153600 3153600 3153600 3153600 31536000 31536000 31536000 31536000 31536000 31536000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 103280400
-step_divide 1
-runge_kutta 5
-bad_step_max 4000
-cvode true
-cvode_steps 500
-cvode_order 5

USER_GRAPH 1
    -chart_title "Chemical interaction of cement leachate with soil minerals"
    -headings               time Quartz Kaolinite K-Feldspar Muscovite Calcite pH
    -axis_titles            "Time, seconds" "SI" "pH"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 graph_x TOTAL_TIME / 3.1536e7 # time in years on x-axis
20 GRAPH_Y SI("Quartz"), SI("Kaolinite"), SI("K-Feldspar"), SI("Muscovite"), SI("Calcite")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true

END


USER_GRAPH 2
 -chart_title "Key Chemicals"
 -axis_titles Years "mmol / L"
 -axis_scale x_axis 0 5
 -initial_solutions  true
 -start
 40 graph_x total_time / 3.1536e7  # time in years on x-axis
 50 graph_y tot("Si") * 1e3        # parameter on y-axis
 -end
END

[dlparkhurst]

Muscovite must be defined in a PHASES data block in the database that you use, or in the input file.

I gave you a working file, you will have to look at the input and output to formulate a more specific question.