Computational Models for Evaluating Long Term CO2 Storage in Saline Aquifers

(funded by NSF and KAUST through the Academic Excellence Alliance Program)

Geologic sequestration is a proven means of permanent CO2 greenhouse gas storage, but it is difficult to design and manage such efforts. Predictive computational simulation may be the only means to account for the lack of complete characterization of the subsurface environment, the multiple scales of the various interacting processes, the large areal extent of saline aquifers, and the need for long time predictions. This project investigates high fidelity multiscale and multiphysics algorithms necessary for simulation of multiphase flow and transport coupled with geochemical reactions and related mineralogy, and geomechanical deformation in porous media to predict changes in rock properties during sequestration. The work will result in a prototypical computational framework with advanced numerical algorithms and underlying technology for research in CO2 applications, which has been validated and verified against field-scale experimental tests.

BIGDATA: Collaborative Research: IA: F: Fractured Subsurface Characterization Using High Performance Computing & Guided by Big Data (Funded by NSF)

Natural fractures act as major heterogeneity in the subsurface that control flow and transport of subsurface fluids and chemical species. Their importance cannot be underestimated, because their transmissivity may result in undesired migration during geologic sequestration of CO2, they strongly control heat recovery from geothermal reservoirs, and they may lead to induced seismicity due to fluid injection into the subsurface. Advanced computational methods are critical to design subsurface processes in fractured media for successful environmental and energy applications. Continue reading

Ensemble-based Uncertainty Quantification Method for Fractured Reservoirs

Omar Al-Hinai, Jing Ping, Mary F. Wheeler

Reservoir production management and optimization requires the characterization of the uncertainty in reservoir description. For fractured reservoirs, the connectivity of fracture distributions is crucial for predicting production characteristics. In this case, since the rock property fields are highly non-Gaussian, a method that combines vector-based level-set parameterization technique and ensemble Kalman filter (EnKF) or estimating fracture distributions is developed. The mimetic finite differences approach is utilized as forward model.

Mimetic finite differences approach

Modeling fluid flow through fracture networks is challenging due to the geometric characteristics of fractures. Subsurface fractures are thin inclusions with large surface areas that can intersect in random configurations. Using traditional hexahedral or tetrahedral mesh generation is not a tenable option, as it is difficult to maintain mesh quality and reasonable number elements. In addition, the context of uncertainty quantification adds further challenges. Correctly sampling realizations requires the ability to select uniformly from the space of possible configurations. This means that a very large space of possible fracture configurations is to be modeled without user intervention.

For these reasons, a different approach is needed. Our approach has been to circumvent traditional mesh generation by using methods that allow for general polyhedral elements. First, a baseline rectangular grid is generated over the computational domain. Then, using an in-house utility, simple polygon division operations incorporate the fracture into the mesh. The utility also allows for mesh quality enhancements through simple vertex shifting operations.  This produces a very fast and robust system for meshing complex fracture networks.

The meshes are then passed into Mimpy for forward modeling of the two-phase fluid system. Mimpy is a open source three dimensional implementation for solving two-phase flow equations through porous media. Mimpy uses the Mimetic Finite Difference method, which is a conservative technique similar to the Mixed Finite Element method and the Finite Volume method.

Reference field for fractures and producers.

3D meshes

Parameterization using vector-based level-set method

For reservoirs with complex geology, a good estimation of the geological structures is very important for predicting and optimizing reservoir production. As the property fields of complex reservoirs are usually of bimodal or multimodal distributions, the Gaussian limitation is a major challenge for the application of the EnKF in the estimation of complex reservoirs [1, 2]. In order to improve the performance of the EnKF for highly non-Gaussian problems, we developed vector-based level-set parameterization method. The work is motivated by the method proposed in [3], in the sense that the level set function is applied as an indicator by which the facies fields can be transformed into smoother parameter fields. However, in the vector-based level set parameterization method, not only is the level set function used to describe the facies fields, a group of proper parameters combined with the level set function, which is called the parameter vector, is applied to depict complex shaped facies domains [4]. Compared to the original multimodal distributed parameters, the transformed parameters are in better agreement with the EnKF Gaussianity limitation, which can be updated using the standard EnKF. This approach is flexible. For different types of complex geology, different parameter vectors can be used to describe the features of the reservoirs appropriately.


We apply this method on a synthetic two-dimensional two-phase fractured reservoir. Figure 1 is the reference field with three main fractures. Figure 2 shows fracture realizations and water saturation profiles before and after updating. Figure 3 shows the matches of production data. After updating, the features of fracture distribution in reference field could be captured, and the matches of production data are more reliable.

In conclusion, the combination of mimetic finite differences approach, level-set parameterization and EnKF provides an effective solution to address the challenges in the history matching problem of highly non-Gaussian fractured reservoirs.

Figure 1. Reference field.

Figure 1. Reference field. Black lines denote fracture. Red points are producers. Green point denotes injector.

Fracture realizations and water saturation profiles

Figure 2. Fracture realizations and water saturation profiles (A) before and (B) after updating.

Oil production rate and water cut from well 1 before and after updating.

Figure 3. (A) oil production rate and (B) water cut from well 1 before and after updating. The red curves denote the responses computed with the reference field. The blue and black curves are the responses from fracture realizations before and after updating, and the green curves denote the average.


[1] Aanonsen, S. I., Nævdal, G., Oliver, D. S., Reynolds, A. C., and Vallès, B., 2009. The ensemble Kalman filter in reservoir engineering–a Review. SPE Journal, 14 (3): 393–412. doi:10.2118/117274-PA.
[2] Oliver, D. S., & Chen, Y., 2011. Recent progress on reservoir history matching: a review. Computational Geosciences, 15(1), 185-221.
[3] Chang, H., Zhang, D., and Lu, Z., 2010. History matching of facies distribution with the EnKF and level set parameterization. Journal of Computational Physics 229 (20), 8011–8030.
[4] Ping, J., & Zhang, D., 2013. History matching of fracture distributions by ensemble Kalman filter combined with vector based level set parameterization. Journal of Petroleum Science and Engineering, 108, 288-303.
This research has been presented at several conferences and is currently being prepared for publication.

Posted in Uncategorized | Tagged

Phase Field Fracture Propagation Model

The design and evaluation of hydraulic fracturing jobs is critical for efficient production from shale oil and gas fields. The efficiency of a fracturing job depends on the interaction between hydraulic (induced) and naturally occurring discrete fractures. A rigorous fracture propagation model is therefore necessary to predict fracture growth pattern in a heterogeneous, anisotropic poroelastic medium. The propped fracture length and aperture is dependent on the poroelastic behavior of rock matrix and proppants under consideration. Additionally, we need to account for complex reservoir geometries and features including non-planar fractures, faults and barriers. A long-term production assessment of these reservoirs therefore requires a coupled reservoir-fracture flow and geo-mechanics model which takes into account above considerations. In this report, we present updates for the phase field fracture propagation and the coupled flow and geomechanics model developments for fractured poroelastic reservoirs.

Two phase field fracture propagation models were developed: (1) pressurized and (2) fluid-filled models. The pressurized approach assumes fracture pressure distribution is known a priori (given) and increases linearly with time. The fluid-filled approach, on the other hand, solves a coupled flow problem for the reservoir and fracture domains to determine pressure distribution along the fracture. The latter approach is an active area of research to determine the significance of coupled flow and poro-mechanics on crack growth. The fracture pressure is then assumed to be in equilibrium with the normal component of the reservoir stresses at the fracture interface for both approaches. A brittle fracture theory, as originally presented by Griffith, is invoked along with its underlying assumptions to determine a fracture growth rate and failure criterion.

A substantial research effort is being directed to address the modeling and computational challenges involved herein. Three different numerical solution algorithms were developed to solve the coupled system differing in the treatment of an irreversible crack growth criterion: (1) simple penalization (2) augmented Lagrangian and (3) primal-dual active set method. The active set approach has faster Newton iteration convergence rates and is therefore the method of choice for 3D simulations. A detailed study describing each of these algorithms and the associated numerical results can be found in ICES Reports 13-15, 13-25 and 14-27. The fracture propagation models using phase field approach have the following advantages:

  • The models are easy to implement and use fixed-grid topology. This avoids remeshing required to resolve an exact fracture location.
  • The fracture nucleation, propagation, kinking, and path are intrinsically determined. This reduces computational cost associated with post-processing of quantities such as stress intensity factors when compared to the widely used displacement discontinuity methods.
  • The models can handle joining and branching of multiple fractures without requiring tracking of fracture interfaces.
  • Crack growth in heterogeneous media does not require special treatment. Furthermore, the fracture aperture can be calculated using the phase-field function.
Growth of three parallel penny-shaped fractures in 3D using the active set method.

Figure 1: Growth of three parallel penny-shaped fractures in 3D using the active set method. The middle fracture growth is shunned by stress shadowing effect similar to those studied in 2D fracture growth (ICES Report 14-20).

Fractures propagating at each time in 3D heterogeneous media.

Figure 2: Fractures propagating at each time in 3D heterogeneous media. Both initial fractures grow non-planar then they join and branch. We observe the adaptive mesh around the fractures.

This work is published in the following papers:

Mikelić, Andro, Mary F. Wheeler, and Thomas Wick. “A phase field approach to the fluid filled fracture surrounded by a poroelastic medium.“ ICES Preprint 13-15, Jun 2013, link

Wheeler, M. F., T. Wick, and W. Wollner. “An augmented-Lagrangian method for the phase-field approach for pressurized fractures.” Computer Methods in Applied Mechanics and Engineering 271 (2014): 69-85.

Mikelić, Andro, Mary F. Wheeler, and Thomas Wick. “A quasi-static phase-field approach to pressurized fractures.” Nonlinearity 28, no. 5 (2015): 1371-1399

Heister, Timo, Mary F. Wheeler, and Thomas Wick. “A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach.” Computer Methods in Applied Mechanics and Engineering 290 (2015): 466-495

Wick, Thomas, Sanghyun Lee, and Mary F. Wheeler. “3D Phase-Field for Pressurized Fracture Propagation in Heterogeneous Media.”, ECCOMAS and IACM Coupled Problems, San Servolo, Venice/Italy, May 18-20 2015

Posted in Uncategorized | Tagged

Simulation of the Cranfield CO2 Injection Site with a Drucker-Prager Plasticity Model

Coupled fluid flow and geomechanics simulations have strongly supported CO2 injection planning and operations. Linear elasticity has been the popular material model in CO2 simulation for addressing rock solid material behaviors. On the other hand, nonlinear constitutive models can take into account more realistic rock formation behaviors to model complex, chemically active, and fast injecting operations. For example, failure or damage may occur for rock formation near wellbores due to high fluid injection pressures or flow rates. The damaged formation near wellbores results in the changes in rock porosity or permeability, which impacts fluid flow behaviors. Such failure or damage of rock formations can be well described by the Drucker-Prager plasticity theory.

The Drucker-Prager plasticity solid mechanics module has been implemented into IPARS (Integrated Parallel Accurate Reservoir Simulators developed at the Center for Subsurface Modeling, The University of Texas at Austin). The coupled poro-plasticity system is solved using an iterative coupling scheme: the nonlinear flow and mechanics systems are solved sequentially using the fixed-stress splitting, and iterates until convergence is obtained in the fluid fraction. To the best of our knowledge, the application of this algorithm is new for poro-plasticity. To achieve fast convergence rates for solving the nonlinear solid mechanics problems, a material integrator is consistently formulated and implemented in the IPARS geomechanics module. An enhanced parallel module for general hexahedral finite elements is also developed for IPARS for solving large-scale problems in parallel. A driver for the direct solver SuperLU has been implemented in cases when the linear systems are difficult to converge. A Cranfield CO2 injection model is set up according to the reservoir geological field data and rock plasticity parameters based on Sandia national lab experimental results. This Cranfield model is solved using IPARS and the prediction on CO2 flow and formation deformation is presented.

Figure 1 shows a comparison between poro-elastic and poro-plastic results in an injection well case with homogeneous Cranfield data and rectangular geometry. Figure 2 shows poro-elastic results in an injection well case with heterogeneous Cranfield properties and geometry. Future work includes the incorporation of plastic yielding into the initial condition so that the heterogeneous Cranfield properties and geometry can be used in the poro-plastic case. Together with stress dependent permeability and fluid fraction, this will give us an accurate predictive geomechanics model for matching CO2 injection results at the Cranfield site.

Figure 1. Comparison of elasticity (left) and plasticity (right) models with homogeneous parameters and rectangular geometry. Fluid pressure, vertical displacement, and plastic strain are shown (below).

Comparison of elasticity (left) and plasticity (right) models

Figure 2. Elasticity model with heterogeneous Cranfield properties and geometry. 3D (left) and 2D (right) plots of the vertical displacement component at final simulation time (below).

Elasticity model with heterogeneous Cranfield properties and geometry