A climate driven water resources model

Authors:

1. David Yates, National Center for Atmospheric Research
2. Kathleen Miler, NCAR
3. David Purkey, Stockholm Environment Institute-US Center

A climate–driven water resource planning model of California’s Sacramento River Basin (SACB) is presented, based on the Water Evaluation and Planning model Version 21 (WEAP21). The model’s calibration, validation, limitations, and results are presented. The major contribution includes a seamless integration of the watershed’s surface and sub-surface hydrology, consumptive and non-consumptive demands, and the water management infrastructure and governance that determine how water naturally flows and is managed. Both supply and demand side interactions are addressed simultaneously, allowing for analysis of alternative and/or future climate scenarios that are unbounded by a reliance on historical hydrologic patterns, with an example shown for the Sacramento river basin.

top of page

A Comparison of Computational Issues for Parameter Estimation using Finite Differences and the Adjoint State Method with MODFLOW and MT3DMS

Authors:

1. Tom Clemo, Boise State University
2. Michael Fienen, U.S. Geological Survey

This presentation reports our experiences using both finite differences and the adjoint state method to calculate parameter sensitivities. Neither method is perfectly accurate nor is the finite difference method always continuous with respect to small parameter changes. We explore the influence of those defects on the estimation of hydraulic conductivity of an isthmus separating Crystal and Muskellunge Lakes at the Trout Lake Field Station, Wisconsin. We also explore the influence of intentional manipulation of parameter sensitivity on the parameter estimates. The parameter estimation problem is formulated using a Bayesian approach similar that proposed by Kitanidis (1995). The approach treats hydraulic conductivity as a distributed property that may vary at the scale of the individual model cells. The sensitivity of simulated measurements to nearly every cell in the model is required. Measurements used to constrain the model include hydraulic head, Oxygen 18 fraction, and tritium concentration. Simulation of flow and transport is performed using the MODFLOW-2005 and MT3DMS programs sequentially. An implementation of the adjoint state method for the combined MODFLOW/MT3DMS model has been developed following the approach presented by Samper and Neuman (1986). For a specific change in the hydraulic conductivity, the finite difference calculation is perfectly accurate. The change in simulated measurements is what it is. However, the parameter estimation problem is non-linear making the finite difference sensitivity an approximation to both the parameter change of an update step and to very small parameter changes. The adjoint state method estimates the sensitivity as a continuous derivative of a measurement with respect to the conductivity of a cell. The accuracy of the adjoint state method is limited by the accuracy of the adjoint state solution. The most significant aspect of using the adjoint state method instead of the finite difference approach is over two orders of magnitude reduction in computation time of the sensitivities for our problem. More subtle aspects of the accuracy issues influence the problem through the matrix of sensitivities for the measurements. While the sensitivity matrix is not directly used to establish an acceptable solution of the parameter estimation problem, it does influence the solution through selection of update direction and the final stopping point due to allowances for measurement error. The matrix of sensitivities is also used in the calculation of posterior covariance, and more precision is expected at this point in the algorithm. Nonetheless, significant computational savings may be achieved by limiting the number of times that the most highly accurate sensitivities must be calculated. We explore creative debasement of parameter sensitivities in two ways. One is to reduce very large sensitivities which can lead to improved results. The second is to truncate the adjoint state calculation. This cause a trade-off between computation time and accuracy of the sensitivity calculation.

top of page

A Coupled Biogeochemical Reactive Transport Model for a Riverine Water-Column Benthic sediments

Authors:

1. Arash Massoudieh, Civil and Environmental Engineering, University of California, Davis
2. Fabian A Bombardelli, University of California, Davis
3. Timothy R Ginn, University of California, Davis

Sediments play a major role in fate and transport of highly sorbing metallic contaminants in water bodies. Sediment-associated contaminants can be released gradually due to diffusive processes, or may be re-entrained in the watercolumn as a result of high flows or anthropogenic activities. In addition biotic and abiotic geochemical transformations in the bed sediments control both the mobility and toxicity of sediment-associated contaminants. Therefore to quantify the fate and transport of contaminants with high affinity to solid particles in water bodies, various processes ranging from hydrodynamic and sediment transport to biogeochemical transformation of contaminants in the sediments need to be considered. In the research being presented, a multi-scale, quasi-two-dimensional, biogeochemical reactive transport model is presented, with description of both theory and numerical implementation. The model includes aqueous water column, suspended sediment, aqueous bed porewater, and sediment bed phases, with associated contaminant transformation and transport, as well as sediment resuspension, deposition and burial. The sediment bed domain is modeled using a set of vertical, Lagrangean (moving with the bed surface) one-dimensional sub-models which take into account burial and compaction as well as diffusive transport and biogeochemical reactions affecting the contaminants. Reactions include speciation, sorption/desorption, and microbially-mediated multiple terminal-electron accepting processes as kinetically-controlled redox reactions. The results of application of the model to simulate the fate and transport of Mercury and its methylation due to activities of sulfate and iron reducing bacteria in Colusa Basin Drain, California will also be presented.

top of page

A Coupled Land Surface-Subsurface Biogeochemical Model for Aqueous and Gaseous Nitrogen Losses Induced by Fertilizer Application

Authors:

1. Chuanhui Gu, Berkeley water center, UC berkeley
2. Federico Maggi, Berkeley Water Center, UC Berkeley
3. William Riley, Lawrence Berkeley National Laboratory
4. Curt Oldenburg, Lawrence Berkeley National Laboratory
5. Norm Miller, Lawrence Berkeley National Laboratory

A Coupled Land Surface-Subsurface Biogeochemical Model for Aqueous and Gaseous Nitrogen Losses Induced by Fertilizer Application Chuanhui GuTP PT, F. Maggi,P1P W.J. RileyTP PT, N.L.Miller2, C.M. OldenburgP2P Abstract In recent years concern has grown over the contribution of nitrogen (N) fertilizers to nitrate (NOB3PB-P) water pollution and atmospheric pollution of nitrous oxide (NB2BO), nitric oxide (NO), and ammonia (NHB3B). Characterizing the amount and species of N losses is therefore essential in developing a strategy to estimate and mitigate N leaching and emission to the atmosphere. Indeed, transformations of nitrogen depend strongly on water content, soil temperature, and nitrogen concentration. Land surface processes therefore have to be taken into account to properly characterize N biogeochemical cycling. However, most current nitrogen biogeochemical models take the land surface as the upper boundary by lumping the complex processes above the surface as known boundary condtions. In this study, an extant subsurface mechanistic N cycle model (TOUGHREACT-N) was coupled with the community land model (CLM). The resulting coupled model extends the modeling capability of TOUGHREACT-N to include the important energy, momentum, and moisture dynamics provided by CLM. The coupled model showed a significant impact of land-surface diurnal forcing on soil temperature and moisture and on nitrogen fluxes. We also discuss field applications of the model and discuss how temporal dynamics of nitrogen fluxes are affected

top of page

A COUPLED MODELING APPROACH FOR INCORPORATING VARIABLY SATURATED WATER FLOW AND SOLUTE TRANSPORT IN GROUND WATER MODELS

Authors:

1. Navin Kumar, UC Riverside
2. Jirka Simunek, UC Riverside

Owing to the constraints on computational resources, incorporation of variably saturated water flow and solute transport in to ground water models has posed many challenges for ground water modelers in the past. Simplified representation of vadose zone processes in ground water models is a very attractive approach but leads to poor modeling efficiency. A balance between modeling efficiency and computational demand of the adopted approach is thus needed to reasonably represent the effects of vadose zone processes. Recently, the HYDRUS package for MODFLOW has been developed that has been shown to represent vadose zone flow processes in ground water models at different spatial and temporal scales. Here, we extend the capabilities of the HYDRUS package to additionally represent the solute transport in vadose zone. Along with modeling water flow, the new HYDRUS package simulates solute transport such that the MODFLOW-HYDRUS code produce concentrations as a function of time that can be incorporated into the source function for MT3D. The capabilities of the HYDRUS package is tested using examples that represent a variety of spatial and temporal scales.

top of page

A Double Continuum Approach for Two-Phase Flow in Porous Media

Authors:

1. Leopold Stadler, Chair of Water Resources Management and Modeling of Hydrosystems, TU Berlin, Germany
2. Julia Mayer, Chair of Water Resources Management and Modeling of Hydrosystems, TU Berlin, Germany
3. Reinhard Hinkelmann, Chair of Water Resources Management and Modeling of Hydrosystems, TU Berlin, Germany
4. Rainer Helmig, Department of Hydromechanics and Modeling of Hydrosystems, Universität Stuttgart, Germany

The top soil layers of natural slopes are often highly porous when they consist of numerous macropores. During heavy rainfall, such structures lead to a fast by-passing of water through low permeable layers. Tracer experiments in the field at the natural slope Heumöser Hang in Austria showed the effect of such structures in a fast response of springs on the natural slope. Further investigations figured out that the fast response of the springs was even generated by a complex network of macropores. Such fast infiltration processes play a crucial role for the creeping and failure of hillslopes. For a simulation as close as possible to reality, it is required to describe the complex process interactions of the infiltration with the highly heterogeneous porous structures by accounting for the soil matrix as well as for the macropores. Consequently, this contribution deals with the development of a double continuum model for two-phase flow in porous media. The soil is subdivided into two interacting continua, one representing the soil matrix and the other the macropores. The flow in every continuum is determined by the solution of the two-phase flow equations for two immiscible fluids in a porous medium. Thus, the model concept accounts for the very different flow processes in the soil matrix as well as in the macropores. The exchange fluxes between the continua are described along the macropore surfaces based on a Darcy-type formulation. Therefore, the flow over the macropore surfaces is included in the mass balance equations as a further boundary. Among other things, it is possible to describe the above-metnioned by-passing via macropores during the infiltration with the presented model concept. The double continuum approach is specially suitable for simulations on large scales. The approach leads to a complex system of four coupled non-linear equations. It has been implemented in the modeling system MUFTE-UG. MUFTE is the abbreviation for Multiphase Flow Transport and Energy model and UG for unstructured grids. MUFTE contains physical concepts and discretization methods for multiphase flow and transport processes in porous and fractured-porous media, whereas UG provides data structures and fast solvers based on parallel, adaptive Multigrid methods. The idea of the new approach can be easily transfered into other existing simulation software for flow and transport processes in porous media including modifications for the specific problems. To test the newly developed double continuum approach, several examples have been analysed. As a first example, the displacement of water by a DNAPL (Dense Non-Aqueous Phase Liquid) is investigated because the density difference of the two fluids is small and therefore, this example is numerically \'easier\' compared to other two-phase systems such as water-air. The second example deals with a simplified test case of water infiltration and by-passsing via macropores in an unsaturatedsoil. A third example will treat water infiltration in a natural slope. In the range of possibilities, a comparion of numerical simulations and field measurement will be carried out.

top of page

A Fracture-Only Reservoir Simulator with Physically-Based Transfer Functions

Authors:

1. Evren Unsal, Imperial College London
2. Huiyun Lu, Imperial College London
3. Stephan K Matthai, Imperial College London
4. Martin J Blunt, Imperial College London

We propose a simulation methodology that combines the strengths of discrete fracture models with conventional dual porosity simulation. Constructing a grid and solving for flow in both fracture and matrix in a discrete fracture model is frequently so computationally demanding that only small systems can be studied. In contrast, while dual porosity models are more computationally efficient and can be applied at the field scale, they average the fracture properties and the transfer of fluids between fracture and matrix. In our approach we capture the complex geometry and connectivity of the fractures through explicit gridding of the fracture network. However, to avoid the prohibitive computational cost associated with gridding both the fracture and the matrix, we apply transfer functions to accommodate the flow of fluids between these two domains. We use a physically-based approach to modeling the transfer that overcomes many of the limitations of current formulations. The model is based on CSMP, an object-oriented discrete fracture simulator. Only transport in the fractures is considered with fracture/matrix transfer accommodated by source/sink terms associated with each fracture grid block. We discuss how to assign the transfer to account for the physical properties and geometry of the system. We validate the method through comparison with one-dimensional analytical solutions and comparison with experiments and simulations where both fracture and matrix are represented. We then present three-dimensional simulations of multiphase flow in a geologically realistic fracture network. To predict recovery it is necessary to represent the fracture network and the fluid displacement accurately. We discuss the strengths and limitations of our approach and the circumstances under which it is valid.

top of page

A global pressure approach for modelling compressible multiphase flow in heterogenous porous media

Authors:

1. Raphaël DI CHIARA, Institut de Mécanique des Fluides et des Solides 2,rue Boussingault 67000 Strasbourg FRANCE
2. Gerhard SCHÄFER, Institut de Mécanique des Fluides et des Solides 2,rue Boussingault 67000 Strasbourg FRANCE
3. Michel QUINTARD, Institut de Mécanique des Fluides de Toulouse Av. du professeur Camille Soula 31400 Toulouse FRANCE
4. Philippe ACKERER, Institut de Mécanique des Fluides et des Solides 2,rue Boussingault 67000 Strasbourg FRANCE
5. Guy CHAVENT, Institut National de Recherche en Informatique et Automatique Rocquencourt, Le Chesnay, France
6. Jean-Marie CÔME, BURGÉAP, 19, rue de la Villette, 69425 Lyon, FRANCE

A method for the simulation of compressible three-phase flows in heterogeneous porous formations is proposed taking into account gravity, capillary effects and heterogeneity. Governing equations are written in a fractional flow formulation in terms of a global pressure equation and two saturation equations. A global pressure-saturation formulation is much more efficient from the computational point of view than classical approaches in the case where three-phase relative permeability and capillary pressure curves satisfy a total differential (TD) condition as introduced by Chavent & Jaffré (1985). Hence, approximations using current approaches have to deal with two capillary pressure gradients coming from the saturation equations and would therefore require the resolution of four linear systems. However, due to the TD condition, the global formulation requires only the solution of two systems. In the developed approach, the solution of the global pressure equation requires the knowledge of the corresponding three-phase relative permeabilities and the global capillary pressure satisfying the TD condition. In absence of experimental three-phase data, regularity, monotonicity and bound constraints ensure that relative permeabilities are “physical”. Thus, a constrained optimization procedure is used to determine the preliminary secondary variables of the fractional flow from the effective saturation ternary diagram. The obtained three-phase relative permeabilities are compared to those obtained by Stone\'s model. Two efficient numerical methods are used to solve the global pressure equation and the two saturation equations for the water and oil phase. Discontinuous finite elements are used to approximate the convective term of the two saturation equations while the mixed finite element method (with mass lumping) is chosen to solve the global pressure equation and the diffusive part of the saturation equations. In a first test series, we will compare the fully global approach to the pseudo-global approach of Chen and Ewing (1996). Furthermore, the numerical results are compared with analytical two-phase solutions and runs obtained with the multiphase multicomponent simulator SIMUSCOPP (Le Thiez and Ducreux, 1997).

top of page

A level set method for non-zero contact angle drainage and imbibition in realistic porous media

Authors:

1. Masa Prodanovic, The University of Texas at Austin
2. Steven L. Bryant, The University of Texas at Austin

An accurate description of the mechanics of pore level displacement of immiscible fluids could significantly improve the macroscopic parameter predictions such as capillary pressure saturation - interfacial area relationships in real porous media. Slow displacement can be modeled as a quasi-static, capillarity-controlled process. At constant pressure and interfacial tension, pore scale fluid-fluid interfaces are modeled as constant mean curvature surfaces. For a given location in the pore space geometry, such surfaces need not be stable (in which case fluid typically moves until it finds an appropriate location) nor easy to calculate. Further, tracking the topological changes of the interface, such as splitting or merging, is nontrivial. The irregular pore spaces in natural porous media admit a multiplicity of such changes. Assuming quasi-static displacement, we recently developed ([1,2,3]) a simple but robust model based on the level set method for determining fluid interface position during both drainage and imbibition. The original method assumed zero-contact angle and arrived at geometrically correct interfaces. At the same time, the level set based description allows for robust handling of topology changes and independence from the pore space complexity. We now show the extension of the method to handle non-zero contact angles. While slightly more computationally involved than the original, the method preserves the benefits of the level set based description of phases: our simulations establish the exact position and shape of the interface in porous geometries, from which fluid volumes, contact areas and interface curvatures can be readily obtained. We show many examples of 2D and 3D displacements in individual pores (openings) and throat (narrow constrictions), as well as simulations in imaged porous samples and comparison to experimental results. Haines jumps at drainage and pore body imbibition events are readily identified, as well as snap-off and trapping of the non-wetting phase. [1] M. Prodanovic and S. L. Bryant. A level set method for determining critical curvatures for drainage and imbibition. Journal of Colloid and Interface Science, 304 (2006) 442--458. [2] M. Prodanovic and S.L. Bryant. In Focus on Water Resource Research ed. E. Heikkinen, Nova Science Publishers, Hauppage, New York chap. Resolving Meniscus Movement Within Rough Confining Surfaces Via the Level Set Method. (in press) [3] M. Prodanovic and S.L. Bryant. Physics-driven Interface Modeling for Drainage and Imbibition in Fractures. SPE paper 110448, Proceedings of SPE Annual Technical Conference and Exhibition, Anaheim, California, U.S.A., November 2007.

top of page

A mass conservative numerical scheme for reactive solute transport in soil

Authors:

1. Florin Adrian Radu, UFZ-Centre for environmental research
2. Sabine Attinger, UFZ-Centre for environmental research

Environmental pollution of groundwater and soil by organic compounds is nowadays a serious and widespread problem. A comprehensive active remediation often is not feasible from technical or financial reasons. Alternatively, in situ bioremediation (natural attenuation) has been recognized as a promising approach to restore sites contaminated with organic pollutants because it is less costly than active remediation strategies, the contaminants can ultimately be transformed to innoucous by-products with the help of microorganisms (not just transferred to another phase or place) and it can operate in situ. However, the decision to apply natural attenuation at a specific site depends essentially on the reliable prediction of the fate of the contaminant plume. Together with laboratory and field experiments, mathematical models can be used to predict the evolution of a site over long time periods. We present a mass conservative numerical scheme for reactive solute transport in porous media. The transport is modelled by a convection-diffusion-reaction equation, including equilibrium or nonequilibrium sorption. The scheme is based on the mixed finite element method, more precisely the lowest order Raviart-Thomas elements and one-step Euler implicit. The underlying fluid flow is described by the Richards equation, a possible degenerate parabolic equation. A numerical scheme for it was already presented and analysed in [1, 2]. Algorithmic aspects and convergence results will be shown. The scheme is analysed in a general framework, especially taking in account the low regularity of the solution of the Richards equation. An order of convergence estimate is derived which depends on the accuracy of the approximation of the fluid flow. The theoretical results are sustained by relevant numerical studies. REFERENCES [1] F. A. Radu, I. S. Pop, and P. Knabner, Order of convergence estimates for an Euler implicit, mixed finite element discretization of Richards\' equation, SIAM J. Numer. Anal. 42, No. 4, pp.~1452-1478, (2004). [2] F. A. Radu, I. S. Pop, and P. Knabner, Error estimates for a mixed finite element discretization of some degenerate parabolic equations, Numerische Mathematik, to appear, (2008).

top of page

A multicomponent transport and enhanced anaerobic dechlorination model of a single fracture – clay matrix system

Authors:

1. Julie Chambon, Departement of Environmental engineering, Technical University of Denmark (DTU)
2. Philip John Binning, Departement of Environmental engineering, Technical University of Denmark (DTU)
3. Poul Løgstrup Bjerg, Departement of Environmental engineering, Technical University of Denmark (DTU)

TCE is a wide spread subsurface contaminant and an important threat to groundwater quality. Many TCE contaminated sites occur in fractured clay systems, and the remediation of these sites is challenging. At such sites, TCE can flow preferentially along fast pathways, formed by the vertical fracture network, and diffuse into the clay matrix. Counter diffusion of TCE to the fracture can take place for hundreds of years after the removal of the contamination source, causing long-term contamination of an underlying aquifer. Recent laboratory and field experiments have shown that bioremediation may be an attractive method for TCE decontamination. Chlorinated solvents can be anaerobically degraded through sequential reactions to a non toxic end product, ethene. Bioremediation, where an electron donor is injected into the fracture system, is a promising remediation strategy that may be able to reduce clean-up times. In order to evaluate remediation potential, a numerical model was constructed to simulate the transport and degradation of TCE in a fracture – clay matrix system. The model is based on field and experimental work including batch and diffusion experiments. The model is an advance on the state of the art as it includes the most up to date knowledge of the sequential degradation processes and involves full coupling between electron donor, biomass and substrate concentrations. The numerical model couples the transport of fluid in the fracture network with diffusive transport in the clay matrix. The model is constructed using Comsol Multiphysics, a generic finite-element partial differential equation solver, and is composed of a single 1D vertical fracture model coupled with a 2D clay matrix model. Equations for TCE and the biodegradation daughter products are written for both the fracture and matrix. The remediation process is simulated by a constant concentration of electron donor in the fracture. The electron donor diffuses slowly into the clay matrix, in order to enhance biodegradation. The sequential dechlorination process, which occurs in the matrix, is modeled with modified Monod’s kinetics including limiting substrate and competitive inhibition between the chlorinated ethenes. Monod’s kinetics is also used to model growth and decay of biomass in the clay matrix. The numerical model is tested by comparison with batch experimental data for TCE dechlorination and field and experimental data on clay transport and decay. First, site specific results from the batch experiments of reductive TCE dechlorination, in an enrichment culture, using lactate or propionate as electron donor, allow testing the degradation part of the model, including only competitive inhibition and growth and decay of biomass. In a second step, the dechlorination model is tested against field experimental data from a contaminated aquifer. The degradation process was monitored under three different conditions, current site conditions, addition of electron donor (lactate) and addition of both electron donor and dechlorinating biomass. With these data sets, it was possible to test the different components of the degradation model. Finally, the model is tested by comparison with existing approaches, including analytical solutions obtained under simplifying assumptions. The numerical model is employed to simulate the entire history of a contaminated site, including both the contamination phase and subsequent remediation. In the remediation phase, the concentration of TCE and its daughter products is modeled in the clay matrix and at the outlet of the fracture. The amount of substrate, the duration of the injection in the fracture, and the time required to meet the drinking water guidelines at the fracture outlet are assessed. The numerical model allows the formulation and assessment of a complete remediation strategy for TCE contamination in fractured clay based on enhanced anaerobic dechlorination.

top of page

A Multilevel Multiscale Mimetic Method for Two-Phase Flows in Porous Media

Authors:

1. Daniil Svyatskiy, LANL
2. Konstantin Lipnikov, LANL
3. David Moulton, LANL

Flow simulations in porous media involve a wide range of strongly coupled scales. The length scale of short and narrow channels is on the order of millimeters, while the size of a simulation domain may be several kilometers. The permeability of rock formations is highly heterogeneous and may span several orders of magnitude, from nearly impermeable barriers to high-permeable flow channels. For such complex systems fully resolved simulations become computationally intractable. To address this problem we developed the new Multilevel Multiscale Mimetic (M^3) method. Using the same mathematical model with averaged parameters to perform simulations at a much coarser scale does not adequately capture the influence of the fine-scale structure. In contrast, the M^3 method constructs a hierarchical sequence of coarse-scale models, which provides a framework to capture fine-scale effects more accurately. Most of existing model upscaling approaches are based on a two-level structure. Using a two-level structure most multiscale methods achieve a coarsening factor of approximately 10 in each coordinate direction, while the trends in finescale realizations of large reservoirs require a coarsening factor of 100 or more. Using the multilevel hierarchy of the M^3 method we achieve large coarsening factors with small computational cost. The M^3 method provides locally conservative velocity fields on all scales, which guarantees local mass conservation. This is a crucial requirement for modeling two-phase flows. Due to its algebraic nature, the method can be adapted to other fine-scale discretizations, such as the Mixed Finite Element and Finite Volume methods, and can handle full permeability tensors and general polygonal meshes. The conservative coarsening procedure is defined by velocity coarsening parameters. These parameters play a critical role in the accuracy of the M^3 method. We implement a black box, problem-dependent, and computationally inexpensive strategy to estimate them. In most multiscale methods the specific parameters that define the coarsening procedure are computed at the initial time step, with high accuracy, and are not changed during the simulation. Our numerical experiments demonstrate that it is important to update these parameters in time, even with moderate accuracy. Thus, we propose to update our velocity coarsening parameters a few times during the simulation (e.g., every 500 time steps) using an efficient algebraic multigrid algorithm with a modest convergence tolerance. With this update strategy the error in the M^3 solution is comparable to the error in the fine-scale solution. The M^3 method has been applied to the upscaling benchmark from the 10th SPE Comparative Solution project. The permeability field has large channelized structures, which is a challenging problem for multiscale methods. To discretize this system in time we use the IMPES time discretization scheme (IMplicit Pressure and Explicit Saturation). The saturation is updated using the Darcy velocities provided by the pressure solver. The numerical results demonstrate that with a large coarsening factor, such as 30, the M^3 solution is close to the fine-scale solution. In other numerical tests for larger problems we implement more aggressive coarsening with a factor of 64, and also observe good agreement with the finescale solution. The M^3 method speeds up the pressure solver up to 80 times, and the overall simulation 8 times, with respect to the fine-scale simulation.

top of page

A multiscale continuum formulation for reactive transport in hierarchical porous media

Authors:

1. Qinjun Kang, Los Alamos National Laboratory
2. Peter C. Lichtner, Los Alamos National Laboratory
3. Partha P. Mukherjee, Los Alamos National Laboratory

A multiscale continuum formulation for reactive transport in a hierarchical porous medium is developed. The porous medium consists of a primary continuum coupled to multiple secondary continua either in series or parallel. Multiscale phenomena are captured by taking into account isolated pores described by the secondary continua which are coupled to the primary continuum fluid by diffusion. An analytical solution is derived for the stationary state of a single component system with linear reaction kinetics. This solution is shown mathematically to be equivalent to a single-continuum formulation, with an effective surface area that is a function of the various parameters defining the reactive system. As a consequence it is not possible to directly measure experimentally the effective surface area for the system. One implication of this result is that it provides a possible explanation for the observed discrepancy between field and laboratory measured kinetic rate constants which are based on a single-continuum model interpretation. An example is presented on the basis of a multiscale, synthetic, structured porous medium describing transport of a tracer and linear reaction kinetics.

top of page

A multiscale finite element framework for simulating reactive contaminant transport in heterogeneous porous media

Authors:

1. Kalyana Babu Nakshatrala, University of Illinois at Urbana-Champaign
2. Albert Valocchi, University of Illinois at Urbana-Champaign

Mixing of chemical species across plume boundaries has a major influence upon reactive pollutant fate in the subsurface. Small-scale heterogeneity leads to irregular plume boundaries which enhances mixing-controlled reactions through increasing the interfacial area of the plume. Therefore, it is crucial to capture this small-scale heterogeneity in order to properly model reactive transport. Unfortunately, computational limitations do not permit full resolution of the smallest scales of heterogeneity, and thus it is necessary to use a coarse numerical grid, particularly for cases with a large number of reactive species. In order to capture the sub-grid scale heterogeneity effects, we present two multiscale methods for reactive contaminant transport. Multiscale finite element methods have been proven successful for diffusion and Darcy flow problems, and we extend two popular methods to advective-diffusive reactive system. We report our preliminary studies for steady state fast bi-molecular reactions. The first method is based on the multi-scale decomposition of the solution, and is an extension of the numerical subgrid upscaling method NSUM proposed by Arbogast [1] (which has been applied to Darcy flow). The governing system is divided into two sub-problems – coarse-scale and finescale. The concentration of each chemical species is split into two components – coarse-scale (which is defined at the grid scale) and fine-scale (which is defined at the smallest modeled scale). The fine-scale sub-problems are solved locally by constructing numerical Green’s functions, which are independent of the coarse-scale problem. The localization of fine-scale sub-problems is achieved by the closure assumption, which in our case is the prescription of homogeneous Dirichlet boundary condition for the fine-scale sub-problem. The second multiscale method constructs special shape functions that capture the underlying small-scale information. This method is an extension of the multiscale finite element method (MsFEM) proposed by Hou and Wu [2]. We generate basis functions by solving homogeneous advection-diffusion equation with heterogeneous velocity and permeability in each element. With proper choice of boundary conditions, these basis functions form a partition of unity. One uses these special functions as shape functions. Finally, we compare the numerical solutions obtained using the proposed multiscale finite element methods on both model and realistic problems. References [1] T. Arbogast. Numerical subgrid upscaling of two-phase flow in porous media. In Z. Chen, R. E. Ewing, and Z.-C. Shi, editors, Numerical Treatment of Multiphase Flows in Porous Media, volume 552 of Lecture Notes in Physics, pages 35–49, Springer-Verlag, Berlin, 2000. [2] T. Y. Hou and X. H. Wu. A multiscale finite element method for elliptic problems in composite materials and porous media. Journal of Computational Physics, 134:169–1189, 1997.

top of page

A new approach for hydrogeological modeling in discretely-fractured media: from 3D geological representation to numerical simulations

Authors:

1. Blessent Daniela, Laval University, Québec, Canada
2. Therrien René, Laval University, Québec, Canada

A new approach to simulate fluid flow and solute transport in discretely-fractured media has been developed by coupling the HydroGeoSphere numerical code (Therrien et al. 2005) with the mesh generation software LaGriT (Los Alamos National Laboratory) and by representing complex fractured geological formation with the 3D geological modeling platform GOCAD (Mira Geoscience). GOCAD is an advanced software designed to visualize and edit 3D geological objects that are discretized with different geometrical entities, such as points, segments and facets. Since a discrete fracture conceptual model is used, each fracture is explicitly defined by its own geometry and then discretized with 2D finite elements. The surrounding rock matrix, which can be porous and contribute to flow or transport, is discretized with 3D finite elements. For the fractures and matrix, 2D triangles and 3D tetrahedra are chosen because they are the most flexible elements to handle complex geometrical domains that are typical of fractured geological formations. The first step of the approach presented here uses the GOCAD platform to represent fractures as 2D non planar triangulated surfaces in the 3D space. If two surfaces intersect each other, a conforming triangulation at their intersection is needed. In the next step, GOCAD surface files are imported in the LaGriT mesh generation software. To integrate triangulated fracture surfaces within 3D tetrahedral elements representing the rock matrix, a background hexahedral mesh is built. This mesh is basically used to develop a regular point distribution inside the simulation domain and it can be easily refined where required, for example around fractures. All nodes belonging to the hexahedral mesh and to the triangulated fracture surfaces are then copied to a new mesh object and connected to form a 3D tetrahedral mesh. Original fracture surfaces can be extracted from the tetrahedral mesh, together with their connectivity information. The method also ensures that each triangle lying on a fracture surface corresponds to the face of a surrounding tetrahedron. The resulting 3D tetrahedral mesh generated by LaGriT, as well as the triangular fracture meshes, is finally directly read by the HydroGeoSphere code, which performs all calculations with tetrahedral and triangular elements. Illustrative examples have been designed to verify and validate the approach. These examples show that the proposed spatial discretization can reproduce solute propagation along fractures more accurately than other fracture discretization configurations, which tend to lengthen fracture pathways. The reason is that fractures, which can be planar or not, are built and discretized with their actual path length. Simulation results for simple scenarios are also compared to existing analytical solutions. Although complex 3D meshes can be generated, the approach allows for easy refinement around fractures, representing a great asset for the accuracy of the numerical solution. In addition, simulations show that the total CPU time is only slightly increased by using a tetrahedral mesh, making the use of this kind of mesh to represent complex geometries a more suitable choice compared to blocks or prisms based meshes. In conclusion, the examples show that this mesh generation approach and the improved HydroGeoSphere version constitute a new suitable and more realistic way to investigate fractured geological media.

top of page

A new mixed finite element method on hexahedra, its equivalent sparse finite volume formulation, and convergence of multipoint flux approximations in 3D

Authors:

1. Sebastien F. Matringe, Stanford University
2. Ruben Juanes, Massachusetts Institute of Technology
3. Hamdi A. Tchelepi, Stanford University

Multipoint flux approximation (MPFA) methods are finite volume methods that approximate the vector variable (the flux) through interfaces as a linear combination of a scalar variable (the pressure) evaluated at the centers of all the surrounding gridcells. These methods are used to discretize of the second-order elliptic problem on distorted or unstructured grids with full-tensor coefficients. The MPFA O-method [1] is implemented in virtually all petroleum reservoir simulators but a proof of its convergence remains elusive for 3D hexahedral grids. Russell and Wheeler [2] established the equivalence between traditional finite differences and mixed finite elements (MFE) based on the lowest order Raviart–Thomas space (RT0) with quadrature. Arbogast et al. [3] extended this work to tensor coefficients in 3D using an expanded MFE method, that discretizes the scalar, its gradient, and the vector variable. Unfortunately, there is severe loss of accuracy for discontinuous coefficients. Klausen and Winther [4] proved convergence of the MPFA O-method in 2D by means of a broken RT0 space. Recently, Wheeler and Yotov [5] used the first order Brezzi-Douglas-Marini space (BDM1) and a quadrature to prove convergence of MPFA for triangular and quadrilateral elements and 3D tetrahedra. 3D hexahedra were not considered. The difficulty in extending this approach to hexahedra is to obtain a vector space conforming in H(div) and that localizes to a finite volume method under appropriate quadrature. Here, we present a new vector space defined on the reference hexahedron that is conforming in H(div) and that has four degrees of freedom per face [7]. The new space is constructed by enriching the space of linear polynomials by divergence-free polynomials of higher degree. The new velocity space is compatible with a constant pressure space to define a consistent and convergent MFE method. In addition, the MFE method is localized into the MPFA O-method by the application of the trapezoidal quadrature. Thus, we prove, for the first time, convergence of the MPFA O-method on 3D rectangular hexahedra [8]. We obtain optimal error estimates for the new MFE and for the MPFA O-method in both the scalar and vector variables. These results are confirmed by numerical convergence studies. References [1] I. Aavatsmark. An introduction to multipoint flux approximations for quadrilateral grids. Comput. Geosci., 6:405–432, 2002. [2] T. F. Russell and M. F. Wheeler. Finite element and finite difference methods for continuous flows in porous media. In R. E. Ewing, editor, The Mathematics of Reservoir Simulation, pages 35–106. SIAM, Philadelphia, PA, 1983. [3] T. Arbogast, M. F. Wheeler, and I. Yotov. Mixed finite elements for elliptic problems with tensor coefficients as cell-centered finite differences. SIAM J. Numer. Anal., 34:828–852, 1997. [4] R. A. Klausen and R. Winther. Convergence of multipoint flux approximations on quadrilateral grids. Numer. Methods Partial Differential Equations, 22(6):1438–1454, 2006. [5] M. F. Wheeler and I. Yotov. A multipoint flux mixed finite element method. SIAM J. Numer. Anal., 44(5):2082–2106, 2006. [6] S. F. Matringe, R. Juanes, and H. A. Tchelepi. Streamline tracing for general triangular and quadrilateral grids. Soc. Pet. Eng. J., 12(2):217—233, 2007. [7] S. F. Matringe, R. Juanes, and H. A. Tchelepi. Mixed finite element and related control volume discretizations for reservoir simulation on three-dimensional unstructured grids. In SPE Reservoir Simulation Symposium, Houston, TX, February 26–28 2007. (SPE 106117 – Submitted to Soc. Pet. Eng. J.) [8] S. F. Matringe, R. Juanes, and H. A. Tchelepi. A new mixed finite element on hexahedra that reduces to a cell-centered finite difference method. Numer. Math., 2007. (Submitted).

top of page

A Parallel and Adaptive Algorithm for Computation of Shallow Water Flows

Authors:

1. Eli Atetjevich, Ca. Dept of Water Resources
2. Peter O. Schwartz, Lawrence Berkeley
3. Dan Graves, Lawrence Berkeley
4. Julie Percelay, Lawrence Berkeley
5. Qiang Shu, Ca. Dept of Water Resources
6. Phil Colella, Lawrence Berkeley

Water resource management in the Sacramento-San Joaquin Delta requires hydrodynamic and water simulation over large, multi-scale domains with complex geometry. Often problems must be solved over a variety of hydrologic conditions, in which case simulations are long and performance-bound. We are developing a 2D hydrodynamic model for a shallow water estuary that is designed for multi-processor computation using a Cartesian mesh, embedded boundaries (Cartesian cut cells) to conform to bathymetry and the natural boundaries of the system, and adaptive mesh refinement (AMR) for efficient handling of multi-scale problems. The Cartesian solver is based on a second order Godunov algorithm with a Corner Transport Upwind correction of fluxes to achieve better resolution in directions oblique to the grid. Due to its shock-capturing properties, this algorithm has the potential to achieve superior results in transient problems such as flood flows. The stability problems of small cut cells are handled by hybridizing the conservative discretization with a stable, non-conservative discretization at irregular cells, with the result that local accuracy and stability depends only on the dimensions of the larger Cartesian grid cells.We present the formulation of our model, testing results, geometry processing techniques and an application in the Bay Delta. We also describe the institutional setting in which we intend to use and further develop the model and how this affects our design. In the future, our work will be extended to moving boundary problems such as intertidal wetting and drying and flooding, particle transport and three dimensional modeling.

top of page

a parallel-computing method for modeling large-scale groundwater flow

Authors:

1. Dong Yanhui, Institute of Geology and Geophysics Chinese Academy of Sciences, China
2. Guomin Li, hydrogeology Prof. at Institute of Geology and Geophysics Chinese Academy of Sciences, China

Though many advancements and refinements of the computational methods have been achieved over the last decade, large-scale simulation of ground water flow in complicated geological formations remains a challenge. Many problems arise from computational efforts required for solving large sparse-matrix systems. High-performance, parallel-computing techniques, being increasingly used for computationally-intensive scientific applications, are a powerful tool to modeling large-scale ground water flow. In the paper, a parallel solver for MODFLOW2000 is developed for solving large-scale groundwater flow. The parallel solver is based on PCG(Preconditioned Conjugate-Gradient) Package of MODFLOW2000 with OpenMP Application Program Interface (API). And then, PGFSA, a 3D parallel simulator based on a code named GFSA(Li Guomin, 1994) is developed combining Message Passing Interface (MPI) and OpenMP. Simulation results show that the parallel-computing technique enhances modeling capability efficiently in speedup of computing times for large-scale modeling studies. Citation: Arlen W. Harbaugh et. al. 2000, Modflow-2000, the U.S. Geological Survey Modular Ground-Water Model-User Guide to Modularization Concepts and the Ground-Water Flow Process, U.S. Geological Survey Open-File Report 00-92;Michael J. Quinn, Parallel Programming in C with MPI and OpenMP, EISBN: 0-07-282256-2; Mary C. Hill 1990, Preconditioned Conjugate-Gradient 2 (PCG2), a Computer Program for Solving Ground-Water Flow Equations, U.S. Geological Survey Water-Resources Investigations Report 90-4048.

top of page

A path forward for biogeochemical reactive transport modeling of in situ uranium bioremediation

Authors:

1. Steven Yabusaki & Philip Long, Pacific Northwest National Laboratory
2. Yilin Fang, Pacific Northwest National Laboratory
3. Mary Lipton, Pacific Northwest National Laboratory
4. Jill Banfield, University of California Berkeley
5. Robert Hettich, Oak Ridge National Laboratory
6. Nathan VerBerkmoes, Oak Ridge National Laboratory
7. Lee Kerkhof, Rutgers University
8. Derek Lovley, University of Massachusetts Amherst
9. Darrell Chandler, Akonni Biosystems
10. Aaron Peacock, Haley and Aldrich

Ongoing bioremediation research has shown that subsurface uranium immobilization can be achieved in situ by stimulating indigenous metal reducing bacteria with electron donor amendments. Reliable, cost-effective uranium bioremediation in the field will require the development of a systematic and quantitative predictive understanding of the complex linkages between hydrology, geochemistry, and microbiology. Yet, identifying appropriate conceptual process models and unraveling interactions between processes under dynamic and variable hydrologic and geochemical conditions is a significant scientific challenge. We present a development path for integrating microbial processes into a field-scale, multicomponent biogeochemical reactive transport simulation capability that is part of a research project on long-term uranium immobilization at a U.S. DOE, Office of Science Integrated Field Challenge (IFC) Site. The characterization of the microbiological processes at this site is based on gene expression, proteogenomic analysis, stable isotope probing, and phospholipid fatty acid profiles. Under the IFC, focused microarrays and analysis systems are being developed to provide near real-time, in-field mRNA and protein assays to identify temporal and spatial changes in microbial community composition and activity before, during, and after biostimulation. Coupling the growth, metabolism, and function of the dominant microbial populations with concomitant changes in the subsurface biogeochemistry provides a basis for assessing and developing mechanistic process models. Information obtained from microbial assessment techniques and their use in the simulation capability are described below. The cellular concentration of expressed genes over time provides a measure of how individual species are adapting their structure and function to the electron donor amendment and the evolving biogeochemical environment. Gene expression data provides a direct measure of the potential rate of Fe and U reduction that can be compared with U removal rates. Information on up regulation or down regulation of specific genes can provide new knowledge on the conditions controlling the onset and termination of TEAPs as well as microbial function kinetics at different stages of the field bioremediation. Proteogenomic approaches, based on the integration of metagenomic and proteomic data, provide a means to define the strain and species makeup of samples and to track temporal and spatial changes in microbial activity levels. The protein complement can be monitored at the strain level as cells proliferate and TEAP is utilized. Proteogenomics methodologies for relatively complex natural consortia, for which only partly representative isolate genomes are available, are only now being developed, but initial results from the IFC are promising and are expected to yield insights that will directly inform the modeling process. Stable isotope probing will be used to follow labeled carbon in the electron donor amendment through the microbial community over time. This will allow the identification of the principal microorganisms responsible for consuming the electron donor in the initial stages of the biostimulation, quantification of biomass yield for these microorganisms, and, in later stages, identification of predation and/or utilization of decaying biomass by succeeding populations. As our ability to assess the metabolic state of microorganisms in the context of changing biogeochemistry advances, we foresee an evolution to more robust, quantitative linkages to the stoichiometry of biologically mediated reactions and associated rate laws. Comprehensive, large-scale simulation of field-scale experiments can then provide a framework to evaluate the scale-up of the fundamental mechanisms, and identify the conditions where these mechanisms will be important under field conditions.

top of page

A PILOT TEST OF CO2 SEQUESTRATION IN A DEPLETED GAS FIELD IN THE PO PLAIN, ITALY: FLUID-DYNAMICAL AND GEO-MECHANICAL ISSUES

Authors:

1. Massimiliano Ferronato, DMMMSA - University of Padova
2. Giuseppe Gambolati, DMMMSA - University of Padova
3. Carlo Janna, DMMMSA - University of Padova
4. Pietro Teatini, DMMMSA - University of Padova

The 1997 Kyoto protocol prescribes for the agreeing nations a reduction of the carbon dioxide emission to the atmosphere that for Italy should be in prospective as much as 100 Mt/y by 2012. Geological CO2 sequestration in deep geological formations can help fulfill this requirement. The potential to store large CO2 quantity is offered by gas reservoirs (either in production or depleted) and saline aquifers. A pilot project of CO2 sequestration in a depleted gas field located in the sedimentary basin underlying the Po river plain, northern Italy, is under planning. The reservoir is formed by a gentle anticline about 1400 m deep and surrounded by an active aquifer. Sealed by a 200 m thick impermeable formation, the geologic unit of interest is made of several sandy/clayey sequences bounded by regional faults and crossed by a number of local fractures. The field development started in 1950, was practically completed in the mid \'70s, and since the early \'90s the reservoir is used for hydrocarbon storage. The pilot test consists of injecting 8000 t/year of CO2 over 3 years. A numerical study is performed to explore the extent of migration of the injected CO2 plume and investigate the geo-mechanical aspects related to the safety of sequestration, i.e. possible pre-existing fault activation, caprock fracturing, and ground surface uplift induced by changes in the effective stress. The simulations are carried out with the aid of finite element / interface element codes for the geo-mechanical simulations. A number of parametric scenarios are analyzed to address the major uncertainties concerning the location of the injection well, the geomechanical properties of the sedimentary sequence, and the natural stress regime. The results show that the pilot test is quite safe from the geo-mechanical viewpoint.

top of page

A Reactive Transport Model for Lactate Stimulated Hexavalent Chromium Reduction at Hanford 100H Site

Authors:

1. Sumit Mukhopadhyay, Lawrence Berkeley National Laboratory
2. Eric L. Sonnenthal, Lawrence Berkeley National Laboratory
3. Boris A. Faybishenko, Lawrence Berkeley National Laboratory
4. Susan S. Hubbard, Lawrence Berkeley National Laboratory

Lawrence Berkeley National Laboratory and Pacific Northwest Laboratory are collaborating on a project to perform field-scale investigations of biostimulation of hexavalent chromium, Cr(VI), reduction at the DOE Hanford 100H site in Washington. Biostimulation for Cr(VI) reduction with HRC® (Hydrogen Release Compound), produced by Regenesis, is based on stimulating the natural microbial population to induce a dynamic process of biogeochemical transformation in the subsurface. The goal is to reduce Cr(VI) to the trivalent state, which precipitates as chromite [Cr(OH)3(s)], thus removing toxic Cr(VI) from groundwater. HRC is a polylactate polymer bonded that slowly releases lactic acid after injection into an aquifer. HRC thus acts as a long-time electron donor and carbon source for native bacteria that stimulate Cr(VI) reduction to the trivalent state. This process occurs either directly, where Cr(VI) acts as terminal electron acceptor for specific subsurface microbes, or indirectly, where the addition of lactate (via HRC) to the aquifer causes the microbial population to remove oxygen, nitrate, sulfate and other competing electron acceptors and depress the redox potential. Species thus produced also include ferrous iron and sulfide. Changes to the subsurface caused by biostimulation may include such processes as dissolution and precipitation of minerals, surface complexation, generation of gas, changes in soil water content and oxygen concentration, sorption, attachment/detachment, oxidation and reduction, biofilm generation, and changes in permeability and effective porosity. We are developing a reactive transport model to investigate the lactate-induced reduction of Cr(VI) in the groundwater of the Hanford 100H aquifer. The ultimate objective of the modeling activity is to investigate some or all of the biogeochemical processes described above with respect to lactate-induced Cr(VI) reduction in groundwater. In this paper, however, we focus on the impact of certain geochemical reactions (such as goethite reduction to Fe(II) and subsequent reduction of CR(VI) by FE(II); impact of bicarbonate on the extent of CR(VI) reduction, etc.) controlling Cr(VI) reduction in groundwater. The numerical model is based on the TOUGHREACT reactive geochemical transport software [Xu et al., 2006]. The model domain included three-dimensional heterogeneous saturated Hanford sediments. Spatial distribution of hydraulic conductivity was obtained through geophysical measurements [Hubbard et al., 2006]. The grid is extremely refined and radial near the injection and pumping wells, but is rectangular and less refined farther away from these wells. HRC was injected through an injection well, and groundwater was pumped from a monitoring well situated about 5 m away from the injection well. The flow model (i.e., without reactive transport) was constrained against Br tracer injection test data. We compare simulated time-series of Cr(VI) in the dissolved phase with measured Cr(VI) field data. Modeling confirms that FE(II) has a strong impact on the extent of Cr(VI) reduction. In addition, since it is assumed that the model domain is located in the saturated zone, no separate gas phase was present. Consequently, carbon dioxide produced during the reaction was accounted for as dissolved bicarbonates (HCO3-). This resulted in an oversaturation with HCO3- (as no carbon dioxide could leave the model domain), ultimately resulting in an underprediction of Cr(VI) reduction. We are thus currently extending the model to include an unsaturated zone and other biogeochemical processes.

top of page

A Scientific Data Repository Based on Ontology Design

Authors:

1. Amanda Hines, US Army Engineer Research and Development Center
2. Chris Bogen, US Army Engineer Research and Development Center

Over the past few years computer simulations have gained wider acceptance as viable alternatives to field testing. Due to this recognition, the numbers of simulations and types of simulations have increased significantly. This increase in computational modeling has produced the need for access to large amounts of data in a short amount of time. Also, it is not uncommon for data requirements to change at any moment, especially in a research environment such as the US Army Engineer Research and Development Center (ERDC). The ERDC recognized the need for a centralized, scientific data repository. It was essential for this repository to be easily accessible to users and flexible enough to handle the range of data required by the modelers. This presentation will discuss the repository, iSciDat, which has been created to address these issues at the ERDC. At the design stage and throughout the implementation of iSciDat, modelers were contacted to provide input about the process. Initially, these modelers were from the groundwater and surface-water fields. Based on their input, the repository’s initial functionality was to store soil properties. After obtaining a list from the modelers, these properties were generated inside iSciDat. However, due to the evolving nature of research, the list has grown and will continue to grow. The best-equipped people to add or change these properties are the modelers themselves. This presented a dilemma since most modelers are not database experts. Therefore, the repository had to be created with enough flexibility that the modelers could easily manage the data definitions and relationships on their own. With this as the major driving force for design, we determined that an ontology-based repository was the best option. Ontologies provide the flexibility and hierarchical structure needed for this environment. This implementation of the repository allows modelers full control of their data storage without the daily need of a database manager and it also provides accessibility to store and retrieve data through the Internet. This presentation will discuss in depth the design rationale and implementation for the ERDC repository iSciDat.

top of page

A series solution for multi-layer aquifers with natural geometry

Authors:

1. Sanders Wong, Dept. of Civil and Environmental Engineering, University of Waterloo
2. James R. Craig, Dept. of Civil and Environmental Engineering, University of Waterloo

An analytical solution has been developed for a two-dimensional steady-state saturated groundwater flow system with multiple layers of arbitrary geometry. Each layer has a unique hydraulic conductivity value, and there is no limitation upon the number of layers included in the model. Similar to Toth’s solution (Toth, J., A theoretical analysis of groundwater flow in small drainage basins, J. Geophys. Res., 67, p4795-4812, 1963), the system is enclosed on the bottom and sides by no-flow boundary conditions; a known configuration of the water table leads to a specified-head condition along the top of model. The analytical series solution is obtained by maintaining flow and head continuity across each interface by minimizing the sum of squared errors at a set of control points. The specified head condition at the water table is treated similarly. Similar to Read (Read, W.W., Series solutions for Laplace’s equation with nonhomogeneous mixed boundary conditions and irregular boundaries, Math. Comput. Modelling, 17(12), p9-19, 1993), the geometry of each layer can be represented using arbitrary continuous geometric functions, e.g., Fourier series or splines. The method will be demonstrated via application to two-layer and multi-layer aquifer systems with different arbitrary shapes of the water table, bottom boundary, and layer interfaces. The solution, implemented in Matlab, is computationally inexpensive, highly accurate, and, unlike numerical solutions, is easily reconfigured for different geometries and performs well for very low height-to-length ratios.

top of page

A Stochastic-Lagrangian Model for Multiphase Flow in Porous Media: Upscaling of Non-Equilibrium Pore Scale Dynamics

Authors:

1. Manav Tyagi, ETH Zurich
2. Patrick Jenny, ETH Zurich
3. Hamdi A Tchelepi, Stanford University

CO2 storage in subsurface formations involves many complex physical processes that are well understood at the pore scale. However, in a real scenarios, it is not feasible to perform pore scale simulations, and one needs a large scale model. Unfortunately, the Darcy based large scale models for multiphase flow and transport are questionable, particularly in the context of CO2 storage. One of the main assumptions in these models is the concept of relative permeability and capillary pressure, which are expressed as functions of saturation. Moreover, these quantities are typically measured under equilibrium conditions and scale independence is assumed. However, under unstable conditions these assumptions may not be valid. We proposed [1] an alternative modeling approach based on the Lagrangian movement of stochastic particles, which represent infinitesimal fluid phase volumes and evolve such that their statistics represent the statistics of the actual fluid volumes. Each particle is labeled by a fluid phase, has its individual velocity, a position in physical space and possibly other properties, e.g. CO2 concentration. As a particle moves in physical space, its properties change according to the stochastic processes such that the specified Lagrangian statistics is honored. The goal of this work is to model these stochastic processes based on the pore scale physics. For the illustration purpose a simplified 3D pore-network consisting of spherical pores and cylindrical throats is considered. The flow through the network is described by simple rules for the fluid movement in the throats and for the pore filling [2]. First, flow with self-similar average saturation profiles are investigated. These cases correspond to multiphase Darcy flow, where the relative permeabilities are uniquely determined by the saturation values. More interesting are those scenarios, where no self-similar profile can observed and classical multiphase Darcy model fails. For such cases non-equilibrium phenomena at the pore scale play an important role and we demonstrate how these can be modeled in our stochastic particle framework. With the help of pore-network simulations first the non-equilibrium pore scale dynamics has been studied. Then, based on the extracted Lagrangian statistics, comparative simulations with stochastic particle method have been performed. References: 1. M. Tyagi, P. Jenny, I. Lunati, H. A. Tchelepi, 2006, Multiscale method for multiphase flow in porous media using stochastic particles, CMWR XVI , 19-22 June 2006, Copenhagen, Denmark. 2. R. Lenormand, E. Touboul, and C. Zarcone, 1988 Numerical models and experiments on immiscible displacements in porous media, J. Fluid Mech. 189, 165-187.

top of page

A study of vortices in shallow water with BGK algorithm

Authors:

1. Man Yue, Lam, Mphil student, Deparment of Civil Engineering, The Hong Kong University of Science and Technology
2. Hongwei, Liu, Research Assist., Deparment of Civil Engineering, The Hong Kong University of Science and Technology
3. Mohamed S. Ghidaoui, Professor, Deparment of Civil Engineering, The Hong Kong University of Science and Technology

Large scale vortices with vertical axis and which extend from the bed to the water surface are common in shallow water flows. For example, such vortices can be found in the wake of islands, in the downstream of headlands and other protrusions, at the confluence of two streams with differing speeds (e.g., channel junction), in compound and composite channel, etc. This paper conducts numerical experiments in order to investigate the origin of the vortices and their growth in scale as well as their sensitivity to inflow perturbation, Froude number and bottom friction. The numerical model used is based on the finite volume approach where the fluxes are determined from the Boltzmann Bhatnagar-Gross-Krook (BGK) (Ghidaoui et al. 2001). The BGK scheme is second order in both time and space, conserves both mass and momentum, satisfies the entropy condition and has been previously applied to various hydraulic problems including bores, oblique hydraulic jumps, tidal waves and shallow mixing layers. References: Ghidaoui, M. S., Deng, J. Q., Gray, W. G., Xu, K. 2001 A Boltzmann based model for open channel flows. Int. J. Numer. Meth. Fluids, 35, 449-494

top of page

A survey on non-negative numerical formulations for tensorial diffusion equations

Authors:

1. Kalyana Babu Nakshatrala, University of Illinois at Urbana-Champaign
2. Albert Valocchi, University of Illinois at Urbana-Champaign

Robustness of numerical methods for flow and transport problems in porous media is important for development of simulators to be used in a wide range for applications in subsurface hydrology and contaminant transport. Obtaining non-negative solutions on unstructured meshes along with local and global mass balance is a particularly essential feature in the simulation of reactive transport of contaminants since the nonlinear reaction term leads to nonsensical unstable solutions for negative concentrations. The heterogeneity in the velocity field will give rise to a non-homogeneous anisotropic dispersion tensor with non-negligible cross terms. Several studies have shown that standard treatment of the cross-dispersion term will result in non-positive solutions on general computational grids. Several ad-hoc methods have been proposed in the literature. For example, a post processing step is performed in which one does some sort of “smoothing”. But this procedure, in many cases, will violate mass balance. Some other methods are limited in their range of applicability (e.g., the method can handle only structured grids, etc.). Herein we consider only the dispersion/diffusion process for steady single-phase flow in heterogeneous anisotropic porous media. Such a flow can be described by the Poisson’s equation with a tensorial diffusion coefficient, which when written in the mixed form is similar to the Darcy equations. We survey various popular numerical formulations from both the finite element and finite volume literature, and comment on their numerical performance with respect to the following metrics: • what are the restrictions on mesh and medium properties to obtain non-negative solutions, • whether the method exhibits local and global mass balance, • what is the accuracy of the method for smooth problems, and • can the method handle unstructured 2D and 3D grids (not just logically rectangular grids)? An ideal formulation is the one that exhibits local and global mass balance, second-order accuracy for smooth problems, and unconditionally produces non-negative solutions on unstructured meshes. From the finite element literature we consider – the classical single-field Galerkin formulation, Raviart-Thomas element, stabilized mixed formulation based on the variational multiscale, stabilized mixed discontinuous Galerkin formulation. From the finite volume literature we consider – multi-point flux approximations (MPFA) proposed by Aavatsmark and co-workers, the formulation by Herrera and Valocchi [1], and a new nonlinear two-point finite volume method proposed by Lipnikov and co-workers [2]. We provide general guidelines based upon the robustness of each method, the ease of programming, flexibility for incorporating advection and transient conditions, and computational efficiency. References [1] P. Herrera and A. Valocchi. Positive solution of two-dimensional solute transport in heterogeneous aquifers. Ground Water, 44:803–813, 2006. [2] K. Lipnikov, M. Shashkov, D. Svyatskiy, and Y. Vassilevski. Monotone finite volume schemes for diffusion equations on unstructured triangular and shape-regular polygonal meshes. Journal of Computational Physics, 227:492–512, 2007.

top of page

A Thermodynamically-Based Model For Predicting Microbial Growth And Community Composition Coupled To System Geochemistry

Authors:

1. Jonathan D. Istok, Oregon State University

‘Bioimmobilization’ of redox-sensitive metals and radionuclides is being investigated as a way to remediate contaminated groundwater and sediments. In this approach, growth-limiting substrates are added to stimulate the activity of targeted groups of indigenous microorganisms and create conditions favorable for the microbially-mediated precipitation (‘bioimmobilization’) of targeted contaminants. We present an approach for modeling this process that couples thermodynamic descriptions for microbial growth with associated geochemical reactions. A synthetic microbial community is defined as a collection of defined microbial groups; each with a growth equation derived from bioenergetic principles. The growth equations and standard-state free energy yields are appended to a thermodynamic database for geochemical reactions and the combined equations are solved simultaneously to predict the effect of added substrates on microbial biomass, community composition, and system geochemistry. This approach, with a single set of thermodynamic parameters (one for each growth equation), was used to predict the results of laboratory and field bioimmobilization experiments at two geochemically diverse research sites. Predicted effects of ethanol or acetate addition on uranium and technetium solubility, major ion geochemistry, mineralogy, microbial biomass and community composition were in general agreement with experimental observations although the available experimental data precluded rigorous model testing. Model simulations provide insight into the long-standing difficulty in transferring experimental results from the laboratory to the field and from one field site to the next, especially if the form, concentration, or delivery of growth substrate is varied from one experiment to the next. Although originally developed for use in better understanding bioimmobilization of uranium and technetium via reductive precipitation, the modeling approach is potentially useful for exploring the coupling of microbial growth and geochemical reactions in a variety of basic and applied biotechnology research settings.

top of page

A variational multiscale high-resolution method for the simulation of unstable multiphase flow in heterogeneous formations

Authors:

1. Francois-Xavier Dub, Massachusetts Institute of Technology
2. Luis Cueto-Felgueroso, Massachusetts Institute of Technology
3. Ruben Juanes, Massachusetts Institute of Technology

Multiscale phenomena are ubiquitous to flow and transport in porous media. They manifest themselves through at least the following three facets: (1) effective parameters in the governing equations are scale dependent; (2) some features of the flow (especially sharp fronts and boundary layers) cannot be resolved on practical computational grids; and (3) dominant physical processes may be different at different scales. Numerical methods should therefore reflect the multiscale character of the solution. We concentrate on the development of simulation techniques that account for the heterogeneity present in realistic reservoirs, and have the ability to solve for coupled pressure-saturation problems (on coarse grids). We express the governing equations of multiphase flow as a pressure equation (of near-elliptic nature) and a saturation equation (quasi-hyperbolic equation). Both are nonlinear, but only weakly coupled. We propose a variational multiscale (VMS) method for the solution of the pressure equation, which splits the original problem rigorously into a coarse-scale problem and a subgrid-scale problem. A high-order finite volume method is used to solve the saturation equation. The proposed VMS method employs a low-order mixed finite element method at the coarse scale, and a finite volume method at the subgrid scale. We identify a weak compatibility condition that allows for subgrid communication across element interfaces, which turns out to be essential for obtaining high-quality solutions. We also introduce an effective, locally conservative formulation of multiscale sources (wells) that relies on a decomposition of fine-scale source terms into coarse and deviatoric components, and on the definition of multiscale “well” basis functions. The saturation equation is then solved on the fine scale by a high-resolution explicit finite difference method. In the course of the simulation, adaptivity is used to minimize the number of shape functions that need to be updated for the solution of the pressure equation. The proposed approach provides highly accurate solutions to flow scenarios that display viscous fingering and channeling, with challenging permeability fields and concentrated fine-scale source terms.

top of page

Adaptive Kronrod-Patterson integration of highly non-linear element matrices

Authors:

1. Hans Janssen, Department of Civil Engineering, Technical University of Denmark

Current adaptive mesh-refinement approaches – be it r-, h- or p-based – all seek to minimise the deviation between a reference solution and the actual discretised solution. The deviation between reference and discretised solution may arise from two sources: deviation due to the inaccurate discretised representation of the continuous fields of independent variables – the discretisation error –, and deviation due to the inaccurate numerical integration of the element storage and transfer matrix coefficients – the integration error. Most adaptive mesh-refinement approaches resolve mesh regions with too large deviations by addition (h/p-based adaptivity) of discretisation nodes, some by relocation (r-based adaptivity). It is hence implicitly assumed that the discretisation error dominates the deviation from the reference solution. For linear PDE’s, with constant storage and transfer coefficients, this assumption is most likely valid. PDE’s describing variably saturated subsurface or building material moisture transfer are however highly non-linear: the moisture capacity and permeability of building materials might vary some ten orders of magnitude from dryness to saturation. For such PDE’s, the actual discretisation may be sufficient to accurately represent the continuous fields of independent variables, but not to accurately integrate the elements’ storage and transport matrices. While such deficient numerical integration may of course be resolved by local densification of discretisation nodes (mesh-refinement), this comes at the price of extra degrees of freedom in the system and a complex mesh-refinement algorithm. An alternative technique, directly targeting the integration error, is presented here. It equally allows to enhance the numerical fidelity of flow and transport simulations, however with a minimal implementational and computational expense, while moreover being complementary to adaptive mesh-refinement approaches. In the first section, the dominance of the integration error in the accuracy of simulation of moisture transfer in building materials will be exemplified. To that goal, one-dimensional and two-dimensional simulations of moisture absorption and drying are performed, with monotonously coarsening discretisations. It will be verified that the deviation from the reference solution is strongly dependent on the order of the used numerical integration scheme, confirming the dominance of the integration error. As it is unnecessary and uneconomical to apply high-order integration for each and every element in the discretisation, an adaptive strategy will be put forward. To that aim, the recursive Kronrod-Patterson numerical integration scheme is implemented with an a posteriori error criterion on the resulting values of the element’s storage and transport matrices. The investigation will show that such adaptive approach provides the right balance between computational accuracy and expense. In a final section a more pragmatic and computationally inexpensive simplification will be presented, applying an a priori estimation of the required integration order, based on the expected variation of the parameters in the element’s storage and transport matrices. It will be demonstrated that, if correctly calibrated, it yields similar results at an even lower computational cost. While the study here is limited to single-phase moisture transfer in porous building materials and the general finite-element method without mesh-adaptation, the presented principle is equally applicable for other strongly non-linear transport phenomena and in more advanced finite element approaches, as well as in other discretisation schemes making use of numerical integration.

top of page

Adaptive Markov Chain Monte Carlo Sampling and High Performance Computing for Estimating Parameters in High-Resolution Three-Dimensional Flow and Transport Models

Authors:

1. Jasper A. Vrugt, Los Alamos National Laboratory
2. Bruce A. Robinson, Los Alamos National Laboratory
3. Cajo ter Braak, Wageningen University

Markov chain Monte Carlo (MCMC) methods have found widespread use in many fields of study to estimate the average properties of complex systems, and for posterior inference in a Bayesian framework. Existing theory and experiments prove convergence of well constructed MCMC schemes to the appropriate limiting distribution under a variety of different conditions. In practice, however this convergence is often observed to be disturbingly slow. This is frequently caused by an inappropriate selection of the proposal distribution used to generate trial moves in the Markov Chain. I this talk, I will show that significant improvements to the efficiency of MCMC simulation can be made by using a self adaptive Differential Evolution learning strategy within a population based evolutionary framework. This scheme, entitled DiffeRential Evolution Adaptive Metropolis or DREAM, runs multiple different chains simultaneously for global exploration, and automatically tunes the scale and orientation of the proposal distribution during the search. I will demonstrate the advantages of DREAM using several synthetic benchmark problems, including a real-world study involving the estimation of permeability and dispersion parameters in a high resolution three-dimensional flow and particle tracking model using Magnetic Resonance Imaging (MRI) data of a conservative tracer moving through a flow cell. This inverse problem is solved using high performance computing with 100 different computational nodes.

top of page

Adaptive multiscale finite element method for subsurface flow simulation.

Authors:

1. J.M. van Esch, Delft university of technology, Deltares
2. F.B.J. Barends, Delft university of technology, Deltares

Natural geological formations generally show multiscale structural and functional heterogeneity evolving over many orders of magnitude in space and time. In subsurface hydrological simulations the geological model focuses on the structural hierarchy of physical subunits and the flow model addresses the functional hierarchy of the flow process. Flow quantities like pressure, flux and dissipation relate to each other by constitutive relations and subunit parameters like porosity and hydraulic permeability. Hydraulic permeability includes the (steady state) intrinsic permeability and the (time dependent) relative permeability for unsaturated conditions. Laboratory experiments provide measurements of the subunit parameters on the fine scale. The intrinsic permeability is the most dominant parameter and the most heterogeneous parameter affecting the fluid flow on the field scale. In general constitutive relations change for different scales. However, for the multiscale problem under investigation, Darcy’s law is supposed to remain valid on both laboratory scale and field scale. If laboratory measurements are treated stochastically within the geological model, then the structural model of the subsurface should be built on the same scale as these supporting measurements. Fully resolved flow simulations on this scale are intractable and a new multiscale technique has been developed. The newly developed adaptive multiscale finite element method extends the original two-level multiscale finite element method to a hierarchy of scales. The multiscale finite element method captures the fine scale behavior on a coarse mesh by multiscale basis functions. Originally the weights of the multiscale basis functions follow from a local flow simulation. The local flow problems are closed by dimensionally reduced flow computation results. In this form the method forms a class of subdomain decomposition techniques and suitable for parallel computation. Variational multigrid methods implicitly construct multiscale basis functions if certain restrictions on the interpolation operator are satisfied. The set of basis functions is defined recursively from the finest level as multigrid operates on a predefined grid hierarchy. The newly developed method applies variational coarsening. At first a multistep coarsening procedure provides a separation of fine grid unknowns and coarse grid unknowns on triangular and tetrahedral elements and constructs a local intergrid transfer operator. The coefficients of the transfer operator, which resemble the weights of the multiscale basis functions, follow from a dimensionally reduced flow problem. Secondly multiscale homogenization extracts the effective permeability tensor from the multiscale basis functions in a pre-calculation phase. The homogenization procedure resembles pressure-dissipation averaging on local flow problems. The effective permeability is not unique and depends on the local boundary conditions. An upscaling criterion detects parts of the domain where averaged permeability tensors may be used on (static) coarse elements, and subsequently the number of integration points can be reduced. Finally adaptive grid refinement concentrates the computational work on localized phenomena like sinks and in infiltration fronts. On the (dynamically) refined parts downscaling of the permeability takes place. A mesh regularization algorithm eliminates non-conforming elements. The parallel implementation of the adaptive multiscale finite element algorithm solves transient partly saturated subsurface flow problem effectively. This will be illustrated by a case study.

top of page

An Accuracy Evaluation of the Parallel Adaptive Hydrology (ADH) Program for Unsaturated Flow Using Analytical Solutions

Authors:

1. Fred T. Tracy, Engineer Research and Development Center, Information Technology Laboratory
2. Stacy E. Howington, Engineer Research and Development Center, Coastal and Hydraulics Laboratory

The parallel Adaptive Hydrology (ADH) finite element computer program has become the standard for many Engineer Research and Development Center (ERDC) flow simulations. One such application is unsaturated flow. Because real-world problems are so complex, testing codes such as ADH for accuracy is often difficult. This is particularly true when flow is in the vadose zone where Richards’ equation is highly nonlinear. Recently, however, Tracy (2006 and 2007) has derived analytical solutions for a box-shaped flow region that is initially very dry until water is applied to the top of the region. This includes two-dimensional (2-D) and three-dimensional (3-D) solutions for steady-state and transient flow with two specified head equations for the top and either specified head or no-flow boundary conditions for the sides of the soil sample. This gives eight scenarios. The bottom of the soil sample remains very dry. Further, this totally unsaturated flow problem can be made more and more nonlinear and thus stress the nonlinear and linear solvers in ADH as much as needed by simply changing a parameter α. The purpose of this work is to test ADH using the above described analytical solutions for various values of α and determine the error in the computations for certain strategic node points. It has been shown using other computational schemes that the errors grow rapidly with the same grid size and time-step size as the problem becomes more nonlinear, and the nonlinear solver often struggles to converge at all. The following are examples of aspects of accuracy that will be considered: (a) effect of grid size, (b) effect of time-step size, (c) effect of nonlinearity as represented by α, (d) robustness of the Newton nonlinear solver and the bi-conjugate gradient stabilized (BiCGSTAB) linear solver with its different preconditioners, (e) effect of “upwinding” the relative hydraulic conductivity computation as compared with using an average value for the element, and (f) effect of increasing the number of processors in the parallel environment. Comparison of results with other groundwater programs will also be made. Both 2- and 3-D analytical solutions will be used in the test. The 2-D solution in the 3-D ADH code will be tested by doing a vertical cross section with a layer of elements in the third direction having the thickness of the grid spacing Δ. In fact, each hexahedron of grid size Δ × Δ × Δ is divided into six tetrahedral elements. - REFERENCES - Tracy, F. T. (2006). “Clean Two- and Three-Dimensional Analytical Solutions of Richards’ Equation for Testing Numerical Solvers.” Water Resources Research, 42, W08503. - Tracy, F. T. (2007). \"Three-Dimensional Analytical Solutions of Richards’ Equation for a Box-Shaped Soil Sample with Piecewise-Constant Head Boundary Conditions on the Top.” Journal of Hydrology, 336, Issues 3-4, 7 April 2007.

top of page

An Adaptive Multiscale Finite Volume Formulation for Two-phase Flow and Transport in Porous Media

Authors:

1. Hui Zhou, Stanford University
2. Seong Lee, Chevron Energy Technology Company
3. Hamdi Tchelepi, Stanford University

Recent advances in multiscale methods have shown great promise in efficiently simulating multiphase flow in heterogeneous media with high resolution. However, existing multiscale methods only solve the flow equations for pressure in a multiscale approach. Saturation in transport equations is usually solved in fine scale, which takes a large portion of CPU cost especially when flow equations are solved with efficient adaptive multiscale approaches. The main difficulty in multiscale modeling of transport equations lies in the nonlocal memory effect in homogenizing a hyperbolic convection equation. The problem is further complicated by the complex interactions between saturation, flow, and spatial heterogeneity in immiscible displacements. We developed an adaptive multiscale finite-volume (MSFV) method for the transport equation of nonlinear two-phase flow in heterogeneous domains. The objective here is not to develop an upscaled (coarse-scale) operator for the transport equation because of the complex interactions between saturation, flow, and spatial heterogeneity in immiscible displacements. Instead, an adaptive reconstruction strategy is developed. The method can be described using prolongation and restriction operators as in a two-level multigrid scheme. The restriction operator is defined as the volume average of the fine-scale saturations in a coarse block. Three adaptive prolongation operators are defined according to the saturation distribution, in which the physical domain is divided into three regions: (1) Region 1 where the injection fluid has not arrived, (2) Region 2, where steep saturation gradients are present, and (3) Region 3 where saturations change slowly in the wake of advancing front. To identify the transition between the regions, specific norm-based criteria are proposed. In Region 1, the transport equation can be completely skipped, whereas in the Regions 2 the local fine-scale transport equations are solved iteratively on a coarse grid (this is referred as Prolongation Operator I). In Regions 3, we developed two approximate prolongation operators: Prolongation Operator II that reconstructs the fine-scale velocity and is locally conservative on the fine grid, and Prolongation Operator III that interpolates saturation changes to yield a locally conservative scheme, but only on the coarse grid. The proposed adaptive multiscale method has been tested with various models, including systems with strong permeability heterogeneity. The results demonstrate clearly that the multi-scale results with adaptive transport calculations are in excellent agreement with the fine-scale reference solutions. Furthermore, the adaptive scheme for coupled flow transport equations yields solutions that are much more computationally efficient than conventional finite difference methods.

top of page

An assessment of the climate change impact on water availability in India using water balance model

Authors:

1. K. Shadananan Nair, Centre for Earth Research and Environment Management

Impact of global climate change, together with the fast increasing population and fast diminishing resources make providing reliable supplies of water a challenging management issue in many parts of the world. This is especially true for a developing country like India with an exploding population, inadequate finance for implementing adaptation measures and several social issues associated with water. Life of millions, especially rural poor has been traditionally linked to agriculture, the largest consumer of water and to related industries. About 18% of world’s population lives in India, but the share of India’s water resources is only 4% of the global. Though water resources development projects have helped in providing adequate supplies of water even in remote villages and in attaining self-sufficiency in food production, fast increasing population and associated needs pose serious threat to water availability. Overuse, pollution and improper management seriously affect the quantity and quality of water. Global climate models predict a possible change in the rainfall and temperature pattern and hence in the water resources of India in the coming decades. An assessment of the change in water availability becomes necessary to prepare a policy and adaptation strategy. In this study, water balance method introduced by Thornthwaite (1948) and modified by Thornthwaite & Mather (1955) has been used to estimate the change in water availability in different States of India. Climate change scenarios for India accepted by the Ministry and the trends in population have been used to compute the change in per capita availability of water from rainfall. Study reveals that the present national per capita water availability of 2150m3 from rainfall will be drastically reduced to 972m3 in three decades. Even today, some states have no surplus water from precipitation annually, after meeting the requirements for evapotranspiration and soil moisture recharge. In almost all parts of India, seasonal water deficiencies show an increasing trend and seasonal surpluses show a decreasing trend. Increase in rainfall predicted by climate models may not be able to compensate for the reduction in soil moisture due to rise in temperature. This will reduce the water level in surface and groundwater resources, and runoff in rivers. This situation will surely worsen the water disputes the economic condition. More investment will be needed in domestic supplies, irrigation, national waterways and for finding alternatives for cheap hydroelectric power. If the present population growth and urbanization continues, the country will find it hard to meet the requirements in domestic and industrial sector and after all in the agricultural sector that may need additional 200MCM water to feed more than 1.5Million population. Guidelines for an appropriate policy and management strategy have been provided, considering socio-economic conditions and natural and man made environmental changes on water.

top of page

An Integrated Hydrology/Hydraulic and Water Quality Model for Watershed-scale Simulations

Authors:

1. Cheng Wang, Department of Civil and Environmental Engeering, Univ. of Central Florida, Orlando, FL, USA
2. Gour-Tsyh Yeh, Department of Civil and Environmental Engeering, Univ. of Central Florida, Orlando, FL, USA
3. Fan Zhang, Environmental Sciences Division, Oak Ridge National Laboratory, Oak Ridge, TN, USA
4. Guobiao Huang, Sutron Corp., West Palm Beach, FL, USA

This paper presents a first principle, physics-based watershed-scale model of integrated hydrology/hydraulic and water quality transport, WASH123D 3.0. This numerical model comprises of three modules: (1) one-dimensional (1-D) simulation module, which is capable of simulating coupled fluid flow, sediment transport and reactive chemical transport and transformation in river/stream networks; (2) two-dimensional (2-D) simulation module, capable of simulating coupled fluid flow, sediment transport, and reactive chemical transport and transformation in two-dimensional overland flow systems; and (3) three-dimensional (3-D) simulation module, capable of simulating coupled fluid flow and reactive chemical transport and transformation in three-dimensional variably saturated subsurface systems. The Saint Venant equation and its simplified versions (diffusion wave and kinematic wave forms) are employed for surface fluid flow simulations and the modified Richards equation is applied for subsurface flow. These governing equations are solved with several physically and mathematically based numerical options. In the transport simulations, fast/equilibrium reactions are decoupled from slow/kinetic reactions with the decomposition of reaction networks, which enables robust numerical integrations. Kinetic variables are adopted as primary dependent variables rather than biogeochemical species to reduce the number of transport equations and simplify reaction terms. In each time step, hydrologic/hydraulic variables are solved in the flow module, kinetic variables are then solved in the transport module. This is followed by solving the reactive chemical system node by node to yield concentrations of all species. In order to obtain accurate, efficient and robust computations, five numerical options are provided to solve the advective-dispersive transport equations and three coupling strategies are employed to deal with the reactive chemistry. Application examples will be presented to demonstrate the design capability of the model. This model may be of interest for the environmental scientists, engineers and decision makers as a comprehensive assessment tool to reliably predict the fluid flow and sediment and contaminant transport on watershed scales so as to evaluate the efficacy of alternative watershed management and remediation techniques prior to incurring expense in the field.

top of page

An integrated model of G. sulfurreducens for investigating bioremediation of U(VI) in subsurface sediments

Authors:

1. Srinath Garg, University of Toronto
2. Radhakrishnan Mahadevan, University of Toronto
3. Timothy Scheibe, Pacific Northwest National Laboratory
4. Yilin Fang, Pacific Northwest National Laboratory
5. Philip Long, Pacific Northwest National Laboratory
6. Derek Lovley, University of Massachusetts

Metals (e.g., chromium, mercury) and radionuclides (e.g., uranium, strontium) are major subsurface contaminants in many US Department of Energy (DOE) facilities. These elements possess high mobility in subsurface environments, and are long lived due to their long half-life periods, and pose a potential risk to humans and the natural environment. In-situ bioremediation is being studied as a potent and powerful technology for remediation of uranium at DOE wastes. This method focuses on the microbial mediated reduction of uranium from soluble U(VI) valence state to the relatively insoluble U(IV) valence state. Since elemental uranium cannot be degraded into daughter products, current research focuses on making use of this radionuclide as an electron acceptor in microbial respiratory processes. In-situ bioremediation can be implemented by introducing a soluble electron donor (e.g., acetate), thereby stimulating the in-situ activity of dissimilatory metal reducing organisms such as G. sulfurreducens. In order to optimize the subsurface bioremediation strategy for maintaining long-term activity of Geobacter species, a framework integrating the genome-scale metabolic model and transport process was developed and numerical simulations of bacterial U(VI) reduction in subsurface sediments were carried out. This model incorporates the flow of groundwater under anaerobic conditions and the dynamics of growth and respiration of G. sulfurreducens with acetate as the electron donor and iron and uranium as the electron acceptors. In contrast to existing models of bioremediation, the new framework accounts for changes in microbial function (respiratory pathways) in response to local geochemical conditions through application of a constraint-based in silico modeling method. Two uses of the in silico model are tested: 1) incorporation of modified microbial growth yield coefficients based on the in silico model, and 2) variation of reaction rates in a reactive transport model based on in silico modeling of a range of local geochemical conditions. Preliminary results from this integrated model will be presented.

top of page

An iterative multiscale finite-volume method converging to the fine-scale solution

Authors:

1. Ivan Lunati, Laboratory of Environmental Fuid Mechanics and Hydrology, EPF Lausanne, Switzerland
2. Seong H. Lee, Chevron Energy Technology Co., San Ramon, CA, USA

The Multiscale Finite-Volume (MSFV) method has been originally developed to efficiently solve homogenous (zero right hand side) elliptic problems on large domains with highly heterogeneous coefficients. The MSFV algorithm has been employed to solve the pressure problem in a pressure-saturation formulation of multiphase flow in porous media. An auxiliary coarse grid (defining the control volumes) is employed, together with its dual, to solve a coarse-scale pressure problem. A set of local numerical solutions on dual cells (basis functions) is used to interpolate the coarse-grid pressure and obtain an approximate fine-scale pressure distribution. Then the saturation problem is solved by a standard Schwarz overlap solution method. Recently, some modifications of the original algorithm have been proposed that allow modeling physically complex flow, e.g., multiphase flow at the presence of gravity and capillary forces (or even the black-oil flow). In this framework, an additional local solution per dual cell is used to correctly capture the effects of the additional physical mechanisms, which appear on the right hand side of the pressure equation. This local solution (correction function) can be regarded as a particular solution of the localized problem, the superimposition of the basis functions being the general solution. Given the localization assumption the definition of basis and correction functions is straightforward and the splitting can be carried out rigorously. For the case of gravity-driven flow, numerical simulations demonstrate excellent agreement between the MSFV solutions for pressure and saturation and the corresponding fine-scale reference solutions. Note, however, that the accuracy of the MSFV solutions depends on the quality of the localization. To improve the boundary conditions of the local problems, hence, the accuracy of the MSFV solution, we develop an iterative scheme and we discuss its similarities with domain decomposition methods. Note that introducing the correction function is essential to obtain an algorithm converging to the fine-scale solution.

top of page

An openMI-coupled surface-subsurface flow simulation model sys-tem to forecast subsurface flood

Authors:

1. Becker, Bernhard, Institute of Hydraulic Engineering and Water Resources Management, RWTH Aachen University
2. Reuter, Christian, Institute of Hydraulic Engineering and Water Resources Management, RWTH Aachen University
3. Schuettrumpf, Holger, Institute of Hydraulic Engineering and Water Resources Management, RWTH Aachen University

The problem of groundwater head rise in urban regions due to river stage rise and inundation during flood water has been reported in several current studies. Groundwater head rise may affect subsurface buildings like underground rail tunnels and cellars and cause flooding of depressed areas or areas be-hind dikes. For some urban regions, site specific groundwater models and inundation models for surface flow already exist. For the case study, we already have a 2D horizontal groundwater model, which has been developed for water supply purpose. This site-specific model is based on the finite element groundwa-ter flow program FEFLOW (DHI-WASY). Also, a 2D storage cell inundation model being used for estimating the danger of flooding in the urban areas exists. For both models, the river stage forms the main boundary conditions. In order to obtain a forecasting tool for operational purpose in „subsurface flood“ protection, these two models shall be coupled iteratively using the openMI standard (www.openmi.org). The challenge in coupling the two models is to make the programs openMI compliant. Usually this requires changes in the source code: the program has to be reorganized to be compiled as a dynamic linked library (dll), and basic program functions must be made available for openMI. Being a commercial program, the source code of FEFLOW is not available for this survey. But FE-FLOW provides a programming interface which allows the user to implement own code into the flow simulation. We show how we made FEFLOW an openMI compliant using this programming interface. The FEFLOW interface and the openMI classes are linked with c++ remote procedure calls and dot net platform invoke routines. The use of existing models in a subsurface flood forecasting modeling system is appropriate due to the following reasons: firstly, no extra site-specific models have to be set up. This reduces the work load for the development of the modeling system, and verification and validation work only needs to be focused on the coupling procedure. Secondly, the forecasting model always will represent the best updated state, because the model components are maintained separately by the authorities who devel-oped them and use them. Thirdly, a high degree of acceptance is achieved when users of the subsur-face modeling system can work with programs and models they already know. In the end, using the openMI standard provides high flexibility for future expansions.

top of page

Analysis of Mixed Explicit-Implicit Method

Authors:

1. Amir Riaz, Stanford University, Stanford CA
2. Romain de Loubens , Schlumberger-Doll Research Center, Cambridge, MA
3. Hamdi Tchelepi, Stanford University, Stanford CA

Efficiency of explicit time integration schemes is limited substantially when the local CFL number varies by orders of magnitude in different regions of the computational domain. Fully implicit schemes allow larger time steps but suffer from excessive numerical dispersion. Higher order spatial discretization can be combined with a fully implicit scheme to minimize the error, however, this approach is also computationally expensive. We present a mixed explicit-implicit finite volume method that retains the benefits of both schemes by reducing the computational cost of higher order schemes through adaptive selection of a small number of fully implicit cells (Adaptive Implicit Method). We show that a $O(\\Delta x)$ error occurs at the boundaries of implicit-explicit cells and propose a new numerical treatment to remove the error. Extensive 2D calculations for miscible and two-phase Darcy flow are carried out to determine the practical feasibility of our adaptive scheme as well as to demonstrate the improvement in accuracy, that is significant in some cases, over conventional low order schemes.

top of page

Analysis of the interplay between contaminant transport scales, aquifer heterogeneity and health related parameters in a risk-driven approach for site characterization

Authors:

1. Felipe P. J. de Barros, UC Berkeley
2. Yoram Rubin, UC Berkeley
3. Reed Maxwell, LLNL

Understanding the role of model parameters in a human health risk assessment due to groundwater contamination may aid decision makers in defining better allocation of resources. Due to our inability to fully characterize the subsurface, it becomes extremely important to develop data acquisition strategies that will best contribute to uncertainty reduction in the final risk. It is of our interest to identify flow and transport conditions in which parametric uncertainty reduction will lead to better estimates of adverse health effects. Furthermore, being able to understand the worth of information from hydrogeological parameters, when counter-balanced with uncertainty and variability in human physiological response and behavioral parameters, may prove useful in defining characterization efforts. By making use of the available computational tools and through a consistent methodology, we investigate the impact of site characterization in the final risk estimates. Simulations are repeated for several different physical conditions such as plume ergodicity, the scale of the capture zone and levels of aquifer heterogeneity. Results indicate that the value of hydrogeological information is dependent on the interplay between length scales of the plume, dimensions of the capture zone and the scale of the screen in which the water is being sampled. The methodology accounts for uncertainty and variability in the risk-related parameters. Sensitivity of model predictions is also investigated for different physiological dose-response models.

top of page

Analysis of the split operator approach for multispecies reactive transport models

Authors:

1. Babak Shafei, PhD student
2. Philippe Van Cappellen, Professor

Analysis of the split operator approach for multispecies reactive transport models: Mixed kinetic and equilibrium reactions Reactive transport models (RTMs) are powerful tools to estimate the fate and transport of chemical constituents. They can be developed to simulate interactions between fluids and solids in porous media through physical and chemical processes. In most cases these interactions involve multi-component chemical systems which introduce multiple unknown species concentrations and nonlinearity in the system. Therefore, multiple advection-dispersion-reaction equations (ADRE) need to be solved simultaneously. Split operator (OS) methods are in widespread use for the numerical solution of coupled nonlinear reactive transport problems. The advantages of OS algorithms come with the cost of the truncation error. Previously developed error analyses have focused on relatively single reaction systems including, for example, unimolecular or bimolecular kinetic reactions, or nonlinear equilibrium adsorption isotherms. An error analysis for OS applied to realistic multispecies and multidimensional reactive transport models considering mixed kinetic and equilibrium reactions is presented. The derived error expression can be applied to a variety of chemical processes taking place in the subsurface such as aqueous complexation, oxidation-reduction reactions, mineral dissolution-precipitation, acid-base equilibria and sorption and ion exchange. The error analysis is illustrated for multi-component reaction networks with particular attention for the errors on concentrations of chemical species participating in different types of reactions (kinetic, equilibrium or both of them). It is shown that when the error related to the effect of the boundary condition (by assuming homogeneous boundary condition) is neglected, the error in the model domain depends on the transport regime (advection or diffusion). Errors in the OS method are compared to those of the fully implicit solution and the differences are discussed.

top of page

Analytical and Numerical Properties of the Diffusive Wave Approximation of the Shallow Water Equations with Applications to Water Flow in Wetlands.

Authors:

1. Mauricio Santillana, Institute for Computational Engineering and Sciences, Univ. of Texas at Austin.
2. Clint Dawson, Institute for Computational Engineering and Sciences, Univ. of Texas at Austin

In this work we address the quantitative and qualitative aspects of modeling water flow driven mainly by gravitational forces and dominated by shear stress, i.e. under uniform and fully developed turbulent flow conditions. These particular flow regimes are present in marshes, wetlands and many estuaries. The effective equation that has been recently used to simulate these particular flow regimes is a doubly nonlinear and degenerate diffusion equation. This equation is often referred to in the literature as the diffusive wave approximation of the shallow water equations (DSW) and it is obtained by approximating the depth averaged continuity equations by empirical laws such as Manning’s or Chezy’s formulas and then combining the resulting expression with the free surface boundary condition. The DSW equation contains as particular cases the Porous Medium equation and the p-Laplacian. In this study we present key characteristics of the solutions of the initial/boundary value problem arising from this equation along with suitable numerical methods aimed at approximating them. We focus our attention in obtaining physical solutions of the DSW equation using the Galerkin method and present a priori error estimates between the true solution and the Galerkin approximate solutions. We present as well some examples where an unstructured-grid finite element code was implemented to solve the DSW equation and used to simulate water flow on real experimental settings as well as real environments.

top of page

APPLICATION OF 2-D COMPUTER MODEL IN RESERVOIR DEVELOPMENT – A CASE STUDY FROM INDIA

Authors:

1. Dr.S.K.Sharma, Geological Research Institute

[Funding dependent abstract - I, being from India need and request for travel grants to participate]. Geothermal exploration carried out in India, so far, has generated valuable data through extensive surface, earth scientific studies backed-up by exploratory drilling down to the depths of about 500 meters at selected geothermal localities concerning structural, geological, geochemical, hydrological and thermal parameters of geothermal systems. One such site is the Tattapani Geothermal Field situated in Surguj