Publications

Abstract not provided.

We consider the influence of ambient groundwater flow on the migration of DNAPL within a fracture network. In context of a modified invasion percolation (MIP) growth algorithm, we formulate a mechanistic model that includes capillary and gravity forces as well as viscous forces within the DNAPL and the ambient groundwater. The MIP model is verified against laboratory experiments, which show good agreement in DNAPL migration path through a two-dimensional fracture network. The results of both simulations and laboratory experiments suggest that ambient groundwater flow can be a significant factor controlling DNAPL migration path, velocity, and channeling pattern in a fracture network.

Abstract not provided.

Abstract not provided.

Fingering, nonmonotonicity, fragmentation, and pulsation within gravity/buoyant destabilized two-phase/unsaturated flow systems has been widely observed with examples in homogeneous to heterogeneous porous media, in single fractures to fracture networks, and for both wetting and nonwetting invasion. To model this phenomena, we consider a mechanistic approach based on forms of modified invasion percolation (MIP) that include gravity, the influence of the local interfacial curvature along the phase-phase interface, and the simultaneous invasion and reinvasion of both wetting and nonwetting fluids. We present example simulations and compare them to experimental data for three very different situations: (1) downward gravity-driven fingering of water into a dry, homogeneous, water-wettable, porous medium; (2) upward buoyancy-driven migration of gas within a water saturated, heterogeneous, water-wettable, porous medium; and (3) downward gravity-driven fingering of water into a dry, water-wettable, rough-walled fracture.

We consider the use of a hypodiffusive governing equation (HDE) for the porous-continuum modeling of gravity-driven fingers (GDF) as occur in initially dry, highly nonlinear, and hysteretic porous media. In addition to the capillary and gravity terms within the traditional Richards equation, the HDE contains a hypodiffusive term that models an experimentally observed hold-back-pile-up (HBPU) effect and thus imparts nonmonotonicity at the wetting front. In its dimensionless form the HDE contains the dimensionless hypodiffusion number, NHD. As NHD increases, one-dimensional (1D) numerical solutions transition from monotonic to nonmonotonic. Considering the experimentally observed controls on GDF occurrence, as either the initial moisture content and applied flux increase or the material nonlinearity decreases, solutions undergo the required transition back to monotonic. Additional tests for horizontal imbibition and capillary rise show the HDE to yield the required monotonie response but display sharper fronts for NHD > 0. Finally, two-dimensional (2D) numerical solutions illustrate that in parameter space where the 1D HDE yields nonmonotonicity, in 2D it forms nonmonotonic GDF.

A simple experiment finds that fracture intersections can act to integrate unsaturated flows, such that regular, low flows entering the intersection from above are transformed into large, less frequent pulses below. At low flows, our simple intersection forms two capillary barriers. Water from above pools at the intersection until sufficient pressure builds to breach the barriers and discharge stored fluid. The barriers then reform and the process is repeated. At low flows, the volume discharged from the intersection remains relatively uniform across a range of flow rates. At higher flows, discharge volume is highly variable, and a viscous stabilized non-pulsating regime occurs at the highest flow rate that we considered.

Abstract not provided.

Methods to determine unsaturated hydraulic properties can exhibit random and nonunique behavior. We assess the causes for these behaviors by visualizing microscale phase displacement processes that occur during equilibrium retention and transient outflow experiments. For both types of experiments we observe the drainage process to be composed of a mixture of fast air fingering and slower air back-filling. The influence of each of these microscale processes is controlled by a combination of the size and the speed of the applied boundary step, the initial saturation and its structure, and small-scale heterogeneities. Because the mixture of these microscale processes yields macroscale effective behavior, measured unsaturated flow properties are also a function of these controls. Such results suggest limitations on the current definitions and uniqueness of unsaturated hydraulic properties.

Hydraulic property measurements often rely on non-linear inversion models whose errors vary between samples. In non-linear physical measurement systems, bias can be directly quantified and removed using calibration standards. In hydrologic systems, field calibration is often infeasible and bias must be quantified indirectly. We use a Monte Carlo error analysis to indirectly quantify spatial bias in the saturated hydraulic conductivity, K{sub s}, and the exponential relative permeability parameter, {alpha}, estimated using a tension infiltrometer. Two types of observation error are considered, along with one inversion-model error resulting from poor contact between the instrument and the medium. Estimates of spatial statistics, including the mean, variance, and variogram-model parameters, show significant bias across a parameter space representative of poorly- to well-sorted silty sand to very coarse sand. When only observation errors are present, spatial statistics for both parameters are best estimated in materials with high hydraulic conductivity, like very coarse sand. When simple contact errors are included, the nature of the bias changes dramatically. Spatial statistics are poorly estimated, even in highly conductive materials. Conditions that permit accurate estimation of the statistics for one of the parameters prevent accurate estimation for the other; accurate regions for the two parameters do not overlap in parameter space. False cross-correlation between estimated parameters is created because estimates of K{sub s} also depend on estimates of {alpha} and both parameters are estimated from the same data.

Double-diffusive finger convection is a hydrodynamic instability that can occur when two components with different diffusivities are oppositely stratified with respect to the fluid density gradient as a critical condition is exceeded. Laboratory experiments were designed using sodium chloride and sucrose solutions in a Hele-Shaw cell. A high resolution, full field, light transmission technique was used to study the development of the instability. The initial buoyancy ratio (R{sub p}), which is a ratio of fluid density contributions by the two solutes, was varied systematically in the experiments so that the range of parameter space spanned conditions that were nearly stable (R{sub p} = 2.8) to those that were moderately unstable (R{sub p} = 1.4). In systems of low R{sub p}, fingers develop within several minutes, merge with adjacent fingers, form conduits, and stall before newer-generated fingers travel through the conduits and continue the process. Solute fluxes in low R{sub p} systems quickly reach steady state and are on the order of 10{sup {minus}6} m{sup 2} sec{sup {minus}1}. In the higher R{sub p} experiments, fingers are slower to evolve and do not interact as dynamically as in the lower R{sub p} systems. Our experiment with initial R{sub p} = 2.8 exhibited flux on the order of that expected for a similar diffusive system (i.e., 10{sup {minus}7} m{sup 2} sec{sup {minus}1}), although the structures were very different than the pattern of transport expected in a diffusing system. Mass flux decayed as t{sup 1/2} in two experiments each with initial R{sub p} = 2.4 and 2.8.

An experimental investigation was conducted to study double-diffusive finger convection in a Hele-Shaw cell by layering a sucrose solution over a more-dense sodium chloride (NaCl) solution. The solutal Rayleigh numbers were on the order of 60,000, based upon the height of the cell (25 cm), and the buoyancy ratio was 1.2. A full-field light transmission technique was used to measure a dye tracer dissolved in the NaCl solution. They analyze the concentration fields to yield the temporal evolution of length scales associated with the vertical and horizontal finger structure as well as the mass flux. These measures show a rapid progression through two early stages to a mature stage and finally a rundown period where mass flux decays rapidly. The data are useful for the development and evaluation of numerical simulators designed to model diffusion and convection of multiple components in porous media. The results are useful for correct formulation at both the process scale (the scale of the experiment) and effective scale (where the lab-scale processes are averaged-up to produce averaged parameters). A fundamental understanding of the fine-scale dynamics of double-diffusive finger convection is necessary in order to successfully parameterize large-scale systems.

The authors performed a simple experiment to elucidate phase structure within a pervasively fractured welded tuff. Dyed water was infiltrated from a surface pond over a 36 minute period while a geophysical array monitored the wetted region within vertical planes directly beneath. They then excavated the rock mass to a depth of {approximately}5 m and mapped the fracture network and extent of dye staining in a series of horizontal pavements. Near the pond the network was fully stained. Below, the phase structure immediately expanded and with depth, the structure became fragmented and complicated exhibiting evidence of preferential flow, fingers, irregular wetting patterns, and varied behavior at fracture intersections. Limited transient geophysical data suggested that strong vertical pathways form first followed by increased horizontal expansion and connection within the network. These rapid pathways are also the first to drain. Estimates also suggest that the excavation captured from {approximately}10% to 1% or less of the volume of rock interrogated by the infiltration slug and thus the penetration depth could have been quite large.

The authors reconceptualize macro modified invasion percolation (MMIP) at the near pore (NP) scale and apply it to simulate the non-wetting phase invasion experiments of Glass et al [in review] conducted in macro-heterogeneous porous media. For experiments where viscous forces were non-negligible, they redefine the total pore filling pressure to include viscous losses within the invading phase as well as the viscous influence to decrease randomness imposed by capillary forces at the front. NP-MMIP exhibits the complex invasion order seen experimentally with characteristic alternations between periods of gravity stabilized and destabilized invasion growth controlled by capillary barriers. The breaching of these barriers and subsequent pore scale fingering of the non-wetting phase is represented extremely well as is the saturation field evolution, and total volume invaded.

Surfactant-enhanced aquifer remediation is an emerging technology for aquifers contaminated with nonaqueous phase liquids (NAPLs). A two-dimensional micromodel and image capture system were applied to observe NAPL mobilization and solubilization phenomena. In each experiment, a common residual NAPL field was established, followed by a series of mobilization and solubilization experiments. Mobilization floods included pure water floods with variable flow rates and surfactant floods with variations in surfactant formulations. At relatively low capillary numbers (N{sub ca}<10{sup {minus}3}), the surfactant mobilization floods resulted in higher NAPL saturations than for the pure water flood, for similar N{sub ca}.These differences in macroscopic saturations are explained by differences in micro-scale mobilization processes. Solubilization of the residual NAPL remaining after the mobilization stage was dominated by the formation of dissolution fingers, which produced nonequilibrium NAPL solubilization. A macroemulsion phase also as observed to form spontaneously and persist during the solubilization stage of the experiments.

The authors consider the ability of the numerical solution of Richards equation to model gravity-driven fingers. Although gravity-driven fingers can be easily simulated using a partial downwind averaging method, they find the fingers are purely artificial, generated by the combined effects of truncation error induced oscillations and capillary hysteresis. Since Richards equation can only yield a monotonic solution for standard constitutive relations and constant flux boundary conditions, it is not the valid governing equation to model gravity-driven fingers, and therefore is also suspect for unsaturated flow in initially dry, highly nonlinear, and hysteretic media where these fingers occur. However, analysis of truncation error at the wetting front for the partial downwind method suggests the required mathematical behavior of a more comprehensive and physically based modeling approach for this region of parameter space.

Recent work demonstrates that phase displacements within horizontal fractures large with respect to the spatial correlation length of the aperture field lead to a satiated condition that constrains the relative permeability to be less than one. The authors use effective media theory to develop a conceptual model for satiated relative permeability, then compare predictions to existing experimental measurements, and numerical solutions of the Reynolds equation on the measured aperture field within the flowing phase. The close agreement among all results and data show that for the experiments considered here, in-plane tortuosity induced by the entrapped phase is the dominant factor controlling satiated relative permeability. They also find that for this data set, each factor in the conceptual model displays an approximate power law dependence on the satiated saturation of the fracture.

Fracture transmissivity and detailed aperture fields are measured in analog fractures specially designed to evaluate the utility of the Reynolds equation. The authors employ a light transmission technique with well-defined accuracy ({approximately}1% error) to measure aperture fields at high spatial resolution ({approximately}0.015 cm). A Hele-Shaw cell is used to confirm the approach by demonstrating agreement between experimental transmissivity, simulated transmissivity on the measured aperture field, and the parallel plate law. In the two rough-walled analog fractures considered, the discrepancy between the experimental and numerical estimates of fracture transmissivity was sufficiently large ({approximately} 22--47%) to exclude numerical and experimental errors (< 2%)as a source. They conclude that the three-dimensional character of the flow field is important for fully describing fluid flow in the two rough-walled fractures considered, and that the approach of depth averaging inherent in the formulation of the Reynolds equation is inadequate. They also explore the effects of spatial resolution, aperture measurement technique, and alternative definitions for link transmissivities in the finite-difference formulation, including some that contain corrections for tortuosity perpendicular to the mean fracture plane and Stokes flow. Various formulations for link transmissivity are shown to converge at high resolution ({approximately} 1/5 the spatial correlation length) in the smoothly varying fracture. At coarser resolutions, the solution becomes increasingly sensitive to definition of link transmissivity and measurement technique. Aperture measurements that integrate over individual grid blocks were less sensitive to measurement scale and definition of link transmissivity than point sampling techniques.

The authors develop and evaluate a modified invasion percolation (MIP) model for quasi-static immiscible displacement in horizontal fractures. The effects of contact angle, local aperture field geometry, and local in-plane interracial curvature between phases are included in the calculation of invasion pressure for individual sites in a discretized aperture field. This pressure controls the choice of which site is invaded during the displacement process and hence the growth of phase saturation structure within the fracture. To focus on the influence of local in-plane curvature on phase invasion structure, they formulate a simplified nondimensional pressure equation containing a dimensionless curvature number (C) that weighs the relative importance of in-plane curvature and aperture-induced curvature. Through systematic variation of C, they find in-plane interracial curvature to greatly affect the phase invasion structure. As C is increased from zero, phase invasion fronts transition from highly complicated (IP results) to microscopically smooth. In addition, measurements of fracture phase saturations and entrapped cluster statistics (number, maximum size, structural complication) show differential response between wetting and nonwetting invasion with respect to C that is independent of contact angle hysteresis. Comparison to experimental data available at this time substantiates predicted behavior.

Dispersion of solutes in a variable aperture fracture results from a combination of molecular diffusion and velocity variations in both the plane of the fracture (macrodispersion) and across the fracture aperture (Taylor dispersion). We use a combination of physical experiments and computational simulations to test a theoretical model in which the effective longitudinal dispersion coefficient D(L) is expressed as a sum of the contributions of these three dispersive mechanisms. The combined influence of Taylor dispersion and macrodispersion results in a nonlinear dependence of D(L) on the Peclet number (Pe = V/D(m), where V is the mean solute velocity,is the mean aperture, and D(m) is the molecular diffusion coefficient). Three distinct dispersion regimes become evident: For small Pe (Pe << 1), molecular diffusion dominates resulting in D(L) proportional to Pe0; for intermediate Pe, macrodispersion dominates (D(L) proportional to Pe); and for large Pe, Taylor dispersion dominates (D(L) proportional to Pe2). The Pe range corresponding to these different regimes is controlled by the statistics of the aperture field. In particular, the upper limit of Pe corresponding to the macrodispersion regime increases as the macrodispersivity increases. Physical experiments in an analog, rough-walled fracture confirm the nonlinear Pe dependence of D(L) predicted by the theoretical model. However, the theoretical model underestimates the magnitude of D(L). Computational simulations, using a particle-tracking algorithm that incorporates all three dispersive mechanisms, agree very closely with the theoretical model predictions. The close agreement between the theoretical model and computational simulations is largely because, in both cases, the Reynolds equation describes the flow field in the fracture. The discrepancy between theoretical model predictions and D(L) estimated from the physical experiments appears to be largely, due to deviations from the local cubic law assumed by the Reynolds equation.

The authors designed and conducted experiments in a heterogeneous sand pack where gravity-destabilized nonwetting phase invasion (CO{sub 2} and TCE) could be recorded using high resolution light transmission methods. The heterogeneity structure was designed to be reminiscent of fluvial channel lag cut-and-fill architecture and contain a series of capillary barriers. As invasion progressed, nonwetting phase structure developed a series of fingers and pools; behind the growing front they found nonwetting phase saturation to pulsate in certain regions when viscous forces were low. Through a scale analysis, they derive a series of length scales that describe finger diameter, pool height and width, and regions where pulsation occurs within a heterogeneous porous medium. In all cases, they find that the intrinsic pore scale nature of the invasion process and resulting structure must be incorporated into the analysis to explain experimental results. The authors propose a simple macro-scale structural growth model that assembles length scales for sub-structures to delineate nonwetting phase migration from a source into a heterogeneous domain. For such a model applied at the field scale for DNAPL migration, they expect capillary and gravity forces within the complex subsurface lithology to play the primary roles with viscous forces forming a perturbation on the inviscid phase structure.