MODFLOW AEM Comparison (NT)

Comparison of a MODFLOW Model with an Analytic Element Model of Mine Dewatering

By Darrel Dunn Ph.D., PG, Consulting hydrogeologist (Professional Synopsis)


Purpose and Scope

The purpose of this webpage is to present a MODFLOW groundwater model that demonstrates the use of dewatering wells for underground mine dewatering, and compare it to an analytic element model (AEM) of the same system. The comparison of the two modeling techniques is limited to this simple dewatering problem. Extrapolation of findings to other modeling scenarios is not intended. The webpage uses relatively non-technical language with some technical footnotes and references. The AEM model of the dewatering system is described on another webpage (mine dewatering). The two models simulate dewatering of a sub-surface mine with two tunnels to ore veins.

MODFLOW and AEM modeling

MODFLOW is a modular three-dimensional finite-difference1 groundwater modeling computer code maintained by the United States Geological Survey (USGS). Finite-difference groundwater modeling was introduced in 1967 (Early Groundwater Model), and more sophisticated models were developed in subsequent years. The technique became widely used after MODFLOW was developed in 1984 (reference 145). It has been evolving into a comprehensive groundwater modeling program as additional modules that expand its capabilities are added. Analytic element modeling began evolving in the 1980s (reference 146) and is not as widely used as finite-difference modeling. A brief description of AEM is provided on the webpage titled Analytic Element Groundwater Modeling. Whereas AEM uses relatively simple mathematical equations to compute the effect of such elements as water wells on groundwater systems; a finite-difference model divides the system into rectangular cells, applies an equation to each cell and solves the equations simultaneously to get the water table elevation in each cell. AEM models are generally used for relatively simple groundwater systems, whereas finite-difference models like MODFLOW can simulate very complex systems. AEM equations apply to extensive areas with no definite boundaries, whereas the extent of a finite-difference model is limited by defined boundaries. Conditions at the boundaries of a finite-difference model must be specified. Boundary conditions may be (1) impermeable boundary, (2) constant head2 boundary, or (3) general head boundary3. Other differences between AEM and finite-difference models exist. Consequently, a comparison of results of AEM and finite-difference models of the same dewatering system is of interest.

Method

FloPy was used to generate the input files for the MODFLOW model. FloPy is a set of Python programs developed by the USGS that can be used to produce input files for various USGS groundwater models, including MODFLOW. Python is a computer programming language that was introduced in 1991 and has evolved into a powerful tool increasingly used in science and technology. It includes matrix manipulation modules that can be used for development of cell by cell MODFLOW input and modules for plotting maps and graphs of model results.

MODFLOW Model Setup - Constant Head Boundary Model

For this MODFLOW model, the gridded area is the same as the area on the maps used to present the results of the AEM model. The grid is shown in Figure 1 superimposed on the pre-mining water table map from the AEM mine dewatering page.

MODFLOW grid for mine dewatering simulation

Figure 1. MODFLOW grid for mine dewatering simulation with AEM initial water table contours.


The grid cells are 200 feet square in the map view. The model contains only one layer and the aquifer is treated as horizontal with the bottom elevation at 8900 feet and top elevation at 10600 feet. So the cells are square columns 1700 feet high.

The boundary conditions are shown in Figure 2. All of the boundary cells (blue) are specified as constant head (fixed water table elevation) taken from the AEM model initial water table map. The northwest boundary is the approximate elevation of the stream. Cells north of the stream are inactive because they are not part of the groundwater flow system being simulated. The hydraulic conductivity (permeability4) assigned to the active cells is also from the AEM model. It is constant, meaning that the aquifer is treated as homogeneous.

MODFLOW Initial Water Table

The initial water table for the MODFLOW model was obtained by running the model without including the drainage wells. It is shown in Figure 2. The initial water table contours of the MODFLOW model differ from those of the AEM model. This difference is because the MODFLOW model includes the effect of the water table elevation on the ability of the aquifer to transmit water toward the stream. As the water table declines toward the stream, the saturated thickness of the aquifer decreases. As the saturated thickness decreases, the ability of the aquifer to transmit water decreases. The AEM model does not simulate this effect. Consequently, the water table contours in the MODFLOW model become more closely spaced near the stream. This closer spacing reflects a greater head gradient needed to counter the reduced saturated thickness of the aquifer. The greater head gradient translates to higher initial water table elevations at the tunnel locations. For example, the initial water table at the northwest end of the lower tunnel is about 20 feet higher in the MODFLOW model than in the AEM model.

MODFLOW boundary conditions and initial head contours.

Figure 2. MODFLOW boundary conditions, initial water table (head) contours, and mine tunnel locations. Constant head boundary cells are blue. Inactive cells are black. The lower tunnel is red. The upper tunnel is blue.


MODFLOW Dewatering Simulation

Figure 3 shows the result of adding the two dewatering wells to the MODFLOW model. This figure corresponds to the fourth figure on the AEM model page (Mine Dewatering). The extraction rates for both wells are the same as in the AEM model. The water table is lowered to below 9900 feet at the well near the southeast end of the lower tunnel, but the well is offset about 100 feet because wells are place in the middle of cells in MODFLOW. Wells can be more accurately positioned in an AEM model. The water table at the northwest end of the lower tunnel is about 9927 feet in the MODFLOW model compared to 9920 feet in the AEM model. The elevation of the tunnel is 9930 feet.

The water table at the north end of the upper tunnel is about 10088 feet in the MODFLOW model compared to about 10035 feet in the AEM model. The water table elevation near the center of the tunnel in the MODFLOW model is 10115 feet. The elevation of the tunnel is 10100 feet. So the MODFLOW model does not show complete dewatering of this tunnel, whereas the AEM model does.

MODFLOW water table produced by two dewatering wells.

Figure 3. MODFLOW water table elevation produced by two dewatering wells.


MODFLOW Mixed Boundary Simulation

A second MODFLOW model was set up with the east and west boundaries changed from constant head to impermeable. The purpose of this model was to see how much effect reducing inheritance from the AEM model would have on the water table near the tunnels. Figure 4 shows the effect on initial water table. The blue contours are the mixed boundary result (impermeable lateral boundaries), and the red contours are from Figure 2 (constant head lateral boundaries). The difference in simulated water table elevations at the boundaries is as much as 100 feet, but the difference in the internal part of the model near the tunnels is much less. Indeed the contours intersect (water table elevation equal) in the central part of the model. The difference at the boundaries results from the fact that groundwater flow near a boundary must be parallel to it and the contours must approach the boundary at a 90 degree angle.

Comparison of initial water tables produced by MODFLOW

Figure 4. Comparison of initial water table produced by MODFLOW for constant head boundaries (red) and mixed boundaries (blue).


Figure 5 shows the result of adding the two dewatering wells to the mixed boundary model. The extraction rates are the same as in the two previous models. The water table is still lowered to below the 9900 feet at the well near the southeast end of the lower tunnel. The water table at the northwest end of the lower tunnel is about 9938 feet (compared 9927 feet in the constant head boundary model). These two water table elevations are calculated as the average of 4 cells whose corners are at the adit.

The water table at the north end of the upper tunnel in the mixed boundary model is about 10074 feet (compared to 10088 feet in the constant head boundary model). The water table near the center of the tunnel in the mixed boundary model is 10099 feet (compared to 10115 in the constant head boundary model). The water table is slightly below the tunnel (10100 feet) in the mixed boundary model.

Pumping water table in mixed boundary MODFLOW model.

Figure 5. MODFLOW water table elevation produced by two wells in the mixed boundary model.


MODFLOW Recharge Simulation

A third MODFLOW model was set up with the east, west, and south boundaries set to impermeable and groundwater recharge added. The south impermeable boundary now simulates a groundwater divide, and simulated recharge of the groundwater by rainfall and snowmelt maintains a water table sloping toward the stream. The concept of groundwater divide is illustrated in Figure 6. Recharge causes the water table to rise between adjacent streams, while groundwater drains toward the streams. This condition results in a water table configuration that follows the land surface. Thus, groundwater divides tend to be below topographic divides. At the divide, groundwater flows downward than toward a stream and does not cross the divide. Hence, a groundwater model can treat the divide as an impermeable boundary because no groundwater crosses it. However, the current model contains only one layer of cells, so the downward flow components are not realistically simulated. Nevertheless, this simulation is more realistic than the previous ones, because the previous ones do not explicitly include groundwater recharge. Instead they cause groundwater to enter the system at the south boundary, which is really a groundwater divide.

Diagram of a groundwater divide.

Figure 6. Diagram illustrating the concept of groundwater divide. The squiggly arrows depict groundwater recharge.


Figure 7 shows the initial water table with the south boundary changed to impermeable and groundwater recharge added. The recharge rate is 4.6 inches per year. This rate causes the initial water table elevation at the tunnels to be about the same as in previous simulations, but the water table is much lower in the south part of the model near the groundwater divide.

MODFLOW initial water table.  Recharge scenario.

Figure 7. MODFLOW initial water table produced by three impermeable boundaries and groundwater recharge.


Figure 8 shows the result of adding the two dewatering wells to the recharge model. In this simulation, the extraction rates of the dewatering wells are changed to produce water table elevations more than 20 feet below the tunnel elevations. The new extraction rates are 52 gpm for the well near the end of the lower (north) tunnel, and 31 gpm for the well near the upper (south) tunnel.

MODFLOW water table elevation.  Recharge model.  2 wells.

Figure 8. MODFLOW water table elevation produced by two wells in the recharge model.


Pre-Mining MODFLOW Three-Dimensional Simulation

The two-dimensional model described in the preceding section does not simulate the effect of such vertical components of groundwater movement as those shown in Figure 6. Consequently, I converted the two-dimensional (one layer) model to a three-dimensional model containing eight layers. All of the layers are 200 feet thick except the bottom layer, which is 300 feet thick. Most of the model cells are 200 X 200 X 200 feet. The vertical discretization of the model is shown in Figure 10. The vertical dimension of the model is the same as the one-layer model described above, but vertical movement of the groundwater is simulated as flow between layers. Vertical hydraulic conductivity controls the rate of this vertical movement. Also, I converted the constant head of 9900 feet representing the stream to cell by cell constant heads that represent the stream gradient across the model. The new stream elevation at the west side of the model is 9805 feet, and at the east side of the model it is 9960 feet. I also extended these constant heads southward from the stream to account for the effect of stream connected alluvium. The result is that all of the constant head boundary cells are in layer 4 of the model, and vertical groundwater flow to and from the stream can occur from cells in layer 5 that are below the constant head cells in layer 4. This three-dimensional model inherits no boundary conditions from the AEM model. It is an independent three-dimensional MODFLOW model of the groundwater system related to the mine.

MODFLOW allows one to use different values for hydraulic conductivity in the horizontal and vertical directions. I used successive trial values for these parameters to find ones that produced a realistic pre-mining water table configuration. The result was a vertical hydraulic conductivity value that is much smaller than than the horizontal hydraulic conductivity. The ratio is 0.0011/0.35, meaning that the horizontal hydraulic conductivity is over 300 times greater than the vertical hydraulic conductivity. Such anisotropy can be produced by horizontal layers of low permeability material. The simulated pre-mining water table configuration is shown in Figure 9. The colors of the contours indicate the model layer that contains the water table at that elevation. Blue, green, red, and magenta indicate layers 1, 2, 3, and 4, respectively. The angular contours are an artifact of the finite-difference grid. The blue cells are constant head cells that represent the stream and stream-connected alluvium.

Pre-mining water table simulated with 3D MODFLOW.

Figure 9. Pre-mining water table simulated with the three-dimensional MODFLOW grid.


Figure 10 is a north-south vertical cross-section through the mine location. It illustrates the configuration of the water table and shows the vertical discretization of the model. The brown curve is the land surface. The blue curve is the simulated water table. The blue cell is constant head representing the stream. Inactive cells are black. Model cells containing the water table are partially saturated. Cells above these are dry cells and cells below them are saturated. The mine tunnels are not shown because this is the pre-mining water table. The lower tunnel is in layer 4 (9930 feet) and the upper tunnel is in layer 3 (10100 feet).

Vertical section through the MODFLOW model.

Figure 10. Vertical cross-section through the MODFLOW model at the mine location.


Dewater MODFLOW Three-Dimensional Simulation

I performed successive trial simulations to find a water well extraction scenario that resulted in the water table being beneath the entire length of both tunnels. The result was a well at the end of the lower tunnel extracting 40 gpm, plus a well at the end of the upper tunnel extracting 5 gpm. The resulting hydraulic head in layer 3, which contains the upper tunnel, is shown in Figure 11. The blank area north of the merged contours contains dry cells. These are model cells that are entirely above the water table. The hydraulic head in layer 4, which contains the lower tunnel, is shown in Figure 12. A side effect of this dewatering is extraction of 3 gpm from the stream, which would be returned to the stream as flow from the wells; and interception of 42 gpm, which would also be returned.

Head in MODFLOW layer 3 with 2 dewatering wells

Figure 11. Hydraulic head in MODFLOW model layer 3 with two dewatering wells.


Head in MODFLOW layer 4 with two dewatering wells.

Figure 12. Hydraulic head in MODFLOW model layer 4 with two dewatering wells.


In this dewatering simulation, layers 1 and 2 are entirely dry, indicating the water table is below 10200 feet throughout the model. Whereas, when the wells are not extracting water, the water table is above 10500 feet in the southeast part of the model, as shown in Figure 9. This result indicates that the dewatering shifts the actual groundwater divide beyond the south edge of the model. Simulating this shift would require the model area to be expanded southward to the next stream. Expanding the model is beyond the scope of this MODFLOW modeling project. An effect of the unchanged impermeable boundary is to cause the dewatering extraction rates to be greater than they would be if the model were expanded to the next stream.

Conclusion

MODFLOW's more realistic treatment of the water table aquifer resulted in differences in water table elevation when compared to the AEM model. MODFLOW provides a smaller estimate of the extraction rates that would be required for the dewatering wells. The combined rate for the two wells in the MODFLOW model is 38 percent of the rate in the AEM model (45 gpm versus 120 gpm). Another application of MODFLOW may be seen on the Surficial Sediment MODFLOW model page.

Footnotes

1A finite-difference model divides the 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.

2Hydraulic head (aka head) at a point in a groundwater system is the elevation to which water would rise in a pipe if its bottom were at that point. A value for head is calculated for each active cell in a finite-difference model. In a water table aquifer, head in a cell that contains the water table is the elevation of the water table. Since this web page is non-technical, I have used the term "water table" instead of "head" to simplify the presentation. Technically, where the head is above the top of the cells (10600 feet), the aquifer is saturated and is acting as a confined aquifer, where "water table" is not the correct term for "head."

3A general head boundary allows flow into the model proportional to the head in the active cells at the boundary.

4Hydraulic conductivity is a term generally used in groundwater applications rather than permeability. Hydraulic conductivity is a measure of the ability of a material to conduct the particular fluid being investigated (usually water). Permeability (technically intrinsic permeability) is a property only of the porous medium and is not affected by the properties of the fluid.


Addendum

If you have a question, send an email: ddunn@dunnhydrogeo.com

If you need a consulting hydrogeologist for groundwater modeling, send an email to the above address.

Posted October 31, 2018

Revised May 30, 2019