A Hybrid FE-FV Discontinuous Method to Solve Flow and Transport Equations in Heterogeneous Porous Media
In heterogeneous porous media material discontinuities affect single phase flow as well as multiphase flow. The finite-element method allows to represent such boundaries using piecewise constant or linear material property variations from finite element to element. However, when a node-centered complementary finite volume method is used to model transport across these material interfaces, this discretization can not represent jump discontinuities in concentration or saturation. To overcome this dilemma there are two options, one can either enrich the nodes at the interfaces by additional degrees of freedom or explode the model along the interfaces. This technique, is another way of resolving this problem. Interface nodes are multiplicated so that they match in number the material domains which they join. Here, we use this latter approach, developing an IMPES transport model that can evolve discontinuities at material interfaces in a heterogeneous porous medium. To achieve pressure continuity across the exploded material interfaces, we implement a new implicit algorithm. For the transport equation, we develop a higher-order-accurate scheme which captures jump discontinuity of the transport variable. We use operator-splitting to solve the transport equation after the pressure equation, with the finite volume method, modeling diffusion implicitly employing the finite element method. We verify the new method by a comparison between its results and those of the continuous one. The main advantage of the discontinuous scheme is the resolution of the effects of the material discontinuities. The first order discretization dependency is removed.
A Comparative Analysis of Series And Finite Element Solutions For Flow In Multi-layer Aquifers With Contiguous Layers
Recently, Wong S. & J.R. Craig [Computational Methods in Water Resources XVII International Conference, 2008] have developed a semi-analytic series (SAS) solution method for simulating 2D steady-state groundwater flow in multi-layer aquifers with natural unconformity. The advantages of the approach include the capability to predict groundwater flow in aquifer systems with a geometrically complex structure, i.e., layers that are contiguous but vary (sometimes dramatically) in thickness across the modelled domain. The SAS approach is unorthodox compared to numerical scheme such as finite element (FE) method when handling an aquifer basin with contiguous layers. This research attempts to compare the robustness of the SAS solution with that of the FE solution under variety of different geometric conditions, including the increasing of system aspect ratios and the inclusion of pinching layers. It was found that the SAS approach is a useful benchmarking tool and that for contiguous layers, both FE methods continue to be highly accurate at even large aspect ratios. Based on the benchmarking experiments, some rules of thumb for mesh generation in FE models of regional aquifer systems with both contiguous and discontiguous layering are identified.
Semi-analytical solutions for flow in heterogeneous media represented using pilot points
The analytic element method (AEM) is a semi-analytical approach for solving a wide variety of linear groundwater problems. In all existing analytic element models, hydraulic conductivity is represented using the zonation approach, i.e., the conductivity field is piecewise constant. Here, a new analytic element method is presented whereby the conductivity field may be represented with a smoothly varying surface interpolated between pilot points. This new approach is based upon a special transformation of the dependent variable in terms of the square root of a term representing the deviation from some average or representative hydraulic conductivity within the domain. The method still has the same benefits as the traditional AEM: solutions are exact solutions to the governing differential equation, and boundary conditions are met with extremely high precision. A number of example applications are presented to demonstrate the efficacy of the method.
The Noordbergum Effect and its Impact on the Estimation of Hydraulic Parameters
For almost 80 years 'reverse water level fluctuations' (known as the Noordbergum Effect) have been witnessed in aquitards adjacent to pumped aquifers. At early time during pumping, water in the aquifer is released from compressive storage and is accompanied by a volume reduction of the aquifer both in the horizontal and vertical directions. Deformation in the pumped aquifer causes deformation in the adjacent layers resulting in an increase in pore pressure and corresponding increase in hydraulic head. The use of type curves to interpret drawdown data when the Noordbergum Effect influences the data is evaluated by generating a synthetic pumping test in a confined aquifer using BIOT2 (Hsieh, 1994), which is a 2D axisymmetric groundwater flow model that considers the effects of linear poroelasticity. The drawdown data from this synthetic pumping test were analyzed using type curve analysis to determine aquifer hydraulic conductivity (K) and specific storage (Ss), as well as aquitard hydraulic conductivity (K'). Results showed that the applicability of traditional techniques varied with radial distance from the pumping well. That is, K estimates were closest to the true value near the pumping well, while K' showed the same at large radial distances. Estimates of Ss were slightly greater than the true value but showed little variation with distance from the pumping well. The results from this analysis were used to develop several guidelines for interpreting pumping test data when the Noordbergum Effect appears in drawdown data. The shear modulus (G) of the aquifer and aquitards are important for determining whether or not the Noordbergum Effect will be present in drawdown data. To estimate G from drawdown data exhibiting the Noordbergum effect, the parameter estimation program (PEST) (Doherty, 2004) was coupled with BIOT2 to analyze drawdown data from a synthetic pumping test. The K, K' and Ss values determined through type curve analysis were used as initial estimates in BIOT2/PEST. Estimation of hydraulic parameters and the shear modulus using BIOT2/PEST revealed that accurate values could be obtained for the synthetic case with known parameter values. To assess the applicability of BIOT2/PEST to estimate hydraulic parameters and shear modulus in a field setting, drawdown data exhibiting the Noordbergum Effect collected at the North Campus Research Site (NCRS) at University of Waterloo were analyzed. Parameter estimation was conducted by: a) matching the drawdown data one at a time from each observation point; and b) matching drawdown data from all monitoring ports simultaneously along a single well. Results show that the estimated parameters (K, K', and G) are variable from one location to the next illustrating the highly heterogeneous nature of the glaciofluvial aquifer.
The Effects of Disconnect Entrapped Air on Hydraulic Conductivity in the Presence of Water Table Fluctuations
The hydraulic conductivity of a natural groundwater system can possess high spatial and temporal variability in the presence of an entrapped air phase (quasi-saturated soils) which is a key factor in controlling hydraulic behavior (Faybishenko, 1995). Research conducted by Faybishenko (1995), Frye et al. (1997) and Zlotnik et al. (2007) provide evidence of reduced hydraulic conductivity caused by entrapped air which can be introduced by water table elevation changes. The hypothesis that the decreases in hydraulic conductivity caused by entrapped air are sensitive to fluctuations in the water table was tested using laboratory experiments. This work investigates the degree of air entrapment by immiscible displacement and the effects of applying increasing confining pressures (water table height) on the quasi-saturated hydraulic conductivity for a range of sands. Results to date indicate that the changes in the volume of entrapped gas caused by 2.5m changes in water table height induced changes in the quasi-saturated hydraulic conductivity by approximately a factor of 2. Such results have implications for site interpretation and simulation given that the behavior of quasi-saturated systems has not been incorporated into current general models of flow and transport. Faybishenko, B.A. 1995. Hydraulic behaviour of quasi-saturated soils in the presence of entrapped air: Laboratory experiments. Water Resources Research 31 (10): 2421-2435. Frye, V.A., Selker, J.S. and S.M. Gorelick. 1997. Experimental investigations for trapping oxygen gas in saturated porous media for in situ bioremediation. Water Resources Research 33 (12): 2687-2696. Zlotnik, V.A., Eisenhauer, D.E., Schlautman, D.J., Zurbuchen, B.R. and D. Van Peursem. 2007. Entrapped air effects on dipole flow test in sand tank experiments: Hydraulic conductivity and head distribution. Journal of Hydrology 339: 193-205.
Evaluating Traditional Hydrogeologic Characterization Approaches in a Highly Heterogeneous Glaciofluvial Aquifer
Hydraulic conductivity (K) and specific storage (Ss) estimates are two of the most essential parameters when designing transient groundwater flow models which are commonly used in contaminant transport and water resource investigations. The purpose of this study was to evaluate the effectiveness of traditional hydrogeologic characterization approaches in a highly heterogeneous glaciofluvial aquifer at the North Campus Research Site (NCRS) situated on the University of Waterloo campus. The site is instrumented with four Continuous Multichannel Tubing (CMT) wells containing a total of 28 monitoring points and a multi-screen well used for pumping at different elevations. Continuous soil cores to a depth of approximately 18 m were collected during the installation of the CMTs and the multi-screen well. The cores were subsequently characterized using the Unified Soil Classification System and grain size analysis. Samples were obtained from the core at approximately 10 cm increments and a falling head permeameter was used to make 471 K estimates. The estimates from the falling head permeameter showed K to vary from 10-4 - 10-10 m/s illustrating the highly heterogeneous nature of the aquifer at the NCRS. A geostatistical analysis performed on the core K dataset yielded a strongly heterogeneous K field for the site. K and Ss estimates were also obtained via slug tests in the CMT ports through type curve analysis. Cross-hole pumping tests were conducted using the center multi-screened well and the 4 CMTs installed in a 5-spot pattern. Pumping was conducted in 7 zones using a straddle packer system and the corresponding drawdown responses were recorded in 28 zones in the CMTs and 3 zones in the center well using pressure transducers. The various K and Ss estimates were then evaluated by simulating the transient drawdown data using a 3D forward numerical model constructed using Hydrogeosphere (Therrien et al., 2005). Simulation was conducted using 3 separate K and Ss fields: 1) a homogeneous case with K and Ss estimates obtained by averaging equivalent K and Ss values from the cross-hole pumping tests, 2) a layered heterogeneous case with K and Ss estimates from the slug tests and 3) a heterogeneous case with the kriged K from the permeameter tests and equivalent Ss from the slug tests. Results showed that, while drawdown predictions, in general, improved as more complexity was introduced into the model, the ability to make accurate drawdown predictions at all of the CMT ports was inconsistent. These results suggest that new techniques may be required to accurately capture subsurface heterogeneity for improved predictions of drawdown responses. Accurate predictions of drawdown responses should lead to improved modeling of contaminant transport and remediation designs.
Effects of Model Layer Simplification Using Composite Hydraulic Properties
The effects of simplified hydraulic property layering within an unconfined aquifer that overlies a leaky confining unit was assessed using radial, axisymmetric two-dimensional (2D) models of aquifer tests and three- dimensional (3D) models. The 3D models simulated time-varying recharge to the unconfined aquifer and pumping stresses from the deeper confined principal water supply aquifer. Hydraulic properties for distinct lithologic units within the unconfined aquifer and leaky confining unit were derived for two sites from the analyses of aquifer-test data. Conceptual flow models were developed by gradually reducing the number of lithologic units in the unconfined aquifer and leaky confining unit by calculating composite hydraulic properties for the simplified lithologic units using either the thickness-weighted average for the hydraulic properties, or from inverse modeling using regression-based parameter estimation techniques. Simulated flow accuracy was assessed using the 2D and 3D models for the simplified lithologic units. The largest residuals in the unconfined aquifer heads, confined hydrogeologic unit heads, and leakage rates to the confined aquifer were simulated when the unconfined aquifer and leaky confining unit were aggregated into a single layer (quasi-3D). The ground-water flow models with one lithologic unit each for the unconfined aquifer and leaky confining unit (fully-3D) resulted in lower head and leakage rate residuals than the quasi-3D model.
Hydrogeology of Regional Valley Fill Aquifers with Mountain System Recharge
Groundwater in the North Okanagan was investigated using an integrated physical, geochemical and numerical approach. The North Okanagan Groundwater Characterization and Assessment (NOGWCA) project began with an investigation of the geology and hydrostratigraphy of the North Okanagan region. The Deep Creek and Fortune Creek watersheds were found to contain multiple valley-fill aquifers which are recharged via mountain system recharge (MSR) and direct recharge to unconfined aquifers in the valley bottom. Detailed hydrometric data indicates groundwater recharge within the alluvial fan of Fortune Creek, and discharge to surface water in the lower reaches of Deep Creek. Valley side recharge from the adjacent mountains generates artesian conditions in the valley center. Physical hydrogeological measurements and groundwater and surface water geochemistry were used to determine the overall groundwater flow regime, inter-aquifer exchange and surface-water groundwater interactions. Conservative elements and deuterium/oxygen isotopes were used in a mixing cell model (MCM) approach to assess groundwater flow between aquifers. Efforts to accurately quantify and understand MSR are hampered by sparse data on the geochemical character of bedrock aquifers. Watershed scale recharge estimates and water balances were derived from a regional integrated climate dataset coupled to FEFLOW simulations. The first stage modeled steady state conditions within the main valley center aquifer. Integrated surface water and groundwater modeling is to be carried out in the future. The groundwater flow modeling will contribute to subsequent water management decisions at the watershed scale. Climate change and economic change scenarios will be considered in the integrated surface water and groundwater modeling.
The Development of Multiple Conceptual Models for the High Risk Saline Water Upconing Area in the Chi River Basin, Thailand
Toward a Distributed Groundwater Model of the Chesapeake Bay Watershed
The Chesapeake Bay is currently listed as an impaired water body and is the object of intensive restoration efforts to mitigate water quality issues, particularly those related to the transport of nutrients and sediments into the Bay. Groundwater discharge contributes significantly to the annual flows of Chesapeake Bay tributaries and is presumed to contribute to the observed lag time between the implementation of management actions and the environmental response in the Chesapeake Bay. The EPA Phase 5 Chesapeake Bay watershed model, based on the HSPF code, simulates river flow and associated transport and fate of nutrients and sediments and is used as a decision support tool. However, this model does not provide a mass-conserving representation of the groundwater fluxes in the watershed in that flow to deep groundwater is removed from model calculations. As an alternative, this study aims at the development of a fully distributed model of the Chesapeake Bay Watershed using PARFLOW, a three-dimensional variably saturated groundwater flow model integrated with an overland flow component. The code can be run on high performance computing facilities, thus allowing for applications at large scale with efficient run times. The CBW encompasses an area of 160,000 sq km and spans five physiographic provinces. We present results of testing that has been carried out to determine optimal gridding and parameterization to represent the topography and relevant subsurface processes across a range scales of the model domain. The development of this distributed large-scale model offers an opportunity of incorporating remote-sensing data such as the MODIS, LDAS and GRACE data products as well the NEXRAD precipitation product. These data sets have spatial resolution relevant to the regional scale of the model and represent an improvement of the modeling performances considering the challenges usually encountered in obtaining temporally and spatially consistent data for large scale hydrologic modeling.
Measuring and understanding total dissolved gas pressure in groundwater
Since dissolved gases are important to a number of aspects of groundwater (e.g. age dating, active or passive bioremediation, greenhouse gas fluxes, understanding biogeochemical processes involving gases, assessing potential impacts of coal bed methane activities), accurate concentration measurements, and understanding of their subsurface behaviour are important. Researchers have recently begun using total dissolved gas pressure (TGP) sensor measurements, more commonly applied for surface water monitoring, in concert with gas composition analyses to estimate more accurate groundwater gas concentrations in wells. We have used hydraulic packers to isolate the well screens where TDP is being measured, and pump tests to indicate that in-well degassing may reduce TDG below background groundwater levels. Thus, in gas-charged groundwater zones, TGPs can be considerably underestimated in the absence of pumping or screen isolation. We have also observed transient decreased TGPs during pumping that are thought to result from ebullition induced when the water table or water level in the well is lowered below a critical hydrostatic pressure.
Measuring Total Dissolved Gas Pressure at a Petroleum Plume Site
Groundwater contamination from petroleum hydrocarbons is ubiquitous across the country, in both urban and rural settings. Natural attenuation of petroleum contaminants may result in the production of gases (e.g. methane, carbon dioxide), in dissolved and potentially gas-phase form, which may affect the extent, persistence and remediation of petroleum hydrocarbon groundwater plumes. Current monitoring practices for gases in groundwater generally involve collecting water samples from wells or gas from gas-water separators during pumping tests, and subsequent analysis in the laboratory. Another potential option is the use of total dissolved gas pressure (TDGP) sensors, which can provide valuable real-time information on dissolved gas conditions while minimizing analytical costs. However, these have not been adequately tested or optimized for use in monitoring petroleum-contaminated groundwater. Preliminary testing of TDGP sensor measurement was performed on a selection of existing wells at a site in Ontario with previously-monitored groundwater contamination by petroleum hydrocarbons. TDGP was measured using a PT4 Tracker (Point Four Systems Inc., B.C.). Other properties such as dissolved oxygen and pH were also measured, and samples were collected and analyzed for major ions, metals, and various petroleum hydrocarbons. Results showed that 3 of the wells had contaminants, as well as elevated methane and dissolved iron. They also had lower nitrate and sulphate concentrations, but so did one uncontaminated well. The TDGP for these wells was elevated compared to background groundwater and compared to that expected for equilibration with the atmosphere. These higher values likely result from the microbial generation of dissolved methane. This data set suggests that natural biodegradation processes are occurring in the petroleum plume. However, some other wells also had elevated TDGP. They could indicate a septic plume, but the relatively low electrical conductivity (EC) is not supportive of this. It was also noted that for some wells, but not all, TDGP increased substantially following pumping, which may indicate that degassed stagnant water in the well needs to be replaced by fresher groundwater prior to TDGP measurement. These preliminary findings suggest that TDGP has the potential to provide real-time insight into where gas-producing reactions (in this case, likely methanogenesis) may be occurring in groundwater, which may be useful in assessing or monitoring natural attenuation of petroleum hydrocarbons. However, there are complicating factors that require further investigation.
Patterns of Fluid Circulation and Steam Generation in Caldera-Hosted Hydrothermal Systems
Steam formation is an important mechanism powering near surface phenomena in active hydrothermal systems (e.g., Yellowstone) and an established ore deposition mechanism in ancient equivalents (e.g., Creede). To gain insights into factors controlling steam formation and distribution in these systems, a series of steady-state numerical models were run on a hypothetical caldera-hosted system based on characteristics of a representative suite of calderas (e.g., Yellowstone, Valles, Creede). Base model conditions consisted of (1) a 10 km-wide caldera with a flat floor and rim height of 800 m; (2) a 500 C intrusion 1.5 km below the caldera centre; (3) a regional conductive heat flux twice continental average (0.10 W/m2); (4) host rock thermal conductivity of 2.5 W/m C, density 2650 kg/m3 and pore fraction 0.05. An impermeable intrusion was modeled with a 500 m wide surrounding region with a permeability (k) 10-3 m2 less than the system meant to represent a ductile region produced by elevated temperature (T > 350 C). The remainder of the system was given homogenous permeability. Cylindrical coordinates were used to represent caldera geometry. For these conditions, a minimum k = 10-15 m2 was required to achieve the target thermal condition of T approximately 220 C at 300 m below ground surface observed in active systems (e.g., Yellowstone). This model also resulted in a continuous steam plume originating at the intrusive contact that reached within 300 m of the surface along the edges of the caldera ~2 km from caldera centre. Models with k < 10-15 m2 produced steam, but at greater depths and failed to match the target conditions. Models with intrusion temperatures reduced by 20% shifted the steam plume toward the caldera centre and reduced its volume, but still achieved target conditions. Increasing intrusion temperature by 10% produced a second distinct plume at the caldera centre that also achieved target conditions. Increasing the rim height for these conditions produced the base (i.e., single plume) conditions. Resurgent doming up to 300 m was also modeled for the caldera floor. Increasing dome height shifted steam towards the caldera centre, narrowing and reducing plume extent at the intrusion edge while forming and increasing the extent of the plume at the centre.
Aquarius - A Modelling Package for Groundwater Flow and Coupled Heat Transport in the Range 0.1 to 100 MPa and 0.1 to 1000 C
Aquarius is a Windows application that models fluid flow and heat transport under conditions in which fluid buoyancy can significantly impact patterns and magnitudes of fluid flow. The package is designed as a visualization tool through which users can examine flow systems in environments, both low temperature aquifers and regions with elevated PT regimes such as deep sedimentary basins, hydrothermal systems, and contact thermal aureoles. The package includes 4 components: (1) A finite-element mesh generator/assembler capable of representing complex geologic structures. Left-hand, right-hand and alternating linear triangles can be mixed within the mesh. Planer horizontal, planer vertical and cylindrical vertical coordinate sections are supported. (2) A menu-selectable system for setting properties and boundary/initial conditions. The design retains mathematical terminology for all input parameters such as scalars (e.g., porosity), tensors (e.g., permeability), and boundary/initial conditions (e.g., fixed potential). This makes the package an effective instructional aide by linking model requirements with the underlying mathematical concepts of partial differential equations and the solution logic of boundary/initial value problems. (3) Solution algorithms for steady-state and time-transient fluid flow/heat transport problems. For all models, the nonlinear global matrix equations are solved sequentially using over-relaxation techniques. Matrix storage design allows for large (e.g., 20000) element models to run efficiently on a typical PC. (4) A plotting system that supports contouring nodal data (e.g., head), vector plots for flux data (e.g., specific discharge), and colour gradient plots for elemental data (e.g., porosity), water properties (e.g., density), and performance measures (e.g., Peclet numbers). Display graphics can be printed or saved in standard graphic formats (e.g., jpeg). This package was developed from procedural codes in C written originally to model the hydrothermal flow system responsible for contact metamorphism of Utah's Alta Stock (Cook et al., AJS 1997). These codes were reprogrammed in Microsoft C# to take advantage of object oriented design and the capabilities of Microsoft's .NET framework. The package is available at no cost by e-mail request from the author.
Evaluation of Zeolite Permeable Treatment Wall for the Removal of Strontium-90 from Groundwater
Experimental and modeling studies have been initiated to evaluate the potential performance of a permeable treatment wall comprised of zeolite-rich rock for the removal of strontium-90 from groundwater. Preliminary column studies were performed using a synthetic groundwater referenced to anticipate field conditions, with a focus on quantifying the competitive ion exchange among five cations (Na+, K+, Ca2+, Mg2+, and Sr2+). Variations of the column configurations addressed the effects of particle size and flow rates on removal efficiency. In general, kinetic effects were not significant for the test conditions. Ongoing studies are focused on the comparison of zeolites obtained from two sources: Teague Mineral Products (Adrian, OR) and Bear River Zeolite (Preston, UT), and on the possible influence of native soil that is mixed with treatment wall material during construction. The results of the performance assessment will support the possible deployment of a full scale treatment wall at a Western New York nuclear facility.
Incorporating Microtopography Into Measures of Wetland Hydroperiod
Hydroperiod is a measure of the temporal variability of land-surface inundation. Frequently, hydrologic analysis uses conceptualizations of hydroperiod that neglect the effects of land surface microtopography and describe land surface as smoothly-varying and deterministic. Wetland topography generally exhibits small-scale (meters to centimeters) fluctuations of land surface (microtopography) that can only practically be described probabilistically. These microtopographic fluctuations produce small-scale variations in land surface inundation and water depth (above or below land surface) that are important considerations for ecologic and hydrologic processes. Alternative forms of hydroperiod are described that convey information about partial areal inundation and water depth variations resulting from microtopography. These enhanced forms of hydroperiod can provide a fuller description of wetland hydrologic conditions compared to traditional approaches and can be useful tools for improved water resources management to meet specific hydroperiod targets for a given species of plant or animal. Using microtopographic datasets derived from wetlands in the Florida Everglades, the enhanced forms of hydroperiod are shown to contain a much greater richness of information than traditional approaches that neglect microtopographic considerations.