basalt dissolution

[onat.yakisik]

Hello,

I am trying to solve basalt in a 0.5 M solution at 150 degrees and 40 bar pressure with a 30-day kinetic reaction. I am using core10.dat. But when i try to simulate the code, it gives me strange results. Is it possible if you can check it and help me if there is an error? I am sharing the code below.

Code: [Select]

#Phases and kinetics from (Zhang et al. 2019
DATABASE C:\phreeqc\database\core10.dat
PHASES
Forsterite
Mg2SiO4 + 4 H+ =  SiO2 + 2 H2O + 2 Mg+2
log_k 27.8626
-delta_H -205.614 kJ/mol
# deltafH -520 kcal/mol
-analytic -7.6195e1 -1.4013e-2 1.4763e4 2.5090e1 -3.0379e5
# Range 0-350
-Vm 43.79
# Extrapol supcrt92
# Ref HDN+78

Fayalite
Fe2SiO4 + 4 H+ = SiO2 + 2 Fe+2 + 2 H2O
log_k 19.1113
-delta_H -152.256 kJ/mol
# deltafH -354.119 kcal/mol
-analytic 1.3853e1 -3.5501e-3 7.1496e3 -6.8710e0 -6.3310e4
# Range 0-350
-Vm 46.39
# Extrapol supcrt92
# Ref HDN+78

Diopside
CaMgSi2O6 + 4 H+ = Ca+2 + Mg+2 + 2 H2O + 2 SiO2
log_k 20.9643
-delta_H -133.775 kJ/mol
# deltafH -765.378 kcal/mol
-analytic 7.1240e1 1.5514e-2 8.1437e3 -3.0672e1 -5.6880e5
# Range 0-350
-Vm 66.09
# Extrapol supcrt92
# Ref HDN+78

Enstatite
MgSiO3 + 2 H+ = H2O + Mg+2 + SiO2
log_k 11.3269
-delta_H -82.7302 kJ/mol
# deltafH -369.686 kcal/mol
-analytic -4.9278e1 -3.2832e-3 9.5205e3 1.4437e1 -5.4324e5
# Range 0-350
-Vm 31.276
# Extrapol supcrt92
# Ref HDN+78

Ferrosilite
FeSiO3 + 2 H+ = Fe+2 + H2O + SiO2
log_k 7.4471
-delta_H -60.6011 kJ/mol
# deltafH -285.658 kcal/mol
-analytic 9.0041 3.7917e-3 5.1625e3 -6.3009 -3.9565e5
# Range 0-350
-Vm 32.952
# Extrapol supcrt92
# Ref HDN+78

Anorthite
CaAl2(SiO4)2 + 8 H+ = Ca+2 + 2 Al+3 + 2 SiO2 + 4 H2O
log_k 26.5780
-delta_H -303.039 kJ/mol
# deltafH -1007.55 kcal/mol
-analytic 3.9717e-1 -1.8751e-2 1.4897e4 -6.3078 -2.3885e5
# Range 0-350
-Vm 100.79
# Extrapol supcrt92
# Ref HDN+78

Hematite
Fe2O3 + 6 H+ = 2 Fe+3 + 3 H2O
log_k 0.1086
-delta_H -129.415 kJ/mol
# deltafH -197.72 kcal/mol
-analytic -2.2015e2 -6.0290e-2 1.1812e4 8.0253e1 1.8438e2
# Range 0-350
-Vm 30.274
# Extrapol supcrt92
# Ref HDN+78

Magnetite
Fe3O4 + 8 H+ = Fe+2 + 2 Fe+3 + 4 H2O
log_k 10.4724
-delta_H -216.597 kJ/mol
# deltafH -267.25 kcal/mol
-analytic -3.0510e2 -7.9919e-2 1.8709e4 1.1178e2 2.9203e2
# Range 0-350
-Vm 44.524
# Extrapol supcrt92
# Ref HDN+78


RATES
Forsterite
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  #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

Fayalite
# from Palandri and Kharaka 2004
# experimental condition range T=8-25C, pH=1.3-11

-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=5.48E+11
12 E1=94400
13 n1=0
20 rem neutral solution parameters
21 a2=5.48E+02
22 E2=94400
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("Fayalite")
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

Diopside
# from Palandri and Kharaka 2004
# experimental condition range T=8-90C, pH=1-6

-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=3.00E+10
12 E1=96100
13 n1=0.710
20 rem neutral solution parameters
21 a2=1.00E-04
22 E2=40600
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("diopside")
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

enstatite
# from Palandri and Kharaka 2004
# experimental condition range T=8-72C, pH=1-12

-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=1.00E+05
12 E1=80000
13 n1=0.600
20 rem neutral solution parameters
21 a2=2.00E+01
22 E2=80000
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("enstatite")
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

Anorthite
# from Palandri and Kharaka 2004
# experimental condition range T=25-95C, pH=2-10.2

-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=2.58E-01
12 E1=16601
13 n1=1.411
20 rem neutral solution parameters
21 a2=1.00E-06
22 E2=17821
30 rem base solution parameters
31 a3=1.00E-22
32 E3=18150
33 n2=-1.767
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("anorthite")
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


magnesite
# from Palandri and Kharaka 2004
# experimental condition range T=25C, pH=0.2-12
# calcite activation energy is assumed
# near equilibrium parameters p=4.0 and q=1.0

-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=1.39E-4
12 E1=14400
13 n1=1.000
20 rem neutral solution parameters
21 a2=5.99E-6
22 E2=23500
30 rem CO2 denpendence parameters
31 a3=6.03E+05
32 E3=62800
33 n2=1.000
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("magnesite")
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)*SR("CO2(g)")^n2    #CO2 rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end


hematite
# from Palandri and Kharaka 2004
# experimental condition range T=25-50C, pH=0-5

-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=161.6
12 E1=66200
13 n1=01
20 rem neutral solution parameters
21 a2=9.96E-4
22 E2=66210
30 rem basic dependence parameters
31 a3=0
32 E3=0
33 n2=0
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("hematite")
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


SOLUTION 1
water 2
pH 7
temp 150
pressure 40
units mol/L

Na 0.5
Cl 0.5

# Si 1.0
# Mg 1.0
# O(0) 4.0

EQUILIBRIUM_PHASES 1 basalt minerals
Forsterite 0 0.2
Diopside   0 0.2757
Anorthite 0 0.4743
Hematite 0 0.025 # Assuming the half is hematite (Fe2O3)
Magnetite 0 0.025 # Assuming the other half of the magmatic oxides are magnetite (Fe3O4)

SAVE SOLUTION 1
SAVE EQUILIBRIUM_PHASES 1

 KINETICS 1
Forsterite
       -formula Mg2SiO4 1.0
     
       -m 2e-1
       -M0 2e-1/30 # moles of solid per kg of water
      -parms 1e3 0.1 # total surface area per kg of water (m2/kgw) and the scaling factor
Diopside
       -formula CaMgSi2O6 1.0
       -m 2.757e-1
       -M0 2.757e-1/30  # moles of solid per kg of water
      -parms 1e3 0.1 # total surface area per kg of water (m2/kgw) and the scaling factor
Anorthite
       -formula CaAl2(SiO4)2 1.0
       -m 4.4743e-1
       -M0 4.4743e-1/30  # moles of solid per kg of water
       -parms 1e3 0.1 # total surface area per kg of water (m2/kgw) and the scaling factor
-time 30.0
-steps  30*86400 # define time steps
-step_divide 1e-2
-cvode true
REACTION
Forsterite 0.0066
Diopside 0.0099
Anorthite 0.01581
INCREMENTAL_REACTIONS TRUE

USE SOLUTION 1
USE EQUILIBRIUM_PHASES 1
USER_GRAPH 1
-chart_title "Forsterite Dissolution"
-axis_titles Day "mol"
-initial_solutions false
-start
10 graph_x total_time/86400
20 graph_y KIN_DELTA("Forsterite")
-end

USER_GRAPH 2
-chart_title "Diopside Dissolution"
-axis_titles Day "mol"
-initial_solutions false
-start
10 graph_x total_time/86400
20 graph_y KIN_DELTA("Diopside")
-end

USER_GRAPH 3
-chart_title "Anorthite Dissolution"
-axis_titles Day "mol"
-initial_solutions false
-start
10 graph_x total_time/86400
20 graph_y KIN_DELTA("Anorthite")
-end
END

[dlparkhurst]

You have include Anorthite, Diopside, and Forsterite in both EQUILIBRIUM_PHASES and KINETICS. If you look at the results, all three minerals react to equilibrium because of the EQUILIBRIUM_PHASES definition and the KINETICS has no effect because the rates are zero at equilibrium. If you want to see a kinetic reaction, you need to set up a water composition that is not in equilibrium with the minerals and use KINETICS for the minerals you want to react kinetically (and exclude these minerals from EQUILIBRIUM_PHASES).

[Divya]

Hello everyone,
I am trying to write a code that include dissolution of wollastonite mineral and precipitation of calcite as we inject CO2 in water  and 0.5M NaHCO3 solution ,with reaction temperature 90C and reaction pressure 200 atm and initial amount of wollastonite is 50 mole. As I'm injecting CO2 into this with unkown amount for 1 and 5 days. As I want to estimate how much amount of CO2 is requied for max precipitation.
I will appreciate anyone who can help me it.

[dlparkhurst]

You can start with the equilibrium case, where you define a NaHCO3 solution and allow wollastonite and calcite to react to equilibrium with EQUILIBRIUM_PHASES as you add CO2 with a REACTION.

[Divya]

As I'm begineer. Can you help me with this?

Code: [Select]

[code]
RATES
    Wollastonite
-start
 10 rem unit should be mol,Liter-1 and second-1
 20 rem parm(1) is surface area in the unit of m2/L
 30 rem calculation of surface area can be found in the note
 40 rem M is current moles of minerals M0 is the initial moles of minerals
 50 rem parm(2) is a scaling factor
 60 rem acid solution parameters
 70 k1=1.60E+04
 80 E1=54651
 90 n1=0.400
100 rem neutral solution parameters
110 k2=5
120 E2=54651
130 rem base solution parameters
140 k3=0
150 E3=0
160 n2=0
170 rem rate=0 if no minerals and undersaturated
180 Si_mineral=SI("wollastonite")
190 if (M<0) then goto 290
200 if (M=0 and Si_mineral<1) then goto 290
210 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
220 if (SA<=0) then SA=1
230 R=8.31451
240 Rate1=k1*EXP(-E1/R/TK)*ACT("H+")^n1
250 Rate2=k2*EXP(-E2/R/TK)
260 Rate3=k3*EXP(-E3/R/TK)*ACT("H+")^n2
270 Rate=(Rate1+Rate2+Rate3)*(1-Si_mineral)*SA*parm(2)
280 moles= rate*Time
290 save moles
-end

SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Na        0.5
    C(4)      0.5
    -water    1 # kg
GAS_PHASE 1
    -fixed_pressure
    -pressure 60
    -volume 1
    -temperature 25
    CO2(g)    60
PHASES
Wollastonite
    CaSiO3 + 2H+ + H2O = Ca+2 + H4SiO4
    log_k     12.996
    delta_h   -19.498 kcal
    -analytical_expression 525.2449 0.2616002 -6403.288 -224.8401 -174483.8 -0.0001135532
    -Vm       39.93 cm3/mol
CO2(g)
    CO2 + H2O = CO3-2 + 2H+
    log_k     -18.16
    delta_h   0.53 kcal
    -T_c      304.2
    -P_c      72.86
    -Omega    0.225
EQUILIBRIUM_PHASES 1
    Wollastonite 0 CaSiO3    50
KINETICS 1
Wollastonite
    -formula  CaSiO3  1
    -m        0
    -m0       50
    -parms    0.1 0.5
    -tol      1e-06
-steps       172800 in 1 steps # seconds
-step_divide 1
-runge_kutta 2
-bad_step_max 500
INCREMENTAL_REACTIONS True
REACTION_TEMPERATURE 1
    90
REACTION_PRESSURE 1
    200
END


[dlparkhurst]

This should get you started. Don't include Wollastonite in KINETICS and EQUILIBRIUM_PHASES, choose one or the other. Give Wollastonite a non-zero M. I've included Calcite in equilibrium phase.

Code: [Select]

RATES
    Wollastonite
-start
 10 rem unit should be mol,Liter-1 and second-1
 20 rem parm(1) is surface area in the unit of m2/L
 30 rem calculation of surface area can be found in the note
 40 rem M is current moles of minerals M0 is the initial moles of minerals
 50 rem parm(2) is a scaling factor
 60 rem acid solution parameters
 70 k1=1.60E+04
 80 E1=54651
 90 n1=0.400
100 rem neutral solution parameters
110 k2=5
120 E2=54651
130 rem base solution parameters
140 k3=0
150 E3=0
160 n2=0
170 rem rate=0 if no minerals and undersaturated
180 Si_mineral=SI("wollastonite")
190 if (M<0) then goto 290
200 if (M=0 and Si_mineral<1) then goto 290
210 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
220 if (SA<=0) then SA=1
230 R=8.31451
240 Rate1=k1*EXP(-E1/R/TK)*ACT("H+")^n1
250 Rate2=k2*EXP(-E2/R/TK)
260 Rate3=k3*EXP(-E3/R/TK)*ACT("H+")^n2
270 Rate=(Rate1+Rate2+Rate3)*(1-Si_mineral)*SA*parm(2)
280 moles= rate*Time
290 save moles
-end

SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Na        0.5
    C(4)      0.5
    -water    1 # kg
GAS_PHASE 1
    -fixed_pressure
    -pressure 60
    -volume 1
    -temperature 25
    CO2(g)    60
PHASES
Wollastonite
    CaSiO3 + 2H+ + H2O = Ca+2 + H4SiO4
    log_k     12.996
    delta_h   -19.498 kcal
    -analytical_expression 525.2449 0.2616002 -6403.288 -224.8401 -174483.8 -0.0001135532
    -Vm       39.93 cm3/mol
CO2(g)
    CO2 + H2O = CO3-2 + 2H+
    log_k     -18.16
    delta_h   0.53 kcal
    -T_c      304.2
    -P_c      72.86
    -Omega    0.225
EQUILIBRIUM_PHASES 1
#    Wollastonite 0 CaSiO3    50
    Calcite 0 10
KINETICS 1
Wollastonite
    -formula  CaSiO3  1
    -m        50
    -m0       50
    -parms    0.1 0.5
    -tol      1e-06
-steps       172800 in 10 steps # seconds
-step_divide 1
-runge_kutta 2
-bad_step_max 500
INCREMENTAL_REACTIONS True
REACTION_TEMPERATURE 1
    90
REACTION_PRESSURE 1
    200
USER_GRAPH 1
    -headings               time Delta_Wollastonite SI_Wollastonite
    -axis_titles            "Time, seconds" "Delta Wollastonite" "SI Wollastonite"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y KIN_DELTA("Wollastonite")
30 GRAPH_SY SI("Wollastonite")
  -end
    -active                 true
END

[Divya]

Hi...
I want to know how much time it will take to dissolve wollastonite completely and form carbonates. at different temperature and CO2 partial pressure of 35?C , 65?C , 90?C  and 100atm, 150atm , 200atm respectively And I also want to make graph of pH evolution,  calcium leaching
profile, calcite precipitation, and  mineral carbonation efficiency at a fixed temperatures with respect to time, temperature and pressure  seperately. Can you help me with the above code by modifying according to the above statement.
I will appriciate your help.

[dlparkhurst]

No, I think you should try first.

[Divya]

Ihave already tried but not getting the result. The dissolution of wollastonite is getting stop after Reaction step 1. It is not getting dissolve completely. Please help me with this.I will truly appreciate your help.

Code: [Select]

RATES
    Wollastonite
-start
 10 rem unit should be mol,Liter-1 and second-1
 20 rem parm(1) is surface area in the unit of m2/g
 30 rem calculation of surface area can be found in the note
 40 rem M is current moles of minerals M0 is the initial moles of minerals
 50 rem parm(2) is a scaling factor
 60 rem acid solution parameters
 70 k1=4.266E-06
 80 E1=54651
 90 n1=0.400
100 rem neutral solution parameters
110 k2=1.32E-09
120 E2=54651
130 rem base solution parameters
140 k3=1.32E-09
150 E3=54651
160 n2=0.5
170 rem rate=0 if no minerals and undersaturated
180 Si_mineral=SI("wollastonite")
190 if (M<0) then goto 290
200 if (M=0 and Si_mineral<1) then goto 290
210 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
220 if (SA<=0) then SA=1
230 R=8.31451
240 Rate1=k1*EXP(-E1/R/TK)*ACT("H+")^n1
250 Rate2=k2*EXP(-E2/R/TK)
260 Rate3=k3*EXP(-E3/R/TK)*ACT("H+")^n2
270 Rate=(Rate1+Rate2+Rate3)*(1-Si_mineral)*SA*parm(2)
280 moles= rate*Time
290 save moles
-end

SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Na        0.5
    C(4)      0.5
    -water    1 # kg
GAS_PHASE 1
    -fixed_pressure
    -pressure 150
    -volume 1
    -temperature 25
    CO2(g)    150
PHASES
Wollastonite
    CaSiO3 + 2H+ + H2O = Ca+2 + H4SiO4
    log_k     12.996
    delta_h   -19.498 kcal
    -analytical_expression 525.2449 0.2616002 -6403.288 -224.8401 -174483.8 -0.0001135532
    -Vm       39.93 cm3/mol
CO2(g)
    CO2 + H2O = CO3-2 + 2H+
    log_k     -18.16
    delta_h   0.53 kcal
    -T_c      304.2
    -P_c      72.86
    -Omega    0.225
Calcite
     CaCO3 = Ca+2 + CO3-2
    log_k     -8.48
    delta_h   -2.297 kcal
    -analytical_expression -171.9065 -0.077993 2839.319 71.595 0 0
    -Vm       36.9 cm3/mol
EQUILIBRIUM_PHASES 1
    Calcite   0 10
KINETICS 1
Wollastonite
    -formula  CaSiO3  1
    -m        50
    -m0       50
    -parms    0.1 0.5
    -tol      1e-08
-steps       864000 in 500 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
REACTION_TEMPERATURE 1
    90
INCREMENTAL_REACTIONS True
USER_GRAPH 1
    -axis_titles            "Time_seconds" "Delta Wollastonite" "SI Wollastonite"
    -chart_title            "Time Delta_Wollastonite SI_Wollastonite"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y KIN_DELTA("Wollastonite")
30 GRAPH_SY SI("Wollastonite")
40 end
  -end
    -active                 true
END


[dlparkhurst]

I think your are misinterpreting this line:

Code: [Select]

-steps       864000 in 500 steps # seconds

It means divide 864000 seconds into 500 steps. The total time integrated will be 10 days.

[dlparkhurst]

Looking a little more closely, there are a couple other problems. Your rate is very slow, requiring more than a billion years to reach equilibrium. In addition, your equilibrium condition is incorrect, you should use SR rather than SI in the following:

Code: [Select]

270 Rate=(Rate1+Rate2+Rate3)*(1-SR("Wollastonite"))*SA*parm(2)


[Divya]

Thank you so much for help.
I tried as you said but I am not understanding one thing that why calcite final moles are decreasing as we are performing wollastonite dissolution and calcite precipitation.Can you help me with code?

[dlparkhurst]

Your Wollanstonite is really soluble--to the tune of several moles per kilogram water. I don't think that is correct.

Wollansonite dissolves, increasing calcium concentration; calcite is then supersaturated and precipitates.  Moles of calcite starts at 10 and increases as shown in the "Phase assemblage" section of the output.

[Divya]

Thank you for such a great help.
I also want to do similar simulation for forsterite mineral

Code: [Select]

DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.7.3-15968\database\phreeqc.dat
SOLUTION_SPECIES
CO3-2 + 2H+ = CO2 + H2O
    log_k     16.681
    delta_h   -5.738 kcal
    analytical_expression 464.1965 0.09344813 -26986.16 -165.75951 2248628.9 0
    dw       1.92e-09
    millero  7.29 0.92 2.07 -1.23 -1.6 0
2CO2 = (CO2)2
    log_k     -1.8
    analytical_expression 8.68 -0.0103 -2190 0 0 0
    millero  14.58 1.84 4.14 -2.46 -3.2 0
2H2O + SiO2 = H4SiO4
    log_k     -2.7
    delta_h   -209.775 kcal
    millero  1.9 1.7 20 -2.7 0.1291 0
H4SiO4 = H3SiO4- + H+
    log_k     -9.93
    delta_h   8.935 kcal
    analytical_expression 6.368 -0.016346 -3405.9 0 0 0
H4SiO4 = H2SiO4-2 + 2H+
    log_k     -21.619
    delta_h   29.714 kcal
    analytical_expression 39.478 -0.065927 -12355.1 0 0 0
RATES
    Forsterite
start
10 rem unit should be mol,kgw-1 and second-1
20 rem parm(1) is surface area in the unit of m2/kgw
30 rem calculation of surface area can be found in the note
40 rem M is current moles of minerals M0 is the initial moles of minerals
50 rem parm(2) is a scaling factor
60 rem acid solution parameters
70 a1=8.38E+04
80 E1=67206
90 n1=0.470
100 rem neutral solution parameters
110 a2=1.58E+03
120 E2=79000
130 rem base solution parameters
140 a3=1.00E-07
150 E3=56637
160 n2=-0.600
170 rem rate=0 if no minerals and undersaturated
180 SR_mineral=SR("forsterite")
190 if (M<0) then goto 310
200 if (M=0 and SR_mineral<1) then goto 310
210 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
220 if (SA<=0) then SA=1
230 R=8.31451
240 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1
250 Rate2=a2*EXP(-E2/R/TK)
260 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
270 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
280 moles= rate*Time
290 rem do not dissolve more minerals than present
300 if (moles>M) then moles=M
310 save moles
end
    Magnesite
start
10 rem unit should be mol,kgw-1 and second-1
20 rem parm(1) is surface area in the unit of m2/kgw
30 rem calculation of surface area can be found in the note
40 rem M is current moles of minerals M0 is the initial moles of minerals
50 rem parm(2) is a scaling factor
60 rem acid solution parameters
70 a1=1.39E-4
80 E1=14400
90 n1=1.000
100 rem neutral solution parameters
110 a2=5.99E-6
120 E2=23500
130 rem CO2 denpendence parameters
140 a3=6.03E+05
150 E3=62800
160 n2=1.000
170 rem rate=0 if no minerals and undersaturated
180 SR_mineral=SR("magnesite")
190 if (M<0) then goto 310
200 if (M=0 and SR_mineral<1) then goto 310
210 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
220 if (SA<=0) then SA=1
230 R=8.31451
240 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1
250 Rate2=a2*EXP(-E2/R/TK)
260 Rate3=a3*EXP(-E3/R/TK)*SR("CO2(g)")^n2
270 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral^4)*SA*parm(2)
280 moles= rate*Time
290 rem do not dissolve more minerals than present
300 if (moles>M) then moles=M
310 save moles
end
PHASES
CO2(g)
    CO2 = CO2
    log_k     -1.468
    delta_h   -4.776 kcal
    analytical_expression 10.5624 -0.023547 -3972.8 0 587460 1.9194e-05
    t_c      304.2
    p_c      72.86
    omega    0.225
Forsterite
    Mg2SiO4 + 4H+ = 2H2O + 2Mg+2 + SiO2
    log_k     27.8626
    delta_h   -520 kcal
    analytical_expression -76.195 -0.014013 14763 25.09 -303790 0
    vm       32.952 cm3/mol
Magnesite
    MgCO3 + H+ = HCO3- + Mg+2
    log_k     2.2936
    delta_h   -265.63 kJ
    analytical_expression -166.65 -0.049469 6434.4 65.506 100.45 0
    vm       28.018 cm3/mol
SOLUTION 1
    temp      35
    pH        7
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    water    1 # kg
GAS_PHASE 1
    fixed_pressure
    pressure 100
    volume 1
    temperature 35
    CO2(g)    100
KINETICS 1
Forsterite
    formula  Mg2SiO4  1
    m        1
    m0       1
    parms    6.639 1
    tol      1e-08
Magnesite
    formula  MgCO3  1
    m        0
    m0       0
    parms    1 0.5
    tol      1e-08
steps       2592000 in 500 steps # seconds
step_divide 1
runge_kutta 3
bad_step_max 500
INCREMENTAL_REACTIONS True
KNOBS
    iterations            3000
    convergence_tolerance 1e-10
    tolerance             1e-15
    step_size             1000
    pe_step_size          10
USER_GRAPH 1
    -axis_titles            "Time_Hours" "Delta_Forsterite" ""
    -chart_title            "Delta Forsterite"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3600
20 GRAPH_Y KIN_Delta("Forsterite")
  -end
    -active                 true
USER_GRAPH 2
    -axis_titles            "Time_Hours" "Mg Concentration (mol/kgw)" ""
    -chart_title            "Magnesium leaching"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3600
20 GRAPH_Y TOT("Mg")
  -end
    -active                 true
USER_GRAPH 3
    -axis_titles            "Time_Hours" "Moles" ""
    -chart_title            "Forsterite moles"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3600
20 GRAPH_Y KIN("Forsterite")
  -end
    -active                 true
SELECTED_OUTPUT 1
    file                 C:\Users\hp\Desktop\Phreeqc\Brucite\selected_output_P(35-100).sel
    reset                false
USER_PUNCH 1
    headings Time pH Mg_conc. SI_Magnesite Delta_brucite Brucite_mol Magnesite_mol
    start
10 PUNCH TOTAL_TIME/3600
20 PUNCH -LA("H+")
30 PUNCH TOT("Mg")
40 PUNCH SI("Magnesite")
50 PUNCH KIN_DELTA("Forsterite")
60 PUNCH KIN("Forsterite")
70 PUNCH KIN("Magnesite")
    end
END
ERROR: Elements in species have not been tabulated, SiO2.
ERROR: Reaction for species has not been defined, SiO2.
ERROR: Calculations terminating due to input errors
but I'm getting error. I don't know what I'm missing...Please guide me with this

[dlparkhurst]

The two species SiO2(aq) and H4SiO4(aq) represent the same species; the two differ by two waters of hydration. You should only use one of them.

If you use an llnl.dat-style database, the master species for Si is SiO2(aq), so H4SiO4 should not be used.

All other databases that include Si use H4SiO4 as the master species. In these databases SiO2(aq) should not be used.

To convert an equation from one to the other, use the following:

Code: [Select]

SiO2 + 2H2O = H4SiO4, log K 0.0


[Divya]

Sorry to bother you. I did what you said but it is showing the same error.

[dlparkhurst]

Code: [Select]

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

RATES
    Forsterite
-start
 10 rem unit should be mol,kgw-1 and second-1
 20 rem parm(1) is surface area in the unit of m2/kgw
 30 rem calculation of surface area can be found in the note
 40 rem M is current moles of minerals M0 is the initial moles of minerals
 50 rem parm(2) is a scaling factor
 60 rem acid solution parameters
 70 a1=8.38E+04
 80 E1=67206
 90 n1=0.470
100 rem neutral solution parameters
110 a2=1.58E+03
120 E2=79000
130 rem base solution parameters
140 a3=1.00E-07
150 E3=56637
160 n2=-0.600
170 rem rate=0 if no minerals and undersaturated
180 SR_mineral=SR("forsterite")
190 if (M<0) then goto 310
200 if (M=0 and SR_mineral<1) then goto 310
210 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
220 if (SA<=0) then SA=1
230 R=8.31451
240 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1
250 Rate2=a2*EXP(-E2/R/TK)
260 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2
270 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
280 moles= rate*Time
290 rem do not dissolve more minerals than present
300 if (moles>M) then moles=M
310 save moles
-end
    Magnesite
-start
 10 rem unit should be mol,kgw-1 and second-1
 20 rem parm(1) is surface area in the unit of m2/kgw
 30 rem calculation of surface area can be found in the note
 40 rem M is current moles of minerals M0 is the initial moles of minerals
 50 rem parm(2) is a scaling factor
 60 rem acid solution parameters
 70 a1=1.39E-4
 80 E1=14400
 90 n1=1.000
100 rem neutral solution parameters
110 a2=5.99E-6
120 E2=23500
130 rem CO2 denpendence parameters
140 a3=6.03E+05
150 E3=62800
160 n2=1.000
170 rem rate=0 if no minerals and undersaturated
180 SR_mineral=SR("magnesite")
190 if (M<0) then goto 310
200 if (M=0 and SR_mineral<1) then goto 310
210 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
220 if (SA<=0) then SA=1
230 R=8.31451
240 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1
250 Rate2=a2*EXP(-E2/R/TK)
260 Rate3=a3*EXP(-E3/R/TK)*SR("CO2(g)")^n2
270 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral^4)*SA*parm(2)
280 moles= rate*Time
290 rem do not dissolve more minerals than present
300 if (moles>M) then moles=M
310 save moles
-end
PHASES
Forsterite
    Mg2SiO4 + 4H+ = H4SiO4 + 2Mg+2
    log_k     27.8626
    delta_h   -520 kcal
    -analytical_expression -76.195 -0.014013 14763 25.09 -303790 0
    -Vm       32.952 cm3/mol
Magnesite
    MgCO3 + H+ = HCO3- + Mg+2
    log_k     2.2936
    delta_h   -265.63 kJ
    -analytical_expression -166.65 -0.049469 6434.4 65.506 100.45 0
    -Vm       28.018 cm3/mol
END
SOLUTION 1
    temp      35
    pH        7
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    -water    1 # kg
GAS_PHASE 1
    -fixed_pressure
    -pressure 100
    -volume 1
    -temperature 35
    CO2(g)    100
KINETICS 1
Forsterite
    -formula  Mg2SiO4  1
    -m        1
    -m0       1
    -parms    6.639 1
    -tol      1e-08
Magnesite
    -formula  MgCO3  1
    -m        0
    -m0       0
    -parms    1 0.5
    -tol      1e-08
-steps       2592000 in 500 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500
INCREMENTAL_REACTIONS True

USER_GRAPH 1
    -axis_titles            "Time_Hours" "Delta_Forsterite" ""
    -chart_title            "Delta Forsterite"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3600
20 GRAPH_Y KIN_Delta("Forsterite")
  -end
    -active                 true
USER_GRAPH 2
    -axis_titles            "Time_Hours" "Mg Concentration (mol/kgw)" ""
    -chart_title            "Magnesium leaching"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3600
20 GRAPH_Y TOT("Mg")
  -end
    -active                 true
USER_GRAPH 3
    -axis_titles            "Time_Hours" "Moles" ""
    -chart_title            "Forsterite moles"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME/3600
20 GRAPH_Y KIN("Forsterite")
  -end
    -active                 true
END


[Divya]

Sorry ...But getting same error again and again

[dlparkhurst]

The file I gave runs for me with version 3.8.2 (https://github.com/usgs-coupled) and either database phreeqc.dat or phreeqc_rates.dat.

If it does not run for you, post your input and output file.

[dlparkhurst]

And it runs faster if you use -cvode in the KINETICS definition.