Flow‐Through Stream Modeling with MODFLOW and MT3D: Certainties and Limitations

This paper aims to assess MODFLOW and MT3D capabilities for simulating the spread of contaminants from a river exhibiting an unusual relationship with an alluvial aquifer, with the groundwater head higher than the river head on one side and lower on the other (flow‐through stream). A series of simulation tests is conducted using a simple hypothetical model so as to characterize and quantify these limitations. Simulation results show that the expected contaminant spread could be achieved with a specific configuration composed of two sets of parameters: (1) modeled object parameters (hydraulic groundwater gradient, hydraulic conductivity values of aquifer and streambed), and (2) modeling parameters (vertical discretization of aquifer, horizontal refinement of stream modeled with River [RIV] package). The influence of these various parameters on simulation results is investigated, and potential complications and errors are identified. Contaminant spread from stream to aquifer is not always reproduced by MT3D due to the RIV package's inability to simulate lateral exchange fluxes between stream and aquifer. This paper identifies the need for a MODFLOW streamflow package allowing lateral stream‐aquifer interactions and streamflow routine calculations. Such developments could be of particular interest for modeling contaminated flow‐through streams.     

Simulating Seepage into Mine Shafts and Tunnels with MODFLOW

In cases when an equivalent porous medium assumption is suitable for simulating groundwater flow in bedrock aquifers, estimation of seepage into underground mine workings (UMWs) can be achieved by specifying MODFLOW drain nodes at the contact between water bearing rock and dewatered mine openings. However, this approach results in significant numerical problems when applied to simulate seepage into an extensive network of UMWs, which often exist at the mine sites. Numerical simulations conducted for individual UMWs, such as a vertical shaft or a horizontal drift, showed that accurate prediction of seepage rates can be achieved by either applying grid spacing that is much finer than the diameter/width of the simulated openings (explicit modeling) or using coarser grid with cell sizes exceeding the characteristic width of shafts or drifts by a factor of 3. Theoretical insight into this phenomenon is presented, based on the so-called well-index theory. It is demonstrated that applying this theory allows to minimize numerical errors associated with MODFLOW simulation of seepage into UMWs on a relatively coarse Cartesian grid. Presented examples include simulated steady-state groundwater flow from homogeneous, heterogeneous, and/or anisotropic rock into a vertical shaft, a horizontal drift/cross-cut, a ramp, two parallel drifts, and a combined system of a vertical shaft connected to a horizontal drift.     

Simulation of Lake-Groundwater Interaction under Steady-State Flow

MODFLOW is one of the most popular groundwater simulation tools available; however, the development of lake modules that can be coupled with MODFLOW is lacking apart from the LAK3 package. This study proposes a new approach for simulating lake - groundwater interaction under steady-state flow, referred to as the sloping lakebed method (SLM). In this new approach, discretization of the lakebed in the vertical direction is independent of the spatial discretization of the aquifer system, which can potentially solve the problem that the lake and groundwater are usually simulated at different scales. The lakebed is generalized by a slant at the bottom of each lake grid cell, which can be classified as fully submerged, dry, and partly submerged. The SLM method accounts for all lake sources and sinks, establishing a governing equation that can be solved using Newton's method. A benchmarking case study was conducted using a modified model setup in the LAK3 user manual. It was found that when there is a sufficient number of layers at the top of the groundwater model, SLM simulates an almost identical groundwater head as the LAK3-based model; when the number of layers decreases, SLM is unaffected while LAK3 may be at a risk of giving unrealistic results. Additionally, the SLM can reflect the relationship between the simulated lake surface area and lake water depth more accurately. Therefore, the SLM method is a promising alternative to the LAK3 package when simulating lake - groundwater interaction.     

Implementation of Solute Transport in the Vadose Zone into the “HYDRUS Package for MODFLOW”

The “HYDRUS package for MODFLOW” is an existing MODFLOW package that allows MODFLOW to simultaneously evaluate transient water flow in both unsaturated and saturated zones. The package is based on incorporating parts of the HYDRUS-1D model (to simulate unsaturated water flow in the vadose zone) into MODFLOW (to simulate saturated groundwater flow). The coupled model is effective in addressing spatially variable saturated-unsaturated hydrological processes at the regional scale. However, one of the major limitations of this coupled model is that it does not have the capability to simulate solute transport along with water flow and therefore, the model cannot be employed for evaluating groundwater contamination. In this work, a modified unsaturated flow and transport package (modified HYDRUS package for MODFLOW and MT3DMS) has been developed and linked to the three-dimensional (3D) groundwater flow model MODFLOW and the 3D groundwater solute transport model MT3DMS. The new package can simulate, in addition to water flow in the vadose zone, also solute transport involving many biogeochemical processes and reactions, including first-order degradation, volatilization, linear or nonlinear sorption, one-site kinetic sorption, two-site sorption, and two-kinetic sites sorption. Due to complex interactions at the groundwater table, certain modifications of the pressure head (compared to the original coupling) and solute concentration profiles were incorporated into the modified HYDRUS package. The performance of the newly developed model is evaluated using HYDRUS (2D/3D), and the results indicate that the new model is effective in simulating the movement of water and contaminants in the saturated-unsaturated flow domains.     

Hydrodynamics of floodplain wetlands in a chalk catchment: The River Lambourn, UK

Variations in floodplain channel water levels and valley floor groundwater levels (measured in piezometers and boreholes) are examined at selected points along the course of the River Lambourn, a chalk river in southern England. A local alluvial gravel aquifer in the valley bottom is associated with numerous small wetlands that extend over much of the river's perennial profile. Variations in hydraulic gradient between local borehole levels and/or floodplain channel water levels are described for three sites in the seasonal section of the channel at Bockhampton, East Garston and West Shefford. The results indicate that observed groundwater levels are closely associated with flows from discrete springs at the margins of the channel and floodplain. However, as the floodplain widens and the alluvial gravel aquifer increases in size, the gravel aquifer accounts for a substantial down-valley component of groundwater flow with a diffuse vertical water flux. In the lower catchment, the exchange of flows between the gravel aquifer and the river enables some attenuation of floodplain water-table variability, providing a stable hydrological regime for valley-bottom wetlands. Catchment controls upon the local, valley-bottom, wetland regime are demonstrated with the application of a simple groundwater model developed using MODFLOW. The model is used to simulate groundwater discharge to the river in the upper and lower catchment, in addition to the water level regime at selected points in the valley bottom in the lower catchment. The results demonstrate the importance of taking catchment-scale water flow into account when managing isolated wetlands in a permeable catchment.     

Embedding isochronous cells overland flow module into MODFLOW

This study integrates a simple overland flow module (isochronous cells model) with the river module of MODFLOW such that temporal and spatial interactions between stream flow and groundwater can be simulated using net rainfall data of a watershed. The isochronous cells model is an efficient travel time runoff approach based on geographic information system (GIS) that considers both spatial and temporal variations of net rainfall through hill slope of the watershed. This overland module is easily coupled with MODFLOW river routing module. Specifically, the stream flow from the isochronous cells model is directly assigned to both sides of river cells of the MODFLOW model. Such an integration of MODFLOW and isochronous cells model is especially useful in watersheds where river flow data are limited. The feasibility of this integrated model was demonstrated using a case study in the middle and downstream regions of the Yitong River watershed, China. Copyright © 2012 John Wiley & Sons, Ltd.     

Vertical Discretization Impact in Numerical Modeling of Matrix Diffusion in Contaminated Groundwater

Understanding the effects of contaminants that can diffuse into low-permeability (“low-k”) zones is crucial for effective groundwater remedial decision-making. Because low-k zones can serve as low-level sources of contamination to more transmissive zones over time, an accurate evaluation of the impacts of matrix diffusion at contaminated sites is vital. This study compared numerical groundwater flow and transport simulations using MODFLOW/RT3D at a hypothetical site using three cases, each with increasing discretization of the vertical 10-m thick domain: (1) a coarse multilayer heterogeneous grid based on one layer for each of four different hydrogeological units, (2) a “low-resolution” discretization approach where the low-k units were divided into several sublayers giving the model 10 layers, and (3) a “high-resolution” numerical model with 199 layers that are a few centimeters thick. When comparing the results of each case, significant differences were observed between the discretizations used, even though all other model input data were identical. The conventional grid models (Cases 1 and 2) appeared to underestimate groundwater plume concentrations by a factor ranging from 1.1 to 36 when compared to the high-resolution grid model (Case 3), and underestimated predicted cleanup times by more than a factor of 10 for some of the hypothetical sampling points in the modeling domain. These results validate the implication of Chapman et al. (2012), that conventional vertical discretization of numerical groundwater flow and transport models at contaminated sites (with layers that are greater than 1 m thick) can lead to significant errors when compared to more accurate high-resolution vertical discretization schemes (layers that are centimeters thick).     

A New Package in MODFLOW to Simulate Unconfined Groundwater Flow in Sloping Aquifers

The nonhorizontal‐model‐layer (NHML) grid system is more accurate than the horizontal‐model‐layer grid system to describe groundwater flow in an unconfined sloping aquifer on the basis of MODFLOW‐2000. However, the finite‐difference scheme of NHML was based on the Dupuit‐Forchheimer assumption that the streamlines were horizontal, which was acceptable for slope less than 0.10. In this study, we presented a new finite‐difference scheme of NHML based on the Boussinesq assumption and developed a new package SLOPE which was incorporated into MODFLOW‐2000 to become the MODFLOW‐SP model. The accuracy of MODFLOW‐SP was tested against solution of Mac Cormack (1969). The differences between the solutions of MODFLOW‐2000 and MODFLOW‐SP were nearly negligible when the slope was less than 0.27, and they were noticeable during the transient flow stage and vanished in steady state when the slope increased above 0.27. We established a model considering the vertical flow using COMSOL Multiphysics to test the robustness of constrains used in MODFLOW‐SP. The results showed that streamlines quickly became parallel with the aquifer base except in the narrow regions near the boundaries when the initial flow was not parallel to the aquifer base. MODFLOW‐SP can be used to predict the hydraulic head of an unconfined aquifer along the profile perpendicular to the aquifer base when the slope was smaller than 0.50. The errors associated with constrains used in MODFLOW‐SP were small but noticeable when the slope increased to 0.75, and became significant for the slope of 1.0.     

Comparing MODFLOW simulation options for predicting intra‐meander flux

During the evolution of meander bends, the intra‐meander groundwater head gradients steepen and generate zones of accelerated water and nutrient intra‐meander fluxes important for ecosystem processes. This paper compares and contrasts three MODFLOW groundwater model packages based on their simulation of intra‐meander flux for two stages of meander evolution observed in a sandbox river table and one level of river bed clogging, where the hydraulic conductivity in the river bed is lower than in the adjacent aquifer. These packages are the Time‐Variant Specified Head package [constant head (CHD)], River package (RIV), and Streamflow‐Routing package (SFR2), each controlling the groundwater or river head bounding the intra‐meander region. The RIV and SFR2 packages fix river stage and allow for variation in groundwater head below the river, which is suggested for simulating intra‐meander flux for all sinuosities with and without river bed clogging whenever river bed parameters are available. The CHD package fixes below river groundwater head and fails to simulate intra‐meander head loss and flux in meanders with high sinuosity or river bed clogging. In low sinuosity meanders and in cases without river bed clogging, there were no significant differences between MODFLOW packages for simulating river intra‐meander head loss and flux. This research demonstrates why MODFLOW users need to consider the limitations of each package when simulating intra‐meander flux in reaches with river bed clogging, high sinuosity, or similarly steep hydraulic gradients. Copyright © 2014 John Wiley & Sons, Ltd.     

Improvement of surface run‐off in the hydrological model ParFlow by a scale‐consistent river parameterization

We propose an improvement of the overland‐flow parameterization in a distributed hydrological model, which uses a constant horizontal grid resolution and employs the kinematic wave approximation for both hillslope and river channel flow. The standard parameterization lacks any channel flow characteristics for rivers, which results in reduced river flow velocities for streams narrower than the horizontal grid resolution. Moreover, the surface areas, through which these wider model rivers may exchange water with the subsurface, are larger than the real river channels potentially leading to unrealistic vertical flows. We propose an approximation of the subscale channel flow by scaling Manning's roughness in the kinematic wave formulation via a relationship between river width and grid cell size, following a simplified version of the Barré de Saint‐Venant equations (Manning–Strickler equations). The too large exchange areas between model rivers and the subsurface are compensated by a grid resolution‐dependent scaling of the infiltration/exfiltration rate across river beds. We test both scaling approaches in the integrated hydrological model ParFlow. An empirical relation is used for estimating the true river width from the mean annual discharge. Our simulations show that the scaling of the roughness coefficient and the hydraulic conductivity effectively corrects overland flow velocities calculated on the coarse grid leading to a better representation of flood waves in the river channels.     

Grid cell distortion and MODFLOW's integrated finite-difference numerical solution

The ground water flow model MODFLOW inherently implements a nongeneralized integrated finite-difference (IFD) numerical scheme. The IFD numerical scheme allows for construction of finite-difference model grids with curvilinear (piecewise linear) rows. The resulting grid comprises model cells in the shape of trapezoids and is distorted in comparison to a traditional MODFLOW finite-difference grid. A version of MODFLOW-88 (herein referred to as MODFLOW IFD) with the code adapted to make the one-dimensional DELR and DELC arrays two dimensional, so that equivalent conductance between distorted grid cells can be calculated, is described. MODFLOW IFD is used to inspect the sensitivity of the numerical head and velocity solutions to the level of distortion in trapezoidal grid cells within a converging radial flow domain. A test problem designed for the analysis implements a grid oriented such that flow is parallel to columns with converging widths. The sensitivity analysis demonstrates MODFLOW IFD's capacity to numerically derive a head solution and resulting intercell volumetric flow when the internal calculation of equivalent conductance accounts for the distortion of the grid cells. The sensitivity of the velocity solution to grid cell distortion indicates criteria for distorted grid design. In the radial flow test problem described, the numerical head solution is not sensitive to grid cell distortion. The accuracy of the velocity solution is sensitive to cell distortion with error <1% if the angle between the nonparallel sides of trapezoidal cells is <12.5 degrees. The error of the velocity solution is related to the degree to which the spatial discretization of a curve is approximated with piecewise linear segments. Curvilinear finite-difference grid construction adds versatility to spatial discretization of the flow domain. MODFLOW-88's inherent IFD numerical scheme and the test problem results imply that more recent versions of MODFLOW 2000, with minor modifications, have the potential to make use of a curvilinear grid.     

Investigating meteorological effect on river flow recession rate in an arid environment

ABSTRACT

This study explores previously unreported linkage between recession rates of rainfall hyetograph and river flow hydrograph in an arid environment in Sultanate of Oman. Ephemeral streams are hydraulically disconnected from groundwater aquifers and depend on rainfall to produce water that flows only for hours or a few days at most. It therefore prompted to hypothesize that the recession rate of the rainfall event controls the corresponding recession of river flow. To test this assumption, 1-h river flow rates and 20-min rainfall rates in Al-Khoud catchment area for the period 1997–2013 were analysed. The river flow recession rate and antecedent river flow were found to be inversely proportional, while their relation improved with increasing time span of cumulating the antecedent river flow. The results further show that the simulation of river flow recession rate can be improved by incorporating the combined effects of rainfall recession rate and antecedent moisture content.     

Numerical modeling groundwater recharge and its implication in water cycles of two interdunal valleys in the Sand Hills of Nebraska

The topography and geomorphology of the sand dunes and interdunal valleys in the Nebraska Sand Hills play important roles in regional water cycle by influencing groundwater recharge and evapotranspiration (ET). In this study, groundwater recharge, associated with precipitation and ET as well as soil hydraulics, and its spatial variations owing to the topography of dunes and valleys are examined. A method is developed to describe the recharge as a function of the storage capacity of dunes of various heights. After the method is tested using observations from a network of wells in the Sand Hills, it is used in the MODFLOW model to simulate and describe recharge effects on groundwater table depth at two different dune-valley sites. Analysis of modeled groundwater budget shows that the groundwater table depth in the interdunal valleys is critically influenced by vertical groundwater flows from surrounding dunes. At the site of higher dunes there are steadier and larger vertical groundwater flows in the dunes from their previous storage of precipitation. These vertical flows change to be horizontal converging groundwater flows and create upwelling in the interdunal valleys, where larger ET loss at the surface further enhances groundwater upwelling. Such interdunal valley is the major concentration area of the surface water and groundwater flow in the Sand Hills. At the site of shallow dunes and a broad interdunal valley the supply of groundwater from the dunes is trivial and inadequate to support upwelling of groundwater in the valley. The groundwater flows downward in the valley, and the valley surface is dry. Weak ET loss at the surface has a smaller effect on the groundwater storage than the precipitation recharge, making such area a source for groundwater.     

A method for estimating spatially variable seepage and hydraulic conductivity in channels with very mild slopes

Infiltration along ephemeral channels plays an important role in groundwater recharge in arid regions. A model is presented for estimating spatial variability of seepage due to streambed heterogeneity along channels based on measurements of streamflow‐front velocities in initially dry channels. The diffusion‐wave approximation to the Saint‐Venant equations, coupled with Philip's equation for infiltration, is connected to the groundwater model MODFLOW and is calibrated by adjusting the saturated hydraulic conductivity of the channel bed. The model is applied to portions of two large water delivery canals, which serve as proxies for natural ephemeral streams. Estimated seepage rates compare well with previously published values. Possible sources of error stem from uncertainty in Manning's roughness coefficients, soil hydraulic properties and channel geometry. Model performance would be most improved through more frequent longitudinal estimates of channel geometry and thalweg elevation, and with measurements of stream stage over time to constrain wave timing and shape. This model is a potentially valuable tool for estimating spatial variability in longitudinal seepage along intermittent and ephemeral channels over a wide range of bed slopes and the influence of seepage rates on groundwater levels. Copyright © 2012 John Wiley & Sons, Ltd.     

Field assessment of surface water–groundwater connectivity in a semi‐arid river basin (Murray–Darling,Australia)

In semi‐arid and arid river basins, understanding the connectivity between rivers and alluvial aquifers is one of the key challenges for the management of groundwater resources. The type of connection present (gaining, losing‐connected, transitional and losing‐disconnected) was assessed at 12 sites along six Murray–Darling Basin river reaches. The assessments were made by measuring the hydraulic head in the riparian zone near the rivers to evaluate if the water tables intersected the riverbeds and by measuring fluid pressure (ψ) in the riverbeds. The rationale for the latter was that ψ will always be greater than or equal to zero under connected conditions (either losing or gaining) and always lesser than or equal to zero under losing‐disconnected conditions. A mixture of losing‐disconnected, losing‐connected and gaining conditions was found among the 12 sites. The losing‐disconnected sites all had a riverbed with a lower hydraulic conductivity than the underlying aquifer, usually in the form of a silty clay or clay unit 0.5–2 m in thickness. The riparian water tables were 6 to 25 m below riverbed level at the losing‐disconnected sites but never lower than 1 m below riverbed level at the losing‐connected ones. The contrast in water table depth between connected and disconnected sites was attributed to the conditions at the time of the study, when a severe regional drought had generated a widespread decline in regional water tables. This decline was apparently compensated near losing‐connected rivers by increased infiltration rates, while the decline could not be compensated at the losing‐disconnected rivers because the infiltration rates were already maximal there. Copyright © 2012 John Wiley & Sons, Ltd.     

Stream water infiltration, bank storage, and storage zone changes due to stream-stage fluctuations

During a flood period, stream-stage increases induce infiltration of stream water into an aquifer; subsequent declines in stream stage cause a reverse motion of the infiltrated water. This paper presents the results of the water exchange rate between a stream and aquifer, the storage volume of the infiltrated stream water in the surrounding aquifer (bank storage), and the storage zone. The storage zone is the part of aquifer where groundwater is replaced by stream water during the flood. MODFLOW was used to simulate stream–aquifer interactions and to quantify rates of stream infiltration and return flow. MODPATH was used to trace the pathlines of the infiltrated stream water and to determine the size of the storage zone. Simulations were focused on the analyses of the effects of the stream-stage fluctuation, aquifer properties, the hydraulic conductivity of streambed sediments, regional hydraulic gradients, and recharge and evapotranspiration (ET) rates on stream–aquifer interactions. Generally, for a given stream–aquifer system, larger flow rates result from larger stream-stage fluctuations; larger storage volumes and storage zones are produced by larger and longer-lasting fluctuations. For a given stream-stage hydrograph, a lower-permeable streambed, an aquitard, or an anisotropic aquifer of low vertical hydraulic conductivity can significantly reduce the rate of infiltration and limit the size of the storage zone. The bank storage solely caused by the stage fluctuation differs slightly between gaining and losing streams. Short-term rainfall recharge and ET loss in the shallow groundwater slightly influence on the flow rate, but their effects on bank storage in a larger area for a longer period can be considerable.     

Estimating streambed parameters for a disconnected river

Evaluation of stream–aquifer interaction and water balance for a catchment often requires specific information on streambed parameters, such as streambed hydraulic conductivity, seepage flux across the streambed and so on. This paper describes a simple, inexpensive instrument that is used to measure these streambed parameters under the condition of a stream disconnected from groundwater. Our method includes a seepage cylinder for simulation of river water depth. The proposed method was applied to estimate the vertical hydraulic conductivity of a streambed and the changes in vertical seepage rate from stream to groundwater with varied stream water depth in the Manasi River of Xinjiang Uygur Autonomous Region, China. The vertical hydraulic conductivities of the streambed determined from 12 sites along the Manasi River vary from 1.01 to 29.m/day where the stream disconnects from the groundwater. The experimental results suggest that there are two kinds of relations between the vertical seepage rate and the simulated stream water depth. One is a linear relation between the two variables with low Reynolds numbers (less than 10); the other is a nonlinear relation (exponential relation) between the two variables with larger Reynolds numbers (greater than 10). This second relationship is quite different from the traditional model that usually calculates the vertical seepage rate from stream to groundwater under the condition of disconnection using a linear relation (Darcy's Law). Our results suggest that a linear relation can only be used for a limited range of river water depth. This method gives a convenient tool for rapidly estimating the streambed hydraulic conductivity and the changes in the vertical seepage rate across streambed with varied stream water depths for the case of a stream disconnected from groundwater. Copyright © 2013 John Wiley & Sons, Ltd.     

Application of an integrated surface water‐groundwater model to multi‐aquifers modeling in Choushui River alluvial fan,Taiwan

This research proposes a combination of SWAT and MODFLOW, MD‐SWAT‐MODFLOW, to address the multi‐aquifers condition in Choushui River alluvial fan, Taiwan. The natural recharge and unidentified pumping/recharge are separately estimated. The model identifies the monthly pumping/recharge rates in multi‐aquifers so that the daily streamflow can be simulated correctly. A multi‐aquifers condition means a subsurface formation composed of at least the unconfined aquifer, the confined aquifer, and an in‐between aquitard. In such a case, the variation of groundwater level is related to pumping/recharge activities in vertically adjacent aquifer and the river‐aquifer interaction. Both factors in turn affect the streamflow performance. Results show that MD‐SWAT‐MODFLOW performs better than SWAT alone in terms of simulated streamflow, especially during low flow period, when pumping/recharge rates are properly estimated. A sensitivity analysis of individual parameter suggests that the vertical leakance may be the most sensitive among all investigated MODFLOW parameters in terms of the estimated pumping/recharge among aquifers, and the Latin‐Hypercube‐One‐factor‐At‐a‐Time sensitivity analysis indicates that the hydraulic conductivity of channel is the most sensitive to the model performance. It also points out the necessity to simultaneously estimate pumping/recharge rates in multi‐aquifers. The estimated net pumping rate can be treated as a lower bound of the actual local pumping rate. As a whole, the model provides the spatio‐temporal groundwater use, which gives the authorities insights to manage groundwater resources. Copyright © 2012 John Wiley & Sons, Ltd.     

Dip and anisotropy effects on flow using a vertically skewed model grid

Darcy flow equations relating vertical and bedding-parallel flow to vertical and bedding-parallel gradient components are derived for a skewed Cartesian grid in a vertical plane, correcting for structural dip given the principal hydraulic conductivities in bedding-parallel and bedding-orthogonal directions. Incorrect-minus-correct flow error results are presented for ranges of structural dip (0 < or = theta < or = 90) and gradient directions (0 < or = phi < or = 360). The equations can be coded into ground water models (e.g., MODFLOW) that can use a skewed Cartesian coordinate system to simulate flow in structural terrain with deformed bedding planes. Models modified with these equations will require input arrays of strike and dip, and a solver that can handle off-diagonal hydraulic conductivity terms.