An assembly model for simulation of large-scale ground water flow and transport

When managing large-scale ground water contamination problems, it is often necessary to model flow and transport using finely discretized domains--for instance (1) to simulate flow and transport near a contamination source area or in the area where a remediation technology is being implemented; (2) to account for small-scale heterogeneities; (3) to represent ground water-surface water interactions; or (4) some combination of these scenarios. A model with a large domain and fine-grid resolution will need extensive computing resources. In this work, a domain decomposition-based assembly model implemented in a parallel computing environment is developed, which will allow efficient simulation of large-scale ground water flow and transport problems using domain-wide grid refinement. The method employs common ground water flow (MODFLOW) and transport (RT3D) simulators, enabling the solution of almost all commonly encountered ground water flow and transport problems. The basic approach partitions a large model domain into any number of subdomains. Parallel processors are used to solve the model equations within each subdomain. Schwarz iteration is applied to match the flow solution at the subdomain boundaries. For the transport model, an extended numerical array is implemented to permit the exchange of dispersive and advective flux information across subdomain boundaries. The model is verified using a conventional single-domain model. Model simulations demonstrate that the proposed model operated in a parallel computing environment can result in considerable savings in computer run times (between 50% and 80%) compared with conventional modeling approaches and may be used to simulate grid discretizations that were formerly intractable.     

MODFLOW/MT3DMS-based simulation of variable-density ground water flow and transport

This paper presents an approach for coupling MODFLOW and MT3DMS for the simulation of variable-density ground water flow. MODFLOW routines were modified to solve a variable-density form of the ground water flow equation in which the density terms are calculated using an equation of state and the simulated MT3DMS solute concentrations. Changes to the MODFLOW and MT3DMS input files were kept to a minimum, and thus existing data files and data files created with most pre- and postprocessors can be used directly with the SEAWAT code. The approach was tested by simulating the Henry problem and two of the saltpool laboratory experiments (low- and high-density cases). For the Henry problem, the simulated results compared well with the steady-state semianalytic solution and also the transient isochlor movement as simulated by a finite-element model. For the saltpool problem, the simulated breakthrough curves compared better with the laboratory measurements for the low-density case than for the high-density case but showed good agreement with the measured salinity isosurfaces for both cases. Results from the test cases presented here indicate that the MODFLOW/MT3DMS approach provides accurate solutions for problems involving variable-density ground water flow and solute transport.     

Comparison of an algebraic multigrid algorithm to two iterative solvers used for modeling ground water flow and transport

Numerical solution of large-scale ground water flow and transport problems is often constrained by the convergence behavior of the iterative solvers used to solve the resulting systems of equations. We demonstrate the ability of an algebraic multigrid algorithm (AMG) to efficiently solve the large, sparse systems of equations that result from computational models of ground water flow and transport in large and complex domains. Unlike geometric multigrid methods, this algorithm is applicable to problems in complex flow geometries, such as those encountered in pore-scale modeling of two-phase flow and transport. We integrated AMG into MODFLOW 2000 to compare two- and three-dimensional flow simulations using AMG to simulations using PCG2, a preconditioned conjugate gradient solver that uses the modified incomplete Cholesky preconditioner and is included with MODFLOW 2000. CPU times required for convergence with AMG were up to 140 times faster than those for PCG2. The cost of this increased speed was up to a nine-fold increase in required random access memory (RAM) for the three-dimensional problems and up to a four-fold increase in required RAM for the two-dimensional problems. We also compared two-dimensional numerical simulations of steady-state transport using AMG and the generalized minimum residual method with an incomplete LU-decomposition preconditioner. For these transport simulations, AMG yielded increased speeds of up to 17 times with only a 20% increase in required RAM. The ability of AMG to solve flow and transport problems in large, complex flow systems and its ready availability make it an ideal solver for use in both field-scale and pore-scale modeling.     

International borders, ground water flow, and hydroschizophrenia

A substantial body of research has been conducted on transboundary water, transboundary water law, and the mitigation of transboundary water conflict. However, most of this work has focused primarily on surface water supplies. While it is well understood that aquifers cross international boundaries and that the base flow of international river systems is often derived in part from ground water, transboundary ground water and surface water systems are usually managed under different regimes, resulting in what has been described as "hydroschizophrenia." Adding to the problem, the hydrologic relationships between surface and ground water supplies are only known at a reconnaissance level in even the most studied international basins, and thus even basic questions regarding the territorial sovereignty of ground water resources often remain unaddressed or even unasked. Despite the tensions inherent in the international setting, riparian nations have shown tremendous creativity in approaching regional development, often through preventive diplomacy, and the creation of "baskets of benefits," which allow for positive-sum, integrative allocations of joint gains. In contrast to the notion of imminent water wars, the history of hydropolitical relations worldwide has been overwhelmingly cooperative. Limited ground water management in the international arena, coupled with the fact that few states or countries regulate the use of ground water, begs the question: will international borders serve as boundaries for increased "flows" of hydrologic information and communication to maintain strategic aquifers, or will increased competition for shared ground water resources lead to the potential loss of strategic aquifers and "no flows" for both ground water users?     

Fracture control of ground water flow and water chemistry in a rock aquitard

There are few studies on the hydrogeology of sedimentary rock aquitards although they are important controls in regional ground water flow systems. We formulate and test a three-dimensional (3D) conceptual model of ground water flow and hydrochemistry in a fractured sedimentary rock aquitard to show that flow dynamics within the aquitard are more complex than previously believed. Similar conceptual models, based on regional observations and recently emerging principles of mechanical stratigraphy in heterogeneous sedimentary rocks, have previously been applied only to aquifers, but we show that they are potentially applicable to aquitards. The major elements of this conceptual model, which is based on detailed information from two sites in the Maquoketa Formation in southeastern Wisconsin, include orders of magnitude contrast between hydraulic diffusivity (K/S(s)) of fractured zones and relatively intact aquitard rock matrix, laterally extensive bedding-plane fracture zones extending over distances of over 10 km, very low vertical hydraulic conductivity of thick shale-rich intervals of the aquitard, and a vertical hydraulic head profile controlled by a lateral boundary at the aquitard subcrop, where numerous surface water bodies dominate the shallow aquifer system. Results from a 3D numerical flow model based on this conceptual model are consistent with field observations, which did not fit the typical conceptual model of strictly vertical flow through an aquitard. The 3D flow through an aquitard has implications for predicting ground water flow and for planning and protecting water supplies.     

Influence of calibration methodology on ground water flow predictions

We constructed a numerical model of transient ground water flow and solute transport for a portion of the Biscayne Aquifer in Florida, and calibrated the model with three different combinations of data from a 193-day period: head (h) data alone, data on h and ground water discharge to a canal (q), and data on h, q, and ground water chloride concentration (C). We used each of the three calibrated models to predict h and q during a 182-day test period separate from the calibration period. All three calibrated models predicted h equally well during the test period (r = 0.95, where r = 1 indicates perfect agreement between measured and simulated values), though the model calibrated on h alone had significantly different parameter values than the other two models. Predictions of q during the test period depended on calibration methodology; models calibrated with multiple targets simulated q more accurately than the model calibrated on h alone (r = 0.79 compared to r = 0.49). Based on the results of these simulations, we conclude: (1) Post-calibration prediction is important in assessing the value of different data types in automated calibration; (2) inverse-solution uniqueness is not a requirement for accurate h predictions; (3) relatively simple models can predict with reasonable accuracy transient ground water flow in a complex aquifer, and parameters governing this prediction can be estimated by nonlinear regression methods that incorporate both h and q data; (4) addition of C data to the calibration did not improve model predictive capacity because the information in the C data was similar to that in the q data, from the perspective of model calibration (the subsurface chemical signal in question was controlled mainly by seepage of high-chloride canal water into the low-chloride ground water system).     

Heat flow and ground water flow in the Great Plains of the United States

Regional groundwater flow in deep aquifers adds advective components to the surface heat flow over extensive areas within the Great Plains province. The regional groundwater flow is driven by topographically controlled piezometric surfaces for confined aquifers that recharge either at high elevations on the western edge of the province or from subcrop contacts. The aquifers discharge at lower elevations to the east. The assymetrical geometry for the Denver and Kennedy Basins is such that the surface areas of aquifer recharge are small compared to the areas of discharge. Consequently, positive advective heat flow occurs over most of the province. The advective component of heat flow in the Denver Basin is on the order of 15 mW m⁻² along a zone about 50 km wide that parallels the structure contours of the Dakota aquifer on the eastern margin of the Basin. The advective component of heat flow in the Kennedy Basin is on the order of 20 mW m⁻² and occurs over an extensive area that coincides with the discharge areas of the Madison (Mississippian) and Dakota (Cretaceous) aquifers. Groundwater flow in Paleozoic and Mesozoic aquifers in the Williston Basin causes thermal anomalies that are seen in geothermal gradient data and in oil well temperature data. The pervasive nature of advective heat flow components in the Great Plains tends to mask the heat flow structure of the crust, and only heat flow data from holes drilled into the crystalline basement can be used for tectonic heat flow studies.     

Implications of ground water chemistry and flow patterns for earthquake studies

Ground water can facilitate earthquake development and respond physically and chemically to tectonism. Thus, an understanding of ground water circulation in seismically active regions is important for earthquake prediction. To investigate the roles of ground water in the development and prediction of earthquakes, geological and hydrogeological monitoring was conducted in a seismogenic area in the Yanhuai Basin, China. This study used isotopic and hydrogeochemical methods to characterize ground water samples from six hot springs and two cold springs. The hydrochemical data and associated geological and geophysical data were used to identify possible relations between ground water circulation and seismically active structural features. The data for delta18O, deltaD, tritium, and 14C indicate ground water from hot springs is of meteoric origin with subsurface residence times of 50 to 30,320 years. The reservoir temperature and circulation depths of the hot ground water are 57 degrees C to 160 degrees C and 1600 to 5000 m, respectively, as estimated by quartz and chalcedony geothermometers and the geothermal gradient. Various possible origins of noble gases dissolved in the ground water also were evaluated, indicating mantle and deep crust sources consistent with tectonically active segments. A hard intercalated stratum, where small to moderate earthquakes frequently originate, is present between a deep (10 to 20 km), high-electrical conductivity layer and the zone of active ground water circulation. The ground water anomalies are closely related to the structural peculiarity of each monitoring point. These results could have implications for ground water and seismic studies in other seismogenic areas.     

The role of hand calculations in ground water flow modeling

Most ground water modeling courses focus on the use of computer models and pay little or no attention to traditional analytic solutions to ground water flow problems. This shift in education seems logical. Why waste time to learn about the method of images, or why study analytic solutions to one-dimensional or radial flow problems? Computer models solve much more realistic problems and offer sophisticated graphical output, such as contour plots of potentiometric levels and ground water path lines. However, analytic solutions to elementary ground water flow problems do have something to offer over computer models: insight. For instance, an analytic one-dimensional or radial flow solution, in terms of a mathematical expression, may reveal which parameters affect the success of calibrating a computer model and what to expect when changing parameter values. Similarly, solutions for periodic forcing of one-dimensional or radial flow systems have resulted in a simple decision criterion to assess whether or not transient flow modeling is needed. Basic water balance calculations may offer a useful check on computer-generated capture zones for wellhead protection or aquifer remediation. An easily calculated "characteristic leakage length" provides critical insight into surface water and ground water interactions and flow in multi-aquifer systems. The list goes on. Familiarity with elementary analytic solutions and the capability of performing some simple hand calculations can promote appropriate (computer) modeling techniques, avoids unnecessary complexity, improves reliability, and is likely to save time and money. Training in basic hand calculations should be an important part of the curriculum of ground water modeling courses.     

Effect of snow structure on water flow and solute transport

Water flow through a melting snow pack modifies its structure and stability and affects the release of water and nutrients into soils and surface waters. Field and laboratory observations indicate a large spatial variability on various scales of the liquid water content and flow, a dominant system feature currently not included in numerical models. We investigated experimentally water and dye tracer movement through microstructurally different snow pack horizons and the persistence of preferential flow paths. Naturally rounded snow of varying grain size was artificially packed to obtain well known conditions by sieving it into rectangular bins. Surface melt was induced with infrared lamps. The flow paths were visualized with tracers and liquid water content was monitored with time domain reflectometry probes. Vertical cuts through the snow pack were imaged. The dye tracer patterns allowed the two flow regimes ‘matrix flow’ and ‘preferential flow’ to be distinguished. Matrix flow is apparently dominated by film and capillary flow in the unsaturated snow matrix. The capillary barrier effect at a boundary between a fine over a coarse textured layer on matrix flow in snow was confirmed. In contrast, preferential flow appears as well‐defined flow fingers that advance from 0·1 to 1 cm s^?1. During a melt phase, the advancing flow fingers enlarge and are only partially time invariant. It remains to be shown whether the continuum concept, including the Darcy–Buckingham law is apt to describe the extremely non‐linear nature of water flow and the travel time of solutes in snow under conditions of melt water percolation. Probably, snow packs that include faceted crystals and large variations in bulk density, feature more pronounced capillary barriers and preferential flow triggering, but also stronger impeding of fingers by lateral dispersion. Further, triggering and persistence of preferential flow is complicated by the usually transient infiltration rate. Copyright © 2004 John Wiley & Sons, Ltd.     

Selenium transport and transformation modelling in soil columns and ground water contamination prediction

Selenium transport and transformation were simulated in a soil column. A one‐dimensional dynamic mathematical and computer model is formulated to simulate, selenate, selenite, selenomethionine, organic selenium, and gaseous selenium. This computer model is based on the mass balance equation, including convective transport, dispersive transport, surface adsorption, oxidation and reduction, volatilization, chemical and biological transformation. The mathematical solution is obtained by the finite difference implicit method. The model was verified by comparison of model results with experimental measurements and also using mass balance calculations in each time step of calculation. For example after 4 days of simulation, the simulated value of adsorbed selenate for depth of 20 cm is 0·2 µmol kg^?1 and the measured value is 0·25 µmol kg^?1. Therefore simulated results are in good agreement with measured values. With this study and its results the distribution of various forms of selenium in soil column to ground water table can be predicted. Copyright © 2008 John Wiley & Sons, Ltd.     

Using hydraulic head measurements in variable-density ground water flow analyses

The use of hydraulic head measurements in ground water of variable density is considerably more complicated than for the case of constant-density ground water. A theoretical framework for dealing with these complications does exist in the current literature but suffers from a lack of awareness among many hydrogeologists. When corrections for density variations are ignored or not properly taken into account, misinterpretation of both ground water flow direction and magnitude may result. This paper summarizes the existing theoretical framework and provides practical guidelines for the interpretation of head measurements in variable-density ground water systems. It will be argued that, provided that the proper corrections are taken into account, fresh water heads can be used to analyze both horizontal and vertical flow components. To avoid potential confusion, it is recommended that the use of the so-called environmental water head, which was initially introduced to facilitate the analysis of vertical ground water flow, be abandoned in favor of properly computed fresh water head analyses. The presented methodology provides a framework for determining quantitatively when variable-density effects on ground water flow need to be taken into account or can be justifiably neglected. Therefore, we recommend that it should become part of all hydrogeologic analyses in which density effects are suspected to play a role.     

An assessment of aquifer storage recovery using ground water flow models

Owing to increased demands on ground water accompanied by increased drawdowns, technologies that use recharge options, such as aquifer storage recovery (ASR), are being used to optimize available water resources and reduce adverse effects of pumping. In this paper, three representative ground water flow models were created to assess the impact of hydrogeologic and operational parameters/factors on recovery efficiency of ASR systems. Flow/particle tracking and solute transport models were used to track the movement of water during injection, storage, and recovery. Results from particle tracking models consistently produced higher recovery efficiency than the solute transport models for the parameters/properties examined because the particle tracking models neglected mixing of the injected and ambient water. Mixing between injected and ambient water affected recovery efficiency. Results from this study demonstrate the interactions between hydrogeologic and operational parameters on predictions of recovery efficiency. These interactions are best simulated using coupled numerical ground water flow and transport models that include the effects of mixing of injected water and ambient ground water.     

Calculating ground water transit time of horizontal flow through leaky aquifers

The calculation of ground water transit times is one important factor in ground water protection. In this paper, we present an analytical solution for the transit time for a Dupuit-type flow system applicable to saturated flow through a horizontal leaky aquifer discharging to a downgradient fixed-head boundary under steady-state conditions. We investigate the influence of leakage when comparing the resulting travel times of our model based on head-dependent leakage with the commonly used model with no leakage and a simplified model with constant leakage. The results show significant differences in the position of the water divide and transit time, suggesting that leakage cannot be ignored.     

Estimating ground water recharge using flow models of perched karstic aquifers

The fraction of rain that is annually recharged to ground water is a function of the transient quantities of precipitation (wet vs. dry years) as well as other meteorological and geologic factors, and thus it is very difficult to estimate. In this study, we have used long records (20 to 30 years) of precipitation and spring discharge to reconstruct the transient character of yearly recharge. These data sets were used to calibrate numerical ground water flow models on the less than 3 km(2) scale for four separate perched karstic aquifers in the Judean and Samarian Mountains of Israel. The stratification and karstic character of the local carbonate rock aquifers cause ground water to flow through discrete dissolution channels and to discharge at isolated springs. An innovative, dual-porosity approach was used where a finite-difference solution simulates flow in the rock matrix, while the karstic channels are simulated using computationally simple drains. Perched conditions are also simulated innovatively using MODFLOW by treating the bottom unsaturated layer as if it is saturated, but by assuming zero pressure head throughout the "unsaturated" layer. Best fitting between measured and computed spring hydrograph data has allowed us to develop a set of empirical functions relating measured precipitation to recharge to the aquifer. The generic methodology presented gives insight into the suspected changes in aquifer recharge rates between particularly wet or dry years.     

Heat flow maps and deep ground water circulation: Examples from Switzerland

The effect of regional and local ground water circulation systems on the Heat Flow Density (HFD) field is demonstrated by two examples from Switzerland, one near St. Gall in an area at the northern border of the Alps, and the other northwest of Zurich along the eastern end of the Jura mountains. Detailed HFD maps of both areas slow pronounced high heat flow zones which are attributed to discharge of subsurface water which has migrated laterally over several 10 km. Seepage velocities on the order of several mm/yr have been calculated. Geothermal information is not available about the infiltration zones where low HFD values are expected. Geochemical and isotopic analysis of water samples from springs and drillholes indicates the recharge zones and demonstrates the effect of extensive regional systems. These results indicate that in regions with significant topographic relief HFD mapping can be seriously biased if drillholes are positioned in valleys which correspond to discharge areas with relatively high HFD, whereas the low heat flow zones remain undetected.     

Effects of model sensitivity and nonlinearity on nonlinear regression of ground water flow

Nonlinear regression is increasingly applied to the calibration of hydrologic models through the use of perturbation methods to compute the Jacobian or sensitivity matrix required by the Gauss-Newton optimization method. Sensitivities obtained by perturbation methods can be less accurate than those obtained by direct differentiation, however, and concern has arisen that the optimal parameter values and the associated parameter covariance matrix computed by perturbation could also be less accurate. Sensitivities computed by both perturbation and direct differentiation were applied in nonlinear regression calibration of seven ground water flow models. The two methods gave virtually identical optimum parameter values and covariances for the three models that were relatively linear and two of the models that were relatively nonlinear, but gave widely differing results for two other nonlinear models. The perturbation method performed better than direct differentiation in some regressions with the nonlinear models, apparently because approximate sensitivities computed for an interval yielded better search directions than did more accurately computed sensitivities for a point. The method selected to avoid overshooting minima on the error surface when updating parameter values with the Gauss-Newton procedure appears for nonlinear models to be more important than the method of sensitivity calculation in controlling regression convergence.     

Use of Picard and Newton iteration for solving nonlinear ground water flow equations

This study examines the use of Picard and Newton iteration to solve the nonlinear, saturated ground water flow equation. Here, a simple three-node problem is used to demonstrate the convergence difficulties that can arise when solving the nonlinear, saturated ground water flow equation in both homogeneous and heterogeneous systems with and without nonlinear boundary conditions. For these cases, the characteristic types of convergence patterns are examined. Viewing these convergence patterns as orbits of an attractor in a dynamical system provides further insight. It is shown that the nonlinearity that arises from nonlinear head-dependent boundary conditions can cause more convergence difficulties than the nonlinearity that arises from flow in an unconfined aquifer. Furthermore, the effects of damping on both convergence and convergence rate are investigated. It is shown that no single strategy is effective for all problems and how understanding pitfalls and merits of several methods can be helpful in overcoming convergence difficulties. Results show that Picard iterations can be a simple and effective method for the solution of nonlinear, saturated ground water flow problems.     

Michigan basin regional ground water flow discharge to three Great Lakes

Ground water discharge to the Great Lakes around the Lower Peninsula of Michigan is primarily from recharge in riparian basins and proximal upland areas that are especially important to the northern half of the Lake Michigan shoreline. A steady-state finite-difference model was developed to simulate ground water flow in four regional aquifers in Michigan's Lower Peninsula: the Glaciofluvial, Saginaw, Parma-Bayport, and Marshall aquifers interlayered with the Till/"red beds," Saginaw, and Michigan confining units, respectively. The model domain was laterally bound by a continuous specified-head boundary, formed from lakes Michigan, Huron, St. Clair, and Erie, with the St. Clair and Detroit River connecting channels. The model was developed to quantify regional ground water flow in the aquifer systems using independently determined recharge estimates. According to the flow model, local stream stages and discharges account for 95% of the overall model water budget; only 50% enters the lakes directly from the ground water system. Direct ground water discharge to the Great Lakes' shorelines was calculated at 36 m3/sec, accounting for 5% of the overall model water budget. Lowland areas contribute far less ground water discharge to the Great Lakes than upland areas. The model indicates that Saginaw Bay receives only approximately 1.13 m3/sec ground water; the southern half of the Lake Michigan shoreline receives only approximately 2.83 m3/sec. In contrast, the northern half of the Lake Michigan shoreline receives more than 17 m3/sec from upland areas.