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:
- 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).
- 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).
- 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).
|BGR/UFZ||OpenGeoSys||Kolditz et al., 2012|
|CNSC-E||COMSOL Multiphysics ®||COMSOL, 2017|
|KAERI||TOUGH-MP FLAC3D||Zhang et al., 2008, Itasca 2018|
|LBNL-D||TOUGH-RBSN||Asahina et al., 2014|
|NCU/TPC||THMC 7.1||Yeh et al., 2013|
|UPC/Andra-H||CODE_BRIGHT||Olivella et al., 1996|
|BGR/UFZ||2D axisymmetric||Triangle|| 2715 elements|
|KAERI||3D||Hexahedral|| 7560 elements|
|LBNL-C||3D||Hexahedral|| 23400 elements|
|LBNL-D||2D||Cells and lattice|| 1401 cells|
3840 lattice elements
|NCU/TPC||2D||Quadrilateral|| 1250 elements|
|Quintessa/RWM||1D finite volume||-||20 cells +1 abstract cell to simulate mass balance in the injector|
|UPC/Andra-H||3D||Hexahedral|| 7168 elements|
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:
- Young's modulus: 𝐸𝐸 = 307 MPa
- Poisson's ratio: 𝜈𝜈 = 0.4
- Porosity: 𝜙𝜙 = 0.44
- 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.
Basic fixed parameters
|E [MPa]||𝝂𝝂 [-]||𝝓𝝓 [-]||k [m2]|
|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|
(𝜙𝜙0 = 0.44)
|3.4 × 10−21|
|Quintessa/RWM||307||0.4||-||3.4 × 10−21|
(𝜙𝜙0 = 0.36–0.38)
(𝜙𝜙0 = 0.36–0.38)
|BGR/UFZ|| * Intrinsic permeability: 1 × 10−21 m2|
* Critical pressure
|CNSC-D|| * Air-entry value|
* Minimum air-entry value
* Maximum intrinsic permeability
* Damage smoothing parameter
|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.
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
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.
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
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
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).
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 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)
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.
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.
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 [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.
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).
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 [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).
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.
- 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.
- COMSOL. 2017. COMSOL Multiphysics version 5.3a. Manual.
- 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.
- Itasca, FLAC3D Version 6.0. 2018. Fast Lagrangian Analysis of Continua in 3 Dimensions. ITASCA Consulting Group Inc. https://www.itascacg.com/software/flac3d .
- 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.
- 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.
- 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.
- Quintessa. 2013. QPAC: Quintessa’s General-Purpose Modelling Software QRS-QPAC-11. http://www.quintessa.org/qpac-overview-report.pdf.
- 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.
- 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.
- 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.
- 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.1918.104.22.1681.
- 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.
- 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.
- 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.
- 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.
- 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.
- Strogatz, S H. 2001. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering, Westview.
- Bani-Yaghoub, M. 2017. Analysis and applications of delay differential equations in biology and medicine, arXiv:1701.04173v1.