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