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:

Descent From Space …

(128) New Feature (03/14/16/JPR): A new option titled “Start With Descent From Outer Space” has been added to the RockWorks17 | Utilities | EarthApps | Flyover | Advanced: Circular program menu.

This new capability will automatically add a virtual introduction to the flyover such that it starts in outer space.  The “camera” position starts out by pointing towards the opposite side of the Earth relative to the project area.  The Earth then spins around such that the project area is directly below.  The “camera” then descends.  As it descends, it gradually tilts to the same angle as that specified for the first circular fly-over.
 Unlike the standard Google Earth fly-in (i.e. clicking on a waypoint), the descent does not appear to exponentially accelerate with the descent.  Instead, a geometric progression is used to slow the descent down as the camera moves closer to the project area.  Then net result is the illusion of a uniform descent thereby providing the viewer with a better spatial orientation.

Click on the link below to view a demonstration of this new capability.


Identifying Spatial Correlations Between Grid Models (Comparing Apples & Oranges)

(127) New Feature (03/04/16/JPR): A new program titled “Correlate Grids” has been added to the Utilities | Grid | Statistics menu.

This program is used to identify spatial correlations between grid models that don’t necessarily represent similar units (e.g. ppm versus occurrences per cell).  Stated euphemistically; it’s a way to compare apples and oranges.

grid_correlations_02Here’s how the program works:

  • For each of the input grids (lead concentrations in soil versus birth defects within the fictitious example shown below), the standard deviation and mean are computed for all of the grid nodes.

  • The grid nodes are then recomputed to represent the number of standard deviations from the mean in order to quantify just how anomalous a node value is relative to the entire node population.

    • The following terminology (from mineral exploration) may help in understanding what these new values represent:

      • Background: Values that are less than 1 standard deviation from the mean.

      • Weakly Anomalous: Values that fall between 1 and 2 standard deviations from the mean.

      • Moderately Anomalous: Values that fall between 2 and 3 standard deviations from the mean.

      • Strongly Anomalous: Values that fall between 3 and 4 standard deviations from the mean.

      • Extremely Anomalous: Values that are greater than 4 standard deviations from the mean.

      • Outrageous Xtreme Tactical Ultra-Anomalous: OK, we make this one up.

    • Note: During this process, negative anomalies are converted to zero.  It is assumed that the user is only interested in positive anomalies.

  • On a node-by-node basis, the absolute value of the differences between the node values are subtracted from the higher value.  For example, let’s say that a lead node has a value of 3.6 while the corresponding birth defects node has a value of 3.1.  The resulting correlation coefficient would be 3.1 ( 3.6 – abs ( 3.6 – 3.1) ).  That node would therefore be considered as a point where both the lead and defects are strongly correlative in a spatial sense.

    • In the diagram below, the spatial correlations,  shown by solid-filled color contour bands, highlight an area where the lead and birth defects are weakly correlative.


Caveats / Disclaimer:

  • This approach assumes that both data sets encompass a large enough area and sample population / distribution to establish the “background” levels for each component.

  • This results do not in any way, shape, or form, establish causality – just spatial correlation.  Indeed, many spatial correlations are often purely coincidental and misuse of this methodology will inevitably result in professional ridicule, humiliation, and sexual dysfunction.