PhreeqcUsers Discussion Forum
Reactive Transport => Reactive Transport Modelling => Topic started by: md.muniruzzaman on October 26, 2016, 11:33:41 AM

Hello everybody,
I am trying to set up a model for oxidation in a simple 1D laboratory sized column containing Pyrite. The system I have conceptualized right now is:
 1D column column (~1cm) containing pyrite and calcite (water saturated)
 Oxic water (equilibrated with atmospheric O2) infiltrates the column by diffusion only (constant concentration boundary at one end of the column)
 Oxidation front propagates through the column by porediffusion and oxidize pyrite
 Formation of Fe(OH)3(a) and Gypsum (maybe?)
I think I should consider kinetic dissolution of pyrite, but as a very first step and to get only the sense if I am producing the right compounds, I have tried pyrite also as an equilibrium_phase.
The overall the simulation result, perhaps, makes sense (attached figure, results after ). However, I am having a major problem of S convergence after a long time! After around (~674 steps = 280.83 days), the simulation stops after the showing that S cannot converge (following error)! ERROR: S has not converged. Total: 1.088142e04 Calculated: 1.432279e04 Residual: 3.441371e05
I have tried with a smaller time_step (even smaller up to 300s, which fulfills Neumann stability criteria) in the TRANSPORT but the problem persists. Just for the curiosity, I tried only the advective transport instead of diffusion and that seems to be working fine (without any error)!
Does anybody have any suggestion on how to solve this? Any comments on the validity of my input file would be greatly appreciated.
Thanks,
Munir

Input File (also attached):
PRINT
reset false
SOLUTION_MASTER_SPECIES
Tracer Tracer 0.0 Tracer 100
SOLUTION_SPECIES
Tracer = Tracer
log_k 0.0;
EQUILIBRIUM_PHASES 120
Calcite 0 7.62 # moles/L_water
Fe(OH)3(a) 0 0
Gypsum 0 0
Pyrite 0 0.85
SOLUTION 120
temp 25
units mol/L
pH 7
water 1
SAVE solution 120
#END
#KINETICS 120
#Pyrite
#m0 0.85 # moles/L_water
#formula FeS2 1.0
#step 9e4 in 100 steps
#INCREMENTAL_REACTIONS true
RATES
Pyrite # rates from data compiled by Williamson and Rimstidt 1994, GCA 58, 5443
start
#1 A = 120 * m0
1 A = 1.43 * m0 # initial surface area in m2/mol_Pyrite indicates 50 um size crystals
10 if SI("Pyrite")>0 then goto 100 # step out when supersaturated...
20 fH = mol("H+")
30 fFe2 = (1 + tot("Fe(2)") / 1e6)
40 if mol("O2") < 1e6 then goto 80
50 rO2 = 10^8.19 * mol("O2")^0.5 * fH^0.11 # ...rate with oxygen
60 rO2_Fe3 = 6.3e4 * tot("Fe(3)")^0.92 * fFe2^0.43 # ...rate with oxygen and Fe3+
70 goto 90
80 rem # rate with Fe3+ without oxygen, and for pH < 3...
81 rFe3 = 1.9e6 * tot("Fe(3)")^0.28 * fFe2^0.52 * fH^0.3
90 rate = A * (m/m0)^0.67 * (rO2 + rO2_Fe3 + rFe3) * (1  SR("Pyrite")) # ...sum terms, zero at equilibrium
100 save rate * time
end
SAVE SOLUTION 120
#END
SOLUTION 0
temp 25
units mol/L
Tracer 2.556e04
pH 7
water 1
EQUILIBRIUM_PHASES 0
O2(g) 0.7 1000
#save solution 0
END
TRANSPORT
cells 20
lengths 5e4 # m
flow_direction diffusion #forward #
boundary_conditions constant flux # at cell 1 and 20
diffusion_coefficient 0.245e9
shifts 1000 # number of diffusive timesteps
punch_frequency 300 # graphing is done after each shift
time_step 36000 # seconds
USER_GRAPH
headings dist pH O2 Tracer SO42 Fe(OH)3 Pyrite, Ca+2 Gypsum
axis_titles "dist (mm)", "pH", "Conc"
initial_solutions false
connect_simulations false
plot_concentration_vs x
start
10 graph_x dist*1e3
20 graph_y LA("H+")
30 graph_sy mol("O2"), tot("Tracer"), mol("SO42"), equi("Fe(OH)3(a)")*1e4, equi("Pyrite")*1e4, mol("Ca+2"), equi("Gypsum")
end
END


That error message occurs occasionally. I think the redox conditions are such that two equations are the same within the machine tolerance; another way to say it is that the species that should cause two equations to differ are very small in concentration. In this case, one unknown is not solved [act(SO42) in your case], and the numerical method fails.
I added the following definition to your input file:
SOLUTION_SPECIES
H2O + 0.01e = H2O0.01
log_k 9
END
This causes a redox aqueous species to be present with a concentration always near 1e9. Sometimes it helps (sometimes not). In your case, the program ran to completion.

Thanks a lot Dr. Parkhurst. This seems to be a very good trick and perhaps solves my issue. I just tried running the simulation for 10 years without having that problem.
I will also post an update after including kinetic oxidation of pyrite and with a refined grid.

Dr. Parkhurst,
I think it has solved my issue. I have just finished running the simulation for 10 years (took almost 80 hours!) including kinetic oxidation of Pyrite.
Thanks again for the help.