Kinetic reaction between groundwater and bentonite

[Daehyun]

Hello,

I am a student researching how bentonite dissolves over time and which secondary minerals precipitate when bentonite and the surrounding groundwater become saturated, following the environmental evolution of a deep geological repository for spent nuclear fuel. Unfortunately, there's nobody around me who I can ask about PHREEQC, so I am writing this post.

I would like to model a simplified reaction scenario initially:

1. Groundwater reaches equilibrium with the CO2 gas present within the repository, then saturates the bentonite.
2. I want to understand how the saturated bentonite dissolves over time and how to input the potential precipitation of secondary minerals.
I look forward to a helpful response. Thank you!


SOLUTION 1 #groundwater composition
    temp      30
    pH        8.985
    pe        -5.976
    units     mol/kgw
    density   1
    Al        2.733e-06
    C(4)      0.001302
    Ca        0.0001425
    Cl(-1)    0.001474
    F         0.0004272
    K         8.456e-06
    Mg        1.195e-05
    Na        0.001649
    Si        0.0001251 as SiO2
    Sr        2.664e-06
    Fe(2)     1.24e-07
    N(5)      1.103e-05
    U(4)      6.335e-09
    S(-2)     6.07e-05
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    CO2(g)    -3.5 0.00012645277459486
END
USE SOLUTION 1 #Bentonite composition
EQUILIBRIUM_PHASES 2
    Albite    0 0.053389576
    Calcite   0 0.009991318
    Cristobalite(alpha) 0 0.216362677
    Montmor-Ca 0 0.191233768
    Quartz    0 0.033286566

RATES

Montmor-Ca
-start
1 eh = 48
2 en = 48
3 eb = 48
4 logkh = -12.71
5 logkn = -14.41
6 logkb = -14.41
7 nh = 0.22
8 nb = -0.13
13 C = 7184.940964
15 R = 8.3144
19 dif_T = 1/TK - 1/303.15
20 pk_H = logkh - (eh*1000)/(2.3025*R)*dif_T
30 rate_H = 10^pk_H*act("H+")^nh
40 pk_w = logkn - (en*1000)/(2.3025*R)*dif_T
50 rate_w = 10^pk_w
60 pk_OH = logkb - (eb*1000)/(2.3025*R)*dif_T
70 rate_OH = 10^pk_OH*act("H+")^nb
80 rate0 = rate_H + rate_w + rate_OH
90 rate = rate0*C*parm(2)*parm(1)*(m/m0)^0.67*(1-SR("Montmor-Ca"))
100 moles = rate*time
110 if SR("Montmor-Ca") > 1 then moles = 0.1*moles
120 save moles
-end

Albite
-start
1 eh = 65
2 en = 69.8
3 eb = 71
4 logkh = -10.16
5 logkn = -12.56
6 logkb = -15.60
7 nh = 0.46
8 nb = -0.572
13 C = 9991.495488
15 R = 8.3144
19 dif_T = 1/TK - 1/303.15
20 pk_H = logkh - (eh*1000)/(2.3025*R)*dif_T
30 rate_H = 10^pk_H*act("H+")^nh
40 pk_w = logkn - (en*1000)/(2.3025*R)*dif_T
50 rate_w = 10^pk_w
60 pk_OH = logkb - (eb*1000)/(2.3025*R)*dif_T
70 rate_OH = 10^pk_OH*act("H+")^nb
80 rate0 = rate_H + rate_w + rate_OH
90 rate = C*rate0*parm(2)*parm(1)*(m/m0)^0.67*(1-SR("Albite"))
100 moles = rate*time
110 if SR("Albite") > 1 then moles = 0.1*moles
120 save moles
-end

Calcite
-start
1 eh = 14.40
2 en = 23.50
3 eb = 35.40
4 logkh = -0.3
5 logkn = -5.8
6 logkb = -3.48
7 nh = 1
8 nb = 1
13 C = 27086.461865
15 R = 8.3144
19 dif_T = 1/TK - 1/303.15
20 pk_H = logkh - (eh*1000)/(2.3025*R)*dif_T
30 rate_H = 10^pk_H*act("H+")^nh
40 pk_w = logkn - (en*1000)/(2.3025*R)*dif_T
50 rate_w = 10^pk_w
60 pk_OH = logkb - (eb*1000)/(2.3025*R)*dif_T
70 rate_OH = 10^pk_OH*act("H+")^nb
80 rate0 = rate_H + rate_w + rate_OH
90 rate = C*rate0*parm(2)*parm(1)*(m/m0)^0.67*(1-SR("Calcite"))
100 moles = rate*time
110 if SR("Calcite") > 1 then moles = 0.1*moles
120 save moles
-end

Quartz
-start
1 eh = 90.9
4 logkh = -13.40
13 C = 44104.699564
15 R = 8.3144
19 dif_T = 1/TK - 1/303.15
20 pk_H = logkh - (eh*1000)/(2.3025*R)*dif_T
30 rate_H = 10^pk_H*act("H+")^nh
80 rate0 = rate_H
90 rate = C*rate0*parm(2)*parm(1)*(m/m0)^0.67*(1-SR("Quartz"))
100 moles = rate*time
110 if SR("Quartz") > 1 then moles = 0.1*moles
120 save moles
-end

Cristobalite(alpha)
-start
1 eh = 65
4 logkh = -12.31
13 C = 43272.535421
15 R = 8.3144
19 dif_T = 1/TK - 1/303.15
20 pk_H = logkh - (eh*1000)/(2.3025*R)*dif_T
30 rate_H = 10^pk_H*act("H+")^nh
80 rate0 = rate_H
90 rate = rate0*C*parm(2)*parm(1)*(m/m0)^0.67*(1-SR("Cristobalite(alpha)"))
110 moles = rate*time
120 if SR("Cristobalite(alpha)") > 1 then moles = 0.1*moles
130 save moles
-end




KINETICS 1
Montmor-Ca
    -formula  Ca0.165Al2.33Si3.67O10(OH)2  1
    -m        0.191233768
    -m0       0.191233768
    -parms    56.38 366.04414
    -tol      1e-07
Quartz
    -formula  SiO2  1
    -m        0.033286566
    -m0       0.033286566
    -parms    0.00226 60.0843
    -tol      1e-07
Albite
    -formula  NaAlSi3O8  1
    -m        0.053389576
    -m0       0.053389576
    -parms    0.00229 262.22347
    -tol      1e-07
Calcite
    -formula  CaCO3  1
    -m        0.009991318
    -m0       0.009991318
    -parms    0.00221 100.0869
    -tol      1e-07
Cristobalite(alpha)
    -formula  SiO2  1
    -m        0.216362677
    -m0       0.216362677
    -parms    0.00264 100.0869
    -tol      1e-06


-steps       31536000000 in 365 steps # seconds 1000 years
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5


SELECTED_OUTPUT 1
    -file                 BENTONITE1.sel
    -pH                   true         
    -Equilibrium phases   Quartz  Albite  Montmor-Ca  Cristobalite(alpha) Calcite
    -kinetic_reactants    Quartz  Albite  Montmor-Ca  Cristobalite(alpha) Calcite

[dlparkhurst]

I don't see a question in your post.

In your first calculation, you need to use SAVE if your want to replace the original SOLUTION 1 with the reacted solution 1. The following breaks the calculation into the two calculations that occur--initial solution, and reaction.

Code: [Select]

SOLUTION 1 #groundwater composition
    temp      30
    pH        8.985
    pe        -5.976
    units     mol/kgw
    density   1
    Al        2.733e-06
    C(4)      0.001302
    Ca        0.0001425
    Cl(-1)    0.001474
    F         0.0004272
    K         8.456e-06
    Mg        1.195e-05
    Na        0.001649
    Si        0.0001251 as SiO2
    Sr        2.664e-06
    Fe(2)     1.24e-07
    N(5)      1.103e-05
    U(4)      6.335e-09
    S(-2)     6.07e-05
    -water    1 # kg
END
USE solution 1
EQUILIBRIUM_PHASES 1
    CO2(g)    -3.5 0.00012645277459486
SAVE solution 1
END


[Daehyun]

Hello Dr.,

Thank you for your advice. I apologize for the lack of content in my question. My question is that I want to model the precipitation of secondary minerals in my desired scenario, but when I add minerals that can precipitate in the Equilibrium_Phases, the model takes too long to run, and the concentrations of elements show negative moles. I believe this is happening because montmorillonite, which should dissolve, continues to precipitate. Do you think my understanding is correct?

Additionally, I'm uploading the modified phreeqc.

Thank you.

SOLUTION 1
    temp      30
    pH        8.985
    pe        -5.976
    units     mol/kgw
    density   1
    Al        2.733e-06
    C(4)      0.001302
    Ca        0.0001425
    Cl(-1)    0.001474
    F         0.0004272
    K         8.456e-06
    Mg        1.195e-05
    Na        0.001649
    Si        0.0001251 as SiO2
    Sr        2.664e-06
    Fe(2)     1.24e-07
    N(5)      1.103e-05
    U(4)      6.335e-09
    S(-2)     6.07e-05
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    CO2(g)    -1 0.04024792723174756
SAVE SOLUTION 2
END
USE SOLUTION 2
EQUILIBRIUM_PHASES 2
Chalcedony 0 0
Clinoptilolite-Ca   0 0
Clinoptilolite-Na   0 0
Dolomite 0 0
Montmor-Na  0 0     

RATES

Montmor-Ca
-start
1 eh = 22.67728
2 en = 22.67728
3 eb = 22.67728
4 logkh = -12.67
5 logkn = -13.67
6 logkb = -14.67
7 nh = 0.5
8 nb = -0.5
13 C = 7184.940964
15 R = 8.3144
19 dif_T = 1/TK - 1/303.15
20 pk_H = logkh - (eh*1000)/(2.3025*R)*dif_T
25  if (M <= 0 and SR("Montmor-Ca") < 0) then goto 120
30 rate_H = 10^pk_H*act("H+")^nh
40 pk_w = logkn - (en*1000)/(2.3025*R)*dif_T
50 rate_w = 10^pk_w
60 pk_OH = logkb - (eb*1000)/(2.3025*R)*dif_T
70 rate_OH = 10^pk_OH*act("H+")^nb
80 rate0 = rate_H + rate_w + rate_OH
90 rate = rate0*C*parm(2)*parm(1)*(m/m0)^0.67*(1-SR("Montmor-Ca"))
100 moles = rate*time
110 if SR("Montmor-Ca") > 1 then moles = 0.1*moles
120 save moles
-end

Albite
-start
1 eh = 65
2 en = 69.8
3 eb = 71
4 logkh = -10.16
5 logkn = -12.56
6 logkb = -15.60
7 nh = 0.46
8 nb = -0.572
13 C = 9991.495488
15 R = 8.3144
19 dif_T = 1/TK - 1/303.15
20 pk_H = logkh - (eh*1000)/(2.3025*R)*dif_T
25  if (M <= 0 and SR("Albite") < 0) then goto 120
30 rate_H = 10^pk_H*act("H+")^nh
40 pk_w = logkn - (en*1000)/(2.3025*R)*dif_T
50 rate_w = 10^pk_w
60 pk_OH = logkb - (eb*1000)/(2.3025*R)*dif_T
70 rate_OH = 10^pk_OH*act("H+")^nb
80 rate0 = rate_H + rate_w + rate_OH
90 rate = C*rate0*parm(2)*parm(1)*(m/m0)^0.67*(1-SR("Albite"))
100 moles = rate*time
110 if SR("Albite") > 1 then moles = 0.1*moles
120 save moles
-end

Calcite
-start
1 eh = 14.40
2 en = 23.50
3 eb = 35.40
4 logkh = -0.3
5 logkn = -5.8
6 logkb = -3.48
7 nh = 1
8 nb = 1
13 C = 27086.461865
15 R = 8.3144
19 dif_T = 1/TK - 1/303.15
20 pk_H = logkh - (eh*1000)/(2.3025*R)*dif_T
25  if (M <= 0 and SR("Calcite") < 0) then goto 120
30 rate_H = 10^pk_H*act("H+")^nh
40 pk_w = logkn - (en*1000)/(2.3025*R)*dif_T
50 rate_w = 10^pk_w
60 pk_OH = logkb - (eb*1000)/(2.3025*R)*dif_T
70 rate_OH = 10^pk_OH*act("H+")^nb
80 rate0 = rate_H + rate_w + rate_OH
90 rate = C*rate0*parm(2)*parm(1)*(m/m0)^0.67*(1-SR("Calcite"))
100 moles = rate*time
110 if SR("Calcite") > 1 then moles = 0.1*moles
120 save moles
-end

Quartz
-start
1 eh = 90.9
4 logkh = -13.40
13 C = 44104.699564
15 R = 8.3144
19 dif_T = 1/TK - 1/303.15
20 pk_H = logkh - (eh*1000)/(2.3025*R)*dif_T
25  if (M <= 0 and SR("Quartz") < 0) then goto 120
30 rate_H = 10^pk_H*act("H+")^nh
80 rate0 = rate_H
90 rate = C*rate0*parm(2)*parm(1)*(m/m0)^0.67*(1-SR("Quartz"))
100 moles = rate*time
110 if SR("Quartz") > 1 then moles = 0.1*moles
120 save moles
-end

Cristobalite(alpha)
-start
1 eh = 65
4 logkh = -12.31
13 C = 43272.535421
15 R = 8.3144
19 dif_T = 1/TK - 1/303.15
20 pk_H = logkh - (eh*1000)/(2.3025*R)*dif_T
25  if (M <= 0 and SR("Cristobalite(alpha)") < 0) then goto 130
30 rate_H = 10^pk_H*act("H+")^nh
80 rate0 = rate_H
90 rate = rate0*C*parm(2)*parm(1)*(m/m0)^0.67*(1-SR("Cristobalite(alpha)"))
110 moles = rate*time
120 if SR("Cristobalite(alpha)") > 1 then moles = 0.1*moles
130 save moles
-end



KINETICS 1
Montmor-Ca
    -formula  Ca0.165Al2.33Si3.67O10(OH)2  1
    -m        0.191233768
    -m0       0.191233768
    -parms    56.38 366.04414
    -tol      1e-07
Quartz
    -formula  SiO2  1
    -m        0.033286566
    -m0       0.033286566
    -parms    0.00226 60.0843
    -tol      1e-07
Albite
    -formula  NaAlSi3O8  1
    -m        0.053389576
    -m0       0.053389576
    -parms    0.00229 262.22347
    -tol      1e-07
Calcite
    -formula  CaCO3  1
    -m        0.009991318
    -m0       0.009991318
    -parms    0.00221 100.0869
    -tol      1e-07
Cristobalite(alpha)
    -formula  SiO2  1
    -m        0.216362677
    -m0       0.216362677
    -parms    0.00264 100.0869
    -tol      1e-06


-steps       7776000 in 90 steps # seconds 90D years
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5
INCREMENTAL_REACTIONS True

SELECTED_OUTPUT 1
    -file                 Sal(F).sel
    -pH                   true
    -equilibrium_phases    Chalcedony Clinoptilolite-Ca  Clinoptilolite-Na  Dolomite  Montmor-Na   
    -saturation_indices   Albite  Montmor-Ca  Calcite  Quartz
                          Cristobalite(alpha)
    -kinetic_reactants    Quartz  Albite  Montmor-Ca  Cristobalite(alpha)  Calcite

[dlparkhurst]

Your rates for Calcite and Montmor-Ca are fast, equilibrating in 1 second. If these rates are correct, move Montmor-Ca and Calcite to EQUILIBRIUM_PHASES and the simulation will run much faster.

You have three SiO2 phases. Do you really need all three?

Code: [Select]

RATES

Montmor-Ca
-start
1 eh = 22.67728
2 en = 22.67728
3 eb = 22.67728
4 logkh = -12.67
5 logkn = -13.67
6 logkb = -14.67
7 nh = 0.5
8 nb = -0.5
13 C = 7184.940964
15 R = 8.3144
19 dif_T = 1/TK - 1/303.15
20 pk_H = logkh - (eh*1000)/(2.3025*R)*dif_T
25  if (M <= 0 and SR("Montmor-Ca") < 0) then goto 120
30 rate_H = 10^pk_H*act("H+")^nh
40 pk_w = logkn - (en*1000)/(2.3025*R)*dif_T
50 rate_w = 10^pk_w
60 pk_OH = logkb - (eb*1000)/(2.3025*R)*dif_T
70 rate_OH = 10^pk_OH*act("H+")^nb
80 rate0 = rate_H + rate_w + rate_OH
90 rate = rate0*C*parm(2)*parm(1)*(m/m0)^0.67*(1-SR("Montmor-Ca"))
100 moles = rate*time
#110 if SR("Montmor-Ca") > 1 then moles = 0.1*moles
120 save moles
-end

Albite
-start
1 eh = 65
2 en = 69.8
3 eb = 71
4 logkh = -10.16
5 logkn = -12.56
6 logkb = -15.60
7 nh = 0.46
8 nb = -0.572
13 C = 9991.495488
15 R = 8.3144
19 dif_T = 1/TK - 1/303.15
20 pk_H = logkh - (eh*1000)/(2.3025*R)*dif_T
25  if (M <= 0 and SR("Albite") < 0) then goto 120
30 rate_H = 10^pk_H*act("H+")^nh
40 pk_w = logkn - (en*1000)/(2.3025*R)*dif_T
50 rate_w = 10^pk_w
60 pk_OH = logkb - (eb*1000)/(2.3025*R)*dif_T
70 rate_OH = 10^pk_OH*act("H+")^nb
80 rate0 = rate_H + rate_w + rate_OH
90 rate = C*rate0*parm(2)*parm(1)*(m/m0)^0.67*(1-SR("Albite"))
100 moles = rate*time
#110 if SR("Albite") > 1 then moles = 0.1*moles
120 save moles
-end

Calcite
-start
1 eh = 14.40
2 en = 23.50
3 eb = 35.40
4 logkh = -0.3
5 logkn = -5.8
6 logkb = -3.48
7 nh = 1
8 nb = 1
13 C = 27086.461865
15 R = 8.3144
19 dif_T = 1/TK - 1/303.15
20 pk_H = logkh - (eh*1000)/(2.3025*R)*dif_T
25  if (M <= 0 and SR("Calcite") < 0) then goto 120
30 rate_H = 10^pk_H*act("H+")^nh
40 pk_w = logkn - (en*1000)/(2.3025*R)*dif_T
50 rate_w = 10^pk_w
60 pk_OH = logkb - (eb*1000)/(2.3025*R)*dif_T
70 rate_OH = 10^pk_OH*act("H+")^nb
80 rate0 = rate_H + rate_w + rate_OH
90 rate = C*rate0*parm(2)*parm(1)*(m/m0)^0.67*(1-SR("Calcite"))
100 moles = rate*time
#110 if SR("Calcite") > 1 then moles = 0.1*moles
120 save moles
-end

Quartz
-start
1 eh = 90.9
4 logkh = -13.40
13 C = 44104.699564
15 R = 8.3144
19 dif_T = 1/TK - 1/303.15
20 pk_H = logkh - (eh*1000)/(2.3025*R)*dif_T
25  if (M <= 0 and SR("Quartz") < 0) then goto 120
30 rate_H = 10^pk_H*act("H+")^nh
80 rate0 = rate_H
90 rate = C*rate0*parm(2)*parm(1)*(m/m0)^0.67*(1-SR("Quartz"))
100 moles = rate*time
#110 if SR("Quartz") > 1 then moles = 0.1*moles
120 save moles
-end

Cristobalite(alpha)
-start
1 eh = 65
4 logkh = -12.31
13 C = 43272.535421
15 R = 8.3144
19 dif_T = 1/TK - 1/303.15
20 pk_H = logkh - (eh*1000)/(2.3025*R)*dif_T
25  if (M <= 0 and SR("Cristobalite(alpha)") < 0) then goto 130
30 rate_H = 10^pk_H*act("H+")^nh
80 rate0 = rate_H
90 rate = rate0*C*parm(2)*parm(1)*(m/m0)^0.67*(1-SR("Cristobalite(alpha)"))
110 moles = rate*time
#120 if SR("Cristobalite(alpha)") > 1 then moles = 0.1*moles
130 save moles
-end
END
SOLUTION 1
    temp      30
    pH        8.985
    pe        -5.976
    units     mol/kgw
    density   1
    Al        2.733e-06
    C(4)      0.001302
    Ca        0.0001425
    Cl(-1)    0.001474
    F         0.0004272
    K         8.456e-06
    Mg        1.195e-05
    Na        0.001649
    Si        0.0001251 as SiO2
    #Sr        2.664e-06
    Fe(2)     1.24e-07
    #N(5)      1.103e-05
    #U(4)      6.335e-09
    S(-2)     6.07e-05
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    CO2(g)    -1 0.04024792723174756
SAVE SOLUTION 2
END
USE SOLUTION 2

EQUILIBRIUM_PHASES 2
Chalcedony 0 0
Clinoptilolite-Ca   0 0
Clinoptilolite-Na   0 0
Dolomite 0 0
Montmor-Na  0 0 
Montmor-Ca  0 0   
Calcite 0 0.009991318
KINETICS 1
#Montmor-Ca
#    -formula  Ca0.165Al2.33Si3.67O10(OH)2  1
#    -m        0.191233768
#    -m0       0.191233768
#    -parms    56.38 366.04414
    #-tol      1e-07
Quartz
    -formula  SiO2  1
    -m        0.033286566
    -m0       0.033286566
    -parms    0.00226 60.0843
    #-tol      1e-07
Albite
    -formula  NaAlSi3O8  1
    -m        0.053389576
    -m0       0.053389576
    -parms    0.00229 262.22347
    #-tol      1e-07
#Calcite
#    -formula  CaCO3  1
#    -m        0.009991318
#    -m0       0.009991318
#    -parms    0.00221 100.0869
#    -tol      1e-07
Cristobalite(alpha)
    -formula  SiO2  1
    -m        0.216362677
    -m0       0.216362677
    -parms    0.00264 100.0869
    #-tol      1e-06

-steps    1 3 1e1 3e1 1e2 3e2 1e3 3e3 1e4 3e4 1e5 3e5 1e6 3e6 1e7 #  7776000 in 90 steps # seconds 90D years
-step_divide 1
-runge_kutta 3
-bad_step_max 500
-cvode true
-cvode_steps 100
-cvode_order 5
INCREMENTAL_REACTIONS True
USER_GRAPH 1
    -headings               time Montmor-Ca Quartz Albite Cristobalite
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y KIN_DELTA("Montmor-Ca"), KIN_DELTA("Quartz")
30 GRAPH_Y KIN_DELTA("Albite"), KIN_DELTA("Cristobalite(alpha)")
  -end
    -active                 true
END

[Daehyun]

Thank you, Doctor!
Thanks to your advice, I believe I can achieve some of the desired results.