OR/18/049 Modelling approaches

Tamayo-Mas, E, Harrington, J F, Brüning, T, Kolditz, O, Shao, H, Dagher, E E, Lee, J, Kim, K, Rutqvist, J, Lai, S H, Chittenden, N, Wang, Y, Damians, I P, Olivella, S. 2018. DECOVALEX-2019 project: Task A - modElliNg Gas INjection ExpERiments (ENGINEER). Nottingham, UK, British geological Survey. (OR/18/049).

Different approaches have been adopted by the participating teams:

1. Classical two-phase flow models, where basic physical principles such as mass and momentum balance apply for each phase. These standard models have been coupled to different mechanical deformation behaviours in order to better represent the experimentally observed flow. In particular, the following deformation models have been considered:

- Model UPC/Andra-H, with a rigid porous medium (developed by UPC/Andra).

- Model LBNL-C, where a linear elastic medium is considered (developed by LBNL).

- Model CNSC-E, where an elastic medium is assumed (developed by CNSC).

- Model CNSC-D, where a damage medium is assumed (also developed by CNSC).

- Model KAERI, where a similar damage medium is assumed (developed by KAERI).

- Model BGR/UFZ, with an elasto-plastic medium (developed by BGR/UFZ).

- Model NCU/TPC, with a viscoelastic model (developed by NCU/TPC).

2. Enriched two-phase flow models, where preferential pathways are considered:
   

   - Model Quintessa/RWM: separate gas and water coexist within an elastic deformation matrix (developed by Quintessa/RWM).

   - Model UPC/Andra-HM1: embedded fracture permeability model where deformation is modelled assuming elasticity with effective stress and suction terms (developed by UPC/Andra).

   - Model UPC/Andra-HM2: embedded fracture permeability model where deformation is modelled assuming elasticity with effective stress and without suction terms (developed by UPC/Andra).

3. Discrete approaches, where fractures are explicitly modelled:

   - Model LBNL-D: a two-phase flow model within a discrete fracture network (developed by LBNL).

   - Model SNL: a conceptual chaotic model developed by SNL.

Differences between the proposed strategies also lie in the software used by the teams (Table 3) and in the assumed geometry to represent the saturated bentonite (Table 4).

Table 3    Codes used by the participating teams.

Model	Software	Reference
BGR/UFZ	OpenGeoSys	Kolditz et al., 2012
CNSC-E	COMSOL Multiphysics ®	COMSOL, 2017
CNSC-D
KAERI	TOUGH-MP FLAC3D	Zhang et al., 2008, Itasca 2018
LBNL-C	TOUGH-FLAC	Rutqvist, 2011
LBNL-D	TOUGH-RBSN	Asahina et al., 2014
NCU/TPC	THMC 7.1	Yeh et al., 2013
Quintessa/RWM	QPAC	Quintessa, 2013
UPC/Andra-H	CODE_BRIGHT	Olivella et al., 1996
UPC/Andra-HM1
UPC/Andra-HM2

Table 4    Test geometries used by the teams.

Model	Geometry	Mesh	Discretisation
BGR/UFZ	2D axisymmetric	Triangle	2715 elements 1447 nodes
CNSC-E	3D	Tetrahedral	22000 elements
CNSC-D	3D	Tetrahedral	22000 elements
KAERI	3D	Hexahedral	7560 elements 8450 nodes
LBNL-C	3D	Hexahedral	23400 elements 24939 nodes
LBNL-D	2D	Cells and lattice	1401 cells 3840 lattice elements
NCU/TPC	2D	Quadrilateral	1250 elements 1326 nodes
Quintessa/RWM	1D finite volume	-	20 cells +1 abstract cell to simulate mass balance in the injector
UPC/Andra-H	3D	Hexahedral	7168 elements 7917 nodes
UPC/Andra-HM1
UPC/Andra-HM2

In order to calibrate and validate these strategies, teams used different parameters and calibration techniques. As seen in Table 5, almost all the teams used for their simulations the following given basic parameters:

1. Young's modulus: 𝐸𝐸 = 307 MPa

1. Poisson's ratio: 𝜈𝜈 = 0.4
2. Porosity: 𝜙𝜙 = 0.44
3. Intrinsic permeability: 𝑘𝑘 = 3.4 × 10⁻²¹ m²

However, in order to improve their calibrated results, some teams (BGR/UFZ, NCU/TPC, UPC/Andra) have used a Young’s modulus and initial permeability and porosity values significantly different from those of the material. Another difference amongst the proposed strategies is the amount of parameters that need to be calibrated to fit experimental results. As seen in Table 6, the complexity of the physical process requires a detailed (and sometimes cumbersome) calibration process.

Table 5    Basic fixed parameters used by the teams.

Model	Basic fixed parameters
	E [MPa]	𝝂𝝂 [-]	𝝓𝝓 [-]	k [m2]
BGR/UFZ	307	0.4	0.44	calibrated
CNSC-E	307	0.4	0.44	3.4 × 10−21
CNSC-D	307	0.4	0.44	3.4 × 10−21
KAERI	307	0.4	0.44	3.4 × 10−21
LBNL-C	307	0.4	0.44	3.4 × 10−21
LBNL-D	307	0.4	Heterogeneous (𝜙𝜙0 = 0.44)	3.4 × 10−21
NCU/TPC	307	0.4	0.43	calibrated
Quintessa/RWM	307	0.4	-	3.4 × 10−21
UPC/Andra-H	-	-	Heterogeneous (𝜙𝜙0 = 0.36–0.38)	calibrated
UPC/Andra-HM1	225	0.125	Heterogeneous (𝜙𝜙0 = 0.36–0.38)	calibrated
UPC/Andra-HM2	307	0.44	0.41	calibrated

Table 6    Calibrated parameters used by the teams.

Model	Calibrated parameters
BGR/UFZ	* Intrinsic permeability: 1 × 10−21 m2 * Critical pressure
CNSC-E	Air-entry value
CNSC-D	* Air-entry value * Minimum air-entry value * Maximum intrinsic permeability * Damage smoothing parameter Swelling coefficient
KAERI	* Damage parameters: tensile strength, residual tensile strength, tensile strain limit, maximum damage value, maximum permeability
LBNL-C	* Hydro-mechanical coupling parameters: swelling coefficient, maximum aperture for stress, reference stress * Effective gas-entry pressure value
LBNL-D	* Mohr-Coulomb failure parameters * Effective gas-entry pressure value * Swelling coefficient
NCU/TPC	* Intrinsic permeability: 3.4 × 10−22 m2 * Air-entry value * Viscous parameters
Quintessa/RWM	* Capillarity compressibility * Capillary spacing * Swelling pressure * Biot coefficient
UPC/Andra-H	* Intrinsic permeability: Het. * 𝑘𝑘0 = 5.59 × 10−18 m2
UPC/Andra-HM1	* Intrinsic permeability: Het. * 𝑘𝑘0 = (2.15–7.10) × 10−18 m2 * Embedded permeability parameters
UPC/Andra-HM2	* Intrinsic permeability: Het. * 𝑘𝑘0 = (2.15–7.10) × 10−19 m2 * Embedded permeability parameters

Classical two-phase flow models

Different approaches have been proposed. All of them assume the standard van Genuchten (1980) model to describe the water retention curve but they differ in the assumed mechanical deformation behaviour.

Model UPC/Andra-H

UPC/Andra developed a two-phase flow calculation on a 3D rigid medium (model UPC/Andra-H). In this model, porosity 𝜙𝜙 [-] is assumed to be space-dependent and it is heterogeneously distributed (randomly generated in the range of 0.36–0.38). The intrinsic permeability 𝑘𝑘 [m²] is porosity-dependent and it is given by the exponential function

where 𝑘𝑘0, 𝜙𝜙0 and 𝑏𝑏𝑘𝑘 are equation definition parameters. Here the calibrated values 𝑘𝑘0 = 5.59´10^-18 m² and 𝑏𝑏𝑘𝑘 = 60 have been used. The retention curve is also porosity-dependent and it is computed as

where both the gas pressure 𝑝𝑝𝑔𝑔 [Pa] and the water pressure 𝑝𝑝𝑙𝑙 [Pa] are computed according to the van Genuchten model, 𝜆𝜆𝑣𝑣𝑣𝑣 is the shape function and 𝑝𝑝𝑣𝑣𝑣𝑣 [Pa] is the capillary pressure

with 𝑝𝑝0 [Pa] and 𝑏𝑏𝑣𝑣𝑣𝑣 two extra equation parameters. In this model

Model LBNL-C

In the poro-elastic model LBNL-C, the effective stress is expressed as

where σ' and σ are the effective and total stress tensors respectively, I is the identity tensor and the pore pressure 𝑝𝑝𝜙𝜙 is defined as

with 𝑝𝑝𝑙𝑙 and 𝑝𝑝𝑔𝑔 liquid and gas phase pressures respectively.

This model assumes that the bentonite behaves elastically, with a volumetric swelling and a swelling stress that depends on the changes in water saturation Δ𝑆𝑆𝑙𝑙, according to

where 𝝈𝝈′ [Pa] is the swelling stress, 𝐊𝐊 [Pa] is the bulk modulus, 𝑆𝑆 [-] is the liquid saturation and

This model assumes a fracture-like behaviour of the flow path. Hence, a pressure dependent permeability function is considered, where 𝑎𝑎 [m] is the element width and 𝑏𝑏ℎ [m] is a non-linear function of the effective minimum compressive stress that reads

with 𝑏𝑏ℎ0 [m] being the (calibrated) maximum aperture for permeability (𝑏𝑏ℎ0 = 2.9 𝜇𝜇m), 𝜎𝜎𝑛𝑛 [Pa] the total stress normal to the fracture and 𝜎𝜎𝑛𝑛, ref [Pa] the adjusted reference stress normal to the fracture (𝜎𝜎𝑛𝑛, ref = 0.2 MPa). The aperture versus pressure relationship of Eq. (8) corresponds to the Bandis et al. (1983) model and its parameters are calibrated by matching pressure and outflow responses observed in the experiments.

Model CNSC-E

The three-dimensional model CNSC-E proposed by CNSC assumes an elastic deformation of the solid material. In this model, the presence of the pore gas phase is accounted by defining the average pore fluid pressure

where 𝜒𝜒 [-] is a function of the matric suction 𝑠𝑠 [Pa] and the air entry suction 𝑠𝑠𝑒𝑒 [Pa] (here calibrated to 𝑠𝑠𝑒𝑒 = 9.81 MPa) determined from the Khalili and Khabbaz (1998) equation

Model CNSC-D

An enhanced version of the above-mentioned model is the damaged-based model CNSC-D, also proposed by CNSC. Model CNSC-D assumes a continuum elastic-damage bulk (Fall et al., 2014), where the rock may progressively degrade due to microcracks. Thus, according to continuum elastic-damage mechanics, the stiffness of the material depends on a damage parameter 𝐷𝐷 (0 ≤ 𝐷𝐷 ≤ 1) and reads

where 𝐸𝐸0 [Pa] is the initial stiffness and 𝑐𝑐 [-] (𝑐𝑐 = 10) is the calibrated damage smoothing parameter (note that both the rock and its damage are assumed isotropic and elastic and thus, 𝐸𝐸, 𝐸𝐸0 and 𝐷𝐷 are all scalar). The intrinsic permeability of this model is also damage dependent and it is computed by means of

where

with 𝑘𝑘max [m²] being a calibrated maximum intrinsic permeability (𝑘𝑘max = 1.00 × 10–18 m²). In this model, the swelling induced volume change d𝜀𝜀vs is related to the change in suction ds by the equation

where 𝐵𝐵𝑠𝑠 [Pa-1] is a calibrated swelling coefficient, Nasir et al. (2016).

Model KAERI

Another damage-based model was proposed by KAERI. Again, bentonite degradation due to microcracks is taken into account by means of the tensile/compressive damage model proposed by Fall et al. (2014). That is, damage variable under tension is described as

where 𝑓𝑓tr [Pa] is the calibrated residual tensile strength (𝑓𝑓tr = 0.2 MPa) and 𝜀𝜀𝑡𝑡0 [-] and 𝜀𝜀𝑡𝑡𝑡𝑡 [-] are the initial and final tensile strain thresholds respectively (𝜀𝜀𝑡𝑡𝑡𝑡 = 5 × 10^-3). At the same time, damage variable under compression reads

where 𝑓𝑓cr [Pa] is the residual compressive strength (𝑓𝑓cr = 3 MPa) and 𝜀𝜀co [-] is the compressive threshold strain. The intrinsic permeability of this model is also characterised by Eq. (12) and (13). In this case, the calibrated maximum intrinsic permeability is assumed to be one order of magnitude smaller than in model CNSC-D (𝑘𝑘max = 1.00 × 10^-19 m²).

Model BGR/UFZ

In model BGR/UFZ, a critical pressure 𝑝𝑝crit [Pa] is introduced and calibrated, which is a sum of the minimal principal stress (confining pressure) and the gas entry pressure. The dilatancy pathway is expressed by permeability change (Xu et al., 2013)

Model NCU/TPC

In model NCU/TPC, a viscoelastic deformation behaviour is assumed. In this model, intrinsic permeability is a non-linear function of porosity and reads

where 𝑘𝑘0 [m²] and 𝜙𝜙0 [-] are the reference intrinsic permeability and the reference porosity respectively and 𝑛𝑛 is the fractional exponent depending on the particle size and packing structure. Note that at this moment in time, see Section 4, this model is not able to correctly reflect the physics of the experiment and further development is needed.

Enriched two-phase flow models with preferential pathways

Three different enriched two-phase flow continuum models have been considered. In these approaches, the bentonite matrix behaves as a linear elastic medium, where preferential gas pathways are included.

Model Quintessa/RWM

Quintessa/RWM developed a model with separate gas and water pathways. The model considers the saturated bentonite as two components: (i) the clay solid ‘grains’ with the non-mobile interlayer water and (ii) the ‘free’ water component. Darcy’s law is assumed to describe the ‘free’ water movement with water permeability 𝑘𝑘𝑤𝑤 [m²]

whereas gas movement is modelled through the Hagen-Poiseuille law. Thus, the opening (and closure) of gas pathways is represented through the capillary radius and its relationship to gas permeability 𝑘𝑘𝑔𝑔 [m²]

where 𝑎𝑎 [m²] is the calibrated capillary spacing (𝑎𝑎 = 5.66 × 10^-6 m²). The capillary radius 𝑟𝑟 [m] is considered to be dependent on the capillarity compressibility 𝛾𝛾 [m2Pa-1] (assumed to be 𝛾𝛾 = 1 × 10^-20 m²Pa^-1) through

where 𝑟𝑟0 [m] is the reference capillary radius (𝑟𝑟0 = 0.05 𝜇𝜇m), 𝜖𝜖𝑐𝑐 [Pa] is the stress for the capillary compressibility, 𝜖𝜖𝑐𝑐0 [Pa] is the reference stress for the capillary compressibility (𝜖𝜖𝑐𝑐0 = 0 MPa) and 𝑝𝑝𝑔𝑔 [Pa] is the gas pressure. As seen by means of Eq. (20) and (21), the coupling of the stresses to the permeability is done through the capillarity radius. In order to model the ceasing of gas flow observed in Figure 8, Quintessa/RWM proposed a simple model of the gas injection system based on the ideal gas law.

Model UPC/Andra-HM1

UPC/Andra-HM1 is the first hydro-mechanical model developed by UPC/Andra. This model is based on the embedded fracture model proposed by Olivella and Alonso, 2008. The basic idea of this model consists in properly representing single fractures embedded in a continuous matrix. These fractures are characterized by their aperture 𝑏𝑏 [m] and spacing 𝑎𝑎 [m] thus leading to an intrinsic permeability

where 𝜙𝜙0 [-] is the initial porosity, 𝑘𝑘0 [m²] is the reference permeability and 𝜙𝜙 [-] is the porosity. As done with the hydraulic model UPC/Andra-H, porosity is assumed to be space-dependent and it is heterogeneously distributed (randomly generated with 𝜙𝜙0 taking values 0.36, 0.37 and 0.38 with a weight of 1/3 each). In this first model, initial intrinsic permeability is also heterogeneous: a random initial intrinsic permeability with 𝑘𝑘0 taking values 2.15 × 10⁻¹⁸ m², 3.90 × 10⁻¹⁸ m² and 7.10 × 10⁻¹⁸ m² with a weight of 1/3 each is considered. This model assumes the same relative permeability for both matrix and fractures. Hence, gas permeability reads

where Eq. (22) has been used.

Model UPC/Andra-HM2

An enhanced version of the previous hydro-mechanical model is also proposed by UPC/Andra. In this second hydro-mechanical model (the so-called UPC/Andra-HM2), different relative permeabilities for the matrix and fractures are assumed. Hence, gas permeability reads

In this second model, a higher fixed initial porosity (𝜙𝜙0 = 0.41) is assumed whereas initial intrinsic permeability is one order of magnitude smaller than in model UPC/Andra-HM1: indeed, 𝑘𝑘0 is also heterogeneous and randomly generated with values 2.15 × 10⁻¹⁹ m² (weight 1/6), 3.90 × 10⁻¹⁹ m² (weight 1/6) and 7.10 × 10⁻¹⁹ m² (weight 2/3).

Discrete approaches

Model LBNL-D

LBNL proposed a discrete technique where a two-phase flow model is coupled to a discrete fracture network (DFN). In particular, the rigid-body-spring network (RBSN) approach, which can be categorized as a lattice model and is based on the rigid-body-spring model developed by Kawai (1978), is assumed to characterize the mechanical and fracture-damage behavior. Thus, the fracture process of a local rigid-body-spring element is realized by degrading the springs. Hence, the stiffness matrix reads

where, similarly as the proposed damage strategies (see Eq. (11)), 𝐷𝐷 is a scalar parameter (0 ≤ 𝐷𝐷 ≤ 1). In LBNL-D model, 𝐷𝐷 is directly switched from 0 to 1 once a fracture event occurs, i.e. once the stress state of an element violates the Mohr-Coulomb criteria. As done by UPC/Andra, see Eq. (22), the permeability is porosity-dependent and reads

where again, 𝜙𝜙0 [-] is the initial porosity, 𝑘𝑘0 [m²] is the reference permeability, 𝜙𝜙 [-] is the porosity,

𝑏𝑏 [m] is the fracture aperture and 𝑎𝑎 [m] is the element width. In this model, as done by CNSC with model CNSC-D, a swelling effect is considered. Here, swelling induced volume change d𝜀𝜀vs is related to the change in liquid saturation dS by the equation

where 𝛼𝛼𝑠𝑠 [-] is the swelling coefficient which is calibrated to match the peak values of total stress responses (𝛼𝛼𝑠𝑠 = 0.1).

Model SNL

Another conceptual model, where special emphasis is placed on the capture of dilatancy was proposed by SNL. This is based on the concept of the delay logistic model (Strogatz, 2001; Bani- Yaghoub, 2017). The model underlying assumption is that given the low permeability of the material, the dominant mechanism for gas migration is first to form a bubble nucleation and then to push the bubble through the clay matrix through matrix dilation and fracturing. Thus, the evolution of mass and pressure within a bubble of a volume V is simply expressed by

where 𝑑𝑑 is the gas mass in the bubble; 𝑃𝑃 is the gas pressure in the bubble; 𝑃𝑃𝑡𝑡 and 𝑃𝑃𝑑𝑑 are the gas pressures in the upstream and the downstream of the bubble movement respectively; 𝑘𝑘𝑡𝑡 and 𝑘𝑘𝑑𝑑 are the permeability of the matrix in the upstream and the downstream of the bubble movement respectively; and t is the time.

For the sake of simplicity, 𝑘𝑘𝑡𝑡 and 𝑘𝑘𝑑𝑑 are assumed to be proportional to the gas pressure 𝑃𝑃. Thus,

where 𝑘𝑘0 and 𝑘𝑘0 are constant. Then, and assuming the ideal gas law, Eq. (28) becomes the continuous logistic equation with

In order to account for the clay ‘memory’ effect, the permeabilities 𝑘𝑘𝑡𝑡 and 𝑘𝑘𝑑𝑑 are assumed to depend not only on the current pressure value (see Eq. (29) and (30)) but on the pressure history. Thus, Eq. (31) becomes

where 𝐺𝐺(𝑑𝑑) is a kernel function.

Due to the preliminary nature of this novel technique, this is not taken into account in the comparison analysis in the Results section.