Strange pressure output from EOS7C

I have been struggling for several months now to understand how pressure is accounted for in some TOUGH2-EOS7C simulations that I’ve been running. The fact that the code outputs a single “Pressure” variable has produced some very problematic results that are very difficult to interpret, and the manual does not provide enough detail to help, so I’m hoping that someone may be able to clarify for me why it is that I’m seeing what I’m seeing. I will explain my issue through an example simulation.

The 1-D system consists of a very low-permeability “shale” unit bounded on the top and bottom by high permeability “sandstone” layers. The constitutive relationships applied to the shale cause extremely high capillary pressures, even at small gas saturations. The system is initialized to a hydrostatic, fully water-saturated condition and then gas phase methane is injected in the middle of the domain for 10,000 years (10 ka). After that, injection is stopped and the system is allowed to continue evolving up to a final time of 1 Million years (1 Ma).

During the injection, gas phase forms and increases in saturation, and the pressure increases everywhere. Then, after injection stops, the pressure goes back to hydrostatic, but gas phase methane continues to evolve. At first glance, this makes sense, but when one tries to understand the evolution of the two phases separately, severe confusion arises. In the TOUGH output, “pressure” is that of the liquid phase when/where only liquid phase is present, but when/where gas phase is present, “pressure” refers to gas phase pressure. Therefore, smoothly varying “pressure” profiles can imply single-phase pressure profiles with extremely sharp gradients that persist for a very long time in the tight rocks that I’m trying to model. Please see the attached plots for results from an example simulation.

In Figure01_P, you can see the initial increase in pressure during the injection, and then a decrease back to hydrostatic afterward. Figure02_Sw, however, shows that gas phase is present in part of the domain throughout the entire simulated timeframe. This means that, in the range where gas phase exists, the hydrostatic P given by the output is in fact gas pressure, as is shown in Figure03_Pg. Where gas phase exists, capillary pressure must also exist, which was calculated from the van Genuchten relationship that was input for the model (see Figure04_Pc). Using the definition of capillary pressure, Pc = Pg – Pw, the water pressure can be calculated in the region where gas phase exists, but it is hydrostatic elsewhere (where there is only water). The results of these calculations are shown in Figure05_Pw.

The water pressure plot is disturbing in the fact that huge gradients exist at the edges of the gas-occupied region, and that these gradients persist for such a long time. Furthermore, the gradients do not translate into the flow patterns that one would expect to see. That is, there should be a very large flow rate where the gradient is really steep, but this is not the case. Likewise, the point where the minimum pressure exists should correspond with a reversal in the water flow direction, but this is not observed either. I realize that the presence of gas phase reduces the relative permeability of the water phase, but by less than an order of magnitude, and this should affect the magnitude but not the direction of water flow, so I don’t think relative permeability alone provides a sufficient explanation for the behavior.

The crux of this problem seems to be rooted in the fact that TOUGH forces the gas phase pressure to attain a hydrostatic pressure profile through time, which doesn’t physically make sense to me. Intuitively, the gas phase pressure should instead achieve a steeper (“gas-static”) profile, and should have a pressure higher than hydrostatic. Can someone please explain to me why TOUGH does this? Thank you very much.