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).
Model | Software | Reference |
BGR/UFZ | OpenGeoSys | Kolditz et al., 2012^{[1]} |
CNSC-E | COMSOL Multiphysics ยฎ | COMSOL, 2017^{[2]} |
CNSC-D | ||
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]} |
UPC/Andra-HM1 | ||
UPC/Andra-HM2 |
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:
- Young's modulus: ๐ธ๐ธ = 307 MPa
- Poisson's ratio: ๐๐ = 0.4
- Porosity: ๐๐ = 0.44
- Intrinsic permeability: ๐๐ = 3.4 ร 10^{โ21} m^{2}
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.
Model | Basic fixed parameters | |||
E [MPa] | ๐๐ [-] | ๐๐ [-] | k [m^{2}] | |
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 |
Model | Calibrated parameters |
BGR/UFZ | * Intrinsic permeability: 1 ร 10^{โ21} m^{2} * 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} m^{2} * 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} m^{2} |
UPC/Andra-HM1 | * Intrinsic permeability: Het. * ๐๐0 = (2.15โ7.10) ร 10^{โ18} m^{2} * Embedded permeability parameters |
UPC/Andra-HM2 | * Intrinsic permeability: Het. * ๐๐0 = (2.15โ7.10) ร 10^{โ19} m^{2} * 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 ๐๐ [m^{2}] 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^{2} 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
where
with ๐๐max [m^{2}] being a calibrated maximum intrinsic permeability (๐๐max = 1.00 ร 10โ18 m^{2}). 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
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} m^{2}).
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^{[15]})
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^{2}] 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^{2}]
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^{2}]
where ๐๐ [m^{2}] is the calibrated capillary spacing (๐๐ = 5.66 ร 10^{-6} m^{2}). The capillary radius ๐๐ [m] is considered to be dependent on the capillarity compressibility ๐พ๐พ [m2Pa-1] (assumed to be ๐พ๐พ = 1 ร 10^{-20} m^{2}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^{[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 [m^{2}] 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} m^{2}, 3.90 ร 10^{โ18} m^{2} and 7.10 ร 10^{โ18} m^{2} 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} m^{2} (weight 1/6), 3.90 ร 10^{โ19} m^{2} (weight 1/6) and 7.10 ร 10^{โ19} m^{2} (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 [m^{2}] 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.
References
- โ 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.1998.48.5.681.
- โ ^{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.
- โ 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.