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.
27 replies
-
Hi Michael,
Only briefly:
(1) While I don't know all the details of the input of your simulations (b.c., properties) nor do I see the detailed output (e.g., fluxes), I don't see anything "problematic".
(2) TOUGH2 prints out the gas pressure and the capillary pressure. That's sufficient. Also printing the liquid pressure would be redundant, as Pl = Pg + Pc (Pc is defined negative in TOUGH2, i.e., Pc = Pl - Pg, not Pc = Pg - Pl, as you write - your plots are correct, though). So again, I don't see anything "problematic" or "difficult to interpret".
(3) Yes, strong pressure gradients exist near phase boundaries; they may induce strong fluxes (don't forget the gravity term when you calculate fluxes!). While I don't see your fluxes, I am very (!) confident that they are calculated correctly. Fortunately, this is very easy to check. Just take two elements in the "problematic" region, extract pressures, relative permeabilities (consider upstream weighting!), densities, viscosities, gravity, nodal distances and cross-sectional area, and plug these values into Darcy's law. You can do this on a piece of paper and compare your number with the flux printed in the TOUGH2 output file for the corresponding connection - the two values will be identical (they may be smaller than what you expected, but they are nevertheless correct).
(4) It is not the case the "TOUGH forces the gas pressure to attain a hydrostatic pressure" - it solves the set of governing equations for the problem at hand (you can easily simulate cases where it gives the "gas-static" gradient you expected and considered more reasonable - but that would be solving a different problem!). Do not expect anything "hydrostatic" if your system has any vertical flow and is not even at steady state! Sorry, I cannot go into more details (due to lack of information and time), but I am again very confident that the TOUGH solution is correct for the problem you specified.
In summary: I don't think this is related to TOUGH at all (check with another code!), but most likely just a difference between expectations and results.
Good luck!
Stefan
-
Stefan,
Thank you for your response. I’ve tried some example calculations as you suggested, but I’m still having trouble interpreting them due to confusion about how the code handles the capillary pressure function. This could be alleviated by checking the capillary pressure output, but I don't see that anywhere. Am I missing something? I searched the output file for “cap,” “pc,” and “cp,” but didn’t see any numbers that could be interpreted as capillary pressure output. This is why I’ve been calculating capillary pressures on my own from the saturation output and the van Genuchten (vG) parameters that I gave as inputs to the model. This brings me to my next question:
How exactly is the parameter CP(3) defined for the vG capillary pressure model in a TOUGH2 input file? In the version of the vG formula that’s given in the TOUGH2 manual (Pcap = -Po([S*]^(-1/lambda)-1)^(1-lambda)), Po should be in units of pressure for the equation to give Pcap in terms of pressure. However, multiple references about TOUGH2 models refer to a “1/alpha” term instead of Po. For instance, the handouts from the TOUGH2 beginner’s class that I took last year show the equation above, except with 1/alpha in place of Po. Rearranging the equation in the CP(3) definition that relates Po and alpha (1/Po = alpha/(rhow*g)), the term 1/alpha = Po/(rhow*g) would be in units of length, not pressure, so the vG formula given in the class would give capillary head instead of capillary pressure. So my question is, should CP(3) be input as entry pressure in units of pressure (Po), entry pressure in units of length (1/alpha), or the reciprocal one of those (1/Po, in units of 1/pressure, or alpha, in units of 1/length)?
For me to set up my simulations correctly, it is imperative that I understand how the vG capillary pressure equation is formulated in the code. Hopefully clarifying this will also help me to correctly post-process and understand my simulation results.
Thank you again in advance for your time,
Mike
-
Mike,
Set KDATA=2 or higher, and you will find PCAP in the output file.
Go by the manual! As described therein, CP(3)=1/P0 [1/Pa].
Regards,
Stefan
-
Stefan,
Thank you for clarifying the CP(3) definition.
The attached files "run.txt" and "out.txt" are the input and output files from the simulation that produced the graphs in my original post. As you can see, I already had KDATA set to 3, but there are no numbers for PCAP in the output file. Is there some other flag that I need to turn on or off or something?
Based on your clarification of CP(3), I know now that the parameters that I input into my model, as well as the postprocessing calculations I performed to make the original graphs, were incorrect (I had been interpreting CP(3) as Po, not 1/Po, based on the way that I had seen parameters reported for the field site that I am trying to model). The hydrostatic gas profile now makes sense, because the way I input the numbers made the capillary pressures very small, when I was in fact trying to make them very large.
I also did some example calculations and found that the outputted fluxes do in fact make sense, when I account for the pressure correctly.
Thanks again,
Mike
-
Mike,
Glad to hear that the results start to make sense.
Please contact the developer of EOS7c regarding the printout of PCAP (in case you have a source code license, it's very easy to add the printout you want). In the meantime, see the attached output file, which shows the PCAP output.
Best,
Stefan
-
Stefan,
I will contact Curt Oldenburg about EOS7C.
Thanks for the example output. I assume that you ran my simulation on your own to produce this file? If so, I am confused about the capillary pressure output it gives.
If what you said about CP(3) is correct, then the capillary pressure-saturation functions that I input into the simulation should look like the ones in the attached file (on the order of 1E-5 Pa for the "shale" rock across a wide range of saturations), but the output you gave me shows PCAP values on the order of 1E+5 or even 1E+7. This discrepancy could be consistent with the code interpreting CP(3) as Po, instead of 1/Po, as I had originally thought.
Could you please verify once more that CP(3) is not in fact supposed to be Po in Pa? If it is, then the water flow rates that I've calculated from the outputted pressures are several orders of magnitude off from the outputted water flow rates.
Thanks again,
Mike
-
Mike,
I used iTOUGH2-EOS7c (I forgot whether I customized the output - which I often do...).
I took the liberty to correct your CP(3) parameters (taking their inverses) before I ran the simulation. The values you had did not make sense. I felt empowered to do so as you acknowledged the mistake in your previous post.
Hope this clarifies it.
Stefan
-
Stefan,
Yes, that clarifies it, and the pressure and flux results now make sense to me. I will look into using iTOUGH2-EOS7c and/or contact Curt about allowing PCAP output from the original EOS7c.
Thank you so much for all your help and for bearing with me through my confusion.
Mike
-
Hello,
I've acquired iTOUGH2-EOS7c but there is no executable, and I am quite lost trying to figure out how to compile the code myself. I installed the files to the suggested directory (C:\iTOUGH2), copied the EOS source code into the "SourceCode" subdirectory, then opened a command prompt and navigated there, but when I type "make," "nmake," or "it2make," it always gives an error message saying that it "is not recongized as an internal or external command, operable program, or batch file." I have the Intel Parallel Studio XE 2017 compiler software installed on my machine, but this is my first time trying to compile TOUGH2, so any help would be greatly appreciated.
Thank you very much in advance,
Mike
-
Michael,
On a PC, you cannot use the unix-style make command (unless you installed, for example, cygwin, an unix-sty;e environment for PCs). Instead, you have to follow the instructions in file read_me_CVF.txt, which was distributed along with the code.
Let me know if you need additional help (but be patient, as I don't have ready access to a PC).
Stefan
-
Mike,
If you decide to go the Cygwin route, I would highly suggest the following blog that explains how to install the components you'll need to compile the codes using the Unix command stucture:
http://fundamentalthinking.blogspot.com/2015/02/installing-fortran-on-windows-with.html
If you'd like, you can skip the last step of installing Eclipse (unless you want to start getting into compiling your own code), but pay special attention the second comment given on the post, which highlights an installation package you'll need that the original author forgot to mention.
-
I am unable to follow the directions in "read_me_cvf" because I do not have the Compaq Visual Fortran compiler. Likewise, using Cygwin seems unnecessarily complicated to me because it involves switching back and forth between Unix and Windows environments. Since I plan to run iTOUGH2 on a PC, and I have a compiler installed on the machine already (the Intel Fortran Compiler), I decided to try compiling the code by adapting project/solution files from a previous compilation of a forward-running version of TOUGH2 that I had done.
I was able to create an executable through the above method, but I still have some general questions about compiling iTOUGH2: (1) does the file called T2 need to be included, like when compiling forward TOUGH2? It was not there when I compiled the code. Also, (2) should I use all of the same settings as when compiling forward-running versions of TOUGH2? If not, (3) what options should be changed?
I went ahead and tried to run a simulation with the executable that I created, but more questions arose from this execution step: (a) what files need to be present in the working directory where I perform execution? Also, (b) can I directly use a forward TOUGH2 input file to execute iTOUGH2 in forward mode? If not, (c) what changes do I need to make to the input file?
Attached is an input file for a simulation that ran successfully in the forward version of EOS7C. I copied "iT2_7c.exe" that I had created above into "C:\iTOUGH2\Executable" and tried running the forward simulation with the same input file (this was the only file in my working directory at the time). As the screenshot shows, the code crashed the first time because it couldn't find a file called "status." Therefore, I created an empty file called "status" in the working directory and re-ran the simulation. This time, a FATAL ERROR occurred again, but no reason was given. The output file is also attached. The code also created a file called "input.txt.msg," but I do not know how to open this file.
Any help would be greatly appreciated.
-
Michael,
Just briefly:
(1) Follow the iTOUGH2 installation read me file I sent you (Digital => Compaq => Intel compilers are just different names for essentially the same compiler, so the instructions are valid and easy to follow if you know how to load a project and create an executable using your particular version of the compiler). Do not use the TOUGH2 installation instructions (e.g., file T2 is obviously not used in iTOUGH2)! The key is that you simply compile and link the list for Fortran source code files given in these instructions.
(2) In the instruction file you will also learn how to run iTOUGH2 in forward mode (no changes needed in the input file). Minor changes and additional features in the forward model are described in http://esd1.lbl.gov/FILES/research/projects/tough/documentation/TOUGH2-In-iTOUGH2_Enhancements.pdf .
(3) Did you look at the output files (not just the screen output you attached) of the run you made for error messages or other information?
Again, read and use the instructions provided (all I could do here is repeat them). If you still have issues, please provide detailed error messages,
Stefan
-
Stefan,
I apologize, I neglected to change the name of the output file so that it would attach to my previous message without zipping itself. I've attached it here, along with the "msg" file that I was able to open this time, after renaming it. I don't see any detailed error messages in either of these files, nor in the the domain geometry/initial condition definition files (GENER, INCON, LINEQ and MESH) that were created when I ran the executable. I've also renamed and attached these for your reference.
I tried again to follow the directions in "read_me_cvf," using my best judgement as to how the steps should be implemented in a different compiler framework that has different words and options, but was unable to successfully create an executable this time. There were two steps in your instructions that I couldn't follow exactly. MFC was not mentioned anywhere in my compiler options, so I was not able to implement the step of "Not Using MFC." Likewise, there was no option for "Maximum optimization," so I picked the one that said "Maximize speed plus Higher Level Optimizations," because this seemed like the closest match.
When I tried to compile after making these changes, I got error messages corresponding to each of the source code files, saying that it "can't open the file for write."
Mike
-
Mike,
output_msg.txt tells you quite precisely what the problem is. Increase these dimension parameters in maxsize.inc and recompile.
It seems you were able to create an executable before - what did you change?
MFC stands for Microsoft Foundation Class Library - as the instructions say: you don't need it. Optimization levels don't matter much, either. Check your file permissions or the path to the directory where you keep your source code. All of this is related to how you set up your project - not the source code itself.
Hope this helps - good luck!
Stefan
-
Stefan,
Thank you for clarifying. After increasing MAXEQ and MAXK to 5, the simulation ran, and capillary pressure values were indeed included in the output.
I apologize for asking so many questions.
Mike
-
Mike,
No worries - glad it works.
Stefan
PS: Generally, MAXEQ=MAXK+1; for EOS7C, MAXEQ=6, otherwise you will have to recompile should you try to do a non-isothermal simulation.
-
Stefan,
Great. Thanks for your advice. I re-compiled with MAXEQ = 6 and made sure that the simulation still ran with the new executable, which it did.
I am very happy to have a working version of iTOUGH2 now, because it really opens up a wide range of new modeling possibilities for me (e.g., incorporating hysteresis, doing inverse modeling, and using the Brooks & Corey formulation).
I have one last question for now if you don't mind: When I first acquired iTOUGH2, I noticed that there was another equation of state for modeling water and methane (eos13). The Licensing and Download website gives a link to the regular TOUGH2 manual next to eos13, but this eos is not mentioned in this manual, nor the iTOUGH2 manual. Could you please let me know where I can learn more about eos13, in particular how it varies from eos7c?
Thank you very much again for all your help,
Mike
-
Mike,
eos13.f has the same capabilities as eos3, with CH4 properties replacing those for air. It does not have two gas components (as eos7c and eos7ca do). This module is no longer supported and will disappear in iTOUGH2 V7.1, released any day (or week/month/year) now!
Stefan
-
Stefan,
That's good to know; I will keep using EOS7c then. I'm sorry, but I thought of one more question already: do you know of any plans to incorporate the functionality of EOS7c-ECBM (specifically, sorption of gases) into the iTOUGH2 framework?
Thanks again,
Mike
-
Mike,
If USGS provides appropriate funding, I'd be happy to incorporate EOS7c-ECBM into iTOUGH2... ;)
(As you may know, I'm no longer at LBNL and I have no general funds that I could use for iTOUGH2 development or support).
Note that you may always use iTOUGH2-PEST and link it to the stand-alone TOUGH2-EOS7c-ECBM executable; it's tedious, but the rest of the world that uses PEST has to deal with the PEST template and instruction files in exactly the same manner.
Stefan
-
Stefan,
Thank you for suggesting that I use iTOUGH2-PEST. Again, I think that this will set me up well to be able to address a wider range of problems, since parameter estimation is something that I have wanted to learn how to do for quite a while anyway. I think I'm close to getting it to work, but I am stuck once again.
Could you please let me know what I'm doing wrong, based on the following summary of the procedure that I tried?
- Created iT2_pest.exe by compiling the following files: eospest.f, it2input.f, it2main.f, it2pest.f, it2stubs.f, it2user.f, it2xxxx.f, mdeppc.f, meshm.f, mio.f90, t2cg22.f, t2f.f, and t2solv.f
- I used the same compiler project as when I sucessfully compiled iT2_7c.exe, but replaced eos7c.f and zevsreal.f with eospest.f, it2pest.f and mio.f90, because the compiler gave error messages when these three files weren't included.
- Successful compilation required me to change the names of the pest subroutines in it2stubs.f
- I got a number of warnings about the mio.f90 file, as shown in the screenshot below, but the executable was still created, so I moved on anyway.
- I copied the executable into the directory C:\iTOUGH2\Executable, because I learned through my experience with iT2_7c that this was necessary.
- Created an inverse input file called testi, to which I've added the .txt extension and attached it here for your reference.
- In this file, I only included the command for the forward execution of TOUGH2_EOS7C_ECBM.exe, using the same input file as we have been discussing throughout this entire thread (TOUGH2.inp, which I've renamed as TOUGH2.txt and attached here as well).
- Created a forward input file called testf, to which I've added the .txt extension and attached it here for your reference.
- This file only includes the characters PEST
- Created an empty file called status (not attached).
- Executed the code by navigating within command prompt to the directory where the aforementioned files were created and typing: itough2 testi testf pest
- This created the output file testi.out, which I've renamed as testi_out.txt and attached here.
Thank you again for all your help!
Mike
- Created iT2_pest.exe by compiling the following files: eospest.f, it2input.f, it2main.f, it2pest.f, it2stubs.f, it2user.f, it2xxxx.f, mdeppc.f, meshm.f, mio.f90, t2cg22.f, t2f.f, and t2solv.f
-
Mike,
Compilation seems ok.
But you need to have at least one template and at least one instruction file (you may define dummy input parameters and output variables to be written to and read from dummy input and output files, respectively).
Stefan
-
Stefan,
Thank you. After spending all of last Friday on the above effort, I realized that running ECBM through PEST won't solve my problem, because I want the added capabilities of iTOUGH2 (Pcap output, hysteresis, etc.). The added capabilities of ECBM (more advanced sorption models) are secondary to the other priorities. Therefore, I went back to using the publicly available iTOUGH2-EOS7C module this week.
Unfortunately, the default version of iTOUGH2-EOS7C does not print out dissolved methane mass fraction, which is a very important variable for my purposes. I tried my best to modify the source code to make it printout XCH4L instead of XTRC (because I'm not trying to model any tracers at the moment), but I don't think I've succeeded, considering that the printed values are still 0 everywhere, while I'm injecting methane over a long period of time and gas phase methane is present. I can share my adjusted source code with you if you'd like, but I'm not sure if it's okay to attach it in this public space, so please let me know if you'd like me to email it to you instead.
Since I have no experience coding in Fortran myself, it is quite difficult, if not impossible, for me to figure out how to solve my issue on my own. Would you mind telling me which line numbers of the eos7c source code I need to adjust in order to simply print out XCH4L values?
Thank you very much again,
Mike
-
Mike,
Go to file eos7c.f subroutine OUT. Add change as indicated below; actually, the code may already be there, so all you need to do is remove/add the comment character C after the line WRITE(66,5040):
IF(NK.EQ.5) XCH4G=PAR(NLOC2+NB+5)
IF(NK.EQ.5) XCH4L=PAR(NLOC2L+NB+5)
SG=PAR(NLOC2+1)
XB=PAR(NLOC2L+NB+2)
XNCG=PAR(NLOC2L+NB+3)
XTRC=PAR(NLOC2L+NB+4)
XNCGG=PAR(NLOC2+NB+3)
XTRCG=PAR(NLOC2+NB+4)
c
IF(SG.EQ.1.) XB=PAR(NLOC2+NB+2)
IF (IHYSTER.EQ.-1) THEN
WRITE(66,5040) ELEM(N)(:IELL),N,PRES,TX,PAR(NLOC2L+1),
& XB,XNCG,XTRC,XCH4G,XNCGG,XTRCG,PERM,
& PHI(N),ICURV(N),SOROLD(N)
ELSE
WRITE(66,5040) ELEM(N)(:IELL),N,PRES,TX,PAR(NLOC2L+1),
C & XB,XNCG,XTRC,XCH4G,XNCGG,XTRCG,PERM
& XB,XNCG,XCH4L,XCH4G,XNCGG,XTRCG,PERM
Stefan
