Hamdi A. Tchelepi

Onset of convection in a gravitationally unstable diffusive boundary layer in porous media

We present a linear stability analysis of density-driven miscible flow in porous media in the context of carbon dioxide sequestration in saline aquifers. Carbon dioxide dissolution into the underlying brine leads to a local density increase that results in a gravitational instability. The physical phenomenon is analogous to the thermal convective instability in a semi-infinite domain, owing to a step change in temperature at the boundary. The critical time for the onset of convection in such problems has not been determined accurately by previous studies. We present a solution, based on the dominant mode of the self-similar diffusion operator, which can accurately predict the critical time and the associated unstable wavenumber. This approach is used to explain the instability mechanisms of the critical time and the long-wave cutoff in a semi-infinite domain. The dominant mode solution, however, is valid only for a small parameter range. We extend the analysis by employing the quasi-steady-state approximation (QSSA) which provides accurate solutions in the self-similar coordinate system. For large times, both the maximum growth rate and the most dangerous mode decay as t 1/4 . The long-wave and the short-wave cutoff modes scale as t 1/5 and t 4/5 , respectively. The instability problem is also analysed in the nonlinear regime by high-accuracy direct numerical simulations. The nonlinear simulations at short times show good agreement with the linear stability predictions. At later times, macroscopic fingers display intense nonlinear interactions that significantly influence both the front propagation speed and the overall mixing rate. A dimensional analysis for typical aquifers shows that for a permeability variation of 1 - 3000 mD, the critical time can vary from 2000 yrs to about 10 days while the critical wavelength can be between 200m and 0.3 m.

Multi-scale finite-volume method for elliptic problems in subsurface flow simulation

In this paper we present a multi-scale finite-volume (MSFV) method to solve elliptic problems with many spatial scales arising from flow in porous media. The method efficiently captures the effects of small scales on a coarse grid, is conservative, and treats tensor permeabilities correctly. The underlying idea is to construct transmissibilities that capture the local properties of the differential operator. This leads to a multi-point discretization scheme for the finite-volume solution algorithm. Transmissibilities for the MSFV have to be constructed only once as a preprocessing step and can be computed locally. Therefore this step is perfectly suited for massively parallel computers. Furthermore, a conservative fine-scale velocity field can be constructed from the coarse-scale pressure solution. Two sets of locally computed basis functions are employed. The first set of basis functions captures the small-scale heterogeneity of the underlying permeability field, and it is computed in order to construct the effective coarse-scale transmissibilities. A second set of basis functions is required to construct a conservative fine-scale velocity field. The accuracy and efficiency of our method is demonstrated by various numerical experiments.

Gravity currents with residual trapping

Motivated by geological carbon dioxide (CO 2 ) storage, we present a vertical-equilibrium sharp-interface model for the migration of immiscible gravity currents with constant residual trapping in a two-dimensional confined aquifer. The residual acts as a loss term that reduces the current volume continuously. In the limit of a horizontal aquifer, the interface shape is self-similar at early and at late times. The spreading of the current and the decay of its volume are governed by power-laws. At early times the exponent of the scaling law is independent of the residual, but at late times it decreases with increasing loss. Owing to the self-similar nature of the current the volume does not become zero, and the current continues to spread. In the hyperbolic limit, the leading edge of the current is given by a rarefaction and the trailing edge by a shock. In the presence of residual trapping, the current volume is reduced to zero in finite time. Expressions for the up-dip migration distance and the final migration time are obtained. Comparison with numerical results shows that the hyperbolic limit is a good approximation for currents with large mobility ratios even far from the hyperbolic limit. In gently sloping aquifers, the current evolution is divided into an initial near-parabolic stage, with power-law decrease of volume, and a later near-hyperbolic stage, characterized by a rapid decay of the plume volume. Our results suggest that the efficient residual trapping in dipping aquifers may allow CO 2 storage in aquifers lacking structural closure, if CO 2 is injected far enough from the outcrop of the aquifer.

Convective dissolution of carbon dioxide in saline aquifers.

[1] Geological carbon dioxide (CO2) storage is a means of reducing anthropogenic emissions. Dissolution of CO2 into the brine, resulting in stable stratification, increases storage security. The dissolution rate is determined by convection in the brine driven by the increase of brine density with CO2 saturation. We present a new analogue fluid system that reproduces the convective behaviour of CO2‐enriched brine. Laboratory experiments and high‐resolution numerical simulations show that the convective flux scales with the Rayleigh number to the 4/5 power, in contrast with a classical linear relationship. A scaling argument for the convective flux incorporating lateral diffusion from downwelling plumes explains this nonlinear relationship for the convective flux, provides a physical picture of high Rayleigh number convection in a porous medium, and predicts the CO2 dissolution rates in CO2 accumulations. These estimates of the dissolution rate show that convective dissolution can play an important role in enhancing storage security. Citation: Neufeld,J.A.,M.A .Hesse,A.Riaz,M. A.H allworth, H. A. Tchelepi, and H. E. Huppert (2010), Convective dissolution of carbon dioxide in saline aquifers, Geophys. Res. Lett., 37, L22404, doi:10.1029/2010GL044728.

Adaptive Multiscale Finite-Volume Method for Multiphase Flow and Transport in Porous Media

We present a multiscale finite-volume (MSFV) method for multiphase flow and transport in heterogeneous porous media. The approach extends our recently developed MSFV method for single-phase flow. We use a sequential scheme that deals with flow (i.e., pressure and total velocity) and transport (i.e., saturation) separately and differently. For the flow problem, we employ two different sets of basis functions for the reconstruction of a conservative fine-scale total velocity field. Our basis functions are designed to have local support, and that allows for adaptive computation of the flow field. We use a criterion based on the time change of the total mobility field to decide when and where to recompute our basis functions. We show that at a given time step, only a small fraction of the basis functions needs to be recomputed. Numerical experiments of difficult two-dimensional and three-dimensional test cases demonstrate the accuracy, computational efficiency, and overall scalability of the method.

Adaptive fully implicit multi-scale finite-volume method for multi-phase flow and transport in heterogeneous porous media

We describe a sequential fully implicit (SFI) multi-scale finite volume (MSFV) algorithm for nonlinear multi-phase flow and transport in heterogeneous porous media. The method extends the recently developed multiscale approach, which is based on an IMPES (IMplicit Pressure, Explicit Saturation) scheme [P. Jenny, S.H. Lee, H.A. Tchelepi, Adaptive multi-scale finite volume method for multi-phase flow and transport, Multiscale, Model. Simul. 3 (2005) 50-64]. That previous method was tested extensively and with a series of difficult test cases, where it was clearly demonstrated that the multiscale results are in excellent agreement with reference fine-scale solutions and that the computational efficiency of the MSFV algorithm is much higher than that of standard reservoir simulators. However, the level of detail and range of property variability included in reservoir characterization models continues to grow. For such models, the explicit treatment of the transport problem (i.e. saturation equations) in the IMPES-based multiscale method imposes severe restrictions on the time step size, and that can become the major computational bottleneck. Here we show how this problem is resolved with our sequential fully implicit (SFI) MSFV algorithm. Simulations of large (million cells) and highly heterogeneous problems show that the results obtained with the implicit multi-scale method are in excellent agreement with reference fine-scale solutions. Moreover, we demonstrate the robustness of the coupling scheme for nonlinear flow and transport, and we show that the MSFV algorithm offers great gains in computational efficiency compared to standard reservoir simulation methods.

Gravity currents in horizontal porous layers : transition from early to late self-similarity

We investigate the evolution of a finite release of fluid into an infinite, two-dimensional, horizontal, porous slab saturated with a fluid of different density and viscosity. The vertical boundaries of the slab are impermeable and the released fluid spreads as a gravity current along a horizontal boundary. At early times the released fluid fills the entire height of the layer, and the governing equation admits a self-similar solution that is a function of the viscosity ratio between the two fluids. This early similarity solution describes a tilting interface with tips propagating as x ∝ t 1/2 . At late times the released fluid has spread along the boundary and the height of the current is much smaller than the thickness of the layer. The governing equation simplifies and admits a different similarity solution that is independent of the viscosity ratio. This late similarity solution describes a point release of fluid in a semi-infinite porous half-space, where the tip of the interface propagates as x ∝ t 1/3 . The same simplification of the governing equation occurs if the viscosity of the released fluid is much higher than the viscosity of the ambient fluid. We have obtained an expression for the time when the solution transitions from the early to the late self-similar regime. The transition time increases monotonically with increasing viscosity ratio. The transition period during which the solution is not self-similar also increases monotonically with increasing viscosity ratio, for mobility ratios larger than unity. Numerical computations describing the full evolution of the governing equation show good agreement with the theoretical results. Estimates of the spreading of injected fluids over long times are important for geological storage of CO 2 , and for the migration of pollutants in aquifers. In all cases it is important to be able to anticipate when the spreading regime transitions from x ∝ t 1/2 to x ∝ t 1/3 .

Stability, Accuracy and Efficiency of Sequential Methods for Coupled Flow and Geomechanics

This paper (SPE 119084) was accepted for presentation at the SPE Reservoir Simulation Symposium, The Woodlands, Texas, USA, 2–4 February 2009, and revised for publication. Original manuscript received for review 14 November 2008. Revised manuscript received for review 24 June 2010. Paper peer approved 19 July 2010. Summary We perform detailed stability and convergence analyses of sequential-implicit solution methods for coupled fluid flow and reservoir geomechanics. We analyze four different sequential-implicit solution strategies, where each subproblem (flow and mechanics) is solved implicitly: two schemes in which the mechanical problem is solved first—namely, the drained and undrained splits—and two schemes in which the flow problem is solved first—namely, the fixed-strain and fixed-stress splits. The von Neumann method is used to obtain the linear-stability criteria of the four sequential schemes, and numerical simulations are used to test the validity and sharpness of these criteria for representative problems. The analysis indicates that the drained and fixed-strain splits, which are commonly used, are conditionally stable and that the stability limits depend only on the strength of coupling between flow and mechanics and are independent of the timestep size. Therefore, the drained and fixed-strain schemes cannot be used when the coupling between flow and mechanics is strong. Moreover, numerical solutions obtained using the drained and fixed-strain sequential schemes suffer from oscillations, even when the stability limit is honored. For problems where the deformation may be plastic (nonlinear) in nature, the drained and fixed-strain sequential schemes become unstable when the system enters the plastic regime. On the other hand, the undrained and fixed-stress sequential schemes are unconditionally stable regardless of the coupling strength, and they do not suffer from oscillations. While both the undrained and fixed-stress schemes are unconditionally stable, for the cases investigated we found that the fixed-stress split converges more rapidly than the undrained split. On the basis of these findings, we strongly recommend the fixed-stress sequential-implicit method for modeling coupled flow and geomechanics in reservoirs.

Parallel Scalable Unstructured CPR-Type Linear Solver for Reservoir Simulation

We describe a multistage parallel linear solver framework developed as part of the Intersect (IX) next-generation reservoir simulation project. The object-oriented framework allows for wide flexibility in the number of stages, methods and preconditioners. Here, we describe the specific components of a two-stage CPR 1 (Constraint Pressure Residual) scheme designed for large-scale parallel, structured and unstructured linear systems. We developed a highly efficient in-house Parallel Algebraic Multigrid (PAMG) solver as the first stage preconditioner. For the second stage, we use a parallel ILU-type scheme. This new and powerful combination of CPR and PAMG was the result of detailed analysis of the linear system of equations associated with reservoir simulation. Using several difficult reservoir simulation problems, we demonstrate the robustness and excellent parallel scalability of the IX linear solver. For the field case studies, the IX linear solver with CPR and PAMG is at least five times faster than an established and widely used industrial linear solver. The performance advantage of the IX linear solver over traditional reservoir simulation linear solvers increases with both problem size and the number of processors.

Dispersion, permeability heterogeneity, and viscous fingering: Acoustic experimental observations and particle‐tracking simulations

Stable and unstable displacement experiments were performed in millstone and limestone cores. Concentration histories at ten locations along the core samples were obtained by acoustic measurements. Particle‐tracking simulations of the displacements were also made utilizing permeability distributions measured with a permeameter. The combination of experimental observations and simulations indicate that superstable (M<1) displacements suppress the influence of heterogeneity; this suppression was reflected in smaller apparent dispersivities as the mobility ratio decreased below unity. In the millstone, which exhibited random heterogeneity, two‐dimensional particle‐tracking simulations reproduce with reasonable accuracy the growth of the fingered region in unstable displacements. In homogeneous porous media, concentration histories obtained in three‐dimensional simulations did not differ significantly from their two‐dimensional counterparts. In the more heterogeneous limestone, unstable displacements accentua...