OR/18/049 Modelling approaches

From Earthwise
Jump to: navigation, search
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).
  1. 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).
  1. 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[1]
CNSC-E COMSOL Multiphysics ยฎ COMSOL, 2017[2]
KAERI TOUGH-MP FLAC3D Zhang et al., 2008[3], Itasca 2018[4]
LBNL-C TOUGH-FLAC Rutqvist, 2011[5]
LBNL-D TOUGH-RBSN Asahina et al., 2014[6]
NCU/TPC THMC 7.1 Yeh et al., 2013[7]
Quintessa/RWM QPAC Quintessa, 2013[8]
UPC/Andra-H CODE_BRIGHT Olivella et al., 1996[9]
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

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โˆ’21 m2

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.

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)
UPC/Andra-HM1 225 0.125 Heterogeneous
(๐œ™๐œ™0 = 0.36โ€“0.38)
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)[10] 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 ๐‘˜๐‘˜ [m2] 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 m2 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)[11] 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)[12] 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[13]), 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




with ๐‘˜๐‘˜max [m2] being a calibrated maximum intrinsic permeability (๐‘˜๐‘˜max = 1.00 ร— 10โ€“18 m2). 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)[14].


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)[13]. 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 m2).


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[15])



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 [m2] 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 ๐‘˜๐‘˜๐‘ค๐‘ค [m2]


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 ๐‘˜๐‘˜๐‘”๐‘” [m2]


where ๐‘Ž๐‘Ž [m2] is the calibrated capillary spacing (๐‘Ž๐‘Ž = 5.66 ร— 10-6 m2). The capillary radius ๐‘Ÿ๐‘Ÿ [m] is considered to be dependent on the capillarity compressibility ๐›พ๐›พ [m2Pa-1] (assumed to be ๐›พ๐›พ = 1 ร— 10-20 m2Pa-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[16]. 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 [m2] 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โˆ’18 m2, 3.90 ร— 10โˆ’18 m2 and 7.10 ร— 10โˆ’18 m2 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โˆ’19 m2 (weight 1/6), 3.90 ร— 10โˆ’19 m2 (weight 1/6) and 7.10 ร— 10โˆ’19 m2 (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)[17], 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 [m2] 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[18]; Bani- Yaghoub, 2017[19]). 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.


  1. โ†‘ Kolditz, O, Bauer, S, Bilke, L, Bรถttcher, N, Delfs, J O, Fischer, T, Gรถrke, U J, Kalbacher, T, Kosakowski, G, McDermott, C I, Park, C H, Radu, F, Rink, K, Shao, H, Shao, H B, Sun, F, Sun, Y Y, Singh, A K, Taron, J, Walther, M, Wang, W, Watanabe, N, Wu, N, Xie, M, Xu, W and Zehner, B. 2012. OpenGeoSys: An open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (THM/C) processes in porous media. Environmental Earth Sciences. 67, 589โ€“599. doi: 10.1007/s12665-012-1546-x.
  2. โ†‘ COMSOL. 2017. COMSOL Multiphysics version 5.3a. Manual.
  3. โ†‘ Zhang, K, Wu, Y S, and Pruess, K. 2008. User's guide for TOUGH2-MP-A massively parallel version of the TOUGH2 code. Earth Sciences Division, Lawrence Berkeley National Laboratory. https://escholarship.org/uc/item/00d9040f.
  4. โ†‘ Itasca, FLAC3D Version 6.0. 2018. Fast Lagrangian Analysis of Continua in 3 Dimensions. ITASCA Consulting Group Inc. https://www.itascacg.com/software/flac3d .
  5. โ†‘ Rutqvist, J. 2011. Status of the TOUGH-FLAC simulator and recent applications related to coupled fluid flow and crustal deformations. Computers and Geosciences. 37, 739โ€“750. doi: 10.1016/j.cageo.2010.08.006.
  6. โ†‘ Asahina, D, Houseworth, J E, Birkholzer, J T, Rutqvist, J, and Bolander, J E. 2014. Hydro-mechanical model for wetting/drying and fracture development in geomaterials. Computers and Geosciences. 65, 13โ€“23. doi: 10.1016/j.cageo.2013.12.009.
  7. โ†‘ Yeh, G T, Tsai, C H, and Liu, I S. 2013. GMech: A Geo-Mechanics Model for Finite Visco-Elastic Materials: Theoretcal Basis and Numerical Approximation. Graduate Institute of Applied Geology, National Central University, Jhongli City, Taoyuan County, Taiwan.
  8. โ†‘ Quintessa. 2013. QPAC: Quintessaโ€™s General-Purpose Modelling Software QRS-QPAC-11. http://www.quintessa.org/qpac-overview-report.pdf.
  9. โ†‘ Olivella, S, Gens, A, Carrera, J, and Alonso, E E. 1996. Numerical formulation for a simulator (CODE_BRIGHT) for the coupled analysis of saline media. Engineering Computations. 13, 87โ€“112. doi: 10.1108/02644409610151575.
  10. โ†‘ van Genuchten, M T. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil science society of America journal. 44, 892โ€“898. doi: 10.2136/sssaj1980.03615995004400050002x.
  11. โ†‘ Bandis, S C, Lumsden, A C, and Barton, N R. 1983. Fundamentals of rock joint deformation. International Journal of Rock Mechanics and Mining Sciences and Geomechanic Abstracts. 20, 249โ€“68. doi: 10.1016/0148-9062(83)90595-8.
  12. โ†‘ Khalili, N, and Khabbaz, M H. 1998. A unique relationship for ฯ‡ for the determination of the shear strength of unsaturated soils. Gรฉotechnique. 48, 1โ€“7. doi: 10.1680/geot.1998.48.5.681.
  13. โ†‘ 13.0 13.1 Fall, M, Nasir, O and Nguyen, T S. 2014. A coupled hydro-mechanical model for simulation of gas migration in host sedimentary rocks for nuclear waste repositories. Engineering Geology. 176, 24โ€“44. doi: 10.1016/j.enggeo.2014.04.003.
  14. โ†‘ Nasir, O, Nguyen, T S, Barnichon, J D, and Millard, A. 2016. Simulation of hydromechanical behaviour of bentonite seals for containment of radioactive wastes. Canadian Geotechnical Journal. 54, 1055โ€“1070. doi: 10.1139/cgj-2016-0102.
  15. โ†‘ Xu, W J, Shao, H, Hesser, J, Wang, W, Schuster, K, and Kolditz, O. 2013. Coupled multiphase flow and elasto-plastic modelling of in-situ gas injection experiments in saturated claystone (Mont Terri Rock Laboratory). Engineering Geology. 157, 55โ€“68. doi: 10.1016/j.enggeo.2013.02.005.
  16. โ†‘ Olivella, S, and Alonso, E E. 2008. Gas flow through clay barriers. Gรฉotechnique. 58, 157โ€“176. doi: 10.1680/geot.2008.58.3.157.
  17. โ†‘ Kawai, T. 1978. New discrete models and their application to seismic response analysis of structures. Nuclear Engineering and Design. 48, 207โ€“229. doi: 10.1016/0029-5493(78)90217-0.
  18. โ†‘ Strogatz, S H. 2001. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering, Westview.
  19. โ†‘ Bani-Yaghoub, M. 2017. Analysis and applications of delay differential equations in biology and medicine, arXiv:1701.04173v1.