Geofluids

Volume 2018 (2018), Article ID 6873298, 25 pages

https://doi.org/10.1155/2018/6873298

## Migration of Gas in Water Saturated Clays by Coupled Hydraulic-Mechanical Model

IRSN, PSE-ENV/SEDRE/LETIS, BP 17, 92260 Fontenay-aux-Roses, France

Correspondence should be addressed to Magdalena Dymitrowska

Received 9 June 2017; Accepted 11 December 2017; Published 11 March 2018

Academic Editor: Ian Clark

Copyright © 2018 Aliaksei Pazdniakou and Magdalena Dymitrowska. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Understanding the gas migration in highly water saturated sedimentary rock formations is of great importance for safety of radioactive waste repositories which may use these host rocks as barrier. Recent experiments on drainage in argillite samples have demonstrated that they cannot be represented in terms of standard two-phase flow Darcy model. It has been suggested that gas flows along highly localized dilatant pathways. Due to very small pore size and the opacity of the material, it is not possible to observe this two-phase flow directly. In order to better understand the gas transport, a numerical coupled hydraulic-mechanical model at the pore scale is proposed. The model is formulated in terms of Smoothed Particle Hydrodynamics (SPH) and is applied to simulate drainage within a sample reconstructed from the Focused Ion Beam (FIB) images of Callovo-Oxfordian claystone. A damage model is incorporated to take into account the degradation of elastic solid properties due to local conditions, which may lead to formation of new pathways and thus to modifications of fluid transport. The influence of the damage model as well as the possible importance of rigid inclusions is demonstrated and discussed.

#### 1. Introduction

The principal objective of this study is to improve the representation and understanding of gas migration in highly water saturated clays. This objective is closely related to several projects that use a deep sedimentary rock formation, such as Callovo-Oxfordian (COx) clay [1], as a host rock and a natural barrier for radioactive waste repository [2]. These sedimentary formations are characterized by a low permeability and are water saturated in their initial state. However, during the postclosure phase, the gases, mainly hydrogen, are being produced mostly by anaerobic corrosion of iron components of the repository. This important gas production may lead to appearance of a free gas phase and to a local increase of pressure, which could influence confining properties of the repository because of desaturation of clay pores and because of the mechanical disturbance caused by a possible microfracturing in the vicinity of the installation.

In classical approaches, several gas transport mechanisms can be active, separately or simultaneously [3, 4]. With increasing rates of gas production, gas should be transported first in dissolved state by advection and diffusion [5] and then by the flow of small gas bubbles inside the water phase and by viscocapillary flow of gas inside the water phase, or by the dilatation of existing and creation of new localized pathways [3, 6, 7]. The ultimate stage appearing at very high gas production rates would be macroscopic fracturing of the host rock; however, it is rather not expected in the case of repositories, where the gas source terms are not sufficient. Experimental evidences accumulating in recent years indicate that in case of indurated clays the gas migration may directly switch from dissolution-diffusion stage to the dilatant pathways mechanism. Indeed, when an increasing gas pressure is applied to a saturated argillite sample, at a certain gas pressure (“entry pressure”), a rapid but intermittent increase of the gas outflow is usually observed [3, 6, 7]. This gas flow is not associated with a significant sample desaturation, as it would be the case in viscocapillary flows. On the other hand, the gas flow causes a macroscopic deformation of the clay sample and the associated permeability tends to increase when repeating gas pulse tests. The prevalence of such a dilatant flow can be justified by the fact that the major part of the porosity of indurated clays exhibits pore size of a few nanometers [4]. Thus, due to associated high capillary pressure in such water saturated nanosize pores, the viscocapillary and bubbly flows would be restricted to largest pores which do not form a connected macroscopic paths and are not suffecient for the macroscopic gas transfer. On the other hand, the low tensile strength of the claystone facilitates the solid matrix cracking and subsequently the gas percolation [8]. Therefore, in order to get a better understanding of dilatant transport mode, it becomes necessary to take into account the coupling of the two phase fluid/gas flow inside the pore space with solid matrix mechanics, which constitute the scope of our study together with modelling of modification of argillite mechanical properties due to damage induced by gas high pressures. The phenomena of gas dissolution and diffusion are not considered in our model, since they are out of scope of this paper. Gas slippage is also excluded from considerations, since at high gas pressures these effects are of limited importance.

The coupling between the fluid flow hydrodynamics and solid matrix deformation has been tackled recently by several authors. In [9–11], a coupled continuous hydromechanical model with embedded fractures is presented where simple planar fractures are characterized by aperture which varies depending on local conditions. This allows taking into account the dependence of medium permeability on pore fluid pressure and improving to a certain extent the predictions of the fluid flow. In [12], the author studied dilatant flows by means of a single, nonpropagating elliptical crack containing a wetting and a nonwetting fluids and embedded in a linearly elastic solid. For this simple geometry it was shown that the onset of dilatant two-phase flow depends on a dimensionless parameter, and analytic solutions were obtained for some physical properties like crack aperture, capillary pressure, and permeabilities. However, the above-mentioned approaches do not take into account neither the fracture propagation process, nor the complicated underlying pore sizes distribution and connectivity. In [13–16], a more sophisticated continuous hydromechanical model is introduced which not only takes into account variation of permeability due to local stress and fluid pressure state, but also considers variation of solid matrix mechanical properties according to an ad hoc damage model which corresponds to emergence of new microfractures. However, in these works, the details of processes taking place at the pore scale are neglected; that is, each elementary volume of the model is characterized by its permeability; thus the actual geometry of the pore space is replaced by a coarsened property together with governing equations (in particular the fluid flow is described in terms of Darcy equation).

In recent years, due to X-ray computer tomography CT-scan and Focused Ion Beam Scanning Electron Microscope (FIB-SEM) imaging, it become possible to acquire models of porous opaque materials with submicronic resolution (from 40 *μ*m down to 5–10 nm). Thus it is possible to conduct flow simulations inside realistic pore geometries. It is to note that, for indurated argillites, the finest possible linear resolution (5–10 nm) is not yet sufficient to capture the connectivity of pores, and that most often the percolating components of pore space are formed by desiccation cracks, which can be assimilated to rough fractures. In [17], the two-phase flow was studied by Lattice Boltzmann Method inside a rough fracture extracted from a X-ray microtomography of Opalinus clay. It has been shown that, under realistic physical conditions and parameters (surface tension, pressure gradient), the nonwetting phase permeability is zero whenever this phase saturation is lower than about , while the nonwetting phase splits into bubbles which become rapidly stuck on fracture roughness. One of the limitations of this work was the lack of deformable character of the pore space, which hindered the possibility of transferring the conclusions to real argillites fissures.

In order to overcome this difficulty, we propose to implement a model taking into account both two-phase flow inside pore space and its interaction with deformable rock matrix. We aim as well at modelling of the appearance of new pore connections, and for this reason we have given up the Lattice Boltzmann Model, which is solved on underlying regular mesh and thus is not well suited for capturing of fracture initiation. It is important to stress here that these “new” connections are seen not only as true microcracks, but also as a way to take into account the underlying pores connectivity that is not visible on the CT-scans. This will allow us to work with samples that are apparently nonpercolating and to ultimately restore visible connectivity during gas migration.

The model presented in this paper is based on the framework of the Smoothed Particle Hydrodynamics (SPH) which is a Lagrangian mesh-free method [18, 19]. Recently it was successfully applied to simulations of miscible, immiscible [20], and unstable [21] fluid flows as well as to drainage [22] at pore scale. In [23], the SPH model was applied to direct simulations of the fluid flow through the samples obtained from the sandstone micro-CT images. A good agreement between the SPH simulations of two-phase flow and micromodel experiments was reported in [24]. It was also demonstrated [25, 26] that SPH is capable of simulating elastic media and that solid-fluid interactions can be relatively easily coupled in the framework of the same approach [27–30]. An interesting feature of the SPH method is its almost local character, which makes it an ideal choice for massively parallel computing, for example, with Compute Unified Device Architecture (CUDA) on Graphics Processing Units (GPU). For all these reasons and combined with the fact that different physical phenomena can be relatively easily expressed in terms of the SPH model, this method can be considered as an attractive choice for our application.

The article is organized as follows. In Section 2, the physical phenomena taken into account in our model are described. The numerical aspects concerning the discretization of the model equations are provided in Section 3. Some of the validation tests are presented in Section 4. The simulations of drainage in the sample obtained from the FIB images of the COx claystone are discussed in Section 5. The paper ends with conclusions in Section 6.

#### 2. Physical Model

##### 2.1. Fluid Phase

Fluid flow inside a porous medium can be described by different equations depending mainly on flow characteristics and the scale chosen for the flow description. At macroscopic scales the Darcy equation and its modifications are widely used [11, 13]; at microscopic scale, down to several nanometers of confinement, fluid flow can be described by the Navier-Stokes equations; when the confinement is limited to a few nanometers, fluid flow can only be described by molecular dynamics and stochastic methods which directly take into account intermolecular interactions. As it was mentioned before, most of the pores inside of a clayey material are of the size of a few nanometers [4]; however, these pores cannot be identified on images of real clay sample considered here with voxel size of 10 nm. Thus, we can employ the Navier-Stokes equations with a good accuracy. In the scope of the Navier-Stokes equations, the fluid is considered as a continuum matter. In the absence of sources and sinks, the conservation of mass is expressed in Lagrangian frame aswhere is the density, is the velocity, is the time, and denotes the material derivative.

Momentum conservation equation in Lagrangian frame can be written aswhere is the pressure, is the gravitational acceleration, and is the dynamic viscosity. The momentum equation for incompressible Newtonian fluid in Lagrangian frame can be stated asIn order to complete the description of fluid state at one point in terms of density , pressure , and velocity , one has to provide the equation of state which usually relates density to pressure. A common choice, when thermodynamics effects are neglected, is the ideal gas equation of statewhere is the speed of sound in the fluid and is the initial fluid density.

###### 2.1.1. Multiphase Flow

The key problem in the modelling of immiscible multiphase flow is the description of interface dynamics. When both of the fluid phases are described by means of the Navier-Stokes equations (1) and (3), one has to assure the continuity of viscous stress tensor at the fluid-fluid interface, where the difference in interactions and/or density between both fluids causes surface tension effects. The resulting force acting on molecules at the fluid interface tends to minimize the surface area. Surface tension generates a pressure step for a curved interface between two fluids, which can be written in terms of the Laplace lawwhere is the surface tension coefficient, and are the principal radii of curvature, and is the unit interface normal.

A widespread approach consists in introducing the continuum surface force [31–35]. We follow the version presented in [36].

In the framework of this model, a volumetric surface tension force is introduced at the fluid-fluid interface. The direction of this force is normal to the interface and the magnitude is proportional to the local curvature of the interface. For the implementation of this method, the main challenge consists in calculating the local curvature of the interface. This can be done by introduction of the color field defined aswhere subscripts and denote liquid and gas, respectively. The values of the color field can be used to calculate the local curvature of the interface as follows:where the unit normal is given byThe volumetric surface tension force is expressed asTo complete the description of the interface phenomena, the contact angle can be introduced by assigning the color function value to the solid phase so thatThese values are then used in calculations of the unit normal near the solid surface.

##### 2.2. Solid Phase

Clay material is usually composed of a mixture of various minerals. At the macroscopic scale this mixture can be considered as a continuum medium with homogeneous mechanical properties. Obviously, this is no longer true at the microscopic scale where the clayey matrix and the mineral inclusions of various sizes can be clearly distinguished. Due to their mechanical properties, the mineral inclusions can be considered as incompressible rigid bodies, while the clayey matrix exhibits mainly elastic response to the applied stress [2, 37] up to a certain point where it starts to fail. Therefore, in our model we are going to implement the possibility to take into account both rigid inclusions and elastic clayey matrix. Generally, clay matrix exhibits elastic properties of transversely isotropic material [1], which, for sake of simplicity, is replaced with isotropic elastic behaviour in our model.

###### 2.2.1. Elastic Material

The deformations inside an elastic material are given by the strain or deformation tensorwhere is the local displacement. The Hooke law linearly relates the strain tensor to the stress tensor aswhere is the stiffness tensor. For an isotropic medium, the components of the stiffness tensor are given bywhere and are the Lamé coefficients and is the Kronecker delta function. Then, (12) can be simplified aswhere is the unit tensor. This expression can be also written aswhere is the bulk modulus and can be associated with pressure inside the elastic medium. The dynamic behaviour of an elastic medium is described byThe volumetric elastic force density components of an isotropic medium can be expressed in terms of partial derivatives of the local displacement field as

###### 2.2.2. Damage Model

If pressure build-up or mechanical stress applied to the clay is high enough, its transport properties may be modified because of the elastic deformation of clayey matrix, as a result of, for example, dilatation of existing pores creating preferential gas paths. However, in the case of nonpercolating pore space the transport of gas phase can only take place while creating new microcracks.

Emergence of such microcracks can be modelled as a modification of mechanical properties of the clayey material under certain conditions. Since the exact models for this kind of features are not available, we propose here to reuse the ideas related to damage and fracture propagation in elastic materials. In the simplest isotropic scalar model, the process of degradation can be quantified by introduction of the damage variable [13, 14, 38, 39] which varies from 0 for intact to 1 for completely damaged material. Young’s modulus of the damaged elastic material is changed as follows [13]:where is Young’s modulus of the intact material and is the total damage variable.

In order to relate the damage variable to mechanical load inside the clayey matrix, we adopt the criterion proposed by [14] with modifications introduced by [13]. According to this criterion, two damage modes are accounted for by introduction of damage variables and responsible for tensile and compressive or shear damage, respectively.

To verify if the elastic material fails under the tensile load, the Rankine criterion is applied [14]:where is the tensile strength and is the smallest principal stress. If the criterion is satisfied, the damage variable associated with it is defined as [14]where is the residual tensile strength, is the strain at the elastic limit, and is the ultimate tensile strain at which the material can be considered as completely damaged.

To verify if the elastic phase fails under the shear or compressive load, the Mohr-Coulomb criterion is applied [14]where is the internal friction angle, is the uniaxial compressive strength, and is the largest principal stress. The damage variable associated with this criterion is given by [14]where is the residual compressive strength and is the compressive strain at the elastic limit. Following [13, 14], in (20) and (22) in case can be replaced with principal strains and , respectively. A schematic shape of the constitutive damage law given by (20) and (22) is presented in Figure 1.