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

Authors: Aliaksei Pazdniakou; Magdalena Dymitrowska
Journal: Geofluids
 (2018)

---

## Abstract

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


---

## Body

## 1. Introduction

The principal objective of this study is to improve the representation and understanding of gas migration in highly water saturated clays. This objective is closely related to several projects that use a deep sedimentary rock formation, such as Callovo-Oxfordian (COx) clay [1], as a host rock and a natural barrier for radioactive waste repository [2]. These sedimentary formations are characterized by a low permeability and are water saturated in their initial state. However, during the postclosure phase, the gases, mainly hydrogen, are being produced mostly by anaerobic corrosion of iron components of the repository. This important gas production may lead to appearance of a free gas phase and to a local increase of pressure, which could influence confining properties of the repository because of desaturation of clay pores and because of the mechanical disturbance caused by a possible microfracturing in the vicinity of the installation.In classical approaches, several gas transport mechanisms can be active, separately or simultaneously [3, 4]. With increasing rates of gas production, gas should be transported first in dissolved state by advection and diffusion [5] and then by the flow of small gas bubbles inside the water phase and by viscocapillary flow of gas inside the water phase, or by the dilatation of existing and creation of new localized pathways [3, 6, 7]. The ultimate stage appearing at very high gas production rates would be macroscopic fracturing of the host rock; however, it is rather not expected in the case of repositories, where the gas source terms are not sufficient. Experimental evidences accumulating in recent years indicate that in case of indurated clays the gas migration may directly switch from dissolution-diffusion stage to the dilatant pathways mechanism. Indeed, when an increasing gas pressure is applied to a saturated argillite sample, at a certain gas pressure (“entry pressure”), a rapid but intermittent increase of the gas outflow is usually observed [3, 6, 7]. This gas flow is not associated with a significant sample desaturation, as it would be the case in viscocapillary flows. On the other hand, the gas flow causes a macroscopic deformation of the clay sample and the associated permeability tends to increase when repeating gas pulse tests. The prevalence of such a dilatant flow can be justified by the fact that the major part of the porosity of indurated clays exhibits pore size of a few nanometers [4]. Thus, due to associated high capillary pressure in such water saturated nanosize pores, the viscocapillary and bubbly flows would be restricted to largest pores which do not form a connected macroscopic paths and are not suffecient for the macroscopic gas transfer. On the other hand, the low tensile strength of the claystone facilitates the solid matrix cracking and subsequently the gas percolation [8]. Therefore, in order to get a better understanding of dilatant transport mode, it becomes necessary to take into account the coupling of the two phase fluid/gas flow inside the pore space with solid matrix mechanics, which constitute the scope of our study together with modelling of modification of argillite mechanical properties due to damage induced by gas high pressures. The phenomena of gas dissolution and diffusion are not considered in our model, since they are out of scope of this paper. Gas slippage is also excluded from considerations, since at high gas pressures these effects are of limited importance.The coupling between the fluid flow hydrodynamics and solid matrix deformation has been tackled recently by several authors. In [9–11], a coupled 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 acquire3D 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 micromodel experiments was reported in [24]. It was also demonstrated [25, 26] that SPH is capable of simulating elastic media and that solid-fluid interactions can be relatively easily coupled in the framework of the same approach [27–30]. An interesting feature of the SPH method is its almost local character, which makes it an ideal choice for massively parallel computing, for example, with Compute Unified Device Architecture (CUDA) on Graphics Processing Units (GPU). For all these reasons and combined with the fact that different physical phenomena can be relatively easily expressed in terms of the SPH model, this method can be considered as an attractive choice for our application.The article is organized as follows. In Section2, 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(1)DρDt=-ρ∇·v,where ρ is the density, v is the velocity, t is the time, and D denotes the material derivative.Momentum conservation equation in Lagrangian frame can be written as(2)ρDvDt=-∇p+∇·μ∇v+∇vT-23∇·vI+ρg,where p is the pressure, g is the gravitational acceleration, and μ is the dynamic viscosity. The momentum equation for incompressible Newtonian fluid in Lagrangian frame can be stated as(3)ρDvDt=-∇p+μ∇2v+ρg.In order to complete the description of fluid state at one point in terms of density ρ, pressure p, and velocity 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(4)p=cs2ρ-ρ0,where cs is the speed of sound in the fluid and ρ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 Δp for a curved interface between two fluids, which can be written in terms of the Laplace law(5)Δp=σ1R1+1R2=σ∇·n^,where σ is the surface tension coefficient, R1 and R2 are the principal radii of curvature, and 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 forceFs(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(6)cr=ρlr-ρgrρlr+ρgr,where subscripts l and g denote liquid and gas, respectively. The values of the color field c(r) can be used to calculate the local curvature κ of the interface as follows:(7)κ=-∇·n^=n·∇nn2-∇·nn,where the unit normal n^ is given by(8)n^=∇cr∇cr.The volumetric surface tension force Fs(r) is expressed as(9)Fsr=12σκn.To complete the description of the interface phenomena, the contact angle θ can be introduced by assigning the color function value to the solid phase so that(10)cr=cos⁡θ.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(11)ε=12∇u+∇uT,where u is the local displacement. The Hooke law linearly relates the strain tensor ε to the stress tensor σ as(12)σ=C:ε,where C is the stiffness tensor. For an isotropic medium, the components of the stiffness tensor are given by(13)Cijkl=λsδijδkl+μsδikδjl+δilδjk,where λs and μs are the Lamé coefficients and δij is the Kronecker delta function. Then, (12) can be simplified as(14)σ=2μsε+λs∇·uI=μs∇u+∇uT+λs∇·uI,where I is the unit tensor. This expression can be also written as(15)σ=μs∇u+∇uT-23∇·uI+Ks∇·uI,where Ks is the bulk modulus and Ks(∇·u) can be associated with pressure inside the elastic medium. The dynamic behaviour of an elastic medium is described by(16)ρs∂2u∂t2=∇·σ.The volumetric elastic force density components ∇·σ=(Fx,Fy,Fz)=Felastic of an isotropic medium can be expressed in terms of partial derivatives of the local displacement field as(17)Fx=2μs+λs∂2ux∂x∂x+μs∂2ux∂y∂y+∂2ux∂z∂z+μs+λs∂2uy∂x∂y+∂2uz∂x∂zFy=2μs+λs∂2uy∂y∂y+μs∂2uy∂x∂x+∂2uy∂z∂z+μs+λs∂2ux∂x∂y+∂2uz∂z∂yFz=2μs+λs∂2uz∂z∂z+μs∂2uz∂x∂x+∂2uz∂y∂y+μs+λs∂2ux∂x∂z+∂2uy∂y∂z.


#### 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 variableD [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]:(18)E=1-DE0,where E0 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 Dt and Dc 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]:(19)σ3≤-ft,where ft is the tensile strength and σ3 is the smallest principal stress. If the criterion is satisfied, the damage variable Dt associated with it is defined as [14](20)Dt=0ifεt0≤ε1-ftrE0εifεtu≤ε≤εt01ifε≤εtu,where ftr is the residual tensile strength, εt0 is the strain at the elastic limit, and ε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](21)σ1-σ31+sin⁡ϕ1-sin⁡ϕ≥fc,where ϕ is the internal friction angle, fc is the uniaxial compressive strength, and σ1 is the largest principal stress. The damage variable Dc associated with this criterion is given by [14](22)Dc=0ifε≤εc01-fcrE0εifεc0≤ε,where fcr is the residual compressive strength and εc0 is the compressive strain at the elastic limit. Following [13, 14], ε in (20) and (22) in 3D case can be replaced with principal strains ε3 and ε1, respectively. A schematic shape of the constitutive damage law given by (20) and (22) is presented in Figure 1.Figure 1
An example of constitutive damage law.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 ft is significantly smaller than their uniaxial compressive strength fc 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.Table 1
COx clay properties.
Parameter Value Young’s modulusE 3 GPa Poisson’s ratioν 0.1 Tensile strengthft 4 MPa Uniaxial compressive strengthfc 20 MPa Internal friction angle 25∘ Strain at elastic limitεt0 0.0013 Ultimate tensile strainεtu 0.0026 Compressive strain at elastic limitεc0 0.015The total damageD is a combination of Dt and Dc and can be calculated as [13](23)D=ζtDt+ζcDc,where the coefficients ζt and ζ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 valueD(t+Δt) based on D value calculated by (23)(24)Dt+Δt=Dt+D-Dtτ,where τ 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≈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,(25a)dvCMdt=FBmB(25b)drCMdt=vCM(25c)dωBdt=IB-1TB-ωB×IBωB(25d)dϕBdt=ωB,where rCM is the position of the center of mass and vCM its velocity, FB is the total external force acting on the body, mB is the mass of the body, ωB is the angular velocity, IB is the tensor of the moment of inertia, TB is the total torque acting on the body, and ϕ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(26)rCM=1mB∫ρrrd3r.The tensor of the moment of inertia is given by(27)IB=∫ρrr-rCM2I-r-rCMr-rCMTd3r.The total torque acting on the body can be calculated as(28)TB=∫r-rCM×Fd3r,where 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 interaction between phases at the solid-fluid interface, the boundary conditions must be defined. In most cases, the velocity of both phases at the contact pointri is supposed to be the same:(29)vfri=vsri.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(30)σs·ns=-σf·nf,where σf and σs are fluid and solid phase stress tensors, respectively, nf is the unit normal pointing out of the fluid phase, and ns=-nf.At the interface between the solid phases, the continuity of normal stresses is imposed:(31)σs·ns=-σe·ne.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.


## 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(1)DρDt=-ρ∇·v,where ρ is the density, v is the velocity, t is the time, and D denotes the material derivative.Momentum conservation equation in Lagrangian frame can be written as(2)ρDvDt=-∇p+∇·μ∇v+∇vT-23∇·vI+ρg,where p is the pressure, g is the gravitational acceleration, and μ is the dynamic viscosity. The momentum equation for incompressible Newtonian fluid in Lagrangian frame can be stated as(3)ρDvDt=-∇p+μ∇2v+ρg.In order to complete the description of fluid state at one point in terms of density ρ, pressure p, and velocity 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(4)p=cs2ρ-ρ0,where cs is the speed of sound in the fluid and ρ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 Δp for a curved interface between two fluids, which can be written in terms of the Laplace law(5)Δp=σ1R1+1R2=σ∇·n^,where σ is the surface tension coefficient, R1 and R2 are the principal radii of curvature, and 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 forceFs(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(6)cr=ρlr-ρgrρlr+ρgr,where subscripts l and g denote liquid and gas, respectively. The values of the color field c(r) can be used to calculate the local curvature κ of the interface as follows:(7)κ=-∇·n^=n·∇nn2-∇·nn,where the unit normal n^ is given by(8)n^=∇cr∇cr.The volumetric surface tension force Fs(r) is expressed as(9)Fsr=12σκn.To complete the description of the interface phenomena, the contact angle θ can be introduced by assigning the color function value to the solid phase so that(10)cr=cos⁡θ.These values are then used in calculations of the unit normal near the solid surface.


## 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 Δp for a curved interface between two fluids, which can be written in terms of the Laplace law(5)Δp=σ1R1+1R2=σ∇·n^,where σ is the surface tension coefficient, R1 and R2 are the principal radii of curvature, and 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 forceFs(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(6)cr=ρlr-ρgrρlr+ρgr,where subscripts l and g denote liquid and gas, respectively. The values of the color field c(r) can be used to calculate the local curvature κ of the interface as follows:(7)κ=-∇·n^=n·∇nn2-∇·nn,where the unit normal n^ is given by(8)n^=∇cr∇cr.The volumetric surface tension force Fs(r) is expressed as(9)Fsr=12σκn.To complete the description of the interface phenomena, the contact angle θ can be introduced by assigning the color function value to the solid phase so that(10)cr=cos⁡θ.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(11)ε=12∇u+∇uT,where u is the local displacement. The Hooke law linearly relates the strain tensor ε to the stress tensor σ as(12)σ=C:ε,where C is the stiffness tensor. For an isotropic medium, the components of the stiffness tensor are given by(13)Cijkl=λsδijδkl+μsδikδjl+δilδjk,where λs and μs are the Lamé coefficients and δij is the Kronecker delta function. Then, (12) can be simplified as(14)σ=2μsε+λs∇·uI=μs∇u+∇uT+λs∇·uI,where I is the unit tensor. This expression can be also written as(15)σ=μs∇u+∇uT-23∇·uI+Ks∇·uI,where Ks is the bulk modulus and Ks(∇·u) can be associated with pressure inside the elastic medium. The dynamic behaviour of an elastic medium is described by(16)ρs∂2u∂t2=∇·σ.The volumetric elastic force density components ∇·σ=(Fx,Fy,Fz)=Felastic of an isotropic medium can be expressed in terms of partial derivatives of the local displacement field as(17)Fx=2μs+λs∂2ux∂x∂x+μs∂2ux∂y∂y+∂2ux∂z∂z+μs+λs∂2uy∂x∂y+∂2uz∂x∂zFy=2μs+λs∂2uy∂y∂y+μs∂2uy∂x∂x+∂2uy∂z∂z+μs+λs∂2ux∂x∂y+∂2uz∂z∂yFz=2μs+λs∂2uz∂z∂z+μs∂2uz∂x∂x+∂2uz∂y∂y+μs+λs∂2ux∂x∂z+∂2uy∂y∂z.


### 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 variableD [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]:(18)E=1-DE0,where E0 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 Dt and Dc 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]:(19)σ3≤-ft,where ft is the tensile strength and σ3 is the smallest principal stress. If the criterion is satisfied, the damage variable Dt associated with it is defined as [14](20)Dt=0ifεt0≤ε1-ftrE0εifεtu≤ε≤εt01ifε≤εtu,where ftr is the residual tensile strength, εt0 is the strain at the elastic limit, and ε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](21)σ1-σ31+sin⁡ϕ1-sin⁡ϕ≥fc,where ϕ is the internal friction angle, fc is the uniaxial compressive strength, and σ1 is the largest principal stress. The damage variable Dc associated with this criterion is given by [14](22)Dc=0ifε≤εc01-fcrE0εifεc0≤ε,where fcr is the residual compressive strength and εc0 is the compressive strain at the elastic limit. Following [13, 14], ε in (20) and (22) in 3D case can be replaced with principal strains ε3 and ε1, respectively. A schematic shape of the constitutive damage law given by (20) and (22) is presented in Figure 1.Figure 1
An example of constitutive damage law.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 ft is significantly smaller than their uniaxial compressive strength fc 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.Table 1
COx clay properties.
Parameter Value Young’s modulusE 3 GPa Poisson’s ratioν 0.1 Tensile strengthft 4 MPa Uniaxial compressive strengthfc 20 MPa Internal friction angle 25∘ Strain at elastic limitεt0 0.0013 Ultimate tensile strainεtu 0.0026 Compressive strain at elastic limitεc0 0.015The total damageD is a combination of Dt and Dc and can be calculated as [13](23)D=ζtDt+ζcDc,where the coefficients ζt and ζ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 valueD(t+Δt) based on D value calculated by (23)(24)Dt+Δt=Dt+D-Dtτ,where τ 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≈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,(25a)dvCMdt=FBmB(25b)drCMdt=vCM(25c)dωBdt=IB-1TB-ωB×IBωB(25d)dϕBdt=ωB,where rCM is the position of the center of mass and vCM its velocity, FB is the total external force acting on the body, mB is the mass of the body, ωB is the angular velocity, IB is the tensor of the moment of inertia, TB is the total torque acting on the body, and ϕ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(26)rCM=1mB∫ρrrd3r.The tensor of the moment of inertia is given by(27)IB=∫ρrr-rCM2I-r-rCMr-rCMTd3r.The total torque acting on the body can be calculated as(28)TB=∫r-rCM×Fd3r,where 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.2.1. Elastic Material

The deformations inside an elastic material are given by the strain or deformation tensor(11)ε=12∇u+∇uT,where u is the local displacement. The Hooke law linearly relates the strain tensor ε to the stress tensor σ as(12)σ=C:ε,where C is the stiffness tensor. For an isotropic medium, the components of the stiffness tensor are given by(13)Cijkl=λsδijδkl+μsδikδjl+δilδjk,where λs and μs are the Lamé coefficients and δij is the Kronecker delta function. Then, (12) can be simplified as(14)σ=2μsε+λs∇·uI=μs∇u+∇uT+λs∇·uI,where I is the unit tensor. This expression can be also written as(15)σ=μs∇u+∇uT-23∇·uI+Ks∇·uI,where Ks is the bulk modulus and Ks(∇·u) can be associated with pressure inside the elastic medium. The dynamic behaviour of an elastic medium is described by(16)ρs∂2u∂t2=∇·σ.The volumetric elastic force density components ∇·σ=(Fx,Fy,Fz)=Felastic of an isotropic medium can be expressed in terms of partial derivatives of the local displacement field as(17)Fx=2μs+λs∂2ux∂x∂x+μs∂2ux∂y∂y+∂2ux∂z∂z+μs+λs∂2uy∂x∂y+∂2uz∂x∂zFy=2μs+λs∂2uy∂y∂y+μs∂2uy∂x∂x+∂2uy∂z∂z+μs+λs∂2ux∂x∂y+∂2uz∂z∂yFz=2μs+λs∂2uz∂z∂z+μs∂2uz∂x∂x+∂2uz∂y∂y+μs+λs∂2ux∂x∂z+∂2uy∂y∂z.


## 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 variableD [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]:(18)E=1-DE0,where E0 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 Dt and Dc 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]:(19)σ3≤-ft,where ft is the tensile strength and σ3 is the smallest principal stress. If the criterion is satisfied, the damage variable Dt associated with it is defined as [14](20)Dt=0ifεt0≤ε1-ftrE0εifεtu≤ε≤εt01ifε≤εtu,where ftr is the residual tensile strength, εt0 is the strain at the elastic limit, and ε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](21)σ1-σ31+sin⁡ϕ1-sin⁡ϕ≥fc,where ϕ is the internal friction angle, fc is the uniaxial compressive strength, and σ1 is the largest principal stress. The damage variable Dc associated with this criterion is given by [14](22)Dc=0ifε≤εc01-fcrE0εifεc0≤ε,where fcr is the residual compressive strength and εc0 is the compressive strain at the elastic limit. Following [13, 14], ε in (20) and (22) in 3D case can be replaced with principal strains ε3 and ε1, respectively. A schematic shape of the constitutive damage law given by (20) and (22) is presented in Figure 1.Figure 1
An example of constitutive damage law.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 ft is significantly smaller than their uniaxial compressive strength fc 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.Table 1
COx clay properties.
Parameter Value Young’s modulusE 3 GPa Poisson’s ratioν 0.1 Tensile strengthft 4 MPa Uniaxial compressive strengthfc 20 MPa Internal friction angle 25∘ Strain at elastic limitεt0 0.0013 Ultimate tensile strainεtu 0.0026 Compressive strain at elastic limitεc0 0.015The total damageD is a combination of Dt and Dc and can be calculated as [13](23)D=ζtDt+ζcDc,where the coefficients ζt and ζ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 valueD(t+Δt) based on D value calculated by (23)(24)Dt+Δt=Dt+D-Dtτ,where τ 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≈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,(25a)dvCMdt=FBmB(25b)drCMdt=vCM(25c)dωBdt=IB-1TB-ωB×IBωB(25d)dϕBdt=ωB,where rCM is the position of the center of mass and vCM its velocity, FB is the total external force acting on the body, mB is the mass of the body, ωB is the angular velocity, IB is the tensor of the moment of inertia, TB is the total torque acting on the body, and ϕ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(26)rCM=1mB∫ρrrd3r.The tensor of the moment of inertia is given by(27)IB=∫ρrr-rCM2I-r-rCMr-rCMTd3r.The total torque acting on the body can be calculated as(28)TB=∫r-rCM×Fd3r,where 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 interaction between phases at the solid-fluid interface, the boundary conditions must be defined. In most cases, the velocity of both phases at the contact pointri is supposed to be the same:(29)vfri=vsri.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(30)σs·ns=-σf·nf,where σf and σs are fluid and solid phase stress tensors, respectively, nf is the unit normal pointing out of the fluid phase, and ns=-nf.At the interface between the solid phases, the continuity of normal stresses is imposed:(31)σs·ns=-σe·ne.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 functionf, its value can be found as [40](32)fx=∫fx′δx-x′dx′,where δ(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](33)fx=∫fx′Wx-x′,hdx′+errsmoothing,where W(x-x′,h) is a smoothing function, h is a smoothing length, and errsmoothing is the approximation error. The integral in (33) represents the smoothed value of the function f [40](34)fSx=∫fx′Wx-x′,hdx′.The approximation error errsmoothing can be estimated as errsmoothing=O(h2) 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](35)∇·fxS=∫fx′·∇Wxx-x′,hdx′.The higher order derivatives can be obtained in the similar way.


#### 3.1.1. Smoothing Function

The smoothing functionW(x-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(x-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](36)Wx-x′,h=16πh312-3r2+3r3if0≤r<0.51-r3if0.5≤r<10ifr≥1,where r=x-x′/h. It is clear that W(x-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 massmi and density ρi associated with it. The particle volume can be defined as(37)Vi=miρi.The continuous integrals (34) and (35) can be approximated by their discrete analogues,(38a)fi=∑jfjWxi-xj,hVj=∑jfjmjρjWxi-xj,h(38b)∇fi=∑jfj∇Wxi-xj,hVj=∑jfjmjρj∇Wxixi-xj,h,where the summation in function term (38a) is carried out for all neighbouring particles within the distance h from the ith particle including j=i; note that in the derivative term j≠i (38b). According to [44], the error introduced by the finite sum approximation of continuous integrals can be estimated as errintegral=O(Δ/h), where Δ is the length scale characterizing the average distance between the particles.Whenf is a density, (38a) is transformed to(39)ρi=∑jmjWxi-xj,h.Usually, the masses of SPH particles representing the same material are equal. Another useful quantity is the number density, which can be calculated as(40)ni=ρimi=∑jWxi-xj,hand 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](41)ρi=∑jmjWxi-xj,h∑jWxi-xj,h/nj.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]:(42)∇f=1F∇Ff-f∇F,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=ρ, and F=1/ρ, substitution of which in (42) yields the following expressions of SPH derivative term:(43a)∇fi=∑jfj-fimjρj∇Wxixi-xj,h(43b)∇fi=1ρi∑jfj-fimj∇Wxixi-xj,h(43c)∇fi=ρi∑jfjρj2+fiρi2mj∇Wxixi-xj,h.The calculation of a Laplacian of some functions is often required for the transport equation. It can be obtained as(44)∇2fi=∑jfjmjρj∇2Wxixi-xj,h.This formula also can be symmetrized just as SPH derivatives terms; however, even symmetrized formulations may suffer from ∇2Wxi(xi-xj,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]:(45)∇·F∇fi=∑jFj+Fifi-fjmjρjxi-xjxi-xj2∇Wxixi-xj,h,where f and F are some scalar functions. It is useful to note that [44](46)∇Wxixi-xj,h=xi-xjxi-xj∂Wrij,h∂rij,where rij=xi-xj.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:(47a)dxidt=vi(47b)dvidt=ai,where vi is the particle velocity and ai 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](48a)vit+Δt2=vit-Δt2+Δtait(48b)xit+Δt=xit+Δtvit+Δt2.It is similar to the velocity Verlet algorithm, which can be expressed as [44](49a)xit+Δt=xit+Δtvit+Δt2ait2(49b)vit+Δt=vit+Δt2ait+ait+Δt.Both methods are second-order integration schemes. For numerical stability the time step Δt must obey the Courant-Friedrich-Levy condition [35](50)Δt≤0.25hvmaxand particle acceleration related condition(51)Δt≤0.25hamax.When viscous dissipation term is present, the following condition must be satisfied:(52)Δt≤0.125ρh2μ.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](53a)DρiDt=ρi∑jmjρjvi-vj·∇Wxixi-xj,h(53b)DρiDt=∑jmjvi-vj·∇Wxixi-xj,h.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 h2) [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:(54a)∇pi=∑jmjρjpi+pj∇Wxixi-xj,h(54b)∇pi=ρi∑jpjρj2+piρi2mj∇Wxixi-xj,h.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:(55)μ∇2vi=2μ∑jmjρjvi-vjxi-xj∂Wrij,h∂rij.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(56)μ∇2vi=∑j4μiμjμi+μjmjρjvi-vjxi-xj∂Wrij,h∂rij.Thus, the discretized momentum equation can be written as follows:(57)DviDt=-∑jpjρj2+piρi2mj∇Wxixi-xj,h+∑j4μiμjμi+μjmjρjρivi-vjxi-xj∂Wrij,h∂rij+g+Fsρi,where Fs 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 soundcs 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 (Section2.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 cf=1 and cg=-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(58)cri=∑jmjρjcjWxi-xj,h.Also, the color function can be calculated based on (6). The color function normal is calculated as(59a)nri=∑jmjρjcrj-cri∇Wxixi-xj,h.The interface local curvature κ given by (7) is calculated as(60)κri=1nri2∑jmjρjnrjnri-nrinrj·∇Wxixi-xj,h.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 Fs(ri) 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]:(61)∂2F∂xα∂xβi=∑jmjρj5xα,i-xα,jxβ,i-xβ,jxi-xj-δαβFi-Fj1xi-xj∂Wrij,h∂rij,where the Greek indexes α and β indicate spatial dimensions and δαβ 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 θ 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 ofN points of mass mp can be found from(62)rCM=1MB∑i=1Nmpri.The total external force acting on the rigid body is given by(63)FB=∑i=1Nmpai.The tensor of the moment of inertia is obtained from(64)IB=∑i=1Nmpri-rCM2I-ri-rCMri-rCMT.The total torque acting on the body can be calculated as(65)TB=∑i=1Nri-rCM×mpai,where ai 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 aCM and the angular acceleration of the rigid body αB are calculated as [27](66a)aCMt=FBtMB(66b)αBt=IBt-1TBt-ωBt-Δt2×IBtωBt-Δt2.Then, the velocity of the center of mass vCM and the angular velocity of the rigid body ωB are updated [27](67a)vCMt+Δt2=vCMt-Δt2+aCMtΔt(67b)ωBt+Δt2=ωBt-Δt2+αBtΔt.Finally, the position of the center of mass is updated by [27](68)rCMt+Δt=rCMt+vCMt+Δt2Δt.The velocities and the positions of the rigid body points can be obtained as follows [27](69a)vit+Δt2=vCMt+Δt2+ωBt+Δt2×rit-rCMt(69b)rit+Δt=rCMt+Δt+RBt+Δt2rit-rCMt,where RB is the rotation matrix. It can be calculated as [27](70)RBt+Δt2=RxωB,xt+Δt2ΔtRyωB,yt+Δt2ΔtRzωB,zt+Δt2Δt,where Rx(θ), Ry(θ), and Rz(θ) are the matrices which yield the rotation about principal coordinate axes(71)Rxθ=1000cos⁡θ-sin⁡θ0sin⁡θcos⁡θ(72)Ryθ=cos⁡θ0sin⁡θ010-sin⁡θ0cos⁡θ(73)Rzθ=cos⁡θ-sin⁡θ0sin⁡θcos⁡θ0001.


### 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 lengthh. 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:(74)csolid=Ksρ0,solid,where Ks is the bulk modulus of the solid phase and ρ0,solid the initial solid phase density. Ks for the elastic phase is defined by its elastic properties, and Ks 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]:(75)fi=-mjmi+mjBrini,where fi is the force acting on particle i, B is the symmetric function which defines the interaction strength, and ni 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 v^i for update of the SPH particles positions in the adopted numerical integration scheme instead of the usual particle velocity vi which is used in the momentum equation. The corrected velocity is calculated by [41](76)v^i=vi+ϵ∑jmj2vj-viρi+ρj,where ϵ 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].


## 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 functionf, its value can be found as [40](32)fx=∫fx′δx-x′dx′,where δ(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](33)fx=∫fx′Wx-x′,hdx′+errsmoothing,where W(x-x′,h) is a smoothing function, h is a smoothing length, and errsmoothing is the approximation error. The integral in (33) represents the smoothed value of the function f [40](34)fSx=∫fx′Wx-x′,hdx′.The approximation error errsmoothing can be estimated as errsmoothing=O(h2) 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](35)∇·fxS=∫fx′·∇Wxx-x′,hdx′.The higher order derivatives can be obtained in the similar way.


### 3.1.1. Smoothing Function

The smoothing functionW(x-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(x-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](36)Wx-x′,h=16πh312-3r2+3r3if0≤r<0.51-r3if0.5≤r<10ifr≥1,where r=x-x′/h. It is clear that W(x-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 massmi and density ρi associated with it. The particle volume can be defined as(37)Vi=miρi.The continuous integrals (34) and (35) can be approximated by their discrete analogues,(38a)fi=∑jfjWxi-xj,hVj=∑jfjmjρjWxi-xj,h(38b)∇fi=∑jfj∇Wxi-xj,hVj=∑jfjmjρj∇Wxixi-xj,h,where the summation in function term (38a) is carried out for all neighbouring particles within the distance h from the ith particle including j=i; note that in the derivative term j≠i (38b). According to [44], the error introduced by the finite sum approximation of continuous integrals can be estimated as errintegral=O(Δ/h), where Δ is the length scale characterizing the average distance between the particles.Whenf is a density, (38a) is transformed to(39)ρi=∑jmjWxi-xj,h.Usually, the masses of SPH particles representing the same material are equal. Another useful quantity is the number density, which can be calculated as(40)ni=ρimi=∑jWxi-xj,hand 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](41)ρi=∑jmjWxi-xj,h∑jWxi-xj,h/nj.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]:(42)∇f=1F∇Ff-f∇F,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=ρ, and F=1/ρ, substitution of which in (42) yields the following expressions of SPH derivative term:(43a)∇fi=∑jfj-fimjρj∇Wxixi-xj,h(43b)∇fi=1ρi∑jfj-fimj∇Wxixi-xj,h(43c)∇fi=ρi∑jfjρj2+fiρi2mj∇Wxixi-xj,h.The calculation of a Laplacian of some functions is often required for the transport equation. It can be obtained as(44)∇2fi=∑jfjmjρj∇2Wxixi-xj,h.This formula also can be symmetrized just as SPH derivatives terms; however, even symmetrized formulations may suffer from ∇2Wxi(xi-xj,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]:(45)∇·F∇fi=∑jFj+Fifi-fjmjρjxi-xjxi-xj2∇Wxixi-xj,h,where f and F are some scalar functions. It is useful to note that [44](46)∇Wxixi-xj,h=xi-xjxi-xj∂Wrij,h∂rij,where rij=xi-xj.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:(47a)dxidt=vi(47b)dvidt=ai,where vi is the particle velocity and ai 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](48a)vit+Δt2=vit-Δt2+Δtait(48b)xit+Δt=xit+Δtvit+Δt2.It is similar to the velocity Verlet algorithm, which can be expressed as [44](49a)xit+Δt=xit+Δtvit+Δt2ait2(49b)vit+Δt=vit+Δt2ait+ait+Δt.Both methods are second-order integration schemes. For numerical stability the time step Δt must obey the Courant-Friedrich-Levy condition [35](50)Δt≤0.25hvmaxand particle acceleration related condition(51)Δt≤0.25hamax.When viscous dissipation term is present, the following condition must be satisfied:(52)Δt≤0.125ρh2μ.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](53a)DρiDt=ρi∑jmjρjvi-vj·∇Wxixi-xj,h(53b)DρiDt=∑jmjvi-vj·∇Wxixi-xj,h.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 h2) [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:(54a)∇pi=∑jmjρjpi+pj∇Wxixi-xj,h(54b)∇pi=ρi∑jpjρj2+piρi2mj∇Wxixi-xj,h.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:(55)μ∇2vi=2μ∑jmjρjvi-vjxi-xj∂Wrij,h∂rij.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(56)μ∇2vi=∑j4μiμjμi+μjmjρjvi-vjxi-xj∂Wrij,h∂rij.Thus, the discretized momentum equation can be written as follows:(57)DviDt=-∑jpjρj2+piρi2mj∇Wxixi-xj,h+∑j4μiμjμi+μjmjρjρivi-vjxi-xj∂Wrij,h∂rij+g+Fsρi,where Fs 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 soundcs 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 (Section2.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 cf=1 and cg=-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(58)cri=∑jmjρjcjWxi-xj,h.Also, the color function can be calculated based on (6). The color function normal is calculated as(59a)nri=∑jmjρjcrj-cri∇Wxixi-xj,h.The interface local curvature κ given by (7) is calculated as(60)κri=1nri2∑jmjρjnrjnri-nrinrj·∇Wxixi-xj,h.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 Fs(ri) 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]:(61)∂2F∂xα∂xβi=∑jmjρj5xα,i-xα,jxβ,i-xβ,jxi-xj-δαβFi-Fj1xi-xj∂Wrij,h∂rij,where the Greek indexes α and β indicate spatial dimensions and δαβ 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 θ in the Verlet algorithm [48]. Both options are implemented in our model.


## 3.1.1. Smoothing Function

The smoothing functionW(x-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(x-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](36)Wx-x′,h=16πh312-3r2+3r3if0≤r<0.51-r3if0.5≤r<10ifr≥1,where r=x-x′/h. It is clear that W(x-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 massmi and density ρi associated with it. The particle volume can be defined as(37)Vi=miρi.The continuous integrals (34) and (35) can be approximated by their discrete analogues,(38a)fi=∑jfjWxi-xj,hVj=∑jfjmjρjWxi-xj,h(38b)∇fi=∑jfj∇Wxi-xj,hVj=∑jfjmjρj∇Wxixi-xj,h,where the summation in function term (38a) is carried out for all neighbouring particles within the distance h from the ith particle including j=i; note that in the derivative term j≠i (38b). According to [44], the error introduced by the finite sum approximation of continuous integrals can be estimated as errintegral=O(Δ/h), where Δ is the length scale characterizing the average distance between the particles.Whenf is a density, (38a) is transformed to(39)ρi=∑jmjWxi-xj,h.Usually, the masses of SPH particles representing the same material are equal. Another useful quantity is the number density, which can be calculated as(40)ni=ρimi=∑jWxi-xj,hand 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](41)ρi=∑jmjWxi-xj,h∑jWxi-xj,h/nj.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]:(42)∇f=1F∇Ff-f∇F,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=ρ, and F=1/ρ, substitution of which in (42) yields the following expressions of SPH derivative term:(43a)∇fi=∑jfj-fimjρj∇Wxixi-xj,h(43b)∇fi=1ρi∑jfj-fimj∇Wxixi-xj,h(43c)∇fi=ρi∑jfjρj2+fiρi2mj∇Wxixi-xj,h.The calculation of a Laplacian of some functions is often required for the transport equation. It can be obtained as(44)∇2fi=∑jfjmjρj∇2Wxixi-xj,h.This formula also can be symmetrized just as SPH derivatives terms; however, even symmetrized formulations may suffer from ∇2Wxi(xi-xj,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]:(45)∇·F∇fi=∑jFj+Fifi-fjmjρjxi-xjxi-xj2∇Wxixi-xj,h,where f and F are some scalar functions. It is useful to note that [44](46)∇Wxixi-xj,h=xi-xjxi-xj∂Wrij,h∂rij,where rij=xi-xj.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:(47a)dxidt=vi(47b)dvidt=ai,where vi is the particle velocity and ai 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](48a)vit+Δt2=vit-Δt2+Δtait(48b)xit+Δt=xit+Δtvit+Δt2.It is similar to the velocity Verlet algorithm, which can be expressed as [44](49a)xit+Δt=xit+Δtvit+Δt2ait2(49b)vit+Δt=vit+Δt2ait+ait+Δt.Both methods are second-order integration schemes. For numerical stability the time step Δt must obey the Courant-Friedrich-Levy condition [35](50)Δt≤0.25hvmaxand particle acceleration related condition(51)Δt≤0.25hamax.When viscous dissipation term is present, the following condition must be satisfied:(52)Δt≤0.125ρh2μ.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](53a)DρiDt=ρi∑jmjρjvi-vj·∇Wxixi-xj,h(53b)DρiDt=∑jmjvi-vj·∇Wxixi-xj,h.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 h2) [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:(54a)∇pi=∑jmjρjpi+pj∇Wxixi-xj,h(54b)∇pi=ρi∑jpjρj2+piρi2mj∇Wxixi-xj,h.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:(55)μ∇2vi=2μ∑jmjρjvi-vjxi-xj∂Wrij,h∂rij.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(56)μ∇2vi=∑j4μiμjμi+μjmjρjvi-vjxi-xj∂Wrij,h∂rij.Thus, the discretized momentum equation can be written as follows:(57)DviDt=-∑jpjρj2+piρi2mj∇Wxixi-xj,h+∑j4μiμjμi+μjmjρjρivi-vjxi-xj∂Wrij,h∂rij+g+Fsρi,where Fs 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 soundcs 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 (Section2.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 cf=1 and cg=-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(58)cri=∑jmjρjcjWxi-xj,h.Also, the color function can be calculated based on (6). The color function normal is calculated as(59a)nri=∑jmjρjcrj-cri∇Wxixi-xj,h.The interface local curvature κ given by (7) is calculated as(60)κri=1nri2∑jmjρjnrjnri-nrinrj·∇Wxixi-xj,h.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 Fs(ri) 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]:(61)∂2F∂xα∂xβi=∑jmjρj5xα,i-xα,jxβ,i-xβ,jxi-xj-δαβFi-Fj1xi-xj∂Wrij,h∂rij,where the Greek indexes α and β indicate spatial dimensions and δαβ 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 θ 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 ofN points of mass mp can be found from(62)rCM=1MB∑i=1Nmpri.The total external force acting on the rigid body is given by(63)FB=∑i=1Nmpai.The tensor of the moment of inertia is obtained from(64)IB=∑i=1Nmpri-rCM2I-ri-rCMri-rCMT.The total torque acting on the body can be calculated as(65)TB=∑i=1Nri-rCM×mpai,where ai 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 aCM and the angular acceleration of the rigid body αB are calculated as [27](66a)aCMt=FBtMB(66b)αBt=IBt-1TBt-ωBt-Δt2×IBtωBt-Δt2.Then, the velocity of the center of mass vCM and the angular velocity of the rigid body ωB are updated [27](67a)vCMt+Δt2=vCMt-Δt2+aCMtΔt(67b)ωBt+Δt2=ωBt-Δt2+αBtΔt.Finally, the position of the center of mass is updated by [27](68)rCMt+Δt=rCMt+vCMt+Δt2Δt.The velocities and the positions of the rigid body points can be obtained as follows [27](69a)vit+Δt2=vCMt+Δt2+ωBt+Δt2×rit-rCMt(69b)rit+Δt=rCMt+Δt+RBt+Δt2rit-rCMt,where RB is the rotation matrix. It can be calculated as [27](70)RBt+Δt2=RxωB,xt+Δt2ΔtRyωB,yt+Δt2ΔtRzωB,zt+Δt2Δt,where Rx(θ), Ry(θ), and Rz(θ) are the matrices which yield the rotation about principal coordinate axes(71)Rxθ=1000cos⁡θ-sin⁡θ0sin⁡θcos⁡θ(72)Ryθ=cos⁡θ0sin⁡θ010-sin⁡θ0cos⁡θ(73)Rzθ=cos⁡θ-sin⁡θ0sin⁡θcos⁡θ0001.


## 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 lengthh. 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:(74)csolid=Ksρ0,solid,where Ks is the bulk modulus of the solid phase and ρ0,solid the initial solid phase density. Ks for the elastic phase is defined by its elastic properties, and Ks 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]:(75)fi=-mjmi+mjBrini,where fi is the force acting on particle i, B is the symmetric function which defines the interaction strength, and ni 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 v^i for update of the SPH particles positions in the adopted numerical integration scheme instead of the usual particle velocity vi which is used in the momentum equation. The corrected velocity is calculated by [41](76)v^i=vi+ϵ∑jmj2vj-viρi+ρj,where ϵ 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 wallsx=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](77)vx,t=ρg2μxL-x+VxL+∑n=1∞ρgμL22-1n-1nπ3-2V-1n+1nπ×exp⁡-nπ2μtρL2sin⁡nπxL.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 with5·10-5m 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-3m, the body force is g=10-3m/s2, and the imposed velocity is V=2·10-4m/s. The fluid density and viscosity are ρ=1000kg/m3 and μ=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 Figure2. 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).Figure 2
Transient Poiseuille-Couette flow. Velocity profile evolution with time. Velocity profiles correspond tot=10-2s, t=4⋅10-2s, t=8⋅10-2s, t=16⋅10-2s, and t=1s.


### 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 permeabilityK of a porous medium can be found from Darcy’s law(78)v=-Kμ∇p-ρg,where μ is the viscosity, v is the mean fluid seepage velocity, ∇p is the pressure gradient, and g is the body force. For a spatially periodic simple cubic array of spheres, the permeability can be found analytically [54, 55] as(79)Ka=τ06πRk∗,where τ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]:(80)k∗=6πμRvF.Values of k∗ can be found in [54, 55] for various values of solid volume fraction ϕ(81)ϕ=VsolidVcell=4πR33L3.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 Knum 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 isL=20h, where h=10-4 m is the smoothing length. The center of the solid sphere rs 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 g=(0,0,5∗10-3)m/s2, the fluid density is ρ=1000kg/m3, and viscosity is μ=10-3 Pa·s. The speed of sound is cs=0.1m/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ϕ values, three values of solid concentration are considered ϕ=0.343, ϕ=0.45, and ϕ=0.5236. Coarse solid spheres are discretized with elementary cubes of size lc=h; that is, each elementary cube contains 8 solid SPH particles. The elementary cube is considered solid if its center rc satisfies the condition rc-rs≤R. In order to test the influence of the sphere surface discretization, spheres are also discretized with solid SPH particles which satisfy the condition rp-rs≤R, where rp is the location of the particle. This way provides a refined discretization and equivalent to discretizing the sphere with elementary solid cubes of size lc=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 ϕ=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.Figure 3
Coarse and refined discretizations of the sphere forϕ=0.45 (a central cut).Numerical errors for the calculated permeability values are summarized in Table2. 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.Table 2
Numerical error for permeability of periodic arrays of spheres for variousϕ values, solid sphere discretizations, and interparticle spacing as compared to Ka.
ϕ Discretization K a, m2 e r r , % (spacing h/2) e r r , % (spacing h/3) 0.343 Coarse 3.1728 ∗ 1 0 - 8 - 37.47 -21.26 0.45 Coarse 1.5904 ∗ 1 0 - 8 - 34.19 -17.22 0.5236 Coarse 1.0107 ∗ 1 0 - 8 - 37.47 -21.31  0.343 Refined 3.1728 ∗ 1 0 - 8 - 31.03 -8.32 0.45 Refined 1.5904 ∗ 1 0 - 8 - 30.11 - 9.24 0.5236 Refined 1.0107 ∗ 1 0 - 8 - 33.66 - 9.80It must be noted (see Table3) 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 ϕ=0.343, L=40h, h/3 setup is only 7.02% accurate compared to 3.97% accuracy for ϕ=0.343, L=20h, h/4 setup, in spite of larger total number of SPH particles and better sphere surface discretization.Table 3
The calculated permeabilityKnum of periodic array of spheres for ϕ=0.343, L=20h with refined discretization.
Interparticle spacing K n u m, m2 err, %. h/2 1.1115 ∗ 1 0 - 8 - 30.11 h/3 1.4434 ∗ 1 0 - 8 - 9.24 h/4 1.5272 ∗ 1 0 - 8 - 3.97Based 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 withh/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 atx=-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:(82a)vgx,t=-Fl22μgx2l2-xlμg-μfμg+μf-2μgμg+μf,-h≤x≤0(82b)vfx,t=-Fl22μfx2l2-xlμg-μfμg+μf-2μfμg+μf,0≤x≤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 alongy- and z-axis. The cubic spline smoothing function is used with smoothing length h=10-4m, the gas/liquid region width is l=5·10-5m, and the body force is equal to F=10-1Pa/m. The liquid density is ρf=103kg/m3 and viscosity μf=10-2 Pa·s; the gas density is ρg=1kg/m3 and viscosity μg=10-4 Pa·s.The numerical velocity profile is compared to the analytical solution in Figure4. 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.Figure 4
Two-phase Poiseuille flow. Analytical solution (black solid line) for velocity profile is compared to the SPH numerical solution (points).


### 4.4. Laplace Law

In order to test the continuum force algorithm described in Section2.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 Lcell=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=L3/4π1/3. After the equilibrium state has been obtained, the pressure difference Δp=pg-pf 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.Figure 5
Laplace’s law for spherical bubbles. Pressure step versush/R.


### 4.5. Reaction of Elastic Phase to Fluid Pressure

The elastic-fluid coupling can be tested by calculation of the relative volume variationdV/V of an elastic cube due to the pressure step dp applied by the surrounding fluid:(83)Ks=-dpdV/V.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 ε~0.1. However, from Table 1 it follows that clayey matrix 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.Figure 6
Elastic cube volume change due to applied fluid pressure.


## 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 wallsx=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](77)vx,t=ρg2μxL-x+VxL+∑n=1∞ρgμL22-1n-1nπ3-2V-1n+1nπ×exp⁡-nπ2μtρL2sin⁡nπxL.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 with5·10-5m 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-3m, the body force is g=10-3m/s2, and the imposed velocity is V=2·10-4m/s. The fluid density and viscosity are ρ=1000kg/m3 and μ=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 Figure2. 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).Figure 2
Transient Poiseuille-Couette flow. Velocity profile evolution with time. Velocity profiles correspond tot=10-2s, t=4⋅10-2s, t=8⋅10-2s, t=16⋅10-2s, and t=1s.


## 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 permeabilityK of a porous medium can be found from Darcy’s law(78)v=-Kμ∇p-ρg,where μ is the viscosity, v is the mean fluid seepage velocity, ∇p is the pressure gradient, and g is the body force. For a spatially periodic simple cubic array of spheres, the permeability can be found analytically [54, 55] as(79)Ka=τ06πRk∗,where τ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]:(80)k∗=6πμRvF.Values of k∗ can be found in [54, 55] for various values of solid volume fraction ϕ(81)ϕ=VsolidVcell=4πR33L3.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 Knum 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 isL=20h, where h=10-4 m is the smoothing length. The center of the solid sphere rs 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 g=(0,0,5∗10-3)m/s2, the fluid density is ρ=1000kg/m3, and viscosity is μ=10-3 Pa·s. The speed of sound is cs=0.1m/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ϕ values, three values of solid concentration are considered ϕ=0.343, ϕ=0.45, and ϕ=0.5236. Coarse solid spheres are discretized with elementary cubes of size lc=h; that is, each elementary cube contains 8 solid SPH particles. The elementary cube is considered solid if its center rc satisfies the condition rc-rs≤R. In order to test the influence of the sphere surface discretization, spheres are also discretized with solid SPH particles which satisfy the condition rp-rs≤R, where rp is the location of the particle. This way provides a refined discretization and equivalent to discretizing the sphere with elementary solid cubes of size lc=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 ϕ=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.Figure 3
Coarse and refined discretizations of the sphere forϕ=0.45 (a central cut).Numerical errors for the calculated permeability values are summarized in Table2. 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.Table 2
Numerical error for permeability of periodic arrays of spheres for variousϕ values, solid sphere discretizations, and interparticle spacing as compared to Ka.
ϕ Discretization K a, m2 e r r , % (spacing h/2) e r r , % (spacing h/3) 0.343 Coarse 3.1728 ∗ 1 0 - 8 - 37.47 -21.26 0.45 Coarse 1.5904 ∗ 1 0 - 8 - 34.19 -17.22 0.5236 Coarse 1.0107 ∗ 1 0 - 8 - 37.47 -21.31  0.343 Refined 3.1728 ∗ 1 0 - 8 - 31.03 -8.32 0.45 Refined 1.5904 ∗ 1 0 - 8 - 30.11 - 9.24 0.5236 Refined 1.0107 ∗ 1 0 - 8 - 33.66 - 9.80It must be noted (see Table3) 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 ϕ=0.343, L=40h, h/3 setup is only 7.02% accurate compared to 3.97% accuracy for ϕ=0.343, L=20h, h/4 setup, in spite of larger total number of SPH particles and better sphere surface discretization.Table 3
The calculated permeabilityKnum of periodic array of spheres for ϕ=0.343, L=20h with refined discretization.
Interparticle spacing K n u m, m2 err, %. h/2 1.1115 ∗ 1 0 - 8 - 30.11 h/3 1.4434 ∗ 1 0 - 8 - 9.24 h/4 1.5272 ∗ 1 0 - 8 - 3.97Based 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 withh/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 atx=-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:(82a)vgx,t=-Fl22μgx2l2-xlμg-μfμg+μf-2μgμg+μf,-h≤x≤0(82b)vfx,t=-Fl22μfx2l2-xlμg-μfμg+μf-2μfμg+μf,0≤x≤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 alongy- and z-axis. The cubic spline smoothing function is used with smoothing length h=10-4m, the gas/liquid region width is l=5·10-5m, and the body force is equal to F=10-1Pa/m. The liquid density is ρf=103kg/m3 and viscosity μf=10-2 Pa·s; the gas density is ρg=1kg/m3 and viscosity μg=10-4 Pa·s.The numerical velocity profile is compared to the analytical solution in Figure4. 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.Figure 4
Two-phase Poiseuille flow. Analytical solution (black solid line) for velocity profile is compared to the SPH numerical solution (points).


## 4.4. Laplace Law

In order to test the continuum force algorithm described in Section2.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 Lcell=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=L3/4π1/3. After the equilibrium state has been obtained, the pressure difference Δp=pg-pf 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.Figure 5
Laplace’s law for spherical bubbles. Pressure step versush/R.


## 4.5. Reaction of Elastic Phase to Fluid Pressure

The elastic-fluid coupling can be tested by calculation of the relative volume variationdV/V of an elastic cube due to the pressure step dp applied by the surrounding fluid:(83)Ks=-dpdV/V.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 ε~0.1. However, from Table 1 it follows that clayey matrix 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.Figure 6
Elastic cube volume change due to applied fluid pressure.


## 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.


### 5.1. Characterization of the Sample

The COx claystone sample EST27405-1 is represented by a series of 103 images of771×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×10nm3. The images are treated based on gray level in order to recognize the pore-space and clayey matrix. The porosity of the sample is ε=0.0307±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 of385×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.Figure 7
Largest percolating pore components inside the Callovo-Oxfordian clay sample.Due to significant computational costs (one simulation for54h×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 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.Figure 8
(a) A disconnected pore componentε=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.
(a) (b)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 Figure8(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 Hr=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,Fz) is applied to the gas in the reservoir outside the sample, which creates a pressure difference between upper z=60 (pg=HrFzρg) and lower z=10 (pg=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×54h×70h size is given in Figure 8(b). One can see that the gas phase (green) penetrates the sample 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 (Δ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:(84a)Re=ρgVgLμg(84b)M=μwμg(84c)Ca=μgVgσ(84d)Bo=ρw-ρggL2σ.The Reynolds number Re which represents the ratio between the inertial and viscous forces is very small Re≪1 for gas flow in clay due to small flow velocity Vg 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≪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.Table 4
Dimensionless numbers characterizing two-phase water-gas flow in clay microfractures.
Number Re M Ca Bo Expected value 4 ∗ 1 0 - 5 1.1 ∗ 1 0 2 7.2 ∗ 1 0 - 10 4.7 ∗ 1 0 - 6 Simulation value ≈ 0.5 1.1 ∗ 1 0 2 ≈ 1.25 ∗ 1 0 - 4 0It can be observed that the Reynolds and capillary numbers observed in our simulations (for mean gas velocityVg≈5∗10-5m/s before the percolation) are several orders of magnitude higher than the expected real ones as a consequence of significantly reduced viscosity (μw=10-7 Pa·s, μg=10-9 Pa·s) and surface tension (σ=4∗10-10 N·m−1) 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 valuecs=0.5 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 ΔP = 0.3 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 Pgas,in = 15 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 Pgas,in and simulation Pgas,in = 0.3 Pa pressure values at the sample entry is 5∗107. 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 Rp between the gas pressure at the sample entry Pgas,in=0.3Pa and Young’s modulus E=60Pa of the elastic phase used in the simulation(85)Rp=Pgas,inEstays the same as in repository conditions where Rp≈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∗107). These rescaled values will be preceded by ~.


#### 5.2.1. Drainage with and without Damage

The pressure differenceΔP~15MPa is imposed between upper and lower sample sides, while the confining stress of σconf~6MPa 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 ratioRp=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.Figure 9
LiquidSw (broken lines) and gas Sg (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 σconf~6MPa with (magenta) and without (black) damage applied to the elastic phase. (b) Elastic matrix with application of damage model under confining stresses of σconf~10MPa (red), σconf~8MPa (green), σconf~6MPa (magenta), and σconf~4MPa (blue).
(a) (b)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 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 (Figure9(a)) increases the drainage rate at all the stages of simulation and the nonwetting fluid percolation time is accelerated (from 8.74∗10-3 s to 7.58∗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σconf~10MPa, 8 MPa, 6 MPa, and 4 MPa applied to lateral sides and for the same pressure difference ΔP~15MPa.Same as for (Figure9(a)), the initial compression phase can be observed in Figure 9(b). It is more pronounced for σconf~4MPa 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.Figure 10
(a) Number of damaged elastic SPH particles. (b) The total damage calculated as a sum over all the damaged points.
(a) (b)The percolation times and gas saturation at percolation are summarized in Table5. 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 tperc 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.Table 5
Percolation timestperc and gas saturation Sg,perc at percolation for various confining stresses.
σ c o n f, ~MPa t p e r c, s S g , p e r c 10 9.18 ∗ 1 0 - 3 0.5158 8 8.52 ∗ 1 0 - 3 0.5074 6 7.58 ∗ 1 0 - 3 0.5043 4 6.56 ∗ 1 0 - 3 0.4971  6, no damage 8.74 ∗ 1 0 - 3 0.5039  Incompressible porous matrix 9.84 ∗ 1 0 - 3 0.5159Evolution of the number of damaged points and of the total damage is presented in Figure10. 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.Figure 11
A cross-section (SPH points inside the45h<z<46h range) of the unit cell under confining stress σconf~8MPa. 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.
(a) (b) (c) (d)The seepage velocities of fluid phases during drainage are plotted in Figure12. 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.Figure 12
Water (a) and gas (b) seepage velocities during drainage as functions of time for various confining stresses.
(a) (b)Figure 13
(a) Water and gas mean pressures inside the sample as functions of time. (b) Capillary pressure as a function of gas saturation.
(a) (b)The waterPw and gas Pg pressures averaged over the percolating pore volume and the capillary pressure Pc=Pg-Pw 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 Sg=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 Figure13 are about 3.6 times smaller than those which can be expected for this pore size in the case of water/air interface with σ=7.2∗10-2N·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 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 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)).Figure 14
A cross-section (SPH points inside the58h<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 Rp=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.
(a) (b) (c) (d)When using SPH model for the elastic phase, a fracturing can be observed even without applying a damage model. For the simulations presented in Sections5.2.1 and 5.2.2 the ratio Rp=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 be no initiation of new fractures. However, for larger values, for example, Rp=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 may result in formation of narrow fractures as displayed in Figure14(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 influence of the rigid inclusions on the flow conditions, two rigid spheres were incorporated into the elastic matrix as displayed in Figure15. 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.Figure 15
A cross-section (a)24h<z<25h and (b) 33h<x<34h of the unit cell under confining stress σconf~6MPa at t = 12 s. Two spherical rigid inclusions are incorporated into initially elastic solid matrix.
(a) (b)


## 5.1. Characterization of the Sample

The COx claystone sample EST27405-1 is represented by a series of 103 images of771×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×10nm3. The images are treated based on gray level in order to recognize the pore-space and clayey matrix. The porosity of the sample is ε=0.0307±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 of385×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.Figure 7
Largest percolating pore components inside the Callovo-Oxfordian clay sample.Due to significant computational costs (one simulation for54h×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 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.Figure 8
(a) A disconnected pore componentε=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.
(a) (b)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 Figure8(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 Hr=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,Fz) is applied to the gas in the reservoir outside the sample, which creates a pressure difference between upper z=60 (pg=HrFzρg) and lower z=10 (pg=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×54h×70h size is given in Figure 8(b). One can see that the gas phase (green) penetrates the sample 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 (Δ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:(84a)Re=ρgVgLμg(84b)M=μwμg(84c)Ca=μgVgσ(84d)Bo=ρw-ρggL2σ.The Reynolds number Re which represents the ratio between the inertial and viscous forces is very small Re≪1 for gas flow in clay due to small flow velocity Vg 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≪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.Table 4
Dimensionless numbers characterizing two-phase water-gas flow in clay microfractures.
Number Re M Ca Bo Expected value 4 ∗ 1 0 - 5 1.1 ∗ 1 0 2 7.2 ∗ 1 0 - 10 4.7 ∗ 1 0 - 6 Simulation value ≈ 0.5 1.1 ∗ 1 0 2 ≈ 1.25 ∗ 1 0 - 4 0It can be observed that the Reynolds and capillary numbers observed in our simulations (for mean gas velocityVg≈5∗10-5m/s before the percolation) are several orders of magnitude higher than the expected real ones as a consequence of significantly reduced viscosity (μw=10-7 Pa·s, μg=10-9 Pa·s) and surface tension (σ=4∗10-10 N·m−1) 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 valuecs=0.5 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 ΔP = 0.3 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 Pgas,in = 15 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 Pgas,in and simulation Pgas,in = 0.3 Pa pressure values at the sample entry is 5∗107. 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 Rp between the gas pressure at the sample entry Pgas,in=0.3Pa and Young’s modulus E=60Pa of the elastic phase used in the simulation(85)Rp=Pgas,inEstays the same as in repository conditions where Rp≈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∗107). These rescaled values will be preceded by ~.


### 5.2.1. Drainage with and without Damage

The pressure differenceΔP~15MPa is imposed between upper and lower sample sides, while the confining stress of σconf~6MPa 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 ratioRp=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.Figure 9
LiquidSw (broken lines) and gas Sg (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 σconf~6MPa with (magenta) and without (black) damage applied to the elastic phase. (b) Elastic matrix with application of damage model under confining stresses of σconf~10MPa (red), σconf~8MPa (green), σconf~6MPa (magenta), and σconf~4MPa (blue).
(a) (b)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 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 (Figure9(a)) increases the drainage rate at all the stages of simulation and the nonwetting fluid percolation time is accelerated (from 8.74∗10-3 s to 7.58∗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σconf~10MPa, 8 MPa, 6 MPa, and 4 MPa applied to lateral sides and for the same pressure difference ΔP~15MPa.Same as for (Figure9(a)), the initial compression phase can be observed in Figure 9(b). It is more pronounced for σconf~4MPa 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.Figure 10
(a) Number of damaged elastic SPH particles. (b) The total damage calculated as a sum over all the damaged points.
(a) (b)The percolation times and gas saturation at percolation are summarized in Table5. 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 tperc 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.Table 5
Percolation timestperc and gas saturation Sg,perc at percolation for various confining stresses.
σ c o n f, ~MPa t p e r c, s S g , p e r c 10 9.18 ∗ 1 0 - 3 0.5158 8 8.52 ∗ 1 0 - 3 0.5074 6 7.58 ∗ 1 0 - 3 0.5043 4 6.56 ∗ 1 0 - 3 0.4971  6, no damage 8.74 ∗ 1 0 - 3 0.5039  Incompressible porous matrix 9.84 ∗ 1 0 - 3 0.5159Evolution of the number of damaged points and of the total damage is presented in Figure10. 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.Figure 11
A cross-section (SPH points inside the45h<z<46h range) of the unit cell under confining stress σconf~8MPa. 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.
(a) (b) (c) (d)The seepage velocities of fluid phases during drainage are plotted in Figure12. 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.Figure 12
Water (a) and gas (b) seepage velocities during drainage as functions of time for various confining stresses.
(a) (b)Figure 13
(a) Water and gas mean pressures inside the sample as functions of time. (b) Capillary pressure as a function of gas saturation.
(a) (b)The waterPw and gas Pg pressures averaged over the percolating pore volume and the capillary pressure Pc=Pg-Pw 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 Sg=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 Figure13 are about 3.6 times smaller than those which can be expected for this pore size in the case of water/air interface with σ=7.2∗10-2N·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 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 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)).Figure 14
A cross-section (SPH points inside the58h<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 Rp=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.
(a) (b) (c) (d)When using SPH model for the elastic phase, a fracturing can be observed even without applying a damage model. For the simulations presented in Sections5.2.1 and 5.2.2 the ratio Rp=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 be no initiation of new fractures. However, for larger values, for example, Rp=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 may result in formation of narrow fractures as displayed in Figure14(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 influence of the rigid inclusions on the flow conditions, two rigid spheres were incorporated into the elastic matrix as displayed in Figure15. 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.Figure 15
A cross-section (a)24h<z<25h and (b) 33h<x<34h of the unit cell under confining stress σconf~6MPa at t = 12 s. Two spherical rigid inclusions are incorporated into initially elastic solid matrix.
(a) (b)


## 5.2.1. Drainage with and without Damage

The pressure differenceΔP~15MPa is imposed between upper and lower sample sides, while the confining stress of σconf~6MPa 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 ratioRp=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.Figure 9
LiquidSw (broken lines) and gas Sg (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 σconf~6MPa with (magenta) and without (black) damage applied to the elastic phase. (b) Elastic matrix with application of damage model under confining stresses of σconf~10MPa (red), σconf~8MPa (green), σconf~6MPa (magenta), and σconf~4MPa (blue).
(a) (b)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 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 (Figure9(a)) increases the drainage rate at all the stages of simulation and the nonwetting fluid percolation time is accelerated (from 8.74∗10-3 s to 7.58∗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σconf~10MPa, 8 MPa, 6 MPa, and 4 MPa applied to lateral sides and for the same pressure difference ΔP~15MPa.Same as for (Figure9(a)), the initial compression phase can be observed in Figure 9(b). It is more pronounced for σconf~4MPa 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.Figure 10
(a) Number of damaged elastic SPH particles. (b) The total damage calculated as a sum over all the damaged points.
(a) (b)The percolation times and gas saturation at percolation are summarized in Table5. 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 tperc 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.Table 5
Percolation timestperc and gas saturation Sg,perc at percolation for various confining stresses.
σ c o n f, ~MPa t p e r c, s S g , p e r c 10 9.18 ∗ 1 0 - 3 0.5158 8 8.52 ∗ 1 0 - 3 0.5074 6 7.58 ∗ 1 0 - 3 0.5043 4 6.56 ∗ 1 0 - 3 0.4971  6, no damage 8.74 ∗ 1 0 - 3 0.5039  Incompressible porous matrix 9.84 ∗ 1 0 - 3 0.5159Evolution of the number of damaged points and of the total damage is presented in Figure10. 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.Figure 11
A cross-section (SPH points inside the45h<z<46h range) of the unit cell under confining stress σconf~8MPa. 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.
(a) (b) (c) (d)The seepage velocities of fluid phases during drainage are plotted in Figure12. 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.Figure 12
Water (a) and gas (b) seepage velocities during drainage as functions of time for various confining stresses.
(a) (b)Figure 13
(a) Water and gas mean pressures inside the sample as functions of time. (b) Capillary pressure as a function of gas saturation.
(a) (b)The waterPw and gas Pg pressures averaged over the percolating pore volume and the capillary pressure Pc=Pg-Pw 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 Sg=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 Figure13 are about 3.6 times smaller than those which can be expected for this pore size in the case of water/air interface with σ=7.2∗10-2N·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 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 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)).Figure 14
A cross-section (SPH points inside the58h<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 Rp=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.
(a) (b) (c) (d)When using SPH model for the elastic phase, a fracturing can be observed even without applying a damage model. For the simulations presented in Sections5.2.1 and 5.2.2 the ratio Rp=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 be no initiation of new fractures. However, for larger values, for example, Rp=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 may result in formation of narrow fractures as displayed in Figure14(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 influence of the rigid inclusions on the flow conditions, two rigid spheres were incorporated into the elastic matrix as displayed in Figure15. 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.Figure 15
A cross-section (a)24h<z<25h and (b) 33h<x<34h of the unit cell under confining stress σconf~6MPa at t = 12 s. Two spherical rigid inclusions are incorporated into initially elastic solid matrix.
(a) (b)


## 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 indurated 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 simulating 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 situations 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 propagation 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 computationally expensive. The main load comes from calculation of neighbouring particles contributions when evaluating densities 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 performed relatively easily, since densities and forces are calculated by summing contributions from neighbouring particles; therefore, each particle can be treated separately by a computational thread.For the present study, the OpenMP and CUDA GPU parallelization 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 permeability values with a reasonable precision (about 30%).The model was then applied to simulate drainage in a small3D 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 number 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 relationships obtained for macroscopic samples. One remarkable feature of capillary pressure curves in Figure13(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 demonstrated 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 inclusions 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 geometries in conditions close to the real ones. To be able to obtain more quantitative information about drainage in argillites 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.



---

*Source: 6873298-2018-03-11.xml*