MODFLOW Model of Groundwater Flow in Sloping Surficial Sand and Gravel
By Darrel Dunn, Ph.D., PG, Hydrogeologist
(View Résumé 🔳)
Purpose and Scope of MODFLOW Webpage
The purpose of this non-technical webpage is to present a MODFLOW model of groundwater moving downslope in surficial sand and gravel over impervious bedrock. This type of groundwater system is of interest because of the common occurrence of contaminants moving downslope from sources on sloping ground toward streams and associated alluvial aquifers. Modeling groundwater flow in such systems is relatively difficult because the saturated zone between the water table and bedrock is often thin. This condition tends to cause failure of the MODFLOW modules that solve the system of equations describing the flow. The tendency to fail is due to repetitive drying and wetting of portions of the model during the iterative solution process. Instability is produced in the MODFLOW solver modules. Nevertheless, a successful model may be developed by experimentation with input variables that control the behavior of the solvers.
This page uses non-technical language plus some technical language that is explained in the text . It also contains endnotes that provide technical details for readers who are interested. MODFLOW is a finite-difference groundwater modeling computer program maintained by the United States Geological Survey (USGS). A finite-difference model divides the groundwater system into rectangular cells, applies an equation that balances inflow, outflow, and change in water content for each cell, and solves the equations simultaneously to get hydraulic head in each cell. Hydraulic head at a point in a groundwater system may be thought of as the water level in a well if it got its water at that point. Groundwater moves in the direction of the head gradient, from higher to lower head. Thus MODFLOW can provide the rate and direction of groundwater movement.
MODFLOW Model
Figure 1 shows the topography of the area simulated in the MODFLOW model and the grid used in the finite-difference computations. The topography is shown as the elevation of the top of each grid cell. This is a sloping one-layer model. The layer represents the surficial sand and gravel. The thickness of the surficial layer is shown in Figure 2. This is also the height of the model cells. The water table is within this layer and the saturated part of the sand and gravel varies in thickness, as shown in Figure 3. Model cells west of the stream are inactive because all of the groundwater enters the stream rather than passing below it. There is no need to extend the active part of the model across the stream.
Figure 1. Elevation of MODFLOW cell tops (land surface above sea level) and model grid1.
Figure 2. Thickness of sand and gravel (cell heights).
Figure 3. Saturated thickness of surficial sand and gravel. Three white cells near the east boundary are dry.
The measure of the capacity of the sand and gravel deposits to conduct water that is used in the model is called hydraulic conductivity. Subsurface material with relatively high hydraulic conductivity, such as sand and gravel, conducts water more easily than silt and clay. Figure 4 shows the distribution of hydraulic conductivity in the model. Typical ranges of hydraulic conductivity for types of sand and gravel deposits are shown next to the colorbar. These ranges are logarithmic. Successive integer powers of 10 are an equal distance apart (1, 10, 100, et cetera; 100, 101, 102, et cetera). The number 1e+05 is common computer notation for 105. The hydraulic conductivities in the model were obtained by calibration. Calibration involved making a model run using hydraulic conductivity estimated from test well data2, then modifying this distribution of hydraulic conductivity via successive model runs until the simulated heads were reasonable3. All simulated heads are within 10 feet of measured heads and are below the land surface.
Figure 4. Hydraulic conductivity distribution in the MODFLOW model.
The model includes other hydrologic input besides hydraulic conductivity. Groundwater recharge due to infiltration of precipitation was set at 3.1 inches per year as the first step in the model calibration4. The climate is semi-arid and the land surface elevation is greater than 5000 feet above sea level. Groundwater inflow at the north (upgradient) edge of the model was set at 5 cubic feet per day. The elevation of the stream in the southwest part of the model was treated as a constant head, and a constant head was applied at the southeast corner of the model to account for water flowing toward the portion of the stream outside of the model boundary. Other boundaries were treated as no-flow (no water crosses the boundary)5.
There is uncertainty regarding reasonable values for groundwater recharge from precipitation in a semi-arid environment. The value of 3.1 inches is in the upper part of the range of published estimates. However. once the model was calibrated at 3.1 inches per year, it became relatively easy to find a lower value providing the same calibration. This possibility results from the fact that there is also uncertainty regarding the hydraulic conductivity values, and the two values are linked. One can reduce the the hydraulic conductivity values for all cells and groundwater recharge rate by the same multiplier and the heads in the model remain unchanged. Therefore, model calibration is not affected. However, the volume of groundwater flowing downslope toward the stream is affected. When the recharge rate is 3.1 inches per year the groundwater discharge to the stream derived from recharge and inflow at the northern boundary is 1,195,682 cubic feet per day. When the recharge rate is 1.55 inches per year, this discharge is halved to 597,341 cubic feet per day.
Conclusion
Figures 1 through 4 illustrate the use of MODFLOW to simulate groundwater flow in surficial sand and gravel on impervious bedrock sloping toward a stream. Solver difficulties encountered in this type of system can be overcome by experimentation with solver input variables. Output files from this model are used to demonstrate groundwater flow path simulation on a separate webpage, which also contains a map of the water table produced by the simulation.
Technical Endnotes
1 This is the numpy array for the model grid. Numpy indices start with zero. MODFLOW grids start with one. The grid size is 1000X1000 feet, which is small enough to capture variation in hydraulic conductivity in the surficial deposits and put observed heads in separate cells. Numpy arrays were used because the model was set up and run using USGS FloPy software to produce input files.
2I used specific capacity data and static water level data to get an initial value for hydraulic conductivity in cells containing the wells. Then I used computerized mapping to get an initial distribution of hydraulic conductivity for the model grid.
3The first set of calibration runs successively changed hydraulic conductivity near the well with the greatest head deviation until all differences between observed and measured heads were less than 10 feet. This amount of deviation allows for steady state modeling of heads that vary with time and were sampled at different times in each well. A second set of calibration runs was required to reduce some heads below the land surface. Again, I started with the greatest head deviation and worked down until all heads were below the land surface. A third set of runs was required to reduce differences between measured and computed heads that were affected by the second set. Calibration statistics are as follows:
Number of observed heads = 29
Sum of head deviations = 68.8 (observed head - calculated head)
Mean error = 2.37
Mean absolute error = 4.38
Root mean square = 5.14
RMS ratio = 0.0098 (Root mean square divided by the maximum difference in observed head across the model.)
There is no obvious clustering of positive or negative head deviations.
4Recharge was initially set at 1.94 inches per year based on literature review of recharge in a semi-arid climate at high elevation. This value resulted in low heads throughout the model and numerical instability related to cells cycling wet/dry during the solution. Increasing recharge to 3.1 inches, representing 19 percent of the mean local precipitation rate, is at the high end of reasonable recharge rates for the climate. Since recharge and hydraulic conductivity are linked in the basic equation, one can simulate lower recharge rates after calibration is achieved by simply multiplying recharge and hydraulic conductivity by the same factor.
5No-flow boundaries were close enough to the expected groundwater gradient and far enough from the central area of interest that they would not significantly affect the objective of the model.
Posted May 17, 2019
Revised September 16, 2019