Ground-water flow was simulated for the Missouri River alluvial aquifer using the three-dimensional finite-difference ground-water flow model MODFLOWARC (Orzol and McGrath, 1992). MODFLOWARC is a modified version of MODFLOW (McDonald and Harbaugh, 1988) that reads and writes files used by ARC/INFO a geographic information system (GIS).
Three-dimensional simulation of ground-water flow in the alluvial aquifer was necessary to accurately determine the hydraulic-head distribution beneath the main rivers and around the multiple well fields in the study area. Discharge from the aquifer to rivers may vary according to river size or depth of incision. Ground-water flow may be divided into smaller flow subsystems because of the degree of interaction of between ground-water and the larger and smaller rivers in the study area. Three-dimensional simulation was also necessary in the analysis of ground-water travel times and the contributing recharge area around each pumping well field because of the vertical flow of ground-water caused by well pumping. Also, pumping from the well fields located near the Missouri River can induce recharge from the river and cause ground-water flow beneath the river.
The modeled area is 73 kilometers by 24 kilometers in size and contains the entire study area. (Click for here for location map [15.4 kb]) The model uses a uniform grid size of 150 meters by 150 meters and contains 310,400 cells in 4 layers, 485 columns, and 160 rows. The irregular shape of the study area reduced the number of active cells in the model to 67,362 with 20,835 active cells in layer 1, 22,198 active cells in layer 2, 17,978 active cells in layer 3, and 6,351 active cells in layer 4. (Click for here for map of all layers [15.0 kb]) The regular grid spacing facilitated data input from the GIS and analysis of model output by the GIS, and the grid size minimized errors in flowpath analysis caused by large grid size (Pollock, 1994; Zheng, 1994).
The model represents the alluvial aquifer using four layers, numbered 1 to 4, of variable thickness with no intervening confining layers. Layer 1 corresponds to the upper part of the aquifer where clays, silts, and finer-grained sands are dominant. The thickness of layer 1 is large enough to account for the anticipated range of water-level variation within the aquifer during ground-water flow simulation and was modeled using unconfined-aquifer hydraulic properties during transient conditions. Layers 2 and 3 correspond to the middle part of the aquifer where sand and gravelly sands predominate. These layers were not anticipated to dewater during the simulations and were modelled using confined aquifer hydraulic properties during transient conditions. Layer 4 corresponds to the deep parts of the aquifer where gravels and sandy gravels are present and was also modeled using confined-aquifer hydraulic properties during transient conditions. Hydrogeologic data from over 1400 locations within the study area was used in the model. Values for hydrogeologic parameters were assigned to each model cell by interpolation from the GIS.
Lithologic descriptions recorded during the installation of wells and boreholes are the most numerous and have the greatest areal extent of all data types in the GIS. The distribution of clay, silt, sand, and gravel within the aquifer was used to distribute the hydraulic conductivity and storage coefficients among model cells and the depth to bedrock was used to determine the geometry of the model.
The Missouri, Kansas, Blue, Little Blue and Fishing Rivers and Cooley Lake are each represented in the model as a head-dependent flow boundary. For each cell in the model where a river or lake affects ground-water flow the altitude of the river or lake stage must be known. Flow into or out of the aquifer at each of the cells where a river or lake is simulated is a function of the river or lake stage with respect to the altitude of the potentiometric surface, the hydraulic conductivity of the riverbed or lakebed material, the cross-sectional area of flow between the river or lake bed and the aquifer, and the altitude of the potentiometric surface with respect to the altitude of the riverbed or lakebed (McDonald and Harbaugh, 1988, p 6-5).
The water table, the upper boundary of the alluvial aquifer, was simulated in the model as a free-surface boundary and was the boundary across which areally distributed recharge entered the aquifer. Recharge to the model was applied to the top-most active cell in each vertical column and was varied areally and temporally as a function of precipitation and the average vertical hydraulic conductivity of class of soils in each model cell (Kelly and Blevins, 1995). Recharge has been estimated to be between 2 and 25 percent of precipitation which was recorded daily at the Kansas City Municipal Airport, between 8-1-93 and 10-31-93 and in Independence, Missouri between 11-1-93 and 3-27-94. Both locations are near the middle of the study area. The adjustment factors for the rate of recharge for each soil vertical hydraulic conductivity class is shown in table 3. Urbanized parts of the study area with unknown vertical hydraulic conductivity of soils were assigned average recharge rates for the study area.
The alluvial valley walls and bedrock were simulated in the model as no-flow boundaries, a form of the specified flow boundary. As discussed previously, the rate of water flow between the alluvial aquifer and the valley walls and bedrock has not been quantified. However, the hydraulic conductivities of the alluvial valley walls and bedrock are several orders of magnitude less than hydraulic conductivities in the alluvial aquifer. Therefore, simulating the alluvial valley walls and bedrock as no-flow boundaries is reasonable because the amount of flow is a negligible percentage of the total flow.
Pumping wells are internal boundaries of the model where water was removed from the model at a specified rate equal to the discharge of each well. The total volume of water withdrawn from the aquifer by pumping of public water-supply wells was recorded monthly for most well fields and were obtained from each water supplier when available. Pumping rates for industrial wells were obtained from owners records. The total annual volume of water removed from the aquifer by other wells in the study area including industrial wells, irrigation wells and small public water suppliers were obtained from the Missouri Department of Natural Resources (1991). Pumping rates for wells with no records were estimated based on type of use and pump rating. The depth of each pumping well in the model was based on the altitude of the screened interval when known or the bottom-most layer when the altitude of the screened interval was unknown.
Hydraulic conductivity or transmissivity data are available for 94 locations within the study area. However, locations for which lithologic data are known are more numerous and have the widest distribution within the study area. Aquifer tests conducted during previous investigations to determine hydraulic conductivity or transmissivity typically were performed in wells where the lithology and altitude of the screened interval were known. Other reported values of hydraulic conductivity were performed in a laboratory on actual samples of aquifer material. Hydraulic data collected during previous investigations were entered into the GIS (Kelly and Blevins, 1995) and used to associate a hydraulic conductivity value with a specific lithology. Typical ranges of the hydraulic conductivities of clays, silts, sands, and gravels (Freeze and Cherry, 1979; Driscoll, 1986) were used where hydraulic conductivity data were unavailable for a specific lithology. Once the relation between hydraulic conductivity and lithologic type was established, hydraulic conductivity was distributed among cells in each layer of the model based on the distribution of lithology among cells in each layer of the model.
The simulated flow of water between model cells in adjacent model layers is controlled by the vertical conductance term. The vertical conductance terms between cells of adjacent layers were calculated using the method presented in McDonald and Harbaugh, 1988, p 5-11 and then multiplied by a factor to simulate the presence of vertical anisotropy in clay, silt, and fine sand deposits. The vertical anisotropy was assumed to decrease with depth because of the increase in particle grain size with depth and the higher probability that fine grained layered depositional features such as overbank and channel fill deposits have been reworked or removed by erosional and depositional processes of the Missouri River. Therefore the vertical conductance terms between adjacent cells of layers 1 and 2 were multiplied by a factor of 0.1 and between layers 2 and 3 by a factor of 0.5. The vertical conductance between layers 3 and 4 was not reduced.
The specific yield is the unconfined storage coefficient for layer 1. Typical specific yield values (Driscoll, 1986) were distributed among model cells in layer 1 based on the distribution of lithology Specific Yield - Layer 1 (11.9 kb). A storage coefficient of 0.001 was used for layers 2, 3, and 4 and represents conditions where water is released from storage due to expansion of the water or compaction of the aquifer material and not actual drainage of the aquifer.
The ground-water flow model was calibrated by adjusting model input data and model geometry to modify model output so that model results matched field observations within an acceptable level of accuracy (Konikow, 1978). Parameters changed during the calibration process include hydraulic conductivity, vertical conductance between model layers, specific yield, river conductance, and recharge rates. After each change in one of these parameters, the simulation was run and simulated ground-water levels were compared to observed ground-water levels. The model accuracy was calculated using the root mean square (RMS) error between actual synoptic measurements of hydraulic head and model generated hydraulic head at the end of each model run. Model accuracy is increased by minimizing the RMS error. The RMS error measures the absolute value of the variation between measured and simulated hydraulic heads.
The strategy for calibration of the ground-water flow model was to use both quasi-steady state hydraulic head data and transient hydraulic head data. The quasi-steady state calibration was used to assess model geometry, confirm the conceptual model of ground-water flow, test the appropriateness of simulated boundary conditions, and obtain approximate transmissivity and recharge arrays in preparation for more rigorous transient calibration. The transient calibration was used to fine-tune the model hydraulic properties through a period of prolonged aquifer drainage from August 1993, just after the peak of the flood of 1993, to February 1994, when river stage and ground-water levels had approached typical conditions for that time of year.
The quasi-steady state hydraulic head data were obtained from the January 1993 synoptic water-level measurement (22.0 kb) conducted during a previous study by Kelly and Blevins (1995). The January 1993 data represent the closest approximation of steady state conditions where water levels, river stage, antecedent precipitation (used to calculate recharge rates), and well pumping data were readily available. However, true steady state conditions were not present hence the name quasi-steady state.
The rationale for using a quasi-steady state calibration was based on the complexity and size of the model and the availability of synoptic water-level data for the study area for January 1993 when hydrologic conditions were relatively stable. Initial calibration simulations used uniform hydraulic conductivity values for each layer, uniform recharge rates, and no well pumping. Errors in model geometry were corrected during this stage of calibration in several locations where data was limited. Large differences between measured and modeled hydraulic head in the area near Lake City and Buckner, Missouri indicated that the buried alluvial channels in the eastern part of the study area needed to be included in the simulation because they are an important outlet for ground-water flow. Other changes to the model geometry in areas where the depth to bedrock was unknown included the extension of layer 4 to the upstream and downstream model boundaries of the Missouri River alluvial aquifer and the lateral extension of layer 3 to more closely approximate the bowl shape of the alluvial aquifer. Numerical instabilities in initial simulations were caused by a large number of model cells in layer 1 going dry. To correct this problem, the thicknesses of layers 2 and 3 were decreased and the thickness of layer 1 was increased to allow more vertical movement of the simulated water table during the iterative solution.
The areal distribution of hydraulic conductivity was based on the lithologic distribution within each layer. Recharge rates were correlated to precipitation and spatially distributed according to soil vertical hydraulic conductivity data (Kelly and Blevins, 1995). Average pumping rates for wells in the study area were determined using pumping records from water suppliers and industries or from the Missouri Department of Natural Resources (1991) when unavailable. The layer, row, column, and withdrawal rate for each well used in the quasi-steady-state calibration simulation are listed for each well field with the results of the particle tracking analysis later in the report. These input parameters were added to the model during the quasi-steady-state calibration simulations. The assignment of hydraulic conductivity based on lithologic distribution, recharge rate based on the vertical permeability of soils, and inclusion of well pumping reduced the RMS error of the quasi-steady-state calibration simulation and more realistically simulated the January 1993 distribution of hydraulic head. However, the RMS error for these simulations was still larger than the maximum measurement error. The effect of the smaller streams and constructed drainage ditches on ground-water flow were not included in the original conceptual model of the flow system. Model simulations indicated that these surface-water bodies had a measurable effect on ground-water flow and the hydraulic head distribution. The addition of these head-dependent flow boundaries to simulate the effects of small streams and ditches on ground-water flow reduced the RMS error to an acceptable value of 1.15 m
The difference between flow into the model and flow out of the model across all model boundaries was 0.13 percent of total flow for the quasi-steady state calibration. The level of accuracy of the simulation in representing the January 1993 hydraulic head distribution was accepted because the hydraulic head distribution for that time was not completely at steady state. Further calibration of the quasi-steady-state model would have resulted in erroneously changing model input parameters to match a hydraulic head distribution that partially resulted from transient ground-water flow.
Transient calibration of the ground-water flow model to hydrologic conditions measured between August 1, 1993 and March 27, 1994 matched the change over time of the simulated hydraulic head distribution with the change over time of the measured hydraulic head distribution. This was done by measuring the changes in various hydrologic stresses that affected the distribution of hydraulic head and simulating those stresses in the model. A stress on the ground-water flow system was any change in river stage, recharge, or well pumping that caused the ground-water flow system, and in particular, the distribution of hydraulic head to change. These changes actually occurred as gradual increases or decreases of river stage, intermittent and varying rates of recharge from precipitation, and intermittent or constant pumping of wells at varying rates. The ground-water-flow model applied areal and temporal changes in stress to the ground-water-flow system by using a series of stress periods. Within each stress period river stage, recharge, and pumping rates of wells were held constant. Each stress period in the transient calibration was 7 days in length. Each stress period was divided into three time steps used by the model to calculate the change of hydraulic head within the stress period. In each stress period the first time step was 1.47 days, the second time step was 2.21 days and the third time step was 3.32 days. Therefore, a series of stress periods each corresponding to a specific time and each representing average stress to the ground-water flow system for that specific time, were used to simulate the effects of changes in stress on ground-water flow. The length of the stress periods was chosen based on the time intervals for data collection and the approximately one week period required to collect the synoptic water level measurements during October 1993 and February 1994.
River stage averaged over each stress period was assigned to each river model cell using interpolation methods previously discussed. Average river stage altitudes for each gaging station in the study area are listed for each stress period in Appendix A in the back of the report. Average recharge for each stress period was assigned to the top most active cell in each vertical column using average precipitation and the vertical hydraulic conductivity class of soil at the cell (Kelly and Blevins, 1995). Precipitation amounts from August 1, 1993 to March 27, 1994 are shown in fig 30. The percentage of precipitation that was supplied to the model as recharge is shown for each soil vertical hydraulic conductivity class (Kelly and Blevins, 1995) in table 3. Average well pumping for each stress period was assigned to each model cell that contained a pumping well. The layer, row, column, and withdrawal rate for each well used in the transient calibration simulation are listed for each stress period in Appendix A in the back of the report.
The hydraulic head data used for transient calibration of the ground-water flow model were obtained from August 1993 flood data and two synoptic water-level measurements from 123 wells between October 18 and 22, 1993 and from 98 wells between February 14 and 18, 1994. The transient calibration simulation was begun right after the flood in the Missouri River basin. At that time the entire alluvial aquifer in the study area was at or near full saturation from flood waters and local precipitation. The initial hydraulic head for each active cell in the model was set equal to the corresponding land-surface altitude. The model was allowed to run 81 simulated days divided into 5 stress periods using the August, 1 1993 river-stage and well-pumping input parameters to produce an initial hydraulic-head distribution at the beginning of the transient simulation that resulted from flood-stage altitudes of rivers in the study area.
Model results from stress periods that corresponded to the October 1993 and February 1994 synoptic water level measurements were compared and the RMS error was calculated for the corresponding stress periods of each transient calibration simulation. For the accepted calibration simulation, the October 1993 RMS error was 0.71 meters and the February RMS error was 0.80 meters. These values are less than the maximum measurement errors previously discussed and indicate the acceptability of the calibrated model.
A sensitivity analysis was performed to assess the response of the model simulation to changes in various input parameter values. The model is sensitive to a parameter when a change of the parameter value changes the distribution of simulated hydraulic head. When the model is sensitive to an input parameter, the value and distribution of that parameter within the model are more accurately determined during model calibration because small changes to the parameter value cause large changes in hydraulic head. If a change of parameter value does not change the simulated hydraulic head distribution, the model is insensitive to that parameter. When the model is insensitive to an input parameter, the value and distribution of that parameter within the model are more difficult to accurately determine from model calibration because large changes to the parameter do not cause large changes in hydraulic head. These values of these parameters may not represent actual values.
Values of the hydraulic conductivity, the vertical conductance between layers, the specific yield of layer 1, the specific storage of layers 2, 3 and 4, the conductance of riverbed material, and the recharge flux were decreased by 50 percent below, and increased by 50 percent above the calibrated parameter value during the sensitivity analysis. The effect of these changes in the value of input parameters on model output were assessed with the transient calibration scenario. In the following discussion, increases or decreases in the RMS error calculated for each sensitivity analysis are with respect to the calibrated RMS errors for October 1993 (0.71 m) and February 1994 (0.80 m).
The sensitivity of the model to changes in hydraulic conductivity is greatest when the hydraulic conductivity is decreased in all layers at the same time. Decreasing the hydraulic conductivity by 50 percent in all layers increased the October 1993 RMS error by 66 percent and the February 1994 RMS error by 18.8 percent. The sum of the individual percent changes in RMS errors caused by decreasing the hydraulic conductivity by 50 percent in each layer is 40.3 percent for October 1993 and 11.9 percent for February 1994. Increasing the hydraulic conductivity by 50 percent in all layers decreased the October 1993 RMS error by 11.8 percent and increased the February 1994 RMS error by 18.4 percent. The sum of the percent changes in RMS errors caused by increasing the hydraulic conductivity by 50 percent in each layer resulted in a 21.5 percent decrease for the October 1993 RMS error and a 17.6 percent increase for the February 1994 RMS error. The model is more sensitive to the sum of individual increases of hydraulic conductivity for the October 1993 data than the February 1994 data.
The model is sensitive to overall changes in the hydraulic conductivity as indicated by the large change of RMS error when the hydraulic conductivity is simultaneously changed in all layers. Also, the effect of changing the hydraulic conductivity for all layers and in each individual layer caused a larger change in the October 1993 RMS error than in the February 1994 RMS error indicating that the model is more sensitive to hydraulic conductivity earlier in the calibration simulation than later. Also, decreasing the hydraulic conductivity by 50 percent in all layers and in each layer caused the October 1993 RMS error to increase and the February 1994 RMS error to decrease. Increasing the hydraulic conductivity by 50 percent in all layers and in each layer caused the October 1993 RMS error to decrease and the February 1994 RMS error to increase.
The insensitivity of the model to changes in the vertical conductance between layers is evident when the small percent changes of the October 1993 and February 1994 RMS errors are compared to the corresponding percent changes in vertical conductance during the sensitivity analysis. The only change in RMS error (-0.7 percent) occurred in the October 1993 RMS error when the vertical conductance was decreased by 50 percent between layers 1 and 2 and when the vertical conductance was decreased by 50 percent between layers 3 and 4.
The sensitivity of the model to changes in the specific yield in layer 1 is evident when the relatively large changes in the October 1993 and February 1994 RMS errors are compared to the corresponding change of specific yield. The small change in October 1993 and February 1994 RMS errors indicate the model is insensitive to changes in the confined storage coefficient and riverbed conductance. The model is somewhat sensitive to changes in recharge rate but large changes in recharge rate result in minor changes in the RMS error.
Back to Missouri River Alluvial Aquifer Ground-Water Protection Project Homepage
Contact address:
Brian P. Kelly, Hydrologist
U.S. Geological Survey
401 NW Capital Drive
Lee's Summit, MO 64086
PHONE: 816-554-2414
FAX: 816-554-9273
EMAIL: bkelly@usgs.gov
U.S. Department of the Interior, U.S. Geological Survey
Maintainer: Rita Choate, Webmaster
Updated: September 4, 2003
Privacy Statement || Disclaimer
|| FOIA || Accessibility
URL: http://missouri.usgs.gov/indep/kelly/mo-alluvial-gw/model/index.htm
![(Click for here for location map [15.4 kb])](../images/locplot.gif){kind=link}
![(Click for here for map of all layers [15.0 kb])](images/modlays.gif){kind=link}
