Back to MODFLOW Tutorial Table of
Contents.
Example Problem 4, A Water Supply
Problem
A fully worked-out version of Example 3 is in the "Examples\example3" directory. The files include fully worked out models for this example named "bound1.mmb" to "bound5.mmb" and an example of an exported grid.
One of the great advantages in using MODFLOW is that it provides an easy way to specify many sorts of boundary conditions. In this exercise, we will go through a series of models that illustrate how to use many of these boundary conditions. This exercise is based on Problem 14 of Andersen (1993). We will use very simple models to illustrate each of the boundary conditions. These boundary conditions are all very similar to one another mathematically and we will illustrate this fact by using different boundary conditions to get the same result. We will have five related models and we will have one of the following boundary conditions in each of the models
In order to have equivalent boundary conditions in all the models, the boundary conditions must apply to exactly the same cells in all the models. This is not usually the case in models of real places. Usually the precise cells used to define a boundary condition are not critical because they are only an approximation of a real boundary. A slight difference in the approximation, doesn't matter.
Before we spend much time setting up the models, lets take a look at how the first four of these boundary conditions differ from one another. In all of them, The amount of flow into or out of the aquifer is a function of the head in the aquifer. The simplest head-dependent boundary is the general head boundary. In it, you specify a boundary head and a conductance. If the cell to which the boundary is connected has a different head than the boundary head, water will flow out of or into the aquifer at a rate that is proportional to the difference in head.
![]() |
Figure 1. Relationship between head and the aquifer and boundary discharge for head-dependent boundary conditions in MODFLOW. (Modified from McDonald and Harbaugh, 1988.) |
The drain package is similar to the general-head boundary except that water can flow out of the aquifer into the drain but can not flow from the drain into the aquifer. This is because the drain will go dry once it no longer receives any water from the aquifer.
Like the drain package, the river package limits flow into the aquifer but unlike the drain package some flow into the aquifer can occur. However, if the head in the aquifer cell attached to the river drops below the base of the river, recharge of the aquifer from the river occurs at a constant rate. Thus, you must provide two elevations; the head in the river and the elevation of the river bottom.
In the evapotranspiration package, when the head in the aquifer cell from which evapotranspiration is occurring is above the evapotranspiration surface, evapotranspiration occurs at the maximum rate. If the head is below the cuttoff elevation evapotranspiration from the aquifer does not occur. However, you do not specify the cuttoff elevation directly. Instead, you give the evapotranspiration extinction depth which is the difference in elevation between the evapotranspiration surface and the cuttoff elevation. When the head in the aquifer cell is in between the evapotranspiration surface and the cutoff elevation, evapotranspiration is proportional to the difference in head between the aquifer cell and the cutoff elevation.
In a two-dimensional model, such as the ones we will use here, it is possible to use a prescribed head to model a general-head boundary condition. We will put the constant head cells beneath the cell for which we wish to have a general-head boundary. Then we will adjust the vertical hydraulic conductivity so that the vertical conductance is the same as the conductance that we would have specified with the general-head boundary.
The models will all have seven rows, and seven columns and they will all be MODFLOW-2000 models. All but the last model will have a single layer. The aquifer will be unconfined and there will be a well in the upper left hand cell.
The following characteristics will apply to all the models. The links are to topics in previous examples that show how to set each characteristic.
Parameter | Argus ONE Layer Name(s) | Numeric Value |
---|---|---|
Initial head | Initial Head Unit1 | 10.0 ft. |
Hydraulic conductivity | Hydraulic Cond Unit1 | 10 ft./day |
Aquifer top | Elevation Top Unit1 | 10 ft. (or higher) |
Aquifer base | Elevation Bottom Unit1 | -50 ft. |
Specific yield | Specific Yield Unit1 | 0.1 |
Grid spacing (uniform) | MODFLOW Domain Outline, MODFLOW FD Grid | 100 ft. |
Pumping rate | Wells Unit1 | -2500 ft.3/day |
Model type | transient | |
Stress period length | 365 days | |
Time steps | 20 | |
Time step multiplier | 1.2 |
When you set up the model, you should make sure that the well screen is completely inside the layer. Thus if you set the well screen from 10 to -50, you must make sure that the top of the model layer is greater than or equal to 10 and the bottom of the layer is less than or equal to -50. If you don't, the MODFLOW-PIE, will treat the top and bottom as if they are at the top and bottom of the geologic unit and will show a warning message. The default value for the top of the top model layer is 0. If you have a well that extends across several layers, you should calculate the pumpage from each layer as described in the documentation for the MODFLOW Wells package (McDonald and Harbaugh, 1988).
To enter most of these values, select "Edit | Layers..." or press the "Layers" button to show the Layers dialog box. Then select the Argus ONE Layer whose value you wish to set. In the lower half of the Layers dialog box, click on value of the parameter until a popup menu appears and then select "Expression" from the popup menu. This will show the Expression Editor. In the upper half of the Expression Editor, type the value that you wish to assign to the layer and click on the "Done" button in the Expression Editor. The value you entered should now be displayed under "Value" in the lower half of the Layers dialog box. You can now either select another Argus ONE layer and assign a value to it, or, if you are done, click on the "Done" button on the Layers dialog box. For more details, see the section on the Expression Editor in the previous example. The links for each parameter are to the locations in previous examples where those parameters were set. You can use the same techniques in this example.
On the "Geology" tab of the "Project Info" dialog box, set the interblock transmissivity to "arithmetic mean". This is one of four possible choices for interblock transmissivity.
The harmonic mean is most appropriate for confined or unconfined aquifers with abrupt discontinuities in transmissivity or hydraulic conductivity. |
The arithmetic mean is most appropriate for unconfined aquifers with uniform hydraulic conductivities. |
The logarithmic mean is most appropriate for confined aquifers with gradually varying transmissivities. |
The "arithmetic and logarithmic" mean is most
appropriate
for unconfined aquifers with gradually varying transmissivities. |
Set up a model with a size of 700 ft x 700 ft. Set the grid density to 100 and make the grid. You can set these in same way you did in Example 1. You can also create the grid manually. If you want a uniform grid, creating it manually acceptable. In general, however, I think creating a domain outline and using the "magic wand" to create the grid is a better method. To manually create a grid, go to the MODFLOW FD Grid layer and select the Grid Tool. Click and drag with the grid tool to create the grid.
In this example, I suggest that you create a domain outline and use the "magic wand" to create the grid. Add a well anywhere in the upper left hand cell and set its pumping rate to -2500 ft3/day. Set the remaining data in the table above using the same methods as in previous examples. Don't forget to select Well on the Packages|Stresses 1 tab of the Project Info Dialog box.
For this model it will be convenient to have precise control on the domain boundary. We want it to be as close to 700 feet long in each direction as we can get. The Utility PIE is useful for this purpose. After drawing a Domain outline that is close to what we want, we can use the EditContours to set the precise node positions.
If you have not already done so, draw a square Domain Outline and
make it 700 feet wide and 700 feet high. Then select "PIEs|Edit...|Edit
Contours..." Next select the MODFLOW Domain Outline layer from the list
of layers. The objects on the Domain Outline layer will be imported
into the EditContours PIE. You can click on any node there to select it
and edit it's position. (You can also copy a contour to the clipboard
and paste it into a text editor to edit it there.)
Once you create the grid, you can save the grid to a file and import it in another model. (You can do this with other types of layer too.) From the MODFLOW FD Grid layer select "File|Export|Export MODFLOW FD Grid". You will then be given a choice of what data to include in the exported grid.
In all the models, we will have a boundary condition running through the middle of the middle column. In the first model, that boundary will be a river.
To create a river boundary, you must first check the "Rivers" checkbox on the Packages|Stresses 1 tab of the Project Info dialog box. You can display the Project Info dialog box by selecting "PIEs | Edit Project Info..."
We will want to add a river boundary with the following characteristics.
Elevation (Stage Stress) | 0.0 ft. |
Width | 20 ft. |
Riverbed hydraulic conductivity | 0.1 ft./d |
Riverbed thickness | l ft. |
River bottom elevation | -2.0 ft. |
To set it up, go to the Line River Unit1 Layer and select the Open Contour Tool. If the Open Contour Tool is not already displayed, select the button beneath the arrow (or hand) button and hold down the mouse button until a pop-up menu appears. (See the illustration above.) Select the Open Contour Tool from the pop-up menu. It is the middle item. Now draw a vertical line through the middle column of the model and assign it a conductance of 2, a bottom of -2 and a stage stress of 0. If you hold down the "Shift" key while drawing the river, it will be easier to draw a vertical line. It is important that the River be in the middle column but it isn't critical exactly where in the middle column it lies. You will assign a conductance multiplier which is the conductance per unit length of the boundary. Thus, conductance multiplier, as used in the dialog box, is not the same as conductance as defined in the MODFLOW documentation. However, it is better to assign the conductance per unit length, rather than conductance, because that allows you to change the grid without having to re-enter data for the boundaries.
In MODFLOW, the riverbed hydraulic conductance is equal to KLW/M where
The MODFLOW-GUI measures the length of the reach (L) and multiplies it by the value you enter so the conductance you should enter is KW/M. For this example, the riverbed hydraulic conductivity is 0.1 ft./day, the river width is 20 ft. and the riverbed thickness is 1 ft. so KW/M is 2.0. The conductances that will be written to the MODFLOW input files will be 200 because the conductance multiplier you enter will be multiplied by the lengths of the cells which are all 100.
You can now run the model. When you are done, use GW_Graph to plot the head in the cell in row 1 column 4. Save the data to a text file an import into a spreadsheet and make a graph of head versus time. Your results should be similar to the ones under "River" in the following table.
It is easy to make most of the remaining models by just making a few simple changes in this model. First, save this model and then save it again with a new name. Copy the river contour to the clipboard. For the General-Head Boundary, go to the project info dialog box and check GHB and uncheck RIV. On the "Line General Head Boundary Unit 1" layer, paste the contour from the clipboard and set its conductance to 2 and the "Head Stress 1" to 0. Clear the remaining values.
To change the boundary condition to a drain, save the model and save it again with a new name. Then uncheck GHB and check DRN in the Project Info dialog box. You can paste the contour on the clipboard to the "Line Drain Unit 1" or just draw an open contour through the middle of the center column on the "Line Drain Unit 1" layer and set "Conductance" to 2, "Elevation" to 0, and "On or Off Stress 1" to True. (In all of these models "Conductance multiplier" is Conductance per unit length. The lengths of the cells are 100 ft. so the conductances in the MODFLOW input files will be 200.).
The model with evapotranspiration is slightly different. In it you must outline the middle column with a closed contour and assign 10 to the evapotranspiration surface, and evapotranspiration extinction depth. Assign 0.2 to EVT Flux Stress1. (In most models, you wouldn't try to outline particular cells but instead would outline the area affected by evapotranspiration. In this case, however, we are trying to use different packages to define equivalent boundary conditions. This requires that exactly the same cells be used to define the boundaries in all the models. Thus, we must be more concerned with the grid than we would be in a more normal situation.) This boundary condition is equivalent to the drain so long as the head is less than 10. MODFLOW will multiply the flux stress of each cell by the cell area to get the maximum volumetric flux rate. In this case that would be 0.2 x 100 x 100 = 2000. You would get this flux if the head was at 10 ft. If the head was at one foot, the flux would be 200. You would get that same flux with the other packages because the conductance calculated by Argus ONE for those other packages is also 200.
It will take a bit more work to set up the last model. First, add a confined layer below the unconfined layer on the Geology tab of the Project Info Dialog box. (Only the top layer can be unconfined. If any of the lower layers are unconfined for part or all of the model, they must be modeled as convertible layers.) Eventually, you will need to unselect evapotranspiration on the Packages|Stresses 1 tab too but don't do that right away because you can copy the closed contour from the evapotranspiration layer and use it in the new model. When you insert a new layer to a model, it is always inserted above the selected layer in the Geology tab of the Project Info dialog box. To add a layer below the bottom layer, click on the "Add" button button on the Geology Tab. The bottom layer should have one column of constant head cells in the middle. There are a couple of ways you can do this both of which are easy. You can use whichever method you prefer.
To get the right conductance between the upper and lower aquifers, you must assign correct values to the tops and bottoms of both aquifers and the vertical hydraulic conductivity (Kz) of both layers. Correct values are listed in the following table.
Geologic Unit | Argus ONE Layer Name | Numeric Value |
---|---|---|
Top, Upper Aquifer | Elevation Top Unit1 | 50 |
Bottom, Upper Aquifer | Elevation Bottom Unit1 | -50 |
Kz, Upper Aquifer | Hydraulic Cond Unit1 | 2 |
Top, Lower Aquifer | Elevation Top Unit2 | -50 |
Bottom, Lower Aquifer | Elevation Bottom Unit2 | -150 |
Kz, Lower Aquifer | Hydraulic Cond Unit2 | 2 |
The reason these values work is that the vertical conductance is
calculated
as .
where
Delta zu = the thickness of the upper layer,
Delta zl = the thickness of the lower layer,
Kzu = the hydraulic conductivity of the upper layer,
Kzl = the hydraulic conductivity of the lower layer,
Delta x = the length of the cell in the x-direction, and
Delta y = the length of the cell in the y-direction.
In this case, Delta x = Delta y = Kzu = Kzl = 100 ft. which makes the vertical conductance 200 just as it is in the general-head boundary.
If you compare the results of the five models, you will see that they are all similar at first. Until the head in Row 1, Column 4 drops below 0 ft., all the boundary conditions are identical in all the models so the head is the same in all the models. When the head drops below 0, the river, general-head boundary, and prescribed head boundary begin supplying water to the model whereas the drain and evapotranspiration boundaries do not. These two groups of models begin to diverge at that point (77.385 days). When the head drops below -2 (207.108 days), the amount of water supplied by the river boundary to the model changes to a constant rate and begins to diverge from the models with the general-head boundary, and the prescribed head boundary.
The similarity between the results of these models illustrates the similarity in the underlying mathematics. In the diagram at the beginning of this example, all of the boundaries had a sloping section. We set up the boundaries so that those sloping sections were the same in all the models. Once the model got out of those sloping sections, the results began to diverge. You can take advantage of that similarity to use the head-dependent flow boundary packages to model boundary conditions different from those implied by their names. The following table from Andersen (1993) illustrates some of the other uses for the head-dependent flow boundaries.
General-Head Boundary | Drain | River | Evapotranspiration |
---|---|---|---|
Rivers | Intermittent streams | Adjacent aquifers | Drains with maximum flow limitation |
Exterior model boundaries | Springs | Wetlands | |
Adjacent aquifers | Ditches |
If you choose the parameters for the head-dependent flow boundaries improperly, you can cause your model to become unstable. Avoid using elevation in the head-dependent boundaries that are lower than the bottom of the cell with the boundary. What can happen is that loss of water may be sufficient to cause one or more cells to go dry. When the cell goes dry, the head-dependent boundary is no longer active. That may allow the cell to rewet (if rewetting is enabled). One or more cells may cycle back and for between a wet and dry state indefinitely. If for some reason, you need to make them elevation of the head-dependent boundary less than the bottom elevation of the boundary, you should carefully examine your model output to make sure that the cells with the head-dependent boundaries did not go dry. Another possible cause of instability would be conductances in the head-dependent boundaries that are much higher than the conductances between cells for groundwater flow. If this is a problem, use a smaller value for the head-dependent boundary conductance.
In this set of examples, we have introduced head-dependent boundary conditions. Such boundaries are used in almost every model. They can be used for real-world boundaries that are different from those implied by their names. For example, a drain boundary can be used to model a spring.
All the head-dependent flow boundaries have some sort of elevation associated with them. If possible, you should try to avoid having that elevation lower than the elevation of the bottom of the cell in which the flow boundary occurs. If you don't, the cell may go dry which can cause problems with your model.
If you wish, you can create a grid manually. This is useful if you just want to create a uniform grid.
Finally, MODFLOW has a number of options for calculating the average hydraulic conductivities between cells. Choosing the correct option will slightly increase the accuracy of the simulation without greatly increasing the time required for a model to run.
The harmonic mean is most appropriate for confined or unconfined aquifers with abrupt discontinuities in transmissivity or hydraulic conductivity. |
The arithmetic mean is most appropriate for unconfined aquifers with uniform hydraulic conductivities. |
The logarithmic mean is most appropriate for confined aquifers with gradually varying transmissivities. |
The "arithmetic and logarithmic" mean is most
appropriate
for unconfined aquifers with gradually varying transmissivities. |
Example Problem 4, A Water Supply
Problem
Back to MODFLOW Tutorial Table of Contents.