Water Rock reaction with Supercritical CO2

[zscowcrofst]

I trying to simulate the interaction between a sandstone and brine at 75?C, followed by exposure to supercritical CO₂ at 85?C and 250 atm pressure. The sandstone consists of 65% Quartz, 20% K-feldspar, and 8% Forsterite (volume fraction). I assume I should first equilibrate the brine with the minerals (Solution 1) and then introduce supercritical CO₂ (Solution 2). The simulation should run for 200 years, tracking mineral volume fractions, saturation indices, and major ion concentration changes over time. The following is my input file. I used the rates in phreeqc and llnl database. I am not sure what the problem is but it looks like rates are not right. The error appears at Quartz though.

Code: [Select]

TITLE Sandstone Equilibrium with Brine and Supercritical CO2 Injection

SOLUTION 1  Brine in equilibrium with minerals
    temp 75
    pH 7
    pe 4
    units mg/L
    Alkalinity 149
    B 16.1
    Ca 2704
    Cl 88409
    K 168
    Mg 889
    Mn 0.81
    Na 49706
    S(6) 10.02
    Si 13.7
    Zn 0.04
    Al 0
    density 1
    -water 1 # kg

EQUILIBRIUM_PHASES 1
    Forsterite 0 8
    K-Feldspar 0 20
    Quartz    0 10

SAVE solution 1

SOLUTION 2  Supercritical CO2 addition
    USE solution 1
    temp 85

GAS_PHASE 1
    -fixed_volume
    CO2(g) 250  # Supercritical CO2 at 250 atm

SAVE solution 2

REACTION 1
    CO2(g)     10
    1 moles in 200 steps

RATES
    Quartz
-start
 1  REM  Specific rate k from Rimstidt and Barnes, 1980, GCA 44,1683
 2  REM  k = 10^-13.7 mol/m2/s (25 C), Ea = 90 kJ/mol
 3  REM  sp. rate * parm(2) due to salts (Dove and Rimstidt, MSA Rev. 29, 259)
 4  REM  PARM(1) = Specific area of Quartz, m^2/mol Quartz
 5  REM  PARM(2) = salt correction: (1 + 1.5 * c_Na (mM)), < 35
10 dif_temp = 1/TK - 1/298
20 pk_w = 13.7 + 4700.4 * dif_temp
40 moles = PARM(1) * M0 * PARM(2) * (M/M0)^0.67 * 10^-pk_w * (1 -  SR("Quartz"))
50 SAVE moles * TIME
-end
    K-feldspar
-start
  1   REM Sverdrup and Warfvinge, 1995, mol m^-2 s^-1
  2   REM PARM(1) = Specific area of Kspar m^2/mol Kspar
  3   REM PARM(2) = Adjusts lab rate to field rate
  4   REM temp corr: from A&P, p. 162. E (kJ/mol) / R / 2.303 = H in H*(1/T-1/281)
  5   REM K-Feldspar parameters
 10  DATA 11.7, 0.5, 4e-6, 0.4, 500e-6, 0.15, 14.5, 0.14, 0.15, 13.1, 0.3
 20  RESTORE 10
 30  READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH
 40  DATA 3500, 2000, 2500, 2000
 50  RESTORE 40
 60  READ e_H, e_H2O, e_OH, e_CO2
 70  pk_CO2 = 13
 80  n_CO2 = 0.6
100 REM Generic rate follows
110 dif_temp = 1/TK - 1/281
120 BC       = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2")
130 REM rate by H+
140 pk_H     = pk_H + e_H * dif_temp
150 rate_H   = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC)
160 REM rate by hydrolysis
170 pk_H2O   = pk_H2O + e_H2O * dif_temp
180 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC)
190 REM rate by OH-
200 pk_OH    = pk_OH + e_OH * dif_temp
210 rate_OH  = 10^-pk_OH * ACT("OH-")^o_OH
220 REM rate by CO2
230 pk_CO2   = pk_CO2 + e_CO2 * dif_temp
240 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2
250 rate     = rate_H + rate_H2O + rate_OH + rate_CO2
260 area     = PARM(1) * M0 *(M/M0)^0.67
270 rate     = PARM(2) * area * rate * (1-SR("K-feldspar"))
280 moles    = rate * TIME
290 SAVE moles
-end
    Forsterite
-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=8.38E+04
 12 E1=67206
 13 n1=0.470
 20 rem neutral solution parameters
 21 a2=1.58E+03
 22 E2=79000
 30 rem base solution parameters
 31 a3=1.00E-07
 32 E3=56637
 33 n2=-0.600
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("forsterite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
 90 Rate=(Rate1+Rate2 + Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

KINETICS 1
Quartz
    -formula  Quartz  1
    -m        65
    -m0       65
    -parms    1e-06
    -tol      1e-08
K-Feldspar
    -formula  K-Feldspar  1
    -m        20
    -m0       20
    -parms    1e-05
    -tol      1e-08
Forsterite
    -formula  Mg2SiO4  1
    -m        8
    -m0       8
    -tol      1e-08
-steps       1
-step_divide 1
-runge_kutta 3
-bad_step_max 500

SELECTED_OUTPUT
    -file output.txt
    -reset false
    -time
    -pH
    -alkalinity
    -si Quartz
    -si K-Feldspar
    -si Albite
    -tot Ca
    -tot Mg
    -tot K
    -tot Si
    -tot Na
    -tot Cl
END


[dlparkhurst]

You are using PARM(2) in the Quartz RATES definition, but only one value is given in
-parms in KINETICS.

[zscowcrofst]

Thank you for the suggestion. I fixed the reaction rates but I am not sure if the way that I am trying to setup the introduction of scCO2 (at 250 atm and 85C) to the brine that is in equilibrium with present minerals are correct. I divided the simulation in two step. in step 1 I put the brine in equilibrium with minerals and saved  the solution:

Code: [Select]

EQUILIBRIUM_PHASES 1
    Forsterite 0 10
    K-Feldspar 0 10
    Quartz    0 10

and Then in second simulation I used Solution 1and added the Gas_Phase. Then I added 10 moles of CO2(g) in 1 step. I don't know if this is enough for 200 years of simulations. The results doesn't seem to be correct. Forsterite and Feldspar start under saturated and pH and every other major ions concertation stays constant nd only Cl concentration goes up.

Also how can I add potential precipitant minerals to the simulations?

Code: [Select]

TITLE Sandstone Equilibrium with Brine and Supercritical CO2 Injection

RATES
    Quartz
-start
 1  REM  Specific rate k from Rimstidt and Barnes, 1980, GCA 44,1683
 2  REM  k = 10^-13.7 mol/m2/s (25 C), Ea = 90 kJ/mol
 3  REM  sp. rate * parm(2) due to salts (Dove and Rimstidt, MSA Rev. 29, 259)
 4  REM  PARM(1) = Specific area of Quartz, m^2/mol Quartz
 5  REM  PARM(2) = salt correction: (1 + 1.5 * c_Na (mM)), < 35
10 dif_temp = 1/TK - 1/298
20 pk_w = 13.7 + 4700.4 * dif_temp
40 moles = PARM(1) * M0 * PARM(2) * (M/M0)^0.67 * 10^-pk_w * (1 -  SR("Quartz"))
50 SAVE moles * TIME
-end
    K-Feldspar
-start
  1   REM Sverdrup and Warfvinge, 1995, mol m^-2 s^-1
  2   REM PARM(1) = Specific area of Kspar m^2/mol Kspar
  3   REM PARM(2) = Adjusts lab rate to field rate
  4   REM temp corr: from A&P, p. 162. E (kJ/mol) / R / 2.303 = H in H*(1/T-1/281)
  5   REM K-Feldspar parameters
 10  DATA 11.7, 0.5, 4e-6, 0.4, 500e-6, 0.15, 14.5, 0.14, 0.15, 13.1, 0.3
 20  RESTORE 10
 30  READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH
 40  DATA 3500, 2000, 2500, 2000
 50  RESTORE 40
 60  READ e_H, e_H2O, e_OH, e_CO2
 70  pk_CO2 = 13
 80  n_CO2 = 0.6
100 REM Generic rate follows
110 dif_temp = 1/TK - 1/281
120 BC       = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2")
130 REM rate by H+
140 pk_H     = pk_H + e_H * dif_temp
150 rate_H   = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC)
160 REM rate by hydrolysis
170 pk_H2O   = pk_H2O + e_H2O * dif_temp
180 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC)
190 REM rate by OH-
200 pk_OH    = pk_OH + e_OH * dif_temp
210 rate_OH  = 10^-pk_OH * ACT("OH-")^o_OH
220 REM rate by CO2
230 pk_CO2   = pk_CO2 + e_CO2 * dif_temp
240 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2
250 rate     = rate_H + rate_H2O + rate_OH + rate_CO2
260 area     = PARM(1) * M0 *(M/M0)^0.67
270 rate     = PARM(2) * area * rate * (1-SR("K-Feldspar"))
280 moles    = rate * TIME
290 SAVE moles
-end
    Forsterite
-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=8.38E+04
 12 E1=67206
 13 n1=0.470
 20 rem neutral solution parameters
 21 a2=1.58E+03
 22 E2=79000
 30 rem base solution parameters
 31 a3=1.00E-07
 32 E3=56637
 33 n2=-0.600
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("forsterite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
 90 Rate=(Rate1+Rate2 + Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

SOLUTION 1 Brine in equilibrium with minerals
    temp      85
    pH        7
    pe        4
    redox     pe
    units     mg/l
    density   1
    Al        0
    Alkalinity 149
    B         16.1
    Ca        2704
    Cl        88409
    K         168
    Mg        889
    Mn        0.81
    Na        49706
    S(6)      10.02
    Si        13.7
    Zn        0.04
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    Forsterite 0 10
    K-Feldspar 0 10
    Quartz    0 10

SAVE solution 1

END

SOLUTION 2 Supercritical CO2 addition
    temp      85
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    1 # kg
USE solution 1
    temp 85

GAS_PHASE 1
    -fixed_volume
    -pressure 1
    -volume 1
    -temperature 85
    CO2(g)    250

SAVE solution 2

REACTION 1
    CO2(g)     1
    10 moles in 1 steps

KINETICS 1
Quartz
    -formula  Quartz  1
    -m        1
    -m0       1
    -parms    1e-06 1
    -tol      1e-08
K-Feldspar
    -formula  K-Feldspar  1
    -m        1
    -m0       1
    -parms    1e-05 1
    -tol      1e-08
Forsterite
    -formula  Mg2SiO4  1
    -m        1
    -m0       1
    -parms    1e-05 1
    -tol      1e-08
-steps       6307200000 in 200 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500

SELECTED_OUTPUT
    -file output.txt
    -reset false
    -time
    -pH
    -alkalinity
    -si Quartz
    -si K-Feldspar
    -si Forsterite
    -tot Ca
    -tot Mg
    -tot K
    -tot Si
    -tot Na
    -tot Cl


USER_GRAPH 1
    -headings               Time pH Ca Mg K Si Na Cl
    -axis_titles            "Time (years)" "Concentration (mg/L)" "pH"
    -chart_title            "Major Ion Concentrations over Time"
    -axis_scale sy_axis     auto 10 1 1
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X TOTAL_TIME/3600/24/365.25
20 GRAPH_Y TOT("Ca"), TOT("Mg"),TOT("K"),TOT("Si"),TOT("Na"),TOT("Cl")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true


USER_GRAPH 2
    -headings               Time Gas_Moles SI_Quartz SI_K-Feldspar SI_Forsterite
    -axis_titles            "Time (years)" "Gas (moles)" "Saturation Index"
    -chart_title            "Gas Moels & Saturation Index"
    -axis_scale sy_axis     auto auto 1 1
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X TOTAL_TIME/3600/24/365.25
20 GRAPH_Y GAS("CO2(g)")
30 GRAPH_SY SI("Quartz"),SI("K-Feldspar"),SI("Forsterite")
  -end
    -active                 true

END



[dlparkhurst]

The headings for your graph were out of order; pH is the last item. So, it is not Cl that is increasing, rather it is Si.

The way you did the GAS_PHASE, the pressure of the gas is about 230 atm. If you want a fixed pressure of 250, you can either make it a fixed-pressure gas phase, where the volume will adjust to maintain the fixed pressure, or, use EQUILIBRIUM_PHASES to fix the fugacity.

I do not think you need the REACTION block.

You can add potential precipitates with EQUILIBRIUM_PHASES by specifying a target SI of zero, and zero moles of the mineral.

If things are not working, simplify. Why do you need 200 steps that you do not understand, when you can not understand 10 or 5? The output file is useful. I tend to look at the results of the final calculation first.

Code: [Select]

TITLE Sandstone Equilibrium with Brine and Supercritical CO2 Injection

RATES
    Quartz
-start
 1  REM  Specific rate k from Rimstidt and Barnes, 1980, GCA 44,1683
 2  REM  k = 10^-13.7 mol/m2/s (25 C), Ea = 90 kJ/mol
 3  REM  sp. rate * parm(2) due to salts (Dove and Rimstidt, MSA Rev. 29, 259)
 4  REM  PARM(1) = Specific area of Quartz, m^2/mol Quartz
 5  REM  PARM(2) = salt correction: (1 + 1.5 * c_Na (mM)), < 35
10 dif_temp = 1/TK - 1/298
20 pk_w = 13.7 + 4700.4 * dif_temp
40 moles = PARM(1) * M0 * PARM(2) * (M/M0)^0.67 * 10^-pk_w * (1 -  SR("Quartz"))
50 SAVE moles * TIME
-end
    K-Feldspar
-start
  1   REM Sverdrup and Warfvinge, 1995, mol m^-2 s^-1
  2   REM PARM(1) = Specific area of Kspar m^2/mol Kspar
  3   REM PARM(2) = Adjusts lab rate to field rate
  4   REM temp corr: from A&P, p. 162. E (kJ/mol) / R / 2.303 = H in H*(1/T-1/281)
  5   REM K-Feldspar parameters
 10  DATA 11.7, 0.5, 4e-6, 0.4, 500e-6, 0.15, 14.5, 0.14, 0.15, 13.1, 0.3
 20  RESTORE 10
 30  READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH
 40  DATA 3500, 2000, 2500, 2000
 50  RESTORE 40
 60  READ e_H, e_H2O, e_OH, e_CO2
 70  pk_CO2 = 13
 80  n_CO2 = 0.6
100 REM Generic rate follows
110 dif_temp = 1/TK - 1/281
120 BC       = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2")
130 REM rate by H+
140 pk_H     = pk_H + e_H * dif_temp
150 rate_H   = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC)
160 REM rate by hydrolysis
170 pk_H2O   = pk_H2O + e_H2O * dif_temp
180 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC)
190 REM rate by OH-
200 pk_OH    = pk_OH + e_OH * dif_temp
210 rate_OH  = 10^-pk_OH * ACT("OH-")^o_OH
220 REM rate by CO2
230 pk_CO2   = pk_CO2 + e_CO2 * dif_temp
240 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2
250 rate     = rate_H + rate_H2O + rate_OH + rate_CO2
260 area     = PARM(1) * M0 *(M/M0)^0.67
270 rate     = PARM(2) * area * rate * (1-SR("K-Feldspar"))
280 moles    = rate * TIME
290 SAVE moles
-end
    Forsterite
-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=8.38E+04
 12 E1=67206
 13 n1=0.470
 20 rem neutral solution parameters
 21 a2=1.58E+03
 22 E2=79000
 30 rem base solution parameters
 31 a3=1.00E-07
 32 E3=56637
 33 n2=-0.600
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("forsterite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
 90 Rate=(Rate1+Rate2 + Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end
END
SOLUTION 1 Brine in equilibrium with minerals
    temp      85
    pH        7
    pe        4
    redox     pe
    units     mg/l
    density   1
    Al        0
    Alkalinity 149
    B         16.1
    Ca        2704
    Cl        88409
    K         168
    Mg        889
    Mn        0.81
    Na        49706
    S(6)      10.02
    Si        13.7
    Zn        0.04
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    Forsterite 0 10
    K-Feldspar 0 10
    Quartz    0 10

SAVE solution 1

END

SOLUTION 2 Supercritical CO2 addition
    temp      85
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    1 # kg
END
USE solution 1
    temp 85

GAS_PHASE 1
    -fixed_volume
    -pressure 1
    -volume 1
    -temperature 85
    CO2(g)    250

SAVE solution 2

#REACTION 1
#    CO2(g)     1
#    10 moles in 1 steps

KINETICS 1
Quartz
    -formula  Quartz  1
    -m        1
    -m0       1
    -parms    1e-06 1
    -tol      1e-08
K-Feldspar
    -formula  K-Feldspar  1
    -m        1
    -m0       1
    -parms    1e-05 1
    -tol      1e-08
Forsterite
    -formula  Mg2SiO4  1
    -m        1
    -m0       1
    -parms    1e-05 1
    -tol      1e-08
-steps       6307200000 in 10 #200 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500


USER_GRAPH 1
    -headings               Time Ca Mg K Si Na Cl pH
    -axis_titles            "Time (years)" "Concentration (mg/L)" "pH"
    -chart_title            "Major Ion Concentrations over Time"
    -axis_scale y_axis      auto auto auto auto log
    #-axis_scale sy_axis     auto 10 1 1
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X TOTAL_TIME/3600/24/365.25
20 GRAPH_Y TOT("Ca"), TOT("Mg"),TOT("K"),TOT("Si"),TOT("Na"),TOT("Cl")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true


USER_GRAPH 2
    -headings               Time Gas_Moles SI_Quartz SI_K-Feldspar SI_Forsterite
    -axis_titles            "Time (years)" "Gas (moles)" "Saturation Index"
    -chart_title            "Gas Moles & Saturation Index"
    -axis_scale sy_axis     auto auto 1 1
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X TOTAL_TIME/3600/24/365.25
20 GRAPH_Y GAS("CO2(g)")
30 GRAPH_SY SI("Quartz"),SI("K-Feldspar"),SI("Forsterite")
  -end
    -active                 true

END

[zscowcrofst]

Thank you for the guidance. I have changed the model and added calcite so that it can precipitate kinetically. I also changed the gas phase from fixed volume to fixed pressure. Additionally, I deleted Solution 2 because it was not being used. Now, Solution 2 is the initial brine in equilibrium with the initial minerals (quartz, K-feldspar, and forsterite), which is saved as Solution 2.

1- I wonder why, despite the solution being in equilibrium with quartz, forsterite, and K-feldspar, the initial SI of these minerals is not starting at 0 (Graph 2).
2- In Graph 3, I used the KIN function to calculate the mineral phase volume. Why is there a change in the phase volume?
 

Code: [Select]

TITLE Sandstone Equilibrium with Brine and Supercritical CO2 Injection

RATES
    Quartz
-start
 1  REM  Specific rate k from Rimstidt and Barnes, 1980, GCA 44,1683
 2  REM  k = 10^-13.7 mol/m2/s (25 C), Ea = 90 kJ/mol
 3  REM  sp. rate * parm(2) due to salts (Dove and Rimstidt, MSA Rev. 29, 259)
 4  REM  PARM(1) = Specific area of Quartz, m^2/mol Quartz
 5  REM  PARM(2) = salt correction: (1 + 1.5 * c_Na (mM)), < 35
10 dif_temp = 1/TK - 1/298
20 pk_w = 13.7 + 4700.4 * dif_temp
40 moles = PARM(1) * M0 * PARM(2) * (M/M0)^0.67 * 10^-pk_w * (1 -  SR("Quartz"))
50 SAVE moles * TIME
-end
    K-Feldspar
-start
  1   REM Sverdrup and Warfvinge, 1995, mol m^-2 s^-1
  2   REM PARM(1) = Specific area of Kspar m^2/mol Kspar
  3   REM PARM(2) = Adjusts lab rate to field rate
  4   REM temp corr: from A&P, p. 162. E (kJ/mol) / R / 2.303 = H in H*(1/T-1/281)
  5   REM K-Feldspar parameters
 10  DATA 11.7, 0.5, 4e-6, 0.4, 500e-6, 0.15, 14.5, 0.14, 0.15, 13.1, 0.3
 20  RESTORE 10
 30  READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH
 40  DATA 3500, 2000, 2500, 2000
 50  RESTORE 40
 60  READ e_H, e_H2O, e_OH, e_CO2
 70  pk_CO2 = 13
 80  n_CO2 = 0.6
100 REM Generic rate follows
110 dif_temp = 1/TK - 1/281
120 BC       = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2")
130 REM rate by H+
140 pk_H     = pk_H + e_H * dif_temp
150 rate_H   = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC)
160 REM rate by hydrolysis
170 pk_H2O   = pk_H2O + e_H2O * dif_temp
180 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC)
190 REM rate by OH-
200 pk_OH    = pk_OH + e_OH * dif_temp
210 rate_OH  = 10^-pk_OH * ACT("OH-")^o_OH
220 REM rate by CO2
230 pk_CO2   = pk_CO2 + e_CO2 * dif_temp
240 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2
250 rate     = rate_H + rate_H2O + rate_OH + rate_CO2
260 area     = PARM(1) * M0 *(M/M0)^0.67
270 rate     = PARM(2) * area * rate * (1-SR("K-Feldspar"))
280 moles    = rate * TIME
290 SAVE moles
-end
    Forsterite
-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=8.38E+04
 12 E1=67206
 13 n1=0.470
 20 rem neutral solution parameters
 21 a2=1.58E+03
 22 E2=79000
 30 rem base solution parameters
 31 a3=1.00E-07
 32 E3=56637
 33 n2=-0.600
 36 rem rate=0 if no minerals and undersaturated
 40 SR_mineral=SR("forsterite")
 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
 80 Rate2=a2*EXP(-E2/R/TK)
 85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
 90 Rate=(Rate1+Rate2 + Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end
END
SOLUTION 1 Brine in equilibrium with minerals
    temp      85
    pH        7
    pe        4
    redox     pe
    units     mg/l
    density   1
    Al        0
    Alkalinity 149
    B         16.1
    Ca        2704
    Cl        88409
    K         168
    Mg        889
    Mn        0.81
    Na        49706
    S(6)      10.02
    Si        13.7
    Zn        0.04
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    Forsterite 0 10
    K-Feldspar 0 10
    Quartz    0 10

SAVE solution 1

END

USE solution 1

GAS_PHASE 1
    -fixed_pressure
    -pressure 1
    -volume 1
    -temperature 85
    CO2(g)    250

SAVE solution 2

KINETICS 1
Quartz
    -formula  Quartz  1
    -m        1
    -m0       1
    -parms    1e-06 1
    -tol      1e-08
K-Feldspar
    -formula  K-Feldspar  1
    -m        1
    -m0       1
    -parms    1e-05 1
    -tol      1e-08
Forsterite
    -formula  Mg2SiO4  1
    -m        1
    -m0       1
    -parms    1e-05 1
    -tol      1e-08
Calcite
    -formula  CaCO3  1
    -m        0
    -m0       0
    -parms    167000 0.6
    -tol      1e-08
-steps       6307200000 in 10 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500

USER_GRAPH 1
    -headings               Time Ca Mg K Si Na Cl pH
    -axis_titles            "Time (years)" "Concentration (mg/L)" "pH"
    -chart_title            "Major Ion Concentrations over Time"
    -axis_scale y_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X TOTAL_TIME/3600/24/365.25
20 GRAPH_Y TOT("Ca"), TOT("Mg"),TOT("K"),TOT("Si"),TOT("Na"),TOT("Cl")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true
USER_GRAPH 2
    -headings               Time Quartz Forsterite K-Feldspar Calcite
    -axis_titles            "Time (years)" "Saturation Index" ""
    -chart_title            "Gas Moles & Saturation Index"
    -axis_scale sy_axis     auto auto 1 1
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X TOTAL_TIME/3600/24/365.25
20 GRAPH_SY SI("Quartz"),SI("Forsterite"),SI("K-Feldspar"), SI("Calcite")
  -end
    -active                 true

USER_GRAPH 3
    -headings               x Quartz Forsterite K-Feldspar Calcite
    -axis_titles            "Phase volumes (cm3/100g)" "" ""
    -axis_scale y_axis      auto auto 10 auto
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 V_Quartz = KIN("Quartz") * 60.08
20 V_Forsterite = KIN("Forsterite") * 140.69
30 V_KFeldspar = KIN("K-Feldspar")*278.33
40 V_Calcite = KIN("Calcite")*100.0869
50 GRAPH_X TOTAL_TIME/3600/24/365.25
60 GRAPH_Y V_Quartz,V_Forsterite,V_KFeldspar, V_Calcite
  -end
    -active                 true

END

[dlparkhurst]

Your first equilibrium calculation for Forsterite, K-feldspar, and Quartz is at 85 C, 1 atm pressure, and no gas phase.

The reaction calculations plotted on your graphs are in the presence of the CO2 gas phase, which would change the saturation indices even before any kinetic reaction occurs.

The SAVE solution 2 has no effect.

The X axis in USER_GRAPH 3 is time; the Y axis is volume. The changes in volume are quite small relative to the total amounts of each kinetic reactant, on the order of 1e-3 or less.

[zscowcrofst]

Thank you David.
I have added a GAS_PHASE to the solution 1 with Ar(g) as an inert gas to bring the pressure to 250atm. But still the starting point of the solution is not looking at equilibrium.

I am also not sure if I defined Calcite kinetic precipitation correctly. I used the Rate expression in llnl database, assuming it is not present initially, m0 = 0 and m = 0.

Code: [Select]

########
#Calcite
########
# Example of KINETICS data block for calcite rate,
#   in mmol/cm2/s, Plummer et al., 1978, AJS 278, 179; Appelo et al., AG 13, 257.
# KINETICS 1
# Calcite
# -tol   1e-8
# -m0    3.e-3
# -m     3.e-3
# -parms 1.67e5   0.6  # cm^2/mol calcite, exp factor
# -time  1 day

Code: [Select]

SOLUTION 1 Brine in equilibrium with minerals
    temp      85
    pH        7
    pe        4
    redox     pe
    units     mg/l
    density   1
    Al        0
    Alkalinity 149
    B         16.1
    Ca        2704
    Cl        88409
    K         168
    Mg        889
    Mn        0.81
    Na        49706
    S(6)      10.02
    Si        13.7
    Zn        0.04
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    Forsterite 0 10
    K-Feldspar 0 10
    Quartz    0 10

GAS_PHASE 1
    -fixed_pressure
    -pressure 250
    -volume 1
    -temperature 85
    Ar(g)     250

SAVE solution 1

END

USE solution 1

GAS_PHASE 2
    -fixed_pressure
    -pressure 250
    -volume 1
    -temperature 85
    CO2(g)    250

KINETICS 1
Quartz
    -formula  Quartz  1
    -m        1
    -m0       1
    -parms    1e-06 1
    -tol      1e-08
K-Feldspar
    -formula  K-Feldspar  1
    -m        1
    -m0       1
    -parms    1e-05 1
    -tol      1e-08
Forsterite
    -formula  Mg2SiO4  1
    -m        1
    -m0       1
    -parms    1e-05 1
    -tol      1e-08
Calcite
    -formula  CaCO3  1
    -m        0
    -m0       0
    -parms    167000 0.6
    -tol      1e-08
-steps       6307200000 in 20 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500

USER_GRAPH 2
    -headings               Time Quartz Forsterite K-Feldspar Calcite
    -axis_titles            "Time (years)" "Saturation Index" ""
    -chart_title            "Saturation Index"
    -axis_scale sy_axis     auto auto 1 1
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X TOTAL_TIME/3600/24/365.25
20 GRAPH_SY SI("Quartz"),SI("Forsterite"),SI("K-Feldspar"), SI("Calcite")
  -end
    -active                 true

USER_GRAPH 3
    -headings               x Quartz Forsterite K-Feldspar Calcite
    -axis_titles            "Phase volumes (cm3/100g)" "" ""
    -chart_title            "Minerals Volum Fraction"
    -axis_scale y_axis      auto auto auto 10 log
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 V_Quartz = KIN("Quartz") * 60.08
20 V_Forsterite = KIN("Forsterite") * 140.69
30 V_KFeldspar = KIN("K-Feldspar")*278.33
40 V_Calcite = KIN("Calcite")*100.0869
50 GRAPH_X TOTAL_TIME/3600/24/365.25
60 GRAPH_Y V_Quartz,V_Forsterite,V_KFeldspar, V_Calcite
  -end
    -active                 true

END


[dlparkhurst]

As soon as you add GAS_PHASE 2 to your calculation, the system is no longer at equilibrium, and that applies to all of your kinetic calculations. If you want the initial equilibrium condition, move your USER_GRAPH definitions up to the place where the initial equilbrium is calculated.

[zscowcrofst]

Thank you that worked.