Effects of fluid flow on inelastic deformation and failure in dilating and compacting rock
W. A. Olsson (505-844-7344) and D. J. Holcomb

This project is collaborative with J. Rudnicki of Northwestern University.

One of the proposed strategies for mitigation of carbon dioxide emissions is to sequester CO₂ in depleted oil or gas reservoirs for periods in excess of 50 to 100 years. The implementation of this strategy safely and efficiently will require an improved understanding of the effects of coupling between fluid flow and the inelastic deformation and failure of porous rock. It is critically important to understand how producing the original contents of the reservoir has changed the stress and deformation states and the properties of the rocks making up the reservoir trap. For example, has inelastic deformation of the reservoir during the original production degraded the porosity and permeability, which make reservoirs attractive for sequestration. Certain deformation structures, such as shear and compaction bands, often compartmentalize a hydrocarbon reservoir. Thus, inelastic deformation during earlier production (or injection of CO₂) may produce structures that interfere with fluid migration, destroying the reservoir as a potential CO₂ repository. Although the importance of coupling between mechanical coupling of fluid flow with inelastic deformation of geomaterials has long been recognized, understanding has progressed slowly because of the difficulty of making incisive measurements and the complexity of the theory.

We will apply coordinated experimental, theoretical, and numerical techniques to the problem, including both the theoretical and numerical aspects (Rudnicki) and the experimental part of the program (SNL). Deformation experiments will be performed on porous rock such as sandstone characteristic of oil and gas reservoirs. Three main objectives of this work will be: (1) to illuminate the phenomenology of dilatant/compactive rock containing pore fluid in the drained and the undrained response, and test predictions of current theory, (2) to collect the appropriate constitutive data constrained and suggested by developing theoretical models, and (3) to construct well-posed experiments that can be used as a validation test against a numerical implementation of the theory to be developed as part of this research. Theoretical work will attempt to extend layer analyses for dilating rock to more general deformation programs. The additional aspects of compaction of porous rock will be included. Porous rock is known to have a cap on the yield surface and this will be addressed in the considerations. The theoretical analytical model will be implemented into a numerical model for comparison to the validation experiment.

Micromechanical processes in porous geomaterials
J. T. Fredrich (505-844-2096)

This project is collaborative with T.-F. Wong of SUNY Stony Brook.

Research has progressed on two main fronts. (1) Microscale Flow Modeling in Porous Media. The first development focused on rigorous formulation of boundary conditions at the edges of the 3d image domain. The numerical lattice Boltzmann (LB) flow simulations require flow definition at the model ingress and exit; however, this condition is not definable a priori for flows in complex porous media. Instead, this requirement is typically handled by "wrapping" flow around the model boundaries so that the exit flow is prescribed at the inlet. Because of the geometric mismatch at the model boundaries for complex porous media, this requires mirroring the model domain so that the exit and inlet geometries are identically matched. This is required on all model boundaries, so that, in effect, the required model domain is 8x the size of the image domain. Because of the high cost associated with the increased model size, alternative formulations were developed to allow the flow conditions at the model boundaries to be defined implicitly, rather than explicitly through formally mirroring the image volume along each of the three primary axes. This resulted in a reduction of the necessary computation by almost an order of magnitude. The second major development related to assessing quantitatively the sensitivity of numerical flow simulations to image segmentation. The 3d image data that are mostly used in this work are gathered at 8-bit resolution equal to 256 grey levels. Approximately 10-15% of the image data fall in the "gray" region between the solid and void phases, meaning that those image voxels contain both solid and part void phases. This aspect is intrinsic to all image data, whether it is gathered by confocal microscopy, nuclear magnetic resonance imaging, or synchrotron computed microtomography. A large 3d data set for sandstone gathered at 1 micron resolution was segmented using three alternative thresholding schemes referred to as a simple segmentation, a local minima segmentation, and a kriging-based segmentation scheme. The first two are global schemes, whereas the third is a local scheme. 3d flow simulations were performed on the three different segmented data sets. As an alternative, we also developed a new numerical framework to allow the raw (unsegmented) image data to be mapped directly onto the LB flow simulation. The formulation, which we term the method of partially saturated computational cells, allows definition of partial flow for image voxels that contain a mixture of solid and void phases. While increasing the numerical computation, our results indicate that the method of partially saturated computational cells provides the most physically realistic approach for utilizing real-world image data. A paper by Fredrich, Noble, and O'Connor describing this work was presented at the fall meeting of the American Geophysical Union. The final work thrust under the first area relates to modeling the flow of melt in partially molten rocks and is being conducted collaboratively with Prof. David L. Kohlstedt at the University of Minnesota. Kohlstedt and his group have supplied "stacks" of optical images by serially sectioning samples that have been deformed at high pressures and temperatures in rock mechanics deformation apparatus. We are performing 3d flow simulations using the binary data, and initial results were presented by Scott, Groebner, Fredrich, O'Connor, and D. L. Kohlstedt at the fall meeting of the American Geophysical Union. This is a unique application of our capabilities as the permeability of the system of interest (partially molten mantle rock) can only be assessed through numerical simulation. (2) Micromechanics of Porous Geomaterials. Optical and scanning electron microscopy studies (SEM) were performed to elucidate the micromechanics of compaction during triaxial compression experiments performed on Castlegate sandstone. The experiments included acoustic emission detection and location, and indicated that initially, acoustic emissions were focused in horizontal bands that initiated at the sample ends (perpendicular to the maximum compressive stress), but with continued loading progressed axially towards the center. The microscopy, that included quantitative stereological measurements, confirmed the existence of fundamentally distinct zones of deformation within the test specimen. In conjunction with the macroscopic stress-strain and acoustic emission data, the microscopy analyses suggest that compaction of this weakly cemented sandstone under the triaxial load path followed during the test was accommodated in two distinct stages. The first stage of compaction was associated with breakage of the minimal cement bonding the framework grains and subsequent rotation of the intact framework grains to yield a more compact grain packing. This stage reduced the macroporosity by an estimated 27%. The second stage of compaction was associated with intense grain comminution and subsequent rotation of the grain fragments and resulted in substantial additional reduction in porosity and bulk volume. A manuscript based on this work by DiGiovanni, Fredrich, Holcomb, and Olsson has been submitted for publication in the Proc. of the Fourth Rock Mechanics Symposium. A second manuscript, by Fossum and Fredrich, was prepared for this meeting that relates generally to the deformation of porous materials (note that Fossum's work was not funded under this project). At low mean stresses, porous geomaterials fail by shear localization, and at higher mean stresses, they undergo strain-hardening behavior. Cap plasticity models attempt to model this behavior using a pressure-dependent shear yield and/or shear limit-state envelope with a hardening or hardening/softening elliptical end cap to define pore collapse. While these traditional models describe compactive yield and ultimate shear failure, difficulties arise when the behavior involves a transition from compactive to dilatant deformation that occurs before the shear failure or limit-state shear stress is reached. In this work, a continuous surface cap plasticity model was used to predict compactive and dilatant pre-failure deformation. During loading the stress point can pass freely through the critical state point separating compactive from dilatant deformation. The predicted volumetric strain goes from compactive to dilatant without the use of a non-associated flow rule. The new model is stable in that Drucker's stability postulates are satisfied. Finally, synchrotron computed microtomography experiments were conducted in collaboration with Dr. Mark L. Rivers and others at the Advanced Photon Source at Argonne National Laboratory in December. Large 3d data sets were gathered on several samples of sandstone, including samples deformed in the laboratory, as well as glass bead packs sintered to differing densities. The data gathered were the highest resolution data gathered at the GeoSoilEnviroCARS sector to date. We are currently analyzing the data sets.

The Role of Fracture Intersections in the Flow and Transport Properties of Rock
H. W. Stockman (702-363-8522)

This project is collaborative with S. Brown of New England Research and A. Caprihan of New Mexico Resonance.

We propose to address several problems concerning flow and transport through fracture intersections, using a combination of quantitative measurements and visual observation in real roughwalled fractures. We will make a variety of observations including channeling in single fractures, channeling at and through single fracture intersections, and channeling through fracture networks. These studies will include transport of solutes and multiple immiscible phases, which will be used to develop quantitative models describing the channeling and mixing behavior in 3-D intersections of rough-walled fractures. Our basic hypothesis - mixing at fracture intersections depends on the interplay between compartmentalized flow in the channels and enhanced dispersion due to fine-scale surface roughness within channels - will be tested directly by experimental observation. This work makes use of our abilities to (1) quantify surface roughness through profilometry, (2) construct high quality replicas of actual fractured rock in transparent plastics, (3) simulate surface roughness with the computer, (4) construct bench-scale fracture systems, (5) observe actual flow and transport with quantitative visualization techniques, and (6) conduct detailed numerical simulations of flow and transport in heterogeneous media using lattice-Boltzmann methods. The proposed experimental observations and numerical studies will provide overlapping and complementary information essential for complete understanding of flow channeling in fractures. The experimental techniques will provide data for building and testing of models. The experiments will allow the complexity of real fracture roughness and intersection geometries to be considered in detail. Numerical studies, on the other hand, will provide insight into processes which are not easily observed in the laboratory (e.g., wide range of Peclet numbers and side-wall porosity and storage in porous rocks). The experimental methods proposed here make extensive use of impermeable, transparent replicas where the study of porous wall effects are not possible.

Results of experiments run by Steve Brown suggested the need to modify the algorithms to allow fluid sinks and sources of arbitrary shapes, so we could simulate injection of dyes near fracture intersections (thereby mimicking experiments). The code was modified to allow such sources, and was tested to verify that accuracy of the flow field was maintained. Major revisions were made to the code to allow very efficient modeling of 3-D fracture junctions with arbitrary offsets. These modifications make it possible to match the scale of the lab experiments. The user interface was improved to simplify the specification and visualization of the intersecting fractures.

In my two years with the Yucca Mountain Waste Package (WP) Design Group, I discovered several applications of LB methods that could help solve significant problems. I am developing a version of the code to model simultaneous heat flow and latent heat, convection, diffusion and corrosion of steel, glass and fuel in the WPs. (The same code will be capable of modeling metallurgical processes and circuit board cooling.) The specific problem I wish to solve is the calculation of redox state and actinide solubility in flooded WPs. Traditionally at YMP, oxidizing conditions in the WP are assumed, leading to grossly conservative estimates of plutonium release and potential for criticality. I recently developed and qualified a version of EQ6, to be used in benchmark comparisons with the LB studies; this code has shown that under reasonable bounding conditions, the competition of fuel degradation and oxygen diffusion maintains very low Pu solubility. It is now left to the LB codes to provide realistic models of convective cell size, oxygen exchange, and corrosion boundary thickness in the WPs.

Three Dimensional Transient Electromagnetic Inversion
G. A. Newman (505-844-8158)

This project is collaborative with D. Alumbaugh at the University of Wisconsin at Madison.

This project has been developing techniques to construct and appraise solutions to the three-dimensional (3D) electromagnetic (EM) inverse problem using full wave equation modeling. We have focused on large scale inverse solutions directly in the frequency domain for both controlled source and magnetotelluric applications. Our philosophy has been to avoid approximations which sacrifice solution accuracy for speed. The need for a full solution to this problem is critical since is serves as an accuracy benchmark on approximate methods now being implemented and tested on workstation platforms. While our work has focused on solving inverse problems in the frequency domain, we have come to realize that there is also a critical need to develop a data inversion capability directly in the time domain, especially given the limitations frequency domain data inversion as outlined in this proposal. Solution of the 3D inverse problem in the time domain has been avoided because the problem is far more difficult compared to the frequency domain problem. Nevertheless our work on the frequency domain problem has provided us with experience and techniques necessary to tackle the problem. The renewal proposal will emphasize the application of 3D transient electromagnetic (TEM) inversion to subsurface site characterization.

This project has submitted several publications on inversion, modeling and appraisal this last year. These works represent a culmination of findings for the BES funding over the period 1996-1999 for frequency domain geophysical methods. With the project renewed for another three-year period, beginning in October 1999, we have begun to focus on the modeling and inverse problems for time domain electromagnetic methods. Currently we are focusing on 3D forward modeling methods that employ explicit time stepping and spectral Lanczos methods. These forward solutions will provide the computational kernel in the 3D inverse solution. Once this later solution has been developed, we will apply it to CO₂ sequestration problems.

Inversion of full waveform seismic data for shallow subsurface properties
D. F. Aldridge (505-844-2823), G. J. Elbring, and D. Womble

High resolution images of the shallow subsurface are of critical importance in the environmental characterization and remediation fields. Full waveform inversion of seismic data offers the potential for providing improved subsurface images. This is particularly true for shallow seismic surveys, where reflected energy may be poor quality or even nonexistent. A full waveform inversion procedure requires robust algorithms for i) realistic modeling of seismic wave propagation in a subsurface geologic environment, and ii) inverting recorded seismic waveforms to recover the three-dimensional (3D) distribution of earth model properties. This research project is focused on developing the requisite forward modeling and inversion algorithms, and siting them in the appropriate computational environments.

A seismic wave propagation algorithm appropriate for 3D isotropic elastic media has been developed. The algorithm is based on the velocity-stress equations of linear elastodynamics, a set of nine, coupled, first-order partial differential equations. Solution yields the three components of the particle velocity vector, and the six independent components of the stress tensor. An explicit, time-domain, finite-difference (FD) numerical scheme is used to solve these equations. Centered spatial and temporal FD operators possess 4th-order and 2nd-order accuracy in the discretization intervals, respectively. The nine dependent variables are stored on uniform, but staggered, spatial and temporal grids.

The algorithm is a direct numerical implementation of the governing partial differential equations of isotropic elasticity. No theoretical approximations, such as far-field distances, high frequencies, weak scattering, or one-way wave propagation, are adopted. Hence, the algorithm generates all seismic arrival types (P-waves, S-waves, reflections, refractions, multiples, mode-conversions, diffractions, head waves, surface waves, etc.) with fidelity, provided spatial and temporal gridding intervals are chosen appropriately. In order to simulate seismic conditions at the earth's surface, an explicit representation of a plane zero-stress boundary is incorporated into the algorithm. Point traction sources, with arbitrary orientations, may be imposed on this surface. Nonplane surface topography will be introduced into the algorithm in the future. However, preliminary testing indicates that a nonplane stress-free surface can be mimicked by assigning the material properties of air to grid nodes above the surface. Algorithm stability is maintained by making the air/earth interface gradational.

Significant computational resources are required to execute this 3D elastic wave propagation algorithm. In order to treat realistic-sized earth models (e.g.., tens of millions of grid nodes) and trace durations (e.g., thousands of timesteps) within reasonable execution times, parallel versions of the algorithm have been developed and sited on appropriate computational platforms. Parallel codes utilized Message Passing Interface (MPI) or Parallel Virtual Machine (PVM) coding protocols to maintain algorithm portability. Synthetic seismograms replicating borehole experiments conducted at Bayou Choctaw Salt Dome, Louisiana have been computed, and are a reasonable match to the field recorded data, albeit with reduced spectral bandwidth. Future algorithm improvements will be directed toward i) understanding numerical dispersion and stability conditions more rigorously, ii) improving the absorbing boundary conditions imposed on the flanks of the 3D grid, iii) improving the frequency bandwidth of the computed signal.

The physics of two-phase immiscible fluid flow in single fractures and fractured rock
R. J. Glass (505-844-5506)

This project is collaborative with H. Rajaram of the University of Colorado at Boulder and M. Nicholl at Oklahoma State University.

Under two-phase, immiscible fluid-flow conditions, phase-geometry within individual fractures and the fracture network (i.e., the geometry saturated with each phase) ultimately controls: the permeability to each phase, fluid pressure/saturation relations, solute dispersion within each phase, matrix block connectivity. Phase geometry within a single fracture is a function of both the aperture field and the two-phase flow processes themselves. Capillary, gravitational, and viscous forces in combination with boundary and initial conditions have all been demonstrated to play roles in the formation of fracture phase-geometry. Smooth to fingered phase-structures form depending on the interplay of these forces emphasizing the importance of understanding the two-phase displacement processes themselves. Exploration of more complicated fracture/matrix and fracture network systems are also required to properly connect two-phase single fracture flow to the behavior of a fractured rock mass. Our approach to building understanding of two-phase flow in fractured rock couples systematic physical experimentation in a variety of specially designed analog systems with concurrent conceptual/numerical model formulation and simulation.

This year we will continue to: develop a generalized form of Reynolds equations which accounts for tortuosity of the fracture plane to better model single phase flow through fractures and to modify our solute transport model to include hydrodynamic boundary layers around entrapped phase. We also plan to design and conduct final experiments, after which we shall compare results to our models for flow and transport. We will conduct a set of experiments varying the bond number for small capillary number and complete our verification of the MIP model; preliminary experiments to study the entrapped phase dissolution process in fractures will begin. We plan to initiate development of a phase dissolution model which combines flow and transport with inter-phase transfer across the phase-phase menisci that uses MIP to properly shrink entrapped ganglia and compare our model results to experiments. We also will continue development and experimentation in our single fracture/matrix and fracture network experiments.

Laboratory Investigation of Permeability Upscaling
V. C. Tidwell (505-844-6025)

This project is collaborative with J. Wilson at New Mexico Institute of Mining and Technology

Parameterization of porous media flow and transport models is often complicated by the inability to make measurements at the desired scale of analysis. This disparity in scales necessitates the use of some averaging or upscaling model to compute the required effective media properties. Although numerous theoretical upscaling models have been proposed, physical data with which to test these models are sparse and limited in scope. To meet this need laboratory investigation of permeability upscaling is being conducted with the use of a novel minipermeameter test system, which we call the Multi-Support Permeameter (MSP). The MSP allows precise, rapid, non-destructive measurement of permeability over a range of different sample supports (i.e., sample volumes). The data are uniquely qualified for investigating permeability upscaling because of the consistency in measurement maintained across the different sample supports and the spatial detail captured in each of the permeability data sets. Systematic investigation of the influence of key porous media characteristics on permeability upscaling is approached through the careful selection of rock samples for analysis (i.e., different depositional and diagenetic environments). Finally, efforts are made to interpret the measured data through comparisons drawn with different theoretical upscaling models.

Key to inverting estimates of permeability from MSP measurements are the geometric factors that account for the divergent flow field geometry around the tip seal. Using a finite element method, new geometric factors for a two-layer permeability system were developed. These factors can be used to determine the permeability of a surface layer when estimates of that layers thickness and the permeability of the basal layer are independently available. This work expands the possible range of minipermeameter application, and yields further insight into the appropriate application of the minipermeameter under classical homogeneous assumptions. Efforts are currently being made to apply these results to determine the importance of surface disturbances from sawing and other rock preparation on measurements of the underlying permeability.

Intuitively we recognize that permeability is in some way related to the textural and structural attributes visible in the host rock; however, establishing such a relationship is problematic. The ability to discern such relationships is of significant value to site characterization programs, as this type of information would aid in predicting permeability at unsampled locations. We explore the potential for such relationships in our data by comparing measured permeabilities from two of our sandstone samples and one tuff sample with digital image data taken from the corresponding block faces. We found direct correlation between the permeability and any textural measure drawn from the digital image data to be very difficult to establish. Alternatively, correlation in the spatial patterns was easily affirmed, beginning with visual comparison of the digital image and permeability data. Quantitative comparison of the spatial patterns was then performed by way of semivariograms calculated from the digital image and permeability data. Results showed consistency in the structure, length-scales, and anisotropy ratios characterizing the semivariograms. Additionally, a Difference-of-Gaussian edge filter was used to segment (delineate) the spatial features in the permeability maps and digital images. Cumulative distribution functions characterizing the size, shape, and orientation of the features were found to be statistically indistinguishable between an image and its corresponding permeability map, but statistically different between rock samples. For visually sharp features, correlation extended to the positioning of the feature in both the digital image and permeability map.

Continuum- and Particle-Level Modeling of Concentrated Suspension Flows
L. A. Mondy (505-844-1755)

This project is collaborative with M. Ingber of the University of New Mexico and A. Graham of Los Alamos National Laboratory.

The purpose of this program is to combine experiments, computations, and theory to make fundamental advances in our ability to predict transport phenomena in concentrated, multiphase, disperse systems, particularly when flowing through geologic media. The ability to accurately model multiphase flow is not only central to petroleum production applications, but also to understanding the transport of natural debris and magma, and the underground disposal of hazardous wastes. We are developing a new particle-level simulation capability, based on a fast multipole method and implemented on massively parallel computers, that will allow the incorporation of an adequate number of suspended particles to determine microstructural information that is either impossible or impractical to obtain through laboratory experiments. This information will be used in the development and refinement of continuum-based approaches. Both levels of models will be validated with actual experimental data where possible.

This year we continued improving the state-of-the-art continuum approach to suspension flow modeling, including effects of non-Newtonian (shear-thinning) suspending fluids, non-neutrally buoyant particles, and flows with significant curvature. For the first time, we believe that the normal forces that become appreciable with curved streamlines are correctly accounted for. Furthermore, the model is frame indifferent. Computational results in two-dimensional flows compare well with laboratory experiments. The model is being extended to three-dimensional flows, with appropriate laboratory experiments being performed concurrently.

Transport Visualization for Studying Mass Transfer and Solute Transport in Permeable Media
L. C. Meigs (505-844-2375)

This project is collaborative with R. Haggerty at Oregon State University and C. Harvey at Massachusetts Institute of Technology.

Laboratory experiments and numerical analyses are being conducted to determine the conditions under which physical mass transfer occurs, and to develop models for predicting the effects of mass transfer over long time scales. Applications of this research are to CO₂ sequestration, nuclear waste disposal, petroleum extraction, and groundwater remediation. Although large bodies of literature describe macro-dispersive transport in multi-Gaussian permeability fields, and mass transfer behavior under dual porosity conditions, little previous work has taken an empirical approach to determine how the spatial structure of permeability in natural rock controls transport behavior. We will be determining the types of permeability fields that lead to transport dominated by mass transfer, and the types of permeability fields that lead to transport dominated by macro-dispersion. We will employ X-ray imaging in actual rock slabs, and transmitted light imaging in constructed media, to gain a detailed understanding of solute transport in heterogeneous material. Such detailed spatial data have been previously unavailable for natural rock. We will construct high-resolution "maps" of flow and transport properties, and use these maps to determine how the spatial structure of small-sale properties control large-scale solute transport.

We also will develop methods to determine the distribution of time-scales (or equivalently, rate coefficients) of mass transfer from the physical and geometric characteristics of rock. A growing body of research demonstrates that mass transfer occurs over multiple time-scales; that conventional, single-rate, models of mass transfer fail to account for the important effects of geometric and chemical complexity found in natural permeable media. We will consider the geometry of flow and transport parameters that lead to multiple rates of mass transfer. Experiments will be conducted at different length- and time-scales, so observations made at one scale may be compared with predictions based on data from a different scale. In particular, we will focus on developing a better understanding of the physical causes of multiple rates of mass transfer, and developing better methods for forecasting the long-term behavior of solutes.

We have made progress in mapping transport regimes based on common dimensionless numbers. For a binary medium (i.e., a medium composed of two hydraulic conductivities) we are able to express the conditions under which a mass transfer formulation is most appropriate and when another formulation (e.g., macrodispersion, advective mass transfer, etc.) is more appropriate. We are in the process conducting numerical simulations to extend our map of transport regimes to non-binary media.

We are beginning to plan all of the experiments for the project, and we have made significant progress in the planning for the next year. First, we will be conducting transmitted-light experiments in sand tanks to test our hypotheses regarding transport regimes. We will be running experiments in a range of media to confirm or reject our hypotheses regarding diffusive mass transfer, advective mass transfer, and macrodispersion regimes. These experiments will be done using a new fluorescing tracer imaging system, which will allow us to obtain concentration maps over a large dynamic range. Second, we will be conducting a series of experiments to test our understanding of multiple time-scales of mass transfer. These experiments will be in media of different length but with the same pore-water velocity, allowing us to test differentiate between non-Fickian dispersion and multiple time-scales of mass transfer.

Atomistic simulation of clay minerals and their interaction with hazardous wastes: Molecular orbital and empirical methods
R. T. Cygan (505-844-7216)

We are incorporating a suite of theoretical tools, that use both empirical and molecular orbital methods, to provide a fundamental understanding of clay mineral bulk structure, surface structure, and fluid-clay sorption interactions. In so doing, we draw upon our previous experience in experimental and spectroscopic examination of sorption, dissolution, and growth of clay phases. Modeling of clay and clay-fluid systems combines energy minimization and molecular dynamics methods using both molecular cluster and periodic systems representations of the clay minerals and the interacting fluids. Substantially more sophisticated calculations that incorporate state-of-the-art quantum chemistry methods complement the empirical models. We use an ab initio approach based on local density functional (LDF) theory to examine the electronic structure of various clay minerals. We are collaborating with leading experts in the structural analysis of clay minerals and in the spectroscopic examination of mineral surfaces. This approach allows us to appropriately modify our simulation approach and/or empirical parameters in order to better simulate the real clay system. The computer simulations also provide a convenient working model for interpreting complex spectroscopic data, that can be continually refined in order to match observed spectral characteristics.

Results to date include the successful analysis of a metal-organic complex within the interlayer of a smectite clay. Molecular modeling of the large tributyl phosphate (TBP) complex of europium provides a test of the sensitivity of force field calculations to the molecules within the interlayer of a trioctahedral hectorite clay. The TBP is used in the processing of actinide materials and often contributes to environmental problems. We have also developed a fully flexible and general set of interatomic potentials for the analysis of simple oxide, hydroxide, oxyhydroxide, and clay phases. This empirical forcefield provides the basis for evaluating numerous minerals common to soils that may potentially interact with environmentally hazardous chemical and radioactive waste. We have successfully modeled the swelling behavior of smectite clays using this forcefield whereby we mimic the observed changes in the basal d-spacing of clay as a function of water content. Not only do we identify the changes in the interlayer hydration but we observe the fine structure in the swelling curve as the water molecules fill in the interlayer between abrupt expansions. We have also examined the intercalation of hydrazine in kaolinite clay. Hydrazine is a component of jet and rocket fuels and is often the primary waste component at government waste sites. We have used empirical and DFT methods to examine the molecular behavior of hydrazine as an intercalating compound. Theoretical structures are in agreement with experimental and spectroscopic observation. Most importantly, the models demonstrate the details of hydrogen bonding within the clay interlayer and how it changes as a function of hydrazine loading. We are working with David Bish (Los Alamos Laboratory) and Cliff Johnston (Purdue University) on the hydrazine project. Recently, we have simulated the surface structure of calcite (dry and wet) using a shell model and molecular dynamics. We are collaborating with Neil Sturchio (Argonne National Laboratory) for comparison of our model results with his experimental X-ray absorption structures. Additionally, we are developing a model for a mixed metal layered hydroxide phase, hydrocalumite, that can be applied in the environmental cleanup of anion species. This work is being done in collaboration with James Kirkpatrick and Andrey Kalinichev at the University of Illinois.