Hydrology [H]

H33H  MS:Exh Hall B   Wednesday
Numerical Simulations of Flow and Transport in Heterogeneous Subsurface Systems Posters
Presiding: A Tompson, Lawrence Livermore National Laboratory; Y Zhang, University of Iowa

H33H-1709 

A Novel Multiple-Relaxation-Time Lattice-Boltzmann Model for Simulating Anisotropic Advective-Diffusive Transport

* Long, W (Wei.Long@bp.com), Pushing Reservoir Limits Team, BP America, Inc., 501 WestLake Park Blvd., Houston, TX 77079, United States Hilpert, M (markus_hilpert@jhu.edu), Johns Hopkins University, 313 Ames Hall, 3400 N. Charles St., Baltimore, MD 21218, United States

We developed a new Lattice-Boltzmann (LB) model to simulate the two-dimensional advection-anisotropic diffusion equation. We approximate the collision term using a multiple-relaxation-time (MRT) method. We validated our model by comparing our simulations to analytical solutions for the cases of pure diffusion and advective-diffusive transport in both isotropic and anisotropic media. Our MRT LB method is more stable than the traditional Bhatnagar-Gross-Krook (BGK) LB model, which uses a single relaxation time.

H33H-1710 

Numerical Simulations of Density-Driven Flow in Heterogeneous Porous Media

* Rapaka, S (saikiran@jhu.edu), Dept of Mechanical Engineering, Johns Hopkins University, 3400 N Charles St, Baltimore, MD 21218, United States Pawar, R J (rajesh@lanl.gov), EES-6, Los Alamos National Laboratory, Los Alamos, NM 87544, United States Stauffer, P H (stauffer@lanl.gov), EES-6, Los Alamos National Laboratory, Los Alamos, NM 87544, United States Zyvoloski, G (gaz@lanl.gov), EES-6, Los Alamos National Laboratory, Los Alamos, NM 87544, United States Chen, S (syc@jhu.edu), Dept of Mechanical Engineering, Johns Hopkins University, 3400 N Charles St, Baltimore, MD 21218, United States Zhang, D (donzhang@ou.edu), Mewbourne School of Petroleum and Geological Engineering, University of Oklahoma, Norman, OK 73019, United States

In the context of geological sequestration of carbon dioxide, it is well known that convective transport is expected to play a crucial role in accelerating the rate of carbon dioxide dissolution into the brine present in the aquifer. Most previous studies of this convective process have considered a homogeneous porous medium. However, the properties of the aquifer are known to be extremely heterogeneous over all length scales. Previous research on convection in heterogeneous media has suggested that global averaged quantities like Rayleigh number etc. may be inadequate for describing convective transport in such systems. In this talk, I will present some recent results by our group on the process of density-driven convection in heterogeneous porous media. We have used high-resolution numerical simulations with a Monte-Carlo approach to understand the role played by permeability heterogeneity. We show that an averaged global Rayleigh number is more useful than previously believed in predicting convective transport in heterogeneous media. We also discuss the variability of the results with increasing variance of the permeability field.

H33H-1711 

Coupled Transport of Magma- and Mantle-Sourced Heat and Helium in Heterogeneous, Fractured Aquifers

Andrews, J L (andre345@umn.edu), University of Minnesota Department of Geology & Geophysics, 310 Pillsbury Drive SE, Minneapolis, MN 55455, United States * Saar, M O (saar@umn.edu), University of Minnesota Department of Geology & Geophysics, 310 Pillsbury Drive SE, Minneapolis, MN 55455, United States

Coupled transport of water and mantle-sourced heat and helium (He) is simulated in order to gain a better understanding of the patterns of temperature, He concentrations, and He isotope ratios (R=3He/4He) observed in groundwater systems where fault structures impact fluid transport, such as grabens, calderas, and volcanic regions. We consider the effects of implementing temperature- and mass-dependent He diffusion coefficients on He signals, as well as permeability, heterogeneity, buoyancy-driven recirculation, and radiogenic heat and 4He production in a variety of one-, two-, and three-dimensional systems with temperature gradients ranging from ~33-90~°C/km. The results of our investigation have applications to geothermal reservoir analyses and studies using heat and/or helium as natural tracers of groundwater flow. We find that even for permeabilities below 10-15~m2, inclusion of the temperature- and mass- dependence on He diffusion can have a large impact on patterns of He concentration and isotope ratios. Further, due to the large difference in the diffusion coefficients of heat and He (~3~orders of magnitude higher for heat) magma-sourced heat and helium signals can be spatially separated, or decoupled, for low-permeability systems in which He transport is dominated by advection, while heat transport is predominately conductive. Even for systems in which heat and He transport are both dominated by advection, temperature and He patterns may not be perfectly synchronized, again due to the wide difference between heat and He diffusion coefficients; while He patterns may be in-line with the dominant flow patterns, temperature profiles tend to be more diffused. Additionally, since He has a low diffusion coefficient, it is prone to entrapment by low-permeability layers, such as crystalline basements. This entrapment of He allows for high concentrations of He in low-permeability layers, such that these layers act as reservoirs of both mantle and crustal He. The depth to these low-permeability layers, therefore, can have a large impact on observed near-surface He signals. Because of the higher diffusion coefficient of 3He relative to 4He, 3He is able to escape more easily, creating lower R-values in these entrapment zones. Buoyancy-driven recirculation cells within a fracture system can lead to decoupling of heat and He isotope ratios, while leaving heat and He concentration relatively well-correlated. Radiogenic processes can obscure these nuances of magmatic heat and He transport by drastically altering the pattern and values of He signals throughout a groundwater system. Our results emphasize the importance of combining temperature, He concentration, and He isotope ratio data towards interpretation of groundwater flow patterns based on these natural tracers.

H33H-1712 

Stratigraphic Controls on Seawater Intrusion, and Implications for Ground-Water Management, Dominguez Gap Area of Los Angeles, California

* Siade, A J (siade@seas.ucla.edu), U.S. Geological Survey, 4165 Spruance Rd., Suite 200, San Diego, CA 92101, United States Nishikawa, T (tnish@usgs.gov), U.S. Geological Survey, 4165 Spruance Rd., Suite 200, San Diego, CA 92101, United States Reichard, E G (egreich@usgs.gov), U.S. Geological Survey, 4165 Spruance Rd., Suite 200, San Diego, CA 92101, United States

Development of ground water in coastal Los Angeles in the 20th century led to extensive water-level declines and associated seawater intrusion. The U.S. Geological Survey (USGS) has developed a solute-transport model to quantitatively test the hydraulic implications of a sequence-stratigraphic model and to assess the possible effects of alternative management strategies. The transport modeling was conducted using SUTRA, a finite-element, density-dependent, ground-water flow and solute-transport model. The SUTRA configuration for this case is two dimensional and considers flow and transport along an approximate flow line extending from the Pacific Ocean through the Dominguez Gap area of coastal Los Angeles. The lithologic representation is based on a stratigraphic cross section developed by Ponti and others (2007) http://pubs.usgs.gov/of/2007/1013/. The transient-state simulation period is from 1850 to 2004. Trial-and-error model calibration was conducted using the measured water levels and chloride (Cl) concentrations at nine wells along the cross section. The results from the calibrated model indicate that faulting can provide the main pathway for downward transport of seawater by juxtaposing low permeability layers with high permeability layers; prior stratigraphic models for the region did not recognize this fault system. Three 20- year management scenarios were considered: (1) status quo, that is, no change in water-management strategies; (2) installation of a slurry wall; and (3) raising inland water levels through increased injection or decreased pumpage. Scenario 1 resulted in increasing Cl concentrations. Scenario 2 slowed Cl migration; however, this did not reverse seawater intrusion. Scenario 3 reversed seawater intrusion, but there remained Cl in the deeper regions that will be removed only by dilution over time.

H33H-1713 

Simulating advective transport through locally refined grids

* Mehl, S), US Geological Survey, 3215 Marine St, Boulder, CO 80303, United States Dickinson, J), US Geological Survey, 520 N Park Suite 221, Tucson, AZ 85719, United States Hanson, R), US Geological Survey, 4165 Spruance Road, San Diego, CA 92101, United States Hill, M C), US Geological Survey, 3215 Marine St, Boulder, CO 80303, United States

Simulation of advective transport with particle tracking is an effective way to approximate contributing areas associated with pumping wells. Using a local-scale model embedded within a regional model improves representation of the hydraulic response due to pumping, and therefore, better representation of the particle trajectories. A consequence of this finer-scale representation near the pumping well is that particles may follow different trajectories than those produced using an unrefined grid. This further implies that failure to represent fine-scale flow features near pumping wells can cause an underestimation of the variability, and therefore, uncertainty, of particle trajectories tracked backward from the pumping well. This work uses new capabilities for tracking particles through locally refined grids to examine the effects of grid refinement and to investigate sensitivity relationships between hydraulic conductivity and porosity parameters and simulated particle termination location and timing. This information is used to assess why some particle trajectories are more sensitive than others and why hydraulic conductivity has greater control over particle trajectories than porosity. These issues are investigated using the shared-node local grid refinement method of MODFLOW-LGR and the particle tracking code, MODPATH. Due to abrupt changes in grid spacing at the interface between grids, a correction of the flows between the locally refined grids is required to simulate accurate particle trajectories. To determine sensitivity and uncertainty in particle termination locations and timing when using embedded models, a universal sensitivity analysis code, UCODE-2005 is used.

H33H-1714 

Simulation of Two Strategies to Enhance Permeable Reactive Barriers in Heterogeneous Aquifer

* Li, L (lin.li@jsums.edu), Jackson State University, 1400 J.R.Lynch Street, Jackson, MS 39217, United States Benson, C (benson@engr.wisc.edu), University of Wisconsin-Madison, 1415 Engineering Drive, Madison, WI 53706, United States

Ground water flow (MODFLOW) and geochemical reactive transport models (RT3D) were used to assess the effectiveness of two strategies in limiting mineral fouling and its impact on hydraulic behavior of continuous-wall permeable reactive barriers (PRBs) employing granular zero valent iron (ZVI). A geochemical algorithm including kinetic expressions of oxidation-reduction and mineral precipitation-dissolution was developed for RT3D. The two strategies that were evaluated are (i) adding pea gravel equalization zones upgradient and down gradient of the reactive zone and (ii) placement of sacrificial pretreatment zones upgradient of the reactive zone. The PRB locates at a three-dimensional heterogeneous sandy aquifer. The sacrificial pretreatment zone contains mixtures of pea gravel and ZVI. Results of simulations show that installation of pea gravel zones provides a more conductive path for ground water flow through the ZVI, which enhances preferential flow and causes greater porosity reductions and shorter residence time in the PRB. After installation of pea gravel zones, the esidence time decreases which is caused by short travel distances in the ZVI due to short circuit of preferential flow. Sacrificial pretreatment zones can be used to elevate the ground water pH and consume many of the mineral forming ions to form secondary minerals in before the reactive zone is reached. The remaining mineral forming ions that pass into the reactive zone cause less mineral fouling. However, mineral fouling by Fe(OH)2 still occurs, and this mineral is formed regardless of the influent mineral forming ions. Addition of the sacrificial pretreatment zone slightly decreases the initial median residence time. However, the pretreatment zone retains higher residence time after 30 yrs due to less mineral fouling in the pure ZVI zone.

H33H-1715 

3-D Numerical Modeling as a Tool for Managing Mineral Water Extraction from a Complex Groundwater Basin in Italy

* Zanini, A (andrea.zanini@unipr.it), Università degli Studi di Parma, Viale G. P. Usberti 181/a, Parma, PR 43100, Italy Tanda, M (mariagiovanna.tanda@unipr.it), Università degli Studi di Parma, Viale G. P. Usberti 181/a, Parma, PR 43100, Italy

The groundwater in Italy plays an important role as drinking water; in fact it covers about the 30% of the national demand (70% in Northern Italy). The mineral water distribution in Italy is an important business with an increasing demand from abroad countries. The mineral water Companies have a great interest in order to increase the water extraction, but for the delicate and complex geology of the subsoil, where such very high quality waters are contained, a particular attention must be paid in order to avoid an excessive lowering of the groundwater reservoirs or great changes in the groundwater flow directions. A big water Company asked our University to set up a numerical model of the groundwater basin, in order to obtain a useful tool which allows to evaluate the strength of the aquifer and to design new extraction wells. The study area is located along Appennini Mountains and it covers a surface of about 18 km2; the topography ranges from 200 to 600 m a.s.l.. In ancient times only a spring with naturally sparkling water was known in the area, but at present the mineral water is extracted from deep pumping wells. The area is characterized by a very complex geology: the subsoil structure is described by a sequence of layers of silt-clay, marl-clay, travertine and alluvial deposit. Different groundwater layers are present and the one with best quality flows in the travertine layer; the natural flow rate seems to be not subjected to seasonal variations. The water age analysis revealed a very old water which means that the mineral aquifers are not directly connected with the meteoric recharge. The Geologists of the Company suggest that the water supply of the mineral aquifers comes from a carbonated unit located in the deep layers of the mountains bordering the spring area. The valley is crossed by a river that does not present connections to the mineral aquifers. Inside the area there are about 30 pumping wells that extract water at different depths. We built a 3-D numerical model of the study area using a finite difference grid describing the surface with 20000 cells (each cell is 30m × 31m). 17 layers represent with high accuracy the boreholes stratigraphy and the available geological cross-sections of the area. The aquifer is described using 6 materials characterized by different hydraulic parameters. Taking into account the particular morphology and the boundary conditions, only 30% of the grid cells were activated in order to improve a better simulation of the physical domain obtaining the reduction of the computation nodes and the speeding up of the numerical processes. MODFLOW 2000 was used to solve the flow problem. After assigning the boundary conditions as fixed heads downstream and input flow discharge upstream, a first calibration of the model has been carried on in steady state condition on the basis of the piezometric heads collected during a measuring campaign. Only a few hydraulic conductivity values were available, so the calibration consisted in varying the hydraulic conductivity of the materials with the aim at reproducing the measured heads in the monitoring points. A refined second calibration in dynamic condition has been carried on because a large dataset of observations during the wells activity were available. The result of the work was a model that allows the study of the flow directions in the aquifers and the analysis of different scenarios of increasing extraction of mineral water with new drilling of pumping wells.

H33H-1716 

Comprehensive 1D Modelling of Reactive Chemical Transport in Unsaturated Soil

* Wissmeier, L (laurin.wissmeier@epfl.ch), Ecole polytechnique federale de Lausanne, Laboratoire de technologie ecologique, Station 2, Lausanne, CH-1015, Switzerland Barry, D A (andrew.barry@epfl.ch), Ecole polytechnique federale de Lausanne, Laboratoire de technologie ecologique, Station 2, Lausanne, CH-1015, Switzerland

Computer models for simulating environmental processes of water flow, solute transport and geochemical reactions have greatly advanced during recent years. However, there is still demand for the development of programs that a capable of simulating the numerous interactions between physical transport processes and biogeochemical reactions in natural soils. We present a new tool for simulating transient vadose zone flow and solute transport according to the moisture- based form of Richards' equation within the widely used geochemical software PHREEQC. The direct implementation into the geochemical framework provides access to comprehensive geochemical models, giving capabilities beyond existing software for coupled unsaturated flow and reaction. Possible reactions include complex aqueous speciation, cation exchange, equilibrium phase dissolution and precipitation, formation of solid solutions, redox reactions, gas phase exchange, surface adsorption considering electrostatics and kinetic reactions with user-defined rate equations, among others. As a result of the close coupling procedure, the influence of geochemical reactions on water content, e.g., through dissolution or precipitation of water-containing phases, can be investigated. For the solution of the partial differential equations of flow and transport, an explicit finite-difference formulation with a second-order space discretization and first-order time discretization was employed. The use of integrated diffusivities transforms Richards' equation into a simple advection-diffusion equation. Changes in water content and solute concentration were conceptualized as local kinetic reactions of individual elements where changes in moisture content result from fluxes of oxygen and hydrogen across cell boundaries. Reactions and chemical element transport are coupled via sequential two-step operator splitting. The scheme was implemented into PHREEQC without any source code modification such that it can be applied by an experienced user within the existing freely available software. In this presentation, we show results from extensive code verification and demonstrate the unique capabilities of the model for simulating surface sorption to variable charge surface sites including the development of a diffuse double layer as well as dissolution reactions with effects on soil moisture. http://ecol.epfl.ch/research/laurin_wissmeier/

H33H-1717 

Macroscopic properties of fractured porous media

Thovert, J (thovert@lcd.ensma.fr), LCD ENSMA, SP2MI, Futuroscope, 86960, France Mourzenko, V V (mourzenk@lcd.ensma.fr), LCD ENSMA, SP2MI, Futuroscope, 86960, France * Adler, P M (padler@ccr.jussieu.fr), Sisyphe-UPMC, 4 place Jussieu, Paris, 75252, France

The determination of the local fields in fractured porous media is a challenging problem, because of the multiple scales that are involved and of the possible nonlinearity of the governing equations. The purpose of this paper is to provide an overall view of the numerical technique which has been used to solve numerous problems. It is based on a three-dimensional discrete description of the fracture network and of the embedding matrix. Any fracture network geometry, any type of boundary condition, and any distribution of the fracture and matrix properties can be addressed, without simplifying approximations. The first step is to mesh the fracture network as it is by triangles of a controlled size. This meshing by an advancing front technique is done successively for each fracture and the intersections between fractures are taken into account. Then, the space in between the fractures is meshed by tetrahedra by the advancing front technique again. The faces of the tetrahedra which are in contact with fractures, coincide with the corresponding triangles in these fractures. The performances of these meshing codes will be illustrated by a few examples. The second step consists in discretizing the conservation equations by the finite volume technique. Specific properties are given to each fracture such as a surface permeability or a joint rigidity. This general technique has been applied to the basic and most important properties of fracture networks and of fractured porous media (1). These properties are single and two phase flows, wether they are accompagnied or not by dispersion of a solute and mechanical properties possibly coupled with flow. These applications will be briefly illustrated by some examples, including when possible comparison with real data. Ref: (1) P.M. Adler, V.V. Mourzenko, J.-F. Thovert, I. Bogdanov, in Dynamics of fluids and transport in fractured rock, ed. B. Faybishenko, Geophysical Monograph Series, 162, 33, 2005.

H33H-1718 

Comprehensive Model for Enhanced Biodegradation of Chlorinated Solvents in Groundwater

Kouznetsova, I (irina.kouznetsova@ed.ac.uk), University of Edinburgh, John Muir building, The King's Buildings, Edinburgh, EH9 3JL, United Kingdom * Gerhard, J I (j.gerhard@ed.ac.uk), University of Edinburgh, John Muir building, The King's Buildings, Edinburgh, EH9 3JL, United Kingdom Mao, X (maoxiaomin@tsinghua.org.cn), China Agricultural University, College of Water Conservancy and Engineering, Beijing, 100083, China Robinson, C (clare.robinson@epfl.ch), Ecole Polytechnique Federale de Lausanne, Laboratoire de technologie Ecologique, Lausanne, CH-1015, Switzerland Barry, A D (andrew.barry@epfl.ch), Ecole Polytechnique Federale de Lausanne, Laboratoire de technologie Ecologique, Lausanne, CH-1015, Switzerland Harkness, M (harkness@crd.ge.com), GE Global Research, One Research Circle, Niskayuna, NY 12309, United States Mack, E E (elizabeth-erin.mack@usa.dupont.com), DuPont Corporate Remediation Group, Glasgow 300, P.O. Box 6300, Newark, DE 19714- 6300, United States Dworatzek, S (sdworatzek@siremlab.com), SiREM, 130 Research Lane Suite 2, Guelph, ON N1G 5G3, Canada

SABRE (Source Area BioREmediation) is a public/private consortium whose charter is to de-termine if enhanced anaerobic bioremediation can result in effective treatment of chlorinated solvent DNAPL source areas. The focus of this 4-year, $5.7 million research and development project is a field site in the United Kingdom containing TCE DNAPL. A comprehensive numerical model for simulating dehalogenation of chlorinated ethenes has been developed. The model considers the kinetic dissolution of DNAPL and nonaqueous organic amendments, bacterial growth and decay, and the interaction of biological and geochemical reactions that might influence biological activity. The model accounts for inhibitory effects of high chlorin-ated solvent concentrations as well as the link between fermentation and dehalogenation due to dynamic hydrogen concentration (the direct electron donor). In addition to the standard biodegradation pathways, sulphate reduction, mineral dissolution and precipitation kinetics are incorporated. These latter processes influence the soil buffering capacity and thus the net acidity generated. One-dimensional simulations were carried out to reproduce the data from columns packed with site soil and groundwater exhibiting both intermediate (250 mg/L) and near solubility (1100 mg/L) TCE concentrations. The modelling aims were to evaluate the key processes underpinning bioremediation success and provide a tool for investigating field sys-tem sensitivity to site data and design variables. This paper will present the model basis and validation and examine sensitivity to key processes including chlorinated ethene partitioning into soybean oil, sulphate reduction, and geochemical influences such as pH and the role of buffering in highly dechlorinating systems.

H33H-1719 

Evaluation of Negatively Correlated Porosity and Permeability on Chemical Migration Through Glacial Outwash Deposits

Morin, R H (rhmorin@usgs.gov), U. S. Geological Survey, Denver Federal Center, Denver, CO 80225, United States * Wellman, T P (twellman@usgs.gov), U. S. Geological Survey, 3215 Marine Street, Boulder, CO 80303, United States

Near surface geophysical logs and hydraulic measurements recorded at decimeter-scale increments along vertical well bores in the glacial outwash deposits of Cape Cod, Massachusetts reveal a negative relation between total porosity and saturated hydraulic conductivity. This finding is somewhat unexpected since many empirical relations predict a positive relation between these parameters for well-sorted sand and gravel deposits with less than one percent fines. To explain the observed negative correlation, we propose a physically-based paradigm that considers the heterogeneity of effective pathways controlling fluid movement. When pathways are conceptualized as idealized conduits with Poiseuille flow the hydraulic conductivity is proportional to the fourth power of their cross-sectional radius while porosity scales to the second power. This disparity in scaling implies that a region may exhibit both lower porosity and greater hydraulic conductivity than at other locations if some or all of its pathways have sufficiently larger effective radii as to offset the conductance loss due to less pore space. The significance of this finding on solute migration is examined using a suite of groundwater models constrained by measurements of hydraulic conductivity and porosity, and observed bedding geometries reported by the U.S. Geological Survey. It is our hypothesis that the negative correlation between hydraulic conductivity and porosity causes localized effects within individual cross beds that eventually diminish at larger length scales, thereby explaining observed tracer behavior. For each model, the relation between porosity and hydraulic conductivity is evaluated as being negatively correlated, positively correlated, and uncorrelated with constant and variable porosity. The range in predicted transport behavior is shown to reflect the variability that could result from assuming alternate parameter relations, as well as the implications of the observed negative correlation.

H33H-1720 

Lattice-Boltzmann Simulations of Heterogeneous Porous Media: A Comparison Between Three Approaches

* Walsh, S D (sdcwalsh@umn.edu), Department of Geology and Geophysics, University of Minnesota, Pilsbury Hall, Minneapolis, MN 55455, Burwinkle, H (hburwin@CLEMSON.EDU), Department of Mathematics and Statistics, Clemson University, O-110 Martin Hall Box 340975, Clemson, SC 29634, O'Grady, R (ogra0014@umn.edu), Department of Geology and Geophysics, University of Minnesota, Pilsbury Hall, Minneapolis, MN 55455, Saar, M O (saar@umn.edu), Department of Geology and Geophysics, University of Minnesota, Pilsbury Hall, Minneapolis, MN 55455,

Lattice-Boltzmann simulations provide a method for modeling fluid flow through complex pore spaces, such as those commonly encountered in many geofluid applications. However, simulating this level of detail comes at a cost -- these explicit, high resolution models often require large amounts of computational resources. An alternative approach, suggested by Dardis and McCloskey (1998), is to introduce a meso-scale Lattice- Boltzmann model that incorporates the porosity of the medium as a model parameter. Rather than a lattice comprising nodes that are either solid or fluid, this approach uses a probabilistic or partial bounce-back model where node properties are varied to reflect the local permeability of the material. These types of models have great potential in a range of geofluids simulations; examples of particular interest to the Geofluids group at the University of Minnesota include the study of pumice permeability, karst formation and well pumping tests. However, there are several different methods for formulating these partial-bounce-back models and little has been done to examine their validity. This presentation compares the predictions of different partial-bounce-back lattice-Boltzmann models under conditions from laminar to turbulent flow. In particular, three models are examined: Dardis and McCloskey (1998), Sukop and Thorne (2004) and an in-house model developed by the Geofluids group at the University of Minnesota. The models' predictions are compared with results from finite element simulations, data from well efficiency experiments and permeability measurements.

H33H-1721 

Characterization of Preferential Flowpaths at the T-Tunnel Complex, Rainier Mesa, Nevada

* Reeves, D M (mreeves@dri.edu), Desert Research Institute, Division of Hydrologic Sciences/Graduate Program of Hydrologic Sciences, 2215 Raggio Parkway, Reno, NV 89512, United States Schultz, R (schultz@mines.unr.edu), University of Nevada, Reno, Department of Geological Sciences and Engineering, 1664 N. Virginia, Reno, NV 89557, United States Bingham, C (kb@umn.edu), University of Minnesota, School of Statistics, 224 Church SE, Minneapolis, NV 55455, United States Pohlmann, K (Karl.Pohlmann@dri.edu), Desert Research Institute, Division of Hydrologic Sciences, 755 E. Flamingo, Las Vegas, NV 89119, United States Russell, C (Chuck.Russell@dri.edu), Desert Research Institute, Division of Hydrologic Sciences, 755 E. Flamingo, Las Vegas, NV 89119, United States Chapman, J (jenny.chapman@dri.edu), Desert Research Institute, Division of Hydrologic Sciences, 755 E. Flamingo, Las Vegas, NV 89119, United States

Rainier Mesa (RM), a tuffaceous plateau on the Nevada Test Site, has been the location of numerous subsurface nuclear tests. The tests were conducted in a series of tunnel complexes located approximately 450 m below the top of the mesa and 1000 m above the regional ground water flow system. The tunnels were constructed near the middle of a 690 m sequence of low-permeability bedded and non-welded vitric and zeolitized tuff units. Though these tuff units are nearly saturated, active ground water flow is restricted to perched lenses occurring within poorly connected normal faults linked to recharge pathways. The perched systems vary from 100 to 150 m above the tunnel complexes which suggests that the now-sealed tunnels could enhance connectivity between otherwise isolated preferential flow pathways. This work represents the first stage of radionuclide transport investigations for the T-tunnel complex at RM: the characterization and stochastic generation of preferential pathways based on fault and joint data collected along tunnel transects. Analysis of fault and joint orientations demonstrate that fractures at RM are steeply-dipping and trend approximately NE-SW. A higher degree of spread about the mean strike relative to the mean dip direction results in an elliptical distribution of fracture orientation. A novel method based on a bivariate normal is used as an alternative to the Fisher distribution to generate fracture orientations. The spatial distribution of fractures along transects indicates fractal clustering (D=0.3) with a power-law distribution of fracture spacing (α=1.1). Fault displacement data are used in conjunction with mechanical models of fault growth to infer both the length and vertical extent of large faults. Major faults are deterministic features in the model domain, while a cascade process is used to govern the spatial distribution of background fractures. Assessment of methods for assigning hydraulic conductivity values to the joints and faults are underway.

H33H-1722 

Development of a Sitewide Groundwater Flow and Radionuclide Transport Model at Idaho National Laboratory

* Huang, H (Hai.Huang@inl.gov), Idaho National Laboratory, P.O. Box 1625, MS 2025, Idaho Falls, ID 83415, United States Magnuson, S (Swen.Magnuson@inl.gov), Idaho National Laboratory, P.O. Box 1625, MS 2025, Idaho Falls, ID 83415, United States Podgorney, R (Robert.Podgorney@inl.gov), Idaho National Laboratory, P.O. Box 1625, MS 2025, Idaho Falls, ID 83415, United States Sondrup, J (Jeffery.Sondrup@inl.gov), Idaho National Laboratory, P.O. Box 1625, MS 2025, Idaho Falls, ID 83415, United States Wood, T (Thomas.Wood@inl.gov), Idaho National Laboratory, P.O. Box 1625, MS 2025, Idaho Falls, ID 83415, United States

A three-dimensional (3D) groundwater flow and radionuclide transport model with a domain of 6400 square kilometers was developed for the Idaho National Laboratory (INL). The model provides a comprehensive evaluation of environmental impacts from operations at the INL on the underlying Snake River Plain Aquifer. The aquifer consists of highly heterogeneous fractured basalt rocks and discontinuous sedimentary interbeds. A multi-objective automated inverse simulation was applied to simultaneously calibrate the model to measured heads, groundwater velocities estimated from isotope studies and derived total water fluxes across model boundaries. A total of 225 wells were used during the calibration of the flow model. A pilot point approach was adopted in the inverse simulations in order to model spatial heterogeneity in the aquifer. This approach resulted in a large inverse problem with more than 350 model parameters to be estimated. A regularization procedure was used to reduce the non-uniqueness of the inverse problem. The transport simulation of tritium transport at the INL over the past 56 years was found to be in agreement with groundwater monitoring data.

H33H-1723 

A Comparison Study of Numerical Simulations in a Borehole Heat Exchanger Field

* Kim, S (kskinc@hanmail.net), Seoul National University, Seoul National University, Seoul, 151-747, Korea, Republic of Bae, G (gokbae@snu.ac.kr), Seoul National University, Seoul National University, Seoul, 151-747, Korea, Republic of Lee, K (kklee@snu.ac.kr), Seoul National University, Seoul National University, Seoul, 151-747, Korea, Republic of Shim, B (boshim@kigam.re.kr), Korea Institute of Geoscience and Mineral Resources, Korea Institute of Geoscience and Mineral Resources, Daejeon, 305-350, Korea, Republic of Song, Y (song@kigam.re.kr), Korea Institute of Geoscience and Mineral Resources, Korea Institute of Geoscience and Mineral Resources, Daejeon, 305-350, Korea, Republic of

Borehole heat exchanger (BHE) field with heat-pump was installed at the building of the Korea Earthquake Research Center in Korea Institute of Geoscience and Mineral (KIGAM). It consists of 28 BHEs equipped with a double U-tube and three monitoring wells. Two BHEs are equipped with optical fiber that measures temperature of the U-tube. A borehole televiewer survey was conducted in four BHEs. The spacing between BHEs is 7 m and the depth of the BHE is 200 m. This geothermal heat-pump system has been operated and monitored since June 2006. We develop a numerical model for BHE field, which can simulate temperature changes in the system with circulating water through the pipe as well as groundwater flow and aquifer temperature changes. This model is based on TOUGH2, a widely accepted three-dimensional numerical simulator for heat and water flow in geothermal systems and verified with the KIGAM data set. Founded on the KIGAM data set, not every BHE need to be operated for air-conditioning of this building at the same time. The spacing in a BHE field is a critical factor for heat exchange efficiency because of thermal interference between BHEs. When some of the BHEs are used, heat exchange efficiency is changed with the selection of BHE array. Using verified model, we find optimal BHE array from six selected cases of the BHE array.

H33H-1724 

Study of Uranium Transport Utilizing Reactive Numerical Modeling and Experimental Data from Heterogeneous Intermediate-Scale Tanks

* Rodriguez, D (drrodrig@mines.edu), Colorado School of Mines, 1500 Illinois St ESE Division, Golden, CO 80401, United States Miller, A (amiller@mines.edu), Colorado School of Mines, 1500 Illinois St ESE Division, Golden, CO 80401, United States Honeyman, B (bhoneyma@minesedu), Colorado School of Mines, 1500 Illinois St ESE Division, Golden, CO 80401, United States

The study of the transport of contaminants in groundwater is critical in order to mitigate risks to downstream receptors from sites where past releases of these contaminants has resulted in the degradation of the water quality of the underlying aquifer. In most cases, the fate and transport of these contaminants occurs in a chemically and physically heterogeneous environment; thereby making the prediction of the ultimate fate of these contaminants difficult. In order to better understand the fundamental processes that have the greatest effect on the transport of these contaminants, careful laboratory study must be completed in a controlled environment. Once the experimental data has been generated, the validation of numerical models may then be achieved. Questions on the management of contaminated sites may center on the long-term release (e.g., desorption, dissolution) behavior of contaminated geomedia. Data on the release of contaminants is often derived from bench-scale experiments or, in rare cases, through field-scale experiments. A central question, however, is how molecular-scale processes (e.g., bond breaking) are expressed at the macroscale. This presentation describes part of a collaborative study between the Colorado School of Mines, the USGS and Lawrence Berkeley National Lab on upscaling pore-scale processes to understanding field-scale observations. In the work described here, two experiments were conducted in two intermediate-scale tanks (2.44 m x 1.22 m x 7.6 cm and 2.44 m x 0.61 m x 7.6 cm) to generate data to quantify the processes of uranium dissolution and transport in fully saturated conditions, and to evaluate the ability of two reactive transport models to capture the relevant processes and predict U behavior at the intermediate scale. Each tank was designed so that spatial samples could be collected from the side of the tank, as well as samples from the effluent end of the tank. The larger tank was packed with a less than 2mm fraction of a composite field material collected from Naturita, Colorado, a Uranium Mill Tailings Remedial Action (UMTRA) Site. The smaller tank was heterogeneously packed into varying layers representing a subdivision of the less than 2mm fraction into two fractions consisting of 0 to 0.250 mm and 0.250 mm to 2 mm. Various physical and chemical parameters were measured in each tank. This paper presents the results from these tank studies as they pertain to a model comparison analysis that was completed. Reactive transport simulations were carried out with the code CrunchFlow and compared with the United States Geologic Service code RATEQ. Both codes were developed to simulate reactive transport, although their code structures are different. RATEQ utilizes MODFLOW to simulate groundwater flow and the framework of MT3DMS to incorporate a reactive transport module. CrunchFlow is self-contained in that the flow and transport portions of the code are solved within the same multicomponent model. The model analysis demonstrated that the incorporation of a kinetic module into RATEQ had the ability to capture the non-equilibrium behavior of the uranium migration in the system as was observed in both intermediate-scale tank experiments. http://laer.mines.edu/projects/uranium.html

H33H-1725 

How does a priori estimated spatial structure of the K field affect hydraulic tomography results?

* Huang, S (syh1019@ntu.edu.tw), Graduate School of Engineering Science and Technology - Doctoral Program,National Yunlin University of Science & Technology, 123, Section 3, University Road, Yunlin, 64002, Taiwan Wen, J (wenjc@yuntech.edu.tw), Graduate School of Engineering Science and Technology - Doctoral Program,National Yunlin University of Science & Technology, 123, Section 3, University Road, Yunlin, 64002, Taiwan Yeh, T J (ybiem@mac.hwr.arizona.edu), Department of Hydrology and Water Resources, The University of Arizona, John Harshbarger Building 1133 E. North Campus Drive, Tucson, Arizona, 85721, United States

Hydraulic tomography (HT) has been widely studied for characterizing the heterogeneity of hydraulic parameters in the subsurface. Although previous researches showed very promising results of hydraulic tomography for identifying hydraulic properties, several issues are still wide open. It is still plagued by a reasonable estimate of the hydraulic conductivity (K) distribution. This article aims at exploring how a priori estimated spatial structure of the K field affects the final result, and assessing the insight of HT. Several synthetic examples with the tomography data are obtained by solving numerically the transient flow equation in a know conductivity field. The sequential successive linear estimator (SSLE) approach is applied to interpret data from the transient hydraulic tomography to an estimated unknown K and specific storage (S) fields of aquifers. Through these numerical examples, the paper demonstrates and answers this question. Keywords: hydraulic tomography, inverse problem, sequential successive linear estimator

H33H-1726 

Three-Dimensional Modeling of Permeability and Porosity Reductions in Saturated Porous Media

Malaguerra, F (flavio.malaguerra@epfl.ch), Ecole polytechnique federale de Lausanne, Laboratoire de technologie ecologique, Institut des sciences et technologies de l'environnement, Lausanne, CH-1015, Switzerland Brovelli, A (alessandro.brovelli@epfl.ch), Ecole polytechnique federale de Lausanne, Laboratoire de technologie ecologique, Institut des sciences et technologies de l'environnement, Lausanne, CH-1015, Switzerland * Barry, D A (andrew.barry@epfl.ch), Ecole polytechnique federale de Lausanne, Laboratoire de technologie ecologique, Institut des sciences et technologies de l'environnement, Lausanne, CH-1015, Switzerland

Changes in hydraulic properties of soils and aquifers as a result of biogeochemical transformations, such as bacteria growth and mineral phase precipitation/dissolution, may lead to significant modifications of the groundwater flow field. This affects in turn the migration pathways and transport rates of solutes, and consequently their spatial distribution. The aim of this work is to investigate different aspects of the field-scale evolution of porosity and hydraulic conductivity in saturated porous media due to bacteria development in the pore space. As a part of this study, a new module was developed for PHWAT to add the capability of modeling clogging. PHWAT is a general flow and multi-component reactive transport computer code based on SEAWAT and PHREEQC-2. The new model incorporates several constitutive equations to convert porosity changes to hydraulic conductivity, as well as biomass attachment/detachment dependent on pore water velocity. Spatial distributions of simulated porosity and hydraulic conductivity changes were compared both against published laboratory data and previous modeling results. We concluded that the model is able to reproduce the clogging process in a reasonably accurate way. Nevertheless, the choice of the constitutive equations and selection of their parameters is problematic. We observed that a single relationship may not be suitable to capture the observed behavior as the hydraulic conductivity decreases. Further research needs to be devoted to understand the pore-scale processes contributing to permeability changes. Following model validation, a synthetic yet realistic contamination scenario was set up to study how porosity and permeability reductions induced by microbial oxidation of contaminants interact with existing geological heterogeneities. We observed strongly nonlinear behavior, resulting in the development of complex spatial contaminant distributions. The original heterogeneous distribution of hydraulic conductivity is modified by the clogging, with the creation of preferential paths that may enhance contaminant spread.

H33H-1727 

Random-walk Particle Tracking for Transport Simulation: Extension to Multi-Component Reactions.

* Scheibe, T D (tim.scheibe@pnl.gov), Pacific Northwest National Laboratory, PO Box 999 MS K9-36, Richland, WA 99354, United States Fang, Y (yilin.fang@pnl.gov), Pacific Northwest National Laboratory, PO Box 999 MS K9-36, Richland, WA 99354, United States Tartakovsky, A M (alexandre.tartakovsky@pnl.gov), Pacific Northwest National Laboratory, PO Box 999 MS K9-36, Richland, WA 99354, United States Redden, G D (george.redden@inl.gov), Idaho National Laboratory, PO Box 1625 MX 2208, Idaho Falls, ID 83404, United States

Random-walk particle tracking methods have been used extensively to simulate transport of conservative solutes in groundwater systems, and are particularly advantageous in cases of heterogeneous flow and advection- dominated transport. However, such models have not been broadly applied to reactive transport simulations, mostly because reactions are formulated in terms of concentrations and conversion back and forth between particle density and concentration leads to inefficiencies and numerical errors. We present an extension of particle tracking methods to include multiple mobile reacting species, in which reactions are formulated in terms of state transition probabilities rather than rates of change of concentration. This approach has previously been successfully applied to cases where there is only one mobile reactant (e.g., sorption/desorption or radioactive decay reactions), but not to multiple interacting mobile species. This particle-tracking method is particularly well- suited to simulation of reactive transport scenarios in which concentrations are poorly-defined because of strong concentration gradients at scales smaller than the REV. As an illustrative example, we apply our approach to an intermediate-scale experiment conducted in a quasi-two-dimensional flow cell in which calcium carbonate precipitation occurs along a narrow interface between two solutions.

H33H-1728 

Pore Scale Simulation of Multiphase Flow using Finite Element Finite Volume Unstructured Mesh Discretization

* Akanji, L (l.akanji06@imperial.ac.uk), Institute of Petroleum Studies, Department of Earth Science and Engineering, Imperial College, Exhibition Road, London, SW7 2AZ, Matthai, S), Institute of Petroleum Studies, Department of Earth Science and Engineering, Imperial College, Exhibition Road, London, SW7 2AZ,

ABSTRACT A full knowledge of the relationship between multiphase flow and pore geometry properties at pore scale will allow a detailed and accurate description of fluid flow on the larger scale using appropriate partial differential equations. However, these constitutive relationships; which are a direct consequence of the complicated geometry of the pore space, are not usually derived from the detailed representation of the pore space but from experiment. The intent of this article is to describe a first principle based numerical simulation method for deriving constitutive relationships governing fluid flow in porous media. The methodology is based on computational fluid dynamics using a finite-element finite-volume (FEMFV) unstructured mesh discretization of the pore geometry; finely resolving individual pores. Steady – state, viscous, laminar flow simulations - using the Reynolds lubrication equation; a simplified form of Navier – Stokes equation for Newtonian fluids, were carried out and benchmarked against the analytical solution. The benchmarked model was subsequently used to compute the integral properties – saturation, capillary pressure, relative permeabilities, pore velocity etc - of multiphase flow. Monitoring of the capillary pressure – saturation relationship within the pores reveals a unique integral curve which can be rationalized in terms of bubbles emerging from the pore throats to the pores and the velocity pattern within the pores shows a stair-step parabolic profile. The relative permeability curves that correspond to these flows are a by-product of the simulations. Expectedly, the porous medium follows a Brooks Corey pattern. Key words: capillary pressure, saturation, pore velocity, relative permeability, channel flow, numerical simulation

H33H-1729 

Investigation of coupled heat and mass transfer in heterogeneous porous media using numerical simulations

Illangasekare, T H (tissa@mines.edu), Environmental Science and Engineering, Colorado School of Mines, 1500 Illinois Street, Golden, CO 80401, United States * Frippiat, C C (cfrippia@mines.edu), Environmental Science and Engineering, Colorado School of Mines, 1500 Illinois Street, Golden, CO 80401, United States * Frippiat, C C (cfrippia@mines.edu), Dept. of Civil and Environmental Engineering, Universite catholique de Louvain, Place du Levant, 1, Louvain-la-Neuve, B-1348, Belgium Zyvoloski, G A (gaz@lanl.gov), Earth and Environmental Sciences, Los Alamos National Laboratory, P.O. box 1663, Los Alamos, NM 87545, United States

A significant body of knowledge exists on separates processes of thermal and mass transport in granular and fractured subsurface formations. However, the need to simulate these processes in a fully coupled way has become necessary to deal with problems associated with long-term-storage of nuclear waste, and the development of new technologies for subsurface remediation. Another emerging area for research is associated with the development of technologies for in situ extraction of underground resources. Numerical models that couple thermal and mass transport processes will play a crucial role in understanding the fundamental processes associated with these new technologies, as well as in making predictions on how complex subsurface systems are expected to behave. It is our hypothesis that heat transport will have a significant impact on distributions of solute concentration, through temperature-dependent dissolution and precipitation, and temperature-dependent rate-limited diffusive transfer of solutes in fractured or highly heterogeneous media. A number of issues related to the validity of existing numerical tools that capture these processes, and their application to field systems through up-scaling need to be investigated. With this overall goal in mind, in this preliminary study, we explore the effect of the variability of subsurface properties on heat and mass transport using simulations conducted using an existing multiphase model. The finite-element code FEHM (Finite-Element Heat and Mass transport code) used in this study was developed at Los Alamos National Laboratory. This code allows for the coupled simulation of flow, heat and mass transport, accounting for density effects and dissolution and/or precipitation reactions. Our analysis is based on two- and three-dimensional simulations using synthetic data sets. Heterogeneous facies distributions are generated according to Markov Chain transition probability models. A distributed source of constant temperature and concentration is located at the inlet model boundary, and the average concentration at the outlet boundary was evaluated. The results that investigated the influence of facies correlation lengths, hydraulic conductivity contrasts, and thermal conductivity contrasts are presented.

H33H-1730 

Lattice Boltzmann Method for Heterogeneous and Anisotropic Advection-Dispersion Equation in Porous Medium Flow

* Servan-Camas, B (bserva1@lsu.edu), Louisiana State University, Department of Civil and Environmental Engineering, 3418G Patrick F. Taylor Hall, Baton Rouge, LA 70803, Tsai, F T (ftsai@lsu.edu), Louisiana State University, Department of Civil and Environmental Engineering, 3418G Patrick F. Taylor Hall, Baton Rouge, LA 70803,

The objective of this study is to solve the heterogeneous, anisotropic advection-dispersion equation (ADE) in porous media using lattice Boltzmann method (LBM). The recovery of anisotropic hydrodynamic dispersion using LBM has not been extensively studied, and using LBM to solve the ADE remains a challenge. In this work, we develop a multidirectional-speed-of-sound (MDSS) method, which are direction-dependent, to directly link the second moment of equilibrium distribution functions (EDFs) with different components in the dispersion tensor. Using MDSS and second-order EDFs, we recovered the ADE with additional numerical diffusion terms, which are not negligible in the general case. We introduced pseudo-velocities to reduce the numerical diffusion terms one order lower and keep the LBM results at second-order accuracy. Two-dimensional mass instantaneous release problems in uniform flow are solved to demonstrate the capability of using MDSS to recover the anisotropic hydrodynamic dispersion. Good agreement is found between the analytical and LBM solutions for high ratio of longitudinal to transverse dispersivities.

H33H-1731 

Numerical Simulation of Biodegradation Processes in Subsurface Systems

* Hernandez-Rendon, M (carmenhr@geofisica.unam.mx), Instituto de Geofisica, Circuito de la Investigacion S/N, Ciudad Universitaria, Mexico, D.F., D.F 04510, Mexico

In this work an evaluation of a numerical scheme to simulate mutispecies reactive transport undergoing biodegradation in porous media is presented. This procedure relies on the application of an operator splitting algorithm introduced by Glowinski [1]. It was used to simulate a system of two species transport with non linear reactions with good results. The main advantage in applying it is that the system of partial differential equations is fully decoupled. Approximation of time derivative is obtained using standard procedures. The spatial advection- diffusion part is solved by means of the standard Galerkin finite element method; nonlinear integrals are evaluated with Gaussian cuadrature. Numerical results are presented for two test cases that take into account the transport of five species; the first one with first order and sequential biotransformation is compared with the analytical solution and with the numerical approximation of the coupled system. In the second, nonlinear and simultaneous reactions are included. [1] Glowinski, R., (2000). Operator-splitting methods for initial value problems: Application to the atmospheric dispersion equations, Lectures 6-8, University of Houston.

H33H-1732 

Numerical investigation of NAPL Source Zone Architecture in Two-Dimensional and Three- Dimensional Unsaturated Porous Media

* Yoon, H (hyoon3@uiuc.edu), University of Illinois at Urbana-Champaign, 205 N Mathews Ave., Urbana, IL 61801, United States Valocchi, A J (valocchi@uiuc.edu), University of Illinois at Urbana-Champaign, 205 N Mathews Ave., Urbana, IL 61801, United States Werth, C J (werth@uiuc.edu), University of Illinois at Urbana-Champaign, 205 N Mathews Ave., Urbana, IL 61801, United States Oostrom, M (mart.oostrom@pnl.gov), Pacific Northwest National Lab, P.O. Box 999, MS K9-33, Richland, WA 99352, United States

The effects of the spatial distribution of soil permeability and water saturation, the NAPL spill scenario, water infiltration events, and vapor transport on NAPL distribution in two-dimensional (2-D) and three-dimensional (3-D) unsaturated porous media were investigated. 3-D homogeneous and heterogeneous fields were considered and a 2-D vertical cross-section along the center of the 3-D field was used for the 2-D simulations. The same NAPL and water infiltration rates over a source zone at the center of the top boundary were used in 2-D and 3-D simulations. The NAPL distribution was most strongly influenced by NAPL evaporation to the atmosphere and NAPL and water infiltration rates. The fraction of total NAPL mass that reached the groundwater table was higher in 2-D than in 3-D. The difference between 2-D and 3-D simulations can be primarily attributed to the following factors. First, water saturation was higher in 2-D than in 3-D because the water plume spread out more evenly due to the additional horizontal direction in the 3D case. Hence, NAPL can migrate vertically faster in 2-D than in 3-D due to the higher NAPL relative permeability in the former. Second, the effect of vapor transport in 3-D was more significant than in 2-D, mainly due to the presence of the additional horizontal direction for vapor transport in 3-D. Hence, more NAPL mass moved out of the NAPL source zone in the vadose zone in the 3-D simulation, resulting in a lower fraction of the total NAPL mass in groundwater. These simulations indicate that the 2-D simulation for organic compounds with high vapor pressure needs to be compared with the 3-D simulation in both homogeneous and heterogeneous unsaturated porous media. The effect of variability in the permeability field and quantitative analysis of dimensionality on NAPL distribution will be further explored through stochastic modeling.

H33H-1733 

A Variational Multiscale – High-Resolution Method for the Simulation of Unstable Multiphase Flow in Heterogeneous Formations

Dub, F (fxdub@mit.edu), Massachusetts Institute of Technology, Civil and Environmental Engineering 77 Massachusetts Ave, Room 48-319, Cambridge, MA 02139, United States * Juanes, R (juanes@mit.edu), Massachusetts Institute of Technology, Civil and Environmental Engineering 77 Massachusetts Ave, Room 48-319, Cambridge, MA 02139, United States

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. In this paper, we concentrate on the development of simulation techniques that account for the heterogeneity present in realistic reservoirs, and have the ability to capture (on coarse grids) the detailed pattern of unstable flows due to viscous fingering and channeling. We express the governing equations of multiphase flow as a pressure equation and a saturation equation. Both are nonlinear but are only weakly coupled. The pressure equation is elliptic, while the saturation equation is quasi-hyperbolic. Traditionally, the large degree of heterogeneity in the coefficients of the pressure equation has been tackled by upscaling the fine-scale properties to coarse-scale effective coefficients. Here, we avoid upscaling and propose a variational multiscale (VMS) method that splits the original problem is (rigorously) into a coarse-scale problem and a subgrid-scale problem. The framework is very flexible with respect to how each of these problems is approximated. 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. The method is therefore locally conservative at both the coarse and fine scales. We pay special attention to the definition of the local boundary conditions for the subgrid problems. In particular, we develop a well model, which accounts for subgrid heterogeneity and radial flow regime in a consistent fashion, without compromising the local mass conservation property. The saturation equation is then solved on the fine scale by a high-resolution explicit finite difference method that captures fingering and channeling when a low-viscosity fluid displaces a more viscous one. We present the application of the method to two-dimensional, highly heterogeneous problems. http://web.mit.edu/juanes/www/