Research Projects


CO2-Water Interfacial Tension

In using geological carbon sequestration (GCS) as a CO2 abatement technology, interfacial tension (IFT) between supercritical CO2 and the aqueous phase in a reservoir is the principal physical property that determines the storage capacity of the reservoir and the storage pressure above which CO2 may flow through the sealing cap rock overlying the storage unit. Very limited high-quality experimental IFT data are available, so alternative methods of predicting IFT over the wide range of storage unit pressures (5-45 MPa) and temperatures (30-110 °C) are needed. Thus we are using MD simulation to investigate the interfacial properties of CO2-H2O and CO2-brine mixtures. Simulations of the CO2-water system have shown that the strength of the short-range part of the CO2-H2O interaction controls the adsorption of CO2 to the fluid-fluid interface, which in turn controls the pressure-dependence of the IFT. Adsorption is apparent in Figure 1, where carbon dioxide density is enhanced adjacent to the water phase. Our simulations also reproduce the experimentally-observed increase in IFT with increasing salinity.

Co2 water

Figure 1. Snapshot of a molecular dynamics simulation cell containing CO2 (EPM2 model) and H2O (SPC/E model) molecules equilibrated at 70 °C and 20 MPa. Oxygen, carbon, and hydrogen atoms are represented by red, cyan, and gray spheres respectively. (Bottom) Time-averaged supercritical CO2 and liquid water density of the same simulation as a function of distance along the cell. The CO2 density is significantly enhanced near the interface.


Magnetic Ordering in Mackinawite

Mackinawite (FeS) is a layer type mineral (P4/nmm space group) composed of sheets of edge-sharing FeS4 tetrahedra stacked along the c-axis and held together through van der Waals forces (Fig. 2a). In nature, this mineral is typically the first solid phase to be precipitated during iron sulfide mineralization mediated by sulfate-reducing bacteria. Mackinawite is known to influence the mobility of environmentally-important trace elements through sorption and redox processes, and it figures importantly in the well-known iron-sulfur hypothesis about the origin of life. Its Mössbauer spectrum exhibits one singlet without quadrupole and hyperfine coupling, indicating that there is no magnetic ordering from room temperature down to 4.2 K. Our DFT calculations, however, show a magnetic moment on the Fe ions and that the ground state of mackinawite should be both metallic and magnetic (Fig. 2b). A similar conflict between Mössbauer spectra and DFT calculations has been reported for Fe-based high-temperature superconductors, such as FeSe (Tc = 12 K), LaO0.89F0.11FeAs (Tc = 26 K), and CeO0.85FeAs (Tc = 46.5 K), which are isostructural with mackinawite and for which it is known that magnetism is suppressed by strong magnetic moment (spin) fluctuations.

If Fe ions in metallic systems like mackinawite have a magnetic moment, multiplet splitting in the Fe 3s core-level photoemission spectrum should occur due to the exchange interaction of core 3s electrons with the net spin of the 3d electrons. The Fe 3s core-level photoemission spectrum of mackinawite, obtained at the Advanced Light Source (Fig. 2c), indeed showed exchange-energy splitting (ΔE3s) consistent with a magnetic moment on the Fe ions. All of these computational and experimental results indicate the existence of strong itinerant spin fluctuations in mackinawite. Because these spin fluctuations are believed to the main mediators of electron-pairing in the Fe-based superconductors, we conjecture that mackinawite may in fact be a superconductor.

Mackinawite structure

Figure 2. (a) Mackinawite structure (Purple = Fe; Orange = S, a = b = 3.674 Å, c = 5.033 Å). (b) DFT-calculated total energy of mackinawite per formula unit vs. the lattice parameter in different magnetic orderings (FM: ferromagnetic; NM: non-magnetic; AFM: antiferromagnetic order). In the AFM ordering, three types were considered. The vertical line represents the experimental lattice parameter. (c) Experimental Fe-3s core-level photoemission spectrum of mackinawite. ΔE3s is the multiplet splitting energy, equal to 2.9 eV.
Back to top

Structure and Dynamics of Water and Solutes Near the Surfaces of Smectite Clay Minerals

The plane of oxygen atoms on the cleavage surface of a 2:1 clay mineral is called a siloxane surface. This plane is characterized by hexagonal symmetry among its constituent oxygen ions. Associated with the siloxane surface is a roughly hexagonal cavity, formed by the bases (triads of oxygen ions) of six corner-sharing silica tetrahedra. This cavity has a diameter of about 0.26 nm. If there is no structural charge localized near a cavity, it can bind water molecules attracted to the proton in the structural OH group nestled inside it, or it can bind hydrophobic molecules, such as methane, that are attracted to its oxygen atoms in preference to those in more polar water molecules. If there is structural charge localized near a cavity, then interlayer cations and water molecules both are attracted to this charge and may compete. In aqueous solutions that contain several types of cations, competition of cations (such as Na+, Ca2+, and CaCl+ ion pairs in Fig. 3) for the structural charge sites results in the well-known cation exchange capacity of clay minerals. The first few statistical monolayers of water on siloxane surfaces (Fig. 3) are structured by their propensity to solvate adsorbed cations while forming hydrogen bonds with surface O atoms.

Density map
Figure 3. Density map of Na(dark blue), Ca2+ (light blue) and Cl- (yellow) ions during a 200 ps interval of a molecular dynamics simulation of a 1.6 M Na-Ca-Cl brine in a 6 nm pore between two smectite clay particles (the clay particle on the left side has a periodic image, not shown, on the right side).  A snapshot of the location of water O atoms in the first two statistical water monolayers on the surface (red spheres) indicates that Na+ and Ca2+ adsorb as outer-sphere surface complexes.

Diffusion in Clay Barriers

Natural and engineered clay barriers are widely used for the isolation of landfills and polluted sites.  These materials are under consideration for similar use as barriers in the subsurface storage of high-level radioactive waste and in the geologic sequestration of CO2.  In pure Na-montmorillonite compacted to > 1 kg dm-3, more than half of the pore space is located in < 1 nm-thick smectite interlayer nanopores.  Extant models of the sealing properties of clay barriers assume that water located in these nanopores diffuses as rapidly as bulk liquid water (i.e., they assume that water diffusion is slowed down only by "geometric" effects such as tortuosity).  Our research on the diffusion of water and ions in clay barriers showed that continuum-scale experimental data can be modeled accurately by considering interlayer nanopores (Fig. 4a) as a distinct "compartment" of the pore space where water and solutes diffuse significantly more slowly than in larger pores.  Our model calculations show that diffusion measurements in cm-scale samples of pure Na-montmorillonite are remarkably consistent with MD simulations predictions of nm-scale diffusion in individual Na-montmorillonite interlayer nanopores (Fig. 4b).

fig4 a & b
Figure 4. (a) MD simulation snapshot of water (O atoms in red, H atoms in gray) and sodium ions (blue) in the two-layer hydrate of Na-montmorillonite (clay Si, Al, Mg atoms are shown in yellow, green, and blue).  (b) Semilog plot of the mean principal value of the apparent diffusion coefficient tensor Da of water tracers, normalized to the self-diffusion coefficient of bulk liquid water (D0) vs. partial montmorillonite dry density (the mass of montmorillonite per volume of montmorillonite and pore space).  Model predictions assume that water diffuses more slowly in smectite interlayer nanopores than in larger pores, as confirmed by MD simulations.

Solute Isotope Fractionation by Diffusion in Liquid Water and Water-Filled Nanopores

The transport of solutes through porous media and their reactivity at nanoparticle surfaces can be kinetically controlled by molecular diffusion in bulk water or near hydrated mineral surfaces. Molecular diffusion in water can also result in significant isotopic fractionation, allowing the inference of geochemical fluxes from observed isotopic distributions. Essential to the interpretation of these distributions is molecular-scale information about the magnitude and mechanisms of solute diffusion and diffusive fractionation. However, the isotopic mass (m) dependence of solute diffusion coefficients (D) remains poorly known despite nearly a century of research on diffusion in liquid water. In 2006, experiments by Frank Richter (U. of Chicago) showed conclusively that lithium and chloride isotopes are fractionated by diffusion in liquid water, whereas magnesium isotopes are not.

In a series of MD simulations designed to provide insight into the results of Frank Richter, we investigated the diffusion of a range of real and hypothetical isotopes (m = 2 to 100 Da) in ambient liquid water (Fig. 5a,b). Our simulations confirmed the experimental findings of the Richter group and predicted the fractionation of potassium isotopes that was later confirmed experimentally by Frank Richter and John Christensen. Our results showed that D has an inverse power-law mass dependence (Fig. 5c) with a power-law exponent in the range between 0 and 0.2 for noble gases, alkali metals, alkaline earth metals, and chloride. Analysis of our simulation results showed that the mass-dependence of D is greatest for solutes whose motions are most weakly coupled to solvent hydrodynamic modes; thus, the power-law exponent of the mass dependence of D increases with the inverse of the average residence time of water molecules in the first solvation shell of the diffusing solute. Our findings for the noble gases—important geochemical indicators of hydrologic transport processes and climatic change—illustrate a key strength of molecular simulations, because aquated noble gases are arduous to probe experimentally but can be routinely probed by MD simulation. A groundwater hydrology application of our MD simulation results is described below (Fig. 6).

Fig5 a-b-c
Figure 5. (a) Snapshot of a MD simulation cell containing one lithium ion (green) and 215 water molecules (O atoms in red, H atoms in gray). (b) Molecular dynamics snapshot of a 40Ar atom (in green) and water molecules located in its first solvation shell (i.e., at Ar-O distances < 0.52 nm). First-shell water molecules preferentially adopt a "straddling" configuration by orienting either a lone pair of electrons (highlighted water molecule) or a hydrogen atom directly away from the solute. (c) Diffusion coefficients of four real or hypothetical isotopes of He, Ne, Ar, and Xe, plotted as a function of the logarithm of isotopic mass. Each symbol shows a D value calculated from a 2 ns "block" of each 8 ns simulation. For all solutes, the results are consistent with an inverse power-law relation between D and m.

Fig6 a-b-c

Figure 6. Illustration of noble gas isotope fractionation by diffusion in liquid water. (a) Schematic view of noble gas fluxes in a Swiss lake where methane ebullition is thought to occur in sediments deposited during the last ~150 years (figures by Matthias Brennwald). (b) Argon concentration in the sediment pore water, normalized to Ar in the overlying lake water and plotted as a function of depth in the sediments, shows a depletion in young sediments that may result from noble gas "stripping" by methane bubbles (diamonds and solid lines: experimental data and model predictions by Brennwald and coworkers; the model accounts for methane ebullition, sediment burial, and diffusion). Right: Faster upward diffusion of light isotopes causes 36Ar/40Ar enrichment near the top of the lake sediments; this is correctly predicted (solid lines) by combining the groundwater hydrology model of Brennwald and coworkers with our MD simulation estimates of the mass dependence of DAr
Back to top

Mn-Vacancy-Induced Photoconductivity of Layer Type Manganese Oxides

Manganese(IV) oxides are known to impact a broad range of biogeochemical processes, mainly through their high capacity for metal sorption and their facile oxidation of organic and inorganic compounds. Mn oxides found in weathering environments and natural waters are produced mainly by bacteria and take on layer structures, i.e., stacks of sheets of edge-sharing MnO6 octahedra. An important structural characteristic of these oxides is the presence of Mn(IV) cation vacancies whose charge deficit is typically compensated by metal cations or protons. These vacancies have long been identified as strong adsorption sites for metals, but they also may play an important role in Mn redox biogeochemistry, particularly photo-induced redox reactions. Because detailed electronic band structure is the key to understanding photo-induced redox reactions, Mn(IV) oxides both vacancy-free and containing cation vacancies charge-compensated by protons were investigated using ab initio quantum mechanics simulations as realized in the code CASTEP.

The ab initio study showed that a Mn(IV) vacancy reduces the band gap energy (Fig. 7a) and separates photo-excited electrons and holes (Fig. 7b). A reduction in band-gap energy generates more pairs of electrons and holes upon illumination, and the distinct separation of charge carriers enhances their transfer before loss by recombination. Therefore, Mn(IV) vacancies enhance effectively the photoconductivity of layer type MnO2, facilitating photo-redox reactions between the mineral and inorganic or organic compounds.

Recent studies in materials science indicate that synthetic layer type Mn(IV) oxide nanoparticles with cation vacancies like those found in the biogenic MnO2minerals are semiconductors that produce photocurrent under visible light stimulation, thus making them very attractive for applications in energy storage, solar cell fabrication, and catalysis. Our prediction of effective band gap energy reduction by cation vacancies indicates that photocurrent production by these layer type MnO2 nanoparticles can be optimized by the control of vacancy concentrations during laboratory synthesis.

Fig 7

Figure 7. Calculated electronic properties of vacancy-free MnO2 (left) and MnO2 containing 12.5 or 3.3 % Mn(IV) vacancies charge-compensated by protons (right). (a) Band structures, showing reduction of the band gap energy (denoted by a red double arrow) from 1.3 to 0.3 eV. (Solid lines denote bands for spin-up electrons while dashed lines denote bands for spin-down electrons.) (b) Charge distributions for valence-band maximum hole states (orange) and conduction-band minimum electron states (blue), showing that a cation vacancy compensated by protons (H) effectively separates photo-induced charge carriers against recombination.


Surface Complexes of Zinc in Birnessite Interlayer

Biogeochemical cycling of zinc (Zn) is strongly influenced by sorption on hexagonal birnessite. Spectroscopic studies have identified that Zn forms both tetrahedral (ZnIV) and octahedral (ZnVI) triple-corner-sharing surface complexes (TCS) at Mn(IV) vacancy sites in hexagonal birnessite (Fig. 8). The octahedral complex is expected because it is similar to that of Zn in the Mn oxide mineral, chalcophanite (ZnMn3O7·3H2O). However, the reason for the occurrence of the four-coordinate Zn surface species remains unclear. We examined this issue using spin-polarized DFT.

Total energy, magnetic moments of vacancy O, and electron-overlap populations between Zn and vacancy O all have proved that ZnIV-TCS is stable in birnessite without a need for Mn(III) substitution in the octahedral sheet and that it is more effective in reducing undersaturation of surface O at a Mn vacancy than is ZnVI-TCS. Geometry-optimization of a hypothetical monohydrate mineral, ZnMn3O7·H2O, which has a similar structure to chalcophanite but contains only tetrahedral Zn, and comparison with chalcophanite structure suggest that occurrence of tetrahedal Zn along with octahedral Zn in hexagonal birnessite is more related to the Zn-birnessite stability mediated by interactions of Zn ions (e.g., through H-bonds) with neighboring Zn ions or with adjacent Mn octahedral sheets and stacking modes of the sheets, rather than to the existence of Mn(III) substitution or local Zn bond strength difference.

Fig 8 a-b-c
Figure 8. Red: O; White: H; Teal: Zn; Purple octahedra: Mn(IV)O6 octahedra. Structure of Zn complexes in (a) chalcophanite (ZnMn3O7·3H2O), viewed perpendicularly to the c-axis. Isolated surface complexes of Zn at a Mn(IV) vacancy of birnessite interlayers: (b) tetrahedral Zn triple-corner-sharing (TCS) complex, with one H2O directly coordinated to Zn and two extra H2O H-bonding to the complex; (c) octahedral Zn TCS complex, with three H2O coordinated to Zn.

Nickel Partitioning at Mn(IV) Vacancies in Birnessite

Besides forming TCS surface complexes, Ni is also able to enter a cation vacancy, in effect substituting for the missing Mn (Ni-Sub). The ratio of Ni-Sub/Ni-TCS is observed to increase with pH and time of equilibration, which can be attributed to a pH-dependent activation energy for the TCS → Sub transformation, but this hypothesis is challenging to test by experiment. Our DFT study of the reaction path (Fig.9) found that: (1) an energy barrier exists between Ni-TCS and Ni-Sub and (2) the height of the barrier increases with the number of protons bound at the vacancy site (i.e., it is pH-dependent). The activation energy explains the predominance of Ni-TCS species over Ni-Sub in short-term experiments. Note also in Fig. 4 the increasing depth of the Ni-Sub energy minimum relative to Ni-TCS as pH increases.
Fig 9 a-b
Figure 9. Reaction path in the Ni-TCS to Ni-Sub transformation at (a) low pH and (b) high pH. The reaction coordinate marked 0.0 corresponds to Ni-TCS while the reaction coordinate marked 1.0 corresponds to Ni-Sub. Low pH and high pH conditions were represented by differing numbers of protons bound at the vacancy site. The energy (ΔE) is expressed relative to the total energy of Ni-TCS (reaction coordinate 0.0).


Proton adsorption on smectite clay minerals

Smectite clay minerals (flake-shaped nanoparticles) are ubiquitous in terrestrial weathering environments where they strongly influence the hydrology and adsorption properties of geologic media. In addition to their well-known permanent structural charge (circa -1 mol kg-1), these minerals carry a variable (pH-dependent) charge as a result of proton adsorption-desorption reactions on their edge surfaces. The protonation of smectite edge surfaces controls smectite dissolution rates, the rheology of smectite-water mixtures, and the adsorption of anions and transition metals on smectite clay minerals. Proton adsorption on smectites is commonly modeled by combining a surface complexation model of oxide-type edge sites (=SiOH, =AlOH2, etc.) with a cation-exchange model of structural charge sites (=X-). Studies of the acid-base surface chemistry of the smectite clay mineral montmorillonite yield conflicting values (ranging from 4.5 to 9.5) for the point of zero net proton charge (p.z.n.p.c.) of montmorillonite edge surfaces. Our research showed that this discrepancy results in part from the erroneous assumption that the net proton surface charge prior to titration is equal to zero. In fact, the best available data on the acid-base titration of Na-montmorillonite can be reconciled (Fig. 10) by discarding this assumption and accounting for the "spillover" of electrostatic potential between basal and edge surfaces.
Fig 10 a-b
Figure 10. Best available data on the acid-base titration on Na-montmorillonite at several ionic strengths, previously interpreted as indicating widely different p.z.n.p.c. values for Na-montmorillonite edge surfaces. A novel simulation scheme (involving a surface chemistry model to simulate both pretreatment and titration) shows that these two data sets can be described using a single model of proton adsorption on montmorillonite.
Back to top