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.     

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.     

Analysis of the thickness of a fresh water lens and of the transition zone between this lens and upwelling saline water

In regions with saline groundwater, fresh water lenses may develop due to rainwater infiltration. The amount of fresh water that is available for e.g. agricultural crops depends on the thickness of the lens and the extent of mixing between fresh and saline water. In this paper, we consider the mixing of fresh water and upward moving saline ground water in low-lying deltaic areas. The parameters that dominate the flow and transport problem are investigated using dimensionless groups and scaled sensitivities. We characterize the numerically simulated thicknesses of the lens and of the mixing zone by spatial moments. Rayleigh number and mass flux ratio, which is the ratio of the salt water seepage and the precipitation, determine the thickness of the fresh water lens. The local thickness of the mixing zone is mainly influenced by the dispersive/diffusive groups and the mass flux ratio. In addition, convergence of streamlines towards an outflow boundary affects the thickness, particularly in the vicinity of this boundary.Analytical and numerical steady state solutions for lens thickness are compared, taking into account upward seepage, for the two cases with and without a density difference between lens and underlying groundwater. Agreement between the numerical and analytical solutions for the lens thickness is good except when the mass flux ratio becomes small. For zero mass flux ratio, it is implicitly assumed in the analytical solution that salt water is stagnant, and that is unrealistic. Relative contributions of longitudinal and transversal hydrodynamic dispersion and diffusion to the thickness of the mixing zone are quantified numerically for different phases of lens formation. Longitudinal dispersion dominates in the early stages of lens formation, while diffusion and transversal dispersion dominate at steady state.     

Grid resolution study of ground water flow and transport

Three-dimensional grids representing a heterogeneous, ground water system are generated at 10 different resolutions in support of a site-scale flow and transport modeling effort. These grids represent hydrostratigraphy near Yucca Mountain, Nevada, consisting of 18 stratigraphic units with contrasting fluid flow and transport properties. The grid generation method allows the stratigraphy to be modeled by numerical grids of different resolution so that comparison studies can be performed to test for grid quality and determine the resolution required to resolve geologic structure and physical processes such as fluid flow and solute transport. The process of generating numerical grids with appropriate property distributions from geologic conceptual models is automated, thus making the entire process easy to implement with fewer user-induced errors. The series of grids of various resolutions are used to assess the level at which increasing resolution no longer influences the flow and solute transport results. Grid resolution is found to be a critical issue for ground water flow and solute transport. The resolution required in a particular instance is a function of the feature size of the model, the intrinsic properties of materials, the specific physics of the problem, and boundary conditions. The asymptotic nature of results related to flow and transport indicate that for a hydrologic model of the heterogeneous hydrostratigraphy under Yucca Mountain, a horizontal grid spacing of 600 m and vertical grid spacing of 40 m resolve the hydrostratigraphic model with sufficient precision to accurately model the hypothetical flow and solute transport to within 5% of the value that would be obtained with much higher resolution.     

Patterns and Age Distribution of Ground-Water Flow to Streams

Simulations of ground-water flow in a generic aquifer system were made to characterize the topology of ground-water flow in the stream subsystem and to evaluate its relation to deeper ground-water flow. The flow models are patterned after hydraulic characteristics of aquifers of the Atlantic Coastal Plain and are based on numerical solutions to three-dimensional, steady-state, unconfined flow. The models were used to evaluate the effects of aquifer horizontal-to-vertical hydraulic conductivity ratios, aquifer thickness, and areal recharge rates on flow in the stream subsystem. A particle tracker was used to determine flow paths in a stream subsystem, to establish the relation between ground-water seepage to points along a simulated stream and its source area of flow, and to determine ground-water residence time in stream subsystems. In a geometrically simple aquifer system with accretion, the source area of flow to streams resembles an elongated ellipse that tapers in the down gradient direction. Increased recharge causes an expansion of the stream subsystem. The source area of flow to the stream expands predominantly toward the stream headwaters. Base flow gain is also increased along the reach of the stream. A thin aquifer restricts ground-water flow and causes the source area of flow to expand near stream headwaters and also shifts the start-of-flow to the drainage basin divide. Increased aquifer anisotropy causes a lateral expansion of the source area of flow to streams. Ground-water seepage to the stream channel originates both from near- and far-recharge locations. The range in the lengths of flow paths that terminate at a point on a stream increase in the downstream direction. Consequently, the age distribution of ground water that seeps into the stream is skewed progressively older with distance downstream. Base flow ia an integration of ground water with varying age and potentially different water quality, depending on the source within the drainage basin. The quantitative results presented indicate that this integration can have a wide and complex residence time range and source distribution.     

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.     

Simulation of submarine ground water discharge to a marine estuary: Biscayne Bay, Florida

Variable density ground water flow models are rarely used to estimate submarine ground water discharge because of limitations in computer speed, data availability, and availability of a simulation tool that can minimize numerical dispersion. This paper presents an application of the SEAWAT code, which is a combined version of MODFLOW and MT3D, to estimate rates of submarine ground water discharge to a coastal marine estuary. Discharge rates were estimated for Biscayne Bay, Florida, for the period from January 1989 to September 1998 using a three-dimensional, variable density ground water flow and transport model. Hydrologic stresses in the 10-layer model include recharge, evapotranspiration, ground water withdrawals from municipal wellfields, interactions with surface water (canals in urban areas and wetlands in the Everglades), boundary fluxes, and submarine ground water discharge to Biscayne Bay. The model was calibrated by matching ground water levels in monitoring wells, baseflow to canals, and the position of the 1995 salt water intrusion line. Results suggest that fresh submarine ground water discharge to Biscayne Bay may have exceeded surface water discharge during the 1989, 1990, and 1991 dry seasons, but the average discharge for the entire simulation period was only approximately 10% of the surface water discharge to the bay. Results from the model also suggest that tidal canals intercept fresh ground water that might otherwise have discharged directly to Biscayne Bay. This application demonstrates that regional scale variable density models are potentially useful tools for estimating rates of submarine ground water discharge.     

Interpreting velocities from heat-based flow sensors by numerical simulation

We have carried out numerical simulations of three-dimensional nonisothermal flow around an in situ heat-based flow sensor to investigate how formation heterogeneities can affect the interpretation of ground water flow velocities from this instrument. The flow sensor operates by constant heating of a 0.75-m-long, 5-cm-diameter cylindrical probe, which contains 30 thermistors in contact with the formation. The temperature evolution at each thermistor can be inverted to obtain an estimate of the ground water flow velocity vector using the standard interpretive method, which assumes that the formation is homogeneous. Analysis of data from heat-based flow sensors installed in a sand aquifer at the Former Fort Ord Army Base near Monterey, California, suggested an unexpected component of downward flow. The magnitudes of the vertical velocities were expected to be much less than those of the horizontal velocities at this site because the sensors were installed just above a clay aquitard. Numerical simulations were conducted to examine how differences in thermal conductivities may lead to spurious indications of vertical flow velocities. We found that a decrease in the thermal conductivity near the bottom of the sensor can perturb the temperature profiles along the instrument in such a manner that analyses assuming homogeneous thermal conductivity could indicate a vertical flow component even though flow is actually horizontal. This work demonstrates how modeling can be used to simulate instrument response to formation heterogeneity and shows that caution must be used in interpreting data from such devices.     

Flow deformation of liquefied sand under constant shear load and its application to analysis of flow slide of infinite slope

This paper intended to evaluate the behavior of saturated sand and sloped ground subjected to flow failure with seepage of pore water in the ground after earthquake and the resultant liquefaction. Triaxial compression tests of sand with constant deviator stress but changing of pore pressure and volume of the specimens were conducted in this study. It was revealed that the relation between the volume change and the amount of shear strain during deformation depended on the initial density of the sand but it did not much depend on shear stress and initial confining stress levels. Based on this test results and numerical analysis of the seepage of pore water in liquefied ground, a methodology was proposed to predict the deformation of inclined ground due to liquefaction.     

A finite element model for simulating surface run‐off and unsaturated seepage flow in the shallow subsurface

This work introduces water–air two‐phase flow into integrated surface–subsurface flow by simulating rainfall infiltration and run‐off production on a soil slope with the finite element method. The numerical model is formulated by partial differential equations for hydrostatic shallow flow and water–air two‐phase flow in the shallow subsurface. Finite element computing formats and solution strategies are presented to obtain a numerical solution for the coupled model. An unsaturated seepage flow process is first simulated by water–air two‐phase flow under the atmospheric pressure boundary condition to obtain the rainfall infiltration rate. Then, the rainfall infiltration rate is used as an input parameter to solve the surface run‐off equations and determine the value of the surface run‐off depth. In the next iteration, the pressure boundary condition of unsaturated seepage flow is adjusted by the surface run‐off depth. The coupling process is achieved by updating the rainfall infiltration rate and surface run‐off depth sequentially until the convergence criteria are reached in a time step. A well‐conducted surface run‐off experiment and traditional surface–subsurface model are used to validate the new model. Comparisons with the traditional surface–subsurface model show that the initiation time of surface run‐off calculated by the proposed model is earlier and that the water depth is larger, thus providing values that are closer to the experimental results.     

Evaluating ground water-sea water interactions via resistivity and seepage meters

We investigated submarine ground water discharge and salt water-fresh water interactions at two locations along the shoreline of the Upper Gulf of Thailand to evaluate mechanisms of water and material transport into the coastal zone. Our data set illustrates the value of using a combined approach consisting of automatic seepage meters to monitor flow rates while assessing the conductivity (salinity) of the subterranean fluids via remote resistivity measurements. Negative correlations between electric conductivities of fluids measured directly inside seepage meter chambers and the remotely assessed resistivities of subsurface pore water show that such measurements may evaluate the spatial distribution of flow rates as well as the subterranean water quality in the coastal zone. Combined seepage and resistivity measurements may thus provide a more complete understanding of coastal ground water dynamics.     

Modeling axisymmetric flow and transport

Unmodified versions of common computer programs such as MODFLOW, MT3DMS, and SEAWAT that use Cartesian geometry can accurately simulate axially symmetric ground water flow and solute transport. Axisymmetric flow and transport are simulated by adjusting several input parameters to account for the increase in flow area with radial distance from the injection or extraction well. Logarithmic weighting of interblock transmissivity, a standard option in MODFLOW, can be used for axisymmetric models to represent the linear change in hydraulic conductance within a single finite-difference cell. Results from three test problems (ground water extraction, an aquifer push-pull test, and upconing of saline water into an extraction well) show good agreement with analytical solutions or with results from other numerical models designed specifically to simulate the axisymmetric geometry. Axisymmetric models are not commonly used but can offer an efficient alternative to full three-dimensional models, provided the assumption of axial symmetry can be justified. For the upconing problem, the axisymmetric model was more than 1000 times faster than an equivalent three-dimensional model. Computational gains with the axisymmetric models may be useful for quickly determining appropriate levels of grid resolution for three-dimensional models and for estimating aquifer parameters from field tests.     

Insights into bedload transport processes of a large regulated gravel‐bed river

A comprehensive monitoring programme focusing on bedload transport behaviour was conducted at a large gravel‐bed river. Innovative monitoring strategies were developed during five years of preconstruction observations accompanying a restoration project. A bedload basket sampler was used to perform 55 cross‐sectional measurements, which cover the entire water discharge spectrum from a 200‐year flood event in 2013 to a rare low flow event. The monitoring activities provide essential knowledge regarding bedload transport processes in large rivers. We have identified the initiation of motion under low flow conditions and a decrease in the rate of bedload discharge with increasing water discharge around bankfull conditions. Bedload flux strongly increases again during high flood events when the entire inundation area is flooded. No bedload hysteresis was observed. The effective discharge for bedload transport was determined to be near mean flow conditions, which is therefore at a lower flow discharge than expected. A numerical sediment transport model was able to reproduce the measured sediment transport patterns. The unique dataset enables the characterisation of bedload transport patterns in a large and regulated gravel‐bed river, evaluation of modern river engineering measures on the Danube, and, as a pilot project has recently been under construction, is able to address ongoing river bed incision, unsatisfactory ecological conditions for the adjacent national park and insufficient water depths for inland navigation. Copyright © 2017 John Wiley & Sons, Ltd.     

Maintenance of headland-associated linear sandbanks: modelling the secondary flows and sediment transport

Linear sandbanks are located globally in areas where there are strong currents and an abundance of sand. In the recent years, these sandbanks have become of strategic interest as a potential source of marine aggregates (sand and gravel) and mineral deposits. They form the seaward boundary of the nearshore zone and therefore are important for the stability of the coastal system. They also commonly reach the sea surface and thus pose a threat to navigation. Headland-associated linear sandbanks are a specific type of sandbanks which are located in the lee of coastal topographic features such as headlands and islands. Interaction between tidal currents and topographic features generate complex three-dimensional circulation patterns that significantly influence the distribution of sediments in the vicinity of the feature. Field and numerical model investigations of the three-dimensional flow structure have been undertaken on the Levillain Shoal, a headland-associated linear sandbank present in the lee of Cape Levillain (Shark Bay, Western Australia). The field data indicated the presence of secondary flows near the tip of the cape and around the bank which were re-produced in the numerical simulations. Numerical results have shown that residual eddies are not representative of the sediment transport and that secondary currents enhance the convergence of sediment towards the sandbank. Maintenance processes have been investigated. Sediment transport paths near the cape and the bank indicate that the sandbank is part of a sand circulation cell where the sand is circulating around the bank with exchanges between the sandbank and the headland.     

Numerical implementation of a backward probabilistic model of ground water contamination

Backward location and travel time probabilities can be used to characterize known and unknown sources or prior positions of ground water contamination. Backward location probability describes the position of the observed contamination at some time in the past; backward travel time probability describes the amount of time prior to observation that the contamination was released from its source or was at a particular upgradient location. The governing equation for backward probabilities is the adjoint of the governing equation for contaminant transport, but with new load terms. Numerical codes that have been written to solve the forward equations of contaminant transport, e.g., the advection-dispersion equation, can also be used to solve the adjoint equation for location and travel time probabilities; however, the interpretation of the results is different and some new approximations must be made for the load terms. We present the governing equations for backward location and travel time probabilities, and provide appropriate numerical approximations for these load terms using the cell-centered finite difference method, one of the most popular numerical methods in ground water hydrology. We discuss some additional numerical considerations for the backward model including boundary conditions, reversal of the flow field, and interpretation of the results. We illustrate the implementation of the backward probability model using hypothetical examples in one- and two-dimensional domains. We also present a three-dimensional application of a pump-and-treat remediation capture zone delineation at the Massachusetts Military Reservation. The illustrations are performed using MODFLOW-96 for flow simulations and MT3DMS for transport simulations.     

Using MODFLOW 2000 to model ET and recharge for shallow ground water problems

In environments with shallow ground water elevation, small changes in the water table can cause significant variations in recharge and evapotranspiration fluxes. Particularly, where ground water is close to the soil surface, both recharge and evapotranspiration are regulated by a thin unsaturated zone and, for accuracy, must be represented using nonconstant and often nonlinear relationships. The most commonly used ground water flow model today, MODFLOW, was originally designed with a modular structure with independent packages representing recharge and evaporation processes. Systems with shallow ground water, however, may be better represented using either a recharge function that varies with ground water depth or a continuous recharge and evapotranspiration function that is dependent on depth to water table. In situations where the boundaries between recharging and nonrecharging cells change with time, such as near a seepage zone, a continuous ground water flux relationship allows recharge rates to change with depth rather than having to calculate them at each stress period. This research article describes the modification of the MODFLOW 2000 recharge and segmented evapotranspiration packages into a continuous recharge-discharge function that allows ground water flux to be represented as a continuous process, dependent on head. The modifications were then used to model long-term recharge and evapotranspiration processes on a saline, semiarid floodplain in order to understand spatial patterns of salinization, and an overview of this process is given.     

Seepage face height, water table position, and well efficiency at steady state

When a fully penetrating well pumps an ideal unconfined aquifer at steady state, the water table usually does not join the water level in the well. There is a seepage face inside the well, which is a key element in evaluating the well performance. This problem is analyzed using the finite-element method, solving the complete equations for saturated and unsaturated flow. The seepage face position is found to be almost independent of the unsaturated zone properties. The numerical results are used to test the validity of several analytic approximations. Equations are proposed to predict the seepage face position at the pumping well for any well drawdown, and the water table position at any distance from the pumping well for any in-well drawdown. Practical hints are provided for installing monitoring wells and evaluating well efficiency.     

Thermal Imagery of Groundwater Seeps: Possibilities and Limitations

Quantifying groundwater flow at seepage faces is crucial because seepage faces influence the hydroecology and water budgets of watersheds, lakes, rivers and oceans, and because measuring groundwater fluxes directly in aquifers is extremely difficult. Seepage faces provide a direct and measurable groundwater flux but there is no existing method to quantitatively image groundwater processes at this boundary. Our objective is to determine the possibilities and limitations of thermal imagery in quantifying groundwater discharge from discrete seeps. We developed a conceptual model of temperature below discrete seeps, observed 20 seeps spectacularly exposed in three dimensions at an unused limestone quarry and conducted field experiments to examine the role of diurnal changes and rock face heterogeneity on thermal imagery. The conceptual model suggests that convective air‐water heat exchange driven by temperature differences is the dominant heat transfer mechanism. Thermal imagery is effective at locating and characterizing the flux of groundwater seeps. Areas of active groundwater flow and ice growth can be identified from thermal images in the winter, and seepage rates can be differentiated in the summer. However, the application of thermal imagery is limited by diverse factors including technical issues of image acquisition, diurnal changes in radiation and temperature, and rock face heterogeneity. Groundwater discharge rates could not be directly quantified from thermal imagery using our observations but our conceptual model and experiments suggest that thermal imagery could quantify groundwater discharge when there are large temperature differences, simple cliff faces, non‐freezing conditions, and no solar radiation.     

A Spatially Periodic Solute Boundary for MT3DMS and PHT3D

The assumption of spatial repetition is commonly made when producing bedform scale models of the hyporheic zone. Two popular solute transport codes, MT3DMS and PHT3D, do not currently provide the necessary boundary condition required to simulate spatial periodicity in hyporheic zone transport problems. In this study, we develop a spatially periodic boundary (SPB) for solutes that is compatible with a SPB that was previously developed for MODFLOW to simulate the flow component of spatially periodic problems. The approach is ideal for simulating groundwater flow and transport patterns under repeating surface features, such as ripples or dunes on the bottom of a lake or stream. The appropriate block‐centered finite‐difference approach to implement the boundary is presented and the necessary source code modifications are discussed. The performance of the solute SPB, operating in conjunction with the groundwater flow SPB, is explored through comparison of a multi‐bedform hyporheic‐zone model with a single bedform variant. The new boundary conditions perform well in situations where both dispersive effects and lateral seepage flux in the underflow regime beneath the hyporheic zone are minimal.