ADVECTION combine with KINETIC

[tongg]

Hi Prof. Parkhurst,

We are trying to investigate the oxygenated water mixes with reduced water, making the mixed fluid reacts with minerals, and we tried to use KINETIC and ADVECTION blocks to see how the minerals would change over geologic time, but I cannot run out the results. Could you help to check where I set it wrong and what I can change to get the results?

And I have one place that confuses me which is when I only run the KINETIC block to see how the minerals changes, the results showed that the illite dissolving completely and k-feldspar precipitating which is not consistent with the mineral stability diagrams, which is if k-feldspar precipitates, illite should precipitate before k-feldspar. What we observe in the research area is also k-feldspar dissolves, can you help to explain what causes this phenomenon?
P.S. You mentioned that the Mg concentration is low in water sample, I tried to increase the Mg concentration, it seems that slow the illite dissolution rate, illite still dissolves completely when the simulation period is extended.

Below is my code:

Code: [Select]

DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.7.3-15968\database\phreeqc.dat
# This model incorporates multiple steps to 1) mix Navajo 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_MASTER_SPECIES
    Acetate       Acetate-         1     59.045          59.045
SOLUTION_SPECIES
Acetate- = Acetate-
    log_k     0
Acetate- + H+ = H(Acetate)
    log_k     4.757
    delta_h   0.41 kJ
    -gamma    0 0
Acetate- + Ca+2 = Ca(Acetate)+
    log_k     1.18
    delta_h   4 kJ
    -gamma    0 0
Acetate- + Na+ = Na(Acetate)
    log_k     -0.18
    delta_h   12 kJ
    -gamma    0 0
Acetate- + K+ = K(Acetate)
    log_k     -0.1955
    delta_h   4.184 kJ
    -gamma    0 0
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 Navajo GW with typical minerals in sandstone.

USE solution 1

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

SAVE solution 2
END


#----------------------------------------------------------------------------------------------------
# Simulation 3: Dissolution of natural gases into the groundwater, leading to a reduing fluid.

USE solution 2
EQUILIBRIUM_PHASES 1
    Albite    0 10
    CH4(g)    0 10
    CO2(g)    1.98 10
    H2(g)     0 10
    Hematite  0 10
    K-feldspar 0 10
    Quartz    0 10
    Siderite  0 10

SAVE solution 3
END



#------------------------------------------------------------------------------------------------------
# Simulation 5: Mixing the Fe-bearing reduced water with a fresh water.
MIX 1
    2    0.7
    3    0.3

SAVE solution 0
END
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
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
Siderite
    -formula  FeCO3  1
    -m        0.44
    -m0       0.44
    -parms    13.55
    -tol      1e-08
K-feldspar
    -formula  KAlSi3O8  1
    -m        2.87
    -m0       2.87
    -parms    87.87
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        0
    -m0       0
    -parms    22.45992
    -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
    Illite    0 0.87
ADVECTION
    -cells 1
    -shifts 5000
    -time_step 315360000000 # seconds
    -warnings true
SELECTED_OUTPUT 1
    -file                 D:\PhD\phreeqc modeling\KINETICS+ADVECTION.sel
    -reset                false
    -pH                   true
    -totals               Fe(2)  Fe  Fe(3)  K
    -equilibrium_phases   Calcite
    -saturation_indices   Calcite  Chalcedony  Dolomite  Hematite
                          Illite  Kaolinite  K-feldspar  Quartz  Siderite
    -kinetic_reactants    Calcite  Illite  Kaolinite  K-feldspar
                          Quartz  Siderite  Hematite
USER_GRAPH 1
    -headings               Quartz Siderite Hematite K-feldspar Anorthite Albite Calcite Illite Kaolinite Goethite
    -axis_titles            "Pore volumes" "Weight % mineral" ""
    -chart_title            "Bleaching of red sandstone by reduced hydrocarbon-bearing fluid"
    -axis_scale x_axis      0.1 5000 auto auto log
    -axis_scale y_axis      0 100 auto auto
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
 10 x = (STEP_NO+0.5)/CELL_NO
 20 weightq = 60.09*KIN("Quartz")/100
 30 weights = 115.85*KIN("Siderite")/100
 40 weighth = 159.68*KIN("Hematite")/100
 50 weightkf = 278.33*KIN("K-feldspar")/100
 60 weightan = 277.41*KIN("Anorthite")/100
 70 weightal = 262.02*KIN("Albite")/100
 80 weightc = 100.09*EQUI("Calcite")/100
 90 weighti = 438.6*EQUI("Illite")/100
100 weightkl = 258.16*KIN("Kaolinite")/100
110 weightg = 88.85*EQUI("Goethite")/100
120 PLOT_XY x, weightq, color = gray, symbol = None
130 PLOT_XY x, weights, color = black, symbol = None
140 PLOT_XY x, weighth, color = red, symbol = None
150 PLOT_XY x, weightkf, color = green, symbol = None
160 PLOT_XY x, weightan, color = cyan, symbol = None
170 PLOT_XY x, weightal, color = magenta, symbol = None
180 PLOT_XY x, weightc, color = blue, symbol = None
190 PLOT_XY x, weighti, color = lime, symbol = None
200 PLOT_XY x, weightkl, color = olive, symbol = None
210 PLOT_XY x, weightkg, color = pink, symbol = None
220 PUT(1, 1)
  -end
    -active                 true
USER_GRAPH 2
    -headings               Fe Si K S Acetate pH
    -axis_titles            "Pore volumes" "Millimoles per kilogram water" "pH"
    -chart_title            "Water composition & pH change"
    -axis_scale x_axis      1 1000 auto auto log
    -axis_scale y_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 x = (STEP_NO+0.5)/CELL_NO
20 PLOT_XY x, TOT("Fe(+2)")*1000, color = red, symbol = None, y_axis = 1
30 PLOT_XY x, TOT("Si")*1000, color = magenta, symbol = None, y_axis = 1
40 PLOT_XY x, TOT("K")*1000, color = green, symbol = None, y_axis = 1
50 PLOT_XY x, TOT("S(-2)")*1000, color = orange, symbol = None, y_axis = 1
60 PLOT_XY x, TOT("C")*1000, color = yellow, symbol = None, y_axis = 1
70 PLOT_XY x, -LA("H+"), color = black, symbol = None, y_axis = 2
80 PUT(1, 1)
  -end
    -active                 true




Thank you in advance!

Best,
Tong

[dlparkhurst]

Let's step back a bit and make sure you are modeling what you think you are modeling.

First, look at the results of your reduction step. You are specifying that you want the reacted solution to be at 1 atm CH4(g), 1 atm H2(g), and 100 atm of CO2(g) plus equilibrium with multiple other phases. There are 10 moles of each reactant in EQUILIBRIUM_PHASES. Here is the result of the reaction.

Code: [Select]

-------------------------------Phase assemblage--------------------------------

                                                      Moles in assemblage
Phase               SI  log IAP  log K(T, P)   Initial       Final       Delta

Albite            0.00   -18.52    -18.52    1.000e+01   1.000e+01   3.450e-04
CH4(g)           -0.00    -2.72     -2.72    1.000e+01   1.143e+01   1.425e+00
CO2(g)           -0.94    -2.31     -1.37    1.000e+01           0  -1.000e+01
H2(g)            -6.51    -9.59     -3.08    1.000e+01           0  -1.000e+01
Hematite         -0.00    -3.39     -3.39    1.000e+01   5.711e+00  -4.289e+00
K-feldspar        0.00   -21.19    -21.19    1.000e+01   1.000e+01   7.405e-07
Quartz           -0.00    -4.10     -4.10    1.000e+01   9.999e+00  -1.048e-03
Siderite          0.00   -10.84    -10.84    1.000e+01   1.858e+01   8.578e+00
All of the H2 and CO2 react so that there is none remaining in equilibrium_phases. P(CO2) after reaction is 10^-0.94 atm, P(H2) is 10^-6.5, and P(CH4) is 1 atm. You have generated methane, and a lot of hematite has been converted to siderite. The solution is a reasonable solution with 1 atm methane, so methane will be the reductant for the bleaching. Maybe that is what you want, but the reactions to make it seem arbitrary.

Perhaps ultimately you want to use Acetate for the reductant, but Acetate is not used in this simulation, so the current simulation is simpler without it.

I guess my biggest problem is related to ADVECTION. This simulation will use EQUILIBRIUM_PHASES 1 in the one cell that is defined, and EQUILIBRIUM_PHASES 4 is not used in the simulation at all (it would be used in cell 4 if ADVECTION had 4 cells). So, you have the same minerals in EQUILIBRIUM_PHASES 1 and KINETICS 1. I'm guessing that you want to change the user number of EQUILIBRIUM_PHASES 4 to EQUILIBRIUM_PHASES 1.

[tongg]

Hi Prof. Parkhurst,

I really appreciate your response so quickly, I changed the EQUILIBRIUM_PHASE4 to EQUILIBRIUM_PHASE1, the code worked, but I got the results after more than 11 hours running...

What we found is that there is a gas reservoir in subsurface, what I think is the minerals of sandstone should equilibrium with these gases too so I put them together in EQUILIBRIUM_PHASE, and also let the fresh water gets equilibrium with them. I know I should set the minerals with specific moles of that area, but when I set the values and ran the codes only with KINETIC, the results showed the changes of minerals precipitation/dissolution amounts instability at the last period.

As you mentioned the produce of methane, we noticed it, so I am wondering if we should consider H2 in our model, and about the unreasonable(?) dissolution of illite and precipitation of k-feldspar, is that related to the produce of methane?

Thank you for answering me these questions!

Best,
Tong 

[dlparkhurst]

I think it is reasonable to specify partial pressures of gases that are consistent with your gas reservoir.

I really can't comment on the rest of your reaction network. If possible, I again think you should work with equilibrium rather than kinetics until you get reactions that you think are reasonable. You can avoid issues with rates of reactions and surface area. If you are not happy with the relative stabilities of K-spar and illite, you may need to consider whether you have consistent log Ks and look at other possible equilibrium constants and other clay minerals.

[tongg]

Hi Prof. Parkhurst,

I got confused when I ran these KINETIC and ADVECTION combined with KINETIC models. In the KINETIC model, I can get the precipitation of goethite and other minerals dissolve/precipitate, but when I run ADVECTION combined with KINETIC (i.e., based on KINETIC, I added ADVECTION), I can't get the precipitation of goethite, only quartz, kaolinite, k-feldspar and calcite have changes, other minerals showed 0 with increasing shifts. Can you help explain why I got such different results? Below is my code.

Thank you in advance.

Best,
Tong

Code: [Select]

# This model incorporates multiple steps to 1) mix Navajo 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

# In solution1 I changed the concentration of Fe and S(-2), either of them increased can get more goethite precipitated but only a little.
SOLUTION 1
    temp      14.3
    pH        5
    pe        7.2
    redox     pe
    units     mmol/l
    density   1
    Acetate   50 mg/l
    Alkalinity 1.52 meq/l
    Ca        0.12
    Cl        0.51
    Fe        0.02
    K         0.08
    Mg        0.44
    Na        0.92
    O(0)      10 mg/l
    S(-2)     10 mg/l
    S(6)      0.64
    Si        0.32
    -water    1 # kg

SAVE solution 1
END


#--------------------------------------------------------------------------------------------------
# Simulation 2: I set all the minerals below as their specific moles I think might be more close to the real situation.

USE solution 1

EQUILIBRIUM_PHASES 2
    Calcite   0 3
    Illite    0 0.87
    K-Feldspar 0 2.87
    Quartz    0 138.15
    Hematite  0 0.44
    Kaolinite    0 0.85

SAVE solution 2
END


#----------------------------------------------------------------------------------------------------
# Simulation 3:Lower the pressure of CO2 and increased CH4 (help dissolve more hematite)

USE solution 2
EQUILIBRIUM_PHASES 3
    CH4(g)    0 10
    CO2(g)    1.5 100
    Calcite   0 3
    Hematite  0 0.44
    Illite    0 0.87
    K-feldspar 0 2.87
    Quartz    0 138.15
    Kaolinite     0 0.85
SAVE solution 3
END



#------------------------------------------------------------------------------------------------------
# Simulation 4: Feel free to change the proportion of these two fluids!

MIX 1
    1    0.6
    3    0.4

SAVE solution 0
END
USE solution 0
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
    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
Quartz
    -formula  SiO2  1
    -m        138.15
    -m0       138.15
    -parms    23.13 0.16
    -tol      1e-008
K-feldspar
    -formula  KAlSi3O8  1
    -m        2.87
    -m0       2.87
    -parms    87.87
    -tol      1e-008
Illite
    -formula  K0.6Mg0.25Al1.8Al0.5Si3.5O10(OH)2  1
    -m        0.87
    -m0       0.87
    -parms    617.23
    -tol      1e-008
Albite
    -formula  NaAlSi3O8  1
    -m        0.29
    -m0       0.29
    -parms    7.6
    -tol      1e-008
Anorthite
    -formula  CaAl2Si2O8  1
    -m        0.27
    -m0       0.27
    -parms    7.5
    -tol      1e-008
-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 1
    Calcite   0 3
    Goethite  0 0
    Kaolinite 0 0.75
ADVECTION
    -cells 1
    -shifts 5000
    -time_step 62712000 # seconds
    -warnings true
SELECTED_OUTPUT 1
    -file                 D:\PhD\phreeqc modeling\Coleman's code combine.sel
    -reset                false
    -pH                   true
    -totals               Fe(2)  Fe  Fe(3)  K
    -equilibrium_phases   Calcite
    -saturation_indices   Calcite  Chalcedony  Dolomite  Hematite
                          Illite  Kaolinite  K-feldspar  Quartz  Siderite
    -kinetic_reactants    Calcite  Illite  Kaolinite  K-feldspar
                          Quartz  Siderite  Hematite
USER_GRAPH 1
    -headings               Quartz Goethite Kaolinite K-feldspar Anorthite Albite Calcite Illite
    -axis_titles            "Pore volumes" "Weight % mineral" ""
    -chart_title            "Bleaching of red sandstone by reduced hydrocarbon-bearing fluid"
    -axis_scale x_axis      0.1 5000 auto auto log
    -axis_scale y_axis      0 100 auto auto
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
 10 x = (STEP_NO+0.5)/CELL_NO
 20 weightq = 60.09*KIN("Quartz")/100
 30 weightg = 88.85*EQUI("Goethite")/100
 40 weightk = 258.16*EQUI("Kaolinite")/100
 50 weightkf = 278.33*KIN("K-feldspar")/100
 60 weightan = 277.41*KIN("Anorthite")/100
 70 weightal = 262.02*KIN("Albite")/100
 80 weightc = 100.09*EQUI("Calcite")/100
 90 weighti = 438.6*KIN("Illite")/100
100 PLOT_XY x, weightq, color = gray, symbol = None
110 PLOT_XY x, weightg, color = brown, symbol = None
120 PLOT_XY x, weightk, color = orange, symbol = None
130 PLOT_XY x, weightkf, color = green, symbol = None
140 PLOT_XY x, weightan, color = cyan, symbol = None
150 PLOT_XY x, weightal, color = magenta, symbol = None
160 PLOT_XY x, weightc, color = blue, symbol = None
170 PLOT_XY x, weighti, color = lime, symbol = None
180 PUT(1, 1)
  -end
    -active                 true





[dlparkhurst]

Do you still have an issue if your replace "Illite" with "Anorthite" in the Anorthite RATES definition?

[tongg]

Hi Prof. Parkhurst,

Thank you for pointing out that issue! Yes, the result doesn't have much difference, I still can't get goethite precipitation.

Best,
Tong 

[dlparkhurst]

Your solution 0 has a lot of Fe, but it is all Fe(2). Goethite is undersaturated in solution 0. There are no oxidizing reactions in the kinetic reactions in KINETICS 1, and no oxidants in equilibrium phases, and further, no sources or sinks of Fe other than goethite. So, the solution simply remains with a lot of Fe(2), but undersaturated with goethite.

What is the point of running 5000 steps? The same reaction will occur at each step.

[tongg]

Hi Prof. Parkhurst

Thanks for your reply. I can get the precipitation of goethite after adding an oxidant in EQUILIBRIUM. The reason we set 5000 shifts is because we want to see how the mineral assemblages change with the fluid advance. So the time I set step in ADVECTION is kinetic time/shifts, but the result figure is really weird, like there are lots of lines connecting the initial point to the different pore volume number... where did I set up wrong?

I really appreciate your help.

Best,
Tong

[dlparkhurst]

I need to see your revised file.

[tongg]

Hi Dr. Parkhurst,

Below is my script.

Code: [Select]

# This model incorporates multiple steps to 1) mix Navajo 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

# In solution1 I changed the concentration of Fe and S(-2), either of them increased can get more goethite precipitated but only a little.
SOLUTION 1
    temp      14.3
    pH        5
    pe        7.2
    redox     pe
    units     mmol/l
    density   1
    Acetate   50 mg/l
    Alkalinity 1.52 meq/l
    Ca        0.12
    Cl        0.51
    Fe        0.02
    K         0.08
    Mg        0.44
    Na        0.92
    O(0)      60 mg/l
    S(6)      0.64
    Si        0.32
    -water    1 # kg

SAVE solution 1
END


#--------------------------------------------------------------------------------------------------
# Simulation 2: I set all the minerals below as their specific moles I think might be more close to the real situation.

USE solution 1

EQUILIBRIUM_PHASES 2
    Calcite   0 3
    Illite    0 0.87
    K-Feldspar 0 2.87
    Quartz    0 138.15
    Hematite  0 0.44
    Kaolinite    0 0.85

SAVE solution 2
END


#----------------------------------------------------------------------------------------------------
# Simulation 3:Lower the pressure of CO2 and increased CH4 (help dissolve more hematite)

USE solution 2
EQUILIBRIUM_PHASES 3
    CH4(g)    0 10
    CO2(g)    1.5 100
    Calcite   0 3
    Hematite  0 0.44
    Illite    0 0.87
    K-feldspar 0 2.87
    Quartz    0 138.15
    Kaolinite     0 0.85
SAVE solution 3
END



#------------------------------------------------------------------------------------------------------
# Simulation 4: Feel free to change the proportion of these two fluids!

MIX 1
    1    0.6
    3    0.4

SAVE solution 0
END
USE solution 0
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
    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("Anorthite") < 1) then goto 200
 40 SA = PARM(1) * M
 50 if (M = 0 and SR("Anorthite") > 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("Anorthite"))
190 moles = r * TIME
200 SAVE moles
-end

KINETICS 1
Quartz
    -formula  SiO2  1
    -m        138.15
    -m0       138.15
    -parms    23.13 0.16
    -tol      1e-008
K-feldspar
    -formula  KAlSi3O8  1
    -m        2.87
    -m0       2.87
    -parms    87.87
    -tol      1e-008
Illite
    -formula  K0.6Mg0.25Al1.8Al0.5Si3.5O10(OH)2  1
    -m        0.87
    -m0       0.87
    -parms    617.23
    -tol      1e-008
Albite
    -formula  NaAlSi3O8  1
    -m        0.29
    -m0       0.29
    -parms    7.6
    -tol      1e-008
Anorthite
    -formula  CaAl2Si2O8  1
    -m        0.27
    -m0       0.27
    -parms    7.5
    -tol      1e-008
-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 1
    Calcite   0 3
    Goethite  0 0
    Kaolinite 0 0.75
    O2(g)     -3 10
ADVECTION
    -cells 1
    -shifts 5000
    -time_step 63072000 # seconds
    -warnings true
SELECTED_OUTPUT 1
    -file                 D:\PhD\phreeqc modeling\Coleman's code combine.sel
    -reset                false
    -pH                   true
    -totals               Fe(2)  Fe  Fe(3)  K
    -equilibrium_phases   Calcite
    -saturation_indices   Calcite  Chalcedony  Dolomite  Hematite
                          Illite  Kaolinite  K-feldspar  Quartz  Siderite
    -kinetic_reactants    Calcite  Illite  Kaolinite  K-feldspar
                          Quartz  Siderite  Hematite
USER_GRAPH 1
    -headings               Quartz Goethite Kaolinite K-feldspar Anorthite Albite Calcite Illite
    -axis_titles            "Pore volumes" "Weight % mineral" ""
    -chart_title            "Bleaching of red sandstone by reduced hydrocarbon-bearing fluid"
    -axis_scale x_axis      0.1 5000 auto auto log
    -axis_scale y_axis      0 100 auto auto
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
 10 x = (STEP_NO+0.5)/CELL_NO
 20 weightq = 60.09*KIN("Quartz")/100
 30 weightg = 88.85*EQUI("Goethite")/100
 40 weightk = 258.16*EQUI("Kaolinite")/100
 50 weightkf = 278.33*KIN("K-feldspar")/100
 60 weightan = 277.41*KIN("Anorthite")/100
 70 weightal = 262.02*KIN("Albite")/100
 80 weightc = 100.09*EQUI("Calcite")/100
 90 weighti = 438.6*KIN("Illite")/100
100 PLOT_XY x, weightq, color = gray, symbol = None
110 PLOT_XY x, weightg, color = brown, symbol = None
120 PLOT_XY x, weightk, color = orange, symbol = None
130 PLOT_XY x, weightkf, color = green, symbol = None
140 PLOT_XY x, weightan, color = cyan, symbol = None
150 PLOT_XY x, weightal, color = magenta, symbol = None
160 PLOT_XY x, weightc, color = blue, symbol = None
170 PLOT_XY x, weighti, color = lime, symbol = None
180 PUT(1, 1)
  -end
    -active                 true





[dlparkhurst]

The script below is reorganized to show the calculations you are making more clearly.

You are running a kinetic simulation with solution 0, equilibrium_phases 1, and kinetics 1 for 1000 time steps.

Then you are running advection for 5000 time steps. The curves connect from the end of the kinetics run to the advection run.

By advecting solution 0 into cell 1, the kinetic reaction that occurs in cell 1 is almost identical at each shift. Until the change in M/M0 is significant, the change in porosity is linear in time. You have most of the information that you need after 1 advection step.

Code: [Select]

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
    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("Anorthite") < 1) then goto 200
 40 SA = PARM(1) * M
 50 if (M = 0 and SR("Anorthite") > 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("Anorthite"))
190 moles = r * TIME
200 SAVE moles
-end
END
# This model incorporates multiple steps to 1) mix Navajo 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

# In solution1 I changed the concentration of Fe and S(-2), either of them increased can get more goethite precipitated but only a little.
SOLUTION 1
    temp      14.3
    pH        5
    pe        7.2
    redox     pe
    units     mmol/l
    density   1
    Acetate   50 mg/l
    Alkalinity 1.52 meq/l
    Ca        0.12
    Cl        0.51
    Fe        0.02
    K         0.08
    Mg        0.44
    Na        0.92
    O(0)      60 mg/l
    S(6)      0.64
    Si        0.32
    -water    1 # kg
END


#--------------------------------------------------------------------------------------------------
# Simulation 2: I set all the minerals below as their specific moles I think might be more close to the real situation.

USE solution 1

EQUILIBRIUM_PHASES 2
    Calcite   0 3
    Illite    0 0.87
    K-Feldspar 0 2.87
    Quartz    0 138.15
    Hematite  0 0.44
    Kaolinite    0 0.85

SAVE solution 2
END


#----------------------------------------------------------------------------------------------------
# Simulation 3:Lower the pressure of CO2 and increased CH4 (help dissolve more hematite)

USE solution 2
EQUILIBRIUM_PHASES 3
    CH4(g)    0 10
    CO2(g)    1.5 100
    Calcite   0 3
    Hematite  0 0.44
    Illite    0 0.87
    K-feldspar 0 2.87
    Quartz    0 138.15
    Kaolinite     0 0.85
SAVE solution 3
END



#------------------------------------------------------------------------------------------------------
# Simulation 4: Feel free to change the proportion of these two fluids!

MIX 1
    1    0.6
    3    0.4

SAVE solution 0
END
EQUILIBRIUM_PHASES 1
    Calcite   0 3
    Goethite  0 0
    Kaolinite 0 0.75
    O2(g)     -3 10
END
KINETICS 1
Quartz
    -formula  SiO2  1
    -m        138.15
    -m0       138.15
    -parms    23.13 0.16
    -tol      1e-008
K-feldspar
    -formula  KAlSi3O8  1
    -m        2.87
    -m0       2.87
    -parms    87.87
    -tol      1e-008
Illite
    -formula  K0.6Mg0.25Al1.8Al0.5Si3.5O10(OH)2  1
    -m        0.87
    -m0       0.87
    -parms    617.23
    -tol      1e-008
Albite
    -formula  NaAlSi3O8  1
    -m        0.29
    -m0       0.29
    -parms    7.6
    -tol      1e-008
Anorthite
    -formula  CaAl2Si2O8  1
    -m        0.27
    -m0       0.27
    -parms    7.5
    -tol      1e-008
-steps       315360000000 in 1000 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 100000
-cvode true
-cvode_steps 100
-cvode_order 5
END
INCREMENTAL_REACTIONS True
END

TITLE Define the primary minerals kinetic equations here, the primary minerals are quartz, kaolinite, k-feldspar, illite, albite, anorthite, siderite and hematite.
USE solution 0
USE equilibrium_phases 1
USE kinetics 1
SELECTED_OUTPUT 1
    #-file                 D:\PhD\phreeqc modeling\Coleman's code combine.sel
    -reset                false
    -pH                   true
    -totals               Fe(2)  Fe  Fe(3)  K
    -equilibrium_phases   Calcite
    -saturation_indices   Calcite  Chalcedony  Dolomite  Hematite
                          Illite  Kaolinite  K-feldspar  Quartz  Siderite
    -kinetic_reactants    Calcite  Illite  Kaolinite  K-feldspar
                          Quartz  Siderite  Hematite
USER_GRAPH 1
    -headings               Quartz Goethite Kaolinite K-feldspar Anorthite Albite Calcite Illite
    -axis_titles            "Time, seconds" "Weight % mineral" ""
    -chart_title            "Bleaching of red sandstone by reduced hydrocarbon-bearing fluid"
    #-axis_scale x_axis      auto auto auto auto log
    #-axis_scale y_axis      auto auto auto auto
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
 10 x = total_time #(STEP_NO+0.5)/CELL_NO
 20 weightq = 60.09*KIN("Quartz")/100
 30 weightg = 88.85*EQUI("Goethite")/100
 40 weightk = 258.16*EQUI("Kaolinite")/100
 50 weightkf = 278.33*KIN("K-feldspar")/100
 60 weightan = 277.41*KIN("Anorthite")/100
 70 weightal = 262.02*KIN("Albite")/100
 80 weightc = 100.09*EQUI("Calcite")/100
 90 weighti = 438.6*KIN("Illite")/100
100 PLOT_XY x, weightq, color = gray, symbol = None
110 PLOT_XY x, weightg, color = brown, symbol = None
120 PLOT_XY x, weightk, color = orange, symbol = None
130 PLOT_XY x, weightkf, color = green, symbol = None
140 PLOT_XY x, weightan, color = cyan, symbol = None
150 PLOT_XY x, weightal, color = magenta, symbol = None
160 PLOT_XY x, weightc, color = blue, symbol = None
170 PLOT_XY x, weighti, color = lime, symbol = None
180 PUT(1, 1)
  -end
    -active                 true
END

ADVECTION
    -cells 1
    -shifts 5000
    -time_step 63072000 # seconds
    -warnings true
END


[tongg]

Hi Dr. Parkhurst,

Thanks for your reply, the figure is much more clear!! But I am wondering if I can let the kinetic happen during the advection? So I can have the figure to show the mineral weight changes with pore volumes/distance? If it can, what settings should I change?

Many thanks for your help.

Best regards,
Tong

[dlparkhurst]

The kinetic reaction is occurring in cell 1 during ADVECTION.

You have only one cell, so there is only one distance. If you have more cells, you can the X variable can be CELL_NO, which would correlate with distance. With TRANSPORT, you can specify cell lengths and have a real distance axis.

I leave it as an exercise to you to develop the plot that you want.