Oxidation and downward leaching in soils

[Tom]

Hi all,

I am looking to simulate the formation of a sulfate enriched layer in soil profiles, due to downward leaching via  infiltration. I think I can begin to simulate the relevant processes using TRANSPORT, considering the geology as a column. The following is my starting conceptual model:

- Pyritic soil, 5 m deep, assumed saturated (but in reality upper part is partially saturated)
- Rain (equilibrated with atmospheric CO2 and O2 - reasonable?) infiltrates profile
- Downward moving solution reacts with mineral phases
- Oxidation front moves down the profile, pyrite oxidises to sulfate
- Presence of Ca (calcite) allows gypsum formation


Ignoring the water table and assuming a simple model as a starting point I have the following input:

--------------------------

SOLUTION 1-50 #Initial solution in column (assumed pure water)
    temp      12
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    1 # kg

EQUILIBRIUM_PHASES 1-50
Pyrite 0 0.05
Calcite 0 10
Quartz 0 10
Kaolinite 0 10

SAVE SOLUTION 1-50

END

SOLUTION 0 #Displacing solution (generic rainwater)
units ppm
pH 5.47
Na 2.05
K 0.35
Ca 1.42
Mg 0.39
Cl 3.47
S(6) 2.19
N(5) 0.27
N(-3) 0.41

EQUILIBRIUM_PHASES 0
CO2(g) -3.5
O2(g) -0.699

SAVE SOLUTION 0

END

TRANSPORT
-cells 50
-length 0.1 # 50 cells @ 0.1m = 5 m total
-dispersivity 0.002 #m
-flow_direction forward
-boundary_conditions flux flux #at column ends
-diffusion_coefficient 1e-9 #m2/s
-punch_cells 1-50 #cells 1-50 in selected output
-time_step 144288 s #based on silty clay k of 0.25 cm/hr or 7e-7 m/s
-shifts 50

SELECTED_OUTPUT
-file sulfatep.xls
-pH true
-distance true
-totals S(6)
-SI gypsum

END

----------------

With this script I can plot the breakthrough curves of sulfate movement down profile (see attachment). At this stage I am wondering if it possible to simulate leaching - i.e. the upper zone becomes depleted in sulfate rather than remaining constant?

I think I am correct in my understanding that PHREEQC considers flow in saturated media and therefore cannot consider the unsaturated zone. Since leaching will occur in the unsaturated zone to the water table, would it be appropriate to model the cells below 2 m as stagnant?

Any other comments on the validity of my input file would be greatly appreciated,

Thanks,

Tom

[dlparkhurst]

The amount of sulfate generated is going to be related to the amount of oxygen in the inflowing water. Pyrite will react to equilibrium, consuming all of the oxygen, possibly modified by sequestering some sulfate in gypsum. Given this equilibrium, the concentration of sulfate will remain constant until the pyrite is removed, normally a long process (relative to the days you are plotting).

Although the flow system is saturated, you can do some things to represent the unsaturated zone. You could fix the partial pressure of O2 in the top cell, but you would not want to keep pyrite equilibrium. You would have to change it to a kinetic reaction (equilibrium between pyrite and O2(g) is not possible). That would allow a greater sulfate concentration to be generated, which would then propogate through the system. I still do not see a mechanism that would decrease sulfate in the UZ until you have removed (or decreased the reactivity) of pyrite.

[Tom]

Thanks David,

I will try your suggestions, particularly the use of kinetic pyrite oxidation and post an update.

Thanks,

Tom

[Tom]

Hi David,

Many thanks for your suggestions, I have proceeded with kinetic pyrite oxidation. As a start, I fixed the partial pressure of O2 in the top cell as that of air, 0.2 using GAS_PHASE.

Using a total time of around 45 years (10000 shifts), the output file was some 7 GB after 30 mins or so, but..

This appears to start depleting the near surface zone of sulfate after around 35-40 years, however sulfate concentrations are still very low (2-3 orders of mag below what I have found in the field).

In a separate file, I increased the ppO2 from 0.2 to 10 and this does seem to deplete the profile after a few years to the background of 20 mgSO4/kgw. SO4 increased to around 185 mgSO4/kgw after about half a year. Perhaps an exaggerated ppO2 is needed to form the profile I am observing in the field, at least over the time span I have so far considered?

Before I examine further changes, would mind taking a look at the input file to check whether it is reasonable please?

It also struck me that the pH doesn't drop below around 9, and begins as high as a very unrealistic 13.

Thanks for any guidance you may have,

Tom

INPUT FILE (and attached)

TITLE Sulfate profile formation (KINETICS with ideal mineralogy, fixing pO2 in top cell)

######INFILTRATING SOLUTION######

SOLUTION 0 #Displacing solution (rainwater)
units ppm
pH 5.47
temp 10
Na 2.05
K 0.35
Ca 1.42
Mg 0.39
Cl 3.47
S(6) 2.19
N(5) 0.27
N(-3) 0.41

EQUILIBRIUM_PHASES 0 #equilibrates with the atmosphere
CO2(g) -3.5
O2(g) -0.699

SAVE SOLUTION 0

END

######INITIAL COLUMN CONDITIONS######

SOLUTION 1 #Initial solution in column (assumed pure water, fixed pO2)
    temp      10
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    1 # kg

EQUILIBRIUM_PHASES 1
Calcite 0 3.59
Quartz 0 17.92
Illite 0 2.34
Ca-Montmorillonite 0 2.45
Kaolinite 0 0.90

GAS_PHASE 1 #Fix the partial-pressure of O2 in the top cell to simulate UZ
-fixed_pressure
-pressure 1.001
-temperature 10
O2(g) 0.2

SAVE SOLUTION 1

SOLUTION 2-50 #Initial solution in column (assumed pure water)
    temp      10
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    1 # kg

EQUILIBRIUM_PHASES 2-50
Calcite 0 3.59
Quartz 0 17.92
Illite 0 2.34
Ca-Montmorillonite 0 2.45
Kaolinite 0 0.90

SAVE SOLUTION 2-50

END

######RATES AND KINETICS######

KINETICS 1-50
Pyrite
-formula FeS2 1.0
-m0 1.06
-m 1.06
-parms 1.98 0.67 0.5 -0.11

RATES #Rate equation from Williamson and Rimstidt (1994), adapted to units of mol/dm3/s
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


######TRANSPORT######

TRANSPORT
-cells 50
-length 0.1 # 50 cells @ 0.1m = 5 m total
-dispersivity 0.005 #m
-flow_direction forward
-boundary_conditions flux flux #at column ends
-diffusion_coefficient 1e-9 #m2/s
-punch_frequency 100 #punch out every 100 shifts
-time_step 144288 s #based on silty clay k of 0.25 cm/hr or 7e-7 m/s
-shifts 10000

######OUTPUT FILE CONTENT######

SELECTED_OUTPUT
-file sulfatep.xls
-pH true
-pe false
-step false
-state false
-simulation true
-solution true
-distance true
-kinetic_reactants Pyrite
-saturation_indices Gypsum
-totals Ca

USER_PUNCH
-headings SO4_mg/kgw
-start
10 PUNCH TOT("S(6)") * GFW("SO4") * 1000
-end

END

[Tom]

Following the last update I have taken a closer look through some batch calculations and can see that the equilibrium dissolution of the clay minerals is causing the high pH.

Equilibrium dissolution is inappropriate in these time-scales and for the system in question. I have modified the file to consider kinetic dissolution of the mineral phases and a much increase total time. I may now also consider oxidation of pyrite by Fe3+ in the presence of dissolved oxygen, should the pH become low enough for this to be reasonable.

Will post an update later with simulations.

[dlparkhurst]

You should use O2(g) in EQUILIBRIUM_PHASES 1 to fix the partial pressure of oxygen (target SI of -0.7 corresponds to atmospheric O2). GAS_PHASE has only a limited amount of O2 that can react, whereas, by using EQUILIBRIUM_PHASES you can have a virtually unlimited amount (10 moles by default). The unlimited reservoir makes sense for the UZ connected to the atmosphere.

You probably did not react very much pyrite with GAS_PHASE, but you will react more with more O2 in EQUILIBRIUM_PHASES. With more oxidation of pyrite, you will want to add a ferric oxyhydroxide phase, like goethite or Fe(OH)3(a). By adding Fe(OH)3(a) you should generate a low pH.

You can decrease the output from TRANSPORT by using -punch_frequency and -print_frequency.

[Tom]

Thanks again David, will give this a go.

Tom

[Tom]

Thanks David, that has resolved the issue and enabled me to simulate the development of a sulfate peak at the base of the "unsaturated" zone (see images).

The simulation clearly has a number of assumptions and workarounds but has achieved what I was after and compares well to my study sites.

Thanks for the help!

[jchen]

Hello,

In your model, where can you express the diameter of the column?

[dlparkhurst]

There is no parameter for the diameter of the column. You need to reference reactants to the volume of solution. Normally, there is 1 kgw (~1 L), although other volumes are possible. Assuming 1 L of water, the other reactants (equilibrium phases, exchangers, etc) are defined in moles, which is effectively mol/L of water.

[jchen]

This is helpful! Thank you very much!

[rfembilejr]

Hi,

I was just wondering if Tom's plots, those with depth on the y-axis and concentration or pH on the x-axis, were made using Phreeqc.
If so, how can I make the same plots? I know with Phreeqc USER_GRAPH we plot ie. concentration vs x, where concentration usually lies on the y-axis and distance (depth) or time on the x-axis.
Thanks.

[Tom]

Hi rfembilejr, sorry for the very late response, but no, these plots were made in Excel from the Phreeqc output data.