Reactive Transport Modelling

[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 think the reaction is just very slow with that set of parameters. Here is a batch calculation with parm(1) increased by several orders of magnitude. You will have to consider what parameters should be changed to suit your environment.

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 you please give your valuable suggestions regarding the selection of time step for the transport an how to incorporate the kinetics with transport. I am considering a column of length 500m with background  solution. An infilling solution with equilibrium phase O2 is injected to it.Rate and Kinetics are defined for column as suggested above. Is that the correct formulation?

[dlparkhurst]

I have not looked carefully, but I think your formulation for transport is OK. I gave a batch kinetics example because it is faster and easier to see if the chemistry is working the way that you think it is. As always, start with a coarse discretization and see if the results are making sense before you refine the grid.