Difference between revisions of "OR/18/049 Modelling approaches"

From Earthwise
Jump to navigation Jump to search
[checked revision][checked revision]
m (1 revision imported)
Line 345: Line 345:
===Model KAERI===
===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)<ref name="Fall 2014">
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)<ref name="Fall 2014"></ref>. That is, damage variable under tension is described as
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.
</ref>. That is, damage variable under tension is described as
[[Image:OR18049equation15.jpg|frameless|center|225px|      ]]
[[Image:OR18049equation15.jpg|frameless|center|225px|      ]]

Latest revision as of 15:23, 3 December 2019

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[edit]

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[edit]

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[edit]

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[edit]

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[edit]

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

Model KAERI[edit]

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

Model BGR/UFZ[edit]

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


Model NCU/TPC[edit]

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[edit]

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[edit]

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[edit]

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[edit]

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[edit]

Model LBNL-D[edit]

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[edit]

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.