Convergence problem: Diffusion and Oxidation

[md.muniruzzaman]

Hello everybody,

I am trying to set up a model for oxidation in a simple 1-D 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 pore-diffusion 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.088142e-04 Calculated: 1.432279e-04 Residual: -3.441371e-05

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 1-20
Calcite  0 7.62     # moles/L_water
Fe(OH)3(a) 0 0
Gypsum     0 0
Pyrite 0 0.85

SOLUTION 1-20
temp 25
units mol/L
pH 7
-water 1

SAVE solution 1-20
#END

#KINETICS 1-20
#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)") / 1e-6)
 40 if mol("O2") < 1e-6 then goto 80
 50 rO2 = 10^-8.19 * mol("O2")^0.5 * fH^-0.11          # ...rate with oxygen
 60 rO2_Fe3 = 6.3e-4 * 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.9e-6 * 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 1-20
#END

SOLUTION 0
temp 25
units mol/L
Tracer 2.556e-04
pH 7
-water 1

EQUILIBRIUM_PHASES 0
O2(g) -0.7 1000

#save solution 0
END

TRANSPORT
-cells 20
-lengths 5e-4   # m
-flow_direction diffusion  #forward #
-boundary_conditions constant flux   # at cell 1 and 20
-diffusion_coefficient 0.245e-9
-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 SO4-2 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("SO4-2"), equi("Fe(OH)3(a)")*1e-4, equi("Pyrite")*1e-4, mol("Ca+2"), equi("Gypsum")
 -end
END
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

[dlparkhurst]

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(SO4-2) in your case], and the numerical method fails.

I added  the following definition to your input file:

SOLUTION_SPECIES
H2O + 0.01e- = H2O-0.01
       log_k -9
END

This causes a redox aqueous species to be present with a concentration always near 1e-9. Sometimes it helps (sometimes not). In your case, the program ran to completion.

[md.muniruzzaman]

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.

[md.muniruzzaman]

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.