Gas handling in reactive transport calculations with PHREEQC

[mpek]

Dear all. I was wondering if anybody could help with the following problem.

I am carrying out simple reactive transport calculations with PHREEQC over a single material domain. In the domain steel is anaerobically corroded (kinetic reaction with a constant rate) and H2 is generated as a result. Boundaries of the domain are kept at constant concentration (redox at boundaries is higher than in the domain). I am intersted in the interplay between H2 production by steel corrosion and O2/H2 diffusion, and the impact on redox potential in the system.

The problem is implemented in the following way:
I define H2(g) as an equilibrium phase with zero initial concentration and a maximum partial pressure (saturation index) of 50 bar (corresponding to hydrostatic pressure at 500 metre depth):

EQUILIBRIUM_PHASES 2-9 
    H2(g)   1.7   0.0

What I want to see is built-up in H2(aq) up to the value defined by the maximum partial pressure of H2(g). And it does what is expected: H2(aq) increases with steel corrosion (linear with time) to a maximum concentration (ca. 0.04 mol/kgw) corresponding to the limiting value of partial pressure of H2(g) = 50 bar.

However, once the maximum concentration is reached H2(g) continues to accumulate in the system as an equilibrium phase. This "precipitated" hydrogen constitutes a large hydrogen sink and provides redox buffering capacity that should not be there. In other words, I would like to see only the redox buffering capacity associated with the H2(aq) and eliminate that associated with the H2(g) equilibrium phase.

I have tried various approaches to achieve this, but with limited success so far:
1) Using GAS_PHASES with fixed gas pressure of 50 bar and initial amount of hydrogen = 0. This does the same as above (dissolved hydrogen and partial pressure are OK), but now hydrogen generated from
anaerobic steel corrosion accumulates as a gas phase (and not as an equilibrium phase) leading to the same redox buffering effect...

2) If no partial pressure limit on H2(g) is implemented, H2(aq) increases linearly with corrosion to unacceptably high values (above the 50 bar H2(g) partial pressure limit).

3) I tested a variety of ways to convert the accumulating H2(g) (either by means of a kinetic or equilibrium reaction) into something else that would not participate in redox calculation. This however gives me lots of convergence issues and never seems to work completely.

4) In a simple batch reaction I managed to use the EQUILIBRIUM_PHASES_MODIFY keyword to eliminate accumulating H2(g). In this approach I would like to automatically correct the #of moles of H2(g) to zero (or some small number) after each chemistry calculation. But it does not work when I use it in conjunction with the TRANSPORT keyword.


I wonder if somebody had an idea how else this can be handled.



mpek

[dlparkhurst]

Have you tried using "precipitate_only" in the definition of H2(g) in EQUILIBRIUM_PHASES? Seems like that might work.

Alternatively, I think a KINETIC reaction might work.

[mpek]

Thank you very much. I think that using "precipitate_only" should fix the issue.

Best regards,

M