Kinetic Dissolution after Initial Equilibrium in Reactive Transport

[Daehyun]

Dear PHREEQC developers and users,

I am trying to set up a long-term reactive transport model to simulate bentonite buffer behavior in a deep geological repository system. My conceptual model includes:

External groundwater inflow (SOLUTION 0).

Bentonite domain (SOLUTION 1-10) initially saturated with groundwater, equilibrated with gases and minerals (including GAS_PHASE and EQUILIBRIUM_PHASES).

During transport, I want minerals to dissolve/precipitate kinetically, not by continuous equilibrium, while diffusion progresses (TRANSPORT with diffusion_only).

I attempted the following approach:

First, I calculated equilibrium of SOLUTION 0 with GAS_PHASE and EQUILIBRIUM_PHASES.

Then I copied the equilibrated solution into SOLUTION 1-10.

I defined KINETICS for minerals, and activated TRANSPORT for 10 cells with diffusion_only flow.

However, I noticed that during TRANSPORT, even with KINETICS defined, PHREEQC keeps calculating EQUILIBRIUM_PHASES (if defined), and essentially overrides the kinetic reactions at every transport shift. As a result, I always get purely equilibrium results, rather than kinetic mineral transformations.

How can I properly set up initial gas-solid-solution equilibrium, but allow only KINETICS to control mineral transformations during TRANSPORT?

I would greatly appreciate any clarification or example models that reflect this kind of long-term coupled system.

Thank you in advance.

Code: [Select]

SOLUTION_MASTER_SPECIES
Ntg Ntg 0 Ntg 28.0134 # N2 gas

SOLUTION_SPECIES
Ntg = Ntg # N2
-dw 1.96e-9
-Vm 7 # Pray et al., 1952, IEC 44. 1146

PHASES
Ntg(g)
    Ntg = Ntg
    log_k     -3.1864
    delta_h   -10.4391 kJ
    -analytical_expression -58.453 0.001818 3199 17.909 -27460 0
    -T_c      126.2
    -P_c      33.5
    -Omega    0.039
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
O2(g)
    O2 = O2
    log_k     -2.8983
    -analytical_expression -7.5001 0.0078981 0 0 200270 0
    -T_c      154.6
    -P_c      49.8
    -Omega    0.021
H2S(g)
    H2S = H+ + HS-
    log_k     -7.93
    delta_h   9.1 kJ
    -analytical_expression -45.07 -0.02418 0 17.9205 0 0
    -T_c      373.2
    -P_c      88.2
    -Omega    0.1
CH4(g)
    CH4 = CH4
    log_k     -2.8
    -analytical_expression 10.44 -0.00765 -6669 0 1014000 0
    -T_c      196.6
    -P_c      45.4
    -Omega    0.008

RATES
Montmor-Ca
10 rem acid solution parameters
11 a1=1.94984E-13
12 E1=48000
13 n1=0.220
20 a2=3.89045E-15
21 E2=48000
30 rem base solution parameters
31 a3=3.89045E-15
32 E3=48000
33 n2=-0.130
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Montmor-Ca")
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 
85 Rate3=a3*EXP(-E3/R/TK)*ACT("OH-")^n2         
90 Rate=(Rate1+ Rate3)*(1-Sr_mineral)*SA*0.1
100 moles= rate*Time
200 save moles
-end


Cristobalite(alpha)
-start
20 rem neutral solution parameters
21 a2=5.88844E-13
22 E2=74500
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Cristobalite(alpha)")
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
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
90 Rate=(Rate2)*(1-Sr_mineral)*SA*0.1
100 moles= rate*Time
200 save moles
-end

Albite
-start
10 rem acid solution parameters
11 a1=6.91831E-11
12 E1=65000
13 n1=0.457
20 rem neutral solution parameters
21 a2=2.75423E-13
22 E2=69800
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Albite")
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
90 Rate=(Rate1+Rate2)*(1-Sr_mineral)*SA*0.1
100 moles= rate*Time
200 save moles
-end

Quartz
-start
20 rem neutral solution parameters
21 a2=3.98107E-14
22 E2=90900
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Quartz")
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
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
90 Rate=(Rate2)*(1-Sr_mineral)*SA*0.1
100 moles= rate*Time
200 save moles
-end


Calcite
-start
10 rem acid solution parameters
11 a1=0.501187
12 E1=14400
13 n1=1
20 rem neutral solution parameters
21 a2=1.54882E-6
22 E2=23500
30 rem basic dependence parameters
31 a3=0.000331
32 E3=35400
33 n2=1
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Calcite")
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("OH")^n2    #base rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*0.1
100 moles= rate*Time
200 save moles
-end


Organic_C
 -start
1   REM      Additive Monod kinetics for SOC (sediment organic carbon)
2   REM      Electron acceptors: O2, NO3, and SO4
10  if (M <= 0) THEN GOTO 200
20  mO2   = MOL("O2")
40  mSO4  = TOT("S(6)")
50  k_O2  = 1.04167E-9    # 1/sec   
70  k_SO4 = 1.e-13     # 1/sec
80  rate  = k_O2 * mO2/(3e-04 + mO2)
100 rate  = rate + k_SO4 * mSO4/(1e-4 + mSO4)
110 moles = rate * M * (M/M0) * TIME
200 SAVE moles
 -end

Pyrite
-start
30 rem Neutral dependence parameters
31 a2=2.81838E-05
32 E2=56900
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Pyrite")
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
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
90 Rate=Rate2*(1-Sr_mineral)*SA*0.1
100 moles= rate*Time
200 save moles
-end


SOLUTION 0
    temp      14.8
    pH        9.05
    pe        -7.665
    redox     pe
    units     mol/kgw
    density   1
    Al        2.733e-06
    C         0.0013 as HCO3
    Ca        0.0001425
    Cl(-1)    0.001474
    F         0.0004272
    Fe        1.24e-07
    K         8.456e-06
    Mg        1.195e-05
    N(5)      1.103e-05
    Na        0.001649
    S         6.05e-05 as SO4
    Si        0.0001251 as SiO2
    Sr        2.664e-06
    -water    1 # kg

SAVE SOLUTION 0

END

SOLUTION 1-10
    temp      14.8
    pH        9.05
    pe        -7.665
    redox     pe
    units     mol/kgw
    density   1
    Al        2.733e-06
    C         0.0013 as HCO3
    Ca        0.0001425
    Cl(-1)    0.001474
    F         0.0004272
    Fe        1.24e-07
    K         8.456e-06
    Mg        1.195e-05
    N(5)      1.103e-05
    Na        0.001649
    S         6.05e-05 as SO4
    Si        0.0001251 as SiO2
    Sr        2.664e-06
    O(0)      0
    -water    1 # kg

GAS_PHASE 1-10
    -fixed_pressure
    -pressure 1
    -volume 1
    -temperature 25
    CO2(g)    0.0004
    Ntg(g)    0.79
    O2(g)     0.2
    CH4(g)    0

EQUILIBRIUM_PHASES 1-10
    Albite    0 0.533924
    Calcite   0 0.098932
    Cristobalite(alpha) 0 2.164
    Dolomite  0 0
    Montmor-Ca 0 1.912036
    Pyrite    0 0
    Quartz    0 0.333
    Beidellite-Ca 0 0

Save solution 1-10


REACTION_TEMPERATURE 1
    60

REACTION_PRESSURE 1
    50

END


KINETICS 1-10
Organic_C
    -formula  CH2O  1
    -m        1
    -m0       1
    -tol      1e-08
Montmor-Ca
    -formula  Ca0.165Al2.33Si3.67O10(OH)2  1
    -m        1.912
    -m0       1.912
    -parms    1
    -tol      1e-08
Quartz
    -formula  SiO2  1
    -m        0.333
    -m0       0.333
    -parms    1
    -tol      1e-08
Albite
    -formula  NaAlSi3O8  1
    -m        0.5338
    -m0       0.5338
    -parms    0.01
    -tol      1e-08
Calcite
    -formula  CaCO3  1
    -m        0.09905
    -m0       0.09905
    -parms    1
    -tol      1e-08
Cristobalite(alpha)
    -formula  SiO2  1
    -m        2.164
    -m0       2.164
    -parms    1
    -tol      1e-08

-steps       1
-step_divide 1
-runge_kutta 3
-bad_step_max 1000


TRANSPORT
    -cells                 10
    -shifts                1
    -time_step             31536000000 # seconds
    -flow_direction        diffusion_only
    -boundary_conditions   constant closed
    -lengths               10*0.1
    -diffusion_coefficient 1e-10
    -print_frequency       1
    -punch_frequency       1


[dlparkhurst]

Once you define EQUILIBRIUM_PHASES 1-10, they remain in the cells for the calculation. You have a couple options.

(1) Use an equilibrium phases definition with a user number that a number above 10, like 100. Then no cell numbered 1-10 will have an equilibrium phase definition.

(2) Before TRANSPORT, delete the equilibrium phases:

Code: [Select]

END
DELETE
-equilibrium_phases 1-10
END
TRANSPORT
...

The GAS_PHASE defiinitions will also be present in the cells. You may want to delete them as well. I'm not sure it makes sense to have gas phases for an aqueous diffusion calculation.