Earth Sciences

TOUGHREACT is a numerical simulation program for chemically reactive non-isothermal flows of multiphase fluids in porous and fractured media. The program was written in Fortran 77 and developed by introducing reactive geochemistry into the multiphase fluid and heat flow simulator TOUGH2. A variety of subsurface thermophysical-chemical processes are considered under a wide range of conditions of pressure, temperature, water saturation, ionic strength, and pH and Eh. Interactions between mineral assemblages and fluids can occur under local equilibrium or kinetic rates. The gas phase can be chemically active. Precipitation and dissolution reactions can change formation porosity and permeability. The program can be applied to many geologic systems and environmental problems, including geothermal systems, diagenetic and weathering processes, subsurface waste disposal, acid mine drainage remediation, contaminant transport, and groundwater quality. Here we present two examples to illustrate applicability of the program: (1) injectivity effects of mineral scaling in a fractured geothermal reservoir and (2) CO 2 disposal in a deep saline aquifer.

A large-scale underground thermal test (Drift Scale Test-DST) in fractured volcanic tuff resulted in changes to water and gas chemistry as well as mineral precipitation and dissolution in fractures. Thermal, hydrological, and chemical (THC) processes in the DST were modeled by Lawrence Berkeley National Laboratory ''LBNL'' and Japan Nuclear Cycle Development Institute ''JNC'' as part of the international working group DECOVALEX. Predictions of THC processes in the DST for the 4-year heating and 4year cooling periods were initially performed by the LBNL group, with the current model reflecting a revised heater operation history and model. JNC used primarily the original data from the prediction and created a new model to evaluate a selected set of data. The approaches taken by the groups differed in several ways and a comparison of the methodologies and results of the simulations allow for a better understanding of modeling coupled processes in unsaturated fractured rock. The LBNL model represented the fractures and rock matrix as a fully interacting dual-continuum (in terms of fluid, heat, and chemical transport) with the local mineral-water-gas reactions treated by kinetic and equilibrium reactions. The JNC model represented the fractures and matrix as a single effective continuum, with equilibrium mineral-water reactions controlling the chemical evolution. Both models considered aqueous species transport, with gas phase CO 2 transport only considered in the LBNL model. Comparisons to data collected from the DST illustrate the behavior of the models and their ability to capture the relevant THC processes. Overall, both models capture the temperature evolution in the rock quite closely, although the JNC model gave a closer match to the initial temperature rise in the rock, likely owing to the use of site-specific thermal data as opposed to average properties used for the LBNL model. Both models showed the contrasting solubility effects of increasing temperature on calcite and silica solubility; yet the dualcontinuum approach better represented the effects of boiling and condensation on aqueous species chemistry and the distribution of mineral precipitation. r

ECO2N is a fluid property module for the TOUGH2 simulator (Version 2.0) that was designed for applications involving geologic storage of CO 2 in saline aquifers. It includes a comprehensive description of the thermodynamics and thermophysical properties of H 2 O-NaCl-CO 2 mixtures, that reproduces fluid properties largely within experimental error for the temperature, pressure and salinity conditions of interest (10°C 6 T 6 110°C; P 6 600 bar; salinity up to full halite saturation). Flow processes can be modeled isothermally or non-isothermally, and phase conditions represented may include a single (aqueous or CO 2 -rich) phase, as well as two-phase mixtures. Fluid phases may appear or disappear in the course of a simulation, and solid salt may precipitate or dissolve. ECO2N can model superas well as sub-critical conditions, but it does not make a distinction between liquid and gaseous CO 2 and hence is not applicable for processes that involve two CO 2 -rich phases. This paper highlights significant features of ECO2N, and presents illustrative applications.

TOUGHREACT is a numerical simulation program for chemically reactive non-isothermal flows of multiphase fluids in porous and fractured media. The program was written in Fortran 77 and developed by introducing reactive geochemistry into the multiphase fluid and heat flow simulator TOUGH2. A variety of subsurface thermo-physical-chemical processes are considered under a wide range of conditions of pressure, temperature, water saturation, ionic strength, and pH and Eh. Interactions between mineral assemblages and fluids can occur under local equilibrium or kinetic rates. The gas phase can be chemically active. Precipitation and dissolution reactions can change formation porosity and permeability. The program can be applied to many geologic systems and environmental problems, including geothermal systems, diagenetic and weathering processes, subsurface waste disposal, acid mine drainage remediation, contaminant transport, and groundwater quality. Here we present two examples to illustrate applicability of the program:

Correlations presented by Spycher et al. (2003) to compute the mutual solubilities of CO 2 and H 2O are extended to include the effect of chloride salts in the aqueous phase. This is accomplished by including, in the original formulation, activity coefficients for aqueous CO 2 derived from several literature sources, primarily for NaCl solutions. Best results are obtained when combining the solubility correlations of Spycher et al. (2003) with the activity coefficient formulation of Rumpf et al. (1994) and Duan and Sun (2003), which can be extended to chloride solutions other than NaCl. This approach allows computing mutual solubilities in a noniterative manner with an accuracy typically within experimental uncertainty for solutions up to 6 molal NaCl and 4 molal CaCl 2.

TOUGHREACT ) is a generalpurpose reactive geochemical transport numerical simulator. It deals with multiphase flow, solute transport and geochemical reactions including aqueous complexation, mineral dissolution/ precipitation and cation exchange. Making use of an extended Debye-Hückel ion activity model, this simulator can handle solutions concentrated to slightly above ~1 molal with caution, and only for NaCl-dominant waters at ionic strengths up to ~4 molal. However, brines produced under natural and artificial conditions are often more concentrated. To handle such brines, a Pitzer activity model was implemented in TOUGHREACT, based on the standard Harvie-Moller-Weare (HMW) formulation that accounts for all binary and ternary combinations of interaction terms. The vapor-pressure-lowering effect caused by the low water activity in brines was also accounted for by this code. The extended version was verified and tested using published results from laboratory experiments and benchmarked against other computer codes. This new version of TOUGHREACT is being applied to the investigation of boiling and evaporation within and around the proposed high-level nuclear waste emplacement tunnels at Yucca Mountain, Nevada. An example application is presented. Processes considered in the example include evaporation of porewater to near dryness, formation of highly concentrated brines, precipitation of deliquescent salts, and generation of acid gases.

Although bioclogging has been used for various applications, its performance has been unpredictable. We carried out a bioclogging experiment and developed a reactive transport model to understand and quantify geophysical responses to the production of a biopolymer (dextran) in a silica sand porous medium. The developed model was used to predict effectiveness of bioclogging and the associated geophysical signature under different treatment conditions. The porosity reduction calculated from the reactive transport modeling matched the concurrent increase in electrical resistivity and confirmed the dominant contribution of pore fluid conductivity to the bulk electrical conductivity. Correlations between measured electrical phase signals and predicted mineral surface area reduction indicated the controlling effect of mineral surface area on polarization. The joint use of the two methods quantitatively estimated the extent of decrease in porosity as well as effective mineral surface area, both of which are challenging to measure nondestructively in biological systems. The joint use of the geophysical methods and reactive transport modeling provided a better understanding of the dominant processes that govern dextran-based bioclogging and the associated geophysical responses. This opens new opportunities for quantifying changes associated with bioclogging processes and can potentially lead to new conceptual models that relate hydraulic property changes to the porosity and mineral surface area, which will enhance prediction, monitoring, and control of applications that use bioclogging.

Correlations presented by Spycher et al. (2003) to compute the mutual solubilities of CO 2 and H 2 O are extended to include the effect of chloride salts in the aqueous phase. This is accomplished by including, in the original formulation, activity coefficients for aqueous CO 2 derived from several literature sources, primarily for NaCl solutions. Best results are obtained when combining the solubility correlations of Spycher et al. with the activity coefficient formulation of and Duan and Sun (2003), which can be extended to chloride solutions other than NaCl. This approach allows computing mutual solubilities in a noniterative manner with an accuracy typically within experimental uncertainty for solutions up to 6 molal NaCl and 4 molal CaCl 2 .

Correlations are presented to compute the mutual solubilities of CO2 and chloride brines at temperatures 12–300°C, pressures 1–600 bar (0.1–60 MPa), and salinities 0–6 m NaCl. The formulation is computationally efficient and primarily intended for numerical simulations of CO2-water flow in carbon sequestration and geothermal studies. The phase-partitioning model relies on experimental data from literature for phase partitioning between CO2 and NaCl brines, and extends the previously published correlations to higher temperatures. The model relies on activity coefficients for the H2O-rich (aqueous) phase and fugacity coefficients for the CO2-rich phase. Activity coefficients are treated using a Margules expression for CO2 in pure water, and a Pitzer expression for salting-out effects. Fugacity coefficients are computed using a modified Redlich–Kwong equation of state and mixing rules that incorporate asymmetric binary interaction parameters. Parameters for the calculation of activity and fugacity coefficients were fitted to published solubility data over the P–T range of interest. In doing so, mutual solubilities and gas-phase volumetric data are typically reproduced within the scatter of the available data. An example of multiphase flow simulation implementing the mutual solubility model is presented for the case of a hypothetical, enhanced geothermal system where CO2 is used as the heat extraction fluid. In this simulation, dry supercritical CO2 at 20°C is injected into a 200°C hot-water reservoir. Results show that the injected CO2 displaces the formation water relatively quickly, but that the produced CO2 contains significant water for long periods of time. The amount of water in the CO2 could have implications for reactivity with reservoir rocks and engineered materials.

Correlations presented by Spycher et al. (2003) to compute the mutual solubilities of CO 2 and H 2 O are extended to include the effect of chloride salts in the aqueous phase. This is accomplished by including, in the original formulation, activity coefficients for aqueous CO 2 derived from several literature sources, primarily for NaCl solutions. Best results are obtained when combining the solubility correlations of Spycher et al. with the activity coefficient formulation of and Duan and Sun (2003), which can be extended to chloride solutions other than NaCl. This approach allows computing mutual solubilities in a noniterative manner with an accuracy typically within experimental uncertainty for solutions up to 6 molal NaCl and 4 molal CaCl 2 .

A field experiment involving the release of carbon dioxide (CO 2 ) into a shallow aquifer was conducted near Bozeman, Montana, during the summer of 2008, to investigate the potential groundwater quality impacts in the case of leakage of CO 2 from deep geological storage. As an essential part of the Montana State University Zero Emission Research and Technology (MSU-ZERT) field program, food-grade CO 2 was injected over a 30 day period into a horizontal perforated pipe a few feet below the water table of a shallow aquifer. The impact of elevated CO 2 concentrations on groundwater quality was investigated by analyzing water samples taken before, during, and following CO 2 injection, from observation wells located in the vicinity of the injection pipe, and from two distant monitoring wells. Field measurements and laboratory analyses showed rapid and systematic changes in pH, alkalinity, and conductance, as well as increases in the aqueous concentrations of naturally occurring major and trace element species.

The Yellowstone geothermal system provides an ideal opportunity to test the ability of reactive transport models to simulate the chemical and hydrological effects of water-rock interaction. Previous studies of the Yellowstone geothermal system have characterized water-rock interaction through analysis of rocks and fluids obtained from both surface and downhole samples. Fluid chemistry, rock mineralogy, permeability, porosity, and thermal data obtained from the Y-8 borehole in Upper Geyser Basin were used to constrain a series of reactive transport simulations of the Yellowstone geothermal system using TOUGHREACT. Three distinct stratigraphic units were encountered in the 153.4 m deep Y-8 drill core: volcaniclastic sandstone, perlitic rhyolitic lava, and nonwelded pumiceous tuff. The main alteration phases identified in the Y-8 core samples include clay minerals, zeolites, silica polymorphs, adularia, and calcite. Temperatures observed in the Y-8 borehole increase with depth from sub-boiling conditions at the surface to a maximum of 169.8°C at a depth of 104.1 m, with near-isothermal conditions persisting down to the well bottom. 1-D models of the Y-8 core hole were constructed to simulate the observed alteration mineral assemblage given the initial rock mineralogy and observed fluid chemistry and temperatures. Preliminary simulations involving the perlitic rhyolitic lava unit are consistent with the observed alteration of rhyolitic glass to form celadonite.

Plugging of flow paths caused by mineral precipitation in fractures above the potential repository at Yucca Mountain, Nevada, would reduce the probability of water seeping into the repository. As part of an ongoing effort to evaluate thermal-hydrologic-chemical (THC) effects on flow in fractured media, we performed a laboratory experiment and numerical simulations to investigate mineral dissolution and precipitation under anticipated temperature and pressure conditions in the repository. To replicate mineral dissolution by vapor condensate in fractured tuff, water was flowed through crushed Yucca Mountain tuff at 94°C. The resulting steady-state fluid composition had a total dissolved solids content of about 140 mg/L; silica was the dominant dissolved constituent. A portion of the steady-state mineralized water was flowed into a vertically oriented planar fracture in a block of welded Topopah Spring Tuff that was maintained at 80°C at the top and 130°C at the bottom. The fracture began to seal with amorphous silica within five days. A 1-D plug-flow numerical model was used to simulate mineral dissolution, and a similar model was developed to simulate the flow of mineralized water through a planar fracture, where boiling conditions led to mineral precipitation. Predicted concentrations of the major dissolved constituents for the tuff dissolution were within a factor of 2 of the measured average steady-1 Correspondence should be addressed to Patrick Dobson (pfdobson@lbl.gov; fax (510) 486-6115) Revised Draft 10/10/02 2 state compositions. The mineral precipitation simulations predicted the precipitation of amorphous silica at the base of the boiling front, leading to a greater than fifty-fold decrease in fracture permeability in 5 days, consistent with the laboratory experiment. These results help validate the use of a numerical model to simulate THC processes at Yucca Mountain. The experiment and simulations indicated that boiling and concomitant precipitation of amorphous silica could cause significant reductions in fracture porosity and permeability on a local scale. However, differences in fluid flow rates and thermal gradients between the experimental setup and anticipated conditions at Yucca Mountain need to be factored into scaling the results of the dissolution/precipitation experiments and associated simulations to THC models for the potential Yucca Mountain repository.

The Yellowstone geothermal system provides an ideal opportunity to test the ability of reactive transport models to accurately simulate water-rock interaction. Previous studies of the Yellowstone geothermal system have characterized water-rock interaction through analysis of rocks and fluids obtained from both surface and downhole samples. Fluid chemistry, rock mineralogy, permeability, porosity, and thermal data obtained from the Y-8 borehole in Upper Geyser Basin were used to constrain a series of reactive transport simulations of the Yellowstone geothermal system using TOUGHREACT. Three distinct stratigraphic units were encountered in the 153.4 m deep Y-8 drill core: volcaniclastic sandstone, perlitic rhyolitic lava, and nonwelded pumiceous tuff. The main alteration phases identified in the Y-8 core samples include clay minerals, zeolites, silica polymorphs, adularia, and calcite. Temperatures observed in the Y-8 borehole increase with depth from sub-boiling conditions at the surface to a maximum of 169.8°C at a depth of 104.1 m, with near-isothermal conditions persisting down to the well bottom. 1-D models of the Y-8 core hole were constructed to determine if TOUGHREACT could accurately predict the observed alteration mineral assemblage given the initial rock mineralogy and observed fluid chemistry and temperatures. Preliminary simulations involving the perlitic rhyolitic lava unit are consistent with the observed alteration of rhyolitic glass to form celadonite.

There is a complex relationship between seismic attributes, including the frequency dependence of reflections and fluid saturation in a reservoir. Observations in both laboratory and field data indicate that reflections from a fluid-saturated layer have an increased amplitude and delayed traveltime at low frequencies, when compared with reflections from a gas-saturated layer. Comparison of laboratory-modeling results with a diffusive-viscous-theory model show that low (<5) values of the quality factor Q can explain the observations of frequency dependence. At the field scale, conventional processing of time-lapse VSP data found minimal changes in seismic response of a gas-storage reservoir when the reservoir fluid changed from gas to water. Lowfrequency analysis found significant seismic-reflectionattribute variation in the range of 15-50 Hz. The field observations agree with effects seen in laboratory data and predicted by the diffusive-viscous theory. One explanation is that very low values of Q are the result of internal diffusive losses caused by fluid flow. This explanation needs further theoretical investigation. The frequencydependent amplitude and phase-reflection properties presented in this paper can be used for detecting and monitoring fluid-saturated layers.

There are numerous laboratory and field examples in which low-frequency components of reflected seismic waves show surprising imaging capabilities. Ironically, such components are often filtered out as useless in conventional data processing. However, as we demonstrate below, this part of the signal contains the most important information about the reservoir.

A complete and exact solution for the problem of an incident P wave scattered by an elastic spherical inclusion is presented and described. The solution can be obtained from either analytical formulae or stable numerical procedures. A method of estimating the number of terms that must be retained in the harmonic series in order to achieve a specified accuracy is given. The results are investigated by calculating synthetic seismograms, scattering diagrams, and scattering cross-sections for a broad frequency band and for both low-velocity and high-velocity inclusions. The fields within the shadow zone are formed primarily from three different types of waves, P waves transmitted through the sphere, P waves diffracted around the sphere, and S waves converted at the boundary of the sphere. The relative contribution from these different waves depends upon the distance of the observation point from the sphere.