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

Aliaksei Pazdniakou, Magdalena Dymitrowska

► To cite this version:

Aliaksei Pazdniakou, Magdalena Dymitrowska. Migration of Gas in Water Saturated Clays by Coupled Hydraulic-Mechanical Model. Geofluids, 2018, Aquitard Fluids and Gases, pp.1-25. 10.1155/2018/6873298. hal-02866936

HAL Id: hal-02866936

https://hal.science/hal-02866936v1

Submitted on 12 Jun 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Research Article

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

**Aliaksei Pazdniakou** and **Magdalena Dymitrowska**

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

Correspondence should be addressed to Magdalena Dymitrowska; magdalena.dymitrowska@irsn.fr

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.

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 sufficient 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 2D 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 3D 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 3D 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 0.4, 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 micro-model 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 as

$$ \frac{D\rho}{Dt} = -\rho\nabla \cdot \mathbf{v}, \qquad (1) $$

where $\rho$ is the density, $\mathbf{v}$ is the velocity, $t$ is the time, and $D$ denotes the material derivative.

Momentum conservation equation in Lagrangian frame can be written as

$$ \rho \frac{D\mathbf{v}}{Dt} = -\nabla p + \nabla \cdot \left[ \mu \left( \nabla \mathbf{v} + (\nabla \mathbf{v})^T - \frac{2}{3} (\nabla \cdot \mathbf{v}) \mathbf{I} \right) \right] + \rho \mathbf{g}, \qquad (2) $$

where $p$ is the pressure, $\mathbf{g}$ is the gravitational acceleration, and $\mu$ is the dynamic viscosity. The momentum equation for incompressible Newtonian fluid in Lagrangian frame can be stated as

$$ \rho \frac{D\mathbf{v}}{Dt} = -\nabla p + \mu \nabla^2 \mathbf{v} + \rho \mathbf{g}. \qquad (3) $$

In order to complete the description of fluid state at one point in terms of density $\rho$, pressure $p$, and velocity $\mathbf{v}$, 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 state

$$ p = c_s^2 (\rho - \rho_0), \qquad (4) $$

where $c_s$ is the speed of sound in the fluid and $\rho_0$ 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 $\Delta p$ for a curved interface between two fluids, which can be written in terms of the Laplace law

$$ \Delta p = \sigma \left( \frac{1}{R_1} + \frac{1}{R_2} \right) = \sigma \nabla \cdot \hat{\mathbf{n}}, \qquad (5) $$

where $\sigma$ is the surface tension coefficient, $R_1$ and $R_2$ are the principal radii of curvature, and $\hat{\mathbf{n}}$ 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 $F_s(\mathbf{r})$ 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 as

$$ c(\mathbf{r}) = \frac{\rho_l(\mathbf{r}) - \rho_g(\mathbf{r})}{\rho_l(\mathbf{r}) + \rho_g(\mathbf{r})}, \qquad (6) $$

where subscripts $l$ and $g$ denote liquid and gas, respectively. The values of the color field $c(\mathbf{r})$ can be used to calculate the local curvature $\kappa$ of the interface as follows:

$$ \kappa = -\nabla \cdot \hat{\mathbf{n}} = \frac{\mathbf{n} \cdot \nabla \|\mathbf{n}\|}{\|\mathbf{n}\|^2} - \frac{\nabla \cdot \mathbf{n}}{\|\mathbf{n}\|}, \qquad (7) $$

where the unit normal $\hat{\mathbf{n}}$ is given by

$$ \hat{\mathbf{n}} = \frac{\nabla c(\mathbf{r})}{\|\nabla c(\mathbf{r})\|}. \qquad (8) $$

The volumetric surface tension force $F_s(\mathbf{r})$ is expressed as

$$ F_s(\mathbf{r}) = \frac{1}{2} \sigma \kappa \mathbf{n}. \qquad (9) $$

To complete the description of the interface phenomena, the contact angle $\theta$ can be introduced by assigning the color function value to the solid phase so that

$$ c(\mathbf{r}) = \cos(\theta). \qquad (10) $$

These 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 tensor

$$ \boldsymbol{\varepsilon} = \frac{1}{2} (\nabla \mathbf{u} + (\nabla \mathbf{u})^T), \quad (11) $$

where **u** is the local displacement. The Hooke law linearly relates the strain tensor **ε** to the stress tensor **σ** as

$$ \boldsymbol{\sigma} = \mathbf{C} : \boldsymbol{\varepsilon}, \quad (12) $$

where **C** is the stiffness tensor. For an isotropic medium, the components of the stiffness tensor are given by

$$ C_{ijkl} = \lambda_s \delta_{ij} \delta_{kl} + \mu_s (\delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}), \quad (13) $$

where $\lambda_s$ and $\mu_s$ are the Lamé coefficients and $\delta_{ij}$ is the Kronecker delta function. Then, (12) can be simplified as

$$ \begin{aligned} \boldsymbol{\sigma} &= 2\mu_s \boldsymbol{\varepsilon} + \lambda_s (\nabla \cdot \mathbf{u}) \mathbf{I} \\ &= \mu_s (\nabla \mathbf{u} + (\nabla \mathbf{u})^T) + \lambda_s (\nabla \cdot \mathbf{u}) \mathbf{I}, \end{aligned} \quad (14) $$

where **I** is the unit tensor. This expression can be also written as

$$ \boldsymbol{\sigma} = \mu_s \left( \nabla \mathbf{u} + (\nabla \mathbf{u})^T - \frac{2}{3} (\nabla \cdot \mathbf{u}) \mathbf{I} \right) + K_s (\nabla \cdot \mathbf{u}) \mathbf{I}, \quad (15) $$

where $K_s$ is the bulk modulus and $K_s(\nabla \cdot \mathbf{u})$ can be associated with pressure inside the elastic medium. The dynamic behaviour of an elastic medium is described by

$$ \rho_s \frac{\partial^2 \mathbf{u}}{\partial t^2} = \nabla \cdot \boldsymbol{\sigma}. \quad (16) $$

The volumetric elastic force density components $\nabla \cdot \boldsymbol{\sigma} = (F_x, F_y, F_z)$ is the **F**_elastic of an isotropic medium can be expressed in terms of partial derivatives of the local displacement field as

$$ \begin{aligned} F_x &= (2\mu_s + \lambda_s) \frac{\partial^2 u_x}{\partial x \partial x} + \mu_s \left[ \frac{\partial^2 u_x}{\partial y \partial y} + \frac{\partial^2 u_x}{\partial z \partial z} \right] \\ &\quad + (\mu_s + \lambda_s) \left[ \frac{\partial^2 u_y}{\partial x \partial y} + \frac{\partial^2 u_y}{\partial x \partial z} \right] \\ F_y &= (2\mu_s + \lambda_s) \frac{\partial^2 u_y}{\partial y \partial y} + \mu_s \left[ \frac{\partial^2 u_y}{\partial x \partial x} + \frac{\partial^2 u_y}{\partial z \partial z} \right] \\ &\quad + (\mu_s + \lambda_s) \left[ \frac{\partial^2 u_x}{\partial x \partial y} + \frac{\partial^2 u_x}{\partial z \partial y} \right] \\ F_z &= (2\mu_s + \lambda_s) \frac{\partial^2 u_z}{\partial z \partial z} + \mu_s \left[ \frac{\partial^2 u_z}{\partial x \partial x} + \frac{\partial^2 u_z}{\partial y \partial y} \right] \\ &\quad + (\mu_s + \lambda_s) \left[ \frac{\partial^2 u_x}{\partial x \partial z} + \frac{\partial^2 u_y}{\partial y \partial z} \right]. \end{aligned} \quad (17) $$

### 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 $D$ [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]:

$$ E = (1 - D) E_0, \quad (18) $$

where $E_0$ is Young's modulus of the intact material and $D$ 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 $D_t$ and $D_c$ 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]:

$$ \sigma_3 \le -f_t, \quad (19) $$

where $f_t$ is the tensile strength and $\sigma_3$ is the smallest principal stress. If the criterion is satisfied, the damage variable $D_t$ associated with it is defined as [14]

$$ D_t = \begin{cases} 0 & \text{if } \varepsilon_{t0} \le \varepsilon \\ 1 - \frac{f_{tr}}{E_0 \varepsilon} & \text{if } \varepsilon_{tu} \le \varepsilon \le \varepsilon_{t0} \\ 1 & \text{if } \varepsilon \le \varepsilon_{tu}, \end{cases} \quad (20) $$

where $f_{tr}$ is the residual tensile strength, $\varepsilon_{t0}$ is the strain at the elastic limit, and $\varepsilon_{tu}$ 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]

$$ \sigma_1 - \sigma_3 \frac{1 + \sin \phi}{1 - \sin \phi} \ge f_c, \quad (21) $$

where $\phi$ is the internal friction angle, $f_c$ is the uniaxial compressive strength, and $\sigma_1$ is the largest principal stress. The damage variable $D_c$ associated with this criterion is given by [14]

$$ D_c = \begin{cases} 0 & \text{if } \varepsilon \le \varepsilon_{c0} \\ 1 - \frac{f_{cr}}{E_0 \varepsilon} & \text{if } \varepsilon_{c0} \le \varepsilon, \end{cases} \quad (22) $$
<figure><img src="image_2.png" /><figcaption>FIGURE 1: An example of constitutive damage law.</figcaption></figure>

<table><caption>TABLE 1: COx clay properties.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Young's modulus <i>E</i></td><td>3 GPa</td></tr><tr><td>Poisson's ratio <i>ν</i></td><td>0.1</td></tr><tr><td>Tensile strength <i>f</i><sub>t</sub></td><td>4 MPa</td></tr><tr><td>Uniaxial compressive strength <i>f</i><sub>c</sub></td><td>20 MPa</td></tr><tr><td>Internal friction angle</td><td>25°</td></tr><tr><td>Strain at elastic limit <i>ε</i><sub>ε0</sub></td><td>0.0013</td></tr><tr><td>Ultimate tensile strain <i>ε</i><sub>tu</sub></td><td>0.0026</td></tr><tr><td>Compressive strain at elastic limit <i>ε</i><sub>ε0</sub></td><td>0.015</td></tr></tbody></table>

where $f_{cr}$ is the residual compressive strength and $\epsilon_{c0}$ is the compressive strain at the elastic limit. Following [13, 14], $\epsilon$ in (20) and (22) in 3D case can be replaced with principal strains $\epsilon_3$ and $\epsilon_1$, respectively. A schematic shape of the constitutive damage law given by (20) and (22) is presented in Figure 1.

It can be assumed that, for a certain load range, the material follows linear elasticity law, and outside this range a nonlinear damage behaviour is observed, which coincides with observations made for clayey materials [2, 37]. It is well known that the Mohr-Coulomb criterion works fine for materials whose tensile strength $f_t$ is significantly smaller than their uniaxial compressive strength $f_c$ which is the case for the COx clay. Although, based on the experimental data [2], other specific damage criteria are available for the COx clay, they can be fitted with the Mohr-Coulomb criterion for the sake of generality, which is done in [2]. Based on these data, the parameters adopted for COx clay are summarized in Table 1.

The total damage $D$ is a combination of $D_t$ and $D_c$ and can be calculated as [13]

$$D = \zeta_t D_t + \zeta_c D_c, \qquad (23)$$

where the coefficients $\zeta_t$ and $\zeta_c$ are calculated from the local stress-strain state as described in [13, 38]. It must be noted that the scalar damage variable $D$ can be replaced with a tensor damage variable [39] which may influence the orientation of microcracks.

In order to avoid sharp restoration of the damage variable which may cause oscillations of elastic properties, the relaxation is applied to obtain a restored damage value $D(t + \Delta t)$ based on $D$ value calculated by (23)

$$D(t + \Delta t) = D(t) + \frac{D - D(t)}{\tau}, \qquad (24)$$

where $\tau$ is the dimensionless relaxation constant which is set to 10 in our simulations.

### 2.2.3. Rigid Bodies.
In their natural state, clayey rocks are highly heterogeneous, especially at the microscale where various rigid inclusions, such as quartz ($E \approx 72$ GPa) or calcite, are observed in a softer clayey matrix. Contact lines between these rigid inclusions and the clayey matrix are considered as probable sites for emergence of the new gas passages. Therefore, this component should be included into our model for further investigation.

Rigid body dynamics is fully described by the Newton-Euler equations,

$$\frac{d\mathbf{v}_{\text{CM}}}{dt} = \frac{\mathbf{F}_B}{m_B} \qquad (25a)$$

$$\frac{d\mathbf{r}_{\text{CM}}}{dt} = \mathbf{v}_{\text{CM}} \qquad (25b)$$

$$\frac{d\boldsymbol{\omega}_B}{dt} = \mathbf{I}_B^{-1} (\mathbf{T}_B - \boldsymbol{\omega}_B \times \mathbf{I}_B \boldsymbol{\omega}_B) \qquad (25c)$$

$$\frac{d\boldsymbol{\phi}_B}{dt} = \boldsymbol{\omega}_B, \qquad (25d)$$

where $\mathbf{r}_{\text{CM}}$ is the position of the center of mass and $\mathbf{v}_{\text{CM}}$ its velocity, $\mathbf{F}_B$ is the total external force acting on the body, $m_B$ is the mass of the body, $\boldsymbol{\omega}_B$ is the angular velocity, $\mathbf{I}_B$ is the tensor of the moment of inertia, $\mathbf{T}_B$ is the total torque acting on the body, and $\boldsymbol{\phi}_B$ are the angles between the body axis and principal coordinate axes. The first two equations describe movement of the center of mass, while the last two describe the rigid body rotation. The position of the center of mass can be found as

$$\mathbf{r}_{\text{CM}} = \frac{1}{m_B} \int \rho(\mathbf{r}) \mathbf{r} \mathrm{d}^3\mathbf{r}. \qquad (26)$$

The tensor of the moment of inertia is given by

$$\mathbf{I}_B = \int \rho(\mathbf{r}) [(\mathbf{r} - \mathbf{r}_{\text{CM}})^2 \mathbf{I} - (\mathbf{r} - \mathbf{r}_{\text{CM}})(\mathbf{r} - \mathbf{r}_{\text{CM}})^T] \mathrm{d}^3\mathbf{r}. \qquad (27)$$

The total torque acting on the body can be calculated as

$$\mathbf{T}_B = \int (\mathbf{r} - \mathbf{r}_{\text{CM}}) \times \mathbf{F} \mathrm{d}^3\mathbf{r}, \qquad (28)$$

where $\mathbf{F}$ is the external force density.

The Newton-Euler equations as described by (25a), (25b), (25c), and (25d) can be integrated by any appropriate numerical method.
**2.3. Boundary Conditions.** In order to provide a correct inter-action between phases at the solid-fluid interface, the boundary conditions must be defined. In most cases, the velocity of both phases at the contact point $\mathbf{r}_i$ is supposed to be the same:

$$ \mathbf{v}_f(\mathbf{r}_i) = \mathbf{v}_s(\mathbf{r}_i). \qquad (29) $$

This condition defines mutual nonpenetration of the phases and the no-slip boundary condition and must be completed with the condition of continuity of normal stresses at the solid-fluid interface

$$ \boldsymbol{\sigma}_s \cdot \mathbf{n}_s = -\boldsymbol{\sigma}_f \cdot \mathbf{n}_f, \qquad (30) $$

where $\boldsymbol{\sigma}_f$ and $\boldsymbol{\sigma}_s$ are fluid and solid phase stress tensors, respectively, $\mathbf{n}_f$ is the unit normal pointing out of the fluid phase, and $\mathbf{n}_s = -\mathbf{n}_f$.

At the interface between the solid phases, the continuity of normal stresses is imposed:

$$ \boldsymbol{\sigma}_s \cdot \mathbf{n}_s = -\boldsymbol{\sigma}_e \cdot \mathbf{n}_e. \qquad (31) $$

For a more realistic model, one should also add the friction force between solid components. In this work the friction force is omitted, but it may be added in the future when major phenomena will be sufficiently well understood.

## 3. Numerical Modeling

**3.1. Smoothed Particle Hydrodynamics.** The model presented in previous section must be solved numerically. In this paper, this is done by means of the Smoothed Particle Hydrodynamics (SPH) which is a Lagrangian mesh-free method originally applied for solution of astrophysical problems [18, 19]. Nowadays the method is widely used in various engineering and scientific areas [40, 41]. For our purposes it is important that the method can be successfully implemented to simulate fluid hydrodynamics, elastic, plastic, and rigid bodies as well as coupled solid-fluid systems. All these simulations can be performed within the same conceptual framework which decreases significantly the amount and the complexity of work to be done. Another considerable advantage is that the method is meshless which simplifies significantly simulations of multiphase flow with moving interfaces and also allows the description of some physical processes (such as fracture propagation) without mesh size dependency. Also the method can be easily parallelized and may profit of CUDA technology for massively parallel computations on GPU.

The SPH model follows the idea that, for a given function $f$, its value can be found as [40]

$$ f(\mathbf{x}) = \int f(\mathbf{x}') \delta(\mathbf{x} - \mathbf{x}') d\mathbf{x}', \qquad (32) $$

where $\delta(\mathbf{x})$ is the Dirac delta function and the integral is taken over the entire space. The Dirac delta function has nonzero value only in one point. However, from (32) one can expect that the value of a function at a point can be approximated by this function values in close proximity of the point by means of a similar expression [40]

$$ f(\mathbf{x}) = \int f(\mathbf{x}') W(\mathbf{x} - \mathbf{x}', h) d\mathbf{x}' + \text{err}_{\text{smoothing}}, \qquad (33) $$

where $W(\mathbf{x} - \mathbf{x}', h)$ is a smoothing function, $h$ is a smoothing length, and $\text{err}_{\text{smoothing}}$ is the approximation error. The integral in (33) represents the smoothed value of the function $f$ [40]

$$ f_S(\mathbf{x}) = \int f(\mathbf{x}') W(\mathbf{x} - \mathbf{x}', h) d\mathbf{x}'. \qquad (34) $$

The approximation error $\text{err}_{\text{smoothing}}$ can be estimated as $\text{err}_{\text{smoothing}} = O(h^2)$ by taking Taylor expansion of the right part of (34). By using the divergence theorem and integration by parts, the smoothed approximation of the spatial derivatives can be obtained as [40]

$$ (\nabla \cdot f(\mathbf{x}))_S = \int f(\mathbf{x}') \cdot \nabla W_x(\mathbf{x} - \mathbf{x}', h) d\mathbf{x}'. \qquad (35) $$

The higher order derivatives can be obtained in the similar way.

**3.1.1. Smoothing Function.** The smoothing function $W(\mathbf{x} - \mathbf{x}', h)$ must satisfy certain conditions. It must have the compact support, be normalized, and converge to the Dirac delta function when $h$ tends to zero. Most of the smoothing functions are positively defined inside the support, which has certain physical meaning. In order to get the same contribution from the points located at the same distance from the center, $W(\mathbf{x} - \mathbf{x}', h)$ must be an even function. It is also observed that smoothing functions with smooth first and higher order derivatives usually provide more stable numerical results, because of smaller sensitivity to particle disorder [40].

A variety of smoothing functions satisfying these and some additional conditions [40] can be constructed, so one can try to find a smoothing function which best fits the precise application needs in terms of accuracy, numerical stability, and performance. One of the most popular smoothing functions, especially in hydrodynamics simulations, is the cubic spline which is adopted in this paper, and in 3D space it can be expressed as [40]

$$ W(\mathbf{x} - \mathbf{x}', h) = \frac{16}{\pi h^3} \begin{cases} \frac{1}{2} - 3r^2 + 3r^3 & \text{if } 0 \le r < 0.5 \\ (1-r)^3 & \text{if } 0.5 \le r < 1 \\ 0 & \text{if } r \ge 1, \end{cases} \qquad (36) $$

where $r = \|\mathbf{x} - \mathbf{x}'\|/h$. It is clear that $W(\mathbf{x} - \mathbf{x}', h)$ has units of inverse of volume.

Smoothing functions proposed by [42] are designed to have some useful properties. Desbrun's spiky smoothing function [43] is appropriate for pressure forces calculation since the spiky shape of its gradient prevents the formation of particle clusters when they get close to each other. In [40] the quadratic smoothing function was proposed which is similar to the cubic spline function but is not a piecewise-defined function which may be advantageous for numerical calculation. Smoothing functions are mostly represented by polynomials for computational convenience; however, one may consider a Gaussian smoothing function [18] when looking for stability and accuracy.
**3.1.2. SPH Model.** In SPH model, the medium is approximated by a set of particles which are initially uniformly or randomly distributed in space. Each SPH particle is a Lagrangian particle and is characterized by its position and velocity. It also has certain physical variables such as mass $m_i$ and density $\rho_i$ associated with it. The particle volume can be defined as

$$V_i = \frac{m_i}{\rho_i}. \qquad (37)$$

The continuous integrals (34) and (35) can be approximated by their discrete analogues,

$$\begin{aligned} f_i &= \sum_j f_j W(\mathbf{x}_i - \mathbf{x}_j, h) V_j \\ &= \sum_j f_j \frac{m_j}{\rho_j} W(\mathbf{x}_i - \mathbf{x}_j, h) \end{aligned} \qquad (38a)$$

$$\begin{aligned} \nabla f_i &= \sum_j f_j \nabla W(\mathbf{x}_i - \mathbf{x}_j, h) V_j \\ &= \sum_j f_j \frac{m_j}{\rho_j} \nabla W_{\mathbf{x}_i} (\mathbf{x}_i - \mathbf{x}_j, h), \end{aligned} \qquad (38b)$$

where the summation in function term (38a) is carried out for all neighbouring particles within the distance $h$ from the $i$th particle including $j = i$; note that in the derivative term $j \neq i$ (38b). According to [44], the error introduced by the finite sum approximation of continuous integrals can be estimated as $\text{err}_{\text{integral}} = O(\Delta/h)$, where $\Delta$ is the length scale characterizing the average distance between the particles.

When $f$ is a density, (38a) is transformed to

$$\rho_i = \sum_j m_j W(\mathbf{x}_i - \mathbf{x}_j, h). \qquad (39)$$

Usually, the masses of SPH particles representing the same material are equal. Another useful quantity is the number density, which can be calculated as

$$n_i = \frac{\rho_i}{m_i} = \sum_j W(\mathbf{x}_i - \mathbf{x}_j, h) \qquad (40)$$

and also represents the inverse of the particle volume. When applied directly, (39) has a drawback to underestimate the density near the material interface, which can cause consistency issues. This can be corrected by using the normalization [40]

$$\rho_i = \frac{\sum_j m_j W(\mathbf{x}_i - \mathbf{x}_j, h)}{\sum_j \left( W(\mathbf{x}_i - \mathbf{x}_j, h) \right) / n_j}. \qquad (41)$$

If the material is in contact with another material represented by SPH particles, the density deficiency near the interface can be corrected by considering the virtual contribution to density from another material particles, which means that summation in (39) is simply extended to neighbouring material particles located within the smoothing function support.

The direct use of (38b) for calculation of derivatives is generally avoided because of low accuracy. For example, the derivative of a constant function calculated by (38b) does not vanish. This drawback can be corrected by using the following identity [41, 44]:

$$\nabla (f) = \frac{1}{F} (\nabla (Ff) - f \nabla F), \qquad (42)$$

where $F$ is a differentiable function. Choosing various $F$ functions, one obtains different formulas for derivative term. The most frequently used functions are $F = 1$, $F = \rho$, and $F = 1/\rho$, substitution of which in (42) yields the following expressions of SPH derivative term:

$$\nabla f_i = \sum_j (f_j - f_i) \frac{m_j}{\rho_j} \nabla W_{\mathbf{x}_i} (\mathbf{x}_i - \mathbf{x}_j, h) \qquad (43a)$$

$$\nabla f_i = \frac{1}{\rho_i} \sum_j (f_j - f_i) m_j \nabla W_{\mathbf{x}_i} (\mathbf{x}_i - \mathbf{x}_j, h) \qquad (43b)$$

$$\nabla f_i = \rho_i \sum_j \left( \frac{f_j}{\rho_j^2} + \frac{f_i}{\rho_i^2} \right) m_j \nabla W_{\mathbf{x}_i} (\mathbf{x}_i - \mathbf{x}_j, h). \qquad (43c)$$

The calculation of a Laplacian of some functions is often required for the transport equation. It can be obtained as

$$\nabla^2 f_i = \sum_j f_j \frac{m_j}{\rho_j} \nabla^2 W_{\mathbf{x}_i} (\mathbf{x}_i - \mathbf{x}_j, h). \qquad (44)$$

This formula also can be symmetrized just as SPH derivatives terms; however, even symmetrized formulations may suffer from $\nabla^2 W_{\mathbf{x}_i}(\mathbf{x}_i - \mathbf{x}_j, h)$ changing a sign for some smoothing functions. This may cause nonphysical behaviour when (44) is applied to calculate forces. This formulation is also sensible to particle disorder [41].

Another way to calculate the divergence of a scalar function can be resumed as follows [41, 44]:

$$\begin{aligned} \nabla \cdot (F \nabla f)_i &= \sum_j (F_j + F_i) (f_i - f_j) \\ & \cdot \frac{m_j}{\rho_j} \frac{\mathbf{x}_i - \mathbf{x}_j}{\|\mathbf{x}_i - \mathbf{x}_j\|^2} \nabla W_{\mathbf{x}_i} (\mathbf{x}_i - \mathbf{x}_j, h), \end{aligned} \qquad (45)$$

where $f$ and $F$ are some scalar functions. It is useful to note that [44]

$$\nabla W_{\mathbf{x}_i} (\mathbf{x}_i - \mathbf{x}_j, h) = \frac{\mathbf{x}_i - \mathbf{x}_j}{\|\mathbf{x}_i - \mathbf{x}_j\|} \frac{\partial W(r_{ij}, h)}{\partial r_{ij}}, \qquad (46)$$

where $r_{ij} = \|\mathbf{x}_i - \mathbf{x}_j\|$.

Using these formulations for scalars and derivative terms, it is possible to discretize various physical equations in order to obtain their SPH analogues which can be solved numerically. For example, the movement of the Lagrangian particles can be described by the following system of equations:

$$\frac{d\mathbf{x}_i}{dt} = \mathbf{v}_i \qquad (47a)$$

$$\frac{d\mathbf{v}_i}{dt} = \mathbf{a}_i, \qquad (47b)$$
where $\mathbf{v}_i$ is the particle velocity and $\mathbf{a}_i$ is the acceleration which results from the sum of the forces acting on the particle. This system can be integrated in time by using one of the explicit time integration schemes. The Leapfrog integration scheme can be stated as [44]

$$ \mathbf{v}_i \left( t + \frac{\Delta t}{2} \right) = \mathbf{v}_i \left( t - \frac{\Delta t}{2} \right) + \Delta t \mathbf{a}_i (t) \quad (48a) $$

$$ \mathbf{x}_i (t + \Delta t) = \mathbf{x}_i (t) + \Delta t \mathbf{v}_i \left( t + \frac{\Delta t}{2} \right). \quad (48b) $$

It is similar to the velocity Verlet algorithm, which can be expressed as [44]

$$ \mathbf{x}_i (t + \Delta t) = \mathbf{x}_i (t) + \Delta t \mathbf{v}_i (t) + \frac{\Delta t^2 \mathbf{a}_i (t)}{2} \quad (49a) $$

$$ \mathbf{v}_i (t + \Delta t) = \mathbf{v}_i (t) + \frac{\Delta t}{2} (\mathbf{a}_i (t) + \mathbf{a}_i (t + \Delta t)). \quad (49b) $$

Both methods are second-order integration schemes. For numerical stability the time step $\Delta t$ must obey the Courant-Friedrich-Levy condition [35]

$$ \Delta t \le 0.25 \frac{h}{\|\mathbf{v}_{\max}\|} \quad (50) $$

and particle acceleration related condition

$$ \Delta t \le 0.25 \sqrt{\frac{h}{\|\mathbf{a}_{\max}\|}}. \quad (51) $$

When viscous dissipation term is present, the following condition must be satisfied:

$$ \Delta t \le 0.125 \frac{\rho h^2}{\mu}. \quad (52) $$

The SPH is usually used to solve systems of partial differential equations which relate the time derivative of variables of interest to spatial derivatives. Since the equations of interest here describing fluid flow, elastic body deformation, and solid body movement satisfy these conditions, let us consider their discretization in detail.

### 3.1.3. SPH Discretization of Navier-Stokes Equations.
In the framework of SPH, the continuity equation (1) can be discretized in various ways depending on the formulation adopted for the spatial derivatives (43a), (43b), and (43c). The most widespread variants are [41]

$$ \frac{D\rho_i}{Dt} = \rho_i \sum_j \frac{m_j}{\rho_j} (\mathbf{v}_i - \mathbf{v}_j) \cdot \nabla W_{\mathbf{x}_i} (\mathbf{x}_i - \mathbf{x}_j, h) \quad (53a) $$

$$ \frac{D\rho_i}{Dt} = \sum_j m_j (\mathbf{v}_i - \mathbf{v}_j) \cdot \nabla W_{\mathbf{x}_i} (\mathbf{x}_i - \mathbf{x}_j, h). \quad (53b) $$

The first one is preferable when simulating a multiphase flow with large density contrasts since it provides better accuracy [41]. It is worth mentioning that these expressions conserve

the mass only approximately (error of the order of $h^2$) [44]. Another possibility to follow the evolution of fluid density consists in using (39) which conserves the mass exactly.

The right part of the momentum conservation equation (3) can be discretized by using (43a), (43b), and (43c). The pressure gradient term can be represented in SPH form by one of the following expressions:

$$ \nabla \rho_i = \sum_j \frac{m_j}{\rho_j} (p_i + p_j) \nabla W_{\mathbf{x}_i} (\mathbf{x}_i - \mathbf{x}_j, h) \quad (54a) $$

$$ \nabla \rho_i = \rho_i \sum_j \left( \frac{p_j}{\rho_j^2} + \frac{p_i}{\rho_i^2} \right) m_j \nabla W_{\mathbf{x}_i} (\mathbf{x}_i - \mathbf{x}_j, h). \quad (54b) $$

Both expressions provide similar results in terms of accuracy; however, the second expression is more stable with respect to particle disorder [41]. The use of (43a) and (43b) must be avoided for discretization of the pressure gradient term, since it becomes unstable for negative pressure values. The viscous dissipation term can be discretized by using (44) or (45). The use of the (45) should be preferred for better accuracy and numerical stability, which leads us to the following expression:

$$ \mu \nabla^2 \mathbf{v}_i = 2\mu \sum_j \frac{m_j}{\rho_j} \frac{\mathbf{v}_i - \mathbf{v}_j}{\|\mathbf{x}_i - \mathbf{x}_j\|} \frac{\partial W(r_{ij}, h)}{\partial r_{ij}}. \quad (55) $$

At the interface between two fluids (of when viscosity term is a spatial variable) one can use the original expression for viscous dissipation term based on (45) which corresponds to the arithmetic mean approximation; however, for better precision it is preferable to use the expression based on the harmonic mean

$$ \mu \nabla^2 \mathbf{v}_i = \sum_j \frac{4\mu_i\mu_j}{\mu_i + \mu_j} \frac{m_j}{\rho_j} \frac{\mathbf{v}_i - \mathbf{v}_j}{\|\mathbf{x}_i - \mathbf{x}_j\|} \frac{\partial W(r_{ij}, h)}{\partial r_{ij}}. \quad (56) $$

Thus, the discretized momentum equation can be written as follows:

$$ \frac{D\mathbf{v}_i}{Dt} = -\sum_j \left( \frac{\rho_j}{\rho_j^2} + \frac{p_i}{\rho_i^2} \right) m_j \nabla W_{\mathbf{x}_i} (\mathbf{x}_i - \mathbf{x}_j, h) \\ + \sum_j \frac{4\mu_i\mu_j}{\mu_i + \mu_j} \frac{m_j}{\rho_j\rho_i} \frac{\mathbf{v}_i - \mathbf{v}_j}{\|\mathbf{x}_i - \mathbf{x}_j\|} \frac{\partial W(r_{ij}, h)}{\partial r_{ij}} + \mathbf{g} \quad (57) \\ + \frac{\mathbf{F}_s}{\rho_i}, $$

where $\mathbf{F}_s$ is the surface tension force.

In weakly compressible formulation adopted in our model, the incompressibility of the fluid is imposed by the proper choice of the speed of sound $c_s$ in the equation of state (4). It must be prescribed based on the desired density variation, keeping in mind that high values of the speed of sound lead to small time steps and noisy pressure distributions [44]. An alternative approach consists in using the projection method [45].
Two principal approaches exist to model the surface tension by using the SPH model. The first one consists in using the pairwise interparticle force which mimics the intermolecular attraction forces [44]. This allows obtaining surface tension effects naturally by adjusting the interaction force constants for different types of particles. However, the introduction of such a force can modify the viscosity of fluids and the equations of state [44].

In our model, the surface tension effects are introduced by discretizing the continuum force model (Section 2.1.1) in terms of SPH. First, a certain value of color function must be attributed to the fluid and gas SPH particles. The simplest way to do this is to use constant values $c_f = 1$ and $c_g = -1$ for fluid and gas SPH particles, respectively. Based on these values one can define a smoothed value of the color function associated with SPH particle

$$c(\mathbf{r}_i) = \sum_j \frac{m_j}{\rho_j} c_j W(\mathbf{x}_i - \mathbf{x}_j, h). \quad (58)$$

Also, the color function can be calculated based on (6). The color function normal is calculated as

$$\mathbf{n}(\mathbf{r}_i) = \sum_j \frac{m_j}{\rho_j} (c(\mathbf{r}_j) - c(\mathbf{r}_i)) \nabla W_{\mathbf{x}_i} (\mathbf{x}_i - \mathbf{x}_j, h). \quad (59a)$$

The interface local curvature $\kappa$ given by (7) is calculated as

$$\kappa(\mathbf{r}_i) = \frac{1}{\|\mathbf{n}(\mathbf{r}_i)\|^2} \cdot \sum_j \frac{m_j}{\rho_j} (\mathbf{n}(\mathbf{r}_j) \|\mathbf{n}(\mathbf{r}_i)\| - \mathbf{n}(\mathbf{r}_i) \|\mathbf{n}(\mathbf{r}_j)\|) \cdot \nabla W_{\mathbf{x}_i} (\mathbf{x}_i - \mathbf{x}_j, h). \quad (60)$$

In order to improve the accuracy of the curvature calculation, we employ the concept of reliable normal introduced in [35]. Then, the continuous surface tension force $F_s(\mathbf{r}_i)$ is calculated by (9). Finally, the contact angle is prescribed by assigning a color function value to the solid and elastic SPH particles based on (10). Integration is performed using the Leapfrog integration scheme (48a) and (48b) while the fluid densities are updated using (39).

### 3.1.4. SPH Discretization of Elastic Bodies Dynamics.
The most frequently used method to model elastic bodies by means of SPH consists in integration of rate of change of the deviatoric stress tensor [25, 26, 28–30] in the framework of the hypoelastic constitutive model. Another approach consists in calculation of elastic forces from expansion of the elastic strain energy [46, 47]. In our model, for the sake of simplicity we discretize the equation governing the dynamic behaviour of the elastic medium (16) by means of SPH.

In order to approximate the force components by SPH, second derivatives of the local displacement can be calculated using the same approach as for the first derivatives. However, this formulation has several disadvantages [41], such as high

sensitivity to particle disorder. A better alternative is to use the following approximation of second derivatives [41]:

$$\left( \frac{\partial^2 F}{\partial x_\alpha \partial x_\beta} \right)_i = \sum_j \frac{m_j}{\rho_j} \left( 5 \frac{(x_{\alpha,i} - x_{\alpha,j})(x_{\beta,i} - x_{\beta,j})}{\|\mathbf{x}_i - \mathbf{x}_j\|} - \delta_{\alpha\beta} \right) \cdot (F_i - F_j) \frac{1}{\|\mathbf{x}_i - \mathbf{x}_j\|} \frac{\partial W(r_{ij}, h)}{\partial r_{ij}}, \quad (61)$$

where the Greek indexes $\alpha$ and $\beta$ indicate spatial dimensions and $\delta_{\alpha\beta}$ is the Kronecker delta function.

When elastic forces acting on elastic particles are calculated, the particle trajectory can be easily obtained by numerical integration of corresponding equations (47a) and (47b) of motion. In order to limit vibrations and dissipate elastic energy, a viscous damping term similar to the one used in the Navier-Stokes equation can be added to (17); another possibility consists in introducing the viscosity-like term $\theta$ in the Verlet algorithm [48]. Both options are implemented in our model.

### 3.2. Discretization of Rigid Body Dynamics.
A rigid body can be discretized with a set of Lagrangian points distributed inside the body volume and at the body surface. This spatial description allows a natural coupling of discretized equations of rigid body dynamics with the other equations in the frame of SPH [27].

In the context of our model, the movement of rigid inclusions results from forces applied on their surfaces by other phases as well as from contact forces between inclusions themselves. In order to decrease the memory use and to save calculation time, one can describe the rigid body by taking into account only the points lying at the body surface and in its close proximity (i.e., the points which may experience the external forces from other phases) without any loss of accuracy since the expected rate of rotation is low, and the final position of the rigid grain is most likely defined by its shape and distribution of external forces at the grain surface rather than by distribution of mass inside this body.

The center of mass of the rigid body represented by a set of N points of mass $m_p$ can be found from

$$\mathbf{r}_{\text{CM}} = \frac{1}{M_B} \sum_{i=1}^{N} m_p \mathbf{r}_i. \quad (62)$$

The total external force acting on the rigid body is given by

$$F_B = \sum_{i=1}^{N} m_p \mathbf{a}_i. \quad (63)$$

The tensor of the moment of inertia is obtained from

$$I_B = \sum_{i=1}^{N} m_p [(\mathbf{r}_i - \mathbf{r}_{\text{CM}})^2 \mathbf{I} - (\mathbf{r}_i - \mathbf{r}_{\text{CM}})(\mathbf{r}_i - \mathbf{r}_{\text{CM}})^T]. \quad (64)$$
The total torque acting on the body can be calculated as

$$T_B = \sum_{i=1}^{N} (\mathbf{r}_i - \mathbf{r}_{CM}) \times (m_p \mathbf{a}_i), \qquad (65)$$

where $\mathbf{a}_i$ is the acceleration of the point due to external forces.

The Newton-Euler equations (25a), (25b), (25c), and (25d) can be integrated using the Leapfrog algorithm as follows [27]. First, the acceleration of the center of mass $\mathbf{a}_{CM}$ and the angular acceleration of the rigid body $\boldsymbol{\alpha}_B$ are calculated as [27]

$$\mathbf{a}_{CM}(t) = \frac{\mathbf{F}_B(t)}{M_B} \qquad (66a)$$

$$\boldsymbol{\alpha}_B(t) = \mathbf{I}_B(t)^{-1} \cdot \left( \mathbf{T}_B(t) - \boldsymbol{\omega}_B \left( t - \frac{\Delta t}{2} \right) \times \mathbf{I}_B(t) \boldsymbol{\omega}_B \left( t - \frac{\Delta t}{2} \right) \right). \qquad (66b)$$

Then, the velocity of the center of mass $\mathbf{v}_{CM}$ and the angular velocity of the rigid body $\boldsymbol{\omega}_B$ are updated [27]

$$\mathbf{v}_{CM} \left( t + \frac{\Delta t}{2} \right) = \mathbf{v}_{CM} \left( t - \frac{\Delta t}{2} \right) + \mathbf{a}_{CM}(t) \Delta t \qquad (67a)$$

$$\boldsymbol{\omega}_B \left( t + \frac{\Delta t}{2} \right) = \boldsymbol{\omega}_B \left( t - \frac{\Delta t}{2} \right) + \boldsymbol{\alpha}_B(t) \Delta t. \qquad (67b)$$

Finally, the position of the center of mass is updated by [27]

$$\mathbf{r}_{CM}(t + \Delta t) = \mathbf{r}_{CM}(t) + \mathbf{v}_{CM}\left(t + \frac{\Delta t}{2}\right) \Delta t. \qquad (68)$$

The velocities and the positions of the rigid body points can be obtained as follows [27]

$$\mathbf{v}_i \left( t + \frac{\Delta t}{2} \right) = \mathbf{v}_{CM} \left( t + \frac{\Delta t}{2} \right) + \boldsymbol{\omega}_B \left( t + \frac{\Delta t}{2} \right) \times \left( \mathbf{r}_i(t) - \mathbf{r}_{CM}(t) \right) \qquad (69a)$$

$$\mathbf{r}_i(t + \Delta t) = \mathbf{r}_{CM}(t + \Delta t) + \mathbf{R}_B \left( t + \frac{\Delta t}{2} \right) \left( \mathbf{r}_i(t) - \mathbf{r}_{CM}(t) \right), \qquad (69b)$$

where $\mathbf{R}_B$ is the rotation matrix. It can be calculated as [27]

$$\mathbf{R}_B \left( t + \frac{\Delta t}{2} \right) = \mathbf{R}_x \left( \boldsymbol{\omega}_{B,x} \left( t + \frac{\Delta t}{2} \right) \Delta t \right) \cdot \mathbf{R}_y \left( \boldsymbol{\omega}_{B,y} \left( t + \frac{\Delta t}{2} \right) \Delta t \right) \qquad (70) \\ \cdot \mathbf{R}_z \left( \boldsymbol{\omega}_{B,z} \left( t + \frac{\Delta t}{2} \right) \Delta t \right),$$

where $\mathbf{R}_x(\theta)$, $\mathbf{R}_y(\theta)$, and $\mathbf{R}_z(\theta)$ are the matrices which yield the rotation about principal coordinate axes

$$\mathbf{R}_x(\theta) = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos(\theta) & -\sin(\theta) \\ 0 & \sin(\theta) & \cos(\theta) \end{pmatrix} \qquad (71)$$

$$\mathbf{R}_y(\theta) = \begin{pmatrix} \cos(\theta) & 0 & \sin(\theta) \\ 0 & 1 & 0 \\ -\sin(\theta) & 0 & \cos(\theta) \end{pmatrix} \qquad (72)$$

$$\mathbf{R}_z(\theta) = \begin{pmatrix} \cos(\theta) & -\sin(\theta) & 0 \\ \sin(\theta) & \cos(\theta) & 0 \\ 0 & 0 & 1 \end{pmatrix}. \qquad (73)$$

### 3.3. Boundary Conditions.
Interaction of the fluid phase described by SPH with solid objects can be obtained in various ways. Boundary forces taking into account viscous stress and fluid pressure are described in [49, 50]; coupling of SPH fluid with rigid solid bodies can be found in [27], while coupling with deformable solids is considered in [51]; coupled models where both the fluid and elastic solid are described in terms of SPH can be found in [28–30].

SPH particles belonging to two different phases or two different objects of the same phase are considered as boundary particles if the distance between them is smaller than smoothing length $h$. Interaction between these particles is defined by the boundary conditions described in Section 2.3.

In our model, we adopted the simplest way to satisfy the conditions (29) and (30) at the solid-viscous fluid interface by extending the summation in momentum conservation equations for all types of interacting particles [29]. For nonviscous fluid (the viscosity is equal to zero) a more complicated approach is necessary [27, 29]. Pressure values for solid phase particles can be obtained similarly to (4) by introduction of the speed of sound:

$$c_{\text{solid}} = \sqrt{\frac{K_s}{\rho_{0,\text{solid}}}}, \qquad (74)$$

where $K_s$ is the bulk modulus of the solid phase and $\rho_{0,\text{solid}}$ the initial solid phase density. $K_s$ for the elastic phase is defined by its elastic properties, and $K_s$ for the rigid phase is set 8 times larger than that for elastic phase. Using this approach, the no-slip boundary condition is imposed at the solid SPH particles.

Interaction between boundary particles belonging to different solid phases or objects can be modelled by introduction of pairwise interparticle forces based on Lennard-Jones potential [52], or repulsive forces modelling conservative full elastic interaction [27, 41]:

$$\mathbf{f}_i = -\frac{m_j}{m_i + m_j} B(\mathbf{r}_i) \mathbf{n}_j, \qquad (75)$$

where $\mathbf{f}_i$ is the force acting on particle $i$, $B$ is the symmetric function which defines the interaction strength, and $\mathbf{n}_i$ is
the local boundary normal. We have also successfully tested the interparticle interaction term based on the repulsive part of the pressure term from the momentum conservation equation.

3.4. XSPH Correction. In our model, we implement the XSPH correction which is a frequently used technique [25, 28, 30, 41] allowing the significant increase in the simulations stability and decrease in nonphysical noisy fluctuations. It consists in using the corrected velocity $\hat{v}_i$ for update of the SPH particles positions in the adopted numerical integration scheme instead of the usual particle velocity $v_i$ which is used in the momentum equation. The corrected velocity is calculated by [41]

$$ \hat{v}_i = v_i + \epsilon \sum_j m_j \frac{2(v_j - v_i)}{\rho_i + \rho_j}, \quad (76) $$

where $\epsilon$ is a constant parameter which varies from 0 to 1 and usually is set to 0.5. The summation in (76) is taken over the particles of the same phase. This modification makes particles move with velocities which are not too different for close particles, while the linear and angular momenta are not modified [41].

## 4. Validation Tests

4.1. Poiseuille-Couette Flow. The capacity of the model to represent correctly important physical phenomena was tested in various validation tests. A simple test to verify the capacity of the model to simulate fluid flow is the transient Poiseuille-Couette flow between two plane walls $x = 0$ and $x = L$. The fluid velocities imposed at these walls are $v(x = 0) = 0$ and $v(x = L) = V$, and the analytical solution for the velocity profile is given by [53]

$$ \begin{aligned} v(x,t) &= \frac{\rho g}{2\mu} x(L-x) + \frac{Vx}{L} \\ &\quad + \sum_{n=1}^{\infty} \left[ \frac{\rho g}{\mu} L^2 \left( \frac{2(-1)^n - 1}{(n\pi)^3} \right) - \frac{2V(-1)^{n+1}}{n\pi} \right] \quad (77) \\ &\quad \times \exp\left(-\frac{(n\pi)^2 \mu t}{\rho L^2}\right) \sin\left(\frac{n\pi x}{L}\right). \end{aligned} $$

For SPH fluid simulations, the cubic spline smoothing function is selected, which is most often used for fluid flow since it provides most accurate results. Fluid particles are initially placed at regular cubic lattice nodes with $5 \cdot 10^{-5}$ m spacing; the solid walls are represented in the same way by immobile SPH particles. The no-slip boundary conditions are applied at these solid phase particles. Spatially periodic boundary conditions are imposed along y- and z-axis. The smoothing length is $h = 10^{-4}$ m, the distance between the walls is $L = 10^{-3}$ m, the body force is $g = 10^{-3}$ m/s², and the imposed velocity is $V = 2 \cdot 10^{-4}$ m/s. The fluid density and

<figure><img src="image_3.png" /><figcaption>FIGURE 2: Transient Poiseuille-Couette flow. Velocity profile evolution with time. Velocity profiles correspond to $t = 10^{-2}$ s, $t = 4 \cdot 10^{-2}$ s, $t = 8 \cdot 10^{-2}$ s, $t = 16 \cdot 10^{-2}$ s, and $t = 1$ s.</figcaption></figure>

viscosity are $\rho = 1000$ kg/m³ and $\mu = 10^{-3}$ Pa·s, respectively. In all cases shown in this paper, small values of body force are imposed following [53].

The obtained velocity profiles at different times are compared with analytical solutions in Figure 2. The agreement between the numerical and analytical solution is very good. The maximal velocity error was smaller than 0.4%. It must be noted that the accuracy of the simulation depends on the smoothing function and the way chosen to discretize the components of the Navier-Stokes equation, especially the viscous dissipation term. The best precision was obtained for the cubic spline and the viscous dissipation term discretization as given in (45).

4.2. Periodic Arrays of Spheres. The capacity of the SPH to simulate fluid flow in porous media can be tested for a nontrivial geometry represented by periodic array of spheres for which an analytical expression of permeability is available [54, 55].

The permeability $K$ of a porous medium can be found from Darcy's law

$$ \langle \mathbf{v} \rangle = -\frac{K}{\mu} (\nabla p - \rho \mathbf{g}), \quad (78) $$

where $\mu$ is the viscosity, $\langle \mathbf{v} \rangle$ is the mean fluid seepage velocity, $\nabla p$ is the pressure gradient, and $\mathbf{g}$ is the body force. For a spatially periodic simple cubic array of spheres, the permeability can be found analytically [54, 55] as

$$ K_a = \frac{\tau_0}{6\pi R} k^*, \quad (79) $$
where $\tau_0$ is the volume of the unit cell, R the radius of the sphere, and $k^*$ the inverse of the dimensionless force exerted by the fluid flow on the array sphere [55]:

$$k^* = \frac{6\pi\mu R\nu}{F}. \qquad (80)$$

Values of $k^*$ can be found in [54, 55] for various values of solid volume fraction $\phi$

$$\phi = \frac{V_{\text{solid}}}{V_{\text{cell}}} = \frac{4\pi R^3}{3L^3}. \qquad (81)$$

The ability of SPH to simulate fluid flow for periodic array of spheres with a good accuracy was demonstrated in [53] for accurately discretized geometries; about 30 SPH particles across the pore throat are necessary in order to obtain the results within about 2% of analytical solution. In our case, we are interested in estimation of possible differences between the analytical solution and the simulation result $K_{\text{num}}$ for relatively coarse geometry and particle (i.e., the number of particles along the smoothing length $h$) discretizations, since for the clay samples images obtained by the 3D imaging methods, the overall sample size is too big to discretize everywhere the pore space with a large number of SPH particles.

For numerical simulation the unit cell size is $L = 20h$, where $h = 10^{-4}$ m is the smoothing length. The center of the solid sphere $\mathbf{r}_s$ is located at the center of the unit cell. Fluid and solid SPH particles are initially placed at regular cubic lattice nodes with $h/2$ interparticle spacing. The no-slip boundary conditions are applied at the middle distance between the solid and fluid particles, while spatially periodic boundary conditions are imposed at the sides of the unit cell. In order to mimic a macroscopic pressure gradient, the flow is driven by the imposed body force $\mathbf{g} = (0, 0, 5 * 10^{-3}) \text{ m/s}^2$, the fluid density is $\rho = 1000 \text{ kg/m}^3$, and viscosity is $\mu = 10^{-3} \text{ Pa.s}$. The speed of sound is $c_s = 0.1 \text{ m/s}$. High values of speed of sound result in very small time steps which slow down the simulation. Therefore, the common approach when using SPH is to adopt smaller values of speed of sound [44] keeping in mind that they should be large enough to ensure the relative fluid density variation inferior to 1% during the simulation.

In order to test the sensitivity of the method to various $\phi$ values, three values of solid concentration are considered $\phi = 0.343$, $\phi = 0.45$, and $\phi = 0.5236$. Coarse solid spheres are discretized with elementary cubes of size $l_c = h$; that is, each elementary cube contains 8 solid SPH particles. The elementary cube is considered solid if its center $\mathbf{r}_c$ satisfies the condition $\|\mathbf{r}_c - \mathbf{r}_s\| \le R$. In order to test the influence of the sphere surface discretization, spheres are also discretized with solid SPH particles which satisfy the condition $\|\mathbf{r}_p - \mathbf{r}_s\| \le R$, where $\mathbf{r}_p$ is the location of the particle. This way provides a refined discretization and equivalent to discretizing the sphere with elementary solid cubes of size $l_c = 0.5h$, each of which contains a single SPH solid particle located in its center. The difference between the coarse and refined discretizations is displayed in Figure 3 for $\phi = 0.45$. The sensitivity of the method to particle discretization is tested by performing the same series of simulations for the $h/3$ interparticle spacing.

Numerical errors for the calculated permeability values are summarized in Table 2. For the $h/2$ interparticle spacing, errors of order of 30% are observed with little influence of the sphere surface discretization, which means that precision is mainly affected by the coarse discretization of flow equations in terms of SPH. For the $h/3$ interparticle spacing, the errors due to particle discretization are less significant, and clear influence of the sphere surface discretization can be observed. The errors of order of 20% are observed for the coarse sphere surface discretization, while for the refined sphere surface discretization the error is less than 10%. As mentioned at the beginning of this section, the pore throat must be very well discretized in order to obtain the permeability value within few percentages of analytical solution when using SPH.

It must be noted (see Table 3) that the precision of the simulation results depends strongly on the number of SPH particles within the smoothing length, which is closely related to the initial interparticle spacing. The larger number of particles provides a better discretization of equations in terms of SPH and consequently more precise results at the expense of significantly increased computational cost. The influence of the number of SPH particles within the smoothing length on the simulation precision can be stronger than accuracy of geometry discretization. Thus, the permeability calculated for $\phi = 0.343$, $L = 40h$, $h/3$ setup is only 7.02% accurate compared to 3.97% accuracy for $\phi = 0.343$, $L = 20h$, $h/4$ setup, in spite of larger total number of SPH particles and better sphere surface discretization.

Based on these simulation results, we can expect for our drainage simulations the precision of about 30%, which is satisfactory for simulations performed for coarse geometries with $h/2$ interparticle spacing.

### 4.3. Two-Phase Poiseuille Flow.
In order to test the capacity of the model to simulate multiphase flow, let us consider a two-phase (liquid-gas) Poiseuille flow between parallel plane walls located at $x = -l$ and $x = l$. The liquid-gas interface is also a plane located at $x = 0$. The gas is located in $x < 0$ and the liquid in $x > 0$ region, respectively. The analytical solution for fluid velocities is the following:

$$v_g(x, t) = -\frac{Fl^2}{2\mu_g} \left[ \frac{x^2}{l^2} - \frac{x}{l} \frac{\mu_g - \mu_f}{\mu_g + \mu_f} - \frac{2\mu_g}{\mu_g + \mu_f} \right], \qquad (82a) \\ -h \le x \le 0$$

$$v_f(x, t) = -\frac{Fl^2}{2\mu_f} \left[ \frac{x^2}{l^2} - \frac{x}{l} \frac{\mu_g - \mu_f}{\mu_g + \mu_f} - \frac{2\mu_f}{\mu_g + \mu_f} \right], \qquad (82b) \\ 0 \le x \le h,$$

where $F$ is a body force or a pressure gradient.

The simulation domain is the same as for the considered above transient Poiseuille-Couette flow simulation. The no-slip boundary conditions are applied at the solid phase SPH particles, while spatially periodic boundary conditions are imposed along $y$- and $z$-axis. The cubic spline smoothing function is used with smoothing length $h = 10^{-4}$ m, the gas/liquid region width is $l = 5 \cdot 10^{-5}$ m, and the body force is equal
<figure><img src="image_4.png" /><figcaption>FIGURE 3: Coarse and refined discretizations of the sphere for $\phi = 0.45$ (a central cut).</figcaption></figure>

<table><caption>TABLE 2: Numerical error for permeability of periodic arrays of spheres for various $\phi$ values, solid sphere discretizations, and interparticle spacing as compared to $K_a$.</caption><thead><tr><th>φ</th><th>Discretization</th><th>K<sub>a</sub>, m<sup>2</sup></th><th>err, % (spacing h/2)</th><th>err, % (spacing h/3)</th></tr></thead><tbody><tr><td>0.343</td><td>Coarse</td><td>3.1728 * 10<sup>-8</sup></td><td>-37.47</td><td>-21.26</td></tr><tr><td>0.45</td><td>Coarse</td><td>1.5904 * 10<sup>-8</sup></td><td>-34.19</td><td>-17.22</td></tr><tr><td>0.5236</td><td>Coarse</td><td>1.0107 * 10<sup>-8</sup></td><td>-37.47</td><td>-21.31</td></tr><tr><td>0.343</td><td>Refined</td><td>3.1728 * 10<sup>-8</sup></td><td>-31.03</td><td>-8.32</td></tr><tr><td>0.45</td><td>Refined</td><td>1.5904 * 10<sup>-8</sup></td><td>-30.11</td><td>-9.24</td></tr><tr><td>0.5236</td><td>Refined</td><td>1.0107 * 10<sup>-8</sup></td><td>-33.66</td><td>-9.80</td></tr></tbody></table>

<table><caption>TABLE 3: The calculated permeability $K_{\text{num}}$ of periodic array of spheres for $\phi = 0.343$, $L = 20h$ with refined discretization.</caption><thead><tr><th>Interparticle spacing</th><th>K<sub>num</sub>, m<sup>2</sup></th><th>err, %.</th></tr></thead><tbody><tr><td>h/2</td><td>1.1115 * 10<sup>-8</sup></td><td>-30.11</td></tr><tr><td>h/3</td><td>1.4434 * 10<sup>-8</sup></td><td>-9.24</td></tr><tr><td>h/4</td><td>1.5272 * 10<sup>-8</sup></td><td>-3.97</td></tr></tbody></table>

to $F = 10^{-1}$ Pa/m. The liquid density is $\rho_f = 10^3$ kg/m$^3$ and viscosity $\mu_f = 10^{-2}$ Pa·s; the gas density is $\rho_g = 1$ kg/m$^3$ and viscosity $\mu_g = 10^{-4}$ Pa·s.

The numerical velocity profile is compared to the analytical solution in Figure 4. A very good agreement is obtained with maximal velocity deviation smaller than 1%. Several simulations performed for even larger ratios of water/gas physical properties also demonstrate a very good accuracy, which proves the capability of the SPH model to simulate accurately multiphase mixtures with highly contrasting physical properties.

**4.4. Laplace Law.** In order to test the continuum force algorithm described in Section 2.1.1, the compliance with Laplace's law for spherical bubbles is tested. Simulation is initiated with a cube of size $L$ consisting of gas SPH particles,

placed at the center of the cubic unit cell of size $L_{\text{cell}} = 46h$ and surrounded by fluid SPH particles. Under the influence of surface tension forces, the gas cube is transformed into a sphere of radius $R = L(3/4\pi)^{1/3}$. After the equilibrium state has been obtained, the pressure difference $\Delta p = p_g - p_f$ across the fluid-gas interface is calculated and plotted for different sphere radii in Figure 5. The difference between the numerical and analytical pressure drop across the interface does not exceed 2% even for the smallest bubbles.

**4.5. Reaction of Elastic Phase to Fluid Pressure.** The elastic-fluid coupling can be tested by calculation of the relative volume variation $dV/V$ of an elastic cube due to the pressure step $dP$ applied by the surrounding fluid:

$$K_s = -\frac{dP}{dV/V} \qquad (83)$$

The volume variations for several pressure steps are compared with analytical solution given by (83) in Figure 6. Relative error varies from 0.33% for small volume variations to 11.4% for relatively large volume variations. This error is mostly due to the fact that the adopted approximation for second-order derivative (61) in elastic force calculation generates errors in the range of 5–10% for large deformations $\epsilon \sim 0.1$. However, from Table 1 it follows that clayey matrix
<figure><img src="image_5.png" /><figcaption>FIGURE 4: Two-phase Poiseuille flow. Analytical solution (black solid line) for velocity profile is compared to the SPH numerical solution (points).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5: Laplace's law for spherical bubbles. Pressure step versus $h/R$.</figcaption></figure>

elastic material fails at much lower strains which limits the influence of second-order derivative errors on the simulation results when applied to an argillite sample.

# 5. Application to Clay Sample

In order to test the ability of our code to simulate gas-water flow inside a realistic pore space in repository conditions, we make a very first application to a small sample of Callovo-Oxfordian (COx) clay.

<figure><img src="image_7.png" /><figcaption>FIGURE 6: Elastic cube volume change due to applied fluid pressure.</figcaption></figure>

## 5.1. Characterization of the Sample.
The COx claystone sample EST27405-1 is represented by a series of 103 images of 771 × 739 pixels size corresponding to consecutive scans orthogonal to z-axis obtained by FIB method [56]. The 3D voxel size is 8.49×10.78×10 nm³. The images are treated based on gray level in order to recognize the pore-space and clayey matrix. The porosity of the sample is $\epsilon = 0.0307 \pm 0.0039$. A more complex morphological analysis was also applied to the sample in order to distinguish various solid inclusions which are most likely small grains of calcite and quartz [56]. However, for this first approach, in the following we will consider the nonporous part of the sample as a continuous elastic phase without rigid inclusions. This simplified approach allows us to use macroscopic mechanical properties of COx rock for solid matrix parameterization which are summarized in Table 1.

In order to decrease the computational costs, a coarsening procedure was applied to the original sample consisting in replacing of 8 neighbouring voxels with a single two times larger voxel whose nature (pore space or clayey matrix) was defined based on the majority rule. Two samples of 385 × 369 × 51 and 192 × 184 × 25 size were obtained by applying this procedure once and twice, respectively. The pore space of the sample consists of several isolated pores. The percolating pores can be assimilated to fissures which connect opposite sides of the sample along the z-axis. They are surrounded by multiple small disconnected pores. The percolating pore components are displayed in Figure 7. As it was mentioned before, most of the pores inside the claystone are of the nanometer size; therefore, some connectivity information can be missing because of insufficient resolution of FIB imaging.

Due to significant computational costs (one simulation for 54h×54h×70h size unit cell requires more than 3 million time steps and lasts about 30 days on our basic GPU card), only one percolating pore component was considered (see
<figure><img src="image_8.png" /><figcaption>FIGURE 7: Largest percolating pore components inside the Callovo-Oxfordian clay sample.</figcaption></figure>

Figure 7). Since the individual shape of the selected percolating component influences the details of the drainage, the obtained results cannot be generalized for all the percolating pore components; however, we expect that they will provide an insight into the phenomena which may occur in these conditions. For simulations we considered a percolating pore (the blue pore in Figure 7) extracted from once (385×369×51) and twice (192 × 184 × 25) coarsened samples. The chosen pore has a fracture like shape with larger entry at the top side and smaller outlet at the bottom side. The results obtained for the percolating pore (see Figure 8(a)) extracted from the sample coarsened twice are very noisy due to insufficient particle resolution; therefore only the results for less coarsened sample will be discussed here.

For the present study, the OpenMP [57] and CUDA GPU parallelization techniques were applied with a significant gain in numerical performance. One drainage simulation as presented further in the chapter took about 30 days on NVIDIA Quadro K5200 standard graphic card, which is about 52 times faster than a single thread processor version. However, this performance could be greatly improved (speed up of about 12 times) by using dedicated GPU units like Nvidia Maxwell P100.

### 5.2. Numerical Simulations.
Several drainage tests are performed in order to study the model behavior under various conditions. A unit cell for simulations is obtained by cutting a part of the sample containing the selected percolating pore component as well as small surrounding disconnected pores. In order to avoid fluid movement across the cell borders, an elastic layer is added to each lateral side of the sample. Finally, the sample is placed between horizontal solid walls in order to prevent its spatial movement during the simulation as displayed in Figure 8(b). Pores inside the sample are initially filled with water phase, while the gas phase is located in the reservoir outside the sample. The height of the reservoir is $H_r = 20h$. The volume of gas is sufficient to completely drain the sample. Spatially periodic boundary conditions are applied to the unit cell along the z-axis. A body force (0, 0, $F_z$) is applied to the gas in the reservoir outside the sample, which creates a pressure difference between upper $z = 60$ ($p_g = H_r F_z \rho_g$) and lower $z = 10$ ($p_g = 0$) sides of the sample (see Figure 8(b)) and causes drainage. Finally, the confining stress is applied to the lateral sides of the sample to study its impact on gas migration and damage of the clay. A cross-section of the unit cell of $54h \times 54h \times 70h$ size is given in Figure 8(b). One can see that the gas phase (green) penetrates the sample

<table><caption>TABLE 4: Dimensionless numbers characterizing two-phase water-gas flow in clay microfractures.</caption><thead><tr><th>Number</th><th>Re</th><th><i>M</i></th><th>Ca</th><th>Bo</th></tr></thead><tbody><tr><td>Expected value</td><td>4 * 10<sup>-5</sup></td><td>1.1 * 10<sup>2</sup></td><td>7.2 * 10<sup>-10</sup></td><td>4.7 * 10<sup>-6</sup></td></tr><tr><td>Simulation value</td><td>≈ 0.5</td><td>1.1 * 10<sup>2</sup></td><td>≈ 1.25 * 10<sup>-4</sup></td><td>0</td></tr></tbody></table>

and liquid phase (yellow) is leaving the sample from through the bottom. The entry and exit pore throats are not visible since they are located in different cross-sections than the one displayed in Figure 8(b).

Usually, it may be possible to prescribe real values of physical parameters of all the phases considered in the model and to solve numerically the corresponding equations. However, in this case the time step would be prohibitively small ($\Delta t < 10^{-12}$ s) due to the explicit integration scheme adopted in the model. Also, the time to completely drain the sample can be very long. For this reason, the physical parameters values are selected based on the dimensionless numbers which define the dominant physical phenomena.

Two-phase flow in porous media can be characterized by several dimensionless numbers:

$$Re = \frac{\rho_g V_g L}{\mu_g} \qquad (84a)$$

$$M = \frac{\mu_w}{\mu_g} \qquad (84b)$$

$$Ca = \frac{\mu_g V_g}{\sigma} \qquad (84c)$$

$$Bo = \frac{(\rho_w - \rho_g) g L^2}{\sigma} \qquad (84d)$$

The Reynolds number $Re$ which represents the ratio between the inertial and viscous forces is very small $Re \ll 1$ for gas flow in clay due to small flow velocity $V_g$ and small characteristic length $L$ which is typically the pore size. The mobility number $M$ gives the ratio between water and gas viscosities. The capillary number $Ca$ measures the ratio of the viscous forces to the surface tension forces and is also usually very small $Ca \ll 1$ for the considered situation. The Bond number $Bo$ characterizes the ratio between external forces (pressure gradient, gravitation, and body forces) and surface tension and describes the shape of bubbles. The dimensionless numbers values expected for the gas-water flow in real clay microfractures [17] and the values applied in our simulations are presented in Table 4. It must be noted that the $Re$ and $Bo$ numbers given in [17] were calculated for micrometer pore sizes. However, the direct extrapolation to nanometer size pores would lead to still smaller numbers, but it does not change the flow regime.

It can be observed that the Reynolds and capillary numbers observed in our simulations (for mean gas velocity $V_g \approx 5 * 10^{-5}$ m/s before the percolation) are several orders of magnitude higher than the expected real ones as a consequence of significantly reduced viscosity ($\mu_w = 10^{-7}$ Pa·s), $\mu_g = 10^{-9}$ Pa·s) and surface tension ($\sigma = 4 * 10^{-10}$ N·m·s⁻¹)
<figure><img src="image_9.png" /><figcaption>FIGURE 8: (a) A disconnected pore component $\epsilon = 0.0017$ extracted from 192 × 184 × 25 sample. (b) A cross-section ($38h < x < 39h$) of the unit cell obtained from 385 × 369 × 51 sample during drainage simulation. Green points correspond to gas particles, yellow points to water particles, and black points to rigid phase particles; other colors reflect relative variation of stress inside the elastic phase.</figcaption></figure>

values; however, they are still sufficiently small and should describe adequately laminar two-phase flow dominated by capillary forces.

As it was previously mentioned, in weakly compressible formulation the speed of sound in the fluid must be high enough to effectively penalize density variations and at the same time small enough to allow reasonably large time steps. The adopted value $c_s = 0.5 \text{ m·s}^{-1}$ for both fluid phases satisfies these conditions with relative fluid density variation of order of 0.12% observed in our simulations for the pressure difference $\Delta P = 0.3 \text{ Pa}$ imposed between sample sides. In order to accelerate the simulations and to more easily observe the effects of damage, the gas pressure at the inlet is set to $P_{\text{gas,in}} = 15 \text{ MPa}$ which is slightly higher than the maximal pressure which could be expected within the repository in COx (the vertical confinement is of about 12 MPa and the gas entry pressure to the sound argillite may be as high as 10 MPa [58]). Simulations with smaller pressure values were conducted but will not be presented in this paper. So, the ratio between $P_{\text{gas,in}}$ and simulation $P_{\text{gas,in}} = 0.3 \text{ Pa}$ pressure values at the sample entry is $5 * 10^7$. Since in these conditions the fluid pressure values calculated by (4) are much smaller than those expected in real conditions, the elastic properties from Table 1 are scaled so that the ratio $R_p$ between the gas pressure at the sample entry $P_{\text{gas,in}} = 0.3 \text{ Pa}$ and Young's modulus $E = 60 \text{ Pa}$ of the elastic phase used in the simulation

$$R_p = \frac{P_{\text{gas,in}}}{E} \qquad (85)$$

stays the same as in repository conditions where $R_p \approx 0.005$. For the sake of better understanding with respect to the aimed application, we are going to comment all simulations by

scaling the numerical results back to the repository conditions values (using the factor $5 * 10^7$). These rescaled values will be preceded by $\sim$.

### 5.2.1. Drainage with and without Damage.
The pressure difference $\Delta P \sim 15 \text{ MPa}$ is imposed between upper and lower sample sides, while the confining stress of $\sigma_{\text{conf}} \sim 6 \text{ MPa}$ is applied to the lateral sides. The drainage tests are performed with and without damage model in the elastic phase. Another test is performed with a completely incompressible rigid porous matrix. Due to the pressure difference between the opposite sides, gas phase enters the sample and pushes the water phase out. Since the gas phase is nonwetting, it propagates mostly inside the wetting water phase.

Since the ratio $R_p = 0.005$ is very small, only small local deformations will be observed in the intact elastic matrix; however, when accumulated through the sample, these deformations will result in faster drainage compared to the completely incompressible rigid porous matrix as it can be observed in Figure 9(a) where the evolution of saturation with time is shown. It can be observed in (Figure 9(a)) that at the very beginning the gas enters the sample with a larger velocity regardless the application of the damage model. This is caused by the establishment of the pressure gradient through the sample, which causes deformation of the elastic matrix increasing the available pore space. The effect is more pronounced for the model with damage due to the softening of the elastic matrix. The compressibility of the fluid phases can also cause the similar effect if the speed of sound is not high enough.

After the initial compression phase, the saturation evolves almost linearly until the percolation, where the curve changes the slope and the drainage continues with a decreasing
<figure><img src="image_10.png" /><figcaption>FIGURE 9: Liquid $S_w$ (broken lines) and gas $S_g$ (solid lines) saturations evolution for simulations varying from rigid matrix to confined elastic matrix with damage. (a) Rigid solid matrix (cyan) and elastic matrix under confining stress $\sigma_{conf} \sim 6$ MPa with (magenta) and without (black) damage applied to the elastic phase. (b) Elastic matrix with application of damage model under confining stresses of $\sigma_{conf} \sim 10$ MPa (red), $\sigma_{conf} \sim 8$ MPa (green), $\sigma_{conf} \sim 6$ MPa (magenta), and $\sigma_{conf} \sim 4$ MPa (blue).</figcaption></figure>

velocity due to the fact that water phase is now mainly located near the walls. It can be expected that the sample will be completely drained after a longer time due to viscous drag force. Incorporation of the damage model (Figure 9(a)) increases the drainage rate at all the stages of simulation and the nonwetting fluid percolation time is accelerated (from $8.74 \times 10^{-3}$ s to $7.58 \times 10^{-3}$ s). However, for this particular geometry of pores, which are percolating already in the initial sample, this is the most important change in behaviour, which stays qualitatively similar with or without damage.

### 5.2.2. Drainage under Different Confining Pressures.
The next series of drainage tests is performed for various confining stresses $\sigma_{conf} \sim 10$ MPa, 8 MPa, 6 MPa, and 4 MPa applied to lateral sides and for the same pressure difference $\Delta P \sim 15$ MPa.

Same as for (Figure 9(a)), the initial compression phase can be observed in Figure 9(b). It is more pronounced for $\sigma_{conf} \sim 4$ MPa where the largest deformation of the elastic matrix by the entering gas is expected due to the smallest applied confining stress. What is more, in this case we note the largest number of damaged particles (Figure 10) causing the elastic matrix softening. As can be seen in Figure 9(b), these effects are largely attenuated with increasing confining stress.

The percolation times and gas saturation at percolation are summarized in Table 5. As expected, the percolation occurs earlier for smaller confining stresses due to the larger deformations and as a consequence due to increased permeability that facilitates drainage. The largest $t_{perc}$ corresponds

to the incompressible porous matrix. The gas saturation at percolation is close to 0.5 for all the cases and increases slightly with confining stress.

Evolution of the number of damaged points and of the total damage is presented in Figure 10. After an initial rapid increase, number of damaged points continues to grow with a reduced rate for all considered confining stresses. A slightly faster increase can be observed just before percolation which is due to the increase of the water pressure (see Figure 13). After percolation the number of damaged points stays almost constant with a very slow grow rate. The number of damaged points is larger for small confining stresses than for high confining stresses since it is more difficult to fracture compressed elastic matrix which results in slower drainage in the last case. Initially damage occurs at the edgy parts of the elastic-fluid interface where the most of the traction stress is accumulated; then it propagates through the elastic phase in the direction of the sample sides (see Figure 11). Also a lot of damaged elastic points were observed in a region where the pore space approaches the sample lateral sides. The total accumulated damage presented in Figure 10(b) follows the same pattern as the number of damaged particles before the percolation, after which it decreases rapidly to a new relatively stable level as a consequence of fast relaxation of gas and water pressures after the percolation (see Figure 13(a)). As it can be observed in Figures 11(c) and 11(d), this reduction of pressure causes also a decrease of number of damaged elastic particles: first particles to restore their original elastic properties are those located the furthest from the damage departure location.
<table><caption>TABLE 5: Percolation times $t_{\text{perc}}$ and gas saturation $S_{g,\text{perc}}$ at percolation for various confining stresses.</caption><thead><tr><th>σ<sub>conf</sub>, ~MPa</th><th>t<sub>perc</sub>, s</th><th>S<sub>g,perc</sub></th></tr></thead><tbody><tr><td>10</td><td>9.18 * 10<sup>-3</sup></td><td>0.5158</td></tr><tr><td>8</td><td>8.52 * 10<sup>-3</sup></td><td>0.5074</td></tr><tr><td>6</td><td>7.58 * 10<sup>-3</sup></td><td>0.5043</td></tr><tr><td>4</td><td>6.56 * 10<sup>-3</sup></td><td>0.4971</td></tr><tr><td>6, no damage</td><td>8.74 * 10<sup>-3</sup></td><td>0.5039</td></tr><tr><td>Incompressible porous matrix</td><td>9.84 * 10<sup>-3</sup></td><td>0.5159</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>FIGURE 10: (a) Number of damaged elastic SPH particles. (b) The total damage calculated as a sum over all the damaged points.</figcaption></figure>

The seepage velocities of fluid phases during drainage are plotted in Figure 12. At the initial stage, the velocities are high due to multiple factors mentioned above. During the linear stage of drainage, the velocities are almost constant and close to each other for all the considered cases. Water velocity gradually decreases before the percolation, since the water which is left in the sample is located close to the pore walls where its velocity tends to zero. After the percolation it stabilizes at a significantly lower level. As expected, gas velocity increases several times after the percolation.

The water $\langle P_w \rangle$ and gas $\langle P_g \rangle$ pressures averaged over the percolating pore volume and the capillary pressure $P_c = \langle P_g \rangle - \langle P_w \rangle$ are presented in Figure 13. The water phase pressure increases rapidly until the establishment of the pressure gradient through the sample; then it grows gradually until the percolation after which it drops and stabilizes at a new level. The gas phase pressure decreases at the initial stage of drainage; then it stays relatively constant till the percolation after which it drops and stabilizes at a new level similarly to the water phase. In some cases, the drop of the pressure can lead to reclosure of the percolating pore passes opened during the gas passage. Capillary pressure decrease with growing gas saturation is related to the geometry of the percolating

pore, in which aperture increases from entry to the central part where it reaches the maximum value. Capillary pressure curves are almost identical for the considered cases starting from $S_g = 0.1$. Just before the percolation the capillary pressure decreases faster than at the earlier stage and stabilizes after the percolation. This decrease is related to the growth of the water pressure which is caused by the gas pushing the water under low pressure from the bottom part of the sample. The results for the smallest confining stress are quite noisy due to important number of damaged elastic points which can cause significant pore surface displacements influencing the stability of the numerical method.

It must be noted that the numerically obtained capillary pressure values displayed in Figure 13 are about 3.6 times smaller than those which can be expected for this pore size in the case of water/air interface with $\sigma = 7.2 * 10^{-2} \text{ N} \cdot \text{m}^{-1}$, since water/air surface tension value was decreased 3.6 times in simulations, in order to speed them up.

### 5.2.3. Damage Model versus Numerical Fracturing.
From the presented simulations it can be concluded that the mechanical interaction between the solid matrix and the fluid phases can influence the drainage process in clays in repository
<figure><img src="image_12.png" /><figcaption>FIGURE 11: A cross-section (SPH points inside the $45h < z < 46h$ range) of the unit cell under confining stress $\sigma_{conf} \sim 8$ MPa. Green points correspond to gas particles, yellow points to water particles, black points to rigid phase particles, cyan points to intact elastic particles, and other colors to elastic particles damaged to various degree with magenta corresponding to completely damaged particles. Images correspond to (a) $t = 4$ s, (b) $t = 52$ s, (c) $t = 80$ s, and (d) $t = 92$ s.</figcaption></figure>

conditions. It has been seen as well that the resulting damaged area was relatively wide and did not localize into new fractures. We have tested several approaches in order to better localize damaged particles. We observed that the details of damage model greatly influenced simulation outcome. The results presented in previous chapters were obtained supposing the reversibility of damage variable; that is, elastic properties of damaged particles can be restored if this is indicated by the damage criterion. Also, in order to preserve the volume of the elastic phase, damaged elastic particles

continue to interact by pressure gradient term (second term in (15)). Under these conditions, the number of damaged elastic particles was bigger for smaller confining stresses, and these particles formed connected regions as it can be seen in Figure 11. These regions can be considered as containing microfractures, which present apertures too small to appear at the model scale. For these reasons, the fluids particles cannot penetrate in these microfractures and create new connections between pores. However, these passages can be created if the fluid pressure force acting on the pore space
<figure><img src="image_13.png" /><figcaption>FIGURE 12: Water (a) and gas (b) seepage velocities during drainage as functions of time for various confining stresses.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 13: (a) Water and gas mean pressures inside the sample as functions of time. (b) Capillary pressure as a function of gas saturation.</figcaption></figure>

surface exceeds the confining stress as illustrated in Figure 14(b) where a formation of new connections visible at the model scale is shown. This crack is big enough to be penetrated by gas particles. It can be observed that the pore space is noticeably larger than at the initial time (see Figure 14(a)).

When using SPH model for the elastic phase, a fracturing can be observed even without applying a damage model. For the simulations presented in Sections 5.2.1 and 5.2.2 the ratio $R_p = 0.005$ is very small, which means that very small local deformations will be observed in the intact elastic matrix. Therefore, without application of damage model, there would
<figure><img src="image_15.png" /><figcaption>FIGURE 14: A cross-section (SPH points inside the $58h < z < 59h$ range) of the unit cell. Color notation is the same as for previous figure showing a cross-section. Images correspond to (a) Initial sample surface, $t = 0$ s. (b) No confining stress is applied. $t = 4$ s. (c) Constant volume boundary conditions. No damage model is applied, and the $R_p = 0.375$. Time $t = 3.6$ s. (d) Constant volume boundary conditions. Damaged particles are considered as ghost particles which do not interact with any other particles. Blue particles correspond to gas phase. $t = 6.8$ s.</figcaption></figure>

be no initiation of new fractures. However, for larger values, for example, $R_p = 0.375$, the deformations may be big enough to separate the SPH particles at distance larger than $h$, which can be interpreted as a creation of a new fracture at the model scale. An example of such a simulation, without application of a damage model can be seen in Figure 14(c). This phenomenon is usually called the numerical fracturing.

It allows in some conditions to localize nicely new fractures; however, it is impossible to parameterize this process in order to fit it to a particular porous material.

A still different behaviour can be observed if damaged elastic particles are considered as ghost particles, which no longer interact with any other particles. In this case, the release of the stress accumulated inside the elastic matrix
<figure><img src="image_16.png" /><figcaption>FIGURE 15: A cross-section (a) $24h < z < 25h$ and (b) $33h < x < 34h$ of the unit cell under confining stress $\sigma_{conf} \sim 6$ MPa at $t = 12$ s. Two spherical rigid inclusions are incorporated into initially elastic solid matrix.</figcaption></figure>

may result in formation of narrow fractures as displayed in Figure 14(d). However, the total volume of the elastic phase is no longer conserved in these conditions.

5.2.4. Solid Inclusions. In order to consider the possible influ-
ence of the rigid inclusions on the flow conditions, two rigid
spheres were incorporated into the elastic matrix as displayed
in Figure 15. The only implemented interaction between solid
and elastic particles was a pressure gradient repulsion. All
other materials parameters were the same as in Section 5.2.2.
It can be observed (see Figure 15(a)) that if the sphere is not
deeply incorporated into the elastic matrix, its position can
be altered by the flow damaging the surrounding elastic
particles, and creating new narrow passages between the rigid
inclusion and the elastic matrix which can be penetrated by
fluid phase particles. For deeply incorporated rigid inclusion
(see Figure 15(b)), very little influence on the surrounding
elastic solid is observed. It can be concluded that the presence
of rigid inclusions in the model provides, as expected, a
different way of creating new paths for gas migration. In
order to obtain a better description of the interaction between
rigid inclusions and elastic matrix, friction forces must be
included at the interface which will decrease to a certain
extent the tangential movement of rigid inclusions in contact
with elastic matrix.

**6. Conclusions**

The aim of this paper was to propose a modelling approach
for the dilatant two-phase flow which is assumed to occur
during gas migration within initially water saturated indurat-
ed clay rock. This type of flow, unlike classical viscoelastic
flows, require an interaction of flow with a deformable solid
matrix.

We developed a numerical model which allows simulat-
ing a two-phase fluid flow inside a porous medium and taking
into account the mechanics of the solid matrix, possibly
composed of two different media: a linearly elastic continuous
medium and isolated rigid inclusions. Then a damage model
is introduced in order to take into account in a controlled
manner the degradation of elastic material properties due to
local stress-strain state. The model can be applied to the FIB-
SEM images of the real clay samples in order not only to
study various phenomena at pore scale but also to calculate
effective transport properties, that is, permeability for varying
pressure gradients and confining stresses taking into account
deformation and degradation of the elastic matrix.

The model was solved in the framework of Smoothed
Particle Hydrodynamics, which is a Lagrangian mesh-free
method. Classical methods to numerically simulate fluid flow
are mostly based on Finite Volume Methods (FVM), while
in solid mechanics the Finite Element Methods (FEM) are
widely applied. These methods allow solving several problems
more efficiently than SPH; however, in some complicated sit-
uations they may suffer from drawbacks related to their mesh
based nature. This includes the cases with deformable and
moving interfaces, dealing with large deformations and prop-
agation of fractures. Use of SPH model is very appropriate
in this context since it provides the opportunity to describe
various types of interactions between different materials
within the same framework. The presented SPH model was
coded using C++ programming language implementing all
the phenomena mentioned above. The SPH is computation-
ally expensive. The main load comes from calculation of
neighbouring particles contributions when evaluating densi-
ties and forces. Its cost grows significantly with the number
of particles located within the smoothing length. Therefore,
the original code was parallelized using OpenMP [57] and
CUDA GPU techniques. Parallelization of SPH can be per-
formed relatively easily, since densities and forces are calcu-
lated by summing contributions from neighbouring particles;
therefore, each particle can be treated separately by a compu-
tational thread.

For the present study, the OpenMP and CUDA GPU par-
allelization techniques were applied with a significant gain in
numerical performance when run on an ordinary graphical
card.

With several elementary tests we have shown that the
model is sufficiently stable and robust for use in conditions
representative for application to gas-water flow in argillites
and that a limited resolution was sufficient to obtain perme-
ability values with a reasonable precision (about 30%).

The model was then applied to simulate drainage in a
small 3D sample obtained from FIB reconstruction of COx
claystone. The simulation parameters were set to approximate
as close as possible the experimental conditions while keeping
in mind numerical constraints. The Re, M, Ca, and Bo num-
ber values of our simulations describe qualitatively two-phase
flow regime in clay microfractures. The elastic properties are
scaled so that proportions between them and Pgas,in are the
same as in repository conditions. The values of Pgas,in are
slightly larger than expected in repository conditions in order
to speed up the simulations.

For the particular geometry of our sample containing a
percolating pore, the overall drainage process stayed similar,
though the softening of clay though damaging resulted in
accelerated drainage velocity. Also, the obtained relation for
the suction curve of an isolated pore is decreasing with gas
saturation, which is the opposite from the experimental rela-
tionships obtained for macroscopic samples. One remarkable
feature of capillary pressure curves in Figure 13(b) is that they
are almost identical for saturation of gas larger than 0.1.

The capability of the model to create new pathways due
to localized damage was also investigated. We have demon-
strated that for a small sample under applied conditions it
is difficult to create macroscopic fractures and damage zone
stays diffuse. However, taking into account the rigid inclu-
sions resulted, as expected, in creation of new pores between
them and the elastic matrix.

We conclude that the developed SPH code is capable of
representing all required phenomena for FIB-SEM based geo-
metries in conditions close to the real ones. To be able to
obtain more quantitative information about drainage in argil-
lites some further work needs to be conducted. In particular,
it is not completely clear whether the macroscopic damage
model is appropriate for use if the solid matrix is divided into
clayey matrix and rigid inclusions (for COx they constitute
about 40% of volume). At the nanometric scale also the mass
exchange phenomena can be very important, especially the
possibility for gas to dissolve and for water to evaporate, since
in tortuous pore space the continuous phase flow may be
slowed down or even stopped. Finally, our model should be
applied to larger samples in order to minimize the boundary
effects and to provide data independent of individual pore
space geometry.

Disclosure

A part of this work was presented during JEMP 2016 confer-
ence, held on 12–14 October 2016 in Anglet, France.

**Conflicts of Interest**

The authors declare that they have no conflicts of interest.

Acknowledgments

The financial support from NEEDS MIPOR is gratefully
acknowledged. The help of Alain Genty (CEA) who provided
a GPU cluster for simulations is greatly appreciated.

**References**

[1] S. Cariou, Z. Duan, C. Davy, F. Skoczylas, and L. Dormieux, "Poromechanics of partially saturated COx argillite," *Applied Clay Science*, vol. 56, pp. 36–47, 2012.

[2] ANDRA Dossier 2005 Argile. Synthesis. Evaluation of the feasibility of a geological repository in an argillaceous formation. ANDRA Report Series, 782 p., ISBN 2-916162-00-3, (2005).

[3] R. Cuss, J. Harrington, R. Giot, and C. Auvray, "Experimental observations of mechanical dilation at the onset of gas flow in Callovo-Oxfordian claystone," *Geological Society, London, Special Publications*, vol. 400, no. 1, pp. 507–519, 2014.

[4] P. Marschall, S. Horseman, and T. Gimmi, “Characterisation of gas transport properties of the Opalinus Clay, a potential host rock formation for radioactive waste disposal,” *Oil & Gas Science and Technology - Revue d’IFP Energies nouvelles*, vol. 60, no. 1, pp. 121–139, 2005.

[5] F. Bardelli, C. Mondelli, M. Didier et al., “Hydrogen uptake and diffusion in Callovo-Oxfordian clay rock for nuclear waste disposal technology,” *Applied Geochemistry*, vol. 49, pp. 168–177, 2014.

[6] C. Davy, S. M'Jahad, F. Skoczylas, and J. Talandier, *Migration de gaz au travers de l'argilite de Bure: mise en evidence d'un passage discontinu*, Rencontres AUGC, 2013.

[7] J. F. Harrington, R. De La Vaissière, D. J. Noy, R. J. Cuss, and J. Talandier, “Gas flow in Callovo-Oxfordian claystone (COx): Results from laboratory and field-scale measurements,” *Miner- alogical Magazine*, vol. 76, no. 8, pp. 3303–3318, 2012.

[8] L. M. Keller, P. Schuetz, R. Erni et al., "Characterization of multi-scale microstructural features in opalinus clay," *Microporous and Mesoporous Materials*, vol. 170, pp. 83–94, 2013.

[9] D. Arnedo, E. E. Alonso, and S. Olivella, "Gas flow in anisotropic claystone: Modelling triaxial experiments," *International Journal for Numerical and Analytical Methods in Geomechanics*, vol. 37, no. 14, pp. 2239–2256, 2013.

[10] P. Gerard, J. Harrington, R. Charlier, and F. Collin, "Modelling of localised gas preferential pathways in claystone," *International Journal of Rock Mechanics and Mining Sciences*, vol. 67, pp. 104–114, 2014.

[11] S. Olivella and E. E. Alonso, "Gas flow through clay barriers," *Géotechnique*, vol. 58, no. 3, pp. 157–176, 2008.

[12] A. Y. Rozhko, "Two-phase fluid-flow modeling in a dilatant crack-like pathway," *Journal of Petroleum Science and Engineering*, vol. 146, pp. 1158–1172, 2016.

[13] M. Fall, O. Nasir, and T. S. Nguyen, “A coupled hydro-mechanical model for simulation of gas migration in host sedimentary rocks for nuclear waste repositories,” *Engineering Geology*, vol. 176, no. 5, pp. 24–44, 2014.
[14] C. A. Tang, L. G. Tham, P. K. K. Lee, T. H. Yanng, and L. C. Li, "Coupled analysis of flow, stress and damage (FSD) in rock failure," *International Journal of Rock Mechanics & Mining Sciences*, vol. 39, pp. 477-489, 2002.

[15] S. Y. Wang, S. W. Sloan, S. G. Fityus, D. V. Griffiths, and C. A. Tang, "Numerical modeling of pore pressure influence on fracture evolution in brittle heterogeneous rocks," *Rock Mechanics and Rock Engineering*, vol. 46, no. 5, pp. 1165-1182, 2013.

[16] T. H. Yang, J. Liu, W. C. Zhu, D. Elsworth, L. G. Tham, and C. A. Tang, "A coupled flow-stress-damage model for groundwater outbursts from an underlying aquifer into mining excavations," *International Journal of Rock Mechanics and Mining Sciences*, vol. 44, no. 1, pp. 87-97, 2007.

[17] S. B. Hassine, M. Dymitrowska, V. Pot, and A. Genty, "Gas migration in highly water-saturated Opalinus clay microfractures using a two-phase TRT LBM," *Transport in Porous Media*, vol. 116, no. 3, pp. 975-1003, 2017.

[18] R. A. Gingold and J. J. Monaghan, "Smoothed particle hydrodynamics: theory and application to non-spherical stars," *Monthly Notices of the Royal Astronomical Society*, vol. 181, pp. 375-398, 1997.

[19] L. B. Lucy, "A numerical approach to the testing of the fission hypothesis," *The Astronomical Journal*, vol. 82, pp. 1013-1024, 1977.

[20] A. M. Tartakovskiy and P. Meakin, "Pore scale modeling of immiscible and miscible fluid flows using smoothed particle hydrodynamics," *Advances in Water Resources*, vol. 29, no. 10, pp. 1464-1478, 2006.

[21] U. C. Bandara, A. M. Tartakovskiy, M. Oostrom, B. J. Palmer, J. Grate, and C. Zhang, "Smoothed particle hydrodynamics pore-scale simulations of unstable immiscible flow in porous media," *Advances in Water Resources*, vol. 62, pp. 356-369, 2013.

[22] A. M. Tartakovskiy, A. L. Ward, and P. Meakin, "Pore-scale simulations of drainage of heterogeneous and anisotropic porous media," *Physics of Fluids*, vol. 19, no. 10, Article ID 103301, 2007.

[23] S. Ovaysi and M. Piri, "Direct pore-level modeling of incompressible fluid flow in porous media," *Journal of Computational Physics*, vol. 229, no. 19, pp. 7456-7476, 2010.

[24] P. Kunz, I. M. Zarikos, N. K. Karadimitriou, M. Huber, U. Nieken, and S. M. Hassanizadeh, "Study of Multi-phase Flow in Porous Media: Comparison of SPH Simulations with Micro-model Experiments," *Transport in Porous Media*, vol. 114, no. 2, pp. 581-600, 2016.

[25] L. Lobovsky and J. Kren, *Applied and Computational Mechanics*, vol. 1, 2007.

[26] J. J. Monaghan, "SPH without a tensile instability," *Journal of Computational Physics*, vol. 159, no. 2, pp. 290-311, 2000.

[27] A. Amicarelli, R. Albano, D. Mirauda, G. Agate, A. Sole, and R. Guandalini, "A smoothed particle hydrodynamics model for 3D solid body transport in free surface flows," *Computers & Fluids. An International Journal*, vol. 116, pp. 205-228, 2015.

[28] Y. Amini, H. Emdad, and M. Farid, "A new model to solve fluid-hypo-elastic solid interaction using the smoothed particle hydrodynamics (SPH) method," *European Journal of Mechanics - B/Fluids*, vol. 30, no. 2, pp. 184-194, 2011.

[29] C. Antoci, M. Gallati, and S. Sibilla, "Numerical simulation of fluid-structure interaction by SPH," *Computers & Structures*, vol. 85, no. 11-14, pp. 879-890, 2007.

[30] M. Lahooti, A. Pisheva, and M. S. Saidi, "Novel 2D algorithm for fluid solid interaction based on the smoothed particle hydrodynamics (SPH) method," *Scientia Iranica*, vol. 18, no. 3 B, pp. 358-367, 2011.

[31] S. Adami, X. Y. Hu, and N. A. Adams, "A conservative SPH method for surfactant dynamics," *Journal of Computational Physics*, vol. 229, no. 5, pp. 1909-1926, 2010.

[32] S. Adami, X. Y. Hu, and N. A. Adams, "A new surface-tension formulation for multi-phase SPH using a reproducing divergence approximation," *Journal of Computational Physics*, vol. 229, no. 13, pp. 5011-5021, 2010.

[33] T. Breinlinger, P. Polfer, A. Hashibon, and T. Kraft, "Surface tension and wetting effects with smoothed particle hydrodynamics," *Journal of Computational Physics*, vol. 243, pp. 14-27, 2013.

[34] M. B. Liu and G. R. Liu, "Meshfree particle simulation of micro channel flows with surface tension," *Computational Mechanics*, vol. 35, no. 5, pp. 332-341, 2005.

[35] J. P. Morris, "Simulating surface tension with smoothed particle hydrodynamics," *International Journal for Numerical Methods in Fluids*, vol. 33, no. 3, pp. 333-353, 2000.

[36] J. U. Brackbill, D. B. Kothe, and C. Zemach, "A continuum method for modeling surface tension," *Journal of Computational Physics*, vol. 100, no. 2, pp. 335-354, 1992.

[37] A. G. Corkum and C. D. Martin, "The mechanical behaviour of weak mudstone (Opalinus Clay) at low stresses," *International Journal of Rock Mechanics and Mining Sciences*, vol. 44, no. 2, pp. 196-209, 2007.

[38] J. Mazars, "A description of micro- and macroscale damage of concrete structures," *Engineering Fracture Mechanics*, vol. 25, no. 5-6, pp. 729-737, 1986.

[39] P. W. Randles and L. D. Libersky, "Smoothed particle hydrodynamics: some recent improvements and applications," *Computer Methods Applied Mechanics and Engineering*, vol. 139, no. 1-4, pp. 375-408, 1996.

[40] M. B. Liu and G. R. Liu, "Smoothed particle hydrodynamics (SPH): an overview and recent developments," *Archives of Computational Methods in Engineering: State-of-the-Art Reviews*, vol. 17, no. 1, pp. 25-76, 2010.

[41] J. J. Monaghan, "Smoothed particle hydrodynamics," *Reports on Progress in Physics*, vol. 68, no. 8, pp. 1703-1759, 2005.

[42] M. Muller, D. Charypar, and M. Gross, "Particle-Based Fluid Simulation for Interactive Applications," *Euragraphics/SIGGRAPH Symp. on Comp. Anim*, 2003.

[43] M. Desbrun and M. P. Gascuel, "Smoothed particles: A new paradigm for animating highly deformable bodies. Computer Animation and Simulation," in *Proceedings of the EG Workshop on Animation and Simulation*, pp. 61-76, Springer-Verlag, 1996.

[44] A. M. Tartakovskiy, N. Trask, K. Pan, B. Jones, W. Pan, and J. R. Williams, "Smoothed particle hydrodynamics and its applications for multiphase flow and reactive transport in porous media," *Computational Geosciences*, vol. 20, no. 4, pp. 807-834, 2016.

[45] X. Y. Hu and N. A. Adams, "An incompressible multi-phase SPH method," *Journal of Computational Physics*, vol. 227, no. 1, pp. 264-278, 2007.

[46] M. Muller, R. Keiser, A. Nealen, M. Pauly, M. Gross, and M. Alexa, "Point based animation of elastic, plastic and melting objects," in *Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation (SCA '04)*, pp. 141-151, Eurographics Association, Aire-la-Ville, Switzerland, 2004.

[47] B. Solenthaler, J. Schläfli, and R. Pajarola, "A unified particle model for fluid-solid interactions," *Computer Animation and Virtual Worlds*, vol. 18, no. 1, pp. 69-82, 2007.

[48] A. Pazdniakou and P. M. Adler, "Lattice spring models," *Transport in Porous Media*, vol. 93, no. 2, pp. 243-262, 2012.
[49] N. Akinci, M. Ihmsen, G. Akinci, B. Solenthaler, and M. Teschner, "Versatile rigid-fluid coupling for incompressible SPH," *ACM Transactions on Graphics*, vol. 31, no. 4, article 62, 2012.

[50] X. Shao, Z. Zhou, N. Magnenat-Thalmann, and W. Wu, "Stable and fast fluid-solid coupling for incompressible SPH," *Computer Graphics Forum*, vol. 34, no. 1, pp. 191–204, 2015.

[51] M. Müller, S. Schirm, M. Teschner, B. Heidelberger, and M. Gross, "Interaction of fluids with deformable solids," *Computer Animation and Virtual Worlds*, vol. 15, no. 3-4, pp. 159–171, 2004.

[52] J. J. Monaghan, "Simulating Free Surface Flows with SPH," *Journal of Computational Physics*, vol. 110, no. 2, pp. 399–406, 1994.

[53] D. W. Holmes, J. R. Williams, and P. Tilke, "Smooth particle hydrodynamics simulations of low Reynolds number flows through porous media," *International Journal for Numerical and Analytical Methods in Geomechanics*, vol. 35, no. 4, pp. 419–437, 2011.

[54] R. E. Larson and J. J. L. Higdon, "A periodic grain consolidation model of porous media," *Physics of Fluids*, vol. 1, no. 1, pp. 38–46, 1989.

[55] A. S. Sangani and A. Acrivos, "Slow flow through a periodic array of spheres," *International Journal of Multiphase Flow*, vol. 8, no. 4, pp. 343–360, 1982.

[56] Y. Song, C. A. Davy, D. Troadec et al., "Multi-scale pore structure of COx claystone: Towards the prediction of fluid transport," *Marine and Petroleum Geology*, vol. 65, pp. 63–82, 2015.

[57] OpenMP Application Programming Interface. Version 4.5 November 2015. http://www.openmp.org/wp-content/uploads/openmp-4.5.pdf.

[58] D. Yang, *Caractérisation par la mesure de perméabilité au gaz de l'endommagement mécanique et hydrique dans l'EDZ des argilites du Callovo-Oxfordien [Ph.D. thesis]*, 2008.
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /></figure>

Hindawi
Submit your manuscripts at
www.hindawi.com