Hydrostatic Pressure Conditions in PetraSim/TOUGH

This is the first in a series of posts discussing the initialization of a TOUGH model to represent realistic natural state conditions.  In almost all cases, it is beneficial to run a model to steady state before adding transient boundary conditions, injection or pumping wells, etc.  We’ll give a very simple example first, where we’ve used EOS3 to establish a hydrostatic pressure gradient in a single phase model.

We’ve created a simple three dimensional model with dimensions of 1000x1000x100 meters.  We use a rectangular grid with spacing of 100m x 100m x 10m in the x, y and z directions.  The top of the model is 900 meters below ground surface. We’ll assume that all of the material above the top of the model is saturated with fluid of a constant density, so the top boundary of the model will be at a hydrostatic pressure based on a water column that is 900 feet thick.

The equation to calculate hydrostatic pressure is:

P = gρh,


P = Hydrostatic Pressure (Pa)
g = Gravity (m/s2)
ρ = Density of Water (kg/m3)
h = height (m)

The hydrostatic pressure calculated at the center of the top layer of cells in the model (a depth of 905 meters below ground surface) is 8.97E6 Pa.  We’ve also estimated a temperature of 47.625 °C based on a geothermal gradient of 25 °C/km with a surface temperature of 25 °C.  The global initial conditions for this model look like this:

Hydrostatic Pressure Global Initial Conditions

Figure 1: Global Initial Conditions based on Pressure and Temperature at a saturated depth of 905 meters.

We’ve made the top cell layer in the model “Fixed State”. By setting the top layer to “Fixed State”, these cells will remain at the pressure and temperature conditions defined in the Initial Conditions Dialog. The cells below the top layer will be at a higher pressure which will be calculated by TOUGH. Since this water will compress, the top layer of cells will act as a source to supply the required water.

Hydrostatic Pressure "Fixed State" cells

Figure 2: The top layer of cells are set as “Fixed State” cells and colored in red. These cells will have a constant pressure and temperature equal to the global initial conditions throughout the simulation.

We set the end time of the model to 100K years with a maximum time step of 100 years. TOUGH does a transient calculation. By setting a long time, the transient will reach steady state.

Hydrostatic Pressure Solution Controls

Figure 3: Solution Controls accessible through the Analysis menu.

The model reaches steady state conditions after around 1400 years and the resulting pressure gradient looks like the figure below.  Note that the resulting minimum pressure in the model is equal to the pressure at the top of the model.  The model cells below the top cell layer have equilibrated to higher pressure values based on hydrostatic pressure.

Hydrostatic pressure final pressure iso-surfaces

Figure 4: Final pressure iso-surfaces show a linear pressure gradient.

To load this pressure gradient as initial conditions in a new model, follow these steps:

  1. On the File menu, click Save and save a new .sim file in a different directory.
  2. With the new .sim file open, on the File menu, click Load Initial Conditions and choose the SAVE file created by the original simulation (stored in the same directory as the original SIM file).

Send comments or suggestions on this post to Alison Alcott at RockWare.

Download the PetraSim input files here:

Tips for Creating Contour Maps of Sparse Groundwater Elevation Data

Contouring sparse data in any mapping program can be challenging.  We’ve found a few tools in the RockWorks15 to be particularly helpful when creating contour maps of sparse groundwater elevation data.

First, let’s take a look at a contour map created using the EZ-Map tool in the RockWorks Utilities.

The EZ-Map tool uses simple triangulation to create contours.  For some data sets, this may be all you need to create a reasonable looking map.  However, it will often be necessary to create grid-based maps with smoother contour lines that extend to the edge of the project.  Note the increase in the groundwater elevation on the eastern edge of the map.  This is something that can be resolved by switching to a grid-based map.  We’ll explore grid-based mapping tools next.

Here is a map created using the default Inverse Distance Weighting settings.  I should note that I had both the “Smoothing” and “High-Fidelity” options turned on during grid creation.

The “bull’s eyes” that you see around the high and low points in the map are typical of the Inverse Distance interpolation method.  One way to resolve this is to decrease the number of points used during interpolation.

To do this, change the Number of Points used for the Inverse Distance Algorithm from 8 (the default setting) to 4 (which is more appropriate for a data set of this size).

Here is a map created with the modified “Number of Points” value:

The bull’s eye effect has been muted somewhat, but notice that the contours don’t honor the data extremely well.  Let’s move on to Kriging.

In the map below, I let the RockWorks program choose the appropriate variogram settings.  With Kriging especially, which can create fairly blocky models, I highly recommend that you turn on both the grid “Smoothing” and “High-Fidelity” options when
creating a contour map. 

This may be a little bit more to your liking, but the general groundwater flow direction could still be better represented along the borders of the map.  Just to cover all of our
bases, here is another map created using Triangulation gridding.  Unfortunately, there are some problems with edge effects in the resulting map as well.

None of these are really doing it for me.  At this point, I think that a lot of people would probably resort to hand drawing their contour maps, or adding additional control points to the data set to force the contours into the shape they have in mind.  Before you resort to these tedious and time-consuming options, I would recommend you look at the Densify and Polyenhancement options available in RockWorks15.

Here is a diagram showing the contour maps created with the Inverse Distance interpolation algorithm, with and without Densify turned on.  As you can see, the densification process (which adds additional control points to the data set before
interpolation using triangulation) straightens out the contour lines quite a bit.

I did the same test using the Kriging algorithm and got the following results.

Last but not least, here is a contour map created with the Polyenhancement option turned on.  When this option is activated, the program fits the data to a polynomial surface and then warps that surface to honor the data points (in this case, I choose a 2nd order polynomial surface).  I think I have my map!


In real life, I’ve found first, second and third order polynomials useful when creating groundwater contour maps.  If the groundwater flow direction is fairly constant through the area, go with a 1st order polynomial (which is a planar surface).  If it is variable because of topography or a feature such as a river or stream, then a 2nd or 3rd order polynomial is the way to go.



Computing Aggregate Reserves for a Site with Two Isolated Carbonate Units

This paper describes how to use RockWorks to compute total economic reserves for a site that includes two carbonate units: an upper limestone and a lower dolomite, separated by a shale unit. It involves creating separate I-Data models using the Stratabound filter, combining the models, and checking them against the observed log data.

Link to original paper: http://www.rockware.com/assets/products/165/casestudies/6/9/computing_aggregate_reserves.pdf



 The purpose of this study is to compute the total economic reserves for a site that includes two carbonate units; an upper limestone and a lower dolomite separated by a shale unit. Quality analyses have been obtained at one-foot intervals within the carbonates. The following diagram depicts a typical log showing the lithology, stratigraphy, and aggregate quality.

Figure 1: Typical log depicting aggregate quality (bargraph on left), stratigraphy (patterns in center), and lithology (subdivisions within stratigraphy)

Step 1. The Problem

Modeling the rock quality en-masse is problematic because the node values would include the quality values for both the limestone and the dolomite. The following diagrams depict a solid model based on the rock quality and a stratigraphic block model. Notice how the rock quality (I-Data) model interpolates quality values where there is no corresponding carbonate.

Figure 2: Problematic “Bulk” Rock Quality Model
Compare the rock quality model with stratigraphy model below and note how quality values are interpolated where there is no carbonate.
Figure 3: Stratigraphic Model

 Compare this stratigraphic model with bulk rock quality model above and note how quality values were interpreted within overburden (light yellow) and interburden.

Step 2. The Solution

The solution to this problem is to use the “Stratabound” option within the I-Data / Model menu. Two rock-quality models were created; one for the upper limestone and another for the lower dolomite.

In the example below, the I-data model is confined to points and nodes within the Hanford Limestone unit.

Figure4: Hanford Limestone Rock-quality Model

 In this example, the I-Data model is confined to points and nodes within Shuller Dolomite.

Figure5: Shuller Dolomite Rock Quality Model

Step 3. Combining the Models

The next step involved adding the two models together and removing all voxels with a quality value less than 50 (the minimum acceptable quality).

Figure6: Fence diagram depicting combined rock-quality models for upper limestone and lower dolomite.

Figure 7: Block Model depicting voxels with a quality value greater than 50.

Figure 8: Block model depicting zones from previous model in which the thickness for any single contiguous ore zone is more than 6 feet thick for any given column.

Figure 10: Block model depicting zones from previous model in which the stripping ratio is less than 1.2. This area represents a good place to start mining in order to gain the highest return on investment.

 Step 4. Checking the Model

The final, and most important step, is to create a 3D log diagram, combine it with the final ore model, and examine the data to see if it make sense.

Figure 11: 3-Dimensional Lithology/Quality Logs Combined With Final Ore Model.

Figure 12: Enlargement of area around highest-ROI ore depicting lithology and quality logs.

Step 5. Conclusion

By combining the preceding approach with increasingly more tolerant filter cutoffs, it is possible to create a mining strategy that will yield the highest return on investment from the onset.


Welcome to the RockWare Blog

In this space we’ll be posting occasional user tips, news, and information relating to RockWare, Inc., the Earth Science Software Company in Golden, Colorado, USA.  We welcome your comments and invite you to stay tuned.

We hope you will find this new resource useful. Thank you for your interest in the Rockware, Inc. Blog.

Kind Regards,
The RockWare Team