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

  1. 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.
  2. 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
  1. 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.
  2. 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.
  3. Select PHAST Options|Free Surface… and check the Use free surface check box.  Click the OK button to close the dialog box.

Data Sets

  1. 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
  1. 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.
Object used to set inactive area on top view of the model.
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

  1. 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
  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.
  2. 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.
  3. 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.)”.
  4. 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.)”.
  5. Try coloring the grid (Data|Color Grid…) with the River Heads (fig. 45).  (In this illustration, the rivers objects have been colored red.)
The top view of the model with grid cells colored to reflect the river head in example 4.
Figure 45. The top view of the model with grid cells colored to reflect the river head in example 4.

Specified Flux Boundaries

  1. 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.
Polygon used to define specified flux boundary in example 4.
Figure 46. Polygon used to define specified flux boundary in example 4.
  1. 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.
Polygon used to specify the specified head boundary condition in example 4.
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.

  1. 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.
Polygon used to specify the leaky boundary on the front view of the model in example 4.
Figure 48. Polygon used to specify the leaky boundary on the front view of the model in example 4.
  1. 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.
  2. 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.
Point object used to specify well boundary in example 4.
Figure 49. Point object used to specify well boundary in example 4.

PHAST Options

  1. 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.
  2. 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.
  3. 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.
  4. 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.
  5. 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.