3D Simulations of Fracture Dissolution Using the Openfoam Toolkit

Although fractures take up only a small fraction (1-2%) of the available pore space in rocks they account for the majority of the subsurface transport of minerals and contaminants. Water percolating through carbonate formations will tend to dissolve the surrounding rock, due to the acidity of the dissolved CO₂. The feedback between the flow of reactant and the increase in permeability due to dissolution of the fracture faces gives rise to a number of interesting and important questions, relating to oil and gas recovery, sequestration, dam stability and the development of caves.

Numerical models of fracture dissolution typically take averages of the flow and reactant concentration across the fracture aperture, giving rise to simplified field equations that can be readily solved on a one or two-dimensional grid. However, fracture dissolution is inherently unstable and highly localized flow paths or wormholes develop [1]. Figure 1 shows a developing wormhole from an initially flat fracture, with a small region of enhanced aperture near the center of the inlet. As the wormhole grows the aperture-averaging approximation breaks down. To study the later development of wormholes, the effects of inertia and turbulence, and the competition between different flow paths a fully three-dimensional model is desirable. Here we report some preliminary developments using the OpenFOAM [2] toolkit.

The simulations solve the Stokes equation for flow in the fracture (although it is straightforward to include fluid inertia) and the convection-diffusion equation for the reactant transport. Reactions at the fracture surfaces are modeled by linear or power-law kinetics and the rates of dissolution computed. The dissolution flux is then used to modify the positions of the surfaces. The additional libraries and modifications to the OpenFOAM source needed to perform these simulations will be outlined. The major difficulty is preserving sufficient mesh quality as the fracture opens; to accomplish this we have implemented a customized source code to relax the mesh on the fracture surfaces.

Results for uniformly dissolving fractures show an excellent agreement between 2D and 3D models. At the same time first we find that the quality of the mesh and the implementation of the boundary conditions are extremely important for the resulting solution. At the meeting we will present the results of ongoing simulations of more complex fracture systems.

Fig 1. (a) A three-dimensional simulation of fracture dissolution based on the OpenFOAM libraries. A flat fracture was seeded with a small region of enhanced aperture at the center of the inlet. The geometry and concentration are shown after breakthrough, when the reactant has reached the outlet in significant quantities. The blue color indicates regions of low concentration and red the regions of high concentration. (b) 2D slices through the fracture, perpendicular (left) and parallel (right) to the flow.

[1] Szymczak and A. J. C. Ladd. The initial stages of cave formation: Beyond the one-dimensional paradigm, Earth Planet. Sci. Lett., 2011, 301, 424-432.

[2] http://www.openfoam.com/