Example 4
Example 4 is described beginning on page 103 of the PHAST
documentation (Parkhurst
and others, 2004). The flow and
transport data file itself is reproduced beginning on page 107.
If PHAST is installed on a Windows computer at the default location,
this file can be found at
C:\Program; Files\USGS\phast-1.1\examples\ex4\ex4.trans.dat.
Initial Set Up
- Start GoPhast as before. In the initial
grid dialog box, set the number of nodes in the X, Y, and Z directions
to be 16, 9, and 5 respectively. Set the distance between nodes
in the X, Y, and Z directions to be 6000, 6000, and 100
respectively. Click the Finish button to close the dialog box.
- In the main GoPhast window select PHAST
Options|Title and Units… Copy the title from ex4.trans.dat to
the
clipboard and paste it in the title section. Assign the following
units. Then click the OK button.
Option
|
Units
|
Time units
|
years
|
Horizontal grid units
|
meters
|
Vertical grid units
|
meters
|
Head units
|
meters
|
Hydraulic conductivity units
|
meters per second
|
Specific storage units
|
1/meters
|
Dispersivity units
|
meters
|
Flux units
|
meters per year
|
Leaky hydraulic conductivity
units
|
meters per second
|
Leaky thickness units
|
meters
|
Well diameter units
|
inches
|
Well flow rate units
|
liters/day
|
River bed hydraulic conductivity
units
|
meters per second
|
River bed thickness units
|
meters
|
- The next step is to activate the chemistry
options. Select PHAST Options|Chemistry Options…
Check the Use solute transport check box.
Uncheck the Use
gas phases, Use
solid solution, and Use kinetics check
boxes. Click the OK button.
- Set the Steady Flow options by selecting PHAST
Options|Steady Flow… and checking the Steady flow check
box.
Change the Head tolerance to 1E-6 and click the OK button.
- Select PHAST Options|Free Surface… and check
the Use free surface check box. Click the OK
button to close the
dialog box.
Data Sets
- Select Data|Edit Data Sets… and enter the
following values for the default formulas.
Data Set
|
Default Formula
|
Kx
|
1.5E-5
|
Ky
|
1.5E-5
|
Kz
|
1.5E-7
|
Porosity
|
0.22
|
Specific_Storage
|
0
|
Longitudinal_Dispersivity
|
2000
|
Horizontal_Transverse_Dispersivity
|
50
|
Vertical_Transverse_Dispersivity
|
50
|
Initial_Head
|
380
|
Chemistry_Initial_Solution
|
2
|
Chemistry_Initial_Equilibrium_Phases
|
2
|
Chemistry_Initial_Surface
|
2
|
Chemistry_Initial_Exchange
|
2
|
- In this example, the bottom layer to the right of
an X coordinate of 48000 is inactive. To set this, first draw an
object on the top view of the model enclosing this area (fig.
44). In the Object Properties dialog box, leave the Elements
radio button checked, set the object’s name to “Inactive_Area”, check
the Set values of enclosed elements check box and set
the Higher
Z-coordinate to 100. On the Data Sets tab, check the
check box
for Active and set the formula to “False”.
After closing the
Object Properties dialog box, color the grid with the Active
data set
and set the selected layer to the bottom layer. Check that the
inactive area is restricted to the bottom layer (fig. 44). Using
the mouse cursor, right click on the top view of the model and select Hide
to hide the “Inactive_Area” object.
 |
Figure 44. Object
used to set inactive area on top view of the model. The grid has
been
colored with the Active data set showing
that the inactive area is restricted to the bottom layer. |
Rivers
- This model has three rivers: the Little River,
the North Fork River and the North Canadian River. The Little
River is defined by a polyline object with vertices at (44000. 15000.),
(44000. 0.), and (90000. 0.). Draw a polyline object on the top
view of the model at approximately those locations. (When drawing
the object, it is possible to see the current location of the cursor in
the status bar.) In the Object Properties dialog box, select the Nodes
radio
button. Change the name of the object to “Little_River”. Check
the Set values of intersected nodes check box.
Then go to the
Vertices tab and correct the locations of the vertices to be exactly
correct.
On the Boundary Conditions tab, select a
River boundary, and set the
following properties.
River Property
|
Value
|
Name
|
Little River
|
Hydraulic Conductivity
|
1
|
Width
|
200
|
Depth
|
1
|
Bed Thickness
|
1
|
- In the table for the river boundary set the Associated
solution to 1. The head must vary
along the length of
the river. A formula can be used to accomplish this. Click
on the cell in the table for head and observe the F() button that
appears. Click it. This will open the Formula Editor.
- The head in Little River should vary from 335
at the upstream end to 275 at the downstream end. The
“Interpolate” or “MultiInterpolate functions can be used to achieve
this goal. On the right side of the Formula Editor, click on the
“+” next to Functions to display the types of functions that are
available. Then click on the “+” next to “Math” to display all
the math functions. Click on MultiInterpolate in this list.
Then click the Function Help button. The
default web browser will
open and display help on the MultiInterpolate function. Read the
description of the MultiInterpolate function to learn what it
does. Go back to the Formula Editor and double-click on
“MultiInterpolate”. Note that MultiInterpolate appears in the
formula text box along with its parameters. Go back the Data Set
and Function List and click on the “-” next to Math and then click on
the “+” next to Object. Select “FractionOfObjectLength” and click
on the Function Help button again. Read the description of what
FractionOfObjectLength does. FractionOfObjectLength will be used
as the first parameter in MultiInterpolate in the Formula text box of
the Formula editor. Select “Position” in the Formula text box and
double-click on FractionOfObjectLength. Next substitute 335 and
275 for Value1 and Value2 respectively. Substitute 0 and 1 for
Distance1 and Distance2 respectively. Remove the extra square
brackets and dots so that the final formula is
“MultiInterpolate(FractionOfObjectLength, 335., 0., 275., 1.)”.
Click OK to close the formula and then click OK again to close the
Object Properties dialog box.
- The North Fork River is entered in the model
in a similar fashion. The coordinates of its vertices are (30000,
36000), (30000, 48000), and (90000, 48000). Its Hydraulic
Conductivity, Width, Depth, Bed Thickness, and Associated Solution are
all the same as in
the Little River. The heads at the upstream and downstream ends
are 335 and 280 respectively so the formula for head is
“MultiInterpolate(FractionOfObjectLength, 335., 0., 280., 1.)”.
- The North Canadian River is entered in the
model in a similar fashion. The coordinates of its vertices are
(60000, 30000), and (90000, 20000). Its Hydraulic Conductivity,
Width, Depth, Bed Thickness, and Associated Solution are all the same
as in the Little
River. The heads at the upstream and downstream ends are 350 and
305 respectively so the formula for head is
“MultiInterpolate(FractionOfObjectLength, 350., 0., 305., 1.)”.
- Try coloring the grid (Data|Color Grid…) with
the River Heads (fig. 45). (In this illustration, the rivers
objects have been colored red.)

|
Figure 45. The top view of the
model with grid cells colored to
reflect the river head in example 4. |
Specified Flux
Boundaries
- The next boundary condition is a specified
flux boundary. It occupies the entire top layer to the right of
30000 except for the first and last rows. Draw a polygon that
encloses this area (fig. 46). In the Object Properties dialog
box, select the Nodes radio
button. Set its name to “Flux_Boundary”. Check the Set
values of enclosed nodes check box. Set the number of
associated
third dimension formulas to one and set Z-coordinate to 400. On
the Boundary Conditions tab, choose a Flux boundary. Set the Flux for
time 0 to –0.055 and set the associated solution to 1. Click OK
to close the dialog box. It is possible to try coloring the grid
with Top_Flux_Boundary_Flux to make sure the boundary condition was
applied properly.

|
Figure 46. Polygon used to
define specified flux boundary in example
4. |
- The next boundary condition is a specified
head boundary (Lake Stanley Draper, figure 47). On the top view
of the model, it looks like the lake only affects one node. However, it
actually extends from an elevation of 300 to 400. Thus, it
affects two nodes: one in the top layer and one in the layer just
beneath it. Draw a polygon surrounding the node at (30,000,
18000) on the top view of the model. In the Object Properties
dialog box, select the Nodes radio
button. Set its name to “Lake_Stanley_Draper”. Check the Set
values of enclosed nodes check box. Set
the Higher
Z-coordinate to 400 and the Lower Z-coordinate to 300. On the
Boundary Conditions tab, select a specified head boundary and set the
head to 348 and the associated solution to 1. Click OK to close
the dialog box.

|
Figure 47. Polygon used to
specify the specified head boundary
condition in example 4. |
Leaky Boundaries
There are two leaky boundaries in the model. Both are in the
XZ plane so they are drawn on the front view of the model. One is at
the front of the model and extends from the left edge of the model to
an X-coordinate of 39000. The other is at the back of the model
and extends from the left of the model to an X-coordinate of
29000. Both leaky boundaries extend over the whole extent of the
model in the Z direction.
- For the front boundary draw a rectangle on the
front view of the model that goes from the left edge to approximately
39000 and that covers the entire height of the model (fig. 48).
In the Object Properties dialog box, select the Nodes radio
button. Set its name to
“Front_Leaky_Boundary”. Check the Set values of enclosed
nodes
check box. Set the Associated third-dimension formulas
to one.
Set the Y-coordinate to 0. On the Boundary Conditions tab, select
the leaky boundary, set the hydraulic conductivity to 1.5e-5, the
thickness to 20000, the head to 320 and the associated solution to
2. Click OK to close the dialog box.

|
Figure 48. Polygon used to
specify the leaky boundary on the front
view of the model in example 4. |
- For the back boundary draw a rectangle on the
front view of the model that goes from the left edge to approximately
29000 and that covers the entire height of the model. In the
Object Properties dialog box, select the
Nodes radio button. Set its name to “Back_Leaky_Boundary”.
Check the Set values of enclosed nodes check
box. Set the Associated third-dimension formulas
to one. Set
the Y-coordinate to
48000. On the Boundary Conditions tab, select the leaky boundary,
set the hydraulic conductivity to 1.5e-5, the thickness to 30000, the
head to 305, and the associated solution to 1. Click OK to close
the dialog box.
- The next boundary condition is a well located
at (X, Y) = (12000, 36000). Make a point object on the top view
of the model at that location
(fig. 49). In the Object Properties dialog box, select the Nodes
radio button and name the object
“Well”. Check the Set values of intersected nodes
check
box. On the Boundary Conditions tab, select the well boundary,
set the Diameter to 2, the First elevation
to 90, the Second elevation
to 110, the Pumping rate to 1, the and the Solution
to 1. Make
sure the Allocate by pressure and mobility check box
is
unchecked. Click OK to close the dialog box.

|
Figure 49. Point object used to
specify well boundary in example 4. |
PHAST Options
- The Solution Method options will be left at
the default values. Select PHAST Options|Solution Method…
Verify that an Iterative solver is selected, Space
differencing is 0, Time differencing is 1,
and the Tolerance
is 1e-10. Click OK to
close the dialog box.
- Next, set the Print Initial options.
Select PHAST Options|Print Initial... Check the
Steady flow
velocities, XYZ heads, and XYZ
steady flow velocities check
boxes. Click OK to close the dialog box.
- Next set the Print Frequency. Select PHAST
Options|Print Frequency... Change XYZ
chemistry to “50000
years”, HDF chemistry to “2000 years”, and XYZ
wells to “2000
years”. Change the rest to “0 default” Check the Save
final head
in *.head.dat check box at the lower left. Click OK to
close the
dialog box.
- To set the time control, select PHAST
Options|Time Control... and set the Time step length
to 2000 and the Ending time to 100000. Click OK
to close the
dialog box.
- The last step is to turn off the printing of
the chemistry for the bottom layer. To do this, draw a polygon
entirely surrounding the grid on the top view of the model. In
the Object Properties dialog box, select the Nodes radio button.
Set the Name to “Deactivate_Printing”. Check the Set
values of
enclosed nodes check box. Set the number of associated
third
dimension formulas to one and set Z-coordinate to 0. On the Data
Sets tab, check the check box for Print_XYZ_Chemistry and change the
formula for it to “0”. Click OK to close the dialog box.
It is now possible to export the transport data file as
ex4.trans.dat,
run the model and view the results as with the previous examples.