Reactive Transport Modelling involving Kinetic oxidation

[gayathri]

Respected Sir,
Being a starter in PHREEQC, I am having some doubts regarding reactive transport modelling. I am working on geochemical transport modelling for simulating acid mine drainage due to oxidation of pyrite. A column with background solution which is in equilibrium with pyrite(SI=-231 and amount =10 moles) is considered and infilling solution with Equillibrium phase oxygen(SI=-0.67) is taken. Kinetics and rates are defined as per PHREEQC.DAT. Transport datablock containing 100 cells of 5m and dispersivity=0.5 is taken. The simulation is running. But I am unable to obtain the reduction in pH.    I am incorporating my input file with  the mail. Kindly see to it and requesting you to provide with necessary corrections wherever required. Awaiting for the reply.
Thanking You,
Gayathri Krishna K

Input File:   
SOLUTION 1-100 Dakshin Clay Mines,Sasthavattom
    temp      25
    pH        6.4
    pe        4
    redox     pe
    units     mg/l
    density   1
    Ca        3.2
    Cl        16
    F         0.02
    K         0.3
    Mg        2
    N(5)      2.8
    Na        4.8
    S(6)      2.5
    -water    1 # kg
EQUILIBRIUM_PHASES 1-100
    Pyrite    -10 10
SOLUTION 101 Dakshin Clay Mines,Sasthavattom
    temp      25
    pH        6.4
    pe        4
    redox     pe
    units     mg/l
    density   1
    Ca        3.2
    Cl        16
    F         0.02
    K         0.3
    Mg        2
    N(5)      2.8
    Na        4.8
    S(6)      2.5
    -water    1 # kg
 
EQUILIBRIUM_PHASES 101
    O2(g)     -0.67 10
 
SAVE solution 101
RATES
    Pyrite
-start
 10 REM PARM(1)=log10(A/V,1/dm) PARM(2)=exp for (M/M0)
 20 REM PARM(3)=exp for O2 PARM(4)=exp for H+
 30 IF(M <= 0) THEN GOTO 100
 40 IF(SI("Pyrite") >= 0) THEN GOTO 100
 50 lograte = -10.19+PARM(1)+PARM(2)*LOG10(M/M0)
 60 lograte = lograte+PARM(3)*LM("O2")+PARM(4)*LM("H+")
 70 moles = (10^lograte)*TIME
 80 IF(moles > M) THEN moles = M
 90 IF(moles >= MOL("O2")/3.5) THEN moles = MOL("O2")/3.5
100 SAVE moles
-end
KINETICS 1-100
Pyrite
    -formula  FeS2  1
    -m        0.162
    -m0       0.162
    -parms    -5 0.1 0.5 -0.11
    -tol      1e-008
-steps       1 in 1 steps # seconds
-step_divide 1
-runge_kutta 6
-bad_step_max 500
INCREMENTAL_REACTIONS True
USER_GRAPH 1
    -headings               Moles  pH
    -axis_titles            "Time in days" "Millimoles oxidised/kgw" "pH"
    -chart_title            "Kinetic oxidation of pyrite"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 x = TOTAL_TIME/3600/24
20 PLOT_XY x, KIN("Pyrite")*1e3
30 PLOT_XY x, -LA("H+")
  -end
    -active                 true
SELECTED_OUTPUT 1
    -file                 F:\IITD\SEM 4\May\Phreeqc Model\Trials\selected_output_1.sel
    -time                 true
    -pH                   true
    -saturation_indices   Pyrite
    -kinetic_reactants    Pyrite
TRANSPORT
    -cells                 100
    -shifts                3250
    -time_step             86400 # seconds
    -flow_direction        back
    -lengths               100*5
    -dispersivities        100*0.5
    -correct_disp          true
    -print_cells           1 50 100
    -punch_cells           1 50 100
USER_GRAPH 2
    -headings               Moles  pH
    -axis_titles            "Distance" "Millimoles oxidized" "pH"
    -chart_title            "Reactive transport"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 x=DIST
20 PLOT_XY x,KIN("Pyrite")*1e3
30 PLOT_XY x,-LA("H+")
  -end
    -active                 true
SELECTED_OUTPUT 2
    -file                 F:\IITD\SEM 4\May\Phreeqc Model\Trials\selected_output_2.sel
    -distance             true
    -pH                   true
    -totals               Fe(2)
    -saturation_indices   Pyrite  Fe(OH)3(a)
    -kinetic_reactants    Pyrite
    -active               true
    -user_punch           true
 

[dlparkhurst]

I suspect the reaction is just very slow with that set of parameters. You may want to work with a batch reaction first to get a feel for the reaction. Here is a simulation that lowers the pH over the course of about a third of a year by increasing parm(1) by 6 orders of magnitude.

SOLUTION
END
EQUILIBRIUM_PHASES
O2(g) -0.7  10
Goethite 0  0
END
Pyrite
  -start
   1 rem   parm(1) = log10(A/V, 1/dm)   parm(2) = exp for (m/m0)
   2 rem   parm(3) = exp for O2      parm(4) = exp for H+

   10 if (m <= 0) then goto 200
   20 if (si("Pyrite") >= 0) then goto 200
   30  rate = -10.19 + parm(1) + parm(3)*lm("O2") + parm(4)*lm("H+") + parm(2)*log10(m/m0)
   40  moles = 10^rate * time
   50 if (moles > m) then moles = m
   60 if (moles >= (mol("O2")/3.5)) then moles = mol("O2")/3.5
   200 save moles
  -end
END
INCREMENTAL_REACTIONS
KINETICS 1
Pyrite
    -formula  FeS2  1
    -m        0.162
    -m0       0.162
    #-parms    -5 0.1 0.5 -0.11
    -parms    1 0.1 0.5 -0.11
    -tol      1e-008
-steps      1 10 100 1000 10000 1e5 1e6 1e7
-step_divide 1
-runge_kutta 6
-bad_step_max 500
END
USE solution 1
USE equilibrium_phases 1
USE kinetics 1
USER_GRAPH 1
    -headings               time pH
    -axis_titles            "Time, in seconds" "pH" ""
    -chart_title            "PARM(1) = 1, fixed P(O2)"
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 graph_x time
20 graph_y -LA("H+")
  -end
    -active                 true
END

[gayathri]

Thank you for the reply. Can I also have suggestions regarding the transport part. I am planning to consider a column with pyrite and infilling solution with equilibrium phase oxygen.Kinetics and rates are specified for all the cells.Is this way feasible?