Ngoc-Tuyen Tran

Duc-Phi Do

Dashnor Hoxha

Minh-Ngoc Vu

Gilles Armand

Kriging-based reliability analysis of the long-term stability of a deep drift constructed in the Callovo-Oxfordian claystone

Keywords:

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The French National Radioactive Waste Management Agency (Andra) is in charge of studying the disposal of high-level and intermediate-level long-lived waste (HLW and ILW-LL) in deep geological disposal (Cigéo project). Cigéo will be located at the border of the Meuse and Haute-Marne departments, nearly 300 km east of Paris. The host formation consists of claystone (Callovo-Oxfordian argillaceous rock -COx) lying between 420 m and 550 m in depth and exhibiting very favorable conditions for a repository of radioactive waste, as it generally has a very low hydraulic conductivity, small molecular diffusion, and significant retention capacity for radionuclide. The underground disposal facility will be divided into two sections depending on the type of waste (ILW-LL or HLW) and the type of excavation used to emplace the waste disposal packages (tunnel for ILW-LL and steel cased micro-tunnel for HLW). The ILW-LL zone is constituted of alveoli/galleries of about 10 m in diameter, which is planned to be excavated by tunnel boring machine (TBM) with an immediate installation of the liner. The ILW-LL cell's liner is the concrete arch with a thickness of 0.5 m covered by 0.2 m of compressible material. The installation of the compressible material in the support system is currently regarded as an innovative and promising construction of the underground structure within the time-dependent geological formation thanks to its high compressibility (deformation more than 50% under stress lower than 3 MPa) to absorb the high convergence of the host rock.

This study aims at investigating the stability of an ILW-LL drift of Cigéo during the exploitation period by accounting for the uncertainty of the host rock behavior. The maximal stress evolution *Corresponding author. Email-address: duc-phi.do@univ-orleans.fr of the concrete liner is considered as the time-dependent indicator of the ILW-LL alveolus' stability.

Actually, the time-dependent response of the ILW-LL cell is governed by the thermo-hydro-mechanical (THM) coupling. Indeed, the heat released from ILW-LL packages emplaced within the alveolus makes the temperature rise in the host rock and the structural elements (≤ 65 °C in the concrete structure), which can induce the thermo-mechanical stress and the thermal pressurization phenomenon [START_REF] Conil | In situ investigation of the THM behavior of the Callovo-Oxfordian claystone[END_REF][START_REF] Vu | Thermal pressurization coefficient of anisotropic elastic porous media[END_REF]. Besides, the timedependent behavior of ILW-LL drift is also controlled by the drainage from the intact rock to the wall, the desaturation of the rock close to the wall due to the ventilation during the exploitation phase, and shrinkage/swelling of COx claystone due to the desaturation/resaturation [START_REF] Wang | Numerical verification of selfsealing of fractured COx during resaturation[END_REF]. However, it has been shown in different contributions [START_REF] Camusso | CIGEO radioactive waste repository project-An observation-based model of claystone behavior for thermomechanical FLAC3D simulations[END_REF][START_REF] Alonso | Effect of thermohydro-mechanical coupling on the evolution of stress in the concrete liner of un underground drift in the Cigéo project[END_REF] that the creep behavior of host rock has a preponderant role in the stress state of the concrete liner, while the effects of the heat loading and the pore water pressure change are marginal.

Therefore, this study focuses only on the uncertainty effect of the creep behavior of the host rock on the exceedance probability of critical stress in the concrete liner of the ILW-LL drift. Both laboratory tests at the sample scale and the convergence measurement in the Andra Underground Research Laboratory (URL) show that the time-dependent response of COx claystone exhibits a significant initial strain rate that decreases with time [START_REF] Armand | Short-and long-term behaviors of drifts in the Callovo-Oxfordian claystone at the Meuse/Haute-Marne Underground Research Laboratory[END_REF][START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of Callovo-Oxfordian claystone: From experimental studies to model calibration and validation[END_REF]. Advances in numerical models have been proposed to reproduce the observation of short-and long-term behaviors of COx rock, e.g. non-local anisotropic elasto-viscoplasticity (Manica et al., 2021a, b), second gradient anisotropic elastoplasticity [START_REF] Pardoen | Modelling the influence of strain localisation and viscosity on the behaviour of underground drifts drilled in claystone[END_REF], double phase-field elastoplastic damage model [START_REF] Yu | Numerical study of thermo-hydromechanical responses of in situ heating test with phase-field model[END_REF], and elastoplastic model combined with weakness planes [START_REF] Souley | 3D anisotropic modelling of deep drifts at the Meuse/Haute-Marne URL[END_REF]. Characterized by various parameters (e.g. larger than 20 parameters), whose physical meaning and calibration are not simple to be evaluated, these sophisticated models are not feasible for the reliability analysis. In this work, the use of the well-known Lemaitre's model is motivated, on one hand, by its simplicity (characterized by only three parameters easily calibrated from creep tests) and, on the other hand, by its accuracy to describe the longterm behavior of host rock (Souley et al., 2017a). Available laboratory triaxial creep tests under different deviatoric stresses are taken to estimate the uncertainty of viscoplastic properties of COx claystone.

The stability of the concrete support of underground structure during the period of exploitation (~100 years) is crucial to ensure the functionality of the mechanical system within the tunnel (i.e. the movement and deposit of waste packages). In this context, a delay effect due to the time-dependent behavior of host rock on the stability of structures is critical and must be taken into account in the optimization design of these constructions. The widely used deterministic (or semi-probabilistic) design methods in geotechnics based on several partial safety factors for considering different uncertainty sources (e.g. geomaterial property uncertainty, loading uncertainties, etc.) have pointed out several limits in rock engineering [START_REF] Spross | On the need for a risk-based framework in Eurocode 7 to facilitate design of underground openings in rock[END_REF]. The difficulty lies in the uncertainty quantification of the delay effect and its propagation on structural stability. This is mainly because the experimental tests in the laboratory and/or in situ can only be conducted within a short period as compared with the service life of almost all structures, particularly with those related to the deep geological nuclear waste disposal.

In the last decades, the reliability analysis has been received much attention, and its progress toward practical applications in geotechnical rock engineering has been largely remarked [START_REF] Wang | LS-SVM and Monte Carlo methods based reliability analysis for settlement of soft clayey foundation[END_REF][START_REF] Contreras | Slope reliability and back analysis of failure with geotechnical parameters estimated using Bayesian inference[END_REF][START_REF] Siacara | Reliability analysis of earth dams using direct coupling[END_REF]. Especially, in many contributions [START_REF] Laso | A level II reliability approach to tunnel support design[END_REF][START_REF] Schweiger | Reliability analysis in geotechnics with deterministic finite elements -Theoretical concepts and practical application[END_REF][START_REF] Mollon | Probabilistic analysis and design of circular tunnels against face stability[END_REF][START_REF] Langford | Reliability based approach to tunnel lining design using a modified point estimate method[END_REF][START_REF] Lü | System reliability assessment for a rock tunnel with multiple failure modes[END_REF], the probabilistic-based method has successfully been applied to the stability analysis and the optimization design of tunnel support in the rock mass. Following that, the Monte Carlo simulation (MCS), the first-or second-order reliability method (FORM/SORM), and response surface method (RSM) are among the most chosen probabilistic methods. However, these contributions dealt almost exclusively with the short-term behavior of underground excavation by accounting for the uncertainty of the elastoplastic parameters of the host rock. Recently, Do et al. (2020) evaluated the failure probability in time of a deep tunnel excavated in a viscoelastic rock by using the MCS. Their study revealed a critical effect of the uncertainty of the viscoelastic parameters on the long-term stability of the tunnel. Ignoring the uncertainty of the time-dependent characteristic of the rock mass can thus lead to an underestimation of the long-term structural stability of underground structures.

Due to the time-demanding computation of the deterministic problem as well as the non-linear behavior of the viscoplastic COx rock, the Kriging metamodeling technique is chosen in this study for the reliability analysis of the ILW-LL alveolus. An extension of the well-known AK-MCS (active learning reliability method combining Kriging and MCS) method is proposed. In comparison with the massive number of trials in the MCS method, this Kriging metamodel allows a fast computation by using a very reasonable number of direct numerical simulation to solve the deterministic problem. This method provides hence a useful and performant tool in the reliability analysis and optimization design, especially when the response of complex structures is only assessable by the numerical simulation, such as the present case of deep drift in the viscoplastic rock.

The paper is organized as follows. In the next section, the geometry, loading (in situ stress), and the constitutive mechanical behavior of all components (host rock, concrete liner, and compressible material) as well as the uncertainty quantification of the long-term mechanical properties of COx claystone are presented. In Section 3, the Kriging metamodeling technique based on the modified AK-MCS method is described. In Section 4, numerical investigations on both the deterministic and the reliability analysis problems are provided. In addition to the uncertainty effect of COx rock, sensitivity analysis on the thickness of both liners (concrete and compressible material) and the compressibility of the compressible layer is also carried out in the reliability analysis. Discussion and some critical remarks derived from these numerical results are then pointed out. Conclusions and perspectives are drawn in the last section. Notably, throughout the paper, the light symbols denote scalars, and the bold images indicate the second-order tensors (as a vector or a matrix).

Description of the studied problem

Geometry and loading

At the position of the Meuse/Haute-Marne (M/HM) URL, the initial stress state is characterized by three principal stresses: vertical stress σv corresponding to the ground weight, minor horizontal stress very close to the vertical stress (σv ≈ σh ≈ 12.5 MPa), and major horizontal stress (σH ≈ 1.3σh) [START_REF] Armand | Short-and long-term behaviors of drifts in the Callovo-Oxfordian claystone at the Meuse/Haute-Marne Underground Research Laboratory[END_REF]. According to the design of the Cigéo project, the alveoli dedicating to store the nuclear waste (ILW-LL and HLW) will be planned to excavate following the major horizontal stress (i.e. the stress in the tunnel section is quasiisotropic). This choice is undertaken due to the long-term safety of nuclear waste management. The length of an ILLW-LL alveolus (about 500 m) is much greater than its diameter (about 10 m); thus, a plane strain model is adopted in this study.

It is assumed that the parallel drifts are separated enough to ignore their mechanical interaction, i.e. this study limits to an individual drift behavior. Fig. 1 illustrates one-quarter of the adopted twodimensional (2D) plane strain model of drift, considering the symmetric conditions. Following that, to simulate the mechanical behavior of the circular drift with 5.05 m of radius, the width and the height of 55 m are taken. Concerning the symmetric boundaries, the normal displacement is fixed. At the top and the right interfaces, normal stress σ0 = 12.5 MPa, which is approximated the measured vertical (σv) and minor horizontal (σh) stresses at the M/HM URL, is imposed. It is worth noting that our proposed methodology can be applied in the general case of anisotropic initial stress state without any additional complexity.

The considered drift is supported by a concrete lining (the class of concrete is C60/75 according to the European Code classification), which is separated from the rock mass by a compressible layer. The drift excavation is assumed to be instantaneous whilst the simultaneous installation of the compressible liner and the final concrete support corresponds to a deconfinement ratio of 85%. The role of the compressible support is to absorb the convergence and to reduce the transmitted stress to the concrete lining which is supposed to be linear elastic.

Time-dependent behavior of COx claystone

In this work, Lemaitre's model is considered to express the timedependent behavior of host rock. Most of the authors dealing with the time-dependent behavior of COx claystone have proposed various (elasto-) viscoplastic models [START_REF] Seyedi | Transverse Action" -A model benchmark exercise for numerical analysis of the Callovo-Oxfordian claystone hydromechanical response to excavation operations[END_REF][START_REF] Stavropoulou | Time-dependent behaviour of the Callovo-Oxfordian claystone-concrete interface[END_REF]. Even if other sophisticated models could be considered for specific features of COx rock behavior (e.g. [START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of Callovo-Oxfordian claystone: From experimental studies to model calibration and validation[END_REF][START_REF] Mánica | A time-dependent anisotropic model for argillaceous rocks. Application to an underground excavation in Callovo-Oxfordian claystone[END_REF][START_REF] Seyedi | Transverse Action" -A model benchmark exercise for numerical analysis of the Callovo-Oxfordian claystone hydromechanical response to excavation operations[END_REF][START_REF] Souley | 3D anisotropic modelling of deep drifts at the Meuse/Haute-Marne URL[END_REF][START_REF] Stavropoulou | Time-dependent behaviour of the Callovo-Oxfordian claystone-concrete interface[END_REF], the viscoplastic model of Lemaitre is sufficiently accurate for long-term design purposes and with a limited number of constitutive parameters.

By neglecting the (thermo-) hydro-mechanical coupling effect, the total strain tensor t ot ε reads

t o t e v = + ε ε ε (1)
where εe and εv indicate the elastic and viscous strain tensors, respectively. The viscous strain rate v

( ) &

ε is expressed as a function of deviatoric stress tensor thanks to Lemaitre's law:

eq v 1/ eq 3 1 2 n m K σ σ λ   =     % & σ ε (2) where eq v v 1/ 2 1 : 3 σ λ λ   = =     & & & n m K ε ε (3) 1 tr( ) 3 = - % I σ σ σ (4)
where σ and % σ are the stress tensor and its corresponding deviatoric part, respectively; and I indicates the second-order identity tensor. The over-dot denotes the time derivative while tr( )

σ
designates the trace of the stress tensor. Also, λ and σeq mean the accumulated viscous strain (i.e. the time-dependent distortion strain) and the equivalent (or the Von-Mises) stress, respectively. The positive triple (K, n, m) is the set of mechanical parameters characterizing the long-term behavior of the material following the Lemaitre model.

Posing α = n/m and after some developments, the accumulated viscous strain can be explicitly expressed as a function of equivalent stress and time:

1/(1 ) /(1 ) 1/(1 ) eq n 1 n t K α α α α λ σ + + + +   =     (5)
Equivalently, the more compact form is written as

e q B C A t λ σ = (6)
where

1 1 n 1 1 , , 1 1 n A B C K α α α α + +   = = =   + +   (7)
Besides, the positive triple is then inferred based on the parameters A, B, and C as follows:

1 , , 1 C B B B B n m K A C C C - - = = = - (8)

Tri-linear elastic model of the compressible liner

The compressible material used as the outer liner to reduce the overstress of the inner lining of drift has already been widely used in various research projects [START_REF] Lombardi | Bau von Tunneln bei grossen Verformungen des Gebirges[END_REF][START_REF] Strohhäusl | TBM tunnelling under high overburden with yielding segmental linings; Eureka Project EU[END_REF][START_REF] Schneider | Compex support system -Komprimierbarer Ringspaltmörtel zur Bewältigung druckhafter Gebirgsbereiche bei TBM-Vortrieben mit starrem Tübbingausbau[END_REF][START_REF] Billig | DeCo Grout -Innovative grout to cope with rock deformations in TBM tunnelling[END_REF][START_REF] Ly | Contribution à l'étude des soutènements compressibles pour les ouvrages souterrains: Validation du modèle rhéologique SC2D pour l'argilite du COx et construction d'un modèle rhéologique Coqcéram pour le matériau compressible[END_REF][START_REF] Gasbarrone | Use of expanded clay as annular gap filling[END_REF]. At the M/HM URL, different types of compressible materials have been tested, and the primary feedback from the experimental program reveals their immense benefits (e.g. reducing the transmitted stress to the concrete liner of drift). Owing to a high void ratio, the mechanical behavior of the compressible materials usually presents a non-linear characteristic due to the pore collapse mechanism under compressive loading. For example, Fig. 2a shows the results of a compression test performed on a specific compressible grout. The material is constituted by expanded polystyrene beads and foam to achieve compressibility up to more than 50% of its initial volume [START_REF] Billig | DeCo Grout -Innovative grout to cope with rock deformations in TBM tunnelling[END_REF]. The stress-strain curve from these experiments can be captured in three stages. Following that, the curve shows an elastic behavior when the volumetric deformation of the compressible material is small. Then a plastic behavior represented by a strain hardening can be observed, while the last stage characterizes the stiffening of the material after the depletion of its porosity.

Many attempts to reproduce the mechanical behavior of the high porous compressible materials have been conducted. For example, in Souley et al. (2017b), an elastoplastic constitutive model was proposed. In this model, the Drucker-Prager criterion and a strain hardening based on an exponential function of the plastic distortion have been applied. Besides, the yield function, the hardening, and the densification mechanism are explicitly related to the pore collapse mechanism (e.g. volumetric plastic strain due to the hydrostatic failure).

For reasons of simplification, a tri-linear elastic behavior was adopted here for the compressible support layer (Fig. 2b). This simplification is acceptable when in our analyses of drift, the compressive stress state in this compressible liner increases monotonically with time. It is worth mentioning that, under the monotonic loading condition, the stress-strain curve provided by our tri-linear elastic model can well fit the experimental data (Fig. 2a).

As a function of the volumetric deformation, the tri-linear elastic model is characterized by the triple values of elastic modulus (Ec1, Ec2, and Ec3) and is expressed in the following form: where εv1 and εv2 present the first and second elastic limits of the trilinear elastic model, respectively. Furthermore, εV and ed are the volumetric deformation (strain) and deviatoric strain tensor, respectively, while the bulk and shear moduli are calculated from the corresponding Young's modulus and Poisson's ratio of each stage:

c1 V c1 d V v1 c1 c2 v1 c2 V c2 d v1 V v 2 c1 c2 v1 c2 c3 v 2 c3 V c3 d v 2 V 2 ( ) ( ) 2 ( ) ( ) ( ) 2 ( ) K G K K K G K K K K K G ε ε ε ε ε ε ε ε ε ε ε ε ε + ≤   = - + < ≤   - + - + < 
c c c c c c / (1 2 ) ( 1, 2, 3) / [2(1 )] i i i i K E i G E ν ν = - =    = +   (10)
The Poisson's ratio of the compressible material is so small and presents a neglected impact on the drift response. Thus, it can be taken equal to zero (νc = 0). For the numerical simulations, this trilinear elastic model has been implemented in the Code_Aster (EDF, 2017).

Choosing the failure criterion

The reliability analysis aims at estimating the failure probability of structure that its responses exceed one or multiple failure criteria by accounting for the propagation of uncertainty of different input parameters. A failure criterion can be mathematically expressed by a limit-state function (LSF) g(x) (also called the performance function). The structure fails when the value of LSF is less than zero (g(x) ≤ 0). In the recent contribution of Do et al. (2020), two time-dependent performance functions were proposed to define the failure state of a deep tunnel excavated in a rheological rock. Following that, the failure of the tunnel can be investigated for the allowable convergence and/or the allowable stress in the system of liners.

Since there is the presence of the outer compressible liner, the high convergence on the tunnel wall can be accepted, and our work focuses mainly on the prediction of failure probability of the concrete liner of drift. We consider that the concrete liner fails when the maximum equivalent stress l ( ) q x overpasses the allowable stress cl σ of the constitutive material. The corresponding performance function is then defined as follows:

cl l ( ) ( ) g q σ = - x x ( 11
)
where the vector x denotes the set of random input variables in the reliability analysis. In comparison with the artificial constitutive materials of the compressible and concrete liners, the uncertainty of the mechanical properties of the natural COx rock is much higher. Consequently, only the uncertainty of the host rock properties is counted in this work.

According to the discussion in Do et al. (2020), for the case that the tunnel response over time is monotonic, the time-dependent reliability analysis can be done similarly to the time-invariant problem by calculating only the probability of failure at a chosen design lifetime of the tunnel. The failure probability of the concrete lining of drift is measured at a period of 100 years, which corresponds to the period of exploitation of the deep geological repository, and hence in Eq. ( 11), the time parameter t is omitted.

It is important to point out that the non-linear behavior of the viscoplastic COx rock is the reason that the performance function in Eq. ( 11) cannot be expressed explicitly in the spaces of random variables but only be numerically evaluated at each selected subset of x.

Uncertainty quantification of mechanical properties of COx rock

In general, the sources of uncertainty of geomaterials can be classified into two groups. The first one, namely the epistemic uncertainty, is associated with the scarcity of data (quality, quantity, specification), and it comes from different stages of sampling preparation, transport, and testing procedure. The second one, called aleatory uncertainty, refers to the inherent spatial variability of rock. Regarding the COx claystone, the inherent variability of the mineralogical composition of this host rock, in addition to the scarcity of high-quality data, can contribute to a significant uncertainty of its mechanical properties.

In the present work, the uncertainty of COx rock properties is quantified from data measured in the laboratory on natural samples obtained from the main level of M/HM URL. According to the contribution of [START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of Callovo-Oxfordian claystone: From experimental studies to model calibration and validation[END_REF], we assume that the mean value and standard deviation are 6 GPa and 3.5 GPa for Young's modulus, respectively. Note that the influence of the Poisson's ratio of rock mass on the behavior of drift is quite moderate as compared with the impact of the other parameters. Thus, its mean value (ν = 0.29) is used in this work as a deterministic parameter by ignoring its uncertainty.

Regarding the uncertainty quantification of the long-term mechanical properties of COx formation, this topic has not yet been discussed in the literature. So far, a universal constitutive model to accurately characterize the time-dependent behavior of COx rock does not exist. Additionally, the number of creep tests carried out on the samples of the same level and with the same condition seems very limited from the statistical point of view.

As the first attempt to quantify the uncertainty of the long-term mechanical properties of COx claystone, in this study, we use the data of seven triaxial creep tests corresponding to seven applied deviatoric stresses performed on the samples taken from the same geological horizon [START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of Callovo-Oxfordian claystone: From experimental studies to model calibration and validation[END_REF][START_REF] Seyedi | Transverse Action" -A model benchmark exercise for numerical analysis of the Callovo-Oxfordian claystone hydromechanical response to excavation operations[END_REF]. The information about the confining pressure and deviatoric stress of each triaxial creep test is summarized in Table 1, while for more details of these tests, the interested readers can refer to [START_REF] Armand | Fundamental aspects of the hydromechanical behaviour of Callovo-Oxfordian claystone: From experimental studies to model calibration and validation[END_REF] and different references cited therein. The calibration of the viscoplastic parameters of Lemaitre's model is conducted by using the compact form expressed in Eq. ( 6). Following this equation, the inverse analysis to determine the three parameters A, B, and C can be carried out by fitting the data of at least two creep tests under two different applied deviatoric stresses [START_REF] Debernardi | New viscoplastic model for design analysis of tunnels in squeezing conditions[END_REF]. Since two creep tests are performed with the same deviatoric stress (see Table 1), thus totally, the statistical analysis of the obtained results of three parameters A, B, and C is conducted with 88 samples generated from at least two tests. The mean and standard deviation of the viscoplastic Lemaitre's parameters are summarized in Table 2. As a function of the applied deviatoric stress and time, the exponential parameters (i.e. the parameters B and C), which characterize the non-linearity of COx rock, present the corresponding mean values of about 2.05 and 0.17, respectively. The coefficient A presents the most considerable uncertainty represented by the highest coefficient of variation (COV) (about 87.5%). Our analysis also exhibits the correlation between these three parameters (see Table 3). More precisely, a relatively significant anti-correlation between the coefficients A and B characterized by a negative value of -0.72 is stated. The other correlations (for example, between A and C or between B and C) seem weaker and are not taken into account in our reliability analysis. Following the histograms captured in Fig. 3, a log-normal distribution is adopted for these parameters. This chosen distribution avoids the negative value of Lemaitre's model parameters (notably in the case of the coefficient A with high COV) generated from the sampling process in the reliability analysis. As an example, Fig. 4 presents the results of a triaxial creep test with the applied deviatoric stress of 27.75 MPa and confining pressure of 12 MPa, which are simulated from the mean values and 1000 random sets of the three parameters A, B, and C of Lemaitre's model. Note here that the number of 1000 random sets of the Lemaitre's model parameters is chosen just for the illustration purpose of the simulated triaxial creep test. But for the reliability analysis (by using the MCS for example), much higher random samples (about one million) will be taken as detailed in the next section.

The comparison of the numerical result using the mean values of A, B, and C, and the experimental data measured up to about 90 d (Fig. 4a) illustrates the correct tendency. However, due to the uncertainty of the long-term mechanical properties, the viscous deformation in COx rock may be higher than 10 times the ones measured with the mean values (Fig. 4a andb). It is important to note here that the calibration and the corresponding statistical analysis of the long-term mechanical properties of COx rock are conducted with the creep tests performed for about 3 months. Consequently, an extrapolation of the viscous deformation to 100 years (Fig. 4b) must also present another source of uncertainty that may underestimate or overestimate the corresponding results of the reliability analysis.

Chosen deterministic parameters

Herein the deterministic parameters like the thickness (l1) and mechanical properties of the tri-linear elastic compressible material (Ec1, Ec2, Ec3, εv1) as well as the elastic moduli of concrete C60/75 (El, νl) are chosen for the numerical applications and are summarized in Table 4. Notably, the mechanical properties of the compressible material of the outer liner are calibrated from the oedometric tests of the compressible grout [START_REF] Billig | DeCo Grout -Innovative grout to cope with rock deformations in TBM tunnelling[END_REF].

The values of the other deterministic parameters, such as the thickness of the final support (l2) and the limitation (εv2) which separates the second and third stages of the tri-linear elastic model of compressible material, will be chosen in the context of a parametric study. Regarding the contrast of Young's modulus between the second stage and the other stages, the parameter εv2 represents, in fact, the compressible potential of the first support layer of drift to adsorb the high time-dependent convergence of host rock.

The methodology of Kriging-based reliability analysis

The application of the probabilistic approach to the stability analysis [START_REF] Li | Reliability analysis of circular tunnel under hydrostatic stress field[END_REF]Lü et al., 2011[START_REF] Lü | System reliability assessment for a rock tunnel with multiple failure modes[END_REF]Lü and Low, 2011;[START_REF] Langford | Reliability based approach to tunnel lining design using a modified point estimate method[END_REF] and optimization design [START_REF] Laso | A level II reliability approach to tunnel support design[END_REF][START_REF] Mollon | Probabilistic analysis and design of circular tunnels against face stability[END_REF][START_REF] Tran | Reliability-based design of deep tunnel excavated in the viscoelastic Burgers rocks[END_REF] of underground structures has been intensively undertaken since the last two decades. The most well-known method may be the MCS or crude MCS, which could provide an accurate estimation of the failure probability and has been regarded as the benchmark to validate the other probabilistic methods. Since the crude MCS requires a massive number of evaluations of structure response, it seems only to be appropriate in case that the closed-form solution is available (e.g. Do et al., 2020). Nevertheless, due to the complex behavior of rock (e.g. the case of drift constructed in the non-linear viscoplastic rock in this study), in most cases, structural behavior is only assessed through numerical calculation. An appropriate approximation method that could reduce the required number of time-consuming numerical simulations is preferred. For example, in [START_REF] Wang | Review of metamodeling techniques in support of engineering design optimization[END_REF], [START_REF] Lü | System reliability assessment for a rock tunnel with multiple failure modes[END_REF][START_REF] Tran | Reliability-based design of deep tunnel excavated in the viscoelastic Burgers rocks[END_REF], the scholars used the RSM to approximate the LSF, while the SORM was chosen to estimate the failure probability. Although the powerful RSM was demonstrated in the short-term reliability of the tunnel, this method can fail when the LSF is highly non-linear, non-smooth, or multimodal (multiple most probable points (MPPs)) [START_REF] Huang | Overview of structural reliability analysis methods -Part III : Global reliability methods[END_REF].

Many other advanced probabilistic methods have been proposed in the literature. Generally, the chosen approach aims at establishing a metamodel (or surrogate model), which approximates the implicit LSF by a mathematical function. The reliability analysis is then conducted on the built metamodel (by using, for example, the classical MCS sampling method), which allows fast computation with high accuracy. The interested reader can refer to [START_REF] Wang | Review of metamodeling techniques in support of engineering design optimization[END_REF], and different references cited therein for the details of various surrogate models (e.g. Kriging, radial basis function (RBF), artificial neural network (ANN), support vector machine (SVM), and polynomial chaos). In this work, the Kriging metamodeling technique is chosen for our reliability analysis of drift thanks to its flexibility in interpolating the sample points and high accuracy, notably for nonlinear problems, according to [START_REF] Wang | Review of metamodeling techniques in support of engineering design optimization[END_REF]. More precisely, a modification of the well-known AK-MCS method [START_REF] Echard | AK-MCS: An active learning reliability method combining Kriging and Monte Carlo Simulation[END_REF][START_REF] Do | Closed-form solution and reliability analysis of deep tunnel supported by a concrete and a covered compressible layer within the viscoelastic Burger rock[END_REF], as detailed in the dissertation of the first author of the present paper [START_REF] Tran | Long-term stability evaluation of underground constructions by considering uncertainties and variability of rock masses[END_REF], is chosen. This modified AK-MCS method will be briefly synthesized below as well as in the Appendix A.

Firstly developed by [START_REF] Krige | A statistical approaches to some basic mine valuation problems on the witwatersrand[END_REF] and then theorized by [START_REF] Matheron | The intrinsic random functions and their applications[END_REF], the Kriging metamodeling technique approximates the implicit LSF by carrying out a Gaussian process:

T ( ) ( ) [ ( )] ( ) ( ) ( ) p i i i g g = + Z k + Z β ≈ = ∑ x x k x x x x β ( 12
)
where the first term at the right-hand side, [k(x)] T β, represents the mean value of the Gaussian process (i.e. the trend of the process) [START_REF] Bichon | Efficient surrogate models for reliability analysis of systems with multiple failure modes[END_REF] with β and k(x) being respectively the vector of regression coefficient and the vector of basis-functions of p elements; and the second term Z(x) is assumed to have the zero-mean stationary Gaussian process and variance

d i i i i R x x θ =   ′ ′ = - -   ∏ x x θ ( 13
)
where d is the dimension of the random variable vector x, and xi is the i-th component of x.

In fact, due to the stationary character of the Gaussian process, for each dimension of the vector x, the correlation functions depend only on the relative positions of two inputs characterized by a scale parameter θi (i = 1, 2, …, d).

The calibration of the unknown parameters 2 , , an d Z σ β θ can be done through an optimization process by using the results evaluated at different observation points (also called training points) obtained from the design of experiment (DoE) step. The resolution of this optimization problem has mainly been discussed in the past and was successfully implemented in some software. In this work, the available MATLAB toolbox DACE [START_REF] Lophaven | DACE -A Matlab Kriging toolbox[END_REF] was used.

Thanks to the built Kriging metamodel ( ) g x , the exceedance probability can be estimated by using one of the well-known direct sampling methods such as the MCS:

MCS ( ) ( ) ( ) f ( ) 1 MCS 1 ( ( ) 0) 1 ( ( )), ( ( )) 0 ( ( ) 0) i N i i i i g x P I g x I g x N g x =  ≤  ≈ =  >   ∑ (14)
where NMCS counts the random samples that are generated for the estimation of failure probability by the MCS.

The probability evaluated by Eq. ( 14) depends on the Kriging metamodel, which is constructed from the results of the direct evaluation of the performance function at the training points of DoE. Thus, to derive the final result of probability, an iterative procedure is carried out. Following that, at each iteration, a subset of new training points can be chosen (based on a so-called learning function) to enrich the DoE. The process is continued until the obtained probability satisfies a stopping criterion. To simplify the presentation, our developed learning function and the chosen convergence criterion are presented in the Appendix A (the interested reader can refer to [START_REF] Tran | Long-term stability evaluation of underground constructions by considering uncertainties and variability of rock masses[END_REF] for much more details). For the sake of clarity, Table 5 summarizes the principal steps to conduct the Kriging-based reliability analysis using the modified AK-MCS method. Construct the Kriging model from S and y (by using the Kriging toolbox) 5

Interpolate g(x (i) ) (i=1, 2, ..., NMCS) from the built Kriging model and calculate the failure probability Pf (using Eq. ( 14)) 6

Verify the stopping criterion (using Eq. ( A7)) and then go to

Numerical applications

Deterministic problems

Before analyzing the results obtained from the reliability analysis, let us discuss the behavior of the concrete liner extracted from the direct numerical simulations. The results of this so-called deterministic problem are of the primary importance for understanding the drift response before accounting for the uncertainty of the mechanical properties of the host rock.

For this aim, the mean values of the short-and long-term mechanical properties of rock mass are taken in this first simulation while the compressible potential of the first support layer is fixed at v2 0 .4 ε = . As the first observation, it is worth pointing out that the numerical simulations present a maximum of equivalent stress at the inner radius of the concrete liner. At the same time, the volumetric deformation is uniform in the whole layer of the compressible support (Fig. 5). This expected result is consistent with different assumptions adopted, notably the hypothesis of isotropic hydrostatic stress at far-field and isotropic homogeneous material of components in the considered problem. By varying the thickness of the concrete liner (l2), the results exhibit a quite significant decrease of the maximum equivalent stress in this last support element (Fig. 6a andb). The mean values of the viscoplastic rock induce in the compressible liner a low volumetric deformation which is mostly smaller than the compressibility v2 0 .4 ε = (see Fig. 6c andd). By changing now the compressible potential v2 ε of the outer liner, we can observe in Fig. 7 that the maximum equivalent stress in the concrete liner seems constant when v2 0 .0 1 5 ε ≥ . The result of equivalent stress presents the highest value with the case v2 0 .0 1 5 ε =

. Notably, this latter value presents in the particular case v2 v1 ( 0 .0 1 5)

ε ε = =
in which the first support layer of drift loses the compressibility and is characterized by only one elastic modulus (Ec = Ec1 = Ec3).

To better understand the effect of rock mechanical properties on the behavior of the final support of drift, we present in Fig. 8 the results of another deterministic analysis in which a random set (noted as set1) of the COx rock properties is taken among the samples generated in the reliability analysis as detailed in the next section. More precisely, the chosen values of parameters A, B, and C, and Young's modulus E of this set are respectively equal to about 1.7, 1.3, 1.1, and 0.7 times of their mean value (i.e. Aset1 = 8.415 × 10 -6 , Bset1 = 2.647, Cset1 = 0.192, and Eset1 = 4.253 GPa). Corresponding to these mechanical properties of the massif, which are more critical than those in the case of using mean values, the results captured in Fig. 8 reveal now the important effect of εv2 on the time-dependent behavior and on the maximum equivalent stress at 100 years in the concrete liner (Fig. 8a andb). Following that, a smaller value of compressible potential εv2 induces a significant increase of compressive stress in the final support element of drift. This is because, in the case of lower compressible potential, the volumetric deformation in the compressible support exceeds this limited value between the second and third elastic stages (i.e. εV > εv2), and the compressible material works with its highest elastic modulus at the rigidification stage (Fig. 8c andd). The results of the maximum equivalent stress in the concrete liner of drift (calculated with two sets of mechanical properties of COx rock: the mean values and set1) exhibit a great difference. Thus, a parametric study is then undertaken to investigate the influence of each parameter of COx rock properties on the response of the concrete layer at 100 years. This study is investigated by varying each parameter of interest while the other mechanical properties are maintained at their corresponding minimum, maximum, or mean values. Instead of using a 95% confidential interval of the distribution, the minimum and maximum values of each mechanical parameter are taken from 10 6 random samples, which are generated for the reliability analysis as detailed in the next section. It is important to note that the anti-correlation of two parameters A and B of Lemaitre's model (see Section 2.4) is counted in these random samples. However, in the present deterministic study and for simplification purposes, we ignore the anti-correlation of these two parameters by accepting the set containing both their maximum and minimum values. In this parametric study, the inner liner's thickness and the outer liner's compressibility are fixed with l2 = 0.5 m and εv2 = 0.4, respectively.

In general, one can state that, with respect to a higher value of one of the long-term mechanical parameters (i.e. either coefficient A or B or C), the maximum equivalent stress in the concrete element tends to increase (see Fig. 9a-c). This phenomenon seems to be consistent with the fact that, under a constant deviatoric stress applied, the increase of one of these parameters yields a higher viscous deformation, according to Lemaitre's model expressed in Eq. ( 6). However, an evident tendency is difficult to be captured regarding the effect of the short-term mechanical properties of COx rock (i.e. Young's modulus), as shown in Fig. 9d. It can be explained by the interaction in a non-linear manner between the convergence of the viscoplastic rock and the reaction from the support system that generates in rock mass a redistribution of stress state represented by a variation of the deviatoric stress with time at a specific point of the host rock. This redistribution depends strongly not only on the shortand long-term mechanical properties of host rock but also on the time-dependent reaction loading of the support system, which in turn depends on the rigidity of both liners. Especially, it can be stated from this parametric study that the maximum equivalent stress in the concrete liner can vastly exceed the elastic limit of the concrete class C60/75. Note that in the typical civil engineering application, a compressible resistance about 60% of the concrete mark is mostly accepted as the elastic limit in the serviceability limit state design of the concrete elements as proposed in European Code 2 (EC2). Beyond this limit, one considers that cracking under compressive loading can occur in the concrete structure. Quantitatively, the effect of each parameter of COx rock properties on the maximum equivalent stress of the concrete liner at 100 years can be evaluated by using the sensitivity analysis. The Cotter method is chosen for this purpose thanks to its applicability for the input variables regardless of their dependence. Known as an efficient and low-cost approach to the sensitivity analysis, the Cotter method allows ranking the input parameters through their indices, called Cotter indices (see Appendix B, [START_REF] Cotter | A screening design for factorial experiments with interactions[END_REF], and [START_REF] Homma | Importance measures in global sensitivity analysis of nonlinear models[END_REF] for more details). As illustrated in Fig. 10, the obtained results of Cotter indices show that the two coefficients A and B affect the most on the concrete liner response. While the Young's modulus presents a moderate effect, and the influence of the coefficient C is lowest, which can be explained by its small COV among the four parameters of the COx rock.

Reliability analysis problem

The exceedance probability of the compressive stress in the concrete liner at t = 100 years with respect to an allowable value (i.e. the chosen elastic limit) is estimated in this part using the modified AK-MCS method. As discussed in a verification test in the Appendix A, the convergence of probability can be attained sooner (i.e. the number of iterations at convergence is lower) if a subset of new training points is used in each iteration. Thus, in the following reliability analysis, a subset of 4 enriched samples is chosen during the iterative construction process of Kriging metamodel thanks to using the parallel calculation in Code_Aster. The initial DoE of 24 training samples of the random vector x (i.e. the mechanical properties of COx rock) is generated by the Latin hypercube sampling (LHS) method. The other necessary parameters for the modified AK-MCS reliability analysis are chosen as γ = 0.01, Nγ = 6, and NMCS = 10 6 , respectively. For the sake of clarity, we remind here that the parameters γ and Nγ indicate the tolerance and the number of iterations used in the stopping criterion (see Eq. (A7) in Appendix A) while NMCS designates the number of random samples generated for the evaluation of the failure probability through the interpolation by MCS.

Fig. 11a illustrates a representative evolution of exceedance probability of the drift's final support during the iterative construction of the modified AK-MCS metamodel. Corresponding to the fixed values of the concrete liner's thickness (l2 = 0.5 m) and the compressible potential of the compressible layer (εv2 = 0.515), the convergence is obtained after 68 iterations with a total number Ncall = 296 of direct evaluations of the performance function in Code_Aster. The exceedance probability is about Pf = 0.13% in this case by adopting the allowable stress σcl = 36 MPa. This stress corresponds to the elastic limit of the concrete class C60/75 as proposed by EC2. An investigation with the smaller values of σcl is also performed. As shown in Fig. 11b, the result highlights an increase as expected of the exceedance probability in the concrete liner when the lower elastic threshold is chosen to limit its maximum equivalent stress. The exceedance probability at 100 years may reach about 98.51% if the threshold is reduced to 16 MPa. However, beyond the allowable stress of 25 MPa, the probability seems smaller than 1%. The numerical results also confirm that the number of direct evaluations of the performance function (Ncall) is more critical for the smaller probability. Now, by fixing the allowable stress of concrete and the compressibility of the compressible liner at σcl = 36 MPa and εv2 = 0.515, respectively, we can investigate the influence of each liner's thickness on the long-term stability of the final support of drift. We can observe in Fig. 12a that the exceedance probability in this concrete support varies from 15.78% (for the case of thickness l2 = 0.25 m) to 0.13% (in the case of l2 = 0.5 m). The exceedance probability is smaller than 1% when the thickness of this final support is higher than 0.33 m. The result of exceedance probability in the concrete liner with the thickness of 0.5 m can increase considerably when the thickness of the outer liner decreases from 0.2 m to 0.1 m, as highlighted in Fig. 12b.

Then by adopting l1 = 0.2 m, l2 = 0.5 m, and σcl = 36 MPa, our last numerical investigations show that the variation of the compressibility of the outer liner can favorably affect the exceedance probability in the final support of drift (Fig. 13). The long-term stability of the concrete liner is higher, representing a lower exceedance probability when the compressible potential in the first support layer increases. For the particular case of an elastic material (i.e. εv2 = εv1 = 0.515) characterized by only one elastic modulus (Ec = 100 MPa), the exceedance probability at 100 years of the concrete liner with 0.5 m of thickness can attain 30.96% instead of 0.13% (as in the case of εv2 = 0.515). This result confirms the evident effect and significant benefit of using the compressible material to ensure the stability of drift support.

Discussion

Based on different simplified assumptions, the study conducted in this work presented our preliminary results of the long-term stability analysis of drift. We can expect that the exceedance probability of the drift support can be higher than the ones obtained in this research when the uncertainty is quantified from the in situ data in which the spatial variability of rock properties and other mechanisms (e.g. damage, THM coupling) are considered.

Indeed, it has been shown in many geotechnical applications that the spatial variability of rock properties can contribute to the most significant source of uncertainty [START_REF] Chen | Reliability assessment on stability of tunnelling perpendicularly beneath an existing tunnel considering spatial variabilities of rock mass properties[END_REF][START_REF] Shokri | A review study about spatial correlation measurement in rock mass[END_REF][START_REF] Gao | Back analysis of displacements for estimating spatial distribution of viscoelastic properties around an unlined rock cavern[END_REF]. However, in the literature, this source has been taken into account in two different manners. The first approach consists of considering all uncertainty sources in the reliability analysis regardless of the heterogeneous distribution in the space of the random parameters. By using this approach in the case of drift, the surrounding rock mass is kept homogeneous, and the problem can be solved similarly to the ones of this work (i.e. by using the modified AK-MCS metamodel).

The second approach aims at accounting for the spatial effect by modeling the heterogeneous medium (also known as the variability problem to distinguish the uncertainty problem in the first approach). For this purpose, the notion of a random field is usually chosen to describe the spatial distribution of a random variable. Mathematically, this random field can be expressed in the form of a correlation function (such as Markovian or Gaussian's function) with an essential characteristic parameter, the spatial correlation length. For numerical modeling, this random field needs to be discretized, and the variability problem becomes a particular case of the uncertainty problem but with a high dimension [START_REF] Tran | Long-term stability evaluation of underground constructions by considering uncertainties and variability of rock masses[END_REF]. Indeed, the dimension of the variability problem can dramatically increase when the spatial correlation length is smaller and/or the adopted geometric model of the structure is large. Correspondingly, the computational demand increases exponentially as the number of variables grows. Solving the great dimensional problem, referred to as the curse of dimensionality, is still an active ongoing research topic in the literature. The efforts of most studies are based on the idea of using the dimension reduction method to construct the appropriate surrogate. Concerning the Kriging-based metamodeling technique, some extensions have been proposed [START_REF] Lelièvre | Développement des méthodes AK pour l'analyse de fiabilité[END_REF][START_REF] Wang | An adaptive high-dimensional model representation method for reliability analysis of geotechnical engineering problems[END_REF] to handle this problem. Although the ability of these improvements was demonstrated, notably in different academic problems, considerable limits always exist in real applications when the structure behavior is complex, and the number of random variables can range from several hundred to several thousand [START_REF] Tran | Long-term stability evaluation of underground constructions by considering uncertainties and variability of rock masses[END_REF].

By neglecting the impact of spatial distribution, the exceedance probability of underground structures provided by the first approach could be overestimated in comparison with the ones calculated by the second approach [START_REF] Griffiths | Influence of spatial variability on slope reliability using 2D random fields[END_REF][START_REF] Lü | Probabilistic assessment of tunnel convergence considering spatial variability in rock mass properties using interpolated autocorrelation and response surface method[END_REF][START_REF] Shokri | A review study about spatial correlation measurement in rock mass[END_REF][START_REF] Tran | Long-term stability evaluation of underground constructions by considering uncertainties and variability of rock masses[END_REF]. However, the computational demand of the first approach is significantly lower whilst the correspondingly overestimated result, represented by a higher probability, seems to be more appreciated for the stability analysis of the structure.

Besides, the estimated long-term exceedance probability of drift support in this work may be highly underestimated due to the ignorance of the damage mechanism as well as the anisotropy and the (thermo-) hydro-mechanical coupling effect. Note that the observation from the measurements at M/HM URL and from the creep tests showed a higher deformation rate linked with the damaged zone whilst anisotropy of convergences was largely observed in all drifts. The consideration of these effects (damage mechanism, anisotropy, coupling loading, and so forth) with the necessary to add the in situ data on the quantification process must contribute to another important uncertainty source of the COx properties, which may increase the risk of long-term stability the of drift's support.

Conclusions

In this work, the Kriging-based reliability analysis was conducted to study the long-term stability of concrete liner of a deep drift by accounting for the uncertainty of the COx claystone. The study was undertaken in the 2D plane strain condition and under the isotropic far-field stresses to model the circular drift oriented in the major horizontal stress direction of COx formation. The time-dependent mechanical behavior of this rock was assumed to be described by the classical viscoplastic model of Lemaitre. The uncertainty of the mechanical properties of COx host rock was quantified from laboratory data obtained from the triaxial creep tests. The tri-linear elastic model was proposed for the compressible liner, which was put between the concrete liner and rock mass. Based on a modification of the well-known AK-MCS method, the efficiency of the Kriging-based reliability analysis was demonstrated. The numerical results of both the deterministic problem and the reliability analysis elucidated the strong dependence of the stability of concrete liner on the uncertainty of rock properties and the compressibility of the outer compressible layer. The exceedance probability with respect to the allowable elastic limit measured at 100 years of this final support of drift decreases when its thickness increases and/or the compressible potential of the outer liner is higher. The numerical investigations confirm the tremendous benefit of compressible material on the longterm stability of drift support from a purely mechanical point of view. However, the results in this preliminary research can be underestimated when the uncertainty related to the spatial variability of COx rock properties, as well as the damage mechanism and/or the coupling effect, is neglected. The consideration of these aspects on the long-term stability analysis and optimization design of drift system support presents an essential topic for future works.

Fig. 1 .

 1 Fig. 1. 2D plane strain model of drift supported by the outer compressible and inner concrete liners in viscoplastic COx rock.

Fig. 2 .

 2 Fig. 2. (a) Stress-strain curves of a compressible material under the oedometric test (Billig et al., 2007); and (b) Tri-linear elastic model of the compressible material versus volumetric deformation.

Fig. 3 .

 3 Fig. 3. Result and histogram of Lemaitre's model parameters calibrated from 88 samples.

Fig. 4 .

 4 Fig. 4. Comparison of the experimental data (red curve) of a creep test with approximation using the mean values (blue curve) and the MCS using 1000 samples (gray curves) of three coefficients A, B, and C of Lemaitre's model (a), and extrapolation to 100 years (b).

σ

  and R indicate the constant process variance and the autocorrelation function (also known as the kernel function), respectively. The following squared exponential autocorrelation (Gaussian) function with respect to the hyperparameter vector θ is used in this work:

  Step 9; otherwise, go to Step 7 7 Identify one or a set of new training points x* using a learning function (see Eqs. (A4) and (A5)) and update the DoE and hence the matrix S 8 Solve the deterministic problem to evaluate the exact performance function of the new training points x* and update y. Then, go to Step 4 9 Calculate the COV of the failure probability: ; otherwise, go to Step 3 and increase NMCS

Fig. 5 .

 5 Fig. 5. (a) Equivalent stress (in Pa) in the concrete liner and (b) volumetric deformation in the outer compressible liner at 100 years.

Fig. 6 .Fig. 7 .

 67 Fig. 6. (a, b) Evolution of the maximum equivalent stress in the concrete liner as a function of time and its thickness; (c) Volumetric stress-strain curve in the compressible liner as a function of the concrete liner thickness; and (d) Evolution of the volumetric strain in the compressible liner as a function of time. These results are obtained by using the mean values of COx mechanical properties and the compressibility εv2=0.4 of the outer layer.

Fig. 8 .

 8 Fig. 8. (a) Evolution of the maximum equivalent stress at 100 years in the concrete liner as a function of the compressible potential εv2 of the compressible layer; (b) Evolution of the maximum equivalent stress in the concrete liner with time; (c) Volumetric strain in the compressible element versus time; and (d) Volumetric stress-strain curves in the compressible liner. A random set of the COx mechanical properties (set1) was used, while the thickness of the concrete liner l2 = 0.5 m is fixed in (b-d).

Fig. 9 .

 9 Fig. 9. Effect of each parameter of COx rock properties on the maximum equivalent stress in the concrete liner at 100 years represented by parametric studies.

Fig. 10 .

 10 Fig. 10. Effect of each parameter of COx rock properties on the maximum equivalent stress in the concrete liner at 100 years represented by the Cotter indices.

Fig. 11 .

 11 Fig. 11. (a) Evolution of exceedance probability of the concrete liner at 100 years during the iterative construction of the modified AK-MCS metamodel (σcl = 36 MPa); and (b) Exceedance probability as a function of allowable stresses σcl. The fixed values of l2 = 0.5 m and εv2 =0.515 were used.

Fig. 12 .

 12 Fig. 12. (a) Exceedance probability of the concrete liner at 100 years as a function of the inner liner thickness (with l1 = 0.2 m), and (b) of the outer liner thickness (with l2 = 0.5 m) using the fixed values of εv2 = 0.515 and σcl = 36 MPa.

Fig. 13 .

 13 Fig. 13. Influence of the compressible potential εv2 of the compressible layer on the exceedance probability of the concrete liner at 100 years (by fixing l1 = 0.2 m, l2 = 0.5 m and σcl = 36 MPa).

Table 1 .

 1 Confining pressure and deviatoric stress of seven triaxial creep tests performed on COx rock.

	Sample No.	Confining pressure (MPa)	Deviatoric stress (MPa)
	1	2	11.5
	2	2	17.25
	3	2	20.7
	4	2	20.7
	5	12	18.5
	6	12	27.75
	7	12	33.3

Table 2 .

 2 Statistical values of the Lemaitre's model coefficients calibrated from 88 samples taken from seven creep tests on COx rock.

	Coefficient	Mean	Standard deviation (std)	COV (std/mean) (%)
	A	4.88 × 10 -6	4.27 × 10 -6	87.5
	B	2.048	0.389	19
	C	0.171	0.013	7.6

Table 3 .

 3 Correlation matrix between the calibrated coefficients of Lemaitre's model.

	Coefficient	Distribution	Coefficient A	Coefficient B	Coefficient C
	A	Normal	1	-0.72	0.35
	B	Normal	-0.72	1	-0.24
	C	Normal	0.35	-0.24	1

Table 4 .

 4 Fixed values of deterministic parameters in the numerical investigations.

	l1 (m)	Ec1 (MPa) Ec2 (MPa) Ec3 (MPa) εv1	El (MPa)	νl
	0.2	100	3	100	0.015	39,100	0.2

Table 5 .

 5 Principal steps of a Kriging-based reliability analysis using the modified AK-MCS method. (i.e. solve the deterministic problem to evaluate the exact performance function at the training points of the initial DoE by using the numerical solution) 3Generate NMCS random samples for the interpolation by MCS 4

	Step No.	Description
	1	Generate NDoE initial training points (matrix S) of the DoE
	2	Determine the vector y

Declaration of competing interest

The authors wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Appendix A: The modified AK-MCS method and verification of its efficiency

(1) The modified AK-MCS method In Eq. ( 14), the performance function's evaluation at each sample of the MCS method is conducted by interpolating the result of the approximated metamodel . Concretely, the performance function for each realization of the random input vector x can be predicted from which follows the normal distribution. The mean and variance functions 2 ( ( ) a n d ( ))

g g µ σ

x

x of the Kriging predictor are derived as [START_REF] Lophaven | DACE -A Matlab Kriging toolbox[END_REF][START_REF] Echard | AK-MCS: An active learning reliability method combining Kriging and Monte Carlo Simulation[END_REF][START_REF] Gaspar | A study on a stopping criterion for active refinement algorithms in Kriging surrogate models[END_REF]:

where

where the vector y gathers the response of the structure (i.e. the exact results of the performance function) evaluated at NDoE training points of DoE defined in the matrix S. Besides, R is the matrix of correlation between each pair of points of DoE while K and r(x) are the regression matrix and the vector of cross-correlation between the prediction point x and each training point of DoE:

As noted in Section 4, in the iterative process to derive the final result of failure probability, new training points must be selected by a learning function to enrich the DoE after each iteration. In the wellknown AK-MCS method, [START_REF] Echard | AK-MCS: An active learning reliability method combining Kriging and Monte Carlo Simulation[END_REF] proposed the following U learning function:

That measures, in fact, the probability of making a mistake on the sign of ( ) g x by replacing ( ) g x by ( ). g x Based on this function, the sample x*, which gives the lowest U value among NMCS samples (i.e. x* = argmin{U(x)}), is selected as a new training point. This point usually corresponds to the nearest point to the limit state. With this learning function, [START_REF] Echard | AK-MCS: An active learning reliability method combining Kriging and Monte Carlo Simulation[END_REF] proposed that the iterative procedure can be stopped when the smallest U value of NMCS samples is higher than 2 (i.e. min{U(x)} > 2).

Although the accuracy of the AK-MCS method has been largely demonstrated [START_REF] Fauriat | AK-SYS: An adaptation of the AK-MCS method for system reliability[END_REF][START_REF] Huang | Assessing small failure probabilities by AK-SS: An active learning method combining Kriging and Subset Simulation[END_REF][START_REF] Liu | Efficient optimization of reliability-constrained structural design problems including interval uncertainty[END_REF], some scholars [START_REF] Gaspar | A study on a stopping criterion for active refinement algorithms in Kriging surrogate models[END_REF][START_REF] Hu | Global sensitivity analysis-enhanced surrogate (GSAS) modeling for reliability analysis[END_REF] showed that stopping criterion (i.e. min{U(x)} > 2) can be too conservative whilst the chosen new training points based on argmin{U(x)} is not optimal. A more comprehensive discussion on the limitation of the original AK-MCS method [START_REF] Echard | AK-MCS: An active learning reliability method combining Kriging and Monte Carlo Simulation[END_REF] can be found in the work of [START_REF] Hu | Global sensitivity analysis-enhanced surrogate (GSAS) modeling for reliability analysis[END_REF].

Besides, it has been discussed in many other contributions [START_REF] Basudhar | Adaptive explicit decision functions for probabilistic design and optimization using support vector machines[END_REF][START_REF] Xiao | A new adaptive sequential sampling method to construct surrogate models for efficient reliability analysis[END_REF][START_REF] Zhao | A novel active learning framework for classification: Using weighted rank aggregation to achieve multiple query criteria[END_REF] that the best candidate point to enrich the DoE should be not only close to the limit state but also far away from the training samples of existing DoE by verifying a distance constraint. In [START_REF] Zhao | A novel active learning framework for classification: Using weighted rank aggregation to achieve multiple query criteria[END_REF], a distance constraint dmin > D was proposed in which the limit distance parameter D varies during the add-sample process, whilst dmin is the minimum distance from the candidate point to training points of DoE.

With the aim to select suitable new training points, a modification of the U learning function is proposed in [START_REF] Tran | Long-term stability evaluation of underground constructions by considering uncertainties and variability of rock masses[END_REF]. Following that, a new training point x* will be chosen from NMCS random samples which has the lowest U value and firstly verifies the following distance constraint:

where x (i) presents the existing training points of the current DoE; and the limit distance parameter D, which can vary after each iteration, is defined as

m in (m ax( )) ( 1,2,..., )

where d denotes the dimension of the vector of hyperparameters θ, and Nint is the current number of iterations.

It can be seen from Eq. (A5) that the sample with the lowest U value (i.e. argmin{U(x)}) may not be chosen as the new training point if it does not verify the distance constraint. Instead, as the principal idea, this modified U learning function aims to find out between the closest points to the limit state, the new training point that is far enough from the existing training points of the current DoE by verifying the distance constraint in Eq. (A5).

In addition to this modified U learning function, an appropriate stopping criterion proposed in the work of [START_REF] Gaspar | A study on a stopping criterion for active refinement algorithms in Kriging surrogate models[END_REF] is also chosen in this work. This stopping criterion is based on a defined stabilization condition of failure probability:

( 2,3,..., )

where (1) f P is the reference failure probability that is used to detect the stabilization, and ( ) f i P is the failure probability in the next Nγ -1 iterations.

(2) Verification test on the efficiency of the modified AK-MCS A numerical test is conducted to verify the efficiency of the modified AK-MCS. The chosen test consists of measuring the failure probability of a deep tunnel supported by two liners installed in sequence in a viscoelastic Burgers rock. In this context, the results ( ) g x ( ) g x obtained from the direct MCS of Do et al. (2020), who investigated the long-term stability of this underground structure by using an analytical solution, will be served as the reference solution. Besides, by comparing with the results provided by the initial AK-MCS method [START_REF] Echard | AK-MCS: An active learning reliability method combining Kriging and Monte Carlo Simulation[END_REF], the improvement and efficiency of our proposed Kriging metamodel are highlighted.

The test is carried out under the same hypothesis in the work of Do et al. (2020) whilst the LSF defined in Eq. ( 11) is considered. The random variable vector x in this LSF now represents four parameters of the Burger model of the rheological rock (x = (GM, ηM, GK, ηK)). All the necessary values of the parameters included in this numerical test are chosen similarly as in Do et al. (2020).

For the Kriging-based reliability analysis, an initial DoE with 24 training points generated quasi-uniformly by the LHS method is used. Further, for the modified AK-MCS metamodel, the values of γ = 0.01 and Nγ = 6 of the stopping criterion are used. The number of random samples for the direct MCS (used as the referent results) and for the interpolation from the constructed metamodels is also taken as NMCS = 10 6 .

The comparison of the obtained results from the Kriging methods with the reference result of the direct MCS in Table A1 confirms the accuracy of these metamodel predictions. A quite similar failure probability measured at 100 years in the second liner but with the lowest number of calls (Ncall) to the performance function (i.e. the resolution of the direct problem) in the modified AK-MCS method demonstrates its improvement with respect to the initial AK-MCS method (see also Fig. A1a).

An investigation on the effect of the chosen number of new training points to enrich the DoE in each iteration on the failure probability by the modified AK-MCS method is also conducted. The results (as illustrated in Fig. A1b) show that, in general, an increase of the new training points in each iteration reduces the total number of iterations at convergence (with the same adopted stopping criterion). However, this last number reaches an asymptotic value when this number of new training points in the added subset of DoE in each iteration becomes significant. It seems that a subset of four new training points in each iteration can be an appropriate choice. This observation is especially interesting in practical engineering, in which many evaluations of the performance function can be numerically conducted through parallel computations. Although the number of calls to the performance function Ncall increases by using a subset of new training points throughout the iterative process, the real time-consuming for the numerical simulation decreases.

Appendix B: Sensitivity analysis by Cotter's method

By ranking input variables, the Cotter method aims at isolating the essential parameters on the structure response. This sensitivity analysis approach is considered as a low-cost method when only 2N + 2 model evaluations are conducted (N is the number of input variables 1 2 ( , ,..., ))

N

x x x [START_REF] Cotter | A screening design for factorial experiments with interactions[END_REF][START_REF] Homma | Importance measures in global sensitivity analysis of nonlinear models[END_REF]. Following that, only the minimum (i.e. low level) and maximum (i.e. high level) values of each input variable are used in these evaluations. Specifically, in the first simulation, all parameters take their minimum value (i.e.