water-rock-gas reaction computing time

[tongg]

Hi David,

I am simulating the water-rock-gas interaction in the deep subsurface. And in my simulation, I first let the fresh water get equilibrium with rocks (e.g., albite, k-feldspar, and quartz) and gases (CO2, CH4 and H2). After equilibrium, the fresh water becomes reduced, then letting reduced water get equilibrium with rocks in that formation. After that, the first fresh water mixes with reduced water, using kinetic to simulate how the minerals change over long time (like 100 years) under the mixing water. But for this simulation, it took me quite a long time to calculate like over 10 hours I still cannot get the results, it is just still running. So I would like to know if is there any wrong sets in the code or what parameters I need to change? Thanks for your help!!

Here is the code:

Code: [Select]

DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.7.3-15968\database\llnl.dat
# This model incorporates multiple steps to 1) mix fresh GW with natural gases to produce reducing fluid,
# 2) equilibrate that fluid with red bed ss to dissolve Fe, and 2) mix that Fe-rich reducing fluid with
# fresh water



SOLUTION 1
    temp      17.1
    pH        7.46
    pe        7.2
    redox     pe
    units     mmol/l
    density   1
    Alkalinity 3.18 meq/l
    Ca        2.03
    Cl        0.076
    Fe        3.44 ug/l
    K         0.07
    Mg        1.81
    Na        0.4
    S(6)      2.28
    -water    1 # kg

END


#--------------------------------------------------------------------------------------------------
# Simulation 2: This equilibrates fresh GW with typical minerals in sandstone.

USE solution 1

EQUILIBRIUM_PHASES 2
    Albite    0 10
    K-Feldspar 0 10
    Quartz    0 10

SAVE solution 2
END


#----------------------------------------------------------------------------------------------------
# Simulation 3: Dissolution of natural gases into the groundwater, leading to a reduing fluid. The pressures can
# be changed depending on how deep this dissolution may be occuring.

USE solution 2
EQUILIBRIUM_PHASES 1
    CH4(g)    -0.12 10
    CO2(g)    2.02 10
    H2(g)     -0.36 10

SAVE solution 3
END


#-----------------------------------------------------------------------------------------------------
# Simulation 4: Letting the reducing fluid equilibrate with red bed sandstone, dissolving Fe.

USE solution 3
EQUILIBRIUM_PHASES 3
    Hematite  0 10
    Albite    0 10
    K-feldspar 0 10
    Quartz    0 10

SAVE solution 4
END

#------------------------------------------------------------------------------------------------------
# Simulation 5: Mixing the Fe-bearing reduced water with a fresh water.

MIX 1
    2    0.5
    4    0.5

SAVE solution 5
END
USE solution 5
EQUILIBRIUM_PHASES 1
    Calcite   0 10
    Chalcedony 0 10
    Fe(OH)3   0 10
    Goethite  0 10
TITLE Define the primary minerals kinetic equations here, the primary minerals are quartz, kaolinite, k-feldspar, illite, albite, anorthite, siderite and hematite. For calcite, define it in equilibrium_phases.

RATES
    Quartz
-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
 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
 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
    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
 60 k_acid = 1.047e-11*EXP((-23600/8.3145)*(1/Tk-1/298.15))
 70 k_neut = 6.918e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
 80 k_base = 8.913e-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
    Illite
-start
 10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
 30 if (M <= 0 and SR("Illite") < 1) then goto 200
 40 SA = PARM(1) * M
 50 if (M = 0 and SR("Illite") > 1) then SA = 1e-05
 60 k_acid = 1.047e-11*EXP(-23600/8.314*(1.0/TK-1.0/298.15))
 70 k_neut = 1.66e-13*EXP(-35000/8.314*(1.0/TK-1.0/298.15))
 80 k_base = 3.02e-17*EXP(-58900/8.314*(1.0/TK-1.0/298.15))
 90 k_rateconst = k_acid*act("H+")^0.34 + k_neut + k_base*act("H+")^(-0.4)
100 r = k_rateconst * SA * (1-SR("Illite"))
190 moles = r * TIME
200 SAVE moles
-end
    Siderite
-start
 10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
 30 if (M <= 0 and SR("Siderite") < 1) then goto 200
 40 SA = PARM(1) * M
 50 if (M = 0 and SR("Siderite") > 1) then SA = 1e-05
 60 k_acid = 6.457e-4*EXP((-36100/8.3145)*(1/Tk-1/298.15))
 70 k_neut = 1.26e-9*EXP((-62760/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("Siderite"))
190 moles = r * TIME
200 SAVE moles
-end
    Hematite
-start
 10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
 30 if (M <= 0 and SR("Hematite") < 1) then goto 200
 40 SA = PARM(1) * M
 50 if (M = 0 and SR("Hematite") > 1) then SA = 1e-05
 60 k_acid = 4.074e-10*EXP((-66200/8.3145)*(1/Tk-1/298.15))
 70 k_neut = 2.512e-15*EXP((-66210/8.3145)*(1/Tk-1/298.15))
 80 k_base = 0
 90 k_rateconst = k_acid*(act("H+")^(0.5))  + k_neut + k_base*act("H+")^0
100 r = (k_rateconst) * SA * (1-SR("Hematite"))
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
 60 k_acid = 1.45*EXP((-58400/8.3145)*(1/Tk-1/298.15))
 70 k_neut = 4.97e-10*EXP((-57000/8.3145)*(1/Tk-1/298.15))
 80 k_base = 7.51e-1*EXP((-55500/8.3145)*(1/Tk-1/298.15))
 90 k_rateconst = k_acid*(act("H+")^(0.335))  + k_neut + k_base*act("OH-")^(0.317)
100 r = (k_rateconst) * SA * (1-SR("Albite"))
190 moles = r * TIME
200 SAVE moles
-end
    Anorthite
-start
 10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
 30 if (M <= 0 and SR("Illite") < 1) then goto 200
 40 SA = PARM(1) * M
 50 if (M = 0 and SR("Illite") > 1) then SA = 1e-05
 60 k_acid = 2.58e-1*EXP(-16601/8.314*(1.0/TK-1.0/298.15))
 70 k_neut = 1.00e-6*EXP(-17821/8.314*(1.0/TK-1.0/298.15))
 80 k_base = 1.00e-22*EXP(-18150/8.314*(1.0/TK-1.0/298.15))
 90 k_rateconst = k_acid*act("H+")^1.411+ k_neut + k_base*act("H+")^(-1.767)
100 r = k_rateconst * SA * (1-SR("Illite"))
190 moles = r * TIME
200 SAVE moles
-end
KINETICS 1
K-feldspar
    -formula  KAlSi3O8  1
    -m        2.87
    -m0       2.87
    -parms    87.87
    -tol      1e-006
Illite
    -formula  K0.6Mg0.25Al2.3Si3.5O10(OH)2  1
    -m        0.87
    -m0       0.87
    -parms    617.23
    -tol      1e-006
Quartz
    -formula  SiO2  1
    -m        138.15
    -m0       138.15
    -parms    166.03
    -tol      1e-006
Siderite
    -formula  FeCO3  1
    -m        0.65
    -m0       0.65
    -parms    13.55
    -tol      1e-006
Albite
    -formula  NaAlSi3O8  1
    -m        0.29
    -m0       0.29
    -parms    7.6
    -tol      1e-006
Hematite
    -formula  Fe2O3  1
    -m        0.44
    -m0       0.44
    -parms    0.7
    -tol      1e-006
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        1
    -m0       1
    -parms    22.45992
    -tol      1e-006
Anorthite
    -formula  CaAl2(SiO4)2  1
    -m        0.27
    -m0       0.27
    -parms    7.51
    -tol      1e-006
-steps       3153600000 in 100 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 5000
-cvode true
-cvode_steps 100
-cvode_order 5
SELECTED_OUTPUT 1
    -file                 D:\PhD\phreeqc modeling\calcite eq.sel
    -reset                false
    -pH                   true
    -totals               Fe(2)  Fe  Fe(3)
    -equilibrium_phases   Chalcedony  Hematite  Illite  Fe(OH)3
    -saturation_indices   Calcite  Chalcedony  Dolomite  Hematite
                          Illite  Kaolinite  K-feldspar  Muscovite
                          Quartz  Siderite
    -kinetic_reactants    Calcite  Illite  Kaolinite  K-feldspar
                          Quartz  Siderite  Hematite
USER_GRAPH 2
    -headings               Years Calcite Chalcedony Hematite Illite Kaolinite K-feldspar Quartz Siderite Goethite
    -axis_titles            "Time" "SI" ""
    -chart_title            "The SI changes of minerals"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/31536000
20 GRAPH_Y SI("Calcite"), SI("Chalcedony"), SI("Hematite"),  SI("Illite"), SI("Kaolinite"), SI("K-feldspar"), SI("Quartz"), SI("Siderite"), SI("Goethite")
  -end
    -active                 true
USER_GRAPH 5
    -headings               Years Chalcedony Goethite Fe(OH)3
    -axis_titles            "Years" "Mole" ""
    -chart_title            "Secondary minerals changes over the years"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/31536000
20 GRAPH_Y EQUI_DELTA("Chalcedony"), EQUI_DELTA("Goethite"), EQUI_DELTA("Fe(OH)3")
  -end
    -active                 true
USER_GRAPH 1
    -headings               Years Calcite Illite Kaolinite K-Feldspar Quartz Sidersite Albite Hematite Anorthite
    -axis_titles            "Years" "Moles" ""
    -chart_title            "Minerals dissolution and precipitation"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/31536000
20 GRAPH_Y EQUI_DELTA("Calcite"), KIN_DELTA("Illite"), KIN_DELTA("Kaolinite"), KIN_DELTA("K-Feldspar"), KIN_DELTA("Quartz"), KIN_DELTA("Siderite"), KIN_DELTA("Albite"), KIN_DELTA("Hematite"),  KIN_DELTA("Anorthite")
  -end
    -active                 true





[dlparkhurst]

A few comments.

(1) I'm not a big fan of llnl.dat. There a lot of unlikely valence states, and a few species that can cause problems. One species is Al13O4(OH)24+7. You might be able to tweak the convergence parameters to get it too converge, but in trying to converge, small changes in the activities of Al+3, OH- or H+, make big changes in molality and cause PHREEQC problems. If you use another database, or remove this species, I think you can get it to run. You can then consider how important the species is to your simulation.

Code: [Select]

SOLUTION_SPECIES
28.0000 H2O + 13.0000 Al+++  =  Al13O4(OH)24+7 +32.0000 H+
#        -llnl_gamma           6.0   
        #log_k           -98.73
      log_k           -300
-delta_H 0       # Not possible to calculate enthalpy of reaction Al13O4(OH)24+7
# Enthalpy of formation: -0 kcal/mol
(2) You should use the following when doing the KINETICS simulation; it avoids repeating the early integration multiple times.

Code: [Select]

INCREMENTAL_REACTIONS true

(3) I think you should reconsider whether you want to include chalcedony, Fe(OH)3, and Goethite in EQUILIBRIUM_PHASES with the KINETICS reactions. First, all of the Fe(OH)3 will immediately dissolve and reprecipitate as goethite because the is the more stable phase. But second, goethite will continually dissolve and precipitate kinetically to hematite, because hematite is more stable than goethite. Similarly, chalcedony will continually dissolve to precipitate kinetically as quartz. Maybe you want to model these transitions, but using 10 moles of each makes me wonder whether you have really thought through these reactions.

[tongg]

Hi Prof. Parkhurst,

Thanks for your kind reply! And your suggestions made the code work! And why I chose chalcedony, goethite, and Fe(OH)3 in equilibrium phase is because I thought they might be the secondary minerals after the water-rock-gas reactions, and want to see how the iron minerals change through evolution. The results showed that Fe(OH)3 dissolved with goethite precipitated at the same time which is consistent with what you said. According to your suggestions, it seems quartz can also be the secondary mineral instead of chalcedony, so I delete these three minerals in equilibrium phase. And after running the codes, the results showed that quartz precipitates at first few decades and reaches equilibrium at 50 years, is that reasonable? In my expectation, the quarts should slowly precipitate over long time.

[dlparkhurst]

It is not my place to say whether it is reasonable or not. It is your investigation and a question of whether you trust your rate equations.

You may want to skip the kinetics calculation and simply consider first the equilibrium end-point that you will approach; that is, just use EQUILIBRIUM_PHASES with sets of minerals with or without chalcedony, quartz, goethite, hematite.  Then you can go back and work on kinetics.

[tongg]

Thank you again. I am not sure if my understanding is right. Do you mean that I just define the mineral sets with EQUILIBRIUM_PHASES, and do not need to define them in KINETIC?

[dlparkhurst]

I am just saying that you can calculate the end point of your reactions if you use EQUILIBRIUM_PHASES for the minerals instead of KINETICS. You will get different end points if you include chalcedony instead of quartz, goethite instead of hematite, and so on. The equilibrium solution composition will be somewhat independent of the masses of minerals. My thought is that you can investigate the end points first to decide which one you like best, and then you can go back and work on the kinetic simulation to reach that point.

[tongg]

Hi Prof. Parkhurst,

Sorry to bother you again. When I tried to see different gas pressures, I set all gases to 0, I met the same problem again. I tried another database (e.g., phreeqc.dat), the code just kept running without results. Can you help with that? Thanks a lot! Below is my new code.

Code: [Select]

DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.7.3-15968\database\phreeqc.dat




SOLUTION 1
    temp      17.1
    pH        7.46
    pe        7.2
    redox     pe
    units     mmol/l
    density   1
    Alkalinity 3.18 meq/l
    Ca        2.03
    Cl        0.076
    Fe        3.44 ug/l
    K         0.07
    Mg        1.81
    Na        0.4
    S(6)      2.28
    -water    1 # kg

END


#--------------------------------------------------------------------------------------------------
# Simulation 2: This equilibrates GW with typical minerals in sandstone.

USE solution 1

EQUILIBRIUM_PHASES 2
    Albite    0 10
    K-Feldspar 0 10
    Quartz    0 10
    Hematite  0 10

SAVE solution 2
END


#----------------------------------------------------------------------------------------------------
# Simulation 3: Dissolution of natural gases into the groundwater, leading to a reduing fluid. The pressures can
# be changed depending on how deep this dissolution may be occuring.

USE solution 2
EQUILIBRIUM_PHASES 1
    CH4(g)    0 10
    CO2(g)    0 10
    H2(g)     0 10

SAVE solution 3
END


#-----------------------------------------------------------------------------------------------------
# Simulation 4: Letting the reducing fluid equilibrate with red bed sandstone, dissolving Fe.

USE solution 3
EQUILIBRIUM_PHASES 3
    Albite    0 10
    Hematite  0 10
    K-feldspar 0 10
    Quartz    0 10

SAVE solution 4
END

#------------------------------------------------------------------------------------------------------
# Simulation 5: Mixing the Fe-bearing reduced water with a fresh water.

MIX 1
    2    0.5
    4    0.5

SAVE solution 5
END
USE solution 5
TITLE Define the primary minerals kinetic equations here, the primary minerals are quartz, kaolinite, k-feldspar, illite, albite, anorthite, siderite and hematite.
RATES
    Quartz
-start
  1 rem unit should be mol,kgw-1 and second-1
  2 rem parm(1) is surface area in the unit of m2/kgw
  3 rem calculation of surface area can be found in the note
  4 rem M is current moles of minerals. M0 is the initial moles of minerals
  5 rem parm(2) is a correction factor
 10 SR_mineral=SR("Quartz")
 11 if (M<0) then goto 200
 12 if (M=0 and SR_mineral<1) then goto 200
 13 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
 20 if (SA<=0) then SA=1
 30 R=8.31451
 40 Rate1=0*EXP(-0/R/TK)*ACT("H+")^0
 50 Rate2=6.64*EXP(-74526/R/TK)
 60 Rate3=0*EXP(-0/R/TK)*ACT("H+")^0
 70 Rate=(Rate1+Rate2+Rate3)*(1-SR_mineral)*SA*parm(2)
100 moles= rate*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
 60 k_acid = 0.5012*EXP(-14400/8.314*(1.0/TK-1.0/298.15))*ACT("H+")
 70 k_neut = 1.5488e-6*EXP(-23500/8.314*(1.0/TK-1.0/298.15))
 80 k_base = 3.31e-4*EXP(-35400/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
    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
 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
    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
 60 k_acid = 1.047e-11*EXP((-23600/8.3145)*(1/Tk-1/298.15))
 70 k_neut = 6.918e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
 80 k_base = 8.913e-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
    Illite
-start
 10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
 30 if (M <= 0 and SR("Illite") < 1) then goto 200
 40 SA = PARM(1) * M
 50 if (M = 0 and SR("Illite") > 1) then SA = 1e-05
 60 k_acid = 1.047e-11*EXP(-23600/8.314*(1.0/TK-1.0/298.15))
 70 k_neut = 1.66e-13*EXP(-35000/8.314*(1.0/TK-1.0/298.15))
 80 k_base = 3.02e-17*EXP(-58900/8.314*(1.0/TK-1.0/298.15))
 90 k_rateconst = k_acid*act("H+")^0.34 + k_neut + k_base*act("H+")^(-0.4)
100 r = k_rateconst * SA * (1-SR("Illite"))
190 moles = r * TIME
200 SAVE moles
-end
    Siderite
-start
 10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
 30 if (M <= 0 and SR("Siderite") < 1) then goto 200
 40 SA = PARM(1) * M
 50 if (M = 0 and SR("Siderite") > 1) then SA = 1e-05
 60 k_acid = 6.457e-4*EXP((-36100/8.3145)*(1/Tk-1/298.15))
 70 k_neut = 1.26e-9*EXP((-62760/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("Siderite"))
190 moles = r * TIME
200 SAVE moles
-end
    Hematite
-start
 10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
 30 if (M <= 0 and SR("Hematite") < 1) then goto 200
 40 SA = PARM(1) * M
 50 if (M = 0 and SR("Hematite") > 1) then SA = 1e-05
 60 k_acid = 4.074e-10*EXP((-66200/8.3145)*(1/Tk-1/298.15))
 70 k_neut = 2.512e-15*EXP((-66210/8.3145)*(1/Tk-1/298.15))
 80 k_base = 0
 90 k_rateconst = k_acid*(act("H+")^(0.5))  + k_neut + k_base*act("H+")^0
100 r = (k_rateconst) * SA * (1-SR("Hematite"))
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
 60 k_acid = 1.45*EXP((-58400/8.3145)*(1/Tk-1/298.15))
 70 k_neut = 4.97e-10*EXP((-57000/8.3145)*(1/Tk-1/298.15))
 80 k_base = 7.51e-1*EXP((-55500/8.3145)*(1/Tk-1/298.15))
 90 k_rateconst = k_acid*(act("H+")^(0.335))  + k_neut + k_base*act("OH-")^(0.317)
100 r = (k_rateconst) * SA * (1-SR("Albite"))
190 moles = r * TIME
200 SAVE moles
-end
    Anorthite
-start
 10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
 30 if (M <= 0 and SR("Illite") < 1) then goto 200
 40 SA = PARM(1) * M
 50 if (M = 0 and SR("Illite") > 1) then SA = 1e-05
 60 k_acid = 2.58e-1*EXP(-16601/8.314*(1.0/TK-1.0/298.15))
 70 k_neut = 1.00e-6*EXP(-17821/8.314*(1.0/TK-1.0/298.15))
 80 k_base = 1.00e-22*EXP(-18150/8.314*(1.0/TK-1.0/298.15))
 90 k_rateconst = k_acid*act("H+")^1.411+ k_neut + k_base*act("H+")^(-1.767)
100 r = k_rateconst * SA * (1-SR("Illite"))
190 moles = r * TIME
200 SAVE moles
-end
    Goethite
-start
 1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters
11 a1=0
12 E1=0
13 n1=0
20 rem neutral solution parameters
21 a2=1.64E+07
22 E2=86500
30 rem base solution parameters
31 a3=0
32 E3=0
33 n2=0
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Goethite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  #acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2    #base rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end
KINETICS 1
Quartz
    -formula  SiO2  1
    -m        138.15
    -m0       138.15
    -parms    23.13 0.16
    -tol      1e-08
Hematite
    -formula  Fe2O3  1
    -m        0.44
    -m0       0.44
    -parms    0.7
    -tol      1e-08
Albite
    -formula  NaAlSi3O8  1
    -m        0.29
    -m0       0.29
    -parms    7.59858
    -tol      1e-08
Anorthite
    -formula  CaAl2(SiO4)2  1
    -m        0.27
    -m0       0.27
    -parms    7.51
    -tol      1e-08
Illite
    -formula  K0.6Mg0.25Al1.8Al0.5Si3.5O10(OH)2  1
    -m        0.87
    -m0       0.87
    -parms    617.23455
    -tol      1e-08
K-feldspar
    -formula  K-feldspar  1
    -m        2.87
    -m0       2.87
    -parms    87.868781
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        0
    -m0       0
    -parms    22.45992
    -tol      1e-08
Siderite
    -formula  FeCO3  1
    -m        0
    -m0       0
    -parms    13.55
    -tol      1e-08
-steps       315360000000 in 1000 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 100000
-cvode true
-cvode_steps 100
-cvode_order 5
INCREMENTAL_REACTIONS True
EQUILIBRIUM_PHASES 4
    Calcite   0 3
    Goethite  0 0
SELECTED_OUTPUT 1
    -file                 D:\PhD\phreeqc modeling\calcite eq.sel
    -reset                false
    -pH                   true
    -totals               Fe(2)  Fe  Fe(3)
    -equilibrium_phases   Calcite
    -saturation_indices   Calcite  Chalcedony  Dolomite  Hematite
                          Illite  Kaolinite  K-feldspar  Muscovite
                          Quartz  Siderite
    -kinetic_reactants    Calcite  Illite  Kaolinite  K-feldspar
                          Quartz  Siderite  Hematite
USER_GRAPH 2
    -headings               Years Calcite Hematite Illite Kaolinite K-feldspar Quartz Siderite Alibite Anorthite Goethite
    -axis_titles            "Time" "SI" ""
    -chart_title            "The SI changes of minerals"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/31536000
20 GRAPH_Y SI("Calcite"), SI("Hematite"),  SI("Illite"), SI("Kaolinite"), SI("K-feldspar"), SI("Quartz"), SI("Siderite")
30 GRAPH_Y SI("Albite"), SI("Anorthite")
40 GRAPH_Y SI("Goethite")
  -end
    -active                 true
USER_GRAPH 1
    -headings               Years Calcite Illite K-Feldspar Kaolinite Quartz Albite Anorthite Hematite Siderite Goethite
    -axis_titles            "Years" "Moles" ""
    -chart_title            "Minerals dissolution and precipitation"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/31536000
20 GRAPH_Y EQUI_DELTA("Calcite")
30 GRAPH_Y KIN_DELTA("Illite")
40 GRAPH_Y KIN_DELTA("K-Feldspar"), KIN_DELTA("Kaolinite")
50 GRAPH_Y KIN_DELTA("Quartz")
60 GRAPH_Y KIN_DELTA("Albite"), KIN_DELTA("Anorthite")
70 GRAPH_Y KIN_DELTA("Hematite")
80 GRAPH_Y KIN_DELTA("Siderite")
90 GRAPH_Y EQUI_DELTA("Goethite")
  -end
    -active                 true
USER_GRAPH 5
    -headings               Quartz Siderite Hematite K-feldspar Anorthite Albite Calcite Illite Kaolinite Goethite
    -axis_titles            "Time" "Mass fraction" ""
    -chart_title            "Mass fraction"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
 10 x = TOTAL_TIME/31536000
 20 massq = 60.09*KIN("Quartz")
 30 masssd = 115.85*KIN("Siderite")
 40 massh = 159.68*KIN("Hematite")
 50 masskf = 278.33*KIN("K-feldspar")
 60 massan = 277.41*KIN("Anorthite")
 70 massal = 262.02*KIN("Albite")
 80 massc = 100.09*EQUI("Calcite")
 90 massi = 438.6*KIN("Illite")
100 massk = 258.16*KIN("Kaolinite")
110 massg = 88.85*EQUI("Goethite")
120 masstot = massq + masssd + massh + masskf + massan + massal + massc + massi + massk + massg
130 mfractionq = massq / masstot
140 mfractions = masssd / masstot
150 mfractionh = massh/ masstot
160 mfractionkf = masskf / masstot
170 mfractionan = massan / masstot
180 mfractional = massal / masstot
190 mfractionc = massc / masstot
200 mfractioni = massi / masstot
210 mfractionk = massk / masstot
220 mfractiong = massg / masstot
230 PLOT_XY x, mfractionq, color = gray, symbol = None
240 PLOT_XY x, mfractions, color = black, symbol = None
250 PLOT_XY x, mfractionh, color = red, symbol = None
260 PLOT_XY x, mfractionkf, color = green, symbol = None
270 PLOT_XY x, mfractionan, color = cyan, symbol = None
280 PLOT_XY x, mfractional, color = magenta, symbol = None
290 PLOT_XY x, mfractionc, color = blue, symbol = None
300 PLOT_XY x, mfractioni, color = purple, symbol = None
310 PLOT_XY x, mfractionk, color = orange, symbol = None
320 PLOT_XY x, mfractionk, color = olive, symbol = None
  -end
    -active                 true





[dlparkhurst]

Illite equilibrates in less than 1 day. If you treat Illite as EQUILIBRIUM_PHASE instead of KINETICS, I think it might work.

[tongg]

Hi Prof. Parkhurst,

Thanks for your reply!! That does make the code work, but the result becomes really weird, how can I know which mineral is suitable to be put in KINETIC or EQUILIBRIUM_PHASE when the code doesn't work?

Thanks a lot!

Best,
Tong

[dlparkhurst]

Illite equilibrated in a few hours, so I assumed it's rate must be fast. Any mineral that you can assume reacts quickly on your time scale of interest is a candidate for EQUILIBRIUM_PHASES.

Don't know about weird results. I think it is weird to have extremely low concentrations of Mg. You'll have to consider how well you know the reactants and the rates.