Search

When a liquid touches a solid surface, it spreads to minimize the system's energy. The classic thin-film model describes the spreading as an interplay between gravity, capillarity, and viscous forces, but it cannot see an end to this process as it does not account for the nonhydrodynamic liquid-solid interactions. While these interactions are important only close to the contact line, where the liquid, solid, and gas meet, they have macroscopic implications: in the partial-wetting regime, a liquid puddle ultimately stops spreading. We show that by incorporating these intermolecular interactions, the free energy of the system at equilibrium can be cast in a Cahn-Hilliard framework with a height-dependent interfacial tension. Using this free energy, we derive a mesoscopic thin-film model that describes the statics and dynamics of liquid spreading in the partial-wetting regime. The height dependence of the interfacial tension introduces a localized apparent slip in the contact-line region and leads to compactly supported spreading states. In our model, the contact-line dynamics emerge naturally as part of the solution and are therefore nonlocally coupled to the bulk flow. Surprisingly, we find that even in the gravity-dominated regime, the dynamic contact angle follows the Cox-Voinov law.

Flow and transport through fractured geologic media often leads to anomalous (non-Fickian) transport behavior, the origin of which remains a matter of debate: whether it arises from variability in fracture permeability (velocity distribution), connectedness in the flow paths through fractures (velocity correlation), or interaction between fractures and matrix. Here we show that this uncertainty of distribution- versus correlation-controlled transport can be resolved by combining convergent and push-pull tracer tests because flow reversibility is strongly dependent on velocity correlation, whereas late-time scaling of breakthrough curves is mainly controlled by velocity distribution. We build on this insight, and propose a Lagrangian statistical model that takes the form of a continuous time random walk (CTRW) with correlated particle velocities. In this framework, velocity distribution and velocity correlation are quantified by a Markov process of particle transition times that is characterized by a distribution function and a transition probability. Our transport model accurately captures the anomalous behavior in the breakthrough curves for both push-pull and convergent flow geometries, with the same set of parameters. Thus, the proposed correlated CTRW modeling approach provides a simple yet powerful framework for characterizing the impact of velocity distribution and correlation on transport in fractured media. © 2014. American Geophysical Union. All Rights Reserved., The authors thank Nicolas Lavenant and Rebecca Hochreutener for their help in setting up and running the field experiments. Funding for this work was provided by the INTERREG IV project CLIMAWAT, the national network of hydrogeological sites H+ (hplus.ore.fr). MD acknowledges the support of the European Research Council (ERC) through the project MHetScale (617511). RJ acknowledges the support of the US Department of Energy (grant DE-SC0003907), and a MISTI Global Seed Funds award.

Groundwater inflows into tunnels constructed in fractured bedrocks not only constitute an important factor controlling the rate of advancement in driving the tunnel but may pose potential hazards. Drawdowns caused by tunnel construction may also induce geotechnical and environmental impacts. Here we present a numerical methodology for the dynamic simulation of the hydrogeological transient conditions induced by the tunnel front advance. The methodology is based on the use of a Cauchy boundary condition at the points lying along the tunnel according to which water discharge, Q , is computed as the product of a leakage coefficient, α , and the head difference, ( H − h ), where H is the prescribed head at the tunnel wall and h is the hydraulic head in the fractured rock in the close vicinity of the tunnel. At a given position of the tunnel, α is zero until the tunnel reaches such position when it is assigned a positive value. The use of step-wise time functions for α allows an efficient and accurate simulation of the transient hydrogeological conditions at and around the tunnel during the excavation process. The methodology has been implemented in TRANMEF-3, a finite element computer code for groundwater flow in 3D fractured media developed at the University of A Coruna, Spain, and has been used to simulate the impact of a tunnel on the groundwater system at the Aspo island (Sweden). This tunnel was constructed to access an underground laboratory for research on radioactive waste disposal. The large amount of available data at this site provides a unique opportunity to test the performance of the numerical model and the proposed methodology for tunnel advance. With just minor calibration, the numerical model is able to reproduce accurately the measurements of inflows into the tunnel at several reaches and hydraulic heads at surface-drilled boreholes. These results obtained at the Aspo site lead us to conclude that accurate predictions of the transient hydrogeological responses induced by tunneling works in fractured bedrocks, can be achieved provided that a sound hydrogeological characterization of large-scale fracture zones is available.

The need to assess quantitatively the safety of waste repositories in deep geological media has fostered the development of efficient numerical models of groundwater flow and contaminant transport in fractured media. These models usually account for water flow through fracture zones embedded in a 3D rock matrix continuum. The first formulation of fractures in groundwater flow finite element models was proposed by Kiraly, and later revisited and generalized by Perrochet. From a mathematical viewpoint, fractures can be considered as m-dimensional manifolds in an n-dimensional Euclidean space (m⩽n). The key step of this formulation lies in an expression relating the hypersurface element dSm to the infinitesimal local co-ordinates dξi(i=1,…,m). Here we present a novel proof for this relation using a different approach to that of Perrochet, and explore the efficiency and accuracy of the formulation. It is shown that the aforementioned relation leads to a general and compact formulation which is not only applicable to elements of any dimension (e.g. 1D, 2D and 3D elements in a 3D domain), but also overcomes the cumbersome and case-specific calculations of traditional approaches. This formulation has been implemented in a versatile finite element program for modelling groundwater flow, solute transport and heat transport in porous and fractured media. The efficiency and accuracy of the proposed formulation has been analysed using a synthetic case dealing with flow and solute transport through a 2D fractured rock block. The proposed formulation, in which fractures are discretized by means of 1D elements is more efficient and accurate than the traditional finite element formulation of discretizing fractures by means of 2D elements. The capability of the proposed formulation to cope with complex systems is illustrated with a case study of groundwater flow induced by the construction of the access tunnel to an underground research laboratory in Aspo (Sweden). The numerical model is able to reproduce the observed records of water levels in boreholes and flow rates into the tunnel. Although the proposed formulation has been implemented and tested within the framework of groundwater flow and solute transport in fractured porous media, it should be of interest for other boundary value problems. Copyright © 2002 John Wiley & Sons, Ltd.

We study experimentally the miscible radial displacement of a more viscous fluid by a less viscous one in a horizontal Hele-Shaw cell. For the range of tested injection rates and viscosity ratios we observe two regimes for the evolution of the fluid-fluid interface. At early times the interface length increases linearly with time, which is typical of the Saffman-Taylor instability for this radial configuration. However, as time increases, the interface growth slows down and scales as ∼t(1/2), as one expects in a stable displacement, indicating that the overall flow instability has shut down. Surprisingly, the crossover time between these two regimes decreases with increasing injection rate. We propose a theoretical model that is consistent with our experimental results, explains the origin of this second regime, and predicts the scaling of the crossover time with injection rate and the mobility ratio. The key determinant of the observed scalings is the competition between advection and diffusion time scales at the displacement front, suggesting that our analysis can be applied to other interfacial-evolution problems such as the Rayleigh-Benard-Darcy instability.

Geological carbon dioxide (CO 2 ) sequestration entails capturing and injecting CO 2 into deep saline aquifers for long-term storage. The injected CO 2 partially dissolves in groundwater to form a mixture that is denser than the initial groundwater. The local increase in density triggers a gravitational instability at the boundary layer that further develops into columnar plumes of CO 2 -rich brine, a process that greatly accelerates solubility trapping of the CO 2 . Here, we investigate the pattern-formation aspects of convective mixing during geological CO 2 sequestration by means of high-resolution three-dimensional simulation. We find that the CO 2 concentration field self-organizes as a cellular network structure in the diffusive boundary layer at the top boundary. By studying the statistics of the cellular network, we identify various regimes of finger coarsening over time, the existence of a non-equilibrium stationary state, and a universal scaling of three-dimensional convective mixing.

The overall objective of this proposal was to develop tools for better understanding, modeling and risk assessment of CO2 permanence in geologic formations at the geologic basin scale.

Oil and natural gas from deep shale formations are transforming the United States economy and its energy outlook. Back in 2005, the US Energy Information Administration published projections of United States natural gas supply that stressed the need to develop an import infrastructure (1): by 2025, imports would account for almost one-third of United States consumption. When we compare those forecasts with the current ones to 2040 (Fig. 1) (2), it is inevitable to feel that a disruptive technology has emerged since. Natural gas consumption is expected to increase significantly over the next three decades, with strong demand growth from the electricity generation and industrial sectors. However, the United States will probably become a net exporter of gas before 2020, increasing domestic production by 44% over the projection period.

Traditionally, seismic interpretation is performed without any account of the flow behavior. Here, we present a methodology to characterize fractured geologic media by integrating flow and seismic data. The key element of the proposed approach is the identification of the intimate relation between acoustic and flow responses of a fractured reservoir through the fracture compliance. By means of synthetic models, we show that: (1) owing to the strong (but highly uncertain) dependence of fracture permeability on fracture compliance, the modeled flow response in a fractured reservoir is highly sensitive to the geophysical interpretation; and (2) by incorporating flow data (well pressures and production curves) into the inversion workflow, we can simultaneously reduce the error in the seismic interpretation and improve predictions of the reservoir flow dynamics.

Public policy and individual incentives determine the patterns of human mobility through transportation networks. In the event of a health emergency, the pursuit of maximum social or individual utility may lead to conflicting objectives in the routing strategies of network users. Individuals tend to avoid exposure so as to minimize the risk of contagion, whereas policymakers aim at coordinated behaviour that maximizes the social welfare. Here, we study agent-driven contagion dynamics through transportation networks, coupled to the adoption of either selfish- or policy-driven rerouting strategies. In analogy with the concept of price of anarchy in transportation networks subject to congestion, we show that maximizing individual utility leads to a loss of welfare for the social group, measured here by the total population infected after an epidemic outbreak.

We study mixing of two fluids of different viscosity in a microfluidic channel or porous medium. We show that the synergetic action of alternating injection and viscous fingering leads to a dramatic increase in mixing efficiency at high P\'eclet numbers. Based on observations from high-resolution simulations, we develop a theoretical model of mixing efficiency that combines a hyperbolic mixing model of the channelized region ahead and a mixing-dissipation model of the pseudosteady region behind. Our macroscopic model quantitatively reproduces the evolution of the average degree of mixing along the flow direction and can be used as a design tool to optimize mixing from viscous fingering in a microfluidic channel.

Convective mixing in porous media is triggered by a Rayleigh-B\'enard-type hydrodynamic instability as a result of an unstable density stratification of fluids. While convective mixing has been studied extensively, the fundamental behavior of the dissolution flux and its dependence on the system parameters are not yet well understood. Here, we show that the dissolution flux and the rate of fluid mixing are determined by the mean scalar dissipation rate. We use this theoretical result to provide computational evidence that the classical model of convective mixing in porous media exhibits, in the regime of high Rayleigh number, a dissolution flux that is constant and independent of the Rayleigh number. Our findings support the universal character of convective mixing and point to the need for alternative explanations for nonlinear scalings of the dissolution flux with the Rayleigh number, recently observed experimentally.

In a recent Perspective (1), Zoback and Gorelick argued that carbon capture and storage (CCS) is likely not a viable strategy for reducing CO2 emissions to the atmosphere. They argued that maps of earthquake epicenters portray earthquakes occurring almost everywhere, suggesting that Earth’s crust is near a critical state, so that increments in fluid pressure from injecting CO2 at 1 to 3 km depth will likely trigger earthquakes within the reservoir and caprock that would be expected to result in leakage of CO2 from the reservoirs to the surface.

We study the displacement of immiscible fluids in deformable, noncohesive granular media. Experimentally, we inject air into a thin bed of water-saturated glass beads and observe the invasion morphology. The control parameters are the injection rate, the bead size, and the confining stress. We identify three invasion regimes: capillary fingering, viscous fingering, and ''capillary fracturing,'' where capillary forces overcome frictional resistance and induce the opening of conduits. We derive two dimensionless numbers that govern the transition among the different regimes: a modified capillary number and a fracturing number. The experiments and analysis predict the emergence of fracturing in fine- grained media under low confining stress, a phenomenon that likely plays a fundamental role in many natural processes such as primary oil migration, methane venting from lake sediments, and the formation of desiccation cracks. DOI: 10.1103/PhysRevLett.108.264504 show that fracturing caused by elastic deformation of particles is the dominant mode of invasion for fine, soft particles under low confining stress. In this Letter, we study fracturing in a system of rigid particles. We show that frictional sliding, rather than elastic deformation, is the key mechanism controlling fracturing. We provide experi- mental evidence for three displacement regimes—capillary fingering, viscous fingering, and capillary fracturing—and derive two dimensionless groups that govern the transitions among these regimes. We inject air into a thin bed of water-saturated glass beads packed in a cylindrical cell. We prepare each granu- lar bed by pouring beads into the cell and then vibrating it to increase the packing homogeneity (25). There are three control parameters: the mean particle diameter, d, the air injection rate, q, and the confining weight, w. The confine- ment is applied by weights placed on a rigid disk that rests on top of the beads. To allow fluids (but not the particles) to leave the cell, the disk is made slightly smaller than the interior of the cell (inner diameter L), leading to a thin gap along the edge of the cell. We inject the fluids with a syringe pump at a fixed rate to the center of the cell, and image the invasion pattern during the experiment with a digital camera (Fig. 1). To analyze the images, we subtract the initial image prior to injection. For most experimental conditions, we perform repeatability experiments. The experiments exhibit three displacement regimes: viscous fingering, capillary fingering, and capillary fractur- ing. We identify the regimes based on the qualitative char- acteristics of the invasion pattern, as has been done previously (2,4,7,18,24). In viscous fingering, the air-water interface grows continuously andat several locations simul- taneously. The resulting pattern is radial and exhibits thin fingers and few trapped water clusters. The pattern appears more space-filling than in two-dimensional experiments

Drops and bubbles are nonspreading, local, compactly supported features. They are also equilibrium configurations in partial wetting phenomena. Yet, current macroscopic theories of capillary-dominated flow are unable to describe these systems. We propose a framework to model multiphase flow in porous media with nonspreading equilibrium configurations. We illustrate our approach with a one-dimensional model of two-phase flow in a capillary tube. Our model allows for the presence of compactons: nonspreading steady-state solutions in the absence of external forces. We show that local rate dependency is not needed to explain globally rate-dependent displacement patterns, and we interpret dynamic wetting transitions as the route from equilibrium, capillary-dominated configurations, towards viscous-dominated flow. Mathematically, these transitions are possible due to nonclassical shock solutions and the role of bistability and higher-order terms in our model.

In carbon capture and storage (CCS), CO 2 is captured at power plants and then injected underground into reservoirs like deep saline aquifers for long-term storage. While CCS may be critical for the continued use of fossil fuels in a carbon-constrained world, the deployment of CCS has been hindered by uncertainty in geologic storage capacities and sustainable injection rates, which has contributed to the absence of concerted government policy. Here, we clarify the potential of CCS to mitigate emissions in the United States by developing a storage-capacity supply curve that, unlike current large-scale capacity estimates, is derived from the fluid mechanics of CO 2 injection and trapping and incorporates injection-rate constraints. We show that storage supply is a dynamic quantity that grows with the duration of CCS, and we interpret the lifetime of CCS as the time for which the storage supply curve exceeds the storage demand curve from CO 2 production. We show that in the United States, if CO 2 production from power generation continues to rise at recent rates, then CCS can store enough CO 2 to stabilize emissions at current levels for at least 100 y. This result suggests that the large-scale implementation of CCS is a geologically viable climate-change mitigation option in the United States over the next century.

In this project we have sought to explain the co-existence of gas and hydrate phases in sediments within the gas hydrate stability zone. We have focused on the gas/brine interface at the scale of individual grains in the sediment. The capillary forces associated with a gas/brine interface play a dominant role in many processes that occur in the pores of sediments and sedimentary rocks. The mechanical forces associated with the same interface can lead to fracture initiation and propagation in hydrate-bearing sediments. Thus the unifying theme of the research reported here is that pore scale phenomena are key to understanding large scale phenomena in hydrate-bearing sediments whenever a free gas phase is present. Our analysis of pore-scale phenomena in this project has delineated three regimes that govern processes in which the gas phase pressure is increasing: fracturing, capillary fingering and viscous fingering. These regimes are characterized by different morphology of the region invaded by the gas. On the other hand when the gas phase pressure is decreasing, the corresponding regimes are capillary fingering and compaction. In this project, we studied all these regimes except compaction. Many processes of interest in hydrate-bearing sediments can be better understood when placed in the context of the appropriate regime. For example, hydrate formation in sub-permafrost sediments falls in the capillary fingering regime, whereas gas invasion into ocean sediments is likely to fall into the fracturing regime. Our research provides insight into the mechanisms by which gas reservoirs are converted to hydrate as the base of the gas hydrate stability zone descends through the reservoir. If the reservoir was no longer being charged, then variation in grain size distribution within the reservoir explain hydrate saturation profiles such as that at Mt. Elbert, where sand-rich intervals containing little hydrate are interspersed between intervals containing large hydrate saturations. Large volumes (of order one pore volume) of gaseous and aqueous phases must be transported into the gas hydrate stability zone. The driver for this transport is the pressure sink induced by a reduction in occupied pore volume that accompanies the formation of hydrate from gas and water. Pore-scale imbibition models and bed-scale multiphase flow models indicate that the rate-limiting step in converting gas to hydrate is the supply of water to the hydrate stability zone. Moreover, the water supply rate is controlled by capillarity-driven flux for conditions typical of the Alaska North Slope. A meter-scale laboratory experiment confirms that significant volumes of fluid phases move into the hydrate stability zone and that capillarity is essential for the water flux. The model shows that without capillarity-driven flux, large saturations of hydrate cannot form. The observations of thick zones of large saturation at Mallik and Mt Elbert thus suggest that the primary control on these systems is the rate of transport of gaseous and aqueous phases, driven by the pressure sink at the base of the gas hydrate stability zone. A key finding of our project is the elucidation of ?capillary fracturing? as a dominant gas transport mechanism in low-permeability media. We initially investigate this phenomenon by means of grain-scale simulations in which we extended a discrete element mechanics code (PFC, by Itasca) to incorporate the dynamics of first single-phase and then multiphase flow. A reductionist model on a square lattice allows us to determine some of the fundamental dependencies of the mode of gas invasion (capillary fingering, viscous fingering, and fracturing) on the parameters of the system. We then show that the morphology of the gas-invaded region exerts a fundamental control on the fabric of methane hydrate formation, and on the overpressures caused by methane hydrate dissociation. We demonstrate the existence of the different invasion regimes by means of controlled laboratory experiments in a radial cell. We collapse the behavior in the form of a phase diagram fully characterized by two dimensionless groups: a modified capillary number and a ?fracturing number? that reflects the balance between the pressure forces that act to open conduits in the granular pack, and frictional forces that resist it. We use all this small-scale knowledge to propose simple mechanistic models of gas migration and hydrate formation at the geologic bed scale. We propose that methane transport in lake and oceanic sediments is controlled by dynamic conduits, which dilate and release gas as the falling hydrostatic pressure reduces the effective stress below the tensile strength of the sediments. We test our model against a four-month record of hydrostatic load and methane flux in Upper Mystic Lake, Mass., USA, and show that it captures the complex episodicity of methane ebullition. Our quantitative conceptualization opens the door to integrated modeling of methane transport to constrain global methane release from lakes and other methane-rich sediment systems, and to assess its climate feedbacks.

The large-scale injection of carbon dioxide (CO2) into saline aquifers is a promising tool for reducing atmospheric CO2 emissions to mitigate climate change. An accurate assessment of the post-injection migration and trapping of the buoyant plume of CO2 is essential for estimates of storage capacity and security, but these physical processes are not fully understood. In Part 1 of this series, we presented a complete solution to a theoretical model for the migration and capillary trapping of a plume of CO2 in a confined, sloping aquifer with a natural groundwater through-flow. Here, we incorporate solubility trapping, where CO2 from the buoyant plume dissolves into the ambient brine via convective mixing. We develop semi-analytical solutions to the model in two limiting cases: when the water beneath the plume saturates with dissolved CO2 very slowly or very quickly (‘instantaneously’) relative to plume motion. We show that solubility trapping can greatly slow the speed at which the plume advances, and we derive an explicit analytical expression for the position of the nose of the plume as a function of time. We then study the competition between capillary and solubility trapping, and the impact of solubility trapping on the storage efficiency, a macroscopic measure of plume migration. We show that solubility trapping can increase the storage efficiency by several-fold, even when the fraction of CO2 trapped by solubility trapping is small.

[1] The strong coupling between multiphase flow and sediment mechanics determines the spatial distribution and migration dynamics of gas percolating through liquid-filled soft granular media. Here, we investigate, by means of controlled experiments and computed tomography (CT) imaging, the preferential mode of gas migration in three-dimensional samples of water-saturated silica-sand and silica-silt sediments. Our experimental system allowed us to independently control radial and axial confining stresses and pore pressure while performing continuous x-ray CT scanning. The CT image analysis of the three-dimensional gas migration provides the first experimental confirmation that capillary invasion preferentially occurs in coarse-grained sediments whereas grain displacement and conduit openings are dominant in fine-grained sediments. Our findings allow us to rationalize prior field observations and pore-scale modeling results, and provide critical experimental evidence to explain the means by which conduits for the transit of methane gas may be established through the gas hydrate stability zone in oceanic sediments, and cause large episodic releases of carbon into the deep ocean.

Flow through lattice networks with quenched disorder exhibits a strong correlation in the velocity field, even if the link transmissivities are uncorrelated. This feature, which is a consequence of the divergence-free constraint, induces anomalous transport of passive particles carried by the flow. We propose a Lagrangian statistical model that takes the form of a continuous time random walk with correlated velocities derived from a genuinely multidimensional Markov process in space. The model captures the anomalous (non-Fickian) longitudinal and transverse spreading, and the tail of the mean first-passage time observed in the Monte Carlo simulations of particle transport. We show that reproducing these fundamental aspects of transport in disordered systems requires honoring the correlation in the Lagrangian velocity. © 2011 American Physical Society.

Viscous fingering is a well-known hydrodynamic instability that sets in when a less viscous fluid displaces a more viscous fluid. When the two fluids are miscible, viscous fingering introduces disorder in the velocity field and exerts a fundamental control on the rate at which the fluids mix. Here we analyze the characteristic signature of the mixing process in viscously unstable flows, by means of high-resolution numerical simulations using a computational strategy that is stable for arbitrary viscosity ratios. We propose a reduced-order model of mixing, which, in the spirit of turbulence modeling and in contrast with previous approaches, recognizes the fundamental role played by the mechanical dissipation rate. The proposed model captures the nontrivial interplay between channeling and creation of interfacial area as a result of viscous fingering.

[1] It has been suggested that volume expansion caused by hydrate dissociation in sediment pores can result in large overpressure, which in turn may destabilize the sediment and trigger massive submarine landslides. Here, we investigate the pressure evolution during thermally-induced dissociation, by means of a pore-scale model that couples dissociation kinetics, multiphase flow and geomechanics. Dissociation is controlled by a self-preservation mechanism: increasing pore pressure reduces the driving force for dissociation. Hence, the overpressure is constrained by the phase equilibrium pressure, regardless of the kinetic rate of dissociation, heat supply, and sediment permeability. Furthermore, we find that the timescale for buildup of pressure by dissociation is typically much larger than that for its dissipation by drainage. Consequently, the overpressure is controlled by the capillary entry thresholds, which depend on the mode of gas invasion. In low-permeability systems, fracturing is the preferred mechanism, occurring at capillary pressures lower than the entry thresholds in the undeformed sediment. Our results suggest that while large overpressures cannot be sustained by rapid dissociation in natural systems, dissociation can induce important geomechanical effects. Gas migration by fracturing provides a possible link between dissociation, sediment deformation and methane venting.

Viscous fingering is a well-known hydrodynamic instability that sets in when a less viscous fluid displaces a more viscous fluid. When the two fluids are miscible, viscous fingering introduces disorder in the velocity field and exerts a fundamental control on the rate at which the fluids mix. We present a fluid dynamics video of the mixing process in a viscously unstable flow, generated from a high-resolution numerical simulation using a computational strategy that is stable for arbitrary viscosity ratios. We develop a two-equation dynamic model of concentration variance and mean dissipation rate to quantify the degree of mixing in such a displacement process. The model reproduces accurately the evolution of these two quantities as observed in high-resolution numerical simulations and captures the nontrivial interplay between channeling and creation of interfacial area as a result of viscous fingering., Comment: Fluid dynamics video is included

[1] Methane is a potent greenhouse gas, but its effects on Earth's climate remain poorly constrained, in part due to uncertainties in global methane fluxes to the atmosphere. An important source of atmospheric methane is the methane generated in organic-rich sediments underlying surface water bodies, including lakes, wetlands, and the ocean. The fraction of the methane that reaches the atmosphere depends critically on the mode and spatiotemporal characteristics of free-gas venting from the underlying sediments. Here we propose that methane transport in lake sediments is controlled by dynamic conduits, which dilate and release gas as the falling hydrostatic pressure reduces the effective stress below the tensile strength of the sediments. We test our model against a four-month record of hydrostatic load and methane flux in Upper Mystic Lake, Mass., USA, and show that it captures the complex episodicity of methane ebullition. Our quantitative conceptualization opens the door to integrated modeling of methane transport to constrain global methane release from lakes and other shallow-water, organic-rich sediment systems, and to assess its climate feedbacks.

Injection of carbon dioxide (CO2) into geological formations is widely regarded as a promising tool for reducing global atmospheric CO2 emissions. To evaluate injection scenarios, estimate reservoir capacity and assess leakage risks, an accurate understanding of the subsurface spreading and migration of the plume of mobile CO2 is essential. Here, we present a complete solution to a theoretical model for the subsurface migration of a plume of CO2 due to natural groundwater flow and aquifer slope, and subject to residual trapping. The results show that the interplay of these effects leads to non-trivial behaviour in terms of trapping efficiency. The analytical nature of the solution offers insight into the physics of CO2 migration, and allows for rapid, basin-specific capacity estimation. We use the solution to explore the parameter space via the storage efficiency, a macroscopic measure of plume migration. In a future study, we shall incorporate CO2 dissolution into the migration model and study the importance of dissolution relative to capillary trapping and the impact of dissolution on the storage efficiency.

We investigate the displacement of one fluid by another in a deformable medium with pore-scale disorder. We develop a model that captures the dynamic pressure redistribution at the invasion front and the feedback between fluid invasion and microstructure rearrangement. Our results suggest how to collapse the transition between invasion percolation and viscous fingering in the presence of quenched disorder. We predict the emergence of a fracturing pattern for sufficiently deformable media, in agreement with observations of drainage in granular material. We identify a dimensionless number that appears to govern the crossover from fingering to fracturing.

The emergence of scaling in transport through interconnected systems is a consequence of the topological structure of the network and the physical mechanisms underlying the transport dynamics. We study transport by advection and diffusion in scale-free and Erd\ifmmode \mbox{\H{o}}\else \H{o}\fi{}s-R\'enyi networks. Velocity distributions derived from a flow potential exhibit power-law scaling with exponent $\ensuremath{\nu}\ensuremath{\approx}\ensuremath{\gamma}+1$, where $\ensuremath{\gamma}$ is the exponent of network connectivity. Using stochastic particle simulations, we find anomalous (nonlinear) scaling of the mean-square displacement with time. We show the connection with existing descriptions of anomalous transport in disordered systems, and explain the mean transport behavior from the coupled nature of particle jump lengths and transition times.

We study a sharp-interface mathematical model of CO[subscript 2] migration in deep saline aquifers, which accounts for gravity override, capillary trapping, natural groundwater flow, and the shape of the plume during the injection period. The model leads to a nonlinear advection–diffusion equation, where the diffusive term is due to buoyancy forces, not physical diffusion. For the case of interest in geological CO[subscript 2] storage, in which the mobility ratio is very unfavorable, the mathematical model can be simplified to a hyperbolic equation. We present a complete analytical solution to the hyperbolic model. The main outcome is a closed-form expression that predicts the ultimate footprint on the CO2 plume, and the time scale required for complete trapping. The capillary trapping coefficient and the mobility ratio between CO[subscript 2] and brine emerge as the key parameters in the assessment of CO[subscript 2] storage in saline aquifers. Despite the many approximations, the model captures the essence of the flow dynamics and therefore reflects proper dependencies on the mobility ratio and the capillary trapping coefficient, which are basin-specific. The expressions derived here have applicability to capacity estimates by capillary trapping at the basin scale., Reed Research (Grant), Massachusetts Institute of Technology. Energy Initiative, ARCO Chair in Energy Studies

[1] We present a phase field model of infiltration that explains the formation of gravity fingers during water infiltration in soil. The model is an extension of the traditional Richards equation, and it introduces a new term, a fourth-order derivative in space, but not a new parameter. We propose a scaling that links the magnitude of the new term to the relative strength of gravity-to-capillary forces already present in Richards' equation. We exploit the thermodynamic framework to design a flow potential that constrains the water saturation to be between 0 and 1, its physically admissible values. The model predicts a saturation overshoot at the wetting front, which is in good agreement with experimental measurements. Two-dimensional numerical simulations predict gravity fingers with the appearance and characteristics observed in visual laboratory experiments. A linear stability analysis of the model shows that there is a direct relation between saturation overshoot and the strength of the front instability. Therefore our theory supports the conjecture that saturation overshoot, a pileup of water at the wetting front, is a prerequisite for gravity fingering.

[1] We present a discrete element model for simulating, at the grain scale, gas migration in brine-saturated deformable media. We rigorously account for the presence of two fluids in the pore space by incorporating forces on grains due to pore fluid pressures and surface tension between fluids. This model, which couples multiphase fluid flow with sediment mechanics, permits investigation of the upward migration of gas through a brine-filled sediment column. We elucidate the ways in which gas migration may take place: (1) by capillary invasion in a rigid-like medium and (2) by initiation and propagation of a fracture. We find that grain size is the main factor controlling the mode of gas transport in the sediment, and we show that coarse-grain sediments favor capillary invasion, whereas fracturing dominates in fine-grain media. The results have important implications for understanding vent sites and pockmarks in the ocean floor, deep subseabed storage of carbon dioxide, and gas hydrate accumulations in ocean sediments and permafrost regions. Our results predict that in fine sediments, hydrate will likely form in veins following a fracture network pattern, and the hydrate concentration will likely be quite low. In coarse sediments, the buoyant methane gas is likely to invade the pore space more uniformly, in a process akin to invasion percolation, and the overall pore occupancy is likely to be much higher than for a fracture-dominated regime. These implications are consistent with laboratory experiments and field observations of methane hydrates in natural systems.

We present a new mixed finite element for hexahedra, conforming in H(div). The new velocity space is similar in construction to that of BDDF of order one, but has four degrees of freedom per face instead of three. The natural degrees of freedom defining the velocity field are the components of the velocity normal to each face at the vertices of that face. We apply the new element to the numerical solution of second order elliptic equations with full tensor coefficients. We prove convergence of the method and provide optimal error estimates for the scalar and vector fields. We use an appropriate numerical quadrature to localize the mixed finite element method into a cell-centered finite difference method. For diagonal tensor coefficients, we recover the usual 7-point stencil. For full tensors, the method leads to a 27-point stencil of multipoint flux approximation (MPFA) methods—specifically the well-known O-method. Our construction allows us to prove, for the first time, convergence of the MFPA O-method on three-dimensional rectangular hexahedra. The results of the error analysis are confirmed by numerical experiments. Although our analysis is restricted to rectangular parallelepipeds, our implementation allows for distorted elements and unstructured grids.

The effect of flow instabilities on capillary trapping mechanisms is a major source of uncertainty in CO2 sequestration in deep saline aquifers. Standard macroscopic models of multiphase flow in porous media are unable to explain and quantitatively predict the onset and structure of viscous-unstable flows, such as the displacement of brine by the injected CO2. We present the first step of a research effort aimed at the experimental characterization and mathematical (continuum) modeling of such flows. Existing continuum models of multiphase flow are unable to explain why preferential flow (fingering) occurs during infiltration into homogeneous, dry soil. We present a macroscopic model that reproduces the experimentally observed features of fingered flows. The proposed model is derived using a phase-field methodology and does not introduce new independent parameters. From a linear stability analysis, we predict that finger velocity and finger width both increase with infiltration rate, and the predictions are in quantitative agreement with experiments., Eni S.p.A. (Firm) (Multiscale Reservoir Science project), Massachusetts Institute of Technology. Department of Civil and Environmental Engineering (Gilbert Winslow Career Development Chair)

Existing continuum models of multiphase flow in porous media are unable to explain why preferential flow (fingering) occurs during infiltration into homogeneous, dry soil. Following a phase-field methodology, we propose a continuum model that accounts for an apparent surface tension at the wetting front and does not introduce new independent parameters. The model reproduces the observed features of fingered flows, in particular, the higher saturation of water at the tip of the fingers, which is believed to be essential for the formation of fingers. From a linear stability analysis, we predict that finger velocity and finger width both increase with infiltration rate, and the predictions are in quantitative agreement with experiments.

The formation of preferential flow paths during infiltration of water into homogeneous, dry soil is an important phenomenon whose explanation and prediction have remained elusive under the standard theories of multiphase flow in porous media. We have recently proposed a macroscopic phase-field model of unsaturated flow in porous media, which explains why such fingering occurs [L. Cueto-Felgueroso and R. Juanes, Phys. Rev. Lett. 101, 244504 (2008)]. Here we present a linear stability analysis of the proposed model for constant-flux infiltration, which allows a quantitative description of the wetting front instability. The present analysis stresses the critical role of the initial water saturation and applied flux ratio in the asymptotic stability of the system, as well as in the transient growth of perturbations, which is consistent with the experimental evidence. The trends in the frequency and growth factor of the most unstable modes predicted by our analysis are also in quantitative agreement with experimental measurements.

We use an analytical model for the post-injection spreading of a plume of CO[subscript 2] in a saline aquifer under the action of buoyancy and capillary trapping to show that the spreading behavior is at all times strongly influenced by the shape of the plume at the end of the injection period. We solve the spreading equation numerically and confirm that, at late times, the volume of mobile CO[subscript 2] is given by existing asymptotic analytical solutions. The key parameters governing plume spreading are the mobility ratio, M, and the capillary trapping number, Gamma—the former sets the shape of the plume at the end of the injection period, and the latter sets the amount of trapping. As a quantitative measure of the dependence of the spreading behavior on the initial shape, we use a volume ratio. That is, we evolve the plume from a true end-of-injection initial shape and also from an idealized “step” initial shape, and we take the ratio of these mobile plume volumes in the asymptotic regime. We find that this volume ratio is a power-law in M, where the exponent is governed exclusively by Gamma. For conditions that are representative of geologic CO[subscript 2] sequestration, the ratio of mobile volumes between “true” and “step” initial plume shapes can be 50% or higher.

We present a sharp-interface mathematical model of CO2 migration in saline aquifers, which accounts for gravity override and capillary trapping. The major differences with respect to previous investigations is that we model the shape of the plume during the injection period, we account for regional groundwater flow during the post-injection period, and we introduce rigorously the reduction in water mobility due to trapping of the CO2. The model leads to a nonlinear advection-diffusion equation, in which the diffusive term is due to buoyancy forces, not physical diffusion. The three key dimensionless groups are the mobility ratio between the injected and initial fluids (M), the gravity number (Ng), and a newly defined trapping coefficient (Г). For the case of interest in geological CO2 storage, in which the mobility ratio is much smaller than 1, the solution is largely insensitive to the value of the gravity number. Under these conditions, the mathematical model can be simplified to a hyperbolic model. We present a complete analytical solution to the hyperbolic model that includes the injection, early post-injection, and late post-injection periods. Despite the fact that the solution involves the interaction of a sharp imbibition front with a drainage rarefaction front, it admits a closed-form expression. The main outcome of the analytical developments presented here is a formula that predicts the ultimate footprint on the CO2 plume, and the time scale required for complete trapping. Both quantities depend strongly on the shape of the plume at the end of the injection period, which must—therefore—be modeled. A second application of the analytical solution is a formulation for upscaling the capillary trapping coefficient from the laboratory scale to the basin scale. Explicit expressions are given for the megascopic and gridblock-effective trapping coefficients, as functions of the local trapping coefficient, the mobility ratio, and the grid resolution. Although the expressions derived are based on a one-dimensional sharp-interface model, we anticipate that they will have broader applicability to injection scenarios with unfavorable mobility ratio and dominated by gravity override.

We propose a novel approach to flash calculation, with particular application to negative flash. The ability to compute a negative flash for any composition state is important in practice, as the construction of analytical solutions for multicomponent systems by the method of characteristics (MOC) relies heavily on the identification of tie lines and tie-line extensions. MOC solutions are at the heart of some techniques for the calculation of the minimum miscibility pressure, and are the key building blocks for fast simulation of multidimensional reservoir flows by the front-tracking / streamline method. The basis of the proposed negative-flash method is a parameterization of the tie-line field. Rather than solving the Rachford Rice equation (or any of its variants) we solve directly for the parameters defining the tie line. For an N-component system, our approach leads to a system of N − 2 quadratic equations, which we solve efficiently using a Newton method. The iterative method is very robust: unlike other negative flash procedures, the solution displays continuous dependence on the overall composition, even in the transition to negative concentrations. We illustrate the properties and behavior of the proposed approach on three-component and four-component systems, and we then generalize the method to systems of N components. From the global triangular structure of the system with constant K-values, it follows that the system of N − 2 quadratic equations can only have two roots. For the important case of three components, the flash calculation is explicit.

[1] In uniform soils that are susceptible to unstable preferential flow, the water saturation may exhibit a nonmonotonic profile upon continuous infiltration. As this nonmonotonicity (also known as saturation overshoot) cannot be described by the conventional Richards equation, there have been proposed possible extensions to the unsaturated flow equations, including a nonmonotonic capillary pressure–saturation curve and a second-order hyperbolic term. Here, we present analytic traveling wave solutions to the extended Richards equation. These new solutions indeed display a nonmonotonic saturation profile, similar to previous simulation results. We show that these extensions need a regularization term to produce a unique solution. We develop complete analytic solutions using a relaxation regularization term, and we discuss the results in terms of recent measurements of saturation overshoot.

Existing continuum models of multiphase flow in porous media are unable to explain why preferential flow (fingering) occurs during infiltration into homogeneous, dry soil. We identify a relevant pattern-forming mechanism in the dynamics of the wetting front, and present a macroscopic model that reproduces the experimentally observed features of fingered flows. The proposed model reveals a scaling between local and nonlocal interface phenomena in imbibition, and does not introduce new independent parameters. The predictions based on this model are consistent with experiments and theories of scaling in porous media.Copyright © 2008 by ASME

Modern reservoir simulation grids are generally composed of distorted hexahedral elements populated with heterogeneous and possibly full-tensor coefficients. The numerical discretization of the reservoir flow equations on such grids is a challenging problem. Finite volume methods based on a two-point flux approximation (TPFA) do not properly account for grid distortion or permeability anisotropy that is misaligned with the grid. Multipoint flux approximation (MPFA) methods have been developed to overcome these shortcomings. Although implemented and used in virtually every commercial reservoir simulator, a proof of convergence for MPFA methods on three-dimensional hexahedral grids has remained elusive. Here, we present a link between MPFA and a new mixed finite element methods (MFEM) on hexahedral grids, which provides a powerful mathematical framework for the analysis of MPFA. First, we introduce a new mixed finite element on 3D hexahedra. The new element defines a velocity field with bilinear normal components through element faces. Thus, the new velocity field is defined by four degrees of freedom per face, which are the normal components of the velocity field at the vertices of each face. The new space is compatible with a piecewise constant pressure discretization and yields a convergent discretization. The application of a vertex-based quadrature rule reduces the new mixed finite element method to a multipoint flux control volume method. For Cartesian grids, this is in fact the classical MPFA O-method. This provides for the first time a direct link between MFEM and MPFA on hexahedral grids, which we use to establish convergence of MPFA for 3D rectangular grids.

The numerical discretization of the reservoir flow equations is a challenging problem on three-dimensional unstructured grids and in the presence of full-tensor permeabilities. Finite volume methods based on a two-point flux approximation (TPFA) do not properly account for grid distortion or permeability anisotropy that is misaligned with the grid. Over the past decade, multipoint flux approximation (MPFA) methods have been devised to address these issues. Recent research has shown, however, that MPFA methods may lead to unphysical nonmonotonic solutions for high anisotropy ratios or strongly distorted grids. In this paper, we work in the mathematical framework of mixed finite element methods (MFEM). We study a low- and high-order accurate MFEM, based on the lowest order Raviart-Thomas space (RT0) and the Brezzi-Douglas-Marini space of order one (BDM1). On tetrahedrons, the BDM1 space can be localized into a compact finite volume method through the use of an appropriate numerical quadrature. However, this is not possible for hexahedrons. The main result of this paper is the development of a new mixed finite element method for 3D hexahedral grids, which leads to a multipoint flux approximation finite volume method through the use of a vertex-based numerical quadrature. This is possible because the new velocity space is designed so that the degrees of freedom are the normal component of the velocity field at the vertices of each face. This establishes for the first time a direct link between MFEM and MPFA methods on 3D hexahedral grids. We compare the performance of classical (RT0 and BDM1) mixed finite element methods with the newly developed method on a number of test cases. We conclude that high-order methods lead to higher-quality velocity fields that are less sensitive to the choice of the grid, and that they have better monotonicity properties for highly anisotropic full-tensor permeability fields. While the added benefit of the new space is that it can be approximated by a cell-centered finite volume scheme, our numerical experiments indicate that the accuracy of the scheme is reduced.

[1] Relative permeabilities are the key descriptors in classical formulations of multiphase flow in porous media. Experimental evidence and an analysis of pore-scale physics demonstrate conclusively that relative permeabilities are not single functions of fluid saturations and that they display strong hysteresis effects. In this paper, we evaluate the relevance of relative permeability hysteresis when modeling geological CO2 sequestration processes. Here we concentrate on CO2 injection in saline aquifers. In this setting the CO2 is the nonwetting phase, and capillary trapping of the CO2 is an essential mechanism after the injection phase during the lateral and upward migration of the CO2 plume. We demonstrate the importance of accounting for CO2 trapping in the relative permeability model for predicting the distribution and mobility of CO2 in the formation. We conclude that modeling of relative permeability hysteresis is required to assess accurately the amount of CO2 that is immobilized by capillary trapping and therefore is not available to leak. We also demonstrate how the mechanism of capillary trapping can be exploited (e.g., by controlling the injection rate or alternating water and CO2 injection) to improve the overall effectiveness of the injection project.

Summary The accuracy of streamline reservoir simulations depends strongly on the quality of the velocity field and the accuracy of the streamline tracing method. For problems described on complex grids (e.g., corner-point geometry or fully unstructured grids) with full-tensor permeabilities, advanced discretization methods, such as the family of multipoint flux approximation (MPFA) schemes, are necessary to obtain an accurate representation of the fluxes across control volume faces. These fluxes are then interpolated to define the velocity field within each control volume, which is then used to trace the streamlines. Existing methods for the interpolation of the velocity field and integration of the streamlines do not preserve the accuracy of the fluxes computed by MPFA discretizations. Here we propose a method for the reconstruction of the velocity field with high-order accuracy from the fluxes provided by MPFA discretization schemes. This reconstruction relies on a correspondence between the MPFA fluxes and the degrees of freedom of a mixed finite-element method (MFEM) based on the first-order Brezzi-Douglas-Marini space. This link between the finite-volume and finite-element methods allows the use of flux reconstruction and streamline tracing techniques developed previously by the authors for mixed finite elements. After a detailed description of our streamline tracing method, we study its accuracy and efficiency using challenging test cases.

Accurate and robust discretization of the fluid flow equations is required to account for the extreme heterogeneity of oil and gas reservoirs, and the combined effect of anisotropy and grid distortion, necessary to adapt the grid to the geology. There is a growing need for handling distorted unstructured grids and full permeability tensors that appear after upscaling of the fine-scale permeability field. The classical cell-centered finite difference method results in a 7-point stencil and is insufficient to account for these effects. A number of closely related methods, coined as multipoint flux approximations (MPFA), have been proposed recently and are currently under active development. The basic principle of MPFA is that the flux across an interface between two gridblocks depends on the state variables (pressure and saturations) of more than two gridblocks. MPFA leads naturally to an enhanced finite volume method with a 27-point stencil. In this paper, we implement a variant of the MPFA method for corner-point geometry hexahedral grids. Motivated by the very high aspect ratio of gridblocks in typical reservoir models (ratios as high as 100:1 are not uncommon), we propose a hybrid method that employs a multipoint flux approximation in the areal direction and a two-point flux approximation in the vertical direction. This restricted MPFA method leads to an 11-point stencil, therefore reducing the computational effort significantly. We discuss the implementation of the method in detail. We evaluate its performance on a number of test cases and show that, for typical applications, this simplification does not greatly compromise the accuracy of the solution.

The complex physics of multiphase flow in porous media are usually modeled at the field scale using Darcy-type formulations. The key descriptors of such models are the relative permeabilities to each of the flowing phases. It is well known that, whenever the fluid saturations undergo a cyclic process, relative permeabilities display hysteresis effects. In this paper we investigate hysteresis in the relative permeability of the hydrocarbon phase in a two-phase system. We propose a new model of trapping and waterflood relative permeability, which is applicable for the entire range of rock wettability conditions. The proposed formulation overcomes some of the limitations of existing trapping and relative permeability models. The new model is validated by means of pore-network simulation of primary drainage and waterflooding. We study the dependence of trapped (residual) hydrocarbon saturation and waterflood relative permeability on several fluid/rock properties, most notably the wettability and the initial water saturation. The relevance of relative permeability hysteresis is then evaluated for modeling geological CO2 sequestration processes. Here we concentrate on CO2 injection in saline aquifers. In this setting, the CO2 is the nonwetting phase, and trapping of the CO2 is an essential mechanism after the injection phase, during the lateral and upward migration of the CO2 plume. We demonstrate the importance of accounting for CO2 trapping in the relative permeability model for predicting the distribution and mobility of CO2 in the formation. We conclude that a proper treatment of the nonwetting phase trapping leads to a higher estimate of the amount of CO2 that it is safe to inject.

This paper presents a front-tracking method for the numerical simulation of first-contact miscible gas injection processes. The method is developed for constructing very accurate (or even exact) solutions to one-dimensional initial-boundary-value problems in the form of a set of evolving discontinuities. The evolution of the discontinuities is given by analytical solutions to Riemann problems. A complete analytical Riemann solver is presented along with methods for simplifying the wave structure for Riemann problems of small amplitude. Several representative examples are used to illustrate the excellent behavior of the front-tracking method. The front-tracking method can be extended to simulate higher-dimensional processes through the use of streamlines. The paper presents an application of this computational framework for the simulation of miscible flooding in a three-dimensional, highly heterogeneous formation, and demonstrate that a miscible water-alternating-gas injection scheme is more efficient than waterflooding or gas injection alone.