H33F-0516 1340h
Evaluating Least Absolute Deviation Regression As An Inverse Model In Groundwater Calibration
A new Least Absolute Deviation and Expectation Maximization (LAD-EM) FORTRAN procedure replaced the MODFLOW groundwater model Least Squares (LSG) parameter estimation inverse procedure. The LSG and LAD-EM MODFLOW inverse models were applied to a multi-layer 20,000 hectare watershed in Iowa. The hydraulic conductivity parameters computed with the LAD-EM MODFLOW inverse model agree with published data. The LSG MODFLOW model predicted a confined aquifer system while the LAD-EM MODFLOW model predicted an unconfined system. Literature documents the area as an unconfined system. The hydraulic conductivity parameters computed by both inverse models were used to create two forward basin models of a larger encompassing 600,000 hectare river basin in Iowa. Application of the LAD-EM MODFLOW model hydraulic conductivities to the river basin produced model heads that were within two meters of observed heads. Application of the LSG MODFLOW model hydraulic conductivities produced model heads that were 15 meters below observed heads. The LAD-EM MODFLOW model is an excellent tool for modeling multi-layer groundwater aquifers.
H33F-0517 1340h
Coupled Geochemical and Reactive Transport Modeling of Organic Contaminants in a Pyrite-Rich Aquifer
Although pH is recognized as a key factor influencing bacterial activity, existing groundwater transport models generally do not directly account for the effect of pH on the biodegradation of organic compounds. The purpose of this study is to develop a coupled reactive transport and geochemical model that explicitly incorporates the effect of spatial and temporal variations of the pH on the biodegradation of organic contaminants. The model consists of two modules: a transport module and a geochemical module. The transport module uses a Crank-Nicholson finite-difference formulation to solve the groundwater flow and transport equations for the hydrocarbon, dissolved oxygen, microbial mass and all reactive groundwater species influencing the hydrocarbon biodegradation and pH distribution. The geochemical module allows for the simulation of both kinetically defined as well as geochemical equilibrium reactions. The governing non-linear system of equations is solved using an iterative multi-step operator-splitting algorithm. Both modules account for heterogeneity in the definition of the hydrogeological and biochemical parameters. For demonstration, the model is applied to a hypothetical pyrite-rich aquifer contaminated with petroleum hydrocarbons. A commonly used practice for the remediation of aquifers contaminated with petroleum hydrocarbons is the delivery of oxygen for the enhanced aerobic biodegradation of the organic contaminant. However, the presence of pyrite may interfere with the intended purpose of the supplied oxygen, leading to undesirable side effects. Specifically, oxygen readily reacts with the sulfide minerals leading to depletion of oxygen and acidification of the subsurface environment and, subsequently, the inadvertent inhibition of the microbial activity. The developed coupled geochemical and reactive transport model is used to quantify these processes and assess the dominance of the various chemical reactions. Both abiotic and biotic pyrite oxidation kinetics are incorporated in the model. The impact of heterogeneity as well as key parameters on the fate and transport of the organic contaminant is also evaluated. The example is used to demonstrate how additional measures such as the injection of alkaline solution with the oxygen may optimize the remediation process.
H33F-0518 1340h
Coupled versus Decoupled Solution Approaches for Non-Linear Fluid Flow Processes in Heterogeneous Porous and Fractured Media
Many fluid-flow processes in the Earth's crust, such as multiphase flow or convection due to temperature and/or concentration gradients, are non-linear in nature. Studying these processes using numerical simulations is challenging. On one hand, numerical methods must be robust and able to deal with the non-linearities efficiently. On the other hand, they must be capable of resolving orders of magnitude variations in permeability and geologically complex structures that often occur in the Earth's subsurface. This usually requires high-resolution meshes, possibly with up to several million degrees of freedom. Traditionally numerical methods have solved such flow processes fully coupled, i.e. solving for the independent variables simultaneously using iterative techniques to account for the non-linearities. While these approaches have solved challenging problems, they have the disadvantage that the global solution matrices are ill conditioned and hence not always suitable for fast matrix solvers such as algebraic multigrid solvers. Furthermore, iterative techniques such as Newton's method may fail to converge. Numerical approaches that are capable of resolving geologically complex structures, for example the finite element method, require upwind-weighting schemes to model advection-dominated fluid flow. Such upwinding techniques, however, may reduce the geometric flexibility of the finite element method, fail to converge if the permeability varies over more than two orders of magnitude, or smear out shock fronts in advection-dominated flows. Here we present the solutions of a decoupled approach, based on a combination of finite volume and finite element methods, and a fully coupled approach, based on an upwind-weighted finite element method, for a variety of non-linear fluid flow problems. The results are compared for accuracy, robustness, and speed. They show that, in general, the decoupled approach is computationally more efficient and robust, because it does not require the use of costly iteration schemes. The decoupled approach can particularly well deal with orders of magnitude variations in permeability. Since higher-order finite volume methods are used in the decoupled approach to solve the advection equations, shock fronts in advection-dominated processes are retained to great accuracy. The exactness of the solutions of the decoupled approach can be improved even further by using internal Picard iterations or improved Euler time-stepping. Results for simulations of convection in sub-seafloor magmatic-hydrothermal systems, where temperature variations are in excess of 1000$^\circ$C, however, show that the gain in accuracy is only small and does not justify the increased computational costs.
H33F-0519 1340h
Mathematical Modeling of Fate and Transport of Aqueous Species in Stormflow Entering Infiltration Basin.
The State of California is evaluating the role of passive stormwater detention facilities for the purpose of attenuating potential dissolved and suspended chemical species that may originate in roadway runoff of rainfall. The engineering design of such infiltration basins requires tools to quantify their performance as recipients of stormwater runoff from roadways, and as filters of aqueous chemical species. For this purpose a one-dimensional unsaturated flow and transport model is developed to estimate the efficiency of storm-water infiltration basins in treating roadway generated metallic and organic pollutants. Kinematic wave approximation is used along with van Genuchten water retention model to simulate water percolation thorough the infiltration basin. For metals a Langmuir type nonlinear competitive sorption isotherm is used for transport of chemicals and a kinetic reversible linear sorption model is considered for organics. The model is applied to known roadway born metallic contaminations such as copper, zinc, lead, chromium, nickel and cadmium, as well as organic species such as diazinon, diuron, ghlyphosate and pyrene, for several representative soil and precipitation condition for California within a period of five years. Representative soil parameters and precipitation patterns are extracted from frequency distributions extracted from a recent study. In addition sensitivity analysis has been done to evaluate the effect of soil property values on the performance of infiltration basins. The results can be used to evaluate the performance of infiltration basins in improving the water quality as well as being used in providing guidelines in design and maintenance of infiltration basins.
H33F-0520 1340h
Calibration of a Complex Three-Dimensional Coastal Aquifer with Density-Dependent Flow
In Los Angeles County hydraulic barriers have been implemented in three coastal areas to reduce saltwater intrusion since the 1950s. Effective evaluation and operation of current barrier facilities is critical for protecting the future basin water supplies. Due to conservation concerns, and because of a general interest in improving the performance of the barriers in Los Angeles County, a modeling study has been initiated for one of the barrier sites. An open-source, fully three-dimensional finite element groundwater model has been selected to solve the coupled flow and mass transport equations. To ensure proper simulation of density-dependent flow and transport, the model was validated using Henry's problem and the output was quantitatively compared with standard published results. A modified kriging method was applied to estimate the heterogeneous intrinsic permeability distribution for each of the five aquifer layers, noting that a comparable level of detail could not be achieved with a more traditional inverse procedure. The model performance for the heterogeneous case was compared with a homogeneous case estimated by traditional inverse procedure. Furthermore, two calibration periods, one short term and one long term, were selected to fine tune the appropriate transport parameters. Typical monitoring intervals were biannual with a smaller subset of observation wells taking measurements of total head and chloride concentration on a monthly basis.
H33F-0521 1340h
A New and Efficient Space-time Sub-discretizaton Methodology for Concurrent Multi-scale Groundwater Modeling
Multi-scale simulations in groundwater often necessitate the use of successive telescopic discretizations, each designed for a given spatial scale. The use of telescopic grids or meshes resulting from such a discretization procedure requires that successive simulations in different scales be carried out, starting with the one covering the largest spatial scale. Using this procedure, boundary conditions for successively smaller-scale models must be derived from larger-scale models. This type of simulation is computationally demanding, labor-intensive, and time-consuming. To overcome these problems, a new methodology for sub-discretization has been developed. This methodology is relatively flexible, allowing a given three-dimensional hexahedral block or element to be sub-discretized with infinite combinations of number of subdivisions in three dimensions. For transient simulations in multi-scale environment, there may be hydrogeologic/anthropogenic features such as fractures or time-varying extraction (or injection) points that cause groundwater pressure or solute concentration to change rapidly with time locally. It is apparent therefore that fine temporal discretization should be confined to the areas where rapid changes are expected. To overcome the problem of non-uniform concurrent location-dependent time discretization requirements, a temporal sub-discretization with virtual nodes to account for local solutions at different time levels where finer time discretization is necessary has been developed. The two methodologies above can be used separately or synergistically combined to provide a powerful solution method for multi-scale simulations. The newly developed methodology of space and time sub-discretization will be presented along with application examples, results, and discussions.
H33F-0522 1340h
Upscaling Relative Permeabilities in a Structured Porous Medium
Upscaling of multi-phase flow problems for a heterogeneous porous medium requires modification of constitutive functions at the grid-block scale. A particular type of heterogeneity that has important environmental consequences involves thin, continuous streaks of high permeability through lower-permeability background rocks. These streaks, which may correspond to features like abandoned wells in mature sedimentary basins, can become preferential flow paths for an invading fluid. Quantification of flow through these types of heterogeneities in deep, geological formations is necessary to estimate the migration and possible leakage of injected fluids such as hazardous liquid wastes, municipal liquid wastes, and, possibly, carbon dioxide. One of the important constitutive functions for proper estimation of flow through these flow paths is the relative permeability function. In the simple case of a single high-permeability streak in a uniform rock matrix, with both materials having identical (local) relative permeability functions, the upscaled relative permeability must be changed significantly to capture the proper leakage due to the nonlinear coupling with phase saturation. Standard petroleum reservoir pseudo functions for relative permeability capture the general features of the upscaled function, but they still produce errors of several hundred percent in the leakage estimation. Detailed three-dimensional numerical simulations and associated upscaled calculations demonstrate the proper form for the upscaled relative permeability, and provide a new derivation of pseudo functions to capture the leakage behavior in larger-scale models.
H33F-0523 1340h
Analytic Element Solutions for Preferential Flow Around and Through Elliptical Inhomogeneities in the Vadose Zone
We present new analytic element solutions for flow around and through elliptical inhomogeneities in the vadose zone. The solutions are analytic, and thus there is no grid; boundary conditions along the elliptical inhomogeneities are met up to machine accuracy, provided that enough terms are used in the series solution. The pressure head, velocity and saturation may be computed analytically at any point in the vadose zone. The approach can handle an arbitrary number of elliptical inhomogeneities, and is limited only by the available computer power. In this presentation we will consider elliptical inhomogeneities that consist of coarser material than the surrounding media; the ellipses are embedded in an otherwise uniform downward flux. The hydraulic conductivity is represented by an exponential function of the pressure head (the Gardner model). Different functions are specified for the inside and the outside of an inhomogeneity; both the saturated hydraulic conductivity and the water retention parameters may differ between the inside and outside. The coarser-soil inhomogeneities may divert flow or contract flow, depending on the magnitude of the downward flux. The resulting preferential flow paths, in case of diversion also called funnel-type flow, may result in large variations of the velocity, much larger than for the equivalent case when flow is saturated. A comparison with a numerical solution will be presented to assess the grid refinement that is needed near inhomogeneity boundaries in finite element solutions of the same problem. Several initial studies of contaminant spreading through multiple elliptical inhomogeneities will be presented, clearly demonstrating the longitudinal spreading caused by advection through fields of ellipses, and the dependency of the longitudinal spreading on the magnitude of the downward flux.
H33F-0524 1340h
Nonlocal Analysis of Mean Non-Reactive Solute Transport in Bounded Randomly Heterogeneous Media
Solute transport in randomly heterogeneous media is described by stochastic transport equations that are typically solved by Monte Carlo simulation. A promising alternative is to solve a corresponding system of statistical moment equations directly. The moment equations are generally integro-differential and include nonlocal parameters depending on more than one point in space-time [Neuman, 1993; Zhang and Neuman, 1996; Guadagnini and Neuman, 2001]. We present recursive approximations, and a numerical algorithm, that allow computing lead ensemble moments of non-reactive solute transport in bounded, randomly heterogeneous media. Our recursive equations are formally valid for mildly heterogeneous aquifers with $\sigma$_{Y}$^{2} <$ 1 where $\sigma$_{Y}$^{2}$ is a measure of log-hydraulic conductivity variance. Our algorithm utilizes a finite element Laplace transform method (FELT) valid for steady state advective velocity fields. Computational results in two spatial dimensions compare well with Monte Carlo results when $\sigma$_{Y}$^{2}$ and the grid Peclet number are small. As these parameters increase, the quality of our moment solution deteriorates. In principle, one should be able to control $\sigma$_{Y}$^{2}$ through conditioning on data and the Peclet number by selecting appropriate space-time discretization intervals.
H33F-0525 1340h
Stochastic Analysis and Numerical Simulations of Hydraulic Head and Baseflow in Heterogeneous Media under Spatial-Temporal Random Recharge
Stochastic analysis and numerical simulations were carried out to study the temporal scaling in the time series of water table fluctuations in a one-dimensional heterogeneous aquifer under spatial-temporal random recharge. It was found in our previous study that scaling of water table fluctuations may exist and the fractal dimensions varies over space based on spectral analyses of the hourly hydraulic head (h) data observed over a four-year period at seven monitoring wells in the Walnut Creek watershed in Iowa. The estimated baseflow in the Walnut Creek and other four watersheds has temporal scaling, but there exits two distinct slopes with a break at about 30 days in the log frequency and log power spectral density plot. It was also found that the hydraulic head in an aquifer may fluctuate as a fractal in time in response to either a white-noise or a fractal recharge process, depending on how quickly the hydraulic head responds to recharge events and the physical parameters of the aquifer (i.e., transmissivity and specific yield). Numerical simulations were conducted to verify whether or not the hydraulic head and flux (baseflow) behave as fractal processes and if their fractal dimensions vary spatially, using a 1-D transient groundwater flow in heterogeneous aquifer subject to temporal and spatial random recharge (white noise in time and exponential covariance in space). The simulation results confirm the previous findings. We also derived the moment equations which were solved to obtain the mean hydraulic head and mean flux (baseflow). The spectrum for mean hydraulic heads and flux (baseflow) are plotted and analyzed to test our hypotheses and the effects of aquifer heterogeneity and spatial-temporal random recharge process on the head fluctuations and spectrum are presented and discussed.
H33F-0526 1340h
Numerical Advances in the Modeling of Highly Heterogeneous Domains Using the Analytic Element Method
The Analytic Element Method (AEM) is an alternative to Finite Element and Finite Difference Methods for solving subsurface flow and transport problems in highly heterogeneous domains. Each hydrogeologic element (e.g. a zone where conductivity differs from the surrounding conductivity, a well, a river, etc.) in the domain of interest is represented by an analytic function. In order to solve the flow problem, the coefficients of the analytic function for each element must be computed so as to satisfy the boundary conditions. For the iterative approach used in this presentation, most of the computational effort is spent on computing the interactions between elements. These interactions are reflected in the element coefficients. Two approaches that can greatly reduce the run-time are examined in this presentation: the superblock approach and the parallel processing approach. The superblock approach groups elements into blocks, and treats each group as a single element while computing the influence of the group's elements on elements that are far from the group. Several superblock formulations exist. For example, the superblocks can be arranged using a Quad-Tree approach to further reduce the computational effort. Two superblock formulations will be examined. Parallel processing also reduces the time needed for solving the flow problem. The computations of the element coefficients are carried out separately for each element during a given step. Interactions are resolved using an iterative algorithm. The elements are hence divided up among many processors, so that a large number of element coefficients can be computed at the same time. The processors only need to communicate the new element coefficients at the end of each iterative step. This enhances the efficiency and the elegance of the approach. The two approaches described above are examined separately and in combination on a set of test problems containing a large number of circular inhomogeneities. Performance and efficiency are evaluated using a set of numerical simulations of flow in highly heterogeneous formations. The results show that as many as 250,000 2D inhomogeneities (and 100,000 3D inhomogeneities) can be simulated regardless of the packing densities and the heterogeneity levels. Variances of log-conductivity up to 10 can be analyzed. In addition to circular (2D) and spherical (3D) shapes the methodology was applied to other shapes.
H33F-0527 1340h
A Sample-based Stochastic Approach For Modeling Of Fluid Flow Through Heterogeneous Unsaturated Fractured Rock
Modeling fluid flow and chemical transport processes in large-scale, three-dimensional, fractured reservoirs is conceptually difficult and computationally demanding. One key issue in the modeling studies of real field problems is how to represent heterogeneous rock properties of fractured media in numerical models. For most studies, large-scale spatial and temporal averaging is employed to represent a heterogeneous fracture and matrix system. An alternative approach is to prescribe the heterogeneous system stochastically, using measurements and calibration data. However, the stochastic approach, compared with the traditional deterministic approach, is computationally more demanding. In this study, we present a modeling approach to examine the effects of rock-property heterogeneity on flow within the unsaturated zone of Yucca Mountain, Nevada. The heterogeneity of each of the geological layers within the unsaturated zone system is represented using a sample-based stochastic distribution scheme. This scheme determines the rock properties of each gridblock based on the statistical information of field-measurement data for the corresponding geological unit. The rock properties are chosen randomly from one of the measured data groups. The frequency for each rock property is conditioned by statistical information from field-measurement data. Simulation results are compared to the traditional stochastic approach, which employs a spherical semivariogram model with empirical log permeability semivariograms to generate a 3-D spatially distributed rock-permeability field. Geostatistical parameters and cumulative distribution functions of rock permeability for the spherical semivariogram model are derived using the measured data.
H33F-0528 1340h
Discrete Analytic Domains: a New Technique for Groundwater Flow Modeling in Layered, Anisotropic, and Heterogeneous Aquifer Systems
A new technique for modeling groundwater flow with analytic solutions is presented. It allows modeling of layered aquifer systems with complex heterogeneity and anisotropy. As with previous AEM techniques, flow in each layer is modeled with two-dimensional analytic solutions, there is high accuracy and resolution, the domain is not discretized into grid blocks or elements, and the modeled area can easily be altered and expanded in the midst of the modeling process. This method differs from previous Analytic Element Method (AEM) techniques by allowing general anisotropy conditions in the aquifer. The flow field is broken into discrete polygonal domains, each with its own definition of isotropic or anisotropic aquifer parameters. An advantage of this approach is that the anisotropy orientation and ratio can differ from one domain to another, a capability not possible with the "infinite domain" of previous AEM formulations. With this approach, the potential and discharge vector functions at a point are the sum of contributions from elements within or on the boundary of the domain containing the point. Unlike previous AEM schemes, elements beyond the domain boundary don't contribute to these functions. Once a solution is in hand, less computation is required to evaluate heads and discharges because fewer elements contribute to the equations. This computational efficiency could prove a significant advantage in large regional models and in solute transport models that use a flow model's velocity field. This technique allows multiple aquifer layers to be stacked vertically, and it has the novel ability to have more layers in the area of interest than in distant areas. This feature can save significant computation and input effort by concentrating layering detail only where it is needed and warranted by data. The boundary conditions at line element boundaries are approximated, including the continuity of flow and head across heterogeneity boundaries. By using high-order line elements with as many as 15 degrees of freedom per element, the boundary condition approximations can be made quite accurate. An example model illustrates the method's capabilities and the accuracy of boundary conditions achieved.
H33F-0529 1340h
Transport Property Modeling in Partially-saturated Rocks Using Pore-scale Simulations
The earth sciences are undergoing a gradual but massive shift from description of the earth and earth systems, toward process modeling and simulation. This shift is very challenging because the underlying physical and chemical processes are often nonlinear and coupled. In addition, we are especially challenged when the processes take place in strongly heterogeneous systems. One example is multiphase fluid flow in rocks, which is a nonlinear, coupled and time-dependent problem and occurs in complex porous systems. To understand and simulate these complex processes, the knowledge of underlying pore-scale processes is essential. To this end, we have initiated computational rock physics to rigorously simulate rock/reservoir properties. The computational rock physics framework is based on digital representations of rocks, which consist of minerals and fluids, and may evolve with time. It also contains modular physical property simulators, with which we directly simulate physical properties of rocks. This computational environment significantly complements the physical laboratory: 1) rigorous prediction of the physical properties, 2) interrelations among the different rock properties using shared digital porous media, and 3) simulation of dynamic problems with multiple physical responses. As a continuing effort of this framework, we investigated transport properties in partially-saturated rocks. Two digital rock structures were used in this study - X-ray tomographic Fontainebleau sandstone and random dense pack of spheres. Partially-saturated rock samples were obtained through two-phase flow simulations. We used both steady-state and unsteady-state simulations to compare static and dynamic properties at different partial saturations. We then performed single-phase and electrical flow simulations to calculate permeability and electrical conductivity. We found that water phase percolates around Sw=20% and air phase loses connectivity around Sw=60-70%. We also observed electrical conductivity hysteresis between drainage and imbibition simulations. The amount of hysteresis shows strong anisotropy, which shows strong dependence on the orientation of fluid flow. However, angular averaged electrical conductivity did not show any significant hysteresis. The relation between peremeability and electrical conductivity varies with different rock geometry. It then is possible to characterize the relation between electrical conductivity and permeability in a given rock formation, and to provide rigorous links to understanding and modeling geological processes.
H33F-0530 1340h
Streamline based analysis of non-Fickian dispersion in porous media with multi-scale heterogeneity
A numerical method to model contaminant transport in heterogeneous geological formations is developed. The present model is based on the framework that takes into account of two different spatial scales of uncertainty, microscopic level and macroscopic level. The microscopic heterogeneity caused by the existence of a soil particle and a crack, which is considered smaller than the scale of the measurement resolution hence unresolved in the real application, is treated by a stochastic method, while the microscopic heterogeneity, which can be treated as a large scale trend and well defined, is solved deterministically by the streamline method. The stochastic method developed is based on the continuous time random walk (CTRW) formalism, which is numerically solved along a streamline. The application of the streamline method to the integration of CTRW in the space domain is based on the fact that the macro-dispersion tensor is usually observed highly anisotropic and that the transverse dispersivity can be negligible compared to the longitudal dispersivity. Since the integration in the space domain is assumed 1-D in the streamline method, the numerical integration can be simplified by applying the Laplace transformation in the time domain. In obtaining a streamline numerically, the mesh dependency observed in the classical finite element approach especially in the case of the unstructured grids has to be eliminated. In the present method, the element free approach is adopted. The accuracy of streamlines obtained by the element free method is evaluated through the comparison with the analytical solution for the repeated spot pattern. The calculations are compared with the experimental results, in which the statistical distribution of the conductivity is well controlled. From the comparison with the experimental result, not only the validity of the present model is discussed but also the effect of the anomalous transport characteristics to the macro-dispersion is clarified.
H33F-0531 1340h
Coupled Fluid-Heat Flow Model Analysis to Explain the Origin and Accumulation of Subsurface Carbon Dioxide: McElmo Dome Case Study
McElmo Dome contains a very large reservoir of natural carbon dioxide. Located within the Paradox Basin, southwestern Colorado, McElmo is among several large carbon dioxide reservoirs in the southwestern U.S. Estimates of total carbon dioxide in the reservoir exceed 2 gigatons, and over 14 megatons are piped annually from McElmo to petroleum fields in the west Texas Permian Basin, for enhanced oil recovery operations. This and other carbon dioxide accumulations in the region have been the subject of many studies since the 1930s, with several hypotheses offered for the origin of McElmo's carbon dioxide. The most cited explanation is thermal degradation of the Mississippian Leadville Formation, a carbonate system that extends over much of the western U.S. The Ute Mountain laccolith is only a few miles from McElmo. This large igneous intrusion penetrates the Leadville carbonates at McElmo, providing the setting for the thermal degradation mechanism. The purpose of this study is to evaluate quantitatively this conceptual model of carbon dioxide origin. We developed mathematical simulation models of coupled heat and fluid flow of the McElmo system. Thermal data suggest that temperatures were probably sufficient for thermal degradation (metamorphism) to generate carbon dioxide. However, uncertainties associated with the intrusion's properties and timing hinder efforts to evaluate whether other mechanisms may play a significant role as well. For instance, coupled fluid-heat flow history model results suggest carbon dioxide may also have been generated as a byproduct of hydrocarbon catagenesis in Mississippian strata.
H33F-0532 1340h
Flow Dimension Analysis of Pumping Tests in Geometric and Random Fields
The generalized radial flow approach is increasingly used to interpret hydraulic tests. The flow dimension of a hydraulic test can be estimated directly from the second derivative of pumping test drawdown versus log time. Previous studies have shown that flow dimension can be described as the change in cross-sectional area of flow and radial distance from the borehole. Using this relationship, we have developed a simple algorithm that generates geometries that correspond to arbitrary flow dimensions. We validated this approach using numerical pumping test simulations in fields mirroring the geometries coupled with curve fitting methods. Our results establish a direct link between physical heterogeneity and flow-dimension. To illustrate the usefulness of flow-dimension for interpreting information about heterogeneity, we simulated pumping tests in spatially correlated random fields. We noted significant changes in flow dimension when the advancing drawdown front encountered strong, contrasting heterogeneities. Using data from randomly located pumping wells, it may be possible to characterize the correlation length of the heterogeneous system. We believe flow dimension can be used with geologically-based conceptual hydrologic models to constrain heterogeneity characteristics in real aquifers.
H33F-0533 1340h
Optimal Mesh Generation for AEM-based Eulerian Transport Simulators
The analytic element method (AEM) is a grid-independent approach for simulating groundwater flow in shallow aquifer systems. Recent advances have enabled AEM flow solutions to be used as the basis for Eulerian (finite difference or finite element) contaminant transport simulators. One of the benefits of such a merger is the removal of the constraints imposed by the flow grid or mesh. The resultant model discretization may be specified to accommodate only the relevant transport constraints (i.e., the Peclet and Courant limitations). Design of the mesh is typically limited by discretization requirements of the flow problem. With the use of AEM flow solutions, grid and mesh geometry may be optimized for a specific transport system without regard for flow system discretization. A two-dimensional mesh generation algorithm is presented that maximizes the required node spacing (governed by Peclet limitations) and therefore the time step (governed by Courant limitations) required for transport models using the spatially continuous velocities and dispersion coefficients obtained from AEM flow solutions. The optimized meshes reduce the computational cost of contaminant transport models by reducing the total number of degrees of freedom. Results from the optimized mesh algorithm illuminate some artifacts of flow discretization that are commonly neglected. The increased computational efficiency of models simulated without these discretization artifacts is quantified, and some non-intuitive results concerning the optimal mesh design for transport simulation are presented.
http://www.groundwater.buffalo.edu
H33F-0534 1340h
A Coupled Systems Approach to Solute Transport Within a Heterogeneous Vadose Zone-Groundwater Environment
A coupled systems approach is presented for the simulation of fluid flow and solute transport within a heterogeneous vadose zone-groundwater environment. Separate model domains were developed for the vadose and saturated zones. The numerical model used is FEMWATER, a three-dimensional finite element model for simulating flow and transport in variably saturated media. A control plane, located at the spatial intersection of the model boundaries, couples the two systems together. The fluid flux and solute mass passing through the control plane during a vadose zone simulation is used to construct the recharge boundary condition in the saturated system. Coupling the two systems allows for the accurate representation of flow and transport within the vadose zone along with the ability to simulate regional plume migration in the saturated zone. The Recharge Boundary Construction Algorithm adaptively reconciles both the spatial and temporal discretization differences between the vadose zone and saturated zone models. The accuracy of the algorithm is determined by a set of user-defined parameters controlling the algorithm's adaptive spatial mesh and adaptive temporal time step routines. The routines are designed to capture the greatest resolution of spatial and temporal changes in fluid flux across the control plane during a vadose zone simulation. A sensitivity analysis performed on the control parameters for the construction of the recharge boundary shows the maximum percent difference in total mass passing through the control plane during a vadose zone simulation as having the greatest effect on the boundaries level of resolution. Numerical simulations are made investigating the impact of effluent recharge at the ground surface on the solute breakthroughs and distributions in the two systems.
H33F-0535 1340h
Simulations and Theory of Density-Dependent Dispersion in Weakly Heterogeneous Porous Media
The effect of density gradients on dispersive mixing of miscible fluids is studied. Density gradients not only affect dispersive transport for unstable flow, but also can significantly affect dispersion under stable conditions. For vertical displacements, laboratory experiments show a reduction of the longitudinal dispersivity under the influence of stabilizing density gradients. Linear Fick's law for the dispersive flux is inadequate to model these experiments. Therefore, alternative theories have been developed that incorporate the effect of density gradients. In this study, accurate numerical simulations of vertical displacements in weakly heterogeneous porous columns are performed. The numerical results confirm experimental observations. Furthermore, the computed concentration profiles are used to validate nonlinear dispersion theories applicable for high density gradients. A comparison is made with the stochastic theory of Welty and Gelhar, with homogenization theory, and with the nonlinear dispersion theory of Hassanizadeh and Leijnse. For small variances in $\log k$, the stochastic and homogenization theory lead to reasonably good predictions of the computed concentration profiles and variances, without any fitting. In the theory of Hassanizadeh and Leijnse a fitting parameter is involved. This nonlinear dispersion parameter is not a true medium parameter, as it is found to be dependent on the flow rate and on the travel distance.
H33F-0536 1340h
Local Discontinuous Galerkin Approximations And Variable Step Size, Variable Order Time Integration For Richards' Equation
Numerical simulation of Richards' equation continues to be difficult. It is highly nonlinear under common constitutive relations and exhibits sharp fronts in both the pressure head and volume fraction for many problems of interest. For a number of multiphase flow problems, the use of variable order and variable step size temporal discretizations has shown some advantages. However, the spatial discretizations commonly used for variably saturated flow are dominated by nonadaptive, low-order finite difference and finite element methods. Discontinuous Galerkin (DG) finite element methods have received significant attention in a number of fields for hyperbolic PDE's and, more recently, for elliptic and parabolic problems. DG approaches like the local discontinuous Galerkin (LDG) method are appealing for modeling subsurface flow since they can lead to velocity fields that are locally mass-conservative without the need for auxiliary variables or alternative meshes. DG discretizations are also inherently local and so better-suited for unstructured meshes and $h$-$p$ adaption strategies than traditional methods. While some work has been done recently for multiphase subsurface flow, there are a range of issues related to the performance of DG methods for highly nonlinear parabolic problems like Richards' equation that have not been investigated fully. In this work, we consider the combination of higher order adaptive time integration with an LDG spatial discretization for Richards' equation. We compare this approach to standard low-order methods for a series of test problems and consider a number of issues including the methods' relative accuracy and computational efficiency.
H33F-0537 1340h
Network Modeling of Anorthite and Kaolinite Reaction Rates in Porous Media
Although subsurface systems consist of porous media with a wide range of physical and chemical properties, reactive transport modeling commonly employs a continuum approach, where each grid block is characterized by spatially-averaged properties and reaction rates are calculated using uniform concentrations. Such spatial averaging can introduce significant error in the representation of geochemical reaction rates. In this work, we use network models to examine the effects of pore-scale heterogeneities on continuum-scale rates of anorthite and kaolinite reactions, and we identify conditions under which the effects of pore-scale heterogeneities are significant in reaction rate up-scaling. The network is constructed to represent consolidated sandstone, with about 10% reactive minerals (anorthite and kaolinite) distributed in reactive pore clusters according to prescribed configurations. The minerals react with acidic brine saturated with high-pressure CO2, representing conditions relevant for geological CO2 sequestration. The reaction rates computed from the network model are compared with those from a continuum model, which simulates the porous medium using spatially-averaged concentrations. Simulation results show that aqueous concentrations at the pore scale vary by orders of magnitude, and their distributions are highly skewed and, in some cases, bimodal. These concentration heterogeneities lead to significant spatial variations in pore-scale reaction rates. As a result, the continuum-scale rates from the network model are significantly different from those from the continuum model. In general, the continuum model overestimates the anorthite dissolution rates; for kaolinite, it either underestimates its precipitation rates, or predicts a different reaction direction from that of the network model. The effects of pore-scale heterogeneities are influenced by hydrodynamic conditions, as well as spatial distributions of reactive minerals. For anorthite dissolution, the continuum model overestimates its rates under all hydrodynamic conditions, with the degree of overestimation reaching a maximum in medium flow conditions. For kaolinite, the continuum model and the network model predict the same reaction direction at slow and fast flow conditions, but opposite reaction directions in medium flow conditions. With regard to the impact of reactive cluster size, larger cluster size leads to larger differences between the rates. These results provide guidelines on the conditions under which the effects of pore-scale heterogeneities are important in reactive transport modeling.
H33F-0538 1340h
Simulating Contaminant Transport In Heterogeneous Media Using The Dual Porosity Method
The dual porosity (DP) formulation is a method originally developed for combined fracture and matrix flow. The typical approach would be to solve the normal system of mass and transport equations on a fully-connected grid representing the fracture domain, with one additional equation at each node which corresponds to the matrix domain. In this formulation, the matrix is not a continuous medium, but provides important fluid and energy storage terms to the fracture node equations. A significant limitation of this approach is that there is no ability to capture gradients within the matrix; the original dual porosity method was sometimes called "quasi-steady" because of this limitation. Since the early development of DP methods, advances have been made in solution techniques and also in increasing the number of matrix modes corresponding to each fracture node. We introduce herein a Generalized Dual Porosity Method (GDPM) in which the number of matrix nodes (per fracture node) can vary spatially within a numerical grid. This allows not only the ability to represent gradients and mass -transfer limited processes within the matrix, but also to place additional matrix nodes efficiently. An algebraic decomposition method is presented in which the added complexity of the additional matrix nodes is efficiently solved. We apply a modification of the method to the problem of incorporating sub-grid scale heterogeneity in a sand and clay aquifer into large-scale flow and transport simulations. Here the sand and clay represent the continuous and discontinuous media, respectively. We assume that the mass of contaminant is primarily contained within the clay layers, which is transported via mass-transfer limited mechanisms to the sand medium. We compare several different numerical formulations of this problem: a very fine grid continuum model (truth), and coarse-grid GDPM solution, and intermediate scale continuum model. The models are compared on the basis of accuracy and computational efficiency.
H33F-0539 1340h
Numerical Investigation of Multiple-, Interacting-Scale Variable-Density Ground Water Flow Systems
The goal of our study is to elucidate the nonlinear processes that are important for multiple-, interacting-scale flow and solute transport in subsurface environments. In particular, we are focusing on the influence of small-scale instability development on variable-density ground water flow behavior in large-scale systems. Convective mixing caused by these instabilities may mix the fluids to a greater extent than would be the case with classical, Fickian dispersion. Most current numerical schemes for interpreting field-scale variable-density flow systems do not explicitly account for the complexities caused by small-scale instabilities and treat such processes as "lumped" Fickian dispersive mixing. Such approaches may greatly underestimate the mixing behavior and misrepresent the overall large-scale flow field dynamics. The specific objectives of our study are: (i) to develop an adaptive (spatial and temporal scales) three-dimensional numerical model that is fully capable of simulating field-scale variable-density flow systems with fine resolution (~1 cm); and (ii) to evaluate the importance of scale-dependent process interactions by performing a series of simulations on different problem scales ranging from laboratory experiments to field settings, including an aquifer storage and freshwater recovery (ASR) system similar to those planned for the Florida Everglades and in-situ contaminant remediation systems. We are examining (1) methods to create instabilities in field-scale systems, (2) porous media heterogeneity effects, and (3) the relation between heterogeneity characteristics (e.g., permeability variance and correlation length scales) and the mixing scales that develop for varying degrees of unstable stratification. Applications of our work include the design of new water supply and conservation measures (e.g., ASR systems), assessment of saltwater intrusion problems in coastal aquifers, and the design of in-situ remediation systems for aquifer restoration. We present preliminary model results for high-resolution simulation of variable-density flow and transport in homogeneous and heterogeneous porous media. We explicitly solve the three-dimensional advection equation using mass-conservative, flux-integral techniques and finite-volume formulations that provide unrestricted time-step capabilities similar to those associated with semi-Lagrangian methods. Our implementation of B. P. Leonard's MACHO (Multidimensional Advective-Conservative Hybrid Operator) and COSMIC (Conservative Operator Splitting for Multidimensions with Inherent Constancy) methods is an Nth-order (e.g., 7th-order or higher) advection scheme that significantly reduces numerical dispersion and can be adapted spatially and temporally as the simulation progresses. The ability of these higher-order methods to yield accurate, nonoscillatory concentration profiles is illustrated and compared to traditional implicit solution methods such as central and upwind differencing, and van Leer flux limiters. We also show preliminary results from our implementation of adaptive mesh refinement (AMR) techniques and discuss the interrelationship between AMR and the Nth-order advection schemes.
H33F-0540 1340h
2-D Numerical Modeling of a Fault Zone Leaking Carbon Dioxide in East Central Utah
The Little Grand Wash fault zone near the city of Green River, Utah, leaks carbon dioxide-rich waters to the surface. Previous studies indicate that the inorganic carbon dioxide travels from depths greater than 1 km to charge a shallower groundwater system that contains at least three aquifers. The aquifers leak fluids to the surface along the fault. We developed a 2-D mathematical model of the system, including coupled heat and multiphase fluid flow, to test and evaluate this conceptual model. Our goals of this effort include estimating quantitative partitioning of carbon dioxide amongst its liquid and gaseous phases with depth, inferring fault and matrix permeabilities, and constraining fluid fluxes in the subsurface system. Boundary conditions assigned in the model include observed surface fluxes and chemical data. Contributing artesian effects on the flow system have also been evaluated.
H33F-0541 1340h
Multicomponent reactive transport in hydrodynamic and geochemical heterogeneous porous media
Methods for analysing groundwater flow and transport in heterogeneous porous media usually focus on the heterogeneity of hydrodynamic parameters. Recent studies have addressed the effects of spatial heterogeneity of transport parameters for a single reactive species. On the other hand, sophisticated deterministic numerical models have been developped for the analysis of multicomponent reactive solute transport. Most of these models account for spatial heterogeneity by parameter zonation. Here we present the results of Montecarlo simulations for the transport of multicomponent chemically-reactive systems through hydrodynamic and geochemical heterogeneous porous media. One-dimensional transport of a set of chemical species suffering cation exchange takes place in a column of porous medium having random log-CEC (cation exchange capacity). This case corresponds to a laboratory column documented in PHREEQM user's manual. The column, initially filled with 1 mM NaNO3 and 0.2 mM KNO3, is flushed by 0.6 mM CaCl2 solution. Simulations of reactive transport have been performed with CORE2D V4. Statistical analyses have been obtained for the spatial distributions of cations and breakthrough curves at points located at increasing distances from the inflow boundary. log-CEC is assumed to be a random Gaussian function with a spherical semivariogram. Here we report the results for three cases: 1) Only heterogeneity in log-CEC with mild, medium and large variances (0.033, 0.1 and 1); 2) Only heterogeneity in log-K with variances of 0.01,0.1 and; and 3) Cross-correlated log-CEC and log-K. The heterogeneity in log-K affects the transport of all reactive and nonreactive species. Its effect increases with increasing log K variance. All species are "retarded" and show longer tails. The variability of CEC has a much smaller effect than the variability of log K. The combined effect of simultaneous heterogeneity in log-K & log-CEC produces much longer tails (especially for Ca) and larger concentration variances. The lack of cross-correlation causes a spread greater than that of possitively correlated K and CEC.
H33F-0542 1340h
Investigations of Observed Water Level Rise in the Culebra Aquifer; WIPP Site, Carlsbad, NM
The Waste Isolation Pilot Plant (WIPP) in southeast New Mexico has been developed for underground disposal of transuranic waste in halite beds of the Permian Salado Formation. Managed by the Department of Energy (DOE), the WIPP has been operational since March 1999. The most important water-bearing unit of the Salado Formation is the Culebra Aquifer, which lies about 200 m below ground surface and 400 m above the repository. In the assessment of compliance monitoring parameters for the year 2000, freshwater heads were compared to trigger value ranges established for 28 monitoring wells in the Culebra. Of these 28 measurements, freshwater heads in 21 wells appeared to be outside the trigger value ranges, with 20 higher and one lower than expected. Head changes in four of the wells could be explained by problems with well casings and/or leaking packers, leaving 17 wells with unexpectedly high freshwater heads. Hydrographs over time show a steady and consistent head level rise in these wells for the last 15 to 20 years. Two possible scenarios for the rise in heads have been formulated: (1) leakage into the Culebra of refining process water discharged onto potash tailings piles, probably through subsidence-induced fractures and/or leaky boreholes; and (2) leakage into the Culebra of water from units above or below the Culebra through poorly plugged and abandoned boreholes. This research examines the plausibility of these scenarios by separately modeling each case using a different mechanism for adding water to the Culebra Aquifer. For leakage from potash tailings piles, the model assumes recharge from a single potash tailings pile known as the Mississippi East site, which is located 10 to 12 km due north and up-gradient of the WIPP site. Storage coefficient and recharge rate from the tailings pile are calibrated to linearized hydrographs from 13 of the 28 monitoring wells using PEST. The second conceptual model identifies 26 boreholes as probable candidates for leakage to the Culebra, and models the boreholes as point sources to the Culebra. Like the first model, flow rates from the boreholes into the Culebra are calibrated to the linearized hydrographs. Calibrations for both models are conducted across 100 realizations of previously calibrated transmissivity fields, allowing for statistical analysis of the results. Results show that calibrations are generally better for the leaky borehole scenario than for the tailings pile scenario. This is most likely due to the distributed nature of the leaky boreholes versus the large point source of the single potash tailings pile. However, each scenario is able to provide very close fits to the hydrographs under certain transmissivity fields. Analyses are not able to exclude one scenario over the other, indicating that the water level rise is most likely due to a combination of sources. Investigations into the relative importance of each source are currently underway.
H33F-0543 1340h
18 Years Later: Revisiting a Groundwater Model of the Cambric Site at NTS
Since its advent in 1974, the Radionuclide Migration Project at the Nevada Test Site has spawned several interesting groundwater modeling ventures. Of interest to this research is the Cambric detonation site, where a tracer test was conducted from 1975 to 1991. Burbey and Wheatcraft (1986) built a groundwater/transport model of the Cambric site and at the time of calibration had achieved a good match to the measured data. Since then the predicted concentrations have diverged from the measured concentrations, which exhibit classic heavy-tailed behavior. It has been hypothesized that the Fractional Advection Dispersion Equation (FADE) will better predict these late-time high concentrations; this research will apply the FADE to the Cambric problem and aims to reach a more complete understanding of the physical significance of the coefficients contained in the FADE. We first built a preliminary groundwater model, employing the traditional Advection Dispersion Equation, in the hopes of duplicating Burbey's predicted concentrations. Burbey used the Deep Well Disposal Model, whereas this investigation used MODFLOW and MT3D. While the new model has produced a breakthrough curve fitting the peak concentration, it too fails to produce the heavy tail seen in the measured data. Also of concern is the nonuniqueness of the new model's solution; the best-fit breakthrough curve can be produced by changing either one of at least two parameters. We believe that both of these shortcomings (under predicted late-time concentrations and non-uniqueness) may be resolved by using the FADE. Not only does fractional theory permit heavy tails, but also it effectively replaces aquifer heterogeneity with fractional derivatives, thereby reducing the probability of a nonunique solution. Future work includes modeling the Cambric problem with Tadjeran and Meerschaert's numerical, fractional, radial-flow transport code (2003) and evaluating the code's applicability to varied flow and transport conditions.