Reactive Transport and Diffusion

[dat]

Hi,

I am new to Reactive transport modelling.
I am planning to do a reactive transport model of hydrogen in reservoir rock and in cap rock. For example, I developed a code for two rock systems of 20 cells. Gas will be initially injected into cells 11-20 and then it will diffuse to cells 1-10.

I want to find the diffusion of hydrogen and mineral distribution and produced gasses in each cell at a given time.

Is it possible to find the above two? Please comment on the code below.
Really appreciate your help.

Code: [Select]

DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.7.3-15968\database\core10.dat
RATES
Quartz
   -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Quartz") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Quartz") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0                 
    70 k_neut = 1.02e-14*EXP((-87700/8.3145)*(1/TK-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-SR("Quartz"))
    190 moles = r * TIME
    200 SAVE moles
     -end

K-Feldspar
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("K-Feldspar") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("K-Feldspar") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 8.71e-11*EXP((-51700/8.3145)*(1/Tk-1/298.15))                 
    70 k_neut = 3.89e-13*EXP((-38000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 6.31e-22*EXP((-94100/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.5 + k_neut + k_base*act("H+")^(-0.823)
    100 r = k_rateconst * SA * (1-SR("K-Feldspar"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Kaolinite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Kaolinite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Kaolinite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 7.7e-14*EXP((-65900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 7.7e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
    80 k_base = 4.70e-18*EXP((-17900/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.777 + k_neut + k_base*act("H+")^(-0.472)
    100 r = k_rateconst * SA * (1-SR("Kaolinite"))
    190 moles = r * TIME
    200 SAVE moles
    -end


SOLUTION 1-10
    temp      60
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1.021
    Cl        17000 mg/L
    Na        8600 mg/L
    Ca        900 mg/L
    Mg        110 mg/L
    K         230 mg/L
    C(4)      440 mg/L
    pressure  100
    -water    1 # kg

 
INCREMENTAL_REACTIONS True

KINETICS 1-10
Quartz
    -formula  SiO2  1
    -m        402.12
    -m0       402.12
    -parms    4.122
    -tol      1e-08

K-Feldspar
    -formula  KAlSi3O8  1
    -m        2.08
    -m0       2.08
    -parms    19.768
    -tol      1e-08

SOLUTION 11-20
    temp      80
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1.021
    Cl        17000 mg/L
    Na        8600 mg/L
    Ca        900 mg/L
    Mg        110 mg/L
    K         230 mg/L
    C(4)      440 mg/L
    pressure  197.4
    -water    1 # kg

 
INCREMENTAL_REACTIONS True

GAS_PHASE 11-20
    -fixed_pressure
    -pressure 197.385
    -volume 1.222
    -temperature 80
    H2(g)     197.385
    CO2(g)   0

KINETICS 11-20
Quartz
    -formula  SiO2  1
    -m        402.12
    -m0       402.12
    -parms    4.122
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        1.46
    -m0       1.46
    -parms    299.027
    -tol      1e-08
K-Feldspar
    -formula  KAlSi3O8  1
    -m        2.08
    -m0       2.08
    -parms    19.768
    -tol      1e-08

-steps       31.536 3153.6 3153.6 31536 31536 31536 31536 31536 315360 315360 3153600 3153600 3153600 3153600 3153600 3153600 3153600 31536000 31536000 31536000 31536000 31536000 31536000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 103280400
-step_divide 1
-runge_kutta 3
-bad_step_max 4000
-cvode true
-cvode_steps 500
-cvode_order 5


TRANSPORT
        -cells                20               
        -length               0.005               
        -shifts               336               
        -time_step            3600             
        -flow_direction       forward
        -boundary_condition   flux flux
        -dispersivity         .05
        -correct_disp         true

END


[dlparkhurst]

Because your kinetic reactions involve non-redox minerals, the diffusion of H2 will have little effect on the rates. You can experiment with and without the GAS_PHASE, but I think the results will be very similar.

PHREEQC allows only advective-dispersive transport or diffusive transport as specified by the -flow_direction identifier in the TRANSPORT block. If you want a diffusive calculation, specify -flow_direction diffusion_only. However, you also need a database that has -dw defined for aqueous species in SOLUTION_SPECIES. Phreeqc.dat, Amm.dat, pitzer.dat, and frezchem.dat have -dw definitions.

There is a -multi_D option in TRANSPORT that allows each aqueous species to have its own diffusion coefficient, so that different species diffuse at different rates, generally subject to no net charge transport (no current).

I don't think flux boundary conditions are allowed in diffusion calculations; you will have to experiment.

The definition of time steps in TRANSPORT overrides the time steps defined in KINETICS.

Code: [Select]

RATES
Quartz
   -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Quartz") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Quartz") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0                 
    70 k_neut = 1.02e-14*EXP((-87700/8.3145)*(1/TK-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-SR("Quartz"))
    190 moles = r * TIME
    200 SAVE moles
     -end

K-Feldspar
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("K-Feldspar") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("K-Feldspar") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 8.71e-11*EXP((-51700/8.3145)*(1/Tk-1/298.15))                 
    70 k_neut = 3.89e-13*EXP((-38000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 6.31e-22*EXP((-94100/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.5 + k_neut + k_base*act("H+")^(-0.823)
    100 r = k_rateconst * SA * (1-SR("K-Feldspar"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Kaolinite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Kaolinite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Kaolinite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 7.7e-14*EXP((-65900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 7.7e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
    80 k_base = 4.70e-18*EXP((-17900/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.777 + k_neut + k_base*act("H+")^(-0.472)
    100 r = k_rateconst * SA * (1-SR("Kaolinite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
END
SOLUTION 1-10
    temp      60
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1.021
    Cl        17000 mg/L
    Na        8600 mg/L
    Ca        900 mg/L
    Mg        110 mg/L
    K         230 mg/L
    C(4)      440 mg/L
    pressure  100
    -water    1 # kg
END

KINETICS 1-10
Quartz
    -formula  SiO2  1
    -m        402.12
    -m0       402.12
    -parms    4.122
    -tol      1e-08

K-Feldspar
    -formula  KAlSi3O8  1
    -m        2.08
    -m0       2.08
    -parms    19.768
    -tol      1e-08
END
SOLUTION 11-20
    temp      80
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1.021
    Cl        17000 mg/L
    Na        8600 mg/L
    Ca        900 mg/L
    Mg        110 mg/L
    K         230 mg/L
    C(4)      440 mg/L
    pressure  197.4
    -water    1 # kg

END

GAS_PHASE 11-20
    -fixed_pressure
    -pressure 197.385
    -volume 1.222
    -temperature 80
    H2(g)     197.385
    CO2(g)   0
END
KINETICS 11-20
Quartz
    -formula  SiO2  1
    -m        402.12
    -m0       402.12
    -parms    4.122
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        1.46
    -m0       1.46
    -parms    299.027
    -tol      1e-08
K-Feldspar
    -formula  KAlSi3O8  1
    -m        2.08
    -m0       2.08
    -parms    19.768
    -tol      1e-08

-steps       31.536 3153.6 3153.6 31536 31536 31536 31536 31536 315360 315360 3153600 3153600 3153600 3153600 3153600 3153600 3153600 31536000 31536000 31536000 31536000 31536000 31536000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 315360000 103280400
-step_divide 1
-runge_kutta 3
-bad_step_max 4000
-cvode true
-cvode_steps 500
-cvode_order 5
END

TRANSPORT
        -cells                20               
        -length               0.005               
        -shifts               100 #336               
        -time_step            3600             
        -flow_direction       diffusion_only
        -boundary_condition   flux flux
  -punch_frequency      100
  -punch_cells          1-20
  -print_frequency      10
        #-dispersivity         .05
        #-correct_disp         true
USER_GRAPH 1
    -headings               dist H2(aq)
    -axis_titles            "Distance, meters" "H2(aq), Molality" "Si, molality"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X DIST
20 GRAPH_Y MOL("H2")
30 GRAPH_SY TOT("Si")
  -end
    -active                 true

END

[dat]

Hi,

Thank you for the comments. I have modified my code.  But I am getting - cvode time step error. Could you please help?

What I wanted to do is I want to define a tracer with the same properties as H2(g)/(aq) and find that diffusivity from that tracer and also with H2(g). Normally hydrogen will inject into cells 183-850 and is supposed to diffuse to 1-182.

I have entered kinetic parameters for both rocks, and hydrogen will be at 182-850.

Here is my full code. Could you please check whether my idea is in the code and why it is giving an error?

Code: [Select]

RATES
Quartz
   -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Quartz") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Quartz") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0                 
    70 k_neut = 1.02e-14*EXP((-87700/8.3145)*(1/TK-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-SR("Quartz"))
    190 moles = r * TIME
    200 SAVE moles
     -end

K-Feldspar
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("K-Feldspar") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("K-Feldspar") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 8.71e-11*EXP((-51700/8.3145)*(1/Tk-1/298.15))                 
    70 k_neut = 3.89e-13*EXP((-38000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 6.31e-22*EXP((-94100/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.5 + k_neut + k_base*act("H+")^(-0.823)
    100 r = k_rateconst * SA * (1-SR("K-Feldspar"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Kaolinite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Kaolinite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Kaolinite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 7.7e-14*EXP((-65900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 7.7e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
    80 k_base = 4.70e-18*EXP((-17900/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.777 + k_neut + k_base*act("H+")^(-0.472)
    100 r = k_rateconst * SA * (1-SR("Kaolinite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Muscovite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Muscovite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Muscovite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 1.41e-12*EXP((-22000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-14*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.82e-15*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.37 + k_neut + k_base*act("H+")^(-0.22)
    100 r = k_rateconst * SA * (1-SR("Muscovite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Albite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Albite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Albite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 6.92e-11*EXP((-65000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.75e-13*EXP((-69800/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.51e-16*EXP((-71000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.457 + k_neut + k_base*act("H+")^(-0.572)
    100 r = k_rateconst * SA * (1-SR("Albite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Pyrite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Pyrite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Pyrite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 3.02e-8*EXP((-56900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-5*EXP((-56900/8.3145)*(1/Tk-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid*act("H+")^(-0.5)  + k_neut*act("O2")^0.5  + k_base*act("H+")^0
    100 r = k_rateconst * SA * (1-SR("Pyrite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Calcite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Calcite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Calcite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0.5012*EXP(-14.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")                 
    70 k_neut = 1.5488e-6*EXP(-23.5e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 3.31e-4*EXP(-35.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("CO2")
    90 k_rateconst = k_acid*act("H+")^1  + k_neut + k_base*act("H+")^1
    100 r = k_rateconst * SA * (1-SR("Calcite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

SOLUTION_MASTER_SPECIES

Tracer Tracer 0.0 1.0 1.0

EXCHANGE_MASTER_SPECIES
X X-

SOLUTION_SPECIES
Tracer = Tracer
-log_k 0.0
#-delta_h -4.184 kJ
#-analytic -9.3114 4.6473e-3 -49.335 1.4341 1.2815e5
#-T_c 33.2
#-P_c 12.80
#-Omega -0.225
-dw 1.0e-9


EXCHANGE_SPECIES
X- = X-
-log_k 0.0

Na+ + X- = NaX
-log_k 0.0
-gamma 4.08 0.082

K+ + X- = KX
-log_k 0.7
-gamma 3.5 0.015
-delta_h  -4.3 # Jardine & Sparks, 1984

Li+ + X- = LiX
-log_k -0.08
-gamma 6.0 0
-delta_h  1.4 # Merriam & Thomas, 1956

# !!!!!
# H+ + X- = HX
# -log_k 1.0
# -gamma 9.0 0


Ca+2 + 2X- = CaX2
-log_k 0.8
-gamma 5.0 0.165
-delta_h  7.2    # Van Bladel & Gheyl, 1980

Mg+2 + 2X- = MgX2
-log_k 0.6
-gamma 5.5 0.2
-delta_h  7.4 # Laudelout et al., 1968


Mn+2 + 2X- = MnX2
-log_k 0.52
-gamma 6.0 0

Fe+2 + 2X- = FeX2
-log_k 0.44
-gamma 6.0 0

Al+3 + 3X- = AlX3
-log_k 0.41
-gamma 9.0 0

AlOH+2 + 2X- = AlOHX2
-log_k 0.89
-gamma 0.0 0


Selected_Output
-file my1_UHS.xls
-temperature true
-totals C(4) C(-4) Tracer H(0) S(6) S(-2) Na Cl N N(-3)
#-kinetic_phases K-feldspar Albite Kaolinite Quartz Calcite Pyrite Muscovite   
-gases CO2(g) CH4(g) H2S(g) H2(g) N2(g)
-water
-charge_balance true
-ionic_strength true


SOLUTION 1-182
    temp      60
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1.021
    Cl        17000 mg/L
    Na        8600 mg/L
    Ca        900 mg/L
    Mg        110 mg/L
    K         230 mg/L
    C(4)      440 mg/L
    pressure  100
    -water    1 # kg

REACTION_TEMPERATURE 1-182
60

REACTION_PRESSURE 1-182
100

GAS_PHASE 1-182
-fixed_pressure
-pressure 0
-temperature 60
CO2(g) 0.0
CH4(g) 0.0
H2S(g) 0.0


KINETICS 1-182

Quartz
    -formula  SiO2  1
    -m        409.387
    -m0       409.387
    -parms    4.122
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        21.52
    -m0       21.52
    -parms    299.027
    -tol      1e-08
K-Feldspar
    -formula  KAlSi3O8  1
    -m        20.32
    -m0       20.32
    -parms    19.768
    -tol      1e-08
Muscovite     
    -formula  KAl3Si3O10(OH)2  1
    -m        11.4
    -m0       11.4
    -parms    25.607
    -tol      1e-08
Pyrite       
    -formula  FeS2 1
    -m        62.3
    -m0       62.3
    -parms    4.354
    -tol      1e-08
Calcite       
    -formula  CaCO3  1
    -m        27.248
    -m0       27.248
    -parms    6.715
    -tol      1e-08

EXCHANGE 1-182
-equilibrate with solution 1-182
X 3.907
END


SOLUTION 183-850 # reservoir rock
    temp      80
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1.021
    Cl        17000 mg/L
    Na        8600 mg/L
    Ca        900 mg/L
    Mg        110 mg/L
    K         230 mg/L
    C(4)      440 mg/L
    pressure  197.4
    -water    1 # kg
    Tracer 0.001


GAS_PHASE 183-850
-fixed_pressure
-pressure 197.4
-temperature 80.0
-volume 3.0
H2(g) 184.56
CH4(g) 3.56
N2(g) 0.4
CO2(g) 1.48
H2S(g) 0.0

REACTION_PRESSURE 183-850
80.0
REACTION_TEMPERATURE 183-850
80.0

EXCHANGE 183-850
-equilibrate with solution 183-850
X 3.722

INCREMENTAL_REACTIONS True

KINETICS 183-850
Quartz
    -formula  SiO2  1
    -m        402.12
    -m0       402.12
    -parms    4.122
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        1.46
    -m0       1.46
    -parms    299.027
    -tol      1e-08
K-Feldspar
    -formula  KAlSi3O8  1
    -m        2.08
    -m0       2.08
    -parms    19.768
    -tol      1e-08

-steps       94608000 94608000 94608000 94608000 94608000 94608000 94608000 94608000 94608000 94608000
-step_divide 1
-runge_kutta 3
-bad_step_max 4000
-cvode true
-cvode_steps 500
-cvode_order 5

END

TRANSPORT
-cells 850
-shifts 10
-flow_direction diffusion_only
-time_step 94608000
-multi_d true
-length 1488*1.0
-dispersivity 1488*0.0
-boundary_conditions flux flux

END




And also, I want to get the results into Excel. The code I have entered for that is not working. Could you please check?

Code: [Select]

-kinetic_phases K-feldspar Albite Kaolinite Quartz Calcite Pyrite Muscovite 


[dlparkhurst]

Go back to a 20 cell model and I will look at it.

[dat]

Hi here is the 20 cell model. Could you please check?

Code: [Select]

DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.7.3-15968\database\core10.dat
RATES
Quartz
   -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Quartz") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Quartz") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0                 
    70 k_neut = 1.02e-14*EXP((-87700/8.3145)*(1/TK-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-SR("Quartz"))
    190 moles = r * TIME
    200 SAVE moles
     -end

K-Feldspar
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("K-Feldspar") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("K-Feldspar") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 8.71e-11*EXP((-51700/8.3145)*(1/Tk-1/298.15))                 
    70 k_neut = 3.89e-13*EXP((-38000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 6.31e-22*EXP((-94100/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.5 + k_neut + k_base*act("H+")^(-0.823)
    100 r = k_rateconst * SA * (1-SR("K-Feldspar"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Kaolinite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Kaolinite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Kaolinite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 7.7e-14*EXP((-65900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 7.7e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
    80 k_base = 4.70e-18*EXP((-17900/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.777 + k_neut + k_base*act("H+")^(-0.472)
    100 r = k_rateconst * SA * (1-SR("Kaolinite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Muscovite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Muscovite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Muscovite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 1.41e-12*EXP((-22000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-14*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.82e-15*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.37 + k_neut + k_base*act("H+")^(-0.22)
    100 r = k_rateconst * SA * (1-SR("Muscovite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Albite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Albite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Albite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 6.92e-11*EXP((-65000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.75e-13*EXP((-69800/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.51e-16*EXP((-71000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.457 + k_neut + k_base*act("H+")^(-0.572)
    100 r = k_rateconst * SA * (1-SR("Albite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Pyrite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Pyrite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Pyrite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 3.02e-8*EXP((-56900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-5*EXP((-56900/8.3145)*(1/Tk-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid*act("H+")^(-0.5)  + k_neut*act("O2")^0.5  + k_base*act("H+")^0
    100 r = k_rateconst * SA * (1-SR("Pyrite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Calcite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Calcite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Calcite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0.5012*EXP(-14.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")                 
    70 k_neut = 1.5488e-6*EXP(-23.5e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 3.31e-4*EXP(-35.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("CO2")
    90 k_rateconst = k_acid*act("H+")^1  + k_neut + k_base*act("H+")^1
    100 r = k_rateconst * SA * (1-SR("Calcite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

SOLUTION_MASTER_SPECIES

Tracer Tracer 0.0 1.0 1.0

EXCHANGE_MASTER_SPECIES
X X-

SOLUTION_SPECIES
Tracer = Tracer
-log_k 0.0
#-delta_h -4.184 kJ
#-analytic -9.3114 4.6473e-3 -49.335 1.4341 1.2815e5
#-T_c 33.2
#-P_c 12.80
#-Omega -0.225
-dw 1.0e-9


EXCHANGE_SPECIES
X- = X-
-log_k 0.0

Na+ + X- = NaX
-log_k 0.0
-gamma 4.08 0.082

K+ + X- = KX
-log_k 0.7
-gamma 3.5 0.015
-delta_h  -4.3 # Jardine & Sparks, 1984

Li+ + X- = LiX
-log_k -0.08
-gamma 6.0 0
-delta_h  1.4 # Merriam & Thomas, 1956

# !!!!!
# H+ + X- = HX
# -log_k 1.0
# -gamma 9.0 0


Ca+2 + 2X- = CaX2
-log_k 0.8
-gamma 5.0 0.165
-delta_h  7.2    # Van Bladel & Gheyl, 1980

Mg+2 + 2X- = MgX2
-log_k 0.6
-gamma 5.5 0.2
-delta_h  7.4 # Laudelout et al., 1968


Mn+2 + 2X- = MnX2
-log_k 0.52
-gamma 6.0 0

Fe+2 + 2X- = FeX2
-log_k 0.44
-gamma 6.0 0

Al+3 + 3X- = AlX3
-log_k 0.41
-gamma 9.0 0

AlOH+2 + 2X- = AlOHX2
-log_k 0.89
-gamma 0.0 0


Selected_Output
-file my1_UHS.xls
-temperature true
-totals C(4) C(-4) Tracer H(0) S(6) S(-2) Na Cl N N(-3)
#-kinetic_phases K-feldspar Albite Kaolinite Quartz Calcite Pyrite Muscovite   
-gases CO2(g) CH4(g) H2S(g) H2(g) N2(g)
-water
-charge_balance true
-ionic_strength true


SOLUTION 1-10
    temp      60
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1.021
    Cl        17000 mg/L
    Na        8600 mg/L
    Ca        900 mg/L
    Mg        110 mg/L
    K         230 mg/L
    C(4)      440 mg/L
    pressure  100
    -water    1 # kg

REACTION_TEMPERATURE 1-10
60

REACTION_PRESSURE 1-10
100

GAS_PHASE 1-10
-fixed_pressure
-pressure 0
-temperature 60
CO2(g) 0.0
CH4(g) 0.0
H2S(g) 0.0


KINETICS 1-10

Quartz
    -formula  SiO2  1
    -m        409.387
    -m0       409.387
    -parms    4.122
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        21.52
    -m0       21.52
    -parms    299.027
    -tol      1e-08
K-Feldspar
    -formula  KAlSi3O8  1
    -m        20.32
    -m0       20.32
    -parms    19.768
    -tol      1e-08
Muscovite     
    -formula  KAl3Si3O10(OH)2  1
    -m        11.4
    -m0       11.4
    -parms    25.607
    -tol      1e-08
Pyrite       
    -formula  FeS2 1
    -m        62.3
    -m0       62.3
    -parms    4.354
    -tol      1e-08
Calcite       
    -formula  CaCO3  1
    -m        27.248
    -m0       27.248
    -parms    6.715
    -tol      1e-08

EXCHANGE 1-10
-equilibrate with solution 1-10
X 3.907
END


SOLUTION 11-20 # reservoir rock
    temp      80
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1.021
    Cl        17000 mg/L
    Na        8600 mg/L
    Ca        900 mg/L
    Mg        110 mg/L
    K         230 mg/L
    C(4)      440 mg/L
    pressure  197.4
    -water    1 # kg
    Tracer 0.001


GAS_PHASE 11-20
-fixed_pressure
-pressure 197.4
-temperature 80.0
-volume 1.222
H2(g) 184.56
CH4(g) 3.56
N2(g) 0.4
CO2(g) 1.48
H2S(g) 0.0

REACTION_PRESSURE 11-20
80.0
REACTION_TEMPERATURE 11-20
80.0

EXCHANGE 11-20
-equilibrate with solution 11-20
X 3.722

INCREMENTAL_REACTIONS True

KINETICS 11-20
Quartz
    -formula  SiO2  1
    -m        402.12
    -m0       402.12
    -parms    4.122
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        1.46
    -m0       1.46
    -parms    299.027
    -tol      1e-08
K-Feldspar
    -formula  KAlSi3O8  1
    -m        2.08
    -m0       2.08
    -parms    19.768
    -tol      1e-08

-steps       94608000 94608000 94608000 94608000 94608000 94608000 94608000 94608000 94608000 94608000
-step_divide 1
-runge_kutta 3
-bad_step_max 4000
-cvode true
-cvode_steps 500
-cvode_order 5

END

TRANSPORT
-cells 20
-shifts 10
-flow_direction diffusion_only
-time_step 94608000
-multi_d true
-length 1488*1.0
-dispersivity 1488*0.0
-boundary_conditions flux flux

END




[dlparkhurst]

It does not look like you even read my previous post, nor put any effort into your script.

Study this.

Code: [Select]

DATABASE Amm.dat
PHASES
Muscovite
# KAl3Si3O10(OH)2 + 10 H+ = K+ + 3 Al+3 + 3 SiO2 + 6 H2O
KAl3Si3O10(OH)2 + 10 H+ = K+ + 3 Al+3 + 3 H4SiO4
log_k 13.5858
-delta_H -243.224 kJ/mol
# deltafH -1427.41 kcal/mol
-analytic 3.3085e1 -1.2425e-2 1.2477e4 -2.0865e1 -5.4692e5
# Range 0-350
-Vm 140.71
# Extrapol supcrt92
# Ref HDN+78
END
RATES
Quartz
   -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Quartz") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Quartz") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0                 
    70 k_neut = 1.02e-14*EXP((-87700/8.3145)*(1/TK-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-SR("Quartz"))
    190 moles = r * TIME
    200 SAVE moles
     -end

K-Feldspar
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("K-Feldspar") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("K-Feldspar") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 8.71e-11*EXP((-51700/8.3145)*(1/Tk-1/298.15))                 
    70 k_neut = 3.89e-13*EXP((-38000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 6.31e-22*EXP((-94100/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.5 + k_neut + k_base*act("H+")^(-0.823)
    100 r = k_rateconst * SA * (1-SR("K-Feldspar"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Kaolinite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Kaolinite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Kaolinite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 7.7e-14*EXP((-65900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 7.7e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
    80 k_base = 4.70e-18*EXP((-17900/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.777 + k_neut + k_base*act("H+")^(-0.472)
    100 r = k_rateconst * SA * (1-SR("Kaolinite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Muscovite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Muscovite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Muscovite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 1.41e-12*EXP((-22000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-14*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.82e-15*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.37 + k_neut + k_base*act("H+")^(-0.22)
    100 r = k_rateconst * SA * (1-SR("Muscovite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Albite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Albite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Albite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 6.92e-11*EXP((-65000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.75e-13*EXP((-69800/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.51e-16*EXP((-71000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.457 + k_neut + k_base*act("H+")^(-0.572)
    100 r = k_rateconst * SA * (1-SR("Albite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Pyrite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Pyrite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Pyrite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 3.02e-8*EXP((-56900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-5*EXP((-56900/8.3145)*(1/Tk-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid*act("H+")^(-0.5)  + k_neut*act("O2")^0.5  + k_base*act("H+")^0
    100 r = k_rateconst * SA * (1-SR("Pyrite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Calcite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Calcite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Calcite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0.5012*EXP(-14.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")                 
    70 k_neut = 1.5488e-6*EXP(-23.5e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 3.31e-4*EXP(-35.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("CO2")
    90 k_rateconst = k_acid*act("H+")^1  + k_neut + k_base*act("H+")^1
    100 r = k_rateconst * SA * (1-SR("Calcite"))
    190 moles = r * TIME
    200 SAVE moles
    -end


EXCHANGE_MASTER_SPECIES
X X-


EXCHANGE_SPECIES
X- = X-
-log_k 0.0

Na+ + X- = NaX
-log_k 0.0
-gamma 4.08 0.082

K+ + X- = KX
-log_k 0.7
-gamma 3.5 0.015
-delta_h  -4.3 # Jardine & Sparks, 1984

Li+ + X- = LiX
-log_k -0.08
-gamma 6.0 0
-delta_h  1.4 # Merriam & Thomas, 1956

# !!!!!
# H+ + X- = HX
# -log_k 1.0
# -gamma 9.0 0


Ca+2 + 2X- = CaX2
-log_k 0.8
-gamma 5.0 0.165
-delta_h  7.2    # Van Bladel & Gheyl, 1980

Mg+2 + 2X- = MgX2
-log_k 0.6
-gamma 5.5 0.2
-delta_h  7.4 # Laudelout et al., 1968


Mn+2 + 2X- = MnX2
-log_k 0.52
-gamma 6.0 0

Fe+2 + 2X- = FeX2
-log_k 0.44
-gamma 6.0 0

Al+3 + 3X- = AlX3
-log_k 0.41
-gamma 9.0 0

AlOH+2 + 2X- = AlOHX2
-log_k 0.89
-gamma 0.0 0
END
SOLUTION 1-10
    temp      60
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1.021
    Cl        17000 mg/L
    Na        8600 mg/L
    Ca        900 mg/L
    Mg        110 mg/L
    K         230 mg/L
    C(4)      440 mg/L
S(6)   1 mg/L
    pressure  100
    -water    1 # kg
END
REACTION_TEMPERATURE 1-10
60
END
REACTION_PRESSURE 1-10
100
END
GAS_PHASE 1-10
-fixed_pressure
-pressure 0
-temperature 60
CO2(g) 0.0
CH4(g) 0.0
H2S(g) 0.0
END
EXCHANGE 1-10
-equilibrate with solution 1
X 3.907
END
RUN_CELL
-cell 1-10
END

KINETICS 1-10
Quartz
    -formula  SiO2  1
    -m        409.387
    -m0       409.387
    -parms    4.122
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        21.52
    -m0       21.52
    -parms    299.027
    -tol      1e-08
K-Feldspar
    -formula  KAlSi3O8  1
    -m        20.32
    -m0       20.32
    -parms    19.768
    -tol      1e-08
Muscovite     
    -formula  KAl3Si3O10(OH)2  1
    -m        11.4
    -m0       11.4
    -parms    25.607
    -tol      1e-08
Pyrite       
    -formula  FeS2 1
    -m        62.3
    -m0       62.3
    -parms    4.354
    -tol      1e-08
Calcite       
    -formula  CaCO3  1
    -m        27.248
    -m0       27.248
    -parms    6.715
    -tol      1e-08
-cvode
END

SOLUTION 11-20 # reservoir rock
    temp      80
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1.021
    Cl        17000 mg/L
    Na        8600 mg/L
    Ca        900 mg/L
    Mg        110 mg/L
    K         230 mg/L
    C(4)      440 mg/L
S(6)   1 mg/L
    pressure  197.4
    -water    1 # kg
    Hdg       0.001
END
GAS_PHASE 11-20
-fixed_pressure
-pressure 197.4
-temperature 80.0
-volume 1.222
H2(g) 184.56
#Hdg(g) 184.56
CH4(g) 3.56
N2(g) 0.4
CO2(g) 1.48
H2S(g) 0.0
END
REACTION_PRESSURE 11-20
80.0
END
REACTION_TEMPERATURE 11-20
80.0
END
EXCHANGE 11-20
-equilibrate with solution 11-20
X 3.722
END
RUN_CELLS
-cell 11-20
END

KINETICS 11-20
Quartz
    -formula  SiO2  1
    -m        402.12
    -m0       402.12
    -parms    4.122
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        1.46
    -m0       1.46
    -parms    299.027
    -tol      1e-08
K-Feldspar
    -formula  KAlSi3O8  1
    -m        2.08
    -m0       2.08
    -parms    19.768
    -tol      1e-08

#-steps       94608000 94608000 94608000 94608000 94608000 94608000 94608000 94608000 94608000 94608000
#-step_divide 1
#-runge_kutta 3
#-bad_step_max 4000
-cvode true
-cvode_steps 500
-cvode_order 5

END

TRANSPORT
-cells 20
-shifts 10
-flow_direction diffusion_only
-time_step 94608000
-multi_d true
-length 1.0
-dispersivity 0.0
-punch_cells      1-20
-punch_frequency  10
-boundary_conditions flux flux
USER_GRAPH
10 GRAPH_X dist
20 GRAPH_Y MOL("H2")
30 GRAPH_SY TOT("Hdg")
END



[dat]

I have looked at the code more. Using Hdg to find the diffusion makes sense (It has a diffusion coefficient).

But If I want to find the diffusion of H(0)-Reactive hydrogen, Is there a way to define that in the solution species with the same diffusion coefficient? I got an error when defining that. I am planning to use core10.dat.

Thank you for the help.

Code: [Select]

SOLUTION_MASTER_SPECIES
Hdg Hdg 0 Hdg 2.016 # H2 gas
H H+ -1 H 1.0079
H(0) H2 0 H
H(+1) H+ -1 0

SOLUTION_SPECIES
Hdg = Hdg # H2
-dw 5.13e-9
-Vm 6.52  0.78  0.12 # supcrt

H2 = H2
log_k -3.1050
-delta_H -4.184 kJ/mol
-analytic -9.3114 4.6473e-3 -49.335 1.4341 1.2815e5
-dw 5.13e-9


PHASES
Hdg(g)
Hdg = Hdg
-analytic   -9.3114    4.6473e-3   -49.335    1.4341    1.2815e5
-T_c  33.2 ; -P_c   12.80 ; -Omega -0.225
H2(g)
H2 = H2
log_k -3.1050
-delta_H -4.184 kJ/mol
# deltafH 0 kcal/mol
-analytic -9.3114 4.6473e-3 -49.335 1.4341 1.2815e5
# Range 0-350
-T_c 33.2 # K
-P_c 12.80
-Omega 0.225 # phreeqc.dat
# Extrapol supcrt92
# Ref WEP+82, Kel60

RATES
Quartz
   -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Quartz") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Quartz") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0                 
    70 k_neut = 1.02e-14*EXP((-87700/8.3145)*(1/TK-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-SR("Quartz"))
    190 moles = r * TIME
    200 SAVE moles
     -end

K-Feldspar
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("K-Feldspar") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("K-Feldspar") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 8.71e-11*EXP((-51700/8.3145)*(1/Tk-1/298.15))                 
    70 k_neut = 3.89e-13*EXP((-38000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 6.31e-22*EXP((-94100/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.5 + k_neut + k_base*act("H+")^(-0.823)
    100 r = k_rateconst * SA * (1-SR("K-Feldspar"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Kaolinite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Kaolinite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Kaolinite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 7.7e-14*EXP((-65900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 7.7e-14*EXP((-22200/8.3145)*(1/Tk-1/298.15))
    80 k_base = 4.70e-18*EXP((-17900/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.777 + k_neut + k_base*act("H+")^(-0.472)
    100 r = k_rateconst * SA * (1-SR("Kaolinite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Muscovite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Muscovite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Muscovite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 1.41e-12*EXP((-22000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-14*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.82e-15*EXP((-22000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.37 + k_neut + k_base*act("H+")^(-0.22)
    100 r = k_rateconst * SA * (1-SR("Muscovite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Albite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Albite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Albite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 6.92e-11*EXP((-65000/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.75e-13*EXP((-69800/8.3145)*(1/Tk-1/298.15))
    80 k_base = 2.51e-16*EXP((-71000/8.3145)*(1/Tk-1/298.15))
    90 k_rateconst = k_acid*act("H+")^0.457 + k_neut + k_base*act("H+")^(-0.572)
    100 r = k_rateconst * SA * (1-SR("Albite"))
    190 moles = r * TIME
    200 SAVE moles
    -end

Pyrite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Pyrite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Pyrite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 3.02e-8*EXP((-56900/8.3145)*(1/Tk-1/298.15))               
    70 k_neut = 2.82e-5*EXP((-56900/8.3145)*(1/Tk-1/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid*act("H+")^(-0.5)  + k_neut*act("O2")^0.5  + k_base*act("H+")^0
    100 r = k_rateconst * SA * (1-SR("Pyrite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
Calcite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    30 if (M <= 0 and SR("Calcite") < 1) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and SR("Calcite") > 1) then SA = 1e-05  #nucleation
    60 k_acid = 0.5012*EXP(-14.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")                 
    70 k_neut = 1.5488e-6*EXP(-23.5e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 3.31e-4*EXP(-35.4e+03/8.314*(1.0/TK-1.0/298.15))*ACT("CO2")
    90 k_rateconst = k_acid*act("H+")^1  + k_neut + k_base*act("H+")^1
    100 r = k_rateconst * SA * (1-SR("Calcite"))
    190 moles = r * TIME
    200 SAVE moles
    -end
END
EXCHANGE_MASTER_SPECIES
X X-
EXCHANGE_SPECIES
X- = X-
-log_k 0.0

Na+ + X- = NaX
-log_k 0.0
-gamma 4.08 0.082

K+ + X- = KX
-log_k 0.7
-gamma 3.5 0.015
-delta_h  -4.3 # Jardine & Sparks, 1984

Li+ + X- = LiX
-log_k -0.08
-gamma 6.0 0
-delta_h  1.4 # Merriam & Thomas, 1956

# !!!!!
# H+ + X- = HX
# -log_k 1.0
# -gamma 9.0 0


Ca+2 + 2X- = CaX2
-log_k 0.8
-gamma 5.0 0.165
-delta_h  7.2    # Van Bladel & Gheyl, 1980

Mg+2 + 2X- = MgX2
-log_k 0.6
-gamma 5.5 0.2
-delta_h  7.4 # Laudelout et al., 1968


Mn+2 + 2X- = MnX2
-log_k 0.52
-gamma 6.0 0

Fe+2 + 2X- = FeX2
-log_k 0.44
-gamma 6.0 0

Al+3 + 3X- = AlX3
-log_k 0.41
-gamma 9.0 0

AlOH+2 + 2X- = AlOHX2
-log_k 0.89
-gamma 0.0 0

Selected_Output
-file Equi,H2,Hdg with modified database.xls
-temperature true
-totals C(4) C(-4) H(0) Hdg H2 S(6) Na Cl N N(-3)
-equilibrium_phases Quartz Kaolinite K-Feldspar Muscovite Pyrite Calcite
#-kinetics H2
-gases CO2(g) CH4(g) H2S(g) H2(g) N2(g)
-water
-charge_balance true
-ionic_strength true
#USER_PUNCH
     #-headings Time(year) KIN(Quartz) KIN_DELTA(Quartz)KIN(Kaolinite) KIN_DELTA(Kaolinite)KIN(K-Feldspar) KIN_DELTA(K-Feldspar)KIN(Muscovite) KIN_DELTA(Muscovite)KIN(Pyrite) KIN_DELTA(Pyrite)KIN(Calcite) KIN_DELTA(Calcite)
     #-start
     #10 punch TOTAL_TIME/31536000 KIN("Quartz") KIN_DELTA("Quartz")KIN("Kaolinite") KIN_DELTA("Kaolinite")KIN("K-Feldspar") KIN_DELTA("K-Feldspar") KIN("Muscovite") KIN_DELTA("Muscovite")KIN("Pyrite") KIN_DELTA("Pyrite")KIN("Calcite") KIN_DELTA("Calcite")
     #-end
END

SOLUTION 1-40
    temp      60
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1.021
    Cl        17000 mg/L
    Na        8600 mg/L
    Ca        900 mg/L
    Mg        110 mg/L
    K         230 mg/L
    C(4)      440 mg/L
    -water    1 # kg

REACTION_TEMPERATURE 1-40
60

REACTION_PRESSURE 1-40
148

EQUILIBRIUM_PHASES 1-40
    Quartz    0 409.387
    Kaolinite 0 21.52
    K-Feldspar 0 20.32
    Muscovite 0  11.4
    Pyrite 0 62.3
    Calcite 0 27.248
    Brucite 0 0
    Diopside 0 0
    Grossular 0 0
    Monticellite 0 0
    Portlandite 0 0
    Troilite 0 0

END


EXCHANGE 1-40
-equilibrate with solution 1-40
X 3.907
END

SOLUTION 41-147 # reservoir rock
    temp      80
    pH        7
    pe        4
    redox     pe
    units     mol/l
    density   1.021
    Cl        17000 mg/L
    Na        8600 mg/L
    Ca        900 mg/L
    Mg        110 mg/L
    K         230 mg/L
    C(4)      440 mg/L
    pressure  197.4
    -water    1 # kg
    Hdg       0.01

REACTION_PRESSURE 41-147
80.0
REACTION_TEMPERATURE 41-147
80.0

GAS_PHASE 41-147
-fixed_pressure
-pressure 197.4
-temperature 80.0
-volume 1.222
H2(g) 197.4

EQUILIBRIUM_PHASES 41-147
    Quartz    0 402.12
    Kaolinite 0 1.46
    K-Feldspar 0 2.08
    Brucite 0 0
    Celadonite 0 0
    Clinochlore-14A 0 0
    Lizardite 0 0
    Mesolite 0 0
    Phlogopite 0 0
    Talc 0 0

EXCHANGE 41-147
-equilibrate with solution 41-147
X 3.722


TRANSPORT
-cells 147
-shifts 100
-flow_direction diffusion_only
-time_step 31536000 # (1 Year)
-multi_d true
-length 147*1.0
-dispersivity 147*0.0
-boundary_conditions flux flux

END






[dlparkhurst]

H2(aq) is a "secondary" master species for H(0). That means that it does not have an identity reaction (H2 = H2), but must be written in terms that can be reduced to "primary" master species. If you look in core10.dat, the equation they use is

Code: [Select]

H2O = H2 + 0.5 O2   
-CO2_llnl_gamma
log_k -46.1066
-delta_H 275.588 kJ/mol
# deltafH -1 kcal/mol
-analytic 6.6835e1 1.7172e-2 -1.8849e4 -2.4092e1 4.2501e5
# Range 0-350
-Vm 5.1427 4.7758 3.8729 -2.9764 -0.209
# Extrapol supcrt92
# Ref SHS89

to which you could add

Code: [Select]

-dw 5.13e-9

You have defined -multi_d, but core10.dat does not have -dw definitions for the aqueous species (SOLUTION_SPECIES). For diffusion calculations, all the aqueous species will be assigned the same default diffusion coefficient, except for your explicit definition of -dw for H2. For the same concentration gradient, all ions would diffuse at the same rate. The calculation is the same as what you would get if you define -diffusion_only without -multi_d (except for H2).

If you use pitzer.dat, Amm.dat, or phreeqc.dat, most aqueous species have -dw defined, in which case -multi_d will have different diffusion coefficients, which allows some ions to diffuse relatively faster or slower; that is, for the same concentration gradient some ions would diffuse faster and some slower.  By default, the combined diffusion of ions will have net zero current.