Constraining the supply of oxygen during water infiltration through tailings

[bschroth]

Hello - I am working on a site in which a tailings impoundment has been receiving irrigation water.  While no groundwater chemistry data are available, ABA data have been collected in several boreholes through the tailings, providing a profile of paste pH, acid-generating potential (AGP), and acid neutralization potential (ANP) with depth.  The data indicate low paste pH (3-5.5) in the top 15 feet and essentially neutral pH below that depth.  All samples show AGP>>ANP, but ANP remains detectable below the upper 15 feet.  Water table is at 50 feet.  The working conceptual model is that oxygen is the limiting factor on the movement of the acid front downward, with pyrite oxidation using up dissolved oxygen in the infiltrating irrigation water and oxygen is diffusion-limited at greater depths in the vadose zone.  I am trying to build a "simple" PHREEQC advection model to illustrate this conceptual process, but am having issues getting the model to converge in the >15-foot depths.  I am simulating kinetics for pyrite and calcite as provided in the PHREEQC documentation along with kinetic expressions for O2 and CO2 diffusion from LeGuern et al. (2003).  The advection block specifies 10 cells (5 feet each) with a time step of 2 years (estimated infiltration velocity from separate model). I have tried using separate kinetic/rate blocks for the upper three cells and the lower seven cells, with the latter containing lower O2 diffusion rate, but the model will not converge.  The closest I have come to illustrating the system is to specify a low SI for O2(g) (-55) for the lower cells.  The model input is shown below (database = llnl.dat).

I would prefer to apply a kinetic model of oxygen diffusion to illustrate the eventual breakthrough of an acid front in the lower cells over time - is there any way I can use the kinetics O2 expression to achieve this without non-convergence issues, or is this oxygen diffusion vs pyrite dissolution kinetics presenting too great a conflict for PHREEQC to provide a reasonable solution?  Thanks in advance for your input

PHASES
# from phreeqc.dat
Al(OH)3(a)
   Al(OH)3 + 3 H+ = Al+3 + 3 H2O
   log_k      10.800
   delta_h -26.500 kcal

Fe(OH)3(a)_pqe
   Fe(OH)3 + 3 H+ = Fe+3 + 3 H2O
   log_k      4.891

# from wateq4f.dat
Basaluminite
        Al4(OH)10SO4 + 10H+ = 4Al+3 + SO4-2 + 10H2O
        log_k           22.7

Witherite_wateq
        BaCO3 = Ba+2 + CO3-2
        log_k           -8.562
        delta_h 0.703 kcal
        -analytical     607.642         0.121098        -20011.25       -236.4948       0.0

Jarosite(ss)        133
        (K0.77Na0.03H0.2)Fe3(SO4)2(OH)6 + 5.8H+ = 0.77K+ + 0.03Na+ + 3Fe+3 + 2SO4-2 + 6H2O
        log_k           -9.83   # WATEQ4F, Alpers and others, 1989

Jarosite-Na(wateq4f)
        NaFe3(SO4)2(OH)6 + 6H+ = Na+ + 3Fe+3 + 2SO4-2 + 6H2O
        log_k           -5.280
        delta_h -36.180 kcal

Jarosite-K(wateq4f)
        KFe3(SO4)2(OH)6 + 6H+ = K+ + 3Fe+3 + 2SO4-2 + 6H2O
        log_k           -9.21
        delta_h -31.28 kcal

Jarosite-H(wateq4f)
        (H3O)Fe3(SO4)2(OH)6 + 5H+ = 3Fe+3 + 2SO4-2 + 7H2O
        log_k           -5.390
        delta_h -55.150 kcal


SELECTED_OUTPUT
        -file Infil_06.sel
        -reset false
USER_PUNCH   
-heading   pH pe  Na K Ca Mg Fe Mn Cl SO4 Alkalinity SiO2 mg_Gypsum
01  PUNCH DESCRIPTION   
01  PUNCH LA("H+")*-1   
02  PUNCH LA("e-")*-1   
05  PUNCH TOT("Na")*RHO*1e3*22.9898   
10  PUNCH TOT("K")*RHO*1e3*39.0983   
15  PUNCH TOT("Ca")*RHO*1e3*40.078   
20  PUNCH TOT("Mg")*RHO*1e3*24.305   
25  PUNCH TOT("Fe")*RHO*1e3*55.847   
30  PUNCH TOT("Mn")*RHO*1e3*54.938   
35  PUNCH TOT("Cl")*RHO*1e3*35.453   
40  PUNCH TOT("S(+6)")*RHO*1e3*96.066   
45  PUNCH ALK*RHO*1e3*61.0173      
70  PUNCH TOT("Si")*RHO*1e3*60.0843
75  PUNCH EQUI_DELTA("Gypsum")*1e3*172


SOLUTION   0   Assumed Irrigation Water
temp   21   
pH   8.5   ##not available-use average of historical samples
pe   12   
redox   pe   
units   mg/L   
#density   1.00011815   
Na   3.9   
K   1.3   
Ca   49    
Mg   16   
Fe   0.05001 #Fe was ND so this is placeholder
Mn   0.003   
Cl   1.9   
S(+6)   9.8   as SO4
Alkalinity   180   as Ca0.5(CO3)0.5
Si   14   as SiO2
#As   0.002   
#O(0) 8.68

PURE 11
Calcite 0.0 0.0
#CO2(g) -3.2 10.0
#O2(g) -0.7 10.0

SAVE solution 0
END

SOLUTION 1-10
pH 7
temp 20


ADVECTION
  -cells 10
  -shifts 20
  -time_step 2 year
  -initial_time 0
  -punch_frequency 1
  -punch_cells 10

INCREMENTAL_REACTIONS true
RATES 1-3
  Pyrite
      -start
   rem  PARM(1) = log10(A/V, 1/dm)
   rem  PARM(2) = exp for (M/M0)
   rem  PARM(3) = exp for O2
   rem  PARM(4) = exp for H+
 10  if (M <= 0) then goto 200
 20  if (SI("Pyrite") >= 0) then goto 200
 30 lograte = -1.0190E+01 + PARM(1) + PARM(2)*LOG10(M/M0)
 40 lograte = lograte + PARM(3)*LM("O2") + PARM(4)*LM("H+")
 50 moles = (10^lograte) * TIME
 60 if (moles > M) then moles = M
 200 SAVE moles
      -end
   Calcite
         -start
          rem M = current number of moles of calcite
          rem M0 = number of moles of calcite initially present
          rem PARM(1) = A/V, cm^2/L
          rem PARM(2) = exponent for M/M0
  10 si_cc = SI("Calcite")
  20 if (M <= 0 and si_cc < 0) then goto 200
  30 k1 = 10^( 0.198 - 444.0 / TK )
  40 k2 = 10^( 2.84 - 2177.0 / TK)
  50 if TC <= 25 then k3 = 10^(-5.86 - 317.0 / TK )
  60 if TC > 25 then k3 = 10^(-1.1 - 1737.0 / TK )
  70 t = 1
  80 if M0 > 0 then t = M/M0
  90 if t = 0 then t = 1
  100 area = PARM(1) * (t)^PARM(2)
  110 rf = k1*ACT("H+")+k2*ACT("CO2")+k3*ACT("H2O")
  120 rem 1e-3 converts mmol to mol
  130 rate = area * 1e-3 * rf * (1 - 10^(2/3*si_cc))
  140 moles = rate * TIME
  200 SAVE moles
        -end
    Oxygen
        -start
  10 si_o2 = SI("O2(g)")
  20 rate =  1.0E-08 * (10^-0.7 - 10^si_o2)
  30 moles = rate * TIME
  40 SAVE moles
        -end
    CO2
        -start
  10 si_co2 = SI("CO2(g)")
  20 rate =  2.0E-09 * (10^-3.5 - 10^si_co2)
  30 moles = rate * TIME
  40 SAVE moles
        -end
 
KINETICS 1-3 Pyrite and calcite formulas from PHREEQC manual example to include dolomite
  Pyrite
    -formula FeS2 1.0
    -m   11.72
    -m0  11.72
    -parms  3.0 0.67 0.5 -0.11
    -tol 1e-9
   Calcite
     -formula CaCO3  0.5  MgCO3  0.5
     -m    0.88
     -m0   0.88
     -parms  0.5 0.3
     -tol 1.e-8
   Oxygen
     -formula O2
     -m  1.0e5
     -m0 1.0e5
     -tol 1.e-8
   CO2
     -formula CO2
     -m  1.0e5
     -m0 1.0e5
     -tol 1.e-8
#-steps  2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 year
-bad_step_max  1000
#-cvode  false
#-runge_kutta  6

EQUILIBRIUM_PHASES 1-3
   Fe(OH)3(a)_pqe  0.0  0.0
   Rhodochrosite  0.0 0.0
   Siderite   0.0 0.0
   Gypsum   0.0 0.0
   SiO2(am) -0.5  0.0
   Basaluminite  0.0  0.0
   Illite  0.0  10.0
 
EQUILIBRIUM_PHASES 4
   Fe(OH)3(a)_pqe  0.0  0.0
   Rhodochrosite  0.0 0.0
   Siderite   0.0 0.0
   Gypsum   0.0 0.0
   SiO2(am) -0.5  0.0
   Basaluminite  0.0  0.0
   Illite  0.0  10.0
   Calcite 0.0 1.38
   Pyrite 0.0 14.28
   O2(g) -55
 
EQUILIBRIUM_PHASES 5
   Fe(OH)3(a)_pqe  0.0  0.0
   Rhodochrosite  0.0 0.0
   Siderite   0.0 0.0
   Gypsum   0.0 0.0
   SiO2(am) -0.5  0.0
   Basaluminite  0.0  0.0
   Illite  0.0  10.0
   Calcite 0.0 1.82
   Pyrite 0.0 14.6
   O2(g) -55
 
EQUILIBRIUM_PHASES 6
   Fe(OH)3(a)_pqe  0.0  0.0
   Rhodochrosite  0.0 0.0
   Siderite   0.0 0.0
   Gypsum   0.0 0.0
   SiO2(am) -0.5  0.0
   Basaluminite  0.0  0.0
   Illite  0.0  10.0
   Calcite 0.0 1.79
   Pyrite 0.0 13.07
   O2(g) -55
 
EQUILIBRIUM_PHASES 7
   Fe(OH)3(a)_pqe  0.0  0.0
   Rhodochrosite  0.0 0.0
   Siderite   0.0 0.0
   Gypsum   0.0 0.0
   SiO2(am) -0.5  0.0
   Basaluminite  0.0  0.0
   Illite  0.0  10.0
   Calcite 0.0 1.57
   Pyrite 0.0 11.72
   O2(g) -55
 
EQUILIBRIUM_PHASES 8
   Fe(OH)3(a)_pqe  0.0  0.0
   Rhodochrosite  0.0 0.0
   Siderite   0.0 0.0
   Gypsum   0.0 0.0
   SiO2(am) -0.5  0.0
   Basaluminite  0.0  0.0
   Illite  0.0  10.0
   Calcite 0.0 1.76
   Pyrite 0.0 10.38
   O2(g) -55
 
EQUILIBRIUM_PHASES 9
   Fe(OH)3(a)_pqe  0.0  0.0
   Rhodochrosite  0.0 0.0
   Siderite   0.0 0.0
   Gypsum   0.0 0.0
   SiO2(am) -0.5  0.0
   Basaluminite  0.0  0.0
   Illite  0.0  10.0
   Calcite 0.0 0.26
   Pyrite 0.0 9.78
   O2(g) -55
 
EQUILIBRIUM_PHASES 10
   Fe(OH)3(a)_pqe  0.0  0.0
   Rhodochrosite  0.0 0.0
   Siderite   0.0 0.0
   Gypsum   0.0 0.0
   SiO2(am) -0.5  0.0
   Basaluminite  0.0  0.0
   Illite  0.0  10.0
   Calcite 0.0 0.23
   Pyrite 0.0 8.63
   O2(g) -55
 
SURFACE 1-10
   Hfo_sOH   Fe(OH)3(a)_pqe  equilibrium_phase   0.001  1e5
   Hfo_wOH   Fe(OH)3(a)_pqe  equilibrium_phase   0.1
 
END

[dlparkhurst]

Use -cvode in KINETICS. I made this change and the file ran to completion, but I did not look at the output. Check to see if the results are reasonable.