Dr Dang-Bao-An Tran

Dr Tat-Thang Nguyen

Dr Anh- Khoa Nguyen

Dr Huy

Keywords: uncertainty, spatial variability, reliability analysis, Kriging metamodel, deep tunnel, compressible material, COx claystone

Most of all, I would like to express my deep gratitude to my wife, Thi-Nhu NGUYEN, for her great patience, endless love, and support through all the ups and downs; to my daughter and my son for making my life so much love, fun

RÉSUMÉ ÉTENDU EN FRANÇAIS

Dans la pratique quotidienne de la science de l'ingénieur, il y a un intérêt croissant pour développer des stratégies et des méthodes de caractérisation et de modélisation de la propagation des incertitudes, de la variabilité et des systèmes aléatoires. En particulier, la l'analyse de fiabilité est devenue un outil puissant dans le processus de conception, d'optimisation de la maintenance et de la prédiction de la durée de vie des structures. Pour les matériaux conçus par l'Homme, actuellement utilisés dans l'industrie, de nombreux efforts sont faits pour assurer leurs propriétés finales, et maitriser les incertitudes sur ces propriétés, si bien que celles-ci, pour un produit donné sont indiquées en même temps que les propriétés ellesmêmes. En pratique de l'ingénierie, l'une des façons de prendre en compte ces incertitudes consiste à utiliser les règles de conception tels que les Eurocodes, l'ASTM, etc. Par exemple, dans diverses approches de conception promues par les Eurocodes, les incertitudes et la variabilité des charges et des propriétés des matériaux sont prise en compte via la méthode des facteurs de sécurité partiels.

Cependant, de telles méthodes ne sont pas imposées à la construction souterraine, malgré leur application dans d'autres domaines variés de la géotechnique (les Eurocodes 7, le code européen de conception géotechnique indiquent explicitement que le domaine de la conception souterraine est hors de ses domaines d'application). Cela s'explique en grande partie par la situation légèrement différente des ouvrages souterraines par rapport à d'autres types de structures. Les sols, et plus souvent les massifs rocheux qui abritent ces ouvrages, contrairement aux matériaux artificiels, sont des matériaux géologiques, naturels et hétérogènes, produits par des processus de la sédimentation, de la diagenèse et du tectonisme. Leur structure interne aléatoire est une combinaison de minéraux constitutifs qui régissent le comportement mécanique local en tout point des massifs. Cette organisation spatiale des minéraux (et pour cette raison, des propriétés) rend la caractérisation des propriétés des géomatériaux plus complexe car les incertitudes épistémiques et aléatoires doivent être qualifiées. L'incertitude aléatoire fait référence à la variabilité naturelle du phénomène considéré et l'incertitude épistémique à un manque de connaissances. L'incertitude épistémique associe tous les paramètres d'entrée de conception géotechnique, tels que le champ de contraintes in situ, les paramètres de résistance du massif rocheux, les modules de déformation, etc., qui n'ont jamais été connus avec précision auparavant. Mais le manque de connaissances ne se limite pas aux incertitudes sur les mesures des propriétés physiques (erreur de mesures et incertitudes inhérentes à la métrologie d'une mesure). En particulier, en géomécanique ce manque est principalement dû au fait qu'il y aura toujours un nombre insuffisant de données pour comprendre et caractériser pleinement les phénomènes en jeu, entre autres (et peut-être surtout) en raison de limitations physiques d'avoir plus qu'un échantillon dans un point d'un massif. En outre, les résultats des analyses déterministes peuvent manquer les mécanismes exacts de défaillance et ignorer la réponse précise de la partie la plus faible de la roche au sens de ses propriétés de caractère aléatoire. De plus, des échantillons provenant d'essais au laboratoire ou in-situ sont plus ou moins perturbés et des erreurs sont introduites dans la mesure et l'estimation de la résistance et de la déformation caractéristiques des roches.

Bien que toutes ces caractéristiques particulières des géomatériaux et des massifs rocheux expliquent à elles seules les besoins d'approches particulières pour évaluer et réduire les incertitudes dans la conception des ouvrages souterraine, elles deviennent plus évidentes lorsque la question est posée pour le ouvrages souterrains de stockage des déchets radioactives. Le stockage souterrain des déchets radioactifs étant une option préférentielle dans de nombreux pays, les questions d'évaluation des incertitudes liées à de telles structures suscitent une attention croissante. Outre les sources d'incertitudes mentionnées ci-dessus, inhérentes à toutes les constructions souterraines, le couplage multiphysique et le comportement à long terme des vii massifs rocheux ainsi que les incertitudes liées au comportement des structures dans de telles conditions sont peut-être les points centraux de la conception des installations de stockage et représente un défi pour les chercheurs et les ingénieurs. D'une part, les études spatiales (quelques dizaines de kilomètres carrés) et temporelles (quelques millions d'années) sur le domaine du stockage ne sont comparables à aucune construction souterraine construite jusqu'à présent. En revanche, la complexité de l'interaction de nombreux phénomènes de nature diverse (mécanique, chimique, hydraulique, thermique, radioactive) et cinétique n'est jamais vue auparavant dans un seul problème. Et pourtant, si l'on considère un seul problème physique (disons mécanique), le problème se ramifie dans d'autres problèmes et mécanismes qui pourraient être actifs à certains moments de la vie de l'excavation souterraine mais qui impactent le comportement les prochaines phases de vie des ouvrages souterrains. Par exemple, le comportement instantané irréversible des roches environnantes contrôlera l'amorçage et la croissance de la zone endommagée d'excavation autour d'une cavité. Ceci, à son tour, aura un impact sur la distribution des contraintes et des déformations ainsi que sur le champ hydraulique autour de la cavité pour, in fine, avoir un impact sur le comportement à long terme de l'excavation. Il est intuitivement clair que toute incertitude sur le comportement irréversible de la roche à tout moment de la vie de l'excavation se propagera aux incertitudes sur les phases futures de la vie de l'excavation.

Les travaux de cette thèse font partie de nombreux autres travaux réalisés par la communauté scientifique concernant le projet du centre de stockage profond de déchets radioactifs en France baptisé par Cigéo, qui sera implémenté dans la couche de l'argilite du Callovo-Oxfordien (COx) aux départements Meuse/Haute-Marne, s'il est autorisé.

Nous n'avons pas l'ambition ici de traiter toutes les sources d'incertitudes dans un si grand projet et même pas toutes celles liées à toute sorte d'excavations souterraines. Notre principale préoccupation dans ce travail est la propagation des incertitudes du comportement mécanique à long terme des massifs rocheux sur la réponse d'une ouvrage souterraine typique.

Les résultats publiés des nombreuses études effectués ces dernières décades ont révélé une incertitude des valeurs des propriétés mécaniques instantanées de l'argilite du Callovo Oxfordien (COx) et la variabilité de ces propriétés en fonction de la profondeur. La même variabilité est également observée sur le comportement à long terme de la roche qui semble être activée à faible contrainte dévia torique. Ces variations des propriétés mécaniques suivent la variation de la composition minéralogique et de l'anisotropie structurale. Par exemple, les levés sismiques des projets de l'Andra montrent également des variations verticales et latérales des propriétés de l'argilite COx. En particulier, ces incertitudes ont un impact sur l'amorçage et l'extension de la zone fracturée due à l'excavation.

La particularité des travaux présentés ici, par rapport aux analyses de fiabilité généralement effectuées sur la mécanique des structures, réside principalement dans la prise en compte du comportement dépendant du temps des roches conduisant à une probabilité de rupture en fonction du temps.

Après une analyse des approches alternatives existant qui pourrait être utilisées d'une part pour l'évaluation des incertitudes liées aux propriétés mécaniques de l'argilite COx et à leur variabilité spatiale et d'autre part, pour quantifier leurs impacts sur la réponse d'une structure de stockage souterraine typique, cette thèse propose une méthodologie optimisée pour les analyses de fiabilité d'un ouvrage souterrain. Cette approche est tout d'abord utilisée pour évaluer l'impact des incertitudes des paramètres mécaniques du COx sur la probabilité de dépassement de la contrainte admissible d'un revêtement d'une galerie. Ensuite, l'impact de la variabilité de ces paramètres sur la probabilité de dépassement de la contrainte admissible à long terme est discuté. Pour la variabilité spatiale, les champs aléatoires gaussiens non corrélés viii et corrélés sont traités. Particulièrement, la méthodologie développée dans cette étude peut être suivie pour traiter l'incertitude dans de nombreux domaines, tels que l'analyse des risques naturels (par exemple, les tremblements de terre, les inondations) et la prise de décision concernant les risques environnementaux. Nous utilisons également le Code-Aster basé sur des scripts Python pour étudier les modèles d'hétérogénéité des matériaux.

La thèse est organisée en quatre chapitres.

Le premier chapitre présente le contexte et l'état de l'art de l'analyse fiabiliste en géotechnique, en particulier pour les ouvrages souterrains dans des massifs rocheux, en tenant compte de l'incertitude et de la variabilité de leurs propriétés mécaniques. Dans ce chapitre bibliographique on réalise une revue de quelques concepts de base de l'analyse de fiabilité utilisés plus tard. Les incertitudes géotechniques sont généralement classées comme épistémiques ou aléatoires. La première incertitude (c'est-à-dire les incertitudes épistémiques) concerne principalement l'erreur de mesure et l'incertitude de transformation. L'erreur de mesure associée à la rareté des données (qualité, quantité, spécification) provient des différentes étapes de l'échantillonnage, de la préparation, du transport et de la procédure de test, liées au manque de connaissances ou de données sur le phénomène observé qui en résulte. Cette source d'incertitude peut être réduite à une valeur minimale en considérant un plus grand nombre d'échantillons ; ce minimum dépende des erreurs systématiques impliquées dans l'ensemble du processus. L'incertitude de transformation est introduite lorsque les mesures sur le terrain ou au laboratoire sont transformées en propriétés du sol de référence à l'aide de modèles de corrélation empiriques ou autres. Cette incertitude peut être réduite en considérant des modèles empiriques ou mathématiques plus raffinés.

L'autre source d'incertitude (c'est-à-dire l'incertitude aléatoire) est connue sous le nom d'incertitude irréductible ou inhérente. La variabilité inhérente au sol qui résulte principalement de processus géologiques naturels est de nature aléatoire. Ce type d'incertitude est associé au caractère aléatoire "inhérent" des processus naturels, qui se manifeste par la variabilité dans le temps pour les phénomènes qui se produisent en un seul endroit (variabilité temporelle), ou par la variabilité dans l'espace pour les phénomènes qui se produisent en différents endroits mais à un seul point dans le temps (variabilité spatiale), ou par variabilité à la fois dans le temps et dans l'espace. Pour prendre en compte l'effet de variabilité spatiale des propriétés des géomatériaux, le concept de champ aléatoire est largement utilisé. Selon ce concept, à n'importe quel endroit de la formation, les propriétés des géomatériaux sont des variables aléatoires caractérisées par des distributions de probabilité et corrélées avec celles des zones adjacentes. Cette corrélation des propriétés entre les zones adjacents (appelée également la fonction de corrélation) est caractérisée par les longueurs de corrélation qui indiquent la plage de fluctuation du paramètre ou décrivent la variabilité spatiale d'une propriété de géomatériaux dans les directions horizontales et verticale. Deux points dans l'espace séparés par une distance supérieure à la longueur de corrélation sont largement non corrélés. La plus petite longueur de corrélation, la plus rapide la fluctuation du paramètre de sol est ; sinon, le paramètre est relativement stable dans une large plage. L'analyse de fiabilité consiste à quantifier la propagation des incertitudes sur la probabilité que la réponse de la structure dépasse un critère de défaillance défini mathématiquement par la fonction d'état limite LSF (Limit State Function ou fonction de performance). Dans les problèmes de fiabilité, les variables d'entrée incertaines, appelées variables aléatoires, et leurs espaces sont divisés en régions de défaillance et de sécurité limitées par le LSF. Le concept de LSF est générique et peut être transféré à une grande famille de problèmes. Il existe également différents types de fonctions de performance. Dans le contexte des ouvrages souterraines, les fonctions de performance bien connues sont généralement liées à la convergence des tunnels, à la capacité de support et au rayon plastique. La quantification de la propagation de ix l'incertitude pour estimer la probabilité de défaillance est le principal objectif de l'analyse de la fiabilité structurelle. Diverses méthodes d'analyse de fiabilité ont été développées et peuvent être grossièrement classées en trois groupes principaux : méthodes de fiabilité locales, méthodes d'échantillonnage et méthodes de fiabilité globale. Les méthodes de fiabilité locales, telles que FORM et SORM, sont bien connues dans les approches analytiques pour approximer la surface de l'état limite autour du point de défaillance le plus probable. Cependant, ces méthodes ne sont pas recommandées pour les fonctions hautement non linéaires et non différenciables car elles doivent calculer le gradient de la LSF. Dans les méthodes de fiabilité d'échantillonnage, la simulation de Monte Carlo (MCS) est peut-être la méthode le plus connue pour estimer la probabilité de défaillance en fonction du nombre de la valeur négative de la fonction de performance sur le nombre total d'essais. Contrairement aux méthodes de fiabilité locales, MCS est une méthode tout à fait robuste et universelle, elle pourrait donc être utilisée pour vérifier d'autres méthodes de fiabilité pour estimer la probabilité de défaillance des structures. Bien qu'il soit robuste et précis, MCS est une méthode généralement couteuse en temps de calcul en raison de l'exigence d'un grand nombre de réalisations. Cela est particulièrement vrai lorsqu'une réalisation de MCS est faite via l'utilisation des méthodes d'éléments finis ou différences finies. Plus la probabilité de défaillance est faible, plus la taille de l'échantillon est grande dans la méthode de Monte Carlo pour garantir la même précision de calcul (c'est-à-dire que l'évaluation des faibles probabilités avec MCS est peu efficace). Ainsi, des approches probabilistes plus avancées nécessitant un nombre d'appels limité au modèle mécanique sont nécessaires. Ces approches peuvent surmonter les inconvénients des méthodes d'échantillonnage MCS et pourraient mieux s'adapter à la non-linéarité élevée et la réponse chronophage d'un modèle original par un métamodèle (c'est-à-dire une approximation fonctionnelle ou une équation analytique). Par conséquent, ces approches permettent de réduire significativement le nombre d'appels au modèle mécanique coûteux en calcul, comme dans la méthodologie MCS brute. Le principe de la méta-modélisation consiste à proposer une fonction purement mathématique, qui remplace le modèle numérique souvent coûteux représentant le comportement du système ou de la structure étudiée. L'idée principale du métamodèle est d'accélérer le calcul et d'améliorer la précision de l'analyse de fiabilité, en particulier dans les problèmes plus complexes et spécifiques (par exemple, non linéaire élevé, la valeur minuscule de la probabilité de défaillance, etc.). Un métamodèle est calibré à partir d'un ensemble de points, appelé plan d'expérience (DoE), pour lesquels le modèle numérique a été évalué. Le principal défi d'une telle analyse de fiabilité basée sur un substitut réside dans la façon de construire un métamodèle adéquat pour approximer le LSF sans gaspillage d'échantillons et d'appels de fonction. Parmi divers métamodèles (par exemple, le Krigeage, les fonctions de base radiale (RBF), le réseau de neurones artificiels (ANN), la machine à vecteur de support (SVM), le chaos polynomial, ...), le métamodèle de Krigeage a l'avantage de préserver la variance de modèle de base et semble mieux adapté à la problématique étudiée ici. De plus, sa flexibilité dans l'interpolation des points d'échantillonnage permet de combiner cette technique avec la méthode d'échantillonnage classique, comme le MCS, pour mesurer la probabilité de défaillance. Nous avons adopté pour ces travaux une version bien connue de cette combinaison sous le nom de la méthode AK-MCS, avec une extension / amélioration récente. Bien que la précision et l'efficacité de la méthode AK-MCS aient été intensivement démontrées dans de nombreux projets de conception structurelle, elle est rarement appliquée en géotechnique, en particulier dans le domaine de la mécanique des roches. Cette étude bibliographique a également révélé qu'une grande partie des études fiabilistes géotechniques concernent les structures géotechniques construites dans les sols. Seuls certains chercheurs ont effectué un calcul probabiliste et utilisé des outils de conception pour les supports de tunnel dans le massif rocheux. De plus, toutes ces études sur la stabilité des ouvrages souterrains profonds ont traité de l'impact de l'incertitude et / ou de la variabilité des propriétés de la roche hôte sur la stabilité à court terme en considérant uniquement le comportement instantané des roches hôtes. De ce x fait les résultats sont d'une pertinence limitée pour la conception des ouvrages de stockage car une large gamme de roches présente une réponse dépendant du temps important, qui dans certains cas peut contribuer à plus de 70% de la convergence totale. Du point de vue déterministe, la prise en compte de l'effet dépendant du temps sur l'analyse et la conception du tunnel peut être réalisée en adoptant une loi de comportement mécanique dépendante du temps approprié du massif rocheux. De nombreux modèles constitutifs, viscoélastiques (linéaire ou non linéaire) ou élasto-viscoplastique, ont été présentés dans la littérature. Ils visent à simuler avec précision les résultats observés en laboratoire (généralement à travers les tests de fluage ou de relaxation uniaxial et / ou triaxial), ou la mesure de convergence de tunnel réalisée in situ. Des modèles de plus en plus sophistiqués permettant de rendre compte non seulement du phénomène élasto-viscoplastique mais aussi de l'effet anisotrope ont été présentés dans la littérature. Si ces modèles correspondent assez bien à différentes observations, la détermination de nombreux paramètres impliqués dans les modèles présente un défi considérable si on tient compte des résultats expérimentaux limités disponibles. Ce défi est plus important encore dans le contexte de l'analyse de fiabilité lorsque la quantification de l'incertitude de l'ensemble des paramètres des modèles complexe devint très difficile, voire impossible. Par conséquent, pour l'analyse de fiabilité dépendante du temps, un modèle simple de roches, capable de décrire les principales caractéristiques du comportement différé, mais impliquant un nombre réduit de paramètres est préférable. En effet, la prise en compte de l'effet du temps sur l'analyse de fiabilité du tunnel dans la masse rocheuse n'a pas encore été menée dans la littérature. La difficulté de ce type de problème est multiple. Premièrement, comme mentionné ci-dessus, la quantification de l'incertitude des paramètres est difficile en raison des résultats expérimentaux limités. Deuxièmement, la durée du fluage-test en laboratoire (quelques heures à plusieurs mois) ou la mesure in situ (plusieurs mois voire plusieurs années) est très courte par rapport à la durée de vie du tunnel (conçu parfois pour plusieurs dizaines voire centaine d'années). Ainsi, la quantification de l'incertitude des propriétés mécaniques en fonction du temps de la masse rocheuse à partir de ces mesures, qui sera utilisée pour extrapoler le résultat à très long terme, doit contribuer à une incertitude (encore plus significative).

Dans le deuxième chapitre, l'analyse de fiabilité d'un tunnel profond creusé dans une roche viscoélastique linéaire a été effectuée. Le tunnel profond étudié se trouve dans une roche viscoélastique homogène, isotrope et incompressible soumise à une contrainte hydrostatique à l'infini. Le tunnel est soutenu par un système de deux revêtements élastiques installés successivement. L'excavation séquentielle et l'effet de l'avancement du front de taille du tunnel sur la section du tunnel considérée sont également pris en compte dans cette étude. Le comportement en fonction du temps du massif est décrit par le modèle viscoélastique de Burgers schématisé par une connexion en série d'un élément Kelvin et d'un élément Maxwell. La procédure pour dériver la solution analytique du comportement du tunnel est basée sur la solution de l'équation intégrale. La validation de la solution développée se fait par comparaison avec les résultats obtenus à partir de la simulation numérique dans le Code Aster. Une étude paramétrique est ensuite menée dans le contexte d'un problème déterministe pour mettre en évidence l'effet de différents paramètres sur la réponse du tunnel. L'accent est mis d'une part sur la vitesse d'avancement du tunnel, le temps d'installation et l'épaisseur des revêtements et, d'autre part, sur les propriétés viscoélastiques du massif rocheux (caractérisées par quatre paramètres du modèle Burger). Ensuite, en utilisant les valeurs moyennes des propriétés viscoélastiques du massif, on étudie tout d'abord la dépendance de la réponse du tunnel de la vitesse d'avancement du tunnel, les temps d'installation de soutènement et de revêtement et l'épaisseur des revêtements. Les résultats montrent qu'une excavation plus rapide conduit à une convergence plus prononcée du tunnel dans le temps. Ce résultat est attendu : un taux d'avancement significatif induit une diminution rapide du confinement et donc du déplacement radial à la surface du tunnel. En termes de contrainte dans les revêtements, la mise en place xi plutôt des éléments de soutènement induit une augmentation de ces contraintes surtout dans le revêtement intérieur. En plus, la convergence du tunnel et des contraintes de compression dans les revêtements diminuent si une épaisseur plus élevée de chaque revêtement est appliquée. En ce qui concerne l'influence des propriétés du massif, le coefficient de viscosité associé à l'amortisseur de Maxwell ne présente que son effet à long terme, comme on peut s'y attendre, lorsque la valeur plus élevée de ce paramètre implique une diminution de la convergence du tunnel et la contrainte équivalente dans les revêtements. L'effet du coefficient de viscosité associé à l'amortisseur de Kelvin semble faible en régime transitoire (à court terme) et à long terme, la contribution de ce paramètre sur le comportement du tunnel et de ses supports n'est pas significative. Le coefficient associé au ressort de Maxwell, le paramètre caractérise le comportement instantané du massif, présente un effet remarquable sur la convergence du tunnel. Une augmentation de ce paramètre provoque dans la roche une déformation initiale plus faible et, par conséquent, une convergence totale plus faible à la surface du tunnel. Les résultats numériques montrent également une variation de la contrainte équivalente dans les revêtements en fonction de ce module élastique associé au ressort de Maxwell. Une valeur plus élevée de ce dernier paramètre peut introduire une augmentation de la contrainte de compression, notamment dans le deuxième revêtement. Parmi les quatre paramètres du modèle Burger, il semble que le module élastique associé au ressort de Kelvin présente l'effet le plus significatif : une variation de ce dernier paramètre peut générer dans le massif ainsi que dans les revêtements une perturbation notable en termes de déplacement et de contrainte. En effet, plus la valeur de ce paramètre est élevée, plus la convergence et la contrainte équivalente dans les deux revêtements sont faibles. L'effet des paramètres de Burgers avec leur incertitude a ensuite été étudié dans l'analyse de fiabilité dépendante du temps. Deux critères d'évaluation sont considérés : le premier correspond à une convergence du tunnel dépassant une certaine valeur tandis que le second présente un critère en contrainte et correspond le cas où l'état de contrainte dépasse une certaine valeur limite de contrainte admissible des matériaux constitutifs. Grâce à la solution analytique développée, la simulation MCS a été utilisée pour évaluer la probabilité de dépassement de seuil de la structure. Cette estimation est utilisée comme solution de référence pour étudier la précision des nouvelles techniques de méta modélisation de Krigeage adoptée par la suite. L'analyse de fiabilité dépendante du temps par MCS montre que la probabilité de dépassement du critère augmente rapidement au cours des premières années avant d'atteindre un taux plus faible à l'étape suivante. Ce résultat, comme attendu, est la conséquence du comportement monotone de la convergence des tunnels dans le temps. En fonction de la vitesse d'avancement du tunnel, des temps d'installation et de l'épaisseur des revêtements, la forte variation de la probabilité de dépassement à chaque instant a été remarquée. Concrètement, correspondant à une excavation rapide du tunnel, les résultats présentent une probabilité de dépassement plus élevée. De manière générale, une augmentation d'épaisseur de chaque revêtement permet de réduire la probabilité de dépassement du déplacement radial au niveau de la paroi du tunnel. On peut également observer qu'une vitesse d'excavation importante du tunnel entraîne une probabilité de dépassement de critère choisi plus élevée dans le second élément de support du tunnel. Par ailleurs, lorsque le premier revêtement est posé immédiatement après l'excavation (cas du béton projeté très souvent utilisé en pratique) on obtient la probabilité maximale de dépassement du critère dans cet élément. Néanmoins, pour une vitesse d'excavation choisie, une installation tôt de chaque revêtement augmente sa probabilité de dépassement des valeurs limites. Cette probabilité est plus faible dans le premier élément de support lorsque le deuxième revêtement est mis en place le plus tôt possible après l'installation du premier. Cependant, il génère dans ce second revêtement une probabilité de dépassement de seuil plus importante. Les applications numériques montrent que, une variation modérée des épaisseurs des revêtements peut générer un changement significatif de cette probabilité dans chaque élément de support. Cette forte diminution de la probabilité est compatible avec une réduction notable de la contrainte en compression dans xii deux revêtements. En considérant les valeurs admissibles plus faibles de la contrainte de compression dans chaque revêtement, nous notons comme attendu une probabilité plus élevée (pour l'épaisseurs des revêtement fixée). Grâce à ces résultats de l'analyse de fiabilité de MCS, la précision et l'efficacité du métamodèle Krigeage choisi ont été vérifiées. Une modification de la méthode bien connue AK-MCS est proposée. Cette modification, consiste à faire un choix optimal les nouveaux échantillons d'apprentissage améliorant et accélérant la procédure de calcul de la probabilité de dépassement avec un nombre d'échantillon plus petite du DoE (indiquant ainsi l'efficacité de la fonction d'apprentissage). Dans ce travail, la fonction d'apprentissage U modifiée est proposée pour sélectionner un nouveau point d'apprentissage approprié qui vérifie deux conditions : le nouveau point sera près de la surface limite et loin des points d'apprentissage existant du DoE. Ainsi, à la différence de la méthode AK-MCS initiale, parmi les échantillons candidats on ne choisit pas forcément celui dont la valeur U est la plus basse. Au lieu de cela, la fonction d'apprentissage U modifiée vise à trouver, parmi les points les plus proches de l'état limite, le nouveau point d'entraînement qui est suffisamment éloigné des points d'entraînement existants du DoE actuel en vérifiant la contrainte condition de distance minimale. Cette distance est calculée comme la distance minimale entre les hyperparamètres du métamodel entre l'itération précédent et l'itération actuelle. Cette distance diminue avec l'augmentation du nombre des itérations. Une valeur élevée de ce paramètre de distance aux premières itérations permet de rechercher le nouveau point d'apprentissage dans l'espace des variables d'entrée dans un domaine suffisamment large, permettant ainsi d'éviter des valeurs minimisant la fonction U et qui de ce fait n'ont donc qu'un effet mineur sur la variation du métamodèle de Krigeage. La validation et l'amélioration significative de cette méthode AK-MCS modifié ont été démontrées dans deux exemples : d'abord dans un problème académique puis dans le cas d'un tunnel creusé dans la roche viscoélastique linéaire de type Burger. Nos vérifications numériques ont également mis en évidence l'avantage d'utiliser un sous-ensemble de nouveaux points pour enrichir le DoE aidant à réduire le nombre d'itérations pour atteindre la convergence de calcul de probabilité de dépassement. Ce fait est utile pour les autres études des chapitres suivants où l'évaluation de la fonction de performance doit être effectuée numériquement avec un calcul parallèle disponible.

Dans le troisième chapitre, la méthode AK-MCS modifiée a été appliquée pour étudier la stabilité à long terme d'un revêtement d'une alvéole profonde dans le contexte du stockage des déchets radioactifs dans l'argilite du Callovo Oxfordien (COx). Dans ce contexte, la stabilité du revêtement en béton de l'alvéole pendant la période d'exploitation d'environ 100 ans est cruciale pour assurer la fonctionnalité du système mécanique à l'intérieur de la galerie. Grâce à ses propriétés favorables (faible conductivité hydraulique, faible diffusion moléculaire, forte capacité de rétention des radionucléides), la roche COx est considérée comme une formation géologique potentielle pour un stockage profond des déchets radioactifs en France. Pour mieux comprendre le comportement de cette roche argileuse, des programmes de recherche intensifs ont été menés combinant des essais au laboratoire, des observations in-situ et des simulations numériques. Parallèlement à de nombreux programmes d'essais au laboratoire (uniaxial/triaxial, fluage et relaxation), des observations in situ ont été réalisées pour caractériser le comportement mécanique et hydromécanique à court et long terme des roches COx. Depuis 2000, l'Agence nationale de gestion des déchets radioactifs (Andra) réalise des études dans un Laboratoire de Recherche Souterrain (URL) à Bure (à près de 300 km à l'Est de Paris) des départements de la Meuse et de la Haute Marne (M/HM). À l'URL M/HM, diverses expériences scientifiques et technologiques ont été menées pour les propriétés/comportements le la roche hôte COx, ainsi que démontrer la faisabilité de la construction/opération d'un stockage profond de déchet radioactive. En particuler, la réponse hydromécanique de la roche COx, vis-à-vis des excavations et des charges de nature différente est exhaustivement étudiée. Les résultats obtenus au cours de ces programmes expérimentaux ont fourni d'énormes xiii connaissances sur le comportement à court et à long terme de l'argilite COx et l'impact de différents types de support sur le comportement des galeries. Du point de vue mécanique, divers essais de laboratoire, confirmés et complétés par les mesures in situ, ont montré un comportement complexe du COx à court et à long terme. Le comportement à court terme de cette roche hôte semble être prédominé par les mécanismes élastoplastiques et d'endommagement tandis que le fluage est plus répandu à long terme. Les résultats des essais de fluage au laboratoire mettent en évidence une phase de fluage transitoire dans la première étape, qui est suivie d'une deuxième phase asymptotique à un taux de déformation constant. En particulier, il est également indiqué que plus de 80% de la convergence des galeries observée a été atteinte au cours des 100 premiers jours et que les taux de convergence diminuent avec le temps. De plus, la surveillance des revêtements des galeries expérimentales confirme les effets significatifs du matériau compressible, qui peut être utilisé pour séparer le revêtement en béton de l'argilite COx (c'est-à-dire que le matériau agit comme revêtement extérieur dans le système de support de galerie). Ce matériau compressible peut être composé de coulis compressible ou d'une couche compressible à l'arrière du revêtement. La raison principale de l'utilisation de ce matériau compressible (avec une compressibilité de 30% à 80% du volume initial) est qu'il permet d'absorber la convergence dans le temps de la roche hôte et de réduire ainsi considérablement la contrainte radiale transmise au revêtement intérieur en béton. Les objectifs dans cette étude sont de quantifier l'effet de l'incertitude des propriétés du COx et de l'impact du matériau compressible sur la stabilité à long terme du revêtement en béton d'une galerie circulaire. Bien que des mécanismes de couplage complexes (tels que les modèles thermohydromécaniques couplés et le mécanisme d'endommagement) puissent régir la réponse en fonction du temps des ouvrages de stockage de déchets nucléaires, cette étude est limitée au seul le comportement mécanique de la roche. Conformément aux choix des galeries de stockage in situ, nous considérons la galerie de notre étude orientée parallèlement à la direction de la contrainte horizontale principale. Pour simuler le comportement mécanique de cette galerie, on adopte un modèle de déformation plane avec une galerie de section circulaire de rayon de 5.05 m et une dimension totale du modèle de de 55 m de chaque côté. Sur les frontières externes du modèle là contrainte in situ (isotrope dans ce cas et égale à 12.5 MPa) est imposée. Cette contrainte isotrope est prise de manière similaire aux contraintes principales suivant la direction verticale et celle de la contrainte mineure horizontale. La galerie considérée est supportée par un revêtement en béton C60/75 qui est séparé du massif par une couche compressible. Trois modèles constitutifs différents sont traités, qui sont respectivement un modèle viscoplastique pour la roche hôte, un modèle élastique non linéaire pour le revêtement compressible et un comportement élastique linéaire pour le revêtement en béton. Le modèle viscoplastique classique de Lemaitre est choisi pour caractériser le comportement à long terme de la roche hôte. Le choix de ce modèle simple est justifié par sa capacité à capter le comportement temporel des roches COx via un nombre limité de paramètres à calibrer, ce qui facilite considérablement le processus de quantification de l'incertitude à partir des essais de fluage au laboratoire. Le comportement du matériau compressible du revêtement extérieur est modélisé par un modèle élastique tri-linéaire avec le module élastique qui varie en fonction de la déformation volumétrique. Comme première tentative de quantification de l'incertitude des propriétés mécaniques à long terme de l'argilite du COx, dans cette étude, nous utilisons les données de sept essais de fluage triaxiaux correspondant à sept contraintes déviatoriques appliquées. L'analyse inverse pour déterminer les paramètres du modèle Lemaitre est réalisée en ajustant les données d'au moins deux tests de fluage sous deux contraintes déviatoriques différentes. Au total, l'analyse statistique des résultats obtenus est réalisée avec 88 échantillons générés à partir d'au moins deux essais. Les applications numériques ont été réalisées à la fois pour l'analyse déterministe et la fiabilité. En particulier, l'analyse de fiabilité basée sur le métamodèle de Krigeage vise à évaluer la probabilité de dépassement de la contrainte admissible dans le revêtement en béton de la galerie à 100 ans. Un sous-ensemble de quatre xiv échantillons enrichis à chaque itération au cours du processus itératif est choisi grâce à la possibilité d'utiliser des calculs parallèles dans Code Aster. Les résultats numériques montrent une diminution, comme prévu, de la probabilité de dépassement de la contrainte admissible en revêtement en béton lorsque cette valeur limite augmente. Pour une valeur fixée de compressibilité de la couche compressible (c'est-à-dire la valeur limite de déformation volumique entre le deuxième et la troisième branche élastique du modèle tri-linéaire) du revêtement compressible, la probabilité de dépassement de la contrainte admissible du revêtement en béton varie considérablement lorsque l'épaisseur du système de support varie. L'effet de la variation de l'épaisseur du revêtement intérieur est même plus remarquable. Le taux de déconfinement au moment de l'installation des revêtements peut également présenter un effet important sur la probabilité de dépassement des valeurs limites du support en béton. La stabilité à long terme du revêtement en béton est plus élevée, lorsque la compressibilité dans la première couche de support augmente. Suite aux résultats obtenus, on observe une faible probabilité de dépassement de la contrainte admissible dans le revêtement de la galerie en adoptant une épaisseur de 0.2(m) dans la couche compressible qui recouvre le support en béton final de 0.5(m) d'épaisseur. La probabilité est inférieure à 1% lorsque la contrainte admissible du béton est supérieure à 26 MPa. Ces résultats numériques confirment l'effet évident et l'avantage très significatif du matériau compressible pour assurer la stabilité de l'élément de support en béton de la galerie ainsi que démontrer la robustesse du dimensionnement actuel de l'alvéole considéré du projet Cigéo.

Le dernier chapitre de la thèse est dédiée à l'étude des incertitudes liées à la variabilité spatiale (incertitude aléatoire) des paramètres du comportement à long terme du COx. Ce type d'incertitude est généralement représenté par des champs aléatoires. Mathématiquement, ce champ aléatoire peut s'exprimer sous la forme d'une fonction de corrélation (telle que la fonction de Markovien ou de Gaussienne) avec un paramètre caractéristique essentiel, la longueur de corrélation spatiale. Pour la modélisation numérique (c'est-à-dire l'évaluation directe de la fonction de performance), ce champ aléatoire doit être discrétisé par une technique appropriée. Après la discrétisation de tous les champs aléatoires, l'évaluation probabiliste peut être effectuée comme un problème d'incertitude abordé dans les chapitres précédents. Néanmoins, l'augmentation significative du nombre de variables aléatoires après la discrétisation de ces champs aléatoires entraîne une difficulté prononcée dans la construction du métamodèle puisque la demande de calcul augmente de façon exponentielle à mesure que le nombre de variables augmente (problème appelé curse of dimensionality). La résolution du problème de grande dimension est toujours un sujet de recherche en cours dans la littérature et très souvent des approches de réduction de modèles sont entreprises. Un outil intéressant appelé représentation de modèle à grande dimension (High-Dimensional Model Representation HDMR) a reçu beaucoup d'attention au cours de ces deux dernières décennies. Dans cette méthode, une procédure de projection systématique entre les entrées et les sorties est prescrite pour révéler la hiérarchie des corrélations entre les variables d'entrée tandis qu'en réalité, pour chaque système physique, seules des corrélations d'ordre relativement faible des variables d'entrée affectent significativement la sortie. Cette propriété permet d'exprimer le HDMR en seulement quelques niveaux hiérarchiques pour représenter le système physique avec une précision satisfaisante. Chaque niveau hiérarchique de HDMR vise à appliquer un opérateur de projection approprié à la fonction de sortie. Entre les deux HDMR bien connus : l'ANOVA (Analyse de variance) -HDMR et le Cut-HDMR, ce dernier a été choisi dans les travaux du quatrième chapitre de cette thèse. L'effet de la variabilité spatiale des propriétés de l'argilite du COx sur la stabilité à long terme du support de la galerie a été étudié en utilisant la méthode AK-MCS modifiée ou la méthode Cut-HDMR en combinaison avec le métamodèle de Krigeage (notée ici AK-HDMR). Les hypothèses adoptées dans le chapitre précédent ont été retenue dans ces études, à ceci près que la variabilité spatiale des propriétés de la roche a été xv prise en compte, ce qui signifie que le massif environnant du tunnel est hétérogène. Chacun des quatre paramètres du modèle Lemaitre, décrivant le comportement dépendant du temps de l'argilite du COx, est décrit par un champ aléatoire éventuellement corrélé avec les champs aléatoires des autres paramètres Chaque champ aléatoire est caractérisé par une longueur de corrélation spatiale appropriée qui doit être établie par des caractérisation spatiale de ces paramètres. Néanmoins, dans le but de vérifier l'efficacité de la méthodologie développée et du fait que toutes les données concernant la longueur de corrélation des propriétés mécaniques de la roche COx ne sont disponibles, dans ces études nous considérons que tous les champs aléatoires ont la même valeur de longueur de corrélation. La méthode bien connue Expansion Optimal Linear Estimation (EOLE) a été adaptée pour discrétiser les champs aléatoires et étendue pour s'attaquer à leur relation de corrélation (c'est-à-dire la discrétisation des champs aléatoires intercorrélés). Cette procédure de discrétisation présente une étape nécessaire pour l'évaluation directe (c'est-à-dire la résolution numérique du problème déterministe) par le code-Aster. Les valeurs de l'ordre d'expansion tronqué (N) et la longueur de corrélation sont deux facteurs principaux pour l'évaluation et le contrôle de la précision des méthodes de discrétisation. Les résultats numériques montrent que l'ordre d'expansion tronqué (N) augmente rapidement lorsque la longueur de corrélation spatiale diminue. La méthode AK-MCS modifiée est d'abord choisie pour estimer la probabilité de dépassement de la contrainte admissible du revêtement en béton de la galerie dans le contexte de la variabilité spatiale de la roche hôte COx. Deux cas de variabilité isotrope et anisotrope dans l'espace des propriétés des roches COx avec le même ordre d'expansion tronqué N (et donc avec les mêmes valeurs de 100 variables aléatoires après la discrétisation) sont pris en considération. Ces études nous permettent de vérifier la possibilité de AK-MCS pour gérer le problème plus complexe, qui, pratiquement, peut représenter jusqu'à 100 variables, selon certaines études dans la littérature. Contrairement à ces dernières contributions dans le contexte des problèmes académiques, notre étude numérique montre que l'AK-MCS ne peut être raisonnablement appliqué que lorsque le nombre maximum de variables est d'environ 50. Au-delà de ce nombre, nous constatons une augmentation drastique des points d'entrainement du DoE et le temps nécessaire pour construire métamodèle de Krigeage. Il ressort de l'analyse de fiabilité que la diminution de la longueur de corrélation réduit la probabilité de dépassement des valeurs limites dans l'élément de support final de la galerie. En fixant la longueur de corrélation selon la direction horizontale et en utilisant la longueur de corrélation verticale plus petite, la probabilité de dépassement de la contrainte admissible du revêtement en béton diminue considérablement. Ces études mettent en évidence l'effet de variabilité spatiale des propriétés des roches COx sur la probabilité de dépassement du revêtement en béton de la galerie. En comparaison avec le problème d'incertitude, une probabilité plus faible a été remarquée dans ce problème de variabilité. Cependant, cette observation doit être renforcée en tenant compte les valeurs de la longueur de corrélation plus petites des propriétés des roches COx. Néanmoins, la limite de l'AK-MCS pour gérer le problème de grande dimension ne permet pas de mener ce type d'étude qui nécessite en fait une extension appropriée de la technique de métamodélisation AK-MCS. Dans une première tentative, la méthode de représentation du modèle à haute dimension (HDMR) a été considérée pour traiter ce problème. L'idée principale de la méthode HDMR est d'approcher une fonction de performance multi-variables par une combinaison des métamodèles de composants d'ordre zéro, premier, second et supérieur. Certaines études dans la littérature montrent que le HDMR peut être tronqué jusqu'au premier ordre (Cut-HDMR1) ou au second ordre (HDMR2) car les interactions de l'ordre supérieure à deux entre les variables d'entrée du HDMR sont négligeables. Pour dériver chaque métamodèle composant de HDMR, différentes méthodes peuvent être utilisées. Dans ce travail, le métamodèle AK-MCS, tel qu'introduit dans les chapitres précédents, est choisi pour construire chaque composant du Cut-HDMR, et nous notons la méthode adoptée comme AK-HDMR1 et AK-HDMR2 correspondant à chaque configuration de la méthode HDMR tronquée au premier ou au second ordre. Les études xvi numériques dans les deux cas de roches viscoélastiques et visoplastiques montrent que la probabilité de dépassement des valeurs limites dans le revêtement en béton prédite par l'AK-HDMR1 n'est pas exacte, et le résultat est largement amélioré par l'AK-HDMR2. Bien que le nombre de métamodèles construits augmente dans l'AK-HDMR2, cette méthode réduit considérablement le nombre de variables aléatoires et donc le temps nécessaire pour construire les métamodèles. Si la précision de cet AK-HDMR2 a été démontrée dans le cas d'un problème d'incertitude, son applicabilité dans le cas d'un problème de variabilité nécessite cependant plus d'améliorations.

Principales publications associées aux travaux

Articles dans des revues internationales à comité de relecture Do D. P., N. T.

Conférences internationales

Tran N. T., D. P. [START_REF] Do | Time-dependent reliability analysis of deep tunnel in the viscoelastic Burger rock with the sequential installation of liners[END_REF]. Time-dependent reliability analysis of lined tunnel within the linear viscoelastic Burgers rocks using the response surface method. In Recent advance in Computational Geotechnics, Nancy, France. http://coggus2.univlorraine.fr/scientific-sessions/.

LIST OF FIGURES

LIST OF TABLES

ABSTRACT

Design an underground structure in the rock mass is traditionally based on the deterministic (or semi-probabilistic) approach, in which the sources of uncertainty are addressed through the partial safety factors. However, the efficiency of this approach to derive an optimal and hence the construction cost has been considered as an open question. Recently, an advanced technique based on reliability analysis is used on the optimization of mechanical structure design from the assessment of the uncertainties propagation on their failure probability. So far, the application of such methods for designing deep tunnels has been only focused on the shortterm behavior of rock masses. In many cases, moreover, these methods lost their effectiveness when a strongly nonlinear behavior is considered. However, this is the case in some specific contexts, such as the deep nuclear waste disposal, where different nonlinear mechanisms (e.g., creep, damage, thermo-hydro-mechanical coupling, etc.) are significant sources of uncertainty on the long-term stability of structures.

This dissertation aims at evaluating the failure probability in the long-term of the concrete lining of deep tunnels by considering the uncertainty of the time-dependent behavior of the host rock. To this purpose, an extension/modification of the reliability analysis technique using the Kriging-based metamodel was carried out. The performance and the efficiency of this metamodel compared to the classical Monte-Carlo Simulation (MCS) method were firstly demonstrated in the context of a sequential construction of a deep tunnel in a linear viscoelastic rock. The metamodel is then applied to analyze the stability of the concrete lining of a gallery for a period of 100-years, which is typically the exploitation period of waste disposal galleries built in Callovo-Oxfordian (COx) claystone -the geological formation chosen for the implementation of the nuclear waste disposal in France if the project is licensed. The long term behavior of COx is described by the viscoplastic model of Lemaitre. The assessment of uncertainties and cross-correlation of Lemaitre's model is carried out using the raw database from creep tests performed by the laboratory test on the sample scale. The reliability analysis exhibited that the long-term stability of concrete inner liner has been affected by the uncertainty of the host rock properties and the great benefit of the compressible layer. According to these studies, the stability can be determined as a function of the compressibility and/or the thickness of the compressible layer. Moreover, the spatial variability of the mechanical properties of the host rock was also undertaken. To deal with the high dimensional problem associated with the discretization of the random fields, an adaptation of the Cut-HDMR method combined with the Kriging-based metamodeling technique was proposed in this study. The applicability of such a method was then tested through numerical applications.

RÉSUMÉ

GENERAL INTRODUCTION Context

In the science and engineering of everyday practice, the mitigation of hazards and risks is of capital importance. Herein there is an increasing interest in developing strategies and methods for characterization and modeling of uncertainties propagation, variability, and random systems. In particular, in the engineering of structures design and on the estimation of the life expectancy of structures, the reliability analysis has become powerful tools for the optimization design process and structure maintenance. For human-made materials currently used in the industry, many efforts are made to assure full control of the final properties, and the uncertainties on the properties of a produced material are given at the same time as the properties themselves. One of the ways these methods are used in practice is through the design rules/codes such as Eurocodes, ASTM, and so on. For example, in various design approaches promoted by Eurocodes, the partial safety factors method is applied to take into consideration the uncertainties and variability of loads and material properties.

However, such methods are not imposed on underground construction, despite their application in other various fields of geotechnics (Eurocodes 7, European geotechnical code of design explicitly indicates the domain of underground designed to be out of its application fields). This is explained, by large, from the slightly different situation of underground excavations as compared with other types of structures. The soils, and more often, the rock masses that host these excavations, in contrast to human-made materials, are natural, heterogeneous geological materials produced by depositional and post-depositional processes, containing tectonism, diagenesis as well as sedimentation [START_REF] Lü | Probabilistic assessment of tunnel convergence considering spatial variability in rock mass properties using interpolated autocorrelation and response surface method[END_REF]. Their random evolutive internal structure is described by a combination of constitutive minerals that govern local mechanical behavior at any point of massifs. This spatial organization of minerals (and for that reason, of properties) makes the characterization of the properties of geomaterials more complex since both epistemic and aleatoric uncertainty should be qualified. Using the terminology used by [START_REF] Sudret | Uncertainty propagation and sensitivity analysis in mechanical models-Contributions to structural reliability and stochastic spectral methods[END_REF] the aleatoric uncertainty makes reference to natural variability in the phenomenon under consideration and epistemic uncertainty to a lack of knowledge. Aleatory uncertainty can be understood as input parameters for geotechnical design, for example, in situ stress field, rock mass strength parameters and modulus of deformation, so forth. Yet, they have never known precisely before. The lack of knowledge is limited to uncertainties in the measurements of the physical properties of the rock masses and lack data to understand and characterize the phenomenon, among other things entirely. Because the physical limitations of a number of measurements could be taken in the region concerned, besides, the results of deterministic analyses can miss the exact failure mechanisms and ignore the precise response that might pass through the weakest part of the rock in the sense of its randomness properties. Moreover, samples from laboratory or field tests are more or less disturbed; yet errors are often introduced in the measurement and estimation of the strength and deformation characteristic of the rocks.

While all these particular features of geomaterials and rock masses explain by themselves the need for specific approaches to mitigate uncertainties or to verify the robustness of underground design, they become more evident when the underground waste storage is in question. Since deep underground radioactive waste disposal has been a preferential option in many countries, the questions of assessment of uncertainties linked with such structures rise with increasing attention. Besides the multiphysics coupling, the long-term behavior of rock masses and the uncertainties linked with the behavior of structures in such conditions are perhaps the central points on the design of storage facilities. Such sources of uncertainties represent a challenge for researchers and engineers. On the one hand, spatial (e.g., several tens of square kilometers) and temporal scale (e.g., some million years) of studies on the storage field are with no comparison with any up-to-now built-up underground construction. On the other hand, the complexity of interaction many phenomena of various nature (mechanical, chemical, hydraulic, thermal, radioactive) and kinetics are never seen before in a single problem. And yet, if one considers just one physical problem (let say mechanical one), it is ramified in other problems and mechanisms that could be active in some moments of the life of underground excavation but constraint behavior of the structure on the next phases of the underground structure. For example, the irreversible instantaneous behavior of surrounding rocks will control the Excavation Damage Zone's initiation and growth around a cavity. This, in turn, will impact the stress and strain distribution as well as the hydraulic field around the cavity to impact the long-term behavior of the excavation. It is intuitively clear that any uncertainty on the irreversible behavior of rock at any moment of the life of the excavation will be propagated to uncertainties in the future phases.

The work of this thesis is a part of the number of studies performed by the scientific community relating to the Cigéo project. This is deep geological disposal for high-level radioactive waste (HLW), and intermediate level long-life radioactive waste (IL-LLW) called Cigéo is planned to be constructed in deep Callovo-Oxfordian claystone (COx) formation if licensed.

We have not to ambition here to deal with all sources of uncertainties in such a big project and even not all those linked to all kinds of underground excavations. Our principal interest in this work is the propagation of uncertainties of long-term mechanical behavior of rock masses on the response of a typical underground excavation.

Because of the above-mentioned reasons being inherent of any underground constructions on the rocks and in long-term analysis, the problem is not straightforward. Certainly, there is a large number of results from laboratory and in situ tests performed on rock masses from Callovo-Oxfordian claystone (COx) [START_REF] Armand | EDZ Characterisation in the Meuse / Haute-Marne Underground Research Laboratory[END_REF]Zhang et al., 2013;[START_REF] Jia | Numerical modeling of the influence of steel corrosion on the thermo-hydro-mechanical response of the Callovo-Oxfordian claystone[END_REF]. These published papers revealed an uncertainty in COx of instantaneous mechanical properties values along with the depth and some random variation with extensions. The same variability is also observed in the long-term behavior of rock that seems to be activated at low differential stress [START_REF] Billig | DeCo Grout -Innovative Grout to cope with rock deformations in TBM tunnelling[END_REF][START_REF] Jia | Hydromechanical modelling of shaft excavation in Meuse/Haute-Marne laboratory[END_REF][START_REF] Martin | Geotechnical feedback of the new galleries excavation at the ANDRA underground research laboratory-Bure (France)[END_REF][START_REF] Mohajerani | Etude expérimentale du comportement thermo-hydromécanique de l'argilite du Callovo-Oxfordien[END_REF][START_REF] Armand | Geometry and properties of the excavation-induced fractures at the meuse/haute-marne URL drifts[END_REF]. These variations of mechanical properties follow the variation of the mineralogical composition and structural anisotropy. For example, the seismic surveys on the Cigeo projects of Andra also show vertical and lateral variations in the properties of COx claystone [START_REF] Armand | Geometry and properties of the excavation-induced fractures at the meuse/haute-marne URL drifts[END_REF]. In particular, these uncertainties impact the initiation and extension of the Excavation Damage Zone [START_REF] Schuëller | A critical appraisal of reliability estimation procedures for high dimensions[END_REF].

For the probabilistic approaches of deep underground tunnel structures, one could cite, among others [START_REF] Pan | Probabilistic evaluation of tunnel face stability in spatially random soils using sparse polynomial chaos expansion with global sensitivity analysis[END_REF][START_REF] Shokri | A review study about spatial correlation measurement in rock mass[END_REF][START_REF] Zhang | Tunneling and deep excavations in spatially variable soil and rock masses: a short review[END_REF] for the tunnel stability. When data on the quality of the rock mass is sufficient, a probabilistic approach allows us either to understand the risks of the project better or to have a more efficient geo-mechanical zoning and a more reliable cost estimate. However, probabilistic calculation and design tools are not yet widely used for tunnel lining [START_REF] Oreste | A probabilistic design approach for tunnel supports[END_REF].

The particularity of the work presented here, compared with reliability analyses generally performed on structural mechanics, resides principally in considering the time-dependent behavior of the rocks leading with probability failure evolving with the time.

After an analysis of alternative approaches, on the one hand, that could be used for the assessment of the uncertainties related to the mechanical properties of COx and their spatial variability; on the other hand, their impacts on the response of a typical underground storage structure, this dissertation proposes an optimized methodology for reliability analyses. This approach is firstly used for the evaluation of the impact of the uncertainties on the mechanical parameters of COx. Then the innovative consideration of the compressible material as the outer liner on the probability failure of the inner liner of a gallery is also investigated. Later the impact of the variability in the space of the rock properties on the long-term failure probability is discussed. For the spatial variability, both the uncorrelated and correlated Gaussian random fields are addressed with nonlinear structural problems. Remarkably, the treatment of uncertainty methodology contained in this study can be followed to treat uncertainty in many fields, such as analysis of natural hazards (e.g., wind, earthquakes, floods) and decision making in environmental risks. We also use the Code-Aster based on Python scripts for studying material heterogeneity models.

Research objectives

The thesis's principal objective is to establish a methodology for estimation in terms of probability to failure of the inner liners, the impact of the uncertainties of mechanical behavior of rock masses on the stability of underground storage excavations. Once established, the methodology is used to evaluate this impact by considering either the epistemic uncertainties linked to the long-term behavior of surrounding rock masses or those linked with these parameters' spatial variability. For that, following auxiliary problems with increasing complexity and a step-by-step research approach is followed consisting to:

(1) Demonstrate and compare diverse reliability analysis methods applied to benchmark tunneling problems. Especially, methods for selecting samples from the meta-model could be applied for the real project to tackle the uncertainty and the inherent spatial soil variability. It is necessary to propose a new approach to deal with these problems base on the viscoelastic/viscoplastic behavior of rock mass.

(2) Provide a comprehensive literature review of state of the art on deep geological disposal of radioactive waste, emphasizing the geomechanical description of COx claystone in terms of rock parameters, mechanical material behavior, and tunnel construction.

(3) Quantify the uncertainty in mechanical properties of COx from the laboratory tests (in a project of Andra) to model a deep tunnel in rocks by considering the uncertainty and variability, and notably regarding time-dependent problems in periods such a hundred years.

(4) Estimate the probability of failure and incorporate the inherent variability in material input parameters into evaluating the stability of tunnels (e.g., liners or tunneling supports). It can be understood as calculating the probability of failure if a prescribed threshold of system response is exceeded. In other words, the computation is based on limited area functions with respect to the uncertainty of the parameter in a quantitative and verifiable manner.

(5) Investigate the tri-linear model adaptation in Aster Code for the compressible material coating used in the project. Then, to perform two-dimensional numerical analyses on the COx rock and support liners responses.

(6) Help the readers understand and apply mathematical frameworks of uncertainty and spatial variability in geotechnical engineering.

The most important objective here is to propose a practical approach to perform a reliabilitybased analysis of the tunnel; yet, the approach may be used in other fields. The nonperformance probability of exceedance of the allowable support strength will be investigated with two main cases: uncertainty problems (as in Chapters 2 and 3) and variability problems (as in Chapter 4).

The layout of the thesis

The dissertation is organized according to the subsequent chapters:

Chapter 1 presents the state-of-the-art of the existing popular models to simulate the viscoelastic and viscoplastic behavior of rock masses by considering uncertainty and variability, which can be divided into four main parts. The first one is a quick review of some basic concepts of reliability analysis and techniques used to solve structure reliability problems.

In which several keywords of the study will be mentioned, namely: uncertainty, variability, reliability analysis, random variables, random fields, performance function, time-dependent reliability analysis. Secondly, some popular reliability analysis methods will be revisited. Thirdly, we will show the state-of-the-art reliability analysis of the underground structure.

Then, the long-term massif behavior, especially viscoplasticity of the rock mass, is also briefly reviewed.

Chapter 2 illustrates the time-dependent reliability analysis of a deep double-lined tunnel in linear visco-elastic rocks. The analytical method will be used to compare and validate the closed-form solution with some suitable reliability methods. That can effectively use in the next sections for the real project and compare with numerical tests in the ideal cases (e.g., incompressible rocks, isotropic and homogenous materials).

Chapter 3 precisely shows the emphasis of this study: Application of tunneling analyses via considering the uncertainty in time-dependent, viscoplastic rocks into a specific project of Andra The method for simulating in 2D the digging of a tunnel that forces on the uncertainty of the geo-mechanical parameters. The expected results will be some failure probabilities of the main-liner corresponding to the prescribed limit states and to considered time (one hundred years). Some parametric case studies are also addressed for obtaining the general and exact behavior of the support system (i.e., the concrete liner).

Chapter 4 characterizes the spatial variability of four input parameters of COx rocks within nonlinear viscoelastic problems. The assessment of stochastic rock characteristics, including spatially correlated properties, is deeply considered. Using similar settings from previous chapters, we can discretize the random fields with expansion optimal linear estimation method (EOLE) and high-dimensional model representation (HDMR). The case studies are different from the ones in Chapter 3 at the evaluation of the spatial variability. This chapter also mentions mathematical frameworks to quantify the variability and the correlation length of rock properties at different scales.

At the end of the thesis, we give conclusions and perspectives.

CHAPTER 1. BACKGROUND OF RELIABILITY ANALYSIS AND ITS APPLICATION FOR UNDERGROUND STRUCTURES

Introduction

This chapter presents the background and state-of-the-art application of the reliability analysis for geotechnical engineering, especially the construction of underground structures, by considering the uncertainty and variability of its surroundings' mechanical properties (rock masses). Firstly, a quick review of some basic reliability analysis concepts, in which the essential keywords (namely: uncertainty, reliability analysis, variability) will be mentioned.

Next, some commonly used reliability analysis methods will be highlighted. Then, we present a quick review of the state-of-the-art reliability analysis of the underground structure. The last part of this chapter aims to synthesize some well-known models that characterize the long term behavior, focusing on the viscoplasticity of rock mass.

Uncertainty concepts

Uncertainties are inherent, unavoidable in the natural geomaterials (rocks, soils), and must be addressed when designing any construction [START_REF] Gulvanesian | Annex C-Calibration Procedure, Leonardo DaVinci Pilot[END_REF]. In geotechnical and geological engineering, uncertainty plays a vital role as a central feature [START_REF] Christian | Geotechnical Engineering Reliability: How Well Do We Know What We Are Doing[END_REF]; thus, it is necessary to extend the deterministic analysis to stochastic (or non-deterministic, or probabilistic) one. Uncertainty in geotechnical response prediction is a complex phenomenon resulting from many different sources, which are typically classified as epistemic or aleatoric [START_REF] Lacasse | Model uncertainty in pile axial capacity calculations[END_REF]. Figure 1-1 demonstrates the types of uncertainty in geotechnical soil properties.

In geotechnics, there are two types of epistemic uncertainties: measurement error and transformation uncertainty. The first type is associated with the scarcity of data (quality, quantity, specification) coming from different sampling stages, preparation, transport, and testing procedure, which is mainly explained by the lack of knowledge or data on the resulting observed phenomenon. This uncertainty can be reduced to a minimum value by considering a larger number of samples, in which the minimum value depends on the systematic errors involved in the whole process. While using empirical or other correlation models when transforming data from laboratories or in situ measurements can be another type of epistemic uncertainty, so-called transformation uncertainty. It is possible to significantly reduce the characteristics of uncertainty by using more refined empirical or mathematical models for the latter type [START_REF] Sanchez-Silva | Physical modeling of soil uncertainty[END_REF]. On the other hand, uncertainties can be classified as aleatoric if the modeler does not predict the possibility of lowering them (Kiureghian and Ditlevsen, 2009). This type is also well-

Uncertainty in soil properties

Epistemic Aleatoric

Measurement error

Transformation uncertainty

Spatial variability

Temporal variability known as inherent uncertainty. The soil variability that mostly results from natural geological processes modifying soil mass in situ is aleatoric (i.e., linked to the inherent randomness). This intrinsic variability of rock is based on the complex history of formation and the continual reworking of geological materials linked to the impossibility of a deterministic description of the phenomenon's intrinsic variability under consideration. The aleatoric uncertainty can be expressed by spatial variability in space for phenomena occurring in different locations but at a single point in time, or by temporal variability over time for phenomena occurring in a single location, or both in time and space.

Spatial variability concepts

In this subsection, some popular theoretical concepts related to spatial variability are presented. Others will be mentioned in chapter 4 (e.g., a method for discretizing random fields, error measures for random field discretization, correlated random fields, and so on).

Spatial variability:

The inherent spatial variability of geomaterials could be simply understood as the variation of rock (soil) properties from one point to another in space [START_REF] Sekhavatian | Application of random set method in a deep excavation: based on a case study in Tehran cemented alluvium[END_REF]. In the reliability analysis, the spatial variability of geomaterials properties is defined through the random field concept.

Random field: Random field theory is a mathematical concept that allows us to model the uncertainties of physical properties in a continuous medium regarding either spatial variability or randomness. In this context, a random field H(x, ) is a collection of random variables indexed by continuous parameters   

x d    ( d    is multidimensional domain).
Here, x also denotes the geographical location (i.e., the coordinates of the system). Note that for a given

x0 ∈ D, H(x0, ) is a random variable. Contrarywise, for a provided outcome  0∈, H(x,  0)
is a realization of the field (the output is a trajectory, i.e., a curve).

The random field is simply a stochastic process that usually takes on values in a Euclidean space and has a parameter space with a dimensionality of at least 1. In geotechnical mechanics, a random field can be in time or distance. However, in most cases in geotechnical mechanics, as well as in our study, random fields are considered in the space; in other words, the spatial variability is modeled by random field theory. Notably, probability distributions and the correlation at adjacent zones can completely describe rock parameters as random variables at any location in a rock layer [START_REF] Lü | Probabilistic assessment of tunnel convergence considering spatial variability in rock mass properties using interpolated autocorrelation and response surface method[END_REF].

Covariance (Coefficients of variation): This is a relatively stable measure of variability, and an essential statistical tool evaluates how much the mean values of two variables change together [START_REF] Cai | Rock mass characterization and rock property variability considerations for tunnel and cavern design[END_REF]. The sign of its value shows they tend to vary together (positive sign) or inversely move (negative one). Let Cov[X, Y] be covariance of two random variables X and Y.

It is analogous to the variance and reads:

[ , ] [( )( )] E( ) E( ) E( )        X Y Cov X Y E X Y XY X Y (1.1)
From the known data, we can use the familiar equation as follow:

1 1 [ , ] ( )( )        n i X i Y i Cov X Y X Y n (1.2)
where X, Y are the mean values of X, Y, respectively; n: number of periods.

In the random fields, these equations can be rewritten as:

   ( ', *) [ ( '), ( *)] ( ') (t') ( *) (t *)            X X X X C t t Cov x t x t t t (1.3)   E ( ') ( *) (t') (t *)     X X X X t t
where X(t'), X(t*) stand for the mean of vector X at the different positions t', t*, respectively. Some examples of covariance for geotechnical properties and in situ tests can be from some authors (e.g., [START_REF] Briaud | Coefficient of variation of in situ tests in sand[END_REF][START_REF] Guharay | Geotechnical Engineering: Reliability Aspects[END_REF]). Readers can find some types of covariance functions in the Appendix.

Correlation coefficient: The random variables could be considered as independent ones, such as in most case studies in literature or correlated together, or both (i.e., some are correlated, others are none). When pair variates X and Y are said correlated, the likelihood of a certain value of Y must depend on X's value. For example, for geotechnical materials, the friction angle and cohesion are often seen as correlated. The main aim of the correlation coefficient , X Y  is to present a dimensionless measure of the degree of correlation between random variables.

, [ , ] X

Y X Y Cov X Y     (1.4)
where X, Y are the standard deviations of X, Y, respectively.

Similarly, in the random field, this equation becomes as:

  ', * ( ', *) ( ') ( *) X X C t t t t t t     (1.5)
In which X(t) is the standard deviation of X at the position t.

Correlation function (or autocorrelation coefficient function):

indicates the spatial correlation of the relevant design material properties (e.g., Young modulus) at distinct places [START_REF] Tian | Bayesian identification of random field model using indirect test data[END_REF].

Correlation length: (or scale of fluctuation) is a characteristic parameter acting in the definition of the autocorrelation coefficient function. The (spatial) correlation length  (or l as in Appendix A1-3), indicates the fluctuation range of the parameter, i.e., it illustrates the variability of material properties in horizontal-/vertical-axes [START_REF] Shen | Non-deterministic analysis of slope stability based on numerical simulation[END_REF]. In the context of random field theory, this length plays the most vital role in measuring variations in space. It can be separated into two points by a distance greater than  are mainly uncorrelated. The smaller the correlation length, the faster fluctuation of soil parameter is; otherwise, the parameter is relatively consistent over a wide range. Mathematically, θ is defined here as the area under the correlation function [START_REF] Vanmarcke | Random fields: Analysis and synthesis[END_REF].

0 ( ) 2 ( )               d d (1.6)
Note that not all correlation functions will satisfy the criterion (Eq. 1.6), and for such random processes, θ =∞. An example of a process with infinite correlation length is a fractal process. Moreover, the correlation length is actually only meaningful for strictly non-negative correlation functions [START_REF] Fenton | Risk assessment in geotechnical engineering[END_REF]. This length can be determined from either experiment, design codes, or other relevant projects.

The consideration of spatial variability in the reliability analysis can be conducted by using a discretization technique that transforms the random field (i.e., infinite dimension problem) into the random variables (i.e., finite dimension problem). This problem can be solved similarly as in the previous case, i.e., the uncertainty problem. Explicitly, an application of an Expansion

Optimal Linear Estimation (EOLE), a popular discretization technique for random fields, will be presented in the last chapter of this thesis to deal with correlated random fields.

Reliability analysis concepts

The simplest understanding of the reliability concept could be a probability of a required function of a system. In other words, the reliability of a system can be understood as the probability of non-failing in meeting demands for a predefined time interval [START_REF] Mahmoudi | Probabilistic analysis of a rock salt cavern with application to energy storage systems[END_REF]. Indeed, one has been developing many different reliability methods to take the uncertainties into account in real engineering structural properties. For example, [START_REF] Baecher | Reliability and statistics in geotechnical engineering[END_REF] showed several primary methodologies to deal with practical reliability analyses. Recently, many researchers, such as [START_REF] Simpson | Reliability in geotechnical design-some fundamentals[END_REF][START_REF] Marelli | UQLab: A framework for uncertainty quantification in Matlab[END_REF][START_REF] Day | Practical Geotechnical Design using Reliability-Based Design Methods[END_REF], found numerous robust tools for the computation of reliability in geotechnical design.

In reliability problems, the uncertain input variables, called random variables, and their space is separated into the failure and safe domains by the so-called limit state function LSF (or limit state surface, LSS). The limit state function g(x) defines a failure criterion that is a function of all random variables x. Note that x could be independent as well as dependent variables.

Besides, x could be identically distributed samples generated according to the probability density function (PDF) in sampling methods. The concept of LSF is generic and can be transferred to a large family of problems. Currently, among various performance functions exist, the well-known ones in the context of underground structures are usually related to tunnel convergence, the support capacity, plastic radius, as brief in Table 1-1 below. The quantification of uncertainty propagation to assess the probability of failure (or exceedance probability, PF) is the main goal of structural reliability analysis. Theoretically, the probability is computed by integrating the joint probability density functions (PDF) of the random variables over the failure region as follows: In many complex structures, the uncertainty of the system response or system performance can be estimated by using specific techniques that can approximate the actual probability of failure value. In this way, many scientific researchers investigated and proposed various reliability analysis methods within different approaches. However, they could be categorized into three main methods according to some recent overviews (Hu and Mahadevan, 2016a;Huang et al., 2016Huang et al., , 2017)), which are captured in Figure 1-3 below. This classification will be detailed in the next sections, especially the well-known Kriging-based method, which is chosen as the principal method throughout this thesis. In many contexts related to geotechnical engineering, the time-dependent behavior of structure induces the variation of exceedance probability as a function of time. That means, in these contexts, the reliability analysis must be addressed to determine the time-dependent failure probability. Accordingly, the first-pass method's essential role in assessing the reliability of the structure under time-variant operating conditions has attracted much attention. [START_REF] Zhang | Time-Dependent Reliability Analysis Through Response Surface Method[END_REF]. The LSF is now written in the form g(x,t), which is a function of x (i.e., a random variable vector) and t (time instance).

( ) 0 ( ( ) 0) ( )      x x x f g P P g f dx (1.
The time-dependent failure probability (also called the cumulative failure probability) (Pf) with the time of interest [0, T] reads:

    (0, ) g( , ) 0, 0,     x f P T P t t T (1.8)
Due to the presence of the time factor, the computation of the probability as in the left-hand term in Eq. (1.8) is often known as a difficult task. The difficulty mainly lies in the monotonic /nonmonotonic behavior of the structure, which includes a time-dependence. For the nonmonotonic case, the maximum failure probability is not easy to find since it often requires either considering the full loading-history or applying the total probability theorem. This case, therefore, has been considerably increased the cost of the calculation. Fortunately, in most cases of underground structures, the structural responses (such as stresses, deformations as well as displacements) mostly vary in monotonic trends. That makes the time-dependent reliability analysis can be performed only at the time-interval border, and the problem becomes the timeinvariant reliability.

Different reliability analysis methods

Local reliability methods

The local reliability methods refer to using the local approximate of actual LSF in computing the probability of failure. These methods often use the norm of "design point (DP)", "most probable failure point," and the "reliability index". For example, the first-, and second-order reliability method (FORM and SORM) are famous in analytical approaches to estimate the LSF around the design point. These methods require handling the gradient of the LSF; thus, it is not recommended to apply for high-nonlinear and non-differentiable functions.

Mean Value First-Order Second Moment (MVFOSM) Method

This method is often called as mean value (MV) method or First-Order Second Moment method (FOSM) originally found by [START_REF] Cornell | A probability-based structural code[END_REF], which is the simplest and least expensive reliability method formed on the approximation of the first-order Taylor series (Huang et al., 2017). The mean  g and variance 2  g of the performance function, reliability index  , and the failure probability are calculated by: ( )

   g x g (1.9) 2 1 1 ( ) ( ) ( , )            n n x x g i j i j i j g g Cov x x x x (1.10) G G     (1.11) ( )     f P (1.12)
Herein, x  stands for the mean of vector x, Cov(xi, xj) is the covariance of xi, xj, and ( )   denotes the standard normal cumulative density function (CDF). The method is satisfactorily exact only when g(x) is almost linear, and variable x is approximately Gaussian. Moreover, it often fails to consistent with various equivalent formulations of the identical problem (Huang et al., 2017).

1.2.1.2. First Order Reliability Method (FORM)
The essence of FORM is to linearly approximate the response function g(x) at the design point (denotes x*). This point is on the limit-state surface (LSF) g(x)=0 and locates nearest to the origin, as indicated in Figure 1 -4 (Silva, 2007;[START_REF] Langford | Application of Reliability Methods To the Design of Underground Structures[END_REF]. One can consider the linear Taylor series expansion around the design point as:

  ( ) ( *) ( *) *      x x x x x T g g g (1.13)
where g(x*)=0 because x* is on the LSF.

Commonly, the gradient vector is replaced by its negative and normalized version, called the alpha-vector:

( ) ( )      x x g g (1.14)
Combination of these two equations, notably all vectors are column vectors, yields:

    ( ) ( *) * ( *) *              x x x x x x x T T T g g g (1.15)
One often rotates and/or transforms the random variable vector from the physical space (xspace or original space) into the standard normal space (U-space). Then, the design point (DP) can be obtained by using the two conditions of the optimization problem: minimize ||U|| and must be subjected to g(U)=0. Once obtaining the DP, the distance from the DP to the origin is denoted as Hasofer-Lind (HL) reliability index HL, that is, HL‖U*‖. Low and Tang (Low and Tang, 1997a;1997b;[START_REF] Christian | Geotechnical Engineering Reliability: How Well Do We Know What We Are Doing[END_REF]2007) illustrated an alternative interpretation of the HL index and the FORM on the basis of an increase in the size of the ellipsoid in the physical space without any rotations or transformations. For the case of correlated or uncorrelated normal variables, the formulation of the HL index HL reads:

    1 min       x x μ C x μ T HL F (1.16) Or, equivalently 1 min                       x R T i i i i HL F i i x x (1.17)
where  are the mean values and standard deviation, respectively. C and R denote the covariance matrix and correlation matrix, respectively. F stands for the failure domain.

Thus, the FORM is capable of approximating the LSF linearly and provides the failure probability as follows:

)

   FORM f HL P (- (1.18)
where Φ(•) denotes the standard normal CDF. Note that U* is not known at the outset and has to be found by iteration. Particularly, when the origin of U-space is already located in the failure zone (F), the failure probability is computed as: Notably, to transforms the random variable vector from other type distribution functions to equivalent standardized normal distribution (U-space), one often applies either Rosenblatt or Nataf transforms [START_REF] Schulze | Reliability as a Public[END_REF]. Besides FORM, one can find a similar method, so-called inverse-FORM, where the different critical points are to argmin{g(U)} and subject to =0 (or Pf=Pf0); one can refer to [START_REF] Youn | Hybrid Analysis Method for Reliability-Based Design Optimization[END_REF][START_REF] Sudret | Uncertainty propagation and sensitivity analysis in mechanical models-Contributions to structural reliability and stochastic spectral methods[END_REF] for more details.

1 ) )        FORM f HL HL P (- ( (1.19)

Second-Order Reliability Method (SORM)

The SORM approximates the failure function by its quadratic form. It provides a more accurate method than the FORM because the expansion of the non-linear LSF also involves the secondorder derivatives of g(x). SORM analysis also needs both the reliability index  and the design point as inputs. That means prior FORM analysis must be required [START_REF] Low | FORM, SORM, and spatial modeling in geotechnical engineering[END_REF].

Because response functions cannot be explicitly expressed in most cases of geotechnical applications, the partial derivatives must be computed numerically utilizing central differences.

The linear estimation of the partial differential coefficient is found by changing each variable by a small amount, while the remaining variables are held at their mean values [START_REF] Langford | Application of Reliability Methods To the Design of Underground Structures[END_REF].

Notably, FORM and SORM methods strongly depend on the mean values of the input variables, where the partial differential coefficients of the safety margin are assessed (invariance problem). Hence these methods are exact only for linear functions, and the error of the approximation could be significant for applying the high nonlinearity of response functions [START_REF] Shen | Non-deterministic analysis of slope stability based on numerical simulation[END_REF].

1.2.1.4. Response Surface Method (RSM)
The RSM is a functional relationship between a dependent output and many related input variables. When the performance function is often implicitly expressed in complicated geotechnical problems, the RSM can be used as a simple yet effective and accurate technique to solve various widespread reliability problems. Box and Wilson (1954) first described this method; then, its extensive applications for reliability analysis can be listed in civil engineering fields as (Low and Tang, 1997a;[START_REF] Low | Reliability analysis using object-oriented constrained optimization[END_REF]Mollon et al., 2009aMollon et al., , 2009b;;Lü and Low, 2011a;[START_REF] Shen | Non-deterministic analysis of slope stability based on numerical simulation[END_REF][START_REF] Low | FORM, SORM, and spatial modeling in geotechnical engineering[END_REF][START_REF] Hamrouni | Reliability analysis of shallow tunnels using the response surface methodology[END_REF]. The fundamental concept of RSM aims to fit the actual LSF g(x) by a closed-form polynomial function,  ( ) x g , by using both deterministic analyses and an iterative algorithm (Lü and Low, 2011a). Among types of functions, the quadratic polynomial function is usually utilized within U-space as:

 0 1 1 1 ( ) ( ) ( )             U U V U a n n n i i ij i j i i j g g a a u a u u (1.20)
where U is the vector of basic random variables ui in U-space and a=[a0, ai, aii, aij] T is the coefficient vector being determined by using the sampling points (DoE). The matrix

V(U) represents [1, ui, ui 2 , uiuj].
A procedure of using RSM for reliability analysis with five main steps in the proposal by (Lü and Low, 2011a). ( 1 Note that, for the correlated variables (with a correlation matrix R), one can use the Cholesky factorization to decompose the matrix R into a triangular matrix L. Let n = LU, the reliability index expression reads, base on Low and Tang algorithm [START_REF] Low | Efficient spreadsheet algorithm for first-order reliability method[END_REF], as follow:

1 1 min min ( ) min R LL U U          T T T T x F x F x F n n n n (1.21)
Once obtaining vector n by performing the five mentioned steps, the original vector x (in Xspace) of d random variables xi (i=1, 2 ..., d) can be obtained from the inverse distribution equation as follow:

  1 ( ) i i x F n    (1.22)
Although the RSM can be seen as one of the most efficient approaches [START_REF] Jiang | Efficient slope reliability analysis at low-probability levels in spatially variable soils[END_REF]Li et al., 2016), it needs a prior assumption on the order and type of polynomials that would be misleading if the true performance response is a multimodal function with several peaks and troughs [START_REF] Liu | Efficient system reliability analysis of soil slopes using multivariate adaptive regression splines-based Monte Carlo simulation[END_REF]. To overcome the assumption, one often uses other methods, such as surrogate models.

Point Estimate Method (PEM)

The PEM proposed by [START_REF] Rosenblueth | Two-point estimates in probabilities[END_REF]) is a numerical procedure that uses Gaussian quadrature to approximate the fundamental statistical moments of design parameters based on other parameters. The basic concept of PEM is to use discrete equivalent distributions (or probability mass functions (PMFs) consisting of a few discrete points) to replace the continuous probability distributions (PDFs) for each input variable. Next, one assesses at each point to estimate the values of the mean and standard deviation of the material response parameter. By using only the first two fundamental statistical moments, the method can be applied without clarifying the form of the distribution [START_REF] Langford | Reliability based approach to tunnel lining design using a modified point estimate method[END_REF].

Notably, this method can be easily applied using the Excavation & Support Design Software, so-called Phase2, a 2D Elastoplastic finite element stress analysis program for underground or surface excavations in rock soil.

1.2.2. Sampling reliability methods

Monte Carlo Simulation (MCS)

MCS (also called "direct MCS", "crude MCS", or "classical MCS") is a popular way to estimate the exceedance probability based on the number of the negative value of performance function over the total number of trials. Unlike the local reliability methods, MCS is an utterly robust and universal method, so it could be used to verify other reliability methods for estimating the failure probability of structures [START_REF] Arab | A fast and robust method for estimating the failure probability of structures[END_REF]. This simulation has many unique properties that make both positive effects and also the opposite. In other words, based on the capability of computing the probability for complex and highly non-linear systems, it is a widely used method that commonly applies to all the benchmark problems. Despite being robust and accurate, it is usually time-consuming because of the requirement of a large number of calls for the computation. In particular, it takes a paramount time when using it with finite element/finite difference methods.

With the MCS method, the failure probability is computed numerically by simulating a sequence of NMC independent realizations:

    ( ) ( ) ( ) ( ) 1 1 ( ) 0 1 ( ) , ( ) 0 ( ) 0          x 1 x 1 x x MC i N i i f i i MC if g P g g N if g (1.23)
where NMC counts the random samples being generated for estimation of probability by the MCS.

Note that, besides the strong capability of this method, the minimum number of realization NMC is found as follows:

2 / 2 2 (1 ) ˆf M C f z P N P     (1.24)
where z/2 is a quantile of the normal distribution, 1- is the confidence level (e.g., 95%),  stands for the relative error of the estimated probability of failure, ˆf P , (e.g., 0.10%). In another approach, the coefficient of variation (Cov) of the estimator failure probability is applied, instead of using Eq. (1.24):

1 ( ) ˆ. f f f MC P Cov P P N   (1.25)
With a predefined Cov (says 0.05), one can identify NMC from Eq. (1.25), then computing the failure probability using Eq. (1.23) according to the acceptable level, according to Cov. The Cov can be used as a convergence criterion to adaptively increase the Monte Carlo sample size until some desired Cov is reached [START_REF] Marelli | UQLab user manual --Structural reliability (Rare event estimation)[END_REF]. For a prescribed Cov, the smaller the probability of failure, the larger the sample size is needed in the MCS, which means dealing with rare events, MCS has low efficiency.

Subset Simulation (SS)

Subset Simulation is one of the advanced sampling methods and currently received great care from researchers to deal with rare events. This method presents an adaptive stochastic simulation, which consists of solving a series of simpler reliability problems. The SS appeared by [START_REF] Au | Estimation of small failure probabilities in high dimensions by subset simulation[END_REF] is also called 'variance reduction techniques' and generates samples that correspond to specified levels of failure probabilities progressively. The core concept starts from a small probability of failure can be stated as a product of larger conditional probabilities.

Let F stand for the target failure region in the space of the variables and let F1  F2  …  Fm = F be a selected sequence of m nested failure regions. Thus the conditional probability is identified as:

1 m 1 1 1 2 ( ) ( ) ( | ) P(F ) ... ( ) ( | )          m m m m i i i P F P F P F F P F P F F (1.26)
Equation ((1.26) shows that by computing the larger probabilities P(F1), P(Fi|Fi-1) (i=2, ..., m) and then producing them, we can determine the small probability (i.e., P(F)).

The procedure of this method is clear with few steps. Firstly, the initial probability P(F1) is concluded by direct MCS. In the next step, the conditional PDF is determined numerically using a Markov Chain Monte Carlo procedure. Besides, the SS makes use of a robust Metropolis-Hastings algorithm, which can efficiently generate random samples with any probability distribution.

In subset simulation, the size of the nested sets is usually chosen regarding the conditional probabilities equal to 0.1 (i.e., 90% of the samples). In other words, it is discarded and needs to be recomputed from the remaining 10% that serve as seeds. This work often leads in general to less than ten nested sets. Having too few nested sets would lead to a similar inefficiency as for direct MCS. However, it is argued that having many more nested sets would lead to an increasing number of samples, which again decreases the efficiency of subset simulation [START_REF] Tran | Reliability-based design of deep tunnel excavated in the viscoelastic Burgers rocks[END_REF].

Latin Hypercube Sampling (LHS)

When considering M independently and identically distributed (idd.) samples, LHS's concept is that the range of each variable input parameter is divided into M equiprobable strata, and one value is randomly selected from each range [START_REF] Helton | Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems[END_REF]. The LHS technique's leading profit is that it offers a boosted sampling number and quicker run-times with a smaller number of interactions. This method can be used to generate the initial design of experiment (DoE) of the approximation methods by metamodeling, which will be revisited later.

The LHS of a random vector x=(x1, x2. . . , xd), denoted (x (1) , x (2) . . . , x (M) ), gives a sample mean () reads:

1 1     M i i Y M (1.27)
According [START_REF] Stein | Large sample properties of simulations using Latin hypercube sampling[END_REF], the output Y = y(x) by using LHS has a smaller variance compared to the sample mean of the MCS. Figure 1-6 below is a particular example of two approaches to create a sample of size n=10, from two variables x = [x1, x2] where both x1 and x2 vary in the range of [0,1]; x1 follows uniform distribution while x2 has a normal distribution. This figure shows that LHS's result is more broadly scattered and does not show the clustering effects found in MCS [START_REF] Iooss | Numerical studies of the metamodel fitting and validation processes[END_REF].

Other Sampling methods

Besides Monte Carlo simulation, Subset simulation, Latin hypercube sampling, some other well-known simulation methods also have gained much attention. The first example is the Importance Sampling (IS), as referred to [START_REF] Bichon | Multimodal reliability assessment for complex engineering applications using efficient global optimization[END_REF][START_REF] Pradlwarter | Application of line sampling simulation method to reliability benchmark problems[END_REF]Ramsey , 2007;[START_REF] Sudret | Uncertainty propagation and sensitivity analysis in mechanical models-Contributions to structural reliability and stochastic spectral methods[END_REF][START_REF] Dubourg | Metamodel-based importance sampling for the simulation of rare events[END_REF][START_REF] Dubourg | Meta-model-based importance sampling for reliability sensitivity analysis[END_REF][START_REF] Papaioannou | Sequential importance sampling for structural reliability analysis[END_REF][START_REF] Al-Bittar | Probabilistic Analysis of Strip Footings Resting on Spatially Varying Soils Using Importance Sampling and Kriging Metamodeling[END_REF][START_REF] Lee | Failure probability calculation method using kriging metamodel-based importance sampling method[END_REF]. Then, the Directional Simulation (DS) that the readers can refer to [START_REF] Bjerager | Probability integration by directional simulation[END_REF][START_REF] Fenrich | A Comparison of Response Surface Methods for Reliability Analysis using Directional Simulation[END_REF][START_REF] Moarefzadeh | Implementation of directional simulation to estimate outcrossing rates in time-variant reliability analysis of structures[END_REF] for more details.

Global reliability methods

Brief introduction of surrogate models

Although the powerful approaches of local reliability methods and sampling methods and are impressive, they remain some disadvantages. The local methods here often fail in situations such as the limit state function (LSF) is highly nonlinear or non-smooth, or the LSF is multimodal (i.e., multiple design points) (Huang et al., 2017a). Therefore, in such situations, they are difficult to be accurately performed the LSF, even though the RSM can perform, but it may not provide encouraging results. Moreover, sampling methods usually require a huge number of performance function evaluations, making them unfeasible if the probabilistic performance is costly to estimate. The sampling methods with a low value of the failure probability cannot be used in all problems.

In engineering, the same program must repeatedly execute response functions due to a ton of random samples to optimize a set of input parameters. This leads to the most satisfying values of the output arguments. Thus, more advanced probabilistic approaches are required, requiring fewer calls to the mechanical model. That requires a so-called meta-model method for using the second level of modeling, then summarizing what has been learned from the simulations and inferring the original computer model's response. This method can cope with the disadvantages of simulation and especially provide structural reliability analysis in complex problems. The technique also could resolve the high non-linearity and time-consuming response of an original model (e.g., a traditional mechanical model) by a metamodel (i.e., a functional approximation or an analytical equation); therefore, significantly cutting down the size of computational calls such as in the classical MCS method. Kriging and kernel regression (i.e., radial basis functions, and support vector regression, etc. ) have proven significantly efficient in this context [START_REF] Vazquez | Sequential search strategies based on kriging[END_REF].

The metamodels are first developed as "surrogates" of the expensive simulation process to enhance the global computation efficiency. They are practical and convergent tools to assist a broad range of interests in modern engineering design [START_REF] Wang | Review of metamodeling techniques in support of engineering design optimization[END_REF][START_REF] Weinmeister | Analysis of a polynomial chaoskriging metamodel for uncertainty quantification in aerodynamics[END_REF]. The principle of meta-modeling (or surrogate modeling) consists of proposing a mathematical function, which replaces the often expensive numerical model representing the behavior of the system or structure. The metamodel's main idea is to fasten the calculation and enhance the accuracy of the analysis, especially in more complex and specific problems (e.g., high nonlinear, the tiny value of the failure probability, etc.). A metamodel is calibrated from a set of points, called experimental design, for which the numerical model has been evaluated. Table 1-2 below shows some recommendations for some common metamodels' right choice based on their characteristics. [START_REF] Hu | Mixed efficient global optimization for time-dependent reliability analysis[END_REF][START_REF] Zhu | Reliability analysis with Monte Carlo simulation and dependent Kriging predictions[END_REF][START_REF] Xiao | A new adaptive sequential sampling method to construct surrogate models for efficient reliability analysis[END_REF][START_REF] Jiang | A general failure-pursuing sampling framework for surrogate-based reliability analysis[END_REF]. It is used as an alternate model ("surrogate model") to replace the expensive performance function. Note that if any sampling method (e.g., as aforementioned) can find the failure probability, the surrogate one, which is even more efficient in most complex cases, also can be applied. Kriging metamodel (or Kriging predictor) is the geostatistical method of prediction by using the Best Linear Unbiased Predictor-BLUP [START_REF] Oliver | Springer Briefs in Agriculture: Basic steps in Geostatistics: The variogram and kriging[END_REF]. This model, as an essential surrogate model, provides better fitting accuracy than other surrogate models [START_REF] Kaymaz | Application of kriging method to structural reliability problems[END_REF][START_REF] Gaspar | Assessment of the efficiency of Kriging surrogate models for structural reliability analysis[END_REF][START_REF] Shi | Kriging response surface reliability analysis of a ship-stiffened plate with initial imperfections[END_REF]. Indeed, this method is more and more attractive due to its numerical accuracy by adding new training samples for each improving iterations [START_REF] Montgomery | Design and analysis of experiments[END_REF]. This model appears in uncertainty propagation [START_REF] Romero | Construction of response surfaces based on progressive-lattice-sampling experimental designs with application to uncertainty propagation[END_REF] and reliability [START_REF] Kaymaz | Application of kriging method to structural reliability problems[END_REF][START_REF] Bichon | Efficient global reliability analysis for nonlinear implicit performance functions[END_REF] studies, and also widespread in optimization. The technique is also known as the active-learning based Kriging method, and its performance function g(x) must satisfy all training points (as in Figure 1234567).

Kriging, developed by Krige [START_REF] Krige | Journal of the Chemical Metallu , Rgical & Mining Society of South Africa[END_REF], was firstly applied in the 1960s for statistical modeling of geospatial data [START_REF] Matheron | Principles of geostatistics[END_REF]. Later, it was applied to deal with structural reliability analysis and to grow up surrogate models [START_REF] Kaymaz | Application of kriging method to structural reliability problems[END_REF].

A training data set ={x (1) , …, x (M) } consists of total M realizations of the input random vector x, and corresponding responses of the model output g(x) are gathered in a vector as y = [g(x (1) ), …,g(x (M) )] T . The M value, the size of initial training points, can be small (e.g., M=12 or 20).

Following [START_REF] Santner | The design and analysis of computer experiments[END_REF][START_REF] Kaymaz | Application of kriging method to structural reliability problems[END_REF], the Kriging surrogate model reads:

( ) ( ) ( ) ( ) ( ) ( )    p T i i i g g = + z k + z x x k x β x β x x (1.28)
where p] T is a regression coefficient vector; k(x) represent explanatory polynomials, with its constants, are selected as the basis functions ki (x) (i = 1,⋯, p). In Eq. (1.28) the first term k(x) T  is the mean value of the Gaussian process (i.e., this is the deterministic part and represents the average trend of the process). While the second term, z(x), is the residual error (or a stochastic process) with zero-mean stationary Gaussian process and variance 2 z

 . The covariance function of two random points is defined as:

    2 z ),z( ) ,   i j i j C (x x x x z R (1.29)
In which,

    2 ' 1 ; , ' exp             θ x x d i i i i R x x (1.30)
where xi is the i-th component of vector x.

Due to this approximation metamodel, the performance function can be predicted for a realization x of the normally distributed random variables:

  2 ( ) ( ), ( )   x x x  g g g  (1.31)
Thus, the dependency of x and its i-th realization is described as:

  ( ) x x,x i i r ( )= R ( 1,..., ) i M  (1.32)
Then the vector represents cross-correlation between the prediction point x and each training point of DoE (M groups), reads:  

1 2 argmin det    R M z  (1.34)
It means that in the general multi-dimensional problem, the correlation function is defined from a set of scale parameters  of which each component represents the characteristic length scale of each input dimension. Contrary to this last case (known as the anisotropic case), one can use a single scan parameter  for all the input dimensions (called an isotropic case). Some widely used correlation functions in the literature are the exponential, squared exponential, and Second-order Markov functions (see (Li et al., 2015)).

The Kriging model can be identified by its mean and variance functions as follows

  1 2 2 1 1 1 ( ) ( ) ( ) ( ) ( ) 1 ( ) ( ) ( ) ( ) ( )            x k x β+r x R y-Kβ x r x R r x u x K R K u x T T g T T T z g (1.35)
where:

1 1 1 ( )     β K R K K R y T T 1 ( ) ( )    u(x) K R r x k x T 2 1 1 1 1 ( - ) ( - )      Y K β R Y K β T z M (1.36)
In Eqs. ((1.35) and (1.36)) the vector y gathers the response of the structure (i.e., the exact results of the performance function) evaluated at M 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:

( ) ( ) ( ) ( ) ( ), 1,..., ( , , ), , 1,..., ( ), 1,..., , 1,...,        x x x x i i i j ij i ij j y g i M R R i j M K k i M j p  (1.37)
By building a Kriging metamodel, the probability of failure can be found by using one of the well-known direct sampling methods, such as the MCS as in Eq (1.23). Table 1-3 below summarizes some general steps to conduct a Kriging-based reliability analysis. Following this iterative procedure, the DoE (or the matrix S) will be enriched by adding one or a set of new training points for each iteration according to a so-called learning function. The Active learning and Kriging-based methods concept is to make a relatively small DOE and sequentially add a new point (or more points) step by step to the DOE. Then the raised Kriging surrogate model progressively approaches the real performance function [START_REF] Tong | A Kriging-Based Active Learning Algorithm for Mechanical Reliability Analysis with Time-Consuming and Nonlinear Response[END_REF]. Many researchers pointed out that the bases of the methods are at selecting a learning function and choosing a new training point corresponding to this function. Following this, several common recent learning functions will be presented in this section.

In practice, AK-based methods can be utilized by MATLAB toolboxes such as MPerK [START_REF] Santner | The design and analysis of computer experiments[END_REF], DACE [START_REF] Lophaven | DACE -A Matlab Kriging Toolbox[END_REF], CODES [START_REF] Lacaze | CODES: A Toolbox for Computational Design, Version 1.0[END_REF], UQLab [START_REF] Lataniotis | UQLAB User Manual -Kriging (Gaussian process modelling[END_REF]. Other toolboxes of these methods can be found in different platforms as in R package: GeoR (Ribeiro and Diggle 2001), GSTAT (Pebesma 2004), RandomFields (Schlather 2012), and DiceKriging [START_REF] Roustant | DiceKriging, DiceOptim: Two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF]. Additionally, one can use a commercial product, so-called Statgraphics Centurion 19, for statistical analysis, data visualization and predictive analytics (http://www.statgraphics.com), that already includes Kriging feature.

U-criterion (AK-MCS)

This criterion was proposed by Echard [START_REF] Echard | AK-MCS: An active learning reliability method combining Kriging and Monte Carlo Simulation[END_REF], by combining the AK-based method and Monte Carlo Simulation (so-called AK-MCS). A candidate sample set  is appeared in the AK-MCS to represent the whole simulation domain for a result x* to update the Kriging model [START_REF] Zhang | REIF: A novel activelearning function toward adaptive Kriging surrogate models for structural reliability analysis[END_REF]:

  (1) (N) * argmin U( ,..., U( )  x x x (i) x    (1.38)
Herein, N is a sizable integer number (e.g., N=10 6 ), and the learning function U is defined as [START_REF] Echard | AK-MCS: An active learning reliability method combining Kriging and Monte Carlo Simulation[END_REF]:

| ( ) | | ( ) | ( ) ( ) ( )      x x x x x g g g g U (1.39)
Note that the ordinary regression model is often sufficient even for strongly nonlinear functions [START_REF] Bichon | Efficient surrogate models for reliability analysis of systems with multiple failure modes[END_REF]. The information then is exclusively carried by the autocorrelation function. The AK-MCS is the crucial method we use in almost all our case studies because of its utilities. The technique offers both the best estimate at an untired point and an approximation of the prediction variance (i.e., the mean squared error (MSE) in DACE toolbox). Thus, the prediction uncertainty can also be estimated by this variance and number as well as the location of training points (Huang et al., 2017a).

The coefficient of variation of the failure probability, Cov(Pf) defined as in Eq. (1.25)(1.25), being lower than 5%, is often used as the convergence criterion of this method.

EFF function (EGRA)

Expected Feasibility Function (EFF) learning function, using in the Efficient Global Reliability Analysis (EGRA) algorithm, gives a suggestion of how well the exact performance value at point x should fulfill the constraint function   .

x g z  In other words, EFF can find the points over the region z   , and defined as:

        ˆˆˆ( ) | ( )| ( ) ( ) x x x x z z EFF g z g f g d g          (1.40) Herein   ˆx g is the predicted response value at point x,     ˆx f g
which is the probability density function (PDF). Let suppose z + and z -denote for z   respectively, Eq.(1.40) can be rewritten as [START_REF] Echard | AK-MCS: An active learning reliability method combining Kriging and Monte Carlo Simulation[END_REF]:

ˆˆ( ) ( ) ( ) ( ) ( ) 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 ( ) ( ) ( ) ( ) ( ) x x x x x x x x x x x x x x x x x g g g g g g g g g g g g g g g g z z z EFF z z ( z z ) z z                                                                                                      ˆ( ) ( ) x x g g                    (1.41)
where z is a constant;  is proportional to the predicted standard deviation. In the case of reliability, set 0,

z  ˆ) 2 (x g   
. Thus, the expected values at the points being closed to zero (i.e., significant uncertainty) will get values of high expected feasibility. EGRA selects the position with the most massive expected value as the new adding sample, and the stopping condition is max (EFF(x)) ≤ 0.001.

The main difference between EGRA and AK-MCS approaches lies in that EGRA estimates the LSF in the whole space; at the same time, AK-MCS addresses itself to the accuracy of the sign of the LSF values among an MCS population (Huang et al., 2017a).

The learning function H

The function H (Lv et al., 2015) can be applied to obtain the uncertainty of g(x) based on the information entropy theory, which can measure the uncertainty. Let ( )

x g  and ( ) x g  denote for ˆ) 2 (x g  
; the function H can be defined as:

        ( ) ( ) ˆˆˆ( ) ( ) ln ( ) ( ) x x x x x x g g H g f g f g d g      (1.42)
The lower information entropy, the more certain the predicted response is. Thus, the candidate point that makes maximum H will be picked as a new adding sample. The above equation can be derived as an analytical form by: ˆˆˆˆˆ2

( ) ( ) 2 ( ) ( ) 1 ln 2 ( ) 2 ( ) ( ) ( ) 2 ( ) ( ) 2 ( ) ( ) 2 ( ) ( ) 2 ( ) ( ) 2 ( ) 2 ( ) x x x x x x x x x x x x x x x x x x g g g g g g g g g g g g g g g g g H                                                                                               (1.43)
The new training point corresponding to argmax{H(x)} is then selected while the convergence is attained when max{H(x)} 0.5 as proposed by (Lv et al. 2015).

Other learning functions

Recently, many learning functions have been appeared and frequently used in AK-based methods. First, the Expected Risk Function (ERF) [START_REF] Yang | Probability and convex set hybrid reliability analysis based on active learning Kriging model[END_REF] is to identify the risk that the sign of a point is not an exact prediction. Second, Least Improvement Function (LIF) is introduced by [START_REF] Sun | LIF: A new Kriging based learning function and its application to structural reliability analysis[END_REF], which takes not only the statistical formation provided by the Kriging model into account but also the joint PDF of random variables. Third, [START_REF] Gaspar | Adaptive surrogate model with active refinement combining Kriging and a trust region method[END_REF] proposed the trust region method for the active learning sampling process. Indeed, over the last two years, numerous researchers have published their proposed functions. For examples, (Wang and Shafieezadeh, 2019a) proposed an error rate-based adaptive Kriging for reliability analysis (REAK); the reliability-based expected improvement function (REIF) by [START_REF] Zhang | REIF: A novel activelearning function toward adaptive Kriging surrogate models for structural reliability analysis[END_REF]; an efficient error-based stopping criterion (ESC) for kriging-based reliability analysis methods by [START_REF] Wang | ESC: an efficient error-based stopping criterion for kriging-based reliability analysis methods[END_REF]. Moreover, a new pointselected learning function, called HPF by [START_REF] Tong | A Kriging-Based Active Learning Algorithm for Mechanical Reliability Analysis with Time-Consuming and Nonlinear Response[END_REF], is applied to the mechanical reliability analysis (involving time-consuming and nonlinear response). Especially for the time-variant hybrid reliability problem under random and interval uncertainties, the timevariant failure probability's upper limit is vastly interested. To adequately estimate it, the adaptive Kriging combined with design point-based importance sampling and meta-model based one are proposed by [START_REF] Ling | Adaptive Kriging coupled with importance sampling strategies for time-variant hybrid reliability analysis[END_REF], etc.

The different problems can be encountered with the diversity of Kriging-based methods or others as Table 1-4 below [START_REF] Lelièvre | La famille des méthodes AK pour l'analyse de fiabilité: synthèse des récentes avancées[END_REF] and Figure A1-2 in appendix A1. et al., 2016). In the context of uncertainty, quantification, PCEs, introduced by [START_REF] Ghanem | Stochastic Finite Element: a Spectral Approach[END_REF], have proven to be well-established tools. For example, its applications in uncertainty propagation [START_REF] Xiu | The Wiener--Askey polynomial chaos for stochastic differential equations[END_REF], sensitivity analysis (Le [START_REF] Gratiet | Metamodel-Based Sensitivity Analysis: Polynomial Chaos Expansions and Gaussian Processes[END_REF], and, to a lesser degree, structural reliability [START_REF] Sudret | Comparison of finite element reliability methods[END_REF].

On geotechnical engineering problems, sparse PCEs are often used to reduce the total computational costs and allow for efficient estimation of failure probabilities [START_REF] Schobi | Application of conditional random fields and sparse polynomial chaos expansions in structural reliability analysis[END_REF]. In which the polynomial functions can be applied as surrogate models to estimate the response of complex systems.

Although PCEs are often regarded as an efficient surrogate modeling technique due to their global convergence behavior, it has been employed only seldom in reliability analysis (see, e.g., [START_REF] Notin | RPCM: a strategy to perform reliability analysis using polynomial chaos and resampling: Application to fatigue design[END_REF]). This situation could result in their lack of accuracy in the model response distribution tails, which are essential in this field [START_REF] Marelli | An active-learning algorithm that combines sparse polynomial chaos expansions and bootstrap for structural reliability analysis[END_REF]. Moreover, most active-learning approaches with surrogates require several local error estimations to adaptively enrich a small set of model evaluations near the LSF. Yet, PCEs cannot offer a variance error such as AK-based methods.

Let consider a finite variance model Y=(X) representing the response of some quantity of interest Y to the random input X   M , modeled by a joint PDF fX. Thus, the functional inner product defined by:

  , ( ) ( ) ( ) ( ) ( )     X X X X x g h g x h x f x dx g h  (1.44)
Herein X represents the input domain. By independent input variables assumption, that is

1 ( ) ( )    X X M i i i f x f x , one can denote (X)
as the following generalized PCEs (see, e.g., [START_REF] Ghanem | Stochastic Finite Element: a Spectral Approach[END_REF][START_REF] Xiu | The Wiener--Askey polynomial chaos for stochastic differential equations[END_REF])

           X X  M Y y  (1.45)
Where y is a real coefficient vector and  is a multi-index that points out the degree of the multivariate polynomial   in each of the input variables Xi:

  ( ) 1        i M i i i X (1.46) Herein ( )   i
i is a polynomial of degree αi that belongs to the orthogonal polynomials family (polynomials and basis), as in Table 1-5 below regarding the marginal PDF fXi. The definitions of the DoE and the polynomial's degrees can be tricky [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF]. Moreover, accuracy evaluation requires cross-validation. Thus, in geotechnical engineering, polynomial chaos expansions are often combined with KARM, for example [START_REF] Schöbi | Combining Polynomial Chaos Expansions and Kriging[END_REF][START_REF] Gao | Fast Uncertainty Quantification in Engine Nacelle Inlet Design Using a Reduced Dimensional Polynomial Chaos Approach[END_REF][START_REF] Weinmeister | Analysis of a polynomial chaoskriging metamodel for uncertainty quantification in aerodynamics[END_REF] or connected with Bootstrap as in [START_REF] Marelli | An active-learning algorithm that combines sparse polynomial chaos expansions and bootstrap for structural reliability analysis[END_REF].

Other surrogate models

Some well-known metamodels can be listed as the artificial neural networks (ANN) [START_REF] Papadrakakis | Structural optimization using evolution strategies and neural networks[END_REF], Radial basis functions (RBF) [START_REF] Krishnamurthy | Comparison of response surface construction methods for derivative estimation using moving least squares, kriging and radial basis functions[END_REF], and extended radial basis functions (ERBF) [START_REF] Mullur | Extended radial basis functions: more flexible and effective metamodeling[END_REF]. In contrast to the traditional response surface methodology (RSM) only approximating a least-squares distance to the data points, the mentioned methods offer an interpolation surface that passes through all training data points. These methods do not assume a unique shape for the approximation (c.f., polynomials in RSM); alternately, they characteristically consist of a group of functions, each incorporated with individual points in the design space [START_REF] Miro | Calibration of numerical models considering uncertainties[END_REF].

Recently, new complex techniques for global reliability methods are investigated in many different approaches. Firstly, the combining Proper orthogonal decomposition (POD) method, also recognized as Karhunen-Loève decomposition [START_REF] Karhunen | Zur spektraltheorie stochasticher[END_REF], with radial basis functions (shortly, POD-RBF) developed by [START_REF] Buljak | Proper orthogonal decomposition and radial basis functions algorithm for diagnostic procedure based on inverse analysis[END_REF]. Secondly, this method's extension is called POD-ERBF by using an extended version of RBFs (instead of the original RBFs) [START_REF] Mullur | Extended radial basis functions: more flexible and effective metamodeling[END_REF]. More recently, a hybrid RNN-GPOD surrogate model for realtime settlement predictions in mechanized tunneling is also taken into account by [START_REF] Cao | A hybrid RNN-GPOD surrogate model for real-time settlement predictions in mechanised tunnelling[END_REF], combining Gappy proper orthogonal decomposition (GPOD) and recurrent neural networks (RNN).

Stop criteria in Meta-modelling

It has been shown that the stopping criterion based on exceedance of the max (or min) value of the chosen learning function concerning an allowable amount can be too conservative. For example, in [START_REF] Gaspar | A study on a stopping criterion for active refinement algorithms in Kriging surrogate models. Safety and reliability of complex engineered systems[END_REF], by using the AK-MCS-IS method, the authors showed that the failure probability prediction could be stabilized much sooner than the stopping criterion defined by min{U(x)}>2 is satisfied. Following the two authors, the additional samples in DoE can have any more significant contribution after the stabilization of the exceedance probability even if their U values are smaller than 2.0. More comprehensive discussions on the limitation of the EGRA and AK-MCS can be found in the work of (Hu and Mahadevan, 2016a). Follow these authors, the variance of the probability prediction consists of two parts: the first part comes from the responses of MCS samples. In contrast, the second part comes from the mutual effects of these individual responses. Thus, the EFF and U learning functions focus only on reducing the different variances in the first part, and their convergence criteria are not defined from the aspect of reliability analysis.

From the previously mentioned discussions, it can be stated that the AK-based models can have sufficient accuracy in terms of failure probability prediction, and hence a stopping criterion based on a defined stabilization condition of the probability is preferred. In our study, we apply a proposed rule by [START_REF] Gaspar | A study on a stopping criterion for active refinement algorithms in Kriging surrogate models. Safety and reliability of complex engineered systems[END_REF] as follows:

  ( ) (1) (1) , 2,...,       i f f f P P i N P (1.47)
where (1) f P is the reference probability of failure using to detect the stabilization, and ( ) i f P (i=2,…, N) is the probability in the i-th iteration. The allowable tolerance =0.015 and N=100 are suggested in [START_REF] Gaspar | A study on a stopping criterion for active refinement algorithms in Kriging surrogate models. Safety and reliability of complex engineered systems[END_REF] for closed-form (analytical) solutions. Besides, these authors proposed that this stopping criterion should be interpreted as an additional criterion to the original proposal by (Echard et al., 2011) (i.e., min{U(x)}>2 in case of the AK-MCS method).

A quite similar stopping criterion, called ks-fold cross-validation, was proposed in [START_REF] Xiao | A new adaptive sequential sampling method to construct surrogate models for efficient reliability analysis[END_REF]. Concerning the rule in Eq. (1.47) these last authors used the failure probability of the current iteration (i.e., the probability estimated from all training samples) as the reference exceedance probability (1) f P . Further, the probability ( ) i f P (i=2,…, N) are the probability of exceeding calculated from the previous (N-1) iterations using the surrogate model that is built with the ith subset left out of the current iteration. The value of =0.01-0.02 was proposed in [START_REF] Xiao | A new adaptive sequential sampling method to construct surrogate models for efficient reliability analysis[END_REF], while a smaller amount could be selected for more accurate results.

Reliability analysis of underground structures: state-of-the-art

Reliability evaluation can evaluate the combined effects of uncertainties in the input parameters and provide a supplement to conventional structural analysis [START_REF] Gharouni-Nik | Reliability analysis of idealized tunnel support system using probability-based methods with case studies[END_REF]. Generally, the ability to tolerate the existing loads of the liner plays a vital role in tunnel excavation. Thus, tunnels must be considered with the applicable structural system to afford these loads. The principles of probability and mathematics for engineering reliability which is now a valuable tool for approximating the appropriate response of systems. This section concentrates on the basics of the reliability in geotechnical engineering to deal with uncertainties, especially with underground structures like tunnels constructed in the rock-mass.

Application of general reliability analysis problems

Recently, reliability analysis has been received much attention from researchers, especially the progress towards practical applications of such methods. However, almost all the specifications and studies concern and are for geotechnical structures in soil problems. Rock engineering principles seem to be neglect in Eurocode 7 (EC7), as highlighted by [START_REF] Harrison | Eurocode 7 and rock engineering: current problems and future opportunities[END_REF]. This reason could be from the fact that the discontinuous, heterogeneous, and anisotropic nature of the rock mass requires a robust deal of empiricismto be involved [START_REF] Lamas | Some key issues regarding application of Eurocode 7 to rock engineering design[END_REF].

The closed-form or semi-closed-form solutions were widely applied to verify various probabilistic methods or surrogate models. Herein, Monte-Carlo simulation (MCS) can be listed as the most application without computational effort, and the MCS results are usually considered as the typical benchmark for comparing other methods, e.g., [START_REF] Lu | Reliability analysis of ground-support interaction in circular tunnels using the response surface method[END_REF][START_REF] Lü | Probabilistic analysis of underground rock excavations using response surface method and SORM[END_REF][START_REF] Idris | Comparison of different probabilistic methods for analyzing stability of underground rock excavations[END_REF]. The FORM and SORM also increasingly received practical attention after a new proposed spreadsheet-based algorithm by (Low and Tang, 1997a;[START_REF] Low | Reliability analysis using object-oriented constrained optimization[END_REF][START_REF] Low | Efficient spreadsheet algorithm for first-order reliability method[END_REF] for both independent standard Gaussian and dependent non-normal variables, such as many types of research by [START_REF] Li | Reliability analysis of circular tunnel under hydrostatic stress field[END_REF][START_REF] Goh | Reliability assessment of stability of underground rock caverns[END_REF]Wang et al., 2013).

During the last decades, many practical procedures for the application of probabilistic stability analyses of underground structures have been provided. One can find the application of FORM/SORM in underground structural analysis by [START_REF] Cai | Refined second-order reliability analysis[END_REF][START_REF] Laso | A level II reliability approach to tunnel support design[END_REF][START_REF] Celestino | Evaluation of tunnel support structure reliability[END_REF][START_REF] Shin | Methodology for quantitative hazard assessment for tunnel collapses based on case histories in Korea[END_REF]. In which almost all problems used classical reliability techniques and relied on a ground-support interaction diagram. Mollon presented a probabilistic analysis and design of circular cross-section tunnels against the stability of the working face by using RSM (Mollon et al., 2009b). Then, [START_REF] Li | Reliability analysis of circular tunnel under hydrostatic stress field[END_REF] followed by [START_REF] Lu | Reliability analysis of ground-support interaction in circular tunnels using the response surface method[END_REF]Lü and Low, 2011a), applied into circular tunnel structures using RSM and SORM. Some researchers concentrated on the probabilistic selection of appropriate rock quality in tunnel surrounding area. [START_REF] Oreste | A probabilistic design approach for tunnel supports[END_REF] performed a probabilistic calculation and used some design tools for tunnel supports. His probabilistic approach has sufficient data on the quality of the rock mass around the tunnel. That leads to a better understanding of the risks, more efficient geomechanical zoning; finally, a more reliable estimation of the costs.

Probabilistic investigations on structural design phenomena have been paid much interest in recent decades. For example, Oste [START_REF] Oreste | A probabilistic design approach for tunnel supports[END_REF] evaluated the geomechanical rock quality's uncertainty to the design of preliminary tunnel supports. [START_REF] Celestino | Evaluation of tunnel support structure reliability[END_REF] evaluated the failure probability by using MCS according to load and resistance factor design principles for the failure modes. Follow these authors, the distribution of tunnel strength capacity was obtained and compared to the load by varying tunnel support parameters. [START_REF] Langford | Reliability based approach to tunnel lining design using a modified point estimate method[END_REF] given a modified PEM for the reliability-based analysis, and then applied it to account for Hoek-Brown parameters' variability when designing the tunnel lining. [START_REF] Lü | System reliability assessment for a rock tunnel with multiple failure modes[END_REF] proposed a new approach for assessing the system reliability of rock tunnels with three failure modes (i.e., inadequate support capacity, excessive tunnel convergence, and insufficient rock bolt length). Note that the reliability of specific components of the tunnel has been studied separately toward nondeterministic conditions. For example, [START_REF] Lü | System reliability assessment for a rock tunnel with multiple failure modes[END_REF] considered rock bolts in the analysis, and [START_REF] Yang | Reliability analysis of shotcrete lining during tunnel construction[END_REF] mainly focused on shotcrete in their reliability analysis. Nevertheless, the question of the reliability of the entire tunnel support system concerning all assembled components is still unresolved in the literature. Briefly, in the past, numerous studies have dealt with the issues of uncertainty of ground conditions and the assessment of rock quality in the tunnel face linking to mechanical properties or geotechnical investigations.

For the-state-of-the-art of Active learning and Kriging-based methods (AK-based methods), one can refer to the summary by [START_REF] Lelièvre | AK-MCSi: A Kriging-based method to deal with small failure probabilities and timeconsuming models[END_REF] as in Appendix A.1. While (Li et al., 2016) utilized the uniform design and support vector machine (SVM) to develop a hybrid methodology for probabilistic tunnel stability analysis.

Spatial variability problems

In geotechnical applications, spatial variability can be considered as the most significant source of uncertainty [START_REF] Shokri | A review study about spatial correlation measurement in rock mass[END_REF]. Thus the characterization of geotechnical variability is the most critical element [START_REF] Phoon | Reliability of geotechnical structures in ISO2394[END_REF]. Because of the inherent existence of inhomogeneity and discontinuity, the local material strength and the mechanical response in rocks are profoundly affected by spatial variability [START_REF] Ji | Effect of 2D spatial variability on slope reliability: A simplified FORM analysis[END_REF].

In geological materials, spatial variability rises under its formation by natural processes acting over an incredibly long time, such as millions of years. It endows Geo-material with some unique statistical features not typically observed in structural material assembled under rigorous quality control. This variability is often characterized by the spatial correlation length, which describes the distance over which the input parameters are identical or correlated. The closer distance (i.e., adjacent locations), the more likely similar parameters are. To characterize the spatial variability through the field, one frequently applies a correlation function such as Markovian or Gaussian. This function shows the connection between the spatial correlation length and the distance separating two points.

The spatial variability can be completely defined by a PDF (i.e., mean and standard deviation), and an autocorrelation function (ACF) with a corresponding autocorrelation distance (i.e., the scale of fluctuation). Generally, at a specific location, rock properties are closer to adjacent locations than those at distant sites. In geotechnical practice, it is difficult to get the real ACF for rock-based on-site investigation data limitation. Yet, since the different CAFs produce insignificant differences in results. Among theoretical ACFs, the squared exponential function is the most well-known ones in the literature (Li et al., 2015;[START_REF] Jiang | Efficient slope reliability analysis at low-probability levels in spatially variable soils[END_REF]; it is also utilized in our study as follow [START_REF] Haj | Probabilistic analysis of an offshore monopile foundation taking into account the soil spatial variability[END_REF].

  2 2 1 2 1 2 1 1 2 2 | | | | ( , ),( , ) exp                                       h v x x
    i i i h v (El
Table A1-4 in the Appendix shows the review of the literature relating to spatial variability problems in the recent decade in the broad context of geotechnical engineering. Indeed, when material properties vary spatially, analyses are often carried out with the help of finite element models (FEM). Thus, in the table, we also mention the list of FEM software that researchers have applied. If the effect of spatial variability was discussed in many problems such as the slope stability, foundation, and dam construction, the contribution in the tunnel structures seems very limited with just only several recent works, see for example (Huang et al., 2017a;[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] Hu | Ground surface settlement analysis of shield tunneling under spatial variability of multiple geotechnical parameters[END_REF]. From these studies, the authors showed that the failure probability could be markedly overestimated if the spatial variability is neglected. In addition, the vertical scale of fluctuation has a much significant effect than the horizontal level of impact on the probabilistic results of tunnel convergence.

The long-term behavior problems

So far, tunnel analysis and tunnel support design have mainly been carried out based on hypotheses of the time-independent behavior of the massif. Typically, the excavation process is assumed to be immediate, yet the linings are supposed to be a quick process installing sometime after or during the full excavation. These statements could cause an underestimation of the attained results since an extensive range of rocks presents a considerable time-dependent reaction, even in many cases that can participate to more than 70% of the total convergence [START_REF] Sulem | Closure analysis in deep tunnels[END_REF]. Several kinds of research in the field [START_REF] Arnau | Longitudinal time-dependent response of segmental tunnel linings[END_REF][START_REF] Barla | Time-dependent modeling of tunnels in squeezing conditions[END_REF][START_REF] Sharifzadeh | Time-dependent behavior of tunnel lining in weak rock mass based on displacement back analysis method[END_REF][START_REF] Vu | Semi-analytical solution for stresses and displacements in a tunnel excavated in transversely isotropic formation with nonlinear behavior[END_REF][START_REF] Bui | Analytical modeling of a deep tunnel inside a poro-viscoplastic rock mass accounting for different stages of its life cycle[END_REF][START_REF] Manh | Anisotropic timedependent modeling of tunnel excavation in squeezing ground[END_REF][START_REF] Zhang | Study on time-dependent behavior and stability assessment of deep-buried tunnels based on internal state variable theory[END_REF] shown that the delay exemplifies a crucial influence on either the final tunnel convergence or the stability of the liners.

From the deterministic point of view, the consideration of the time-dependent effect on the analysis and design of the tunnel can be conducted by adopting an appropriate time-dependent mechanical behavior law of rock mass. Many constitutive models, such as linear viscoelastic, non-linear viscoelastic, and Elasto-Viscoplastic law, were presented in the literature. They aim to simulate accurately the results observed in the laboratory (usually through the uniaxial and/or triaxial creep tests) or the measurement of tunnel convergence conducted in-situ. More and more sophisticated models that allow accounting for not only the Elasto-Visco-plastic phenomenon but also the anisotropic effect were presented in the literature. If these models fit quite well with different observations, the determination of many parameters involved in the models presents a considerable challenge regarding the limited available experimental results. This challenge will be emphasized in the context of reliability analysis when the quantification of uncertainty of the parameters is challenging, even impossible. Consequently, a simple model of rocks involving a lower number of parameters seems preferred for the time-dependent reliability analysis.

The simplest models that could consider the time-dependent effect on the mechanical behavior of rock mass are the linear viscoelastic laws developed from some well-known rheological models. There are several popular viscoelastic models built by combining series or parallels of Maxwell's model and Kelvin's model. A brief synthesis of some commonly used rheological models to characterize the linear viscoelastic behavior of rock mass will be presented in the next section. The advantage of using these models is the possibility to deduce the analytical solutions of deep tunnels excavated in rheological rocks, which allows them to quickly investigate the parametric study and/or validate the developed numerical simulation. This work is also crucial for the reliability analysis problem when one would like to confirm a new method by comparing it with the referent solution of the direct MCS using the analytical solution. It is the case of the present research. In the next chapter, an analytical solution of the deep tunnel excavated in the viscoelastic Burgers rock is firstly presented, while the direct MCS is carried out in the next step to investigate the time-dependent failure probability of underground structure. These results are then utilized to validate a new metamodel that helps to improve the efficiency of the well-known AK-MCS. The validation allows us to apply this new metamodel to study the tunnel in the more appropriated time-dependent constitutive law, such as the nonlinear viscoelastic model of the rock mass. Note that, in this context, the analytical solution is no longer available, and the utilization of the numerical simulation, which is usually expensively time-consuming, is necessary.

From the knowledge of the present author, the consideration of time effect on the reliability analysis of tunnel in rock mass has not yet been conducted in the literature. The difficulty of this kind of problem is multiple. Firstly, as mentioned above, the quantification of parameters' uncertainty is difficult due to the limited experimental results. Secondly, the duration of the creep-test in the laboratory (some hours to several months) or the measurement in situ (several months or even several years) seems to be very short as compared to the life service of the tunnel (designed sometime for hundred years). Thus, the quantification of uncertainty of the time-dependent mechanical properties of rock mass from these measurements will be used to extrapolate the result at the very long-term and contribute to (an even more significant) uncertainty. However, the quantification of this kind of uncertainty (i.e., for the extrapolation) is an open question and will not be discussed in the present research.

Brief syntheses of some constitutive models of rocks

A summary of some commonly used linear-and non-linear-viscoelastic rocks constructed from rheological models is presented in this section. In general, viscoelastic deformation is characterized by a non-linear stress-strain response that travels along a different path upon unloading. Nevertheless, given enough time after the loadings removing, the stress-strain curve approaches the origin [START_REF] Aboudi | Micromechanics of Composite Materials: A Generalized Multiscale Analysis Approach[END_REF]. Note that the well-known linear-viscoelastic, socalled Burgers model, which will be used in the next chapter, belongs to this class of materials.

Then, another famous non-linear viscoelastic one, which will be used as the fourth-and lastchapter model, is also tacked.

The time-dependent constitutive law here can be understood as a stress-strain-time relationship, defined by loading a specimen with constant stress (creep) or strain (relaxation or isometric), illustrated as the curves below (ASTM Standard D2990-09, 2011): In practice, viscoelastic constitutive models are regularly utilized to ascertain the strains and stresses in the rock mass and to find suitable supports. Viscoelasticity is a time-dependent behavior, characterizing materials illustrating both an elastic and a viscous behavior when deformation is imposed [START_REF] Carret | Linear viscoelastic characterization of bituminous mixtures from dynamic tests back analysis[END_REF]. The elastic deformation explains the short-term response to stress, while the viscous one defines the long term response. In contrast to the elastic case, the viscoelastic material response depends on the instantaneous value of the stress as well as its entire history. Besides, viscoelastic behavior presents sensitivity to the speed of the stress. The smaller the amplitude of the response, the higher the speed of the stress is. Such an elastic behavior, viscoelastic behavior, is said to be linear when the response is proportional and additive to the stress. This note can be highlighted from creep or relaxation tests for different values of the stress [START_REF] Salençon | Viscoélasticité pour le calcul des structures[END_REF]. In our study, rheological characteristics are identical in all directions, so-called isotropic materials, and the constitutive equations focus only on primary and/or secondary creep. With the hypothesis of small transformations, the constitutive equations of a general linear viscoelastic solid could be shown as follows [START_REF] Wang | Analytical solutions for tunnels of elliptical cross-section in rheological rock accounting for sequential excavation[END_REF]:

0 0 ( , ) ( ) ( ,0) ( ) , ( , ) ( ) ( , 0) ( ) ( , ) 2 ( , ) 3                                x x x x x x v t ij v v ij ij v t v v kk kk kk de G t e G t dt d d K t K t dt d s t t (1.49)
where x is the position vector; v ij s and v ij e are the stress and strain deviator tensors, respectively (the superscript 'v' stands for viscoelastic) and are defined as (Song et al., 2018a):

1 , 3 v v v ij ij ij kk s      1 3 v v v ij ij ij kk e      (1.50) Herein v ij  and v ij
 are the stress and strain tensors, ij  is the unit tensor; and G(t) and K(t)

represent the shear and bulk relaxation moduli, respectively. Through viscoelastic processes, the relaxation can be integrated for both the short term and long term [START_REF] Nedjar | An approach to the modeling of viscoelastic damage. Application to the long-term creep of gypsum rock materials[END_REF].

Based on analogies to Spring-and-Dashpot models, there are four commonly used physical viscoelastic models, namely: The primary forms of stress-strain-time relationship for these models under uniaxial loading are synthesized based on some formulas in recent articles, e.g., [START_REF] Mogilevskaya | International Journal of Rock Mechanics and Mining Sciences A lined hole in a viscoelastic rock under biaxial far-fi eld stress[END_REF] and [START_REF] Paraskevopoulou | Analysis of time-dependent deformation in tunnels using the Convergence-Confinement Method[END_REF] as follow: No. Medium Forms of stress-strain relationship (a) Maxwell

M M M E          ( ) 3 M M t t E       (1.51) (b) Kelvin K K E        ( ) 1 exp 3 K K K E t t E                   (1.52) (c) Boltzmann (Kelvin-Voigt) M K K K K M M E E E E E             ( ) 1 exp 3 K M K E t t E                    (1.53) (d) Burgers M M K M K M K M M K K M K K E E E E E E                              ( ) 1 exp 3 3 K M K M E t t t E                      
(1.54)

The overdot indicates the time derivative;  stands for the viscosity of the dashpot element (; E denotes the Young modulus of the spring and (E. Notably, in the linear viscoelastic region, the well-known principle called Boltzmann superposition principle (or integral) explaining the response of a material to separate loading histories is vital in forecasting the creep behavior of plastic materials in a variety of test conditions. This principle is that the reaction of a material to a provided load does not depend on the material's response to any pressure already existing in the material. In other words, the creep strain at any time is the sum of the deformations due to each change in stress, σ, that has occurred, where each σ is considered to act continuously as applied stress, starting from the time it happens to any later time [START_REF] Dowling | Mechanical behavior of materials: engineering methods for deformation, fracture, and fatigue[END_REF].

Constitutive models of viscoplastic rocks

The non-linear viscoplasticity model, like linear viscoelasticity, still considers the strains to be infinitesimal. Nevertheless, with the non-linear case, doubling the loading no longer results in a doubling of the response. As a result, the constitutive model is non-linear, and the Boltzmann superposition integral is no longer applicable. In many cases, the stress level at which the material becomes non-linear appears to decrease along with increasing time [START_REF] Aboudi | Micromechanics of Composite Materials: A Generalized Multiscale Analysis Approach[END_REF]. The long-time predictions may require the use of non-linear viscoelastic constitutive equations, even for very low-stress levels. In this section, we revisit a short overview of some common law used for rocks (i.e., Power's law, Lemaitre's model, so on)

The power law is based on the original expression suggested by [START_REF] Obert | Creep in mine pillars: Report of investigation[END_REF]. The relationship between strain, stress, and time can be shown as

1        t ε σ m n a k t (1.55)
where ε a is the axial strain at time t, σ is equivalent stress ( 1 3    ), n is an index of the power function between stress and strain, m denotes the index of power function between strain and time (creep coefficient), and k is constant at the reference time t1 and is related to the modulus of the material.

Other forms of Power's law can be found in the literature; for example, Bailey and Norton [START_REF] Betten | Creep mechanics[END_REF] form:

' ' '   ε σ t ε σ t  n m a n m a A A (1.56)
Herein the parameters A, A', m, and m' remain as the material constants that are generally independent of stress and can be determined in a uniaxial creep test; n and n' are stress exponents.

Among various possibilities, in our thesis (especially in chapter 3 and chapter 4), we choose Lemaitre's model to describe the time-dependent behavior of host rock. In fact, most of the authors dealing with the time-dependent behavior of COx claystone proposed various models based on the elasto-visco-plasticity concept (Seyedi et al., 2017;[START_REF] Stavropoulou | Time-dependent behaviour of the Callovo-Oxfordian claystone-concrete interface[END_REF]. Even if many sophisticated models have been taken into account the unique features of COx rock behavior (see, for example, (Armand et al., 2017c;[START_REF] Mánica | A time-dependent anisotropic model for argillaceous rocks. Application to an underground excavation in Callovo-Oxfordian claystone[END_REF]Seyedi et al., 2017;[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] for more details), the non-linear viscoelastic model of Lemaitre is chosen for its simplicity and limited numbers of constitutive parameters.

By neglecting the thermo-mechanical (or thermo-hydro-mechanical) coupling effect, the total strain tensor tot ε reads:

ε   tot e v   (1.57)
where ε e and ε v indicate the elastic and viscous strain tensors, respectively. The viscous strain rate ( ε  v ) is expressed thanks to Lemaitre's law:

1/ 3 1 2           σ ε   n eq v m eq K (1.58) where 1/ 2 1 : 3            ε ε    n eq v v m K (1.59) 1 tr( ) 3   σ σ σ I  (1.60)
In these equations σ and σ  are the stress tensor and its corresponding deviatoric part, respectively, and I indicate 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 cumulate viscous strain (i.e., the time-dependent distortion strain) and the equivalent stress (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.

Put = n/m, the cumulate viscous strain is explicitly expressed as:

1 1 1 1 1 n 1 (1 )                  n eq t K (1.61)
Equivalently, the more compact form is written as:

   B C eq A t
(1.62)

with:

1 1 n 1 1 (1 ) ; ; 1 1                  n A B C K (1.63)
Besides, the positive triple is then inferred based on the parameters A, B, C as follows:

1 ;

; 1

      C B B B B n m K A C C C (1.64)
Note that creep in the viscoelastic model is generally assumed proportional to the deviatoric stresses, while hydrostatic pressures alone will not produce creep effects. Thus, creep effectively reduces the shear modulus and Young's modulus of a material, while the bulk modulus is not affected by it [START_REF] Fjar | Petroleum related rock mechanics[END_REF].

Summary

This chapter revisited some central concepts of reliability analysis, from the sources of uncertainty to their inclusions, and its applications for underground structures in design problems. Afterwards, the head of uncertainties can be categorized as aleatory and epistemic. Aleatory uncertainty is listed as temporal and spatial. In comparison, parameters and models belong to epistemic uncertainty cases.

Quantitatively, the underground structure's stability analysis is conducted by measuring the failure probability thanks to using an appropriate probabilistic approach. The most famous method is the Monte Carlo Simulation (MCS), which could provide an accurate estimate of the probability due to many trials. This approach has been primarily considered the benchmark to validate the other probabilistic methods (e.g., FORM/SORM or RSM), demonstrating their feasibility and effectiveness in different projects. However, since the sampling methods like the MCS require a massive number of structure response evaluations, they seem unfeasible in the case of rock formations due to complex behavior that can only be solved by numerical tools.

Recently, many advanced probabilistic approaches have been proposed in the literature to overcome these disadvantages. Generally, these approaches aim at establishing a metamodel (or surrogate model), which estimates the implicit LSF by a mathematical function. The reliability analysis is then conducted on the built metamodels (e.g., Kriging, RBF, ANN, SVM, Polynomial Chaos, etc.).

From this bibliographic survey, the Kriging metamodeling technique will be chosen for the reliability analysis in the next chapters thanks to its flexibility in interpolating the sample points and high accuracy, notably for nonlinear problems. Although this metamodeling's approved performance in many structural design projects, it is surprising that it is rarely applied in geotechnics, particularly in the rock engineering field. Here, it is essential to emphasize that almost all studies in the literature considered the uncertainty and/or variability of host rock's short-term mechanical properties on the deep underground structure's stability.

CHAPTER 2. RELIABILITY ANALYSIS OF DEEP TUNNELS EXCAVATED IN THE LINEAR VISCOELASTIC ROCKS

Introduction

In this chapter, the time-dependent reliability analysis of tunnels excavated in the linear viscoelastic rock will be investigated based on analytical approaches. Our analytical solutions are utilized here as a first assessment of the design parameters, offering guidance in the conceptual stage of the design process and a validation step (before tackling numerically).

Precisely, in the early stage, we estimate the time-dependent exceedance probability of tunnel based on the direct MCS, which is conducted by using a new closed-form solution, presented in [START_REF] Do | Time-dependent reliability analysis of deep tunnel in the viscoelastic Burger rock with the sequential installation of liners[END_REF]. Since a common constitutive model employed in rock engineering, the Burgers rock is chosen among linear viscoelastic models mentioned in section 1.4.1. While, in the second stage, an extension of the well-known AK-MCS metamodeling with a novel learning function is proposed, a so-called modified AK-MCS. Then, this modified function's validation work is investigated by comparing it with the results obtained in the first step. This work is especially important because it allows us to apply this technique in the next sections in the specific study-cases with much more complex problems (i.e., the viscoplastic behavior of rock mass, tri-linear elastic of the compressible material such in chapter 3; and additional complicate task with consideration of the spatial variability in chapter 4).

Closed-form solution of double-lined tunnels in linear viscoelastic rock

Theoretically, in linear-viscoelastic rocks, one can consider the time-dependent behavior by numerous analytical approaches [START_REF] Fahimifar | Analytical solution for the excavation of circular tunnels in a viscoelastic Burger's material under hydrostatic stress field[END_REF][START_REF] Nomikos | Supported axisymmetric tunnels within linear viscoelastic burgers rocks[END_REF]Wang et al., 2013[START_REF] Wang | An analytical approach for the sequential excavation of axisymmetric lined tunnels in viscoelastic rock[END_REF]Song et al., 2018b). More specifically, these rocks are often modeled with the combination of Hookean elastic springs and Newtonian viscous dashpots, such as in Table 1-6 and Figure 1-9 (Chapter 1). Especially, in the literature, many closed-form solutions of deep tunnels in the viscoelastic Burgers model, which allows considering the transient and secondary steady-state, are proposed. For instance, regarding [START_REF] Fahimifar | Analytical solution for the excavation of circular tunnels in a viscoelastic Burger's material under hydrostatic stress field[END_REF], these authors investigated the effect of the tunnel face advancement and the stress history before tunnel construction. Nevertheless, their analytical solution was developed under the initial assumption that the behavior of the rock had no effect on the variation of the rock support pressure. Progressively, [START_REF] Nomikos | Supported axisymmetric tunnels within linear viscoelastic burgers rocks[END_REF] have derived a closed-form expression from accounting for the stress state before excavation of the tunnel and the time elapsed between excavation as well as the installation of available support. It should be noted that, in many practical projects, the tunnel is usually driven in complex geological conditions and the tunnel cross-section has a large size. Thus, a sequential excavation and installation of more than one-liner play a vital role in preventing collapse [START_REF] Miura | Study on design and construction method for the New Tomei--Meishin expressway tunnels[END_REF][START_REF] Mason | Support provided to rock excavations by a system of two liners[END_REF][START_REF] Sharifzadeh | Design of sequential excavation tunneling in weak rocks through findings obtained from displacements based back analysis[END_REF]. In other words, the time-dependent behavior of tunnels has crucial effects due to the sequential construction process. Numerous contributions have been proposed to address the problems of installing sequential supports, for example, (Wang et al., 2013[START_REF] Wang | An analytical approach for the sequential excavation of axisymmetric lined tunnels in viscoelastic rock[END_REF]; however, their proposals are still limited for tunnels excavated mainly in generalized Kelvin viscoelastic rock masses.

Viscoelastic of rock mass with double-liners tunnel

A deep double-lined tunnel excavated in a homogenous, isotropic, and incompressible viscoelastic rock mass subjected to hydrostatic far-field stress 0 h p will be investigated in this section. The assumption of incompressible material allows us to obtain analytical solutions rapidly. Notably, as shown in Figure 2-1, the two liners can be installed at separate instants t1 and t2. Furthermore, the first lining is assumed to be firmly bound to the rock mass as well as to the second.

The tunnel excavation can be sequential and represented by a variety of radii over time as follows: The effect of the tunnel face advancement on the interested tunnel section is taken into account.

0 0 ( ) 0 ( ) if ini fin R g t t t R t R t t         (2.
Precisely, a fictitious internal pressure applied to the tunnel's perimeter is mainly considered. The decreasing pressure starts from the initial stress state 0 h p to zero as follows (see also [START_REF] Wang | An analytical approach for the sequential excavation of axisymmetric lined tunnels in viscoelastic rock[END_REF]:

  0 0 ( ) 1 ( ) , h p t p t    (2.3)
Herein the dimensionless parameter (0 1)     can take the following expression:

2 1 ( ) 1 exp . , ( ) l m t m v t R t            (2.4)
In this equation, one assumed that the constant rate l v (i.e., advancement rate) is sufficiently minor to ignore the dynamic effects. The two coefficients t t  respectively. The correspondingly supported pressures 1 ( ) p t 2 ( ) p t at the outer face of these liners read:

1 1 11 1 2 12 2 0 0 ( ) ( ) if ( ) t t p t p t t t t p t t t            (2.5) 2 2 22 2 0 0 ( ) if ( ) t t p t p t t t        (2.6)

Constitutive mechanical behavior law of materials

The time-dependent behavior of rock masses can be explained by the viscoelastic Burgers model, which exhibits the rheological behavior of incompressible rock mass under deviatoric loadings (see Figure 1-9.d). In which, the four-component Burgers model consisted of the connection in series of the Maxwell-and the Kelvin-element. Under instantly applied stress, the initial reaction is determined by the spring of the Maxwell element (marked by parameter GM) while the Kelvin element (GK and K) and the dashpot of Maxwell element (M) control the transient creep and secondary steady-state creep, respectively [START_REF] Jaeger | Fundamentals of rock mechanics, 4th edn Blackwell[END_REF][START_REF] Nomikos | Supported axisymmetric tunnels within linear viscoelastic burgers rocks[END_REF]. Under deviatoric loadings, the behavior of incompressible rock mass can be described by an expression:

( , ) 2 ( )* ( , )

ij ij s r t G t de r t  (2.7)
In which the deviatoric tensors of stress and strain are computed as:

1 1 , 3 3 ij ij ij mm ij ij ij mm s e           (2.8)
The asterisk (*) in Eq. (2.7) introduces the well-known convolution integral being defined as:

0 d ( ) ( ) d ( ) ( ). (0) ( ) d d t g f t g t f t g f t          (2.9)
By considering the following boundary conditions on the tunnel's surface and at the far-field:

1 0 0 ( ( ), ) ( ) ( ) ( , ) r h r R t t p t p t t p           (2.10)
It can be deduced the expressions of radial and hoop stresses in the rock mass (Wang et al. 2013):

    2 2 0 1 0 2 2 2 2 0 1 0 2 2 ( ) ( ) 1 ( ) ( ) ( ) ( ) 1 ( ) ( ) h r h R t R t p p t p t r r R t R t p p t p t r r                                 (2.11)
Correspondingly, the radial displacement in the rock mass reads:

2 0 1 0 0 1 ( , ) ( ) ( ) ( ) ( ) 2 t h r u r t p p p R H t d r                (2.12)
where:

1 1 ( ) . ( ) H t L s G s         (2.13) With ( )
G s is the Laplace transform of the shear relaxation modulus ( ) G t .

The expression ( ) H t for the Burgers viscoelastic model is determined as [START_REF] Wang | An analytical approach for the sequential excavation of axisymmetric lined tunnels in viscoelastic rock[END_REF]:

1 1 1 ( ) exp ( ) K M K K M G H t t t G               (2.14)
For lined tunnels, the isotropic elastic behavior is chosen so that the radial displacement (urL), the radial-and tangential-stresses (rL, L) considering the supports' pressures on their boundaries can be easily deduced. Let t1, t2 are the installed time of the first and second lining, respectively. One can obtain the desired equations for the first lining at instant time t (t  t1):

  2 2 2 2 1 2 1 1 1 2 2 1 1 2 2 2 2 2 1 1 2 1 1 2 1 ( ) ( ) 1 ( , ) . ( ) ( ) . . 2 . L rL L L R R R p t R p t u r t p t p t r G r R R K R R          (2.15) 2 2 2 2 1 2 2 1 1 2 1 2 2 2 2 2 2 1 2 1 2 ( , ) 1 ( ) 1 ( ) rL R R R R r t p t p t r R R r R R                    (2.16) 2 2 2 2 1 2 2 1 1 2 1 2 2 2 2 2 2 1 2 1 2 ( , ) 1 ( ) 1 ( ) L R R R R r t p t p t r R R r R R                     (2.17)
Similarly, for the second liner, the radial displacement and stresses at instant t (t  t2)

2 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2 3 2 2 3 1 ( ) 1 ( , ) . ( ) . . 2 . L rL L L R R R p t u r t p t r G r R R K R R        (2.18) 2 2 3 2 2 2 2 2 2 2 3 ( , ) 1 ( ) rL R R r t p t r R R            (2.19) 2 2 3 2 2 2 2 2 2 2 3 ( , ) 1 ( ) L R R r t p t r R R             (2.20) Herein 1 R (i.e., fin R ) and 2
R indicate the exterior-and interior-radius of the first liner; 3 R means the inner radius of the second liner (as in Figure 2-1). The relation of shear and bulk moduli of each tunnel support and Young's modulus, as well as Poisson ratio, reads:

/(2(1 )), /(1 2 ), ( 1,2)        Li Li Li Li Li Li G E K E i (2.21)

Development of the closed-form solution

This section briefly represents the process to deduce the closed-form solution of the interesting problem, while more details can be found in [START_REF] Do | Time-dependent reliability analysis of deep tunnel in the viscoelastic Burger rock with the sequential installation of liners[END_REF].

The pressures of each liner are acquired from the compatibility conditions of displacement at the interface between the rock mass and the first liner (p1) and at the interface of two support liners (p2). Clarify, the continuity of displacements condition between the first liner and the rock mass is expressed as:

1 1 1 1 1 1 ( , ) ( , ) ( , ), ( ) r r rL u R t u R t u R t t t    (2.22)
Similarly, the compatibility of the displacement at the interface of two liners imposes:

1 2 1 2 2 2 2 2 ( , ) ( , ) ( , ), ( ) rL rL rL u R t u R t u R t t t   
(2.23)

These equations (Eqs. (2.22) and (2.23)) can be rewritten with the form of the standard integral equation (see [START_REF] Do | Time-dependent reliability analysis of deep tunnel in the viscoelastic Burger rock with the sequential installation of liners[END_REF] for this transformation) as follow:

  1 2 ( ) ( ) 1 2 ( ) ( ) ( ), x x t x t a y x E e E e y t dt f x         (2.24)
The solution of this equation can provide the supporting pressure of the first liner as (for the sake of simplicity, see the detailed developments in [START_REF] Do | Time-dependent reliability analysis of deep tunnel in the viscoelastic Burger rock with the sequential installation of liners[END_REF]):

  1 11 1 2 ( )exp , B K K G t t p t t t t             (2.25)   2 12 2 ( )exp , B K K G t p t t t t            (2.26)
as well as the supporting pressure of the second liner: Once obtaining supporting pressure   1 p t , the radial displacement (rock mass) can be identified by using the Eq. (2.12); more specifically, this displacement is indicated according to time periods as follows:

1 0 2 0 1 3 1 2 4 2 ( , ) ( , ) ( , ) if ( , ) ( , ) r r r r r u r t t t u r t t t t u r t u r t t t t u r t t t              
(2.28) For 0 t t  (i.e., the excavation stage), by substituting 1 ( ) 0 p t  and Eq.(2.14) into Eq.(2.12), the displacement reads:

2 0 1 1 1 1 1 1 ( , ) ( )( ) exp ( ) ( ) , 2                         h K r ini r M K K M p G u r t t R v t t D t G t r G
(2.29) where:

  2 2 2 1 1 0 0 ( ) ( ) ( )exp 1 .exp exp t t l K K ini r K ini r K m v G G D t R d m R v d R v                                               (2.30)   2 2 2 1 1 0 0 ( ) ( ) ( ) 1 .exp , t t l ini r ini r m v G t R d m R v d R v                               (2.31)
For 0 1 t t t   , (i.e., after the excavation and before the first liner's installation), the displacement is now derived by:

2 2 0 2 1 0 2 2 1 0 2 1 1 1 ( , ) ( ) exp ( ) exp . ( ) 2 1 1 ( ) ( ) ,                                   h K K r fin fin M K K K K fin M M p G G u r t t R t D t t R D t r G G t R G t
(2.32)

where:

0 0 2 2 1 ( ) ( )exp 1 . xp exp t t l K K t t K fin K m v G G D t d m e d R                                            (2.33) 0 0 2 2 1 ( ) ( ) 1 . xp t t l t t fin m v G t d m e d R                            (2.34)
For 1 2 t t t   (i.e., after the first liner's installation and before the second liner's installation) the displacement reads:

  2 2 0 0 3 0 11 1 0 2 2 2 2 0 0 3 1 0 2 3 1 1 ( , ) ( ) ( ) exp . ( ) exp ( ) 2 1 1 exp . ( ) ( ) ( ) . ( ) h h h K K r fin fin M K K K K h h K fin fin fin K K M M M p p G G u r t p t p t R t D t R t D t r G p p G t R D t G t R G t R G t                                                
(2.35) where:

1 3 11 ( ) ( )exp t K t K G D t p d             (2.36) 1 3 11 ( ) ( ) t t G t p d    
(2.37) For 2 t t  (i.e., finishing the second liner's installation), after some manipulations with using the support pressure acting on the rock 12 ( ) p t (as in Eq. (2.26)), the displacement is deduced: We choose the referent analytical solution introduced in [START_REF] Nomikos | Supported axisymmetric tunnels within linear viscoelastic burgers rocks[END_REF], where the time elapsed between excavation and installation of the support, as well as the effect of support deformation on the time-variant response, were investigated. We keep the hypothesis adopted from Nomikos paper, that is, the excavation of the tunnel is instantaneous, or the tunnel face advancement and the stress history before tunnel construction is not counted. The Nomikos solution can be adopted as a specific case of our advanced solution (by taking the parameter 1 0 m  in Eq. (2.4)). Furthermore, our solution, mainly used for calculating the support pressure 1 ( ) p t and radial displacement 1 ( , ) r u R t ), highly agrees with [START_REF] Nomikos | Supported axisymmetric tunnels within linear viscoelastic burgers rocks[END_REF], as shown in Figure 2-3 below by applying a constant of Ks,, a stiffness of liner used in [START_REF] Nomikos | Supported axisymmetric tunnels within linear viscoelastic burgers rocks[END_REF]) as follows:

  2 2 0 0 4 0 12 1 0 2 2 2 2 0 0 3 2 4 1 0 2 2 3 2 1 1 ( , ) ( ) ( ) exp ( ) exp ( ) 2 1 1 exp ( ) exp ( ) ( ) ( ) 1 ( ) h h h K K r fin fin M K K K K h h K K fin fin fin K K K K M M fin M p p G G u r t p t p t R t D t R t D t r G p p G G R t D t R t D t G t R G t R G t                                                        2 4 1 ( ) fin M R G t     (2.38) where: 2 4 12 ( ) ( )exp t K t K G D t p d             (2.
  2 2 1 2 2 2 1 1 2 1 1 1 2 s L L L R R K R R K G     
(2.41)

Notably, in these studies, the lining is installed at t=t1=30 (days) after the tunnel excavation. Notably, the initial conditions (namely, the boundary conditions and the sequential excavation and installation of liners) are kept as in the mentioned analytical simulation. That means on the top-and right-boundaries. We impose the hydrostatic stress 6.8(MPa) while on the bottom-and left-ones, the fixed normal displacement is applied. The quasi-instantaneous excavation is investigated with the initial time t0=5(minutes) while the two linings are put at t1=3(hours) and t2=2(days), respectively. The tunnel radial displacements in the vicinity and in the rock mass are time-variant functions (as in Figure 2-6 below). These figures show good agreements between the two approaches. One remark also the concordance of stress state calculated in the two liners as highlighted in Figure 234567. Clearly, the creep behavior of time-dependent Burgers rock causes compressive radial stress (r) as well as tangential stress () in both the linings. The figures show that considering the liners' interior to the exterior surface, radial stresses significantly raise while the tangential ones slightly decline. Notably, the maximum equivalent stress (q=|-r|) on every lining is located on its inner surface. In these investigations, our proposed analytical solution differs from numerical ones by less than 3.5%. The influence of various parameters on the tunnel's reactions can be investigated by using parametric studies. These parameters can be a face advancement rate of the tunnel, (or tunnel speed) an installation time, and thickness of liners as well as the viscoelastic Burgers rock properties (i.e., the first four parameters in Table 234567). For this purpose, the interested parameter is varied during the investigation, while all other ones assume to be constant.

In the first part, the tunnel's face advancement rate, the installation time, and the thickness of liners are firstly taken in to account for its effects on the tunnel behavior. These parametric studies are studied along with using the mean values of the four viscoelastic rock properties. Figure 2-8a strongly illustrates a faster excavation (i.e., a higher speed vl) conducts to a more evident tunnel convergence over time. This result is expected when the high rate stimulates a fast decrease of the initial pressure and thus a more critical radial displacement at the tunnel surface. Also, the high excavation rates can increase the equivalent stress in the second liner; yet it can reduce in the first lining (Figure 2-8b, c). However, the effect of the excavation speed on the behavior of coatings is not clear, particularly at the moderate rate of tunnel advancement about vl =1(m/day) and for initial liners installation time. For illustration, when the excavation rate varies from 0.75(m/day) to 4.5(m/day) and with t1=1(day) and t2=5(days), the maximum equivalent stress in the first liner is found with vl=0.75(m/day) (as in Figure 2-8a-center) while the minimum equivalent stress in the second one is obtained with vl=1.5(m/day) (as in Figure 2-8a-right). Thus, the complex evolutions of the stresses in the linings are highlighted. Notably, the liners must withstand at the identical time the increased radial displacement over time of the rock mass and the initial residual stress state on the tunnel surface.

The installation of liners is proven that the earlier time can reduce the convergence of tunnel (Figure 2-9a) by using, without loss of generality, the excavation rate vl=0.75(m/day) and the considered time t=100(years). It is especially more critical with the earlier putting of the first liner generally makes a higher value of equivalent stress. Notably, installing the second liner as early as possible reduces the stress in the first liner but results in higher compressive stress in the second lining.

For considering the thickness of the liners, Figure 2-9b shows a common trend of the tunnel convergence (i.e., radial displacement) and equivalent stress being reduced if a greater depth of each lining is applied. Note that the compressive stresses are analyzed at the inner surface (i.e., the maximum values in the isotropic and homogeneous cases). In the second part, the impact of the four viscoelastic Burgers rock properties on the behavior of the tunnel is assessed. Figure 2-10a shows that the effect of the dashpot of Maxwell element (M) only obviously appears in the long-term. Specifically, the higher value of M denotes a reduction of the tunnel convergence and the equivalent stress in the linings. Figure 2-10b illustrates the effect of dashpot of the Kelvin element (K). Subsequently, a higher viscosity coefficient (K) can slightly reduce the convergence of the tunnel but raise the equivalent stress in the second liner. Figure 2-10c demonstrates that the influence of the spring of the Maxwell element (GM) on the tunnel convergence. This is undoubtedly significant, but the equivalent stress of the first liner is negligible. Particularly, a greater value of this coefficient gives rise to the compressive stress in the second lining. Finally, Figure 2-10d proves that the Kelvin element (GK) effect is the most significant among the other Burgers rock parameters. Undoubtedly, the convergence and the equivalent stresses in the liners can be radically reduced with the higher coefficient (GK). The uncertainty of the four viscoelastic Burgers rock properties (the random vector X = [GM, M, GK, K]) is mainly taken into account for the reliability analysis. In other words, the other uncertainty sources (i.e., the tunnel geometry, hydrostatic stress) are not considered. By means of the reliability analysis, the effects of different parameters (i.e., the advancement rate, the time of installation, and the thickness of the liners) on the exceedance probability of the tunnel are investigated.

The analysis is considered with two failure modes being presumably related: the tunnel convergence exceeds an allowable value and the collapse of liners when the stress state overpasses have been tolerated. Especially, the radial deformation on the boundary of the tunnel is assumed to be limited to a maximum value lim. Therefore, the response function of the first failure mode, GTunnel, reads:

1 lim 1 ( , ) ( , ) , X X    R Tunnel u t G t R (2.42)
Additionally, the LSF of the second mode for each liner GL1, GL2 are expressed in the form: The value for each parameter will be systematically applied in our reliability analysis if the supplementary information is not given herein. The four viscoelastic Burgers rock properties are identical to those in the section mentioned above (see the name and its mean values in Table 2345678). In our current analysis, the uncertainty of these properties regards the mean values with an assumption of Gaussian distribution. Precisely, the parameters are entirely identified by adopting the identical coefficient of variation (COV) equals 0.25 (i.e., its standard deviations are assumed to be 25% of its mean values). Besides, the thickness and the installation time of two liners are chosen respectively l1=46(cm), l2=20(cm), t0=1(day), t1=2(days), t2=5(days), and the excavation rate vl=0.75(m/day). The other parameters, such as the radius of the tunnel, elastic properties of both linings, and the initially applied stresses, are also listed in Table 2-8. It is noteworthy that over the time considered, a monotonic increase in tunnel convergence and the compressive equivalent stresses in the two liners is confirmed in the previous part. Consequently, according to the selected criteria of failure mode, the time-variant reliability analysis can be investigated similarly to the time-independent problem. In other words, any time-independent probabilistic method can easily be applied here to estimate the exceedance probability for a chosen life of the tunnel. Thus, on the basis of the developed closed-form solution, the direct MCS will be conducted to measure the probability of exceedance of the structure. After that, MCS's obtained results will be used as the reference solutions to investigate the accuracy of the new Kriging metamodeling techniques.

1 lim1 1 1 1 1 ( ) ( , ), with ( , ) ( , ) ( , ) X, X X X X         L L L L rL G t q t q t t t (2.43) 2 lim2 2 2 2 2 ( ) ( , ), with ( , ) ( , ) ( , ) X, X X X X         L L L L rL G t q t q t t

Reliability analysis by using MCS

Like the previous parametric study (i.e., section 2.3.5), now we consider the effects of the face advancement rate, installation time as well as the thickness of both linings into reliability analysis problems. For estimating the probability of exceedance, the additional parameters such as the thresholds (i.e., acceptable values) of the convergence and the equivalent stress need to be given. There some arbitrary values are chosen for illustration. Precisely, the allowable convergence of tunnel, the stresses in the first-and the second liner respectively are lim 1.5%   , lim1=40(MPa) and lim2=50(MPa). In some particular parametric studies, these parameters lim lim1 lim2 , ,    are also considered as variants.

Firstly, we take the evaluation of the tunnel deformation exceeds the allowable value lim 1.5%   into account. Following that, Figure 2-12c proves that the higher lim is, the lower the probability can be. Secondly, we consider the exceedance probability of the stresses on the liners (as in Figure 2-13 to Figure 23456789101112131415). These figures obviously show an increasing probability of exceedance overtime and a consistency with the previous parametric analyses. In particular, a crucial excavation tunnel speed leads to a higher exceedance probability in the second lining (Figure 2-13b). Conversely, Figure 2-13a represents a higher probability in the first lining with a lower advancement rate. Besides, for the case of t1=t0=1(day) (i.e., the first liner is put immediately after the excavation), (Figure 2-13c), the exceedance probability reaches the peak at the excavation rate vl=1(m/day) instead of vl=0.5(m/day). In engineering practice, the second lining is usually installed to ensure the long-term stability of the tunnel. At the same time, the first one is placed as soon as the completion of the tunnel excavation to seal the rock. However, with a selected excavation rate, sooner installation of each lining could raise the exceedance probability (Figure 2-13g, h). Figure 2-13g presents the exceedance probability on the first lining may be low if the second support installation time is as short as possible. Nevertheless, this situation induces a higher probability of exceedance on the second lining (Figure 2-13h). The impact of the thicknesses of liners is now investigated by retaining the selected tunnel advancement rate as well as the installation times (i.e., vl=0.75(m/day), t1=2(days), and t2=5(days)). According to the results shown in Figure 2-14 and Figure 2-15, a slight variation in the liners' thickness about a few centimeters can create a surprising change in each lining's exceedance probability. More specifically, an increase of about 3(cm) in the thickness of each liner can reduce the exceedance probability in the first liner (Figure 2-14b) as well as in the second element (Figure 2-15b). This strong decrease in the probability is consistent with a remarkable reduction in compressive stresses in both elements in terms of thicknesses, as illustrated in the previous parametric studies (Figure 2-9b, c). Lastly, considering the smaller allowable values of compressive stress in each liner, a greater probability of exceeding the maximum allowable stress can be determined, as expected, if the liner thicknesses are maintained (Figure 2-14c, d, and Figure 2-15c, d). This means that the thickness of each liner must be increased to achieve the same probability of exceeding the initial limited value of compressive stress in each liner. There may be different combinations of varying liner thicknesses, yet the optimized values can be obtained when the other constraints (e.g., chosen technology and cost-effectiveness) are provided. In this work, the modified AK-MCS will be proposed for the reliability analysis of deep tunnels excavated in the squeezing rocks. This modification of the well-known U learning function in the initial AK-MCS [START_REF] Echard | AK-MCS: An active learning reliability method combining Kriging and Monte Carlo Simulation[END_REF] aims improving more efficiency. The potential gain of energy is considerable for the practical engineering application, especially for the complex structure system such as the stability of drift constructed in the Callovo-Oxfordian claystone (COx) for the deep disposal of radioactive waste. In this framework, the problem of long-term stability analysis accounting for the time-dependent mechanical behavior of rock mass is of most importance, in which solving the deterministic problem is expensively time-consuming by using numerical simulation software.

As mentioned in Chapter 1, the capability to choose the most suitable new training samples to improve the stabilization of the exceedance probability with a smaller size of DoE (NDoE) characterizes the learning function's efficiency. According to (Hu and Mahadevan, 2016a), a new training point is chosen, such that it reduces most significantly the uncertainty of the exceedance probability instead of choosing the lowest U value (in AK-MCS) or the highest EFF value (in EGRA). A new learning function called Global Sensitivity Analysis enhanced Surrogate (GSAS) was presented by (Hu and Mahadevan, 2016a). The efficiency of this method was demonstrated. However, following these authors, the GSAS method may increase the computational overhead required by the algorithm selecting the training points.

Besides, it was highlighted in many 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 both close to the limit-state surface (LSF) and far away from the training samples of existing DoE by verifying a distance constraint. Regarding [START_REF] Zhao | A novel active learning framework for classification: Using weighted rank aggregation to achieve multiple query criteria[END_REF], the following distance constraint dmin>D was proposed in which the limit distance parameter D varies during the add-sample process, whilst dmin is the minimum space from the candidate point to training points of DoE.

In what it follows, the modified U learning function is proposed, so-called Modi-U, which considers the distance constraint. Then, the validation of this modified AK-MCS on the reliability analysis of structures will be demonstrated through a simple academic problem (e.g., a two-dimensional series system with four branches), and for the case of deep tunnel excavated in the linear viscoelastic rock.

Distance constraint U learning function

The modified U learning function proposed in this work aims to gather a new suitable training point that verifies both the conditions: near to the LSF and far away from the training points of existing DoE. More precisely, the chosen new training point x* corresponds to the lowest value of the vector U(x) of the NMCS random samples, which verifies the following distance constraint:

, 1,2,...,

   x* x i DoE D i N
(2.45)

The limit distance parameter, i.e., D in Eq. (2.45), can vary after each iteration and is defined:

      int ( ) 1, min max , 1, 2,...,     i j i N D j d (2.46)
The lowest U value (min{U(x)}) may not be chosen as the new training point if it does not satisfy the distance constraint in Eq.(2.45). Instead, as the principal idea, this modified U learning function aims to find, among the points closest to the limit state, the new training point that is far enough away from the existing training points of the current DoE by checking the distance constraint in Eq.(2.45). Following that, the minimum value among the highest hyperparameter (i.e., among the maximum of vector ) calculated from the initial to the current iteration is chosen as the flexible limit distance parameter D. One can state that as a function of the increased iterative number, this limit distance parameter D decreases. A higher value of D at the first iterations allows us to seek the new training point in the more significant input variables' space instead of using the local value obtained from argmin{U(x)} which may have minor effects on the variation of the Kriging metamodel and hence the exceedance probability.

Note that the distance constraint is also applied for the case that a subset of new training points is taken for each iteration (i.e., the distance between these added points must also verify this condition). In addition to this modified U learning function, the stopping criterion is written in Eq. (1.47) is considered in this work. Our numerical applications in the next subsections showed that a value =0.01 and N =6 is enough to obtain the convergence of exceedance probability. Herein, for all numerical applications, the number of random samples for the direct MCS (used as the referent results) and for the interpolation from the constructed metamodels is taken at 10 6 (NMCS).

Application to a series system with four branches

As the first example, the efficiency of the modified AK-MCS is investigated by considering the following well-known academic performance function (also called a series system with four branches (see [START_REF] Waarts | Structural reliability using finite element methods: an appraisal of directional adaptive response surface sampling[END_REF][START_REF] Echard | AK-MCS: An active learning reliability method combining Kriging and Monte Carlo Simulation[END_REF][START_REF] Vahedi | An adaptive divergence-based method for structural reliability analysis via multiple Kriging models[END_REF][START_REF] Jiang | A general failure-pursuing sampling framework for surrogate-based reliability analysis[END_REF]):

2 1 2 1 2 1 2 1 2 2 1 2 1 2 2 1 6 3 0.1( ) ; ( ) 2 2 ( , ) min 6 3 0.1( ) ; ( ) 2 2                            x x x x x x G x x x x x x x x
(2.47)

where x1 and x2 are two independent random variables.

In order to study the effect of the chosen initial DoE on the evolution of exceedance probability during the iteration process, the initial training points are generated in two different ways: random by Gaussian distribution and quasi-uniform, such as Latin Hypercube Sampling (LHS).

Similarly, as in the study of [START_REF] Echard | AK-MCS: An active learning reliability method combining Kriging and Monte Carlo Simulation[END_REF], in this academic example, an initial DoE of twelve training samples is taken.

With the Gaussian first DoE generation, we present in Figure 2-16 and Figure 2-17 the results of the evolution of exceedance probability Pf (i.e., probability of G(x1,x2) <0) in which the classical AK-MCS and EGRA methods are conducted. These results agree well with the previous discussions of [START_REF] Echard | AK-MCS: An active learning reliability method combining Kriging and Monte Carlo Simulation[END_REF]. Indeed, as illustrated in Figure 2-17, before attaining the convergence of Pf; the evolution of exceedance probability during iteration develops at different levels (four levels). Moreover, the number of repetitions at each level is quite essential, and it seems risky to apply in this case, the stopping criterion defined in Eq.

(1.47). A small number of N may induce a lower exceedance probability evaluated at an early level. With respect to the position of the new learning points relative to the limit state during the iteration (Figure 2-17), it appears that they were adjusted entirely locally as the added samples follow each branch of the four-branch system. Besides, fairly large new training points, which are not close to the LSF, were selected after the convergence of probability.

A remarkable difference can be observed for the case that the initial training points of DoE are generated quasi-uniform by the LHS method (Figure 2-18). Regarding the evolution of exceedance probability during iteration by the AK-MCS method, the number of levels, as well as the number of iterations in each level, reduces importantly with respect to the previous case of randomly initial DoE. Concerning the result obtained from the EGRA method, a notable dispersion of the exceedance probability can be stated before the convergence, and no level was developed in this considered case. From the results highlighted in (Figure 2-18), it is evident that the number of iterations to attain the convergence reduces in this second case by using the quasi-uniform initial DoE. The comparison of the classical AK-MCS and EGRA methods does not highlight an essential difference in terms of iteration number to convergence. The difference seems more pronounced regarding the position of added training points concerning the limit state. While these samples always follow each branch of the four branches system in the EGRA method, the new training points selected by AK-MCS are quite random around four-branches during iteration.

We now examine the efficiency of the modified AK-MCS method by using the distance constraint U learning function. For this purpose, the results obtained from this last method are compared with those of the classical AK-MCS method. These results were conducted on the quasi-uniform initial DoE by accounting for its advantages, as observed previously. Thanks for using the distance constraint, no level of probability Pf during the iteration was generated before convergence in the modified AK-MCS method (Figure 2-18). This can be clarified by the fact that the chosen new training samples in each iteration are far enough with respect to the training points of the existing DoE, which have more effect on the variation of the probability. Note that this academic problem shows that the limit distance parameter D remains constant during iteration (D equals 0.3536). The disappearance of the local convergence of exceedance probability in each level (the case was noted in the original AK-MCS method) allows applying the stopping criterion (see Eq. (1.47)) in our modified AK-MCS method. For example, by taking =0.01 and N =6, the probability of exceedance measured by the modified AK-MCS method is similar to the one provided by the classical AK-MCS method, but the convergence is much more sooner in the former method. Concerning the position of the new training points of DoE during iteration concerning the limit state, we observed that they were selected around the four branches system, similar to the classical AK-MCS method (Figure 2-18b). Finally, Table 2-8 below illustrates the reliability results of the academic four-branches system problem: comparison of EGRA, classical (U), and modified AK-MCS (Modi-U) methods with the direct MCS as well as the final number of iterations for each method. Note that Ncall corresponds to the number of calls to the response function, which is equal to the sum of the number of initial training points and the number of iterations to attain the convergence. The comparison with the referent method (MCS method) confirms the accuracy of these metamodeling techniques. In this section, the results of the previous MCS reliability analysis of tunnels excavated in the Burgers rock can be used as a reference to prove the efficiency of our modified AK-MCS method. For this aim, we adopt the same hypothesis of the previous reliability analysis by considering only the uncertainty of four parameters of the Burgers rocks on the exceedance probability of tunnel. However, to simplify the presentation, we only utilize the failure mode of the second liner when the stress state in this support element exceeds the permissible stress of the constitutive materials. The corresponding performance function of this failure mode is defined in Eq.(2.48).

All the necessary parameters included in the present reliability analysis with the AK-MCS are taken similarly to ones in the previous subsection (2.3.1 and 2.4.1). That means the four Burgers rock properties are considered as random variables with a standard deviation 25% of the mean value (i.e., COV, coefficient of variation is equal to 0.25). The thickness and the installation time of first-, second-lining are chosen respectively l1=46(cm), l2=20(cm), t0=1(day), t1=2(days), t2=5(days), while the excavation rate is vl=0.75(m/day). The allowable stress on the second liner lim2 =50(MPa) is chosen. Note that, for the reliability analysis by metamodeling, an initial DoE with 24 training points generated quasi-uniform by the LHS method is used.

Figure 2-19a shows the evolution of the exceedance probability of the stress on the second liner depending on the number of iterations. These results calculated from three metamodels (i.e., EGRA (EFF), classical (U), and modified AK-MCS (Modi-U) methods) represent in effect the probability measured at t=100(years). Like the previous study cases, before attaining the convergence, the Pf probability estimated from the EGRA and classical AK-MCS methods develops in two levels. Contrary, by using the distance constraint in the modified AK-MCS method, the apparition of these two levels was not observed. Subsequently, the utilization of a dynamic distance parameter D, which is higher in the initial iterations and reduces during the DoE enrichment process , contributes a significant effect to avoid the apparition of the local convergence of exceedance probability. A much smaller number of iterations to attain the global convergence of the Pf probability in the modified AK-MCS method with respect to one of the two other methods confirm an improvement of our proposed Kriging metamodel.

Meanwhile, the accuracy of the three metamodels is also demonstrated in Table 2-9 below. This table shows the reliability results of the tunnel excavated in the viscoelastic Burgers rock with a comparison of EGRA, classical (U), modified AK-MCS (Modi-U) methods, and the direct MCS. Note that the exceedance probability is also measured at t=100 years. Finally, Figure 2-20 illustrates the effect of the number of new training points in a subset to be added in each iteration in the DoE on the convergence of the exceedance probability by the modified AK-MCS method (Modi-U). This investigation is especially interesting for the many projects in which the evaluation of the performance function needs to be conducted through the numerical simulation with the possibility of realizing many calculations at the same time (i.e., parallel computations). The results illustrated in Figure 2-20 admit that, in general, an increase in the number of new training points in each iteration will reduce the total number of iterations at convergence (with the same adopted stopping criterion). However, this last number can reach an asymptotic value when this number of new training points in the added subset of DoE becomes higher. In this study case, a subset of four new training points in each iteration can be an appropriate choice. Indeed, although the number of calls to the response function increases significantly in this last case due to the higher number of additional training points selected, the actual time required for the numerical simulation decreases. The reduction here is that the number of iteration to achieve convergence of exceedance probability is lower.

Concluding Remarks

In this chapter, the reliability analysis of deep tunnels excavated in complex geological conditions was conducted. In this context, a sequential excavation and installation of more than one-liner are necessary to avoid collapse and ensure the safety of structures in both the shortand long-term. To achieve this aim, we proposed in the first stage, a closed-form solution of double-lined tunnels in linear viscoelastic Burgers rock through which the direct Monte Carlo Simulation was then carried out in the next step. The numerical investigation highlighted the critical effects of the different parameters, notably the uncertainty of the rheological parameters of the Burgers rock on the long-term exceedance probability of tunnel. The MCS investigations are helpful when they can be served as a primary tool for the selection of suitable surrogatebased reliability methods, which can be more appropriate in the analysis of complex underground structure system.

The main challenge of the surrogate-based reliability analysis, like on the way to develop an adequate metamodel, is to approximate the LSF without wastes of samples and function calls. Among the various surrogates proposed in the literature, we chose the active learning Kriging metamodel technique. Because this technique not only provides the best estimate value at an untired point but also gives an estimation of the prediction variance. This variance presents the prediction uncertainty, which depends on the number and location of the training data points. This metamodel can be established iteratively through an enrichment process of the DoE. The discussions in the literature showed that the best candidate points to enrich the DoE should not only be close to the LSF but also be far away from the training samples of existing DoE. This challenge is a crucial motivation for us to propose a novel approach: a modified AK-MCS (Modi-U). The novelty of this proposed method lies in consideration of the constraint distance, representing by a parameter D, which can vary during the enrichment process, on the selection of the new training point of DoE. With this condition, the points chosen from the minimum of vector U in the classical AK-MCS may not be added to the DoE. The validation and significant improvement of this modified AK-MCS were demonstrated in two examples: an academic problem and the case of a tunnel excavated in the linear viscoelastic Burgers rock. Our investigations also highlighted the advantage of using a subset of new training points to enrich the DoE since it can lower the number of iterations to attain the convergence of exceedance probability. That is so useful for our studies in the next chapters in which the evaluation of the performance function must be performed numerically in an open-source code (e.g., Code_Aster) with available parallel calculations.

CHAPTER 3. KRIGING-BASED RELIABILITY ANALYSIS OF THE STABILITY OF A DEEP DRIFT SUPPORT FOR NUCLEAR WASTE DISPOSAL IN THE COX CLAYSTONE

Introduction

As demonstrated in the previous chapter, the accuracy of the modified AK-MCS allows us to apply this metamodeling technique to study the long-term stability of deep drift support in the French concept of nuclear waste disposal. In France, ), is in charge of studying the disposal of high-level and intermediate-level long-lived waste (HLW and ILW-LL) in a deep geological repository (Cigéo project). In this context, the stability of the concrete support of the underground structure during the operational period for about 100 years is crucial to ensure the functionality of the mechanical system within the tunnel.

Thanks to its favorable characteristics (incredibly low hydraulic conductivity, small molecular diffusion, and significant retention capacity for radionuclide), COx rock is taken into account as a potential geological formation for a deep geological repository (also disposal of radioactive waste) in France, if the Cigéo project is licensed (see Figure 3-1 below). Specifically, to study the feasibility of constructing and operating the Cigéo project, since 2000, Andra has been conducting an Underground Research Laboratory (URL) in Bure (at nearly 300 km East of Paris) of the Meuse and Haute Marne (M/HM) departments (Armand et al., 2017a;Armand et al., 2017b). To deeply understand the behavior of this argillaceous rock, intensive research programs have been performed combining laboratory tests, in-field observations, and numerical simulations. In parallel to many programs of laboratory tests (uniaxial/triaxial, mono/multi-stage creep and relaxation), in-situ observations have been undertaken for characterizing mechanical and hydromechanical short-, and long-term behavior of COx rock [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]2014;Vu et al., 2020b;[START_REF] Conil | In Situ Investigation of the THM Behavior of the Callovo-Oxfordian Claystone[END_REF]. Thanks to exhaustive laboratory tests and in-situ experiments in the URL, the essential fundamental aspects of the hydro-mechanical behavior of COx have been provided (see, for instance, the work of (Armand et al. 2017a)). Besides, the understanding of drifts' behavior with different excavation/support methods is the primary necessity for the support optimization design. At the M/HM URL, the drift network (as in Figure 3-2 (Armand et al., 2017b;[START_REF] Delay | Synthesis of 20 years Research, Development and Demonstration in Andra's Underground Research Laboratory in Bure for Cigéo Project-France[END_REF]) is constructed at two depths (at -490 m, the main level of the URL, for almost studies, and -445m for an experimental drift). At the -490 m, the anisotropic in-situ stress states have been deeply observed, and among the principal stresses' values, the major is the horizontal one; the vertical and the minor horizontal stresses are identical [START_REF] Wileveau | Complete in situ stress determination in an argillite sedimentary formation[END_REF]. The vertical component σv and the minor horizontal component σh (oriented 065°±10°) are approximately vh12.5MPa. At the same time, the major horizontal component σH (oriented 155°±10°) is larger than the others 30% (i.e., H1.3v 1.3h) [START_REF] Wileveau | Complete in situ stress determination in an argillite sedimentary formation[END_REF]. Regarding in-situ measurements, drifts' convergence data plays a vital role in the orientation excavated following the major and minor horizontal stresses. These convergences (monotonic functions of time) are differently pronounced for most substances. For example, vertical convergence is much more four times than the horizontal one with the drifts parallel to h, (see Fig. 3-2) [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]. Notably, the parallel drifts are often developed following the major in-situ stresses and are separated enough to be considered independent concerning their hydromechanical behavior.

The excavation of galleries at the main level (490 m depth) of the URL has been carried out essentially in the direction of the major and minor horizontal principal stresses. Continuous monitoring around drifts within MHM URL during and after the excavation showed anisotropic hydromechanical responses in term of excavation induced damage zone [START_REF] Armand | Geometry and properties of the excavation-induced fractures at the meuse/haute-marne URL drifts[END_REF], convergence [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], and pore pressure distribution (Vu et al., 2020b). Numerous different numerical models have been developed to explain these phenomena [START_REF] Pardoen | Modelling the influence of strain localisation and viscosity on the behaviour of underground drifts drilled in claystone[END_REF][START_REF] Trivellato | Modélisation en endommagement radoucissant de la zone de rupture autour des ouvrages profonds dans des roches anisotropes quasi-fragiles. Actes des 9èmes Journées Nationales de Géotechnique et de Géologie de l'Ingénieur[END_REF][START_REF] Souley | 3D anisotropic modelling of deep drifts at the Meuse/Haute-Marne URL[END_REF]Manica et al., 2021a, b;[START_REF] Yu | Numerical study of thermo-hydromechanical responses of in situ heating test with phase-field model[END_REF]. The time-dependent behavior of rock mass around drifts is also evidenced by following the evolution of the convergence of drifts' wall and the displacement within the ground [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] Delay | Synthesis of 20 years Research, Development and Demonstration in Andra's Underground Research Laboratory in Bure for Cigéo Project-France[END_REF] In reality, the time-dependent response of the ILW-LL cell can be governed by the thermohydro-mechanical 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 thermomechanical stress and the thermal pressurization phenomenon [START_REF] Vu | Thermo-poro-mechanical coupled processes during thermal pressurization around nuclear waste repository[END_REF]Seyedi et al., 2017;[START_REF] Conil | In Situ Investigation of the THM Behavior of the Callovo-Oxfordian Claystone[END_REF]Vu et al., 2020a;[START_REF] Bumbieler | Feasibility of constructing a full-scale radioactive high-level waste disposal cell and characterization of its thermo-hydro-mechanical behavior[END_REF]Plúa et al., 2021a, b). Moreover, the time-dependent 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, shrinkage/swelling of COx claystone due to the desaturation/resaturation (Wang et al., 2020, a, b). Nevertheless, 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 thermo-hydro-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 on 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 probability of exceedance 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](Armand et al., , 2017)). Advances numerical models have been proposed to reproduce the observation of short and long-term behavior 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] elastoplastic combining with weakness plan model [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), whose physical meaning and calibration are not simple to be evaluated, these sophisticated models are not feasible for the reliability analysis.

Moreover, the lining monitoring confirms the significant effects of the compressible material, which can be used to separate the concrete lining from the COx claystone (i.e., the material acts as the outer liner in the drift support system). The material can be made up of compressible grout, a compressible substitute for pea gravel, or a compressible layer on the backside of the precast segment lining [START_REF] Billig | DeCo Grout -Innovative Grout to cope with rock deformations in TBM tunnelling[END_REF]. The main reason for using this material is that the essential characteristic of a high void ratio (with the compression rates up to 50% of the initial volume under stress lower than 3 MPa are possible [START_REF] Billig | DeCo Grout -Innovative Grout to cope with rock deformations in TBM tunnelling[END_REF]) allows this compressible outer liner to absorb the convergence over-time of the host rock and considerably reduce the radial stress transmitted to the inner liner [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] Delay | Synthesis of 20 years Research, Development and Demonstration in Andra's Underground Research Laboratory in Bure for Cigéo Project-France[END_REF]. In recent times, this benefit provided by this compressible material has been widely tested in the configuration of vast drift excavation in combination with innovative construction methods [START_REF] Bosgiraud | Compressible Arch Segments for the Cigéo Access Ramps, Drifts and Vaults-A Field Test-17011[END_REF][START_REF] Stasiak | Experimental investigation of mode I fracture for brittle tube-shaped particles[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] Zghondi | Qualification, Construction and Analysis of a Precasted Compressible Arch Segments Drift Test in the Andra Meuse/ Haute-Marne Underground Research Laboratory[END_REF][START_REF] Delay | Synthesis of 20 years Research, Development and Demonstration in Andra's Underground Research Laboratory in Bure for Cigéo Project-France[END_REF].

As particular study cases, our objectives focus on studying the effect of the uncertainty of COx rock properties and the consideration of the compressible material on the long-term stability of the concrete liner of the circular drift. Although complex coupling mechanisms (such as coupled thermo-hydro-mechanical models and damage mechanisms) can govern the timedependent response of underground nuclear waste disposal structures, our current study focuses only on the pure mechanical behavior of the drift. Additionally, it is also assumed that drift's axis aligns with the major horizontal stress. The classical viscoplastic model of Lemaitre is chosen to describe the long-term behavior of host rock. This simple model's choice can be explained by its capacity to catch the essentials of time-dependent behavior of COx rock by a limited number of parameters to be calibrated, which significantly facilitates the uncertainty quantification process from the creep tests in the laboratory. The last simplification concerns the behavior of the compressible material of the outer liner, whose non-linear response is modeled as a tri-linear elastic law as a function of the volumetric deformation.

Precisely, to gain our objectives, the modified-Kriging-based approach is chosen to assess the exceedance probability of the final drift support during a period of one hundred years. In the following, the open-source Code_Aster in the 2D plane strain condition is used for the direct evaluations of the performance function (i.e., numerical simulation of the deterministic problem). This task can start from quantifying the uncertainty propagation of the random input variables (i.e., the short-and long-term mechanical properties of COx) on the exceedance probability (Pf) of the final lining of the drift. Throughout the numerical investigations, the compressible layer's role and the uncertainty of COx rock properties on the long-term stability of the final drift lining are highlighted.

Description of the studied problem

Assumptions and preconditions

Although there has been considerable development of three-dimensional (3D) analysis software packages using FEM or Finite Difference algorithms, 2D analysis is still needed as the main engineering tool for practical analysis of tunnel behavior and tunnel support performance in design [START_REF] Vlachopoulos | Appropriate Uses and Practical Limitations of 2D Numerical Analysis of Tunnels and Tunnel Support Response[END_REF]. In which the convergence is usually monitored by boundary displacement or internal pressure relaxation (convergenceconfinement method-CCM). Because the latter allows us to appreciate the nature of the interaction between the rock-mass and the support and the effect of the variation of the rock properties on the support, we decide to use it through the thesis. CCM, also called the stress relief method, is commonly applied for preliminary circular tunnels' assessment of squeezing potential and support requirements in a variety of stress states and geological conditions. A critical key so-called rate of convergence-confinement or deconfinement rate (plays an essential role in this method. Notably, its values applied in each phase are taken from Andra (see Figure 3-4a). This model is often related to the tunnel face progression, and this could be shown as in Figure 3 -3 (El Matarawi, 2016). The readers can figure it out more in many recent articles, such as [START_REF] Alejano | Why the convergence confinement method is not much used in practice[END_REF][START_REF] El Matarawi | Investigating the use of reliability-calibrated partial factors in tunnel design[END_REF][START_REF] Paraskevopoulou | Analysis of time-dependent deformation in tunnels using the Convergence-Confinement Method[END_REF]. The lined tunnel excavation in the rheological rock with the implementation of CCM requires some crucial assumptions, as below:

(1) The tunnel (drift) is of circular section and horizontal axis.

(2) This tunnel is a deeply buried tunnel (H>10*R) so that no linear variation of the stresses as a function of depth must be considered. (3) The rock is isotropic, homogeneous. (4) It is a 2D plane strain mechanical problem and isotropic initial stress state (0=x=y) (5) The excavation speed is low enough (no dynamic stresses). ( 6) It considers quasi-static equilibrium (i.e., not in terms of acceleration).

Moreover, supplementary assumptions are also necessary to address our problems in a more practical and feasible to Andra's project, as follows:

(1) The work is limited to COx claystone, with its rheology suitably described by viscoplasticity.

(2) The tunnel is excavated following the major horizontal stress in the COx.

(3) Only individual drift behavior is taken into account (i.e., parallel drifts are separated enough to ignore their hydro-/mechanical interactions). ( 4) The analyses are limited to available experimental data.

(5) The quantification of uncertainties is limited to the variation in rock properties. ( 6) Normal (Gaussian) distribution is used for the random variables (7) The installation time of double liners is corresponding to the confinement rate =0.85, and they are assumed to be installed immediately after excavation.

Problem statements

Our main considered problem consists of investigating the long-term stability of the drifts' final lining excavated in the major horizontal stress of the COx rock. To clarify step-by-step, it must first model the problem; for instance, Figure 3-4 and Figure 3-5 below illustrate the one-quarter of the adopted 2D plane strain model of drift and is taken the symmetric conditions into account.

To simulate the mechanical behavior of the circular drift with a 5.05m radius, a total dimension of 55m is considered following each direction that is sufficiently large to avoid boundary effects. Correspondingly to the symmetric boundaries, the normal displacement is fixed while isotropic far-field stress of 0=12.5MPa, (i.e., the identical stress on vertical (v) and minor horizontal (h) axis at a distance away from the excavation) is imposed. The considered drift is supported by a concrete lining C60/75 (with the thickness l2), which is separated from the rock mass by a compressible layer (with the thickness l1). This layer must absorb almost the convergence, which increases as a function of time due to the actual behavior of the COx claystone and reduce the radial stress transmitted to the concrete lining, which is supposed to be linear elastic.

The assumptions mentioned above and modeling of tunnel excavation with CCM can be classified into four main phases, as in Figure 3-5a. Phase 1: initialization of the constraints in situ (=0), computation of the nodal reactions on the level of the cavity (wall of the excavation). Phase 2: Deconfinement rock mass to simulate the progressive excavation and the distance of the tunnel face (=0). Phase 3: Installation of two liners at the same time. In this context, the deconfinement rate here values at =0.85. In general, it can be separated the instant time of various liners with the different rate Phase 4: End of the deconfinement with .

To investigate the long-term stability of the concrete lining, we use reliability analysis with the modified Kriging-based method combined Monte Carlo simulation (AK-MCS), detailed in Chapter 2. The performance function described the equivalent stress reaching the allowable pressure of the constitutive material (as in Eq. (2.48)) is recalled as the limit state definition (i.e., the safety and failure state in the concrete liner of drift). Then it will be applied in reliability analysis with the modified AK-MCS with the exceedance probability (Pf) taken from Eq. (1.23).

Constitutive behavior model of materials and uncertainty quantification of COx rock properties

In our study, three different constitutive models are addressed (a viscoplastic model for the geological rock mass, a non-linear elastic model for the compressible liner, and linear elastic for the concrete support liner). For the sake of simplicity, this study is also restricted to the purely mechanical response of drift by ignoring the other phenomena such as damage mechanism and/or coupled (thermo-) hydro-mechanical loadings. The classical viscoplastic model of Lemaitre is chosen to characterize the long-term behavior of host rock. The choice of this simple model can be explained by its capacity to catch the essentials of time-dependent behavior of COx rock by a small number of required input parameters to be calibrated (see in section 1.4.2, Chapter 1), which significantly facilitates the uncertainty quantification process from the creep tests in the laboratory. The last simplification concerns the behavior of the compressible material of the outer liner, whose nonlinear response is modeled as a tri-linear elastic law as a function of the volumetric deformation. Since the behavior of the final tunnel liner is linearly elastic, we just mention here the constitutive models for rock mass and for the compressible zone.

Time-dependent behavior of COx by the non-linear viscoplastic model of Lemaitre

To model the time-dependent behavior of COx claystone, one often uses constitutive models on the basis of elasto-visco-plasticity concept (Seyedi et al., 2017). Especially, the viscous model of Lemaitre is widely chosen to describe the long-term response of host rock. In general, these elasto-visco-plastic models can give the global behavior of the rock mass around drifts. However, to better explain the observed phenomena in-situ, more sophisticated models take the anisotropy effect in the plastic criterion and/or the damage mechanism, such as strain localization, into account. The interested readers can refer to (Seyedi et al., 2017), and different references cited therein for more details.

In our study, Lemaitre's viscoplastic model (see section 1.4.2) is chosen to present the timedependent mechanical behavior of the COx rock around the drift. Most of the authors dealing with the time-dependent behavior of COx claystone have proposed various (elasto-) viscoplastic models (Seyedi et al., 2017;[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 (see, for example, (Armand et al., 2017c;[START_REF] Mánica | A time-dependent anisotropic model for argillaceous rocks. Application to an underground excavation in Callovo-Oxfordian claystone[END_REF]Seyedi et al., 2017;[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. Based on the number of parameters to be calibrated in this model is so small, it can greatly facilitate the uncertainty quantification process and ensures the accuracy of the obtained results of these parameters.

Tri-linear elastic model of the compressible material

Using the compressible material as the outer liner to reduce the overstressing of the inner concrete lining of tunnel excavated in the squeezing rock has been largely considered in many research projects both in laboratory and tunnel construction programs [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. Tunnel boring machines--Trends in design & construction of mechanized tunnelling[END_REF][START_REF] Schneider | Complex support system[END_REF][START_REF] Billig | DeCo Grout -Innovative Grout to cope with rock deformations in TBM tunnelling[END_REF][START_REF] Cucino | Expanded clay in deep mechanised tunnel boring[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. Design and application at the Brenner Base Tunnel[END_REF].

At the M/HM URL, different types of compressible material have been tested, and the feedback from the experiment programs reveals their immense benefit. Thanks to the high compressibility (up to more than 50% of its initial volume as Figure 3-6a [START_REF] Billig | DeCo Grout -Innovative Grout to cope with rock deformations in TBM tunnelling[END_REF]), the compressible materials can drastically reduce the transmitted stress to the final liner of drift. From the experiments, the compressible material's stress-strain curve [START_REF] Billig | DeCo Grout -Innovative Grout to cope with rock deformations in TBM tunnelling[END_REF] can be captured in three stages. Following that, the curve illustrates an elastic behavior when the volumetric deformation of the compressible material is small. Then a plastic behavior representing by a strain hardening can be observed whilst the last stage characterizes the stiffening of the material after the depletion of its porosity. Some attempts to model the mechanical behavior of the high porous compressible materials have been conducted. For example, in [START_REF] Souley | A constitutive model for compressible materials: application to the study of interaction between supports and rock mass[END_REF] an elastoplastic constitutive model was proposed. In this model, the Drucker-Prager's criterion and a strain hardening modeled by an exponential function with respect to the plastic distortion are used. In addition, in this model, 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).

In this work, for the sake of simplicity, a tri-linear-elastic behavior was adopted for the compressible support layer (Figure 3-6b). This simplification is acceptable when in all our analyses of drift, the compressive stress state in this compressible liner presents a monotonic evolution. It is important to note that, under the monotonic loading condition, the stress-strain curve provided by our tri-linear-elastic model can fit quite well the experimental results of Billig [START_REF] Billig | DeCo Grout -Innovative Grout to cope with rock deformations in TBM tunnelling[END_REF].

As a function of the volumetric deformation, this model is characterizing by three values of elastic modulus (Ec1, Ec2, Ec3) and is expressed in the form:

1 1 1 1 2 1 2 2 1 2 1 2 1 2 3 2 3 3 2 2 ; ( ) 2 ; ( ) ( ) 2
+ +                               d d d I e σ I I e I I I e c v c v v c c v c v c v v v c c v c c v c v c v v K G if K K K G if K K K K K G if (3.1)
In Eq.(3.1) two parameters 1 2 ,   v v are the limits at the first/second stages and the second/third stages of the tri-linear elastic model, respectively. Notably, corresponding to the contrast of the mechanical property of the second stage with respect to the other stages, which are calibrated from the compressive tests of the compressible grout in [START_REF] Billig | DeCo Grout -Innovative Grout to cope with rock deformations in TBM tunnelling[END_REF]

  / (1 2 ), ( 1,2,3) / 2(1 )        ci ci c ci ci c K E i G E (3.2)
The Poisson's ratio of the compressible material is so small that it is taken equal to zero (c=0), i.e., ignore the impact of Poisson's ratio on the overall responses. This tri-linear elastic model has been implemented on the Code_Aster (Code_Aster® software) for this work. The results obtained from the tests in the laboratory and the in-situ measurements show a complex medium both in short -, and long-time of COx claystone. The claystone's short-term mechanical behavior is considered a linear behavior under low deviatoric stresses and investigated by means of uniaxial and triaxial compression tests. For low confining pressures, the failure of the samples is brittle and can be idealized by the formation of a shear band inclined with respect to the sample axis. Under high confining pressure, mechanical behavior is ductile [START_REF] Souley | A viscoplastic model including damage for argillaceous rocks[END_REF]. Both at the sample and in-situ scales, sophisticated phenomenological models combining advanced numerical modeling (non-local model, second gradient regularization, phase-field, so on) have been developed to reproduce the observation of short and long-term behavior of COx rock [START_REF] Pardoen | Modelling the influence of strain localisation and viscosity on the behaviour of underground drifts drilled in claystone[END_REF][START_REF] Souley | 3D anisotropic modelling of deep drifts at the Meuse/Haute-Marne URL[END_REF]Manica et al., 2021a;[START_REF] Manica | Numerical simulation of underground excavations in an indurated clay using nonlocal regularisation. Part 2: sensitivity analysis[END_REF]. Ordinarily, these models are characterized by various parameters (e.g., larger than 20), whose physical meaning and calibration are not simple. Consequently, the use of these models in a reliability analysis cannot be feasible.

The long-term one was studied by uniaxial creep and relaxation tests. For example, the results from the creep tests in the laboratory pointed out that the creep rate varies slowly and linearly with stress at low stresses, but it increases and deviates from the linear form under high loading (Chun-Liang [START_REF] Zhang | Thermo-hydro-mechanical behaviour of the Callovo-Oxfordian clay rock[END_REF][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] Zhang | Thermo-hydro-mechanical behavior of clay rock for deep geological disposal of high-level radioactive waste[END_REF]. Furthermore, the lack of measurable creep thresholds at the laboratory means that the creep of COx in their initial anisotropic stress state seems unlikely to be from a physical point of view. Also, depending on the mineralogical composition, the creep is slowed down with respect to the high carbon content.

The noteworthy uncertainty of the COx mechanical properties can relate to the scarcity of highquality data (epistemic uncertainty) and the inherent variability of the mineralogical composition of the host rock (the aleatory uncertainty).

This chapter focuses only on the epistemic uncertainty, while the consideration of the other uncertainty source (the aleatory uncertainty) contributes to another critical issue of Chapter 5. In our study, we take the second-moment statistics for Young's modulus according to the natural samples with UA geological unit in the contribution of Armand (Armand et al., 2017c), that is, the mean and the standard deviation worth of 6.0 GPa and 3.5 GPa, respectively. These statistical values of the short-term mechanical properties of COx are synthesized from the measurements conducted on natural samples obtained from the same depth at the main level of the M / HM URL (i.e., these samples have identical mineralogical composition). Note that the influence of the Poisson's ratio of rock mass on the behavior of drift is quite moderate with respect to the other parameters. Thus its mean value =0.29 is used in this work as a deterministic parameter by ignoring its uncertainty. Regarding the epistemic uncertainty of the long-term mechanical properties of COx, this topic has not yet been discussed in the literature. So far, a universal constitutive model to characterize the time-dependent behavior of COx accurately does not exist, and 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 epistemic uncertainty of the long-term mechanical properties of COx claystone, we use the data of seven tri-axial creep tests corresponding to seven applied deviatoric stresses performed on the samples taken from the same geological horizon (Armand et al., 2017c;Seyedi et al., 2017). The information about the confining pressure and deviatoric stress of each triaxial creep test is summarized in Table 3-10, while for more details of these tests, the interest readers can refer to (Armand et al., 2017c). The calibration of the viscoplastic parameters of Lemaitre's model is conducted by using the compact form expressed in Eqs.

(1.62) and the Best Linear Unbiased Estimator calibration algorithm (BLUE). Following these equations, the inverse analysis to determine the three parameters A, B, C can be carried out by fitting the data of at least two creep tests under two different applied deviatoric stresses. Since there are two creep tests using the same deviatoric stress (see Table 3-10), thus totally, the statistical analysis of the obtained results of three parameters A, B, C is conducted with 88 samples generated from at least two tests. Following the histograms in Figure 3-7, a normal distribution is adopted for these parameters characterizing the viscoplastic Lemaitre's model whilst their corresponding mean and standard deviation are summarized in Table 3-11. As a function of the applied deviatoric stress and time, the exponential parameters (i.e., the parameters B and C), which characterizes the non-linearity of COx rock, present the mean value about 2.05 and 0.17 correspondingly. Parameter A presents the most considerable uncertainty representing by the highest value of COV (about 87.5%). Our analysis also exhibits the correlation between these three parameters (see Table 3456789101112). More precisely, a quite significant anti-correlation between the coefficients A and B characterizing by a negative value -0.72 is stated. The other correlations (for example, between A and C or between B and C) seem weaker, which are, however, taken into account in our reliability analysis. Besides, with regard to Eq. (1.59), we can describe the cumulated viscous strain (as in Figure 3-9, by using the triaxial compressive tests describing the creep process as in Figure 3-8. These tests show different levels of applied axial and lateral stresses (i.e., different q/qpeak levels), as shown in Figure 3-8. The qpeak values were approximated from the triaxial compressive test results. Three q/qpeak levels (50%, 75%, and 90%) were applied and kept for three months. The cumulated viscous strain's evolution during the creep process is also illustrated in Figure 3456789. Obviously, the deviatoric loading level has an important influence on the amplitude of the axial strain. Note that in all creep tests that are carried out, the claystone is always deformed under pressure. As an example, in Figure 3-10, a tri-axial creep test is presented with the applied deviatoric stress 27.75(MPa) and confining pressure 12(MPa), while in Figure 3-10 are gathered all the simulated tri-axial creep tests. These results are simulated from the mean values and 1000 random sets of the three parameters A, B, C of Lemaitre's model. The comparison of the numerical result using the mean values of A, B, C, and the experimental data measured up to about 90 days (Figure 3-10a) 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 ten times the ones measured with the mean values (Figure 3-10). Notice 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 three months. Consequently, an extrapolation of the viscous deformation to 100 years(Figure 3-10b) must also present another source of uncertainty that may underestimate or overestimate the reliability analysis's corresponding results.

Discussions on deterministic problems

This section analyzes the relationship between the maximum equivalent stress in the concrete liner and some affected parameters in the context of the deterministic problem. The concerned parameters can be listed as the thickness of two liners, the deconfinement rates in the CCM, the limitation of deformation of CMC as well as the mechanical properties of the host rock.

The deterministic results are important to understand the drift responses before accounting for the uncertainty of the last parameters (i.e., the mechanical properties of Cox claystone). Our studies focus only on the maximum of the concrete liner's equivalent stresses, hereafter called the interested stresses. For the sake of clarity, we mention here that the equivalent stress present in fact the von Mises stress (the stress related to the distortion energy, as well as the deviatoric strain), which is expressed in the general 3D case as:

        2 2 2 2 2 2 1 6 2                          eq xx yy yy zz zz xx xy yz xz (3.3)
which can also be written as a function of the principal stresses:

      2 2 2 1 2 2 3 3 1 1 2                  eq (3.4)
In the 2D case, it can be reduced to:

2 2 2 3           eq xx xx yy yy xy (3.5)
To illustrate the influence of the mentioned aspects on the stresses, we need to fix some parameters to consider the aspect interest. In other words, except for the variation of the parameter in parametric-study of interest, other parameters of our numerical applications, the thicknesses (l1) and mechanical properties of the tri-linear elastic compressible material (Ec1, Ec2, Ec3, Figure 3-13 above shows the displacements around the tunnel in a specific study case along with three types of the axis being interested. Convergence is about 8.3cm was noted on the tunnel wall. This displacement induces a high volumic deformation of about 0.41, which is uniform in the compressible liner.

1 v  , 2 v  1  ), as in
l1 (m) Ec1 (MPa) Ec2 (MPa) Ec3 (MPa) 1 v  2 v  1  0.2 100

Influence of the thickness of the liners

As the first numerical investigations, we study the time-dependent evolutions of the maximum equivalent stress on the concrete liner with respect to the thickness of two liners. In Figure 3-14, the maximum equivalent stresses quickly increase at the very beginning of the period (i.e., in a few days) and become stable beyond the duration of about ten years. On the contrary, in these cases, the compressible layer's volumetric deformation depends much more on the thickness l1 than l2 (see Figure 3-16). They do not reach the third elastic stage (i.e., the volumetric deformations are still lower than the limit 2  v ). Its compressible potential does not benefit when the volumetric deformations are very far from the limit 2  v .

(a) (b)

Figure 3-16: Volumetric deformations on 1 st liner versus time (a) Variation of 1 st thickness (fixed l2=0.5m), (b) Variation of 2 nd thickness (fixed l1=0.2m)

Influence of the deconfinement rate

The effect of the deconfinement rate variation on the maximum equivalent stresses in the concrete liner is also tested. In Figure 3-17a, the interested stresses slightly turn down from nearly 18.9MPa to about 18 MPa when the deconfinement rate changes in the range of 0.7 to 1.0. Similar to the previous case, the interested stresses become stable after ten years. The results captured in Figure 3-17 prove that the deconfinement rate does not affect the stress in the case that the mean values of COx rock properties are used as input data. Another case study is then conducted by adopting a random set of the input parameters (noted as set-1) with Aset1=8.415e-6; B set1=2.647; C set1=0.192; and Young Modulus: E set1=4.253GPa. These values correspond to about 1.7, 1.3, 1.1, 0.7 times with respect to their mean value. In these simulations, the limit 2 0.4 3-18 illustrates that the deconfinement rate can now effect much more pronounced on the maximum equivalent stress of the concrete liner. More details can be seen in Figure 3-18a; these stresses rapidly reduce from more than 63MPa to lower than 29 MPa when the rate exceeds 0.85

  v . Figure

 

. Although in all five cases of the considered deconfinement rate ([0.7, 0.775, 0.85, 0.925, 1.0]), the volumetric deformations in Figure 3-18c are not much different, the stresses in Figure 3-18b is really dissimilar, especially for two last cases compare to three first ones. This phenomenon (as the deformation in the first three cases) has reached the third branch with 2 0.4

    v v
. To deeply understand it more, let consider these effects as the next section (variation of the limitation 2  v ). The following discussions will focus on the impact of the compressibility (v2) of the compressible liner (CMC), i.e., the limit value of volumic deformation between the second and the third elastic stage (as in Figure 3-6b). Firstly, we use the mean values of the input parameters of COx rock properties. The effects are not obvious, except the cases 0.1   v even for the variation of the thickness l2 (as in Figure 3-19a). The relations of the volumetric-stresses and the deformations (in case of l2=0.5m) are shown in Figure 3-19b. In which no difference occurs with most values of 2  v (except a case study 1 2 0.015

    v v
).

(a) (b) Unlike the first case, in many events with other combinations, the compressibility of the CMC has an extremely significant effect on the distribution of stress state on the concrete liner. In other words, the compressive stress in this inner liner rapidly increases when the volumetric strain in the covered compressible liner v exceeds the limit value 2  v (i.e., the elastic behavior of this outer support element is in the 3-rd stage). Now, to strongly clarify the observed phenomenon, we take the second limitation 2  v varies from 0.015 to 0.515, with the spacing of unity, and the input parameters are taken from the set-1 (mentioned in 3.4.2). Figure 3-20a shows that the effect of 2  v changes significantly in the range from 0.315 to 0.415 (with all three cases of the thickness l2). For the typical case of concrete liner thickness (l2=0.5m), Figure 3-20b illustrates the difference between the four first lines (e.g., 2

 v ranges from 0.015 to 0.315) tendency is difficult to be captured from these results. Following that, with 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 Figure 3-22 a, b, c). This phenomenon seems to be consistent because the increase of one of these parameters yields a highly viscous deformation, according to the Lemaitre model expressed in Eq. (1.58). If it is evident for a constant applied deviatoric stress, it is not the present underground structure. In this last context, the interaction in a non-linear manner between the convergence of the viscousrock and the reaction from the support system generates a redistribution of stress state representing by a variation of the deviatoric stress at a specific point of rock mass in time. This redistribution depends strongly not only on the short and 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 the liners. Concerning the effect of the short-term mechanical properties of COx rock (i.e., Young's modulus), the results from Figure 3-22d present lower equivalent stress in the concrete element when Young's modulus is much higher than its mean value of 6.0 GPa. Maximum equivalent stress in the concrete-liner can especially reach a value of 76 MPa (Figure 3-22), which largely exceeds the compressive resistance of the C60/75. Note that, in the typical civil engineering application by the deterministic (or semi-probabilistic) approach, a resistance about 60% of the concrete mark is mostly accepted in the serviceability limit state design of the concrete elements as proposed in European Code (EC2). Below this limit, one considers that the risk of crack under compressive loading will not happen in the structure. Similarly, the coefficient B also tends to increase the compressive stress in the drift's final support element (i.e., the concrete liner). Notably, in the case of using the Set-1, a rapid change happens when B>2.5). The significant impact of coefficient B on the interested stress can be stated from two figures (Figure 3

Influence of the Young Modulus E

There is little debate that the Young Modulus (E) has a strange impact on interested stress. In Figure 3-29, the parameter has a disproportional effect on the stresses (i.e., E increases from 0.2 to 11.8, along with the variation of ,max  eq reducing from approximately 21.9 MPa to 17.8 MPa. However, Figure 3-30 presents a low maximum equivalent stress when E is very small; then this stress increases to above 65 MPa with the middle range of value E; the next, it suddenly falls to under 30 MPa for the highest range of the modulus (three last values). Although the effect of different parameters was highlighted in the deterministic problem, a clear trend of the COx rock properties, accounting also for the fact of their correlations (i.e., correlation of A, B, C), on the behavior of the final support of drift was not captured. In addition, the obtained results do not allow us to quantify the effect of these parameters on the stability at the long-term of this structural element. For this last aim, the reliability analysis based on the Kriging metamodeling technique will be carried out in the next section.

Discussions on reliability analysis problems

For the sake of clarity, we precise that the probability of exceedance is only measured in the concrete liner when the maximum equivalent stress exceeds the preselected threshold lim2.

This exceedance probability prediction is undertaken for five case-studies by using the modified AK-MCS method; this technique was detailed in Chapter 1 and Chapter 2. As demonstrated previously, the convergence of exceedance probability can be attained sooner (i.e., the number of iterations at convergence is lower) if a set of new training points is used in each iteration. Thus, hereafter, for all the numerical investigations of reliability analysis, a set of four enriched samples in each iteration during the iterative process is chosen thanks to the capability of using parallel calculations 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 of drift support, which relate to the convergence condition and to the interpolation by Monte Carlo Simulation, are chosen with =0.01, N =6, NMCS=10 6 (see their definitions in Sec. 2.5.1 in chapter 2).

Effects of the adopted threshold stresses

The considered concrete liner of drift as adopted in this work has the characteristic compressive strength according to the Eurocode-2 (European code for design of concrete structures) fck= 60MPa (i.e., the concrete class C60/75). Thus, the value of the design compressive strength (fcd) is defined as: / 0.9*60/1.5 36 (MPa)

     cd cc ck C f f (3.6)
where: cc=0.9 denotes the coefficient taking account of long term effects on the compressive strength and of unfavorable impacts resulting from the way the load is applied; C=1.5 is the partial factor for materials for ultimate limit states.

Influence of the deconfinement rate

A discussion of this section reveals how the rate of deconfinement could affect the stability of the final liner with the same conditions (lim2=36MPa, l2=0.5(m)). Figure 3-33a below shows the decrease of Pf and the increasing deconfinement rate, from the peak point at 0.48% (=0.7) to 0.02% (=1), i.e., 24 times the tolerance. Thus, the chosen deconfinement rate to install the liners can strongly affect the probability of the final drift support's exceedance.

Influence of the compressibility of the CMC

As the last numerical investigation, by keeping lim2=36MPa and l2=0.5(m), here, the compressible potential of the outer layer can induce a quite significant variation of the exceedance probability in the drift's final support (Figure 3-33b). The concrete liner's longterm stability is higher, representing a lower probability when the compressible potential in the first support layer is more critical. For a particular case of an elastic material (i.e., 2 1 0.015

v v    
) characterizing by an elastic modulus (Ec=100MPa), the exceedance probability of the concrete liner of 0.5(m) of thickness at 100 years can attain to 51% instead of 0.11% (the case 2 0.515 v

 

). This result confirms the evident effect and incredibly significant benefit of using the compressible material to ensure the stability of drift support. 3.5.5. Summary of the reliability analysis problems

These obtained results from the reliability analysis highlighted a small probability of exceedance of deep drift support by adopting a thickness of 0.2(m) in the compressible layer that covers final concrete support of 0.5m of thickness. The predicted failure probability is smaller than 1% when the allowable stress is higher than 26MPa. However, based on different simplified assumptions, as mentioned in the description section of the studied problem, the reliability analysis in this chapter presents the preliminary results on drift's long-term stability analysis. We can expect that the exceedance probability of the drift support can be higher than the ones obtained in this research when the aleatory uncertainty and other mechanisms (e.g., damage, THM-coupling) are counted.

Concluding Remarks

In this chapter, various numerical applications were conducted to analyze the long-term stability of the deep drift support by taking the epistemic uncertainty of the mechanical properties of the host rock (COx claystone) into account. The support system of drift consists of a compressible material whose mechanical behavior is described by a tri-linear elastic model, contact directly with the host rock, and a linear elastic concrete liner. The study was undertaken in the 2D plane strain problem in the isotropic far-field stresses to model the circular drift oriented in the major horizontal stress of the COx layer. The time-dependent mechanical behavior of this rock was assumed to be described by the classical viscoplastic model of Lemaitre. The epistemic uncertainty of the mechanical properties of COx host rock was quantified from data obtained from the experiments in the laboratory, which were performed on the samples from the same depth.

The numerical results of both the deterministic problem and the reliability analysis, which is based on the modified AK-MCS method, 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 measured at 100 years of this final support of drift decreases when its thickness increases and the compressible potential of the outer liner is higher. The numerical investigations confirm the tremendous benefit of compressible material on the long-term stability of drift support from a purely mechanical point of view. However, the results in this preliminary research can be underestimated when the aleatory uncertainty related to the spatial variability of COx rock properties, as well as the damage mechanism and/or the coupling effect, are neglected. In geotechnical engineering applications, mechanical properties are uncertain, and they naturally vary in space. Such inherent spatial variability in rocks properties can considerably affect the reliability of geotechnical structures (G. [START_REF] Fenton | Risk assessment in geotechnical engineering[END_REF][START_REF] Al-Bittar | Bearing capacity of strip footing on spatially random soils using sparse polynomial chaos expansion[END_REF]Dian Qing Li et al., 2016;[START_REF] Xiao | Full probabilistic design of slopes in spatially variable soils using simplified reliability analysis method[END_REF] and greatly affect the behavior of tunnels during and after construction as well as their long-term responses (Ki Il [START_REF] Song | Effects of spatially variable weathered rock properties on tunnel behavior[END_REF]. For instance, ignoring soil properties' spatial variability often overestimates the probability of exceedance (Pf) of a slope [START_REF] Griffiths | Probabilistic slope stability analysis by finite elements[END_REF]. It could be more dangerous with high coefficients of variation (COVs) of the random parameters or the low factor of safety (FS).

Uncertainty related to a random variable's spatial variability (aleatory uncertainty) is usually represented through the random field concept (see subsection 1.1.2). 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 the numerical modeling (i.e., the direct evaluation of the performance function), this random field needs to be discretized by an appropriate technique. After the discretization of all random fields, the probabilistic assessment can be undertaken as an uncertainty problem tackled in the previous chapters. Nevertheless, the significant increase in the number of random variables after the discretization of these random fields results in pronounced difficulty in the construction of the metamodel since the computational demand increases exponentially as the number of variables grows. This problem is also referred to as the curse of dimensionality.

Solving the great dimensional problem is still an active ongoing research topic in the literature, and almost all research efforts are based on the idea of using dimension reduction methods. An interesting tool termed high-dimensional model representation (HDMR) has been received a lot of attention in the last two decades [START_REF] Rabitz | General foundations of high-dimensional model representations[END_REF]. This method prescribes a systematic mapping procedure between the inputs and outputs to reveal the hierarchy of correlations between the input variables. In reality, only relatively low-order correlations of the input variables significantly affect each physical system's output. This property permits the expression of the HDMR in only a few hierarchical levels to represent the physical system accurately. Each hierarchical level of HDMR aims to apply an appropriate projection operator to the output function. The latter will be chosen in this study between the two well-known HDMRs: the ANOVA (Analysis of variance)-HDMR and Cut-HDMR.

The main aim of this chapter is, therefore, to investigate the effect of the spatial variability of COx rock properties on the stability at the long-term of the deep drift support by using the modified AK-MCS method as well as the applicability of the Cut-HDMR in combination with the Kriging-based reliability analysis (noted herein as AK-HDMR). This chapter is organized as follows. The initial state of the problem with considering spatial variability of the COx rock is firstly introduced. An appropriate discretization technique adopted in the case of crosscorrelated random fields is then detailed. The following section consists of presenting some results of the reliability analysis of deep drift by the modified AK-MCS by accounting for the spatial variability. The comparison of these results (i.e., variability problem) with the ones of the previous chapter (i.e., uncertainty problem) is undertaken, and we highlight in parallel the limit of the AK-MCS to tackle the problem of high dimension. After a brief presentation of the Cut-HDMR method, its combination with the AK-MCS reliability analysis technique to handle COx rock properties' variability is detailed. After several numerical applications, we discuss the applicability of this AK-HMDR method.

4.2. Problem statement and procedure to handle the spatial variability 4.2.1. Set the initial state of the problem

In this chapter, as detailed in the previous chapter, almost the adopted hypothesis will be retained to describe the problem statement of deep drift constructed in the COx host rock. The deep-circular drift excavated in the viscoplastic rocks is supported by a concrete liner (elastic behavior) and a tri-linear compressible material layer. The two-dimensional (2D) plane strain model to represent the circular drift in the hydrostatic far-field stress is also maintained during this work. However, in comparison with the problem tackled in the previous chapter, the spatial variability of COx rock properties is accounted for, which means that the surrounding rock mass of the tunnel can be heterogeneous. The adopted 2D plane-strain model, counting also for the symmetric conditions, seem so to be a strong hypothesis to represent a spatial variability problem that needs to be really modeled in three-dimensional (3D). However, regarding the very high consuming time and computation demand to simulate a time-dependent behavior of deep tunnel in 3D and the fact that our main purpose is to validate the potential applicability of the chosen method (i.e., the modified AK-MCS and the AK-HDMR method as presented below), the adopted hypothesis may be acceptable to simplify the procedure.

Following that, the deep drift support's exceedance probability at the period of exploitation of 100 years will be measured considering the variability in the COx rock properties' space. The reliability analysis bases firstly on the modified AK-MCS metamodel, as mentioned in the last two chapters. Meanwhile, the spatial variability of four random parameters (coefficients A, B, C, and Young modulus E), which present the time-dependent behavior of the COx claystone, will be described by four random fields, and there are cross-correlation relationships between the first three coefficients. Each random field can be characterized by a proper correlation length. Nevertheless, for the methodology verification purpose, and since all the data concerning the correlation length of the mechanical properties of COx rock are available, we assume that all the random fields have the same correlation length.

Below, the well-known Expansion Optimal Linear Estimation (EOLE) method proposed by [START_REF] Li | Optimal discretization of random fields[END_REF] will be adapted to discretize the random fields and extended to tackle their correlation relationship (i.e., discretization of the cross-correlated random fields). This discretization procedure presents a necessary step for the direct evaluation (i.e., solving the deterministic problem numerically) by the open-source Code-Aster.

Discretization of random fields by the EOLE method

For civil structures, the increasing interest of the research community in the probability analysis concerning spatial variability (or space-variant properties) points out the problem of achieving a reliable discretization of random fields. A discretization procedure is required to decrease a continuous random field to a finite set of random variables. The precision of discretization is a crucial point for any subsequent probabilistic investigation of structural response. It is useful to optimally formulate a discretization error estimator to quantify the discretization's precision and to require that the approximate random field fill a prescribed target precision. In which the number of optimization variables is small, however, it could represent the random field, and the objective function must depend on the discretization error estimator and the target precision [START_REF] Allaix | Discretization of 2D random fields: A genetic algorithm approach[END_REF].

In practice, it is not straightforward to sample a Gaussian random field H(x, ) to obtain the realization of H. Consequently, the discretization procedure is appeared to approximate the random field H(x, ) by ˆ( , )



x H

. The methods of discretization can be divided into three groups [START_REF] Sudret | Stochastic Finite Element Methods and Reliability: A State-of-the-Art Report[END_REF], namely: point discretization [START_REF] Der-Kiureghian | The stochastic finite element method in structural reliability[END_REF], average discretization [START_REF] Fenton | Simulation of random fields via local average subdivision[END_REF][START_REF] Griffiths | Probabilistic slope stability analysis by finite elements[END_REF], and series expansion methods [START_REF] Phoon | Implementation of Karhunen-Loeve expansion for simulation using a wavelet-Galerkin scheme[END_REF][START_REF] Sudret | Comparison of finite element reliability methods[END_REF][START_REF] Cho | Probabilistic assessment of slope stability that considers the spatial variability of soil properties[END_REF][START_REF] Jiang | Slope reliability analysis considering spatially variable shear strength parameters using a non-intrusive stochastic finite element method[END_REF]. In which, the series expansion methods, commonly used in the literature, can be listed as Karhunen-Loève expansion (KLE), orthogonal series expansion (OSE), and the expansion optimal linear estimation method (EOLE). The crucial point here is that the random target functions can be adequately simulated by serial expansion methods using a finite number of deterministic functions and random variables-coefficients. These methods are highly appreciated for more complex problems. In our problems, the latter is suitable and chosen among series expansion one, thanks to its accurate, more practical, and efficient than some other series expansion discretization techniques according to [START_REF] Zhang | Time-Dependent Reliability Analysis Through Response Surface Method[END_REF].

The EOLE method can be used to approximate the Gaussian stochastic input parameters. On the basis of the pointwise regression of the original random field wrt the selected values of the field, and compaction of the data by spectral analysis [START_REF] Li | Optimal discretization of random fields[END_REF]; [START_REF] Allaix | Numerical discretization of stationary random processes[END_REF]. The subsequent problem takes into account the derivation of an EOLE illustration of a random field with two aims: (a) limit the computational effort devoted to the application of the resulting representation; and (b) ensuring that the representation meets all accuracy standards. By limiting the highest order of the expansion term and the total number of random variables used to construct a representation, the first objective can be achieved. In contrast, the second goal will be attained by finding the minimum of multiple error estimators.

The EOLE method has the advantage of allowing the error variance of the corresponding discretization scheme to be determined. Thus, one can determine the optimal number of eigenmodes for a specified value of error variance. Notice that the discretization of a random field by this method conducts an expression that gives the random field's value at each point in the space of the rock mass as a function of N random variables (following the standard normal distribution law). This value N equals the number of eigenmodes. For a given value of the error variance on the EOLE, the number N is small for high values of autocorrelation distances. In which, homogeneity is a special case according to the autocorrelation distances gets infinite value. However, N can be an extremely high value when the autocorrelation distances are so small [START_REF] Haj | Probabilistic analysis of an offshore monopile foundation taking into account the soil spatial variability[END_REF].

Mathematically, the EOLE method represents a stochastic field in terms of a linear combination of deterministic functions/vectors ( , )

  , ( ), ( ) ( ) ( ) ( , )           Σ x x x x x x j k j k j k j k Cov H H (4.2)
The EOLE of random variable H(x) onto the random vector χ reads:

1 ( ) ( ) ( ) ( ) ( )          χ χχ x x x Σ x Σ χ μ T H H H (4.3)
where:

( ) ( ) χ Σ x H x
is a vector whose components are given by ( ) ( ) ( ), ( ), ( )

            Σ x x
( ) 1 ( ) ( , ) ( , ) ( )             x x x N j T j H x j j H H   (4.7)
By defining the variance of H(x) is 2 ( )  x , the error variance for EOLE after basic algebra is:

  2 2 ( ) 1 1 [ ( ) ( )] ( )         x x x N T j H x j j
Var H H   (4.8)

In Eq. (4.8), the second term is identical to the variance of ˆ( ) x H . Thus, EOLE always underestimates the exact variance. However, the error decreases monotonically with N, which helps one automatically define the cut-off value of N for a provided tolerance in the variance error ().

Discretization of the cross-correlated random fields

In engineering practice, the parameters of random fields and their cross-correlation, which are often unknown, must be estimated from extensive measurements. Thus, it is a high challenge to either precisely approximate random field parameters or correctly simulates cross-correlated random field samples (RFSs) since the number of measurements is sparse and limited (due to sensor failure, budget limit, etc.) [START_REF] Zhao | Simulation of cross-correlated random field samples from sparse measurements using Bayesian compressive sensing[END_REF].

In our study, the adopted 2D normally distributed random fields are discretized by crosscorrelated in the context of the EOLE method. In which parameter E is independent of others. Thus, it can be directly applied to Eq (4.7). The three coefficients A, B, C being in the Lemaitre model, as mentioned in chapter 3, have cross-correlations between each pair. The crosscorrelation structure between each pair of simulated fields is simply defined by crosscorrelation coefficients. The method requires all cross-correlated fields on the domain to share an identical autocorrelation function. A cross-correlation coefficient easily represents the crosscorrelation structure between each pair of simulated fields. Now, by considering the case of three cross-correlated random fields (the coefficients A, B, C), the cross-correlation matrix reads:

1 0 0 ˆ1 0 0 1 0 0                                                                          C C C C C C C C C C C C C C C C C C C C C C T AB AC AB BC AC BC (4.9) i.e.,                  C C C C Φ Λ Φ T ( 4
                    C C C C C C C C C I I I Φ I I I I I I N N N D N N N N N N (4.11) and,   1 2 3     C C Λ I I I D C N N N diag (4.12)
Then, the cross-correlation block sample matrix is defined as

    1/ 2 3 1 3 3 3 3 ( )                 χ Λ D D D j N N N N N  (4.
( ) 0 0 0 0 ( ) ( ) ( )                                                                                      A C C C C C C C C B C C C C C C C C C C C C C C C C χ N N D j N D A j D B N N N N N N N N j D C j I I I I I I I I       1 3 1 31 1 32 2 33 3 1                                        A C B C C C C C C C N N N N N N N I I I I (4.14) 0
Here, each field is created by using a set of independent random variables, and these sets are then correlated wrt the assumed cross-correlation matrix between three expanded random fields according to the framework presented by Vořechovský [START_REF] Vořechovský | Simulation of simply cross correlated random fields by series expansion methods[END_REF]. Thus, the EOLE representation Gaussian random field in Eq. (4.7) can be rewritten as follows:

, ( ) 1 , ( ) 1 , ( ) 1 ( ) ( , ) ( , ) ( ) ( , ) ( , ) ( ) ( , ) ( , )                                                   x x x x x x x x x D N A j T A A A A j H j j D N B j T B B B B j H j j D N C j T C C C C j H j j H H H H H H       (4.15)
where ,   j j are eigenvalue and eigenvector of each Gaussian auto-correlation matrix (as in the case of non-correlated random fields). Notice that the length of the vector T j  and

( ) x H
 are equivalent to the number Mgrid of the grid points (e.g., if the random fields are distributed with all quadrangle shapes, Mgrid =Mx*My, with Mx and My are the total point numbers in the horizontal (x-axis) and vertical one (y-axis), respectively). In our study, Mgrid is the total number of all cell-centroids.

Application to the COx rock properties

The values of the truncated order of expansion (N) and the correlation length are two main factors for the assessment and control the accuracy of the discretization methods. These factors significantly affect the point-wise estimator for variance error of the discretization. Indeed, the EOLE method always under-represents the true variance of the random field. Thus, the accuracy of this method is strongly correlated with both factors.

We assume that the coefficients A, B, C, and Young modulus E are for Gaussian random fields that share the same autocorrelation length function. Table 4-16 below shows the increasing expansion order (N) corresponding to the reduction of correlation lengths in both cases (isotropic x = y and anisotropic spatial variability x  y, where x and y are correlation lengths in x-and y-axis respectively). In each case, we also consider two options for the limitation of the error of 5% or 20%, which means the maximum variance error can be accepted in the discretization. The current concept here is that whether the higher variance error (i.e., 20%) in some cases is possible to use to replace the lower one without the loss in accuracy estimating of the drift's stability.

It is worth to note that, for the practical simulation, to reduce the EOLE point-wise error variance at the boundaries, we extend the random field mesh with a small value, e.g., 3.0 m (i.e., the boundary of the random filed domain [-3.0, 58] m, more detail in Appendix A3).

Besides, each element size must be sufficiently refined (i.e., LRF/ ≤ 1/6 where LRF and  are the typical element length of random field mesh and the correlation length, respectively [START_REF] Sudret | Stochastic Finite Element Methods and Reliability: A State-of-the-Art Report[END_REF]). Therefore, in the table above, the same conditions, such as the grid, the domain of random fields, the partition size, or the size of cell elements for the discretization are applied.

For the parametric study, the two following chosen cases are considered. More precisely, with respect to a variance error, 5% the isotropic variability x=y=15m (noted as case study 1) and anisotropic variability with x=68m, y=5.05m (case study 2) have the same number N=25. Figures 4.1 and 4.2 are presented, respectively, the results of the eigenvalues and the pointwise estimator for variance errors of these study cases.

In fact, these two study cases (1 and 2), which have an identical value of expansion order term N=25 totally generate 100 random variables in the reliability problem. This large size of random variables corresponding to the theoretically limited value that the AK-MCS method can handle, as discussed in many studies in the literature, which performed, however, in some academic problems. 4.3. Primary numerical results

Results of the deterministic problem

In this part, we present some representative results of the deterministic problem that were conducted with respect to one realization of the discretized random fields of the case studies 1 and 2. to its chosen thickness. Then in Figure 4-8, the effect of the compressibility of the compressible material on the maximum stress in the final support of drift is investigated. The results of both cases 1 and 2, as well as the ones of the uncertainty problem extracted from the previous chapter (Figure 4-8c), are highlighted for comparison purposes. From these results, the spatial variability of the host rock reduces the maximum equivalent stress in the concrete liner. Nevertheless, it is important to clarify here that, in cases 1 and 2, the obtained results are only taken from one realization, and they are not sufficiently representative from a statistical point of view to confirm the observed tendency. The probabilistic results from the following reliability analysis seem much more appropriate for this task.

(a) (b)

Results of the modified AK-MCS reliability analysis

The modified AK-MCS method is firstly chosen to estimate the probability of exceedance of the concrete liner of deep drift in the context of spatial variability of COx host rock. The two cases 1 and 2 are the primary interests because of their significantly high number of random variables (100 variables in each case). These case studies allow us to verify the possibility of AK-MCS to handle the more complex problem, which, theoretically, can account for up to 100 variables, according to some scholars (e.g., see [START_REF] Lelièvre | Development of AK-based method for reliability analyses[END_REF][START_REF] Wang | An adaptive high-dimensional model representation method for reliability analysis of geotechnical engineering problems[END_REF]). Contrary to these last contributions, which limit the academic problems, our numerical study shows that the AK-MCS can only be reasonably applied when the maximum number of variables is about 50. Beyond this number, we state an increase drastic of the training points in the DoE and the time-consuming to construct the Kriging surrogate.

Figures 4-9 and 4-10 below show the exceedance probability in the concrete liner that is predicted as a function of the allowable stress for the two cases 1 and 2. However, in these studies, only the spatial variability of the coefficient A and Young's modulus E are considered to limit the total number of random variables at 50 (w.r.t. the truncated order term N=25). More precisely, by regarding the formula (4.14 and 4.15), only the random variables   4.30), it seems that the spatial variability induces a higher exceedance probability in the concrete support of drift when its allowable stress is chosen at a value inferior to about 30MPa. However, with the limit value lim2=36MPa, in these cases results, a much smaller probability in comparison with the uncertainty problem. But the observations must be verified in the case of more random fields (i.e., the random fields of the two parameters B and C are taken into consideration in these two cases). The effect of all four random fields was then investigated but limited in the case of quite high correlation length representing by a moderate value of the order term (i.e., N value). Concretely, by keeping the correlation length on X-axis at 50m and the thickness of the concrete liner at 50cm, different values of correlation length on Y-axis (e.g., y equals 50m, 45m, 35m, 25m, and 15m) are undertaken.

Table 4-17 and Figure 4-11 below are summarized the exceedance probability at 100 years in the concrete liner by considering the variation of correlation length in the Y-axis. They illustrate that decreasing correlation length reduces the exceedance probability in the final support element of drift. Specifically, in the case of isotropic correlation length x = y = 50m, the exceedance probability value (Pf=0.098%) is slightly smaller than the relevant case of uncertainty problem (Pf=0.11%) evaluated in the previous chapter. By using the anisotropic correlation with a smaller correlation length y, the probability decreases significantly to about Pf=0.02% with y = 15m. Moreover, the results in Table 4-16 below also show a significant increase in the number of call, Ncall, compared to the uncertainty problem. These studies highlight the spatial variability effect of the COx rock properties on the exceedance probability of the concrete liner of drift. In comparison with the uncertainty problem, a smaller probability can be stated in this variability problem. However, this observation must be strengthened in the intensively numerical investigations taking into account different values of correlation length of COx rock properties (e.g., smaller correlation length), and of the concrete liner thickness as well as with a more appropriate geometrical model of drift (e.g., the total 2D plane strain model without symmetric conditions or in the general 3D problems). Nevertheless, the limit of the AK-MCS to handle the curse of dimensionality does not permit to conduct this kind of study, which needs, in fact, an appropriate extension of the AK-MCS metamodeling technique.

HDMR method

The high-dimensional model representation (HDMR) method has been largely considered in recent years to handle the problem of the curse of dimensionality. The principal idea of this method is to approximate a multivariate performance function by a combination of zeroth-, first-, second-, and higher-order component metamodels [START_REF] Li | High dimensional model representations[END_REF]. In other words, HDMR is developed as a set of quantitative models to gain the input-output relationships of large dimensional physical systems with numerous input variables.

HDMR was originally introduced by [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF] for the purpose of conducting sensitivity analysis and was further investigated [START_REF] Sobol | Theorems and examples on high dimensional model representation[END_REF]. HDMR consists in decomposing a performance function G(X) in terms of increasing dimensionality yielding an expression of the form:

1 2 1 2 1 2 1 2 1 2 1 2 0 ... 1 1 1 ... 12... 1 2 ( ) ( ) ( , ) ... ( , ,..., ) ( , ,..., ) l l l m i i i i i i i i i i i i i i i m i i i m m m G g g x g x x g x x x g x x x                   x (4.16)
where g0 indicates a constant term denoting the 0-th order effect. The first-order component function gi(xi) is a univariate function that represents an individual contribution to the output g(x). The second-order one gij(xi, xj) is a bivariate function describing the interactive effect of the input parameters xi and xj upon the output response.

Some experiences in the literature [START_REF] Rabitz | General foundations of high-dimensional model representations[END_REF][START_REF] Sobol | Theorems and examples on high dimensional model representation[END_REF] shows that the form above (Eq. 4.16) can be truncated up to two orders because the high-order interactions among input variable of HDMR are negligible. The function is now written as an accurate approximation of g(x) with respect to the variable vector   1 2 , ,..., ,...,  x

T i m x x x x , namely, 1 2 1 2 1 2 0 1 1 ( ) ( ) ( ) ( , )           x x  m i i i i i i i i i m G G g g x g x x (4.17)
By determining all the summands on the equation above, the HDMR can be used as a computationally efficient meta-model for predicting the response. Practically, two techniques available in the literature to identify the components of decomposition: cut-HDMR and ANOVA-HDMR (also referred to as Random sampling RS-HDMR). The cut-HDMR could perform much better than the other metamodeling methods, especially when handling the function with weak parameter interaction (Chen et al., 2019).

The cut-HDMR is considered in this study. Once a cut-HDMR method is applied, an anchor point (or a reference point)

  1 2 , ,...,  u T m u u u
needs to be defined. The reference point could be selected as either the means of all random variables (e.g., [START_REF] Wang | An adaptive high-dimensional model representation method for reliability analysis of geotechnical engineering problems[END_REF]) or the point among the initial DoE being the closest point to the Limit State Function (LSF) [START_REF] Lelièvre | Development of AK-based method for reliability analyses[END_REF]. For example, the point provides the minimum absolute value of G(u). The component functions of HDMR read as

( )  u o g G (4.18) ( ) ( , )   u i i i i o g x G x g (4.19) 1 2 1 2 1 2 1 2 1 1 2 2 , ( , ) ( , , ) ( ) ( )     u i i i i i i i i i i i i o g x x G x x g x g x g (4.20)
where g0 is the constant value of the performance function measured at the anchor point u is the 0-th order component function. Equations (4.19) and (4.20) define the first-, and secondorder component function, respectively. u denotes the vector being including the remaining element of vector u after element i-th of u is replaced by xi (or xi1, xi2). Specifically, (xi, u i ) = (u1, ..., ui-1, xi, ui+1, …, um) and (xi1, xi2, u i1,i2 ) = (u1, ..., ui1-1, xi1, ui1+1, …, ui2-1, xi2, ui2+1,…, um).

The Cut-HDMR has been proven to be efficient because the functions of higher-order components are usually negligible or have fewer effects on the output of a system and can, therefore, be ignored [START_REF] Wang | An adaptive high-dimensional model representation method for reliability analysis of geotechnical engineering problems[END_REF]. Although the first order HDMR model (noted as HDMR1) is obvious, it may not be accurate for some problems. High-order HDMR models, such as the second-order HDMR model (noted as HMDR2) written in Eq. (4.17), can be applied [START_REF] Chowdhury | High-dimensional model representation for structural reliability analysis[END_REF].

AK-HDMR technique

Recently, different methods can be used to derive each component of HDMR (i.e., surrogate), such as polynomial regressions, Kriging, support vector machines, or RBFs. In this work, the AK-MCS Kriging metamodel, as introduced in the previous chapters, will be chosen to construct each surrogate of the Cut-HDMR, and we note the adopted method as AK-HDMR1 and AK-HDMR2 corresponding to each configuration of the truncated HDMR method at the first or second order.

In the AK-HDMR1 method, a combination of HDMR1 with the Kriging-based reliability analysis, each Kriging metamodels ( , ) u i i G x of the AK-HDMR1 is independently constructed. The final approximate performance functions, results of the summands of these surrogates can be used to interpolate the result of the performance function at each random point x. More precisely, this interpolation leads to the following expression of Kriging prediction and variance [START_REF] Lelièvre | Development of AK-based method for reliability analyses[END_REF]:

    0 1 ( ) (1 ) ( , )        x u   i m k k i G G i m g x (4.21)     2 2 1 ( ) ( , )      x u   i m k k i G G i x (4.22)
with m the dimension of the vector of random variables.

Correspondingly, in the AK-HDMR2 method, the Kriging metamodels ( , )

u i i G x and 1 2 1, 2 ( , , ) u i i i i G x x
of the AK-HDMR2 are independently constructed, which leads to the following expression of Kriging prediction and variance [START_REF] Lelièvre | Development of AK-based method for reliability analyses[END_REF])

        1 2 , 1 2 1 2 1 1 0 1 1 1 2 ( 3) ( ) (2 ) ( , ) ( , , ) 2                  x u u    i i i m m m k k k k i i i G G G i i i i m m g m x x x (4.23)         1 2 , 1 2 1 2 1 1 2 2 2 1 1 1 ( ) (2 ) ( , ) ( , , )               x u u    i i i m m m k k k k i i i G G G i i i i m x x x (4.24)
From these mean (Eq. 4.23) and variance (Eq. 4.24) values of the Kriging prediction, the classical criteria U (i.e., Eq. (1.39)) can be rewritten as follows:

      ( ) ( ) ( )    x x x   k k G k G U (4.25)

Construction procedure of AK-HDMR

The procedure of AK-HDMR can be described as 13 steps as follows (also in Figure 4-12)

(1) Generation Monte-Carlo population X of NMCS samples (points x (k) )

(2) Definition of the initial DOE Random selection of N1 points in X to evaluate on G Drawing of N1 points from the population of Monte-Carlo above. (N1 = 20 seems a good compromise between numerical efforts and efficiency of the method.)

(3) Selection of an anchor point (u) for HDMR decomposition -Assessment of the performance function at these points -Selection of the point corresponding to the lowest absolute value of the performance function. This point is the so-called "anchor point".

(4) Definition of the first DOE by LHS of N2 points -Generation by Latin Hypercube Sampling of a set of points in standard space.

-These points are then projected onto the axes passing through the anchor point.

-Evaluation of the performance function on these projected points

(5) Calibration of Kriging np=do metamodels for HDMR1 or  

1 2 0 0   p d d n C C metamodels for HDMR2. (6) Prediction , G i  , , G i
 of X for each metamodel (7) Evaluation of the response Ĝ of HDMR decomposition using calibrated models.

(8) Estimation of exceedance probability (e.g., by MCS with the estimation of the probability of exceedance by Eq. (1.23), i.e., ˆ( ) 1 1 ( )

   1  p G n k f k p p x n ; where ˆ( ) ( )  1 G k
x is equal to 1 if the kriging prediction is negative for the point u (k) , and equal to 0 otherwise.

(9) Evaluation of the enrichment function U (Eq. 4.25). This requires an identification of the point of the population of Monte-Carlo being enriched. The learning function U is used with the prediction and the variance defined above. The point to enrich is therefore 4.17 below taken from Table 2-9), it can be stated that the AK-HDMR1 cannot be accurately evaluated the exceedance probability of the deep underground structure. More precisely, in comparison with the MCS, the results of AK-HDMR1 and AK-HDMR2 present a difference of about 16.5% and 0.53%, respectively (Table 4.17 below). Although some scholars demonstrated the efficiency of the HDMR1 method (e.g., [START_REF] Lelièvre | Development of AK-based method for reliability analyses[END_REF][START_REF] Wang | An adaptive high-dimensional model representation method for reliability analysis of geotechnical engineering problems[END_REF] by investigating several academic problems, the present numerical investigation shows the inaccuracy of such a method to study in the underground structure in which the interaction of random variables is important. The AK-HDMR2 method provides an appropriate result compared to the modified AK-MCS and MCS methods. However, because the number of metamodels that need to be constructed in this method is essential, the number of calls (Ncall) to the deterministic problem increases significantly regarding modifying one's AK-MCS method.

G G g g G g g G g g G g G G g G g G g g G (
(a) (b) This section proposes two approaches to reduce the size of a metamodel for the drift in the COx rock problem. Since the previous study confirmed that the first-order Cut-HDMR (i.e., AK-HDMR1) is not sufficiently accurate to tackle the problem in which the interaction of the random variable is high. In such a situation, the second-order Cut-HDMR (i.e., AK-HDMR2) is required to improve the accuracy of the obtained results. However, the classical AK-HDMR2 needs a huge number of metamodels: total (1+m+m(m-1)/2) metamodels to be calibrated (i.e., the number of functions need to be identified in Eq. (4.17)) for the vector of m random variables. Therefore, once the AK-HDMR2 is adopted, the number of Kriging metamodels becomes exceptionally large when the size of the problem increases (e.g., the number of random fields, random variables, or DOFs). For example, if m equals 100, one must build 5051 metamodels, which is unaffordable.

Regarding the case of deep drift, the discretization of the four random fields of the COx rock properties (i.e., coefficients A, B, C, and Young's modulus E) induces a remarkably high dimension of vector of the random variables. From the equations (4.14 and 4.15), this random variable vector can be written

        1 1 1 1 , , ,              x A B C E N N N N with dimension 4   m N (i.e.
, four times of the truncated order term of the expansion, N). The previous discussion does not allow us to apply AK-HDMR2 with respect to each variable of each corresponding vector

  1 , A, B, C, E    i N i
. Instead, we propose to construct the metamodels with respect to each

variable component vector   1   i N (i=A, B, C, E) of the global vector x.
Thereupon, the function of AK-HDMR2 in Eq 4.17 can be rewritten as follow:

              0 1 1 3 1 1 3 1 1 3 1 1 3 1 1 1 2 1 1 1 2 1 1 1 2 1 1 1 ( ) ( ) , , , , , , , , , , , , 
                                                    x x u u u u u u u u  A B C E A B A C A E B C A B C E N N N N N N N N A B A C A E N N N N N N N N N B C N N G G g g g g g g g g g       2 1 1 1 2 1 1 1 2 , , , ,                  u u B E C E B E C E N N N N N N N g g (4.30)
Once the spatial variability of the deep drift by using the Cut-HMDR2 method is applied, ten surrogates (w.r.t. ten terms in Eq. 4.30, except g0) is enough to approximate the response. is about the isotropic spatial variability problem with the correlation lengths x = y = 50m, the term order N=6. In these studies, we consider the thickness of concrete liner l2=0.25(m) while the threshold stress is fixed at lim2=36MPa. The comparison of the two methods (e.g., modified AK-MCS and AK-HDMR2 methods) exhibits good accordance, which confirms the accuracy of the AK-HDMR2 method to treat the uncertainty problem (Figures 4-14a). However, in the variability problem , the gap between the two methods is quite large. Notably, the initial number of calls (Ncall) to evaluate the deterministic problem (i.e., the evaluation of initial DoE) in the AK-HDMR2 method is much more important regarding the necessarily constructed surrogates (10 surrogates). Although the more significant number of DoE evaluations, the time-consuming to construct the surrogates of the AK-HDMR2 is significantly decreased, particularly in the variability problem, thanks to the reduced number of random variables being handled in each metamodel. In this chapter, the long-term stability of the underground structure is considered in the more general context by accounting for the spatial variability effect of the host rock properties. The random field concept is adopted to present the variability in the space of each mechanical parameter of the host rock. An extension of the well-known EOLE method is presented to discretize the cross-correlated random fields. The modified AK-MCS method seems appropriate for tackling the variability problem due to the high correlation length represented by a low or moderate number of random variables after the discretization. The highdimensional model representation (HDMR) approach truncated at the first order (Cut-HDMR1) or at the second-order (Cut-HDMR2) is chosen to handle high dimensional problems, known as the curse of dimensionality. In this method, the AK-MCS is combined to construct each component (i.e., surrogate) of the Cut-HDMR. The numerical investigations show that the exceedance probability predicted by the AK-HDMR1 is not accurate, and the AK-HDMR2 dramatically improves the result. Although the number of built surrogates increases in the AK-HDMR2, this method significantly reduces the number of random variables and hence the necessary time to construct each constituted surrogate. Once the accuracy of this AK-HDMR2 is demonstrated in the uncertainty problem, its applicability in the variability problem requires more improvements.

CONCLUSIONS AND PERSPECTIVES

Evaluation of the impact of uncertainties and spatial variability of the time-behavior properties of the host rock owing a high time-dependent behavior on the long-term stability of the constructed underground was the principal objective of this dissertation. Because of the studied period of hundred of years and the high time-cost of a single numerical simulation of such complex problems, an assessment of probability should be performed following a carefully chosen methodology, able to quantify the propagation of uncertainty in the exceedance probability of structure with a reasonable computational demand accurately and time-cost. This work's contribution relates mainly to the selection, adaptation, extension, and validation in different contexts of a chosen probabilistic method based on the metamodeling technique.

This dissertation started with a review of the literature, highlighted the background and the state of the art of the reliability analysis application on the construction of underground structures in rock masses. A revision of essential concepts of reliability analysis (uncertainty, spatial variability, random field, etc.) and the two primary sources of uncertainty, namely aleatory and epistemic, was presented. Aleatory uncertainty relates to inherent variability in the space of the formation properties, while the uncertainties of knowledge with respect to the parameter characterization and/or adopted model belong to the epistemic uncertainty.

Regarding the methodology of reliability analysis, it showed that the direct sampling Monte Carlo Simulation (MCS) had been largely chosen thanks to its simplicity and its capacity to provide an accurate estimate of the probability. The results obtained from this method were usually considered as the benchmark to validate many other probabilistic methods. However, the tremendously required number of evaluations of structure response in the MCS method presents its drawback to being applicable in the case of rock formations owing a complex behavior. The developments of some advanced probabilistic approaches in the last decades allowed us to overcome this disadvantage. In general, these approaches aim at approximating the implicit limit-state function (LSF) by a mathematical function (e.g., Kriging, Radial Basis Functions (RBF), Artificial Neural Network (ANN), Support Vector Machine (SVM), Polynomial chaos,…), the so-called metamodel (or surrogate), through which the prediction of the exceedance probability can be conducted. The main challenge of such surrogate-based reliability analysis lies in the way to build an adequate metamodel to approximate the LSF without wastes of samples and function calls. From this bibliographic survey, it showed that the Kriging metamodeling technique could provide high accuracy, notably for nonlinear problems. In addition, its flexibility in interpolating the sample points allows combining this technique with the classical sampling method, such as the MCS, to measure the exceedance probability. A well-known version of such a combination is the AK-MCS method, with some eventually recent extension/improvement that was then summarized. Surprisingly, although the accuracy and efficiency of the AK-MCS method have been intensively demonstrated in many structural design projects, it is rarely applied in geotechnics, particularly in the rock engineering field. The bibliographic study also revealed that almost all studies on the stability of deep underground structure dealt with the impact of uncertainty and/or variability of host rock properties on the underground's short-term stability.

The adaptation, validation, and efficient investigation of the AK-MCS metamodel were firstly undertaken in the context of a deep tunnel excavated in a rock with a viscoelastic behavior of Burgers. Based on the hypothesis that the deep tunnel is circular, supported by a system of two elastic liners and constructed in a hydrostatic condition stress state at far field, a closed-form solution was derived using an integral equation. In this analytical solution, the sequential tunnel excavation and the installation in a sequence of two liners can be considered. This analytical solution presents an efficient tool to investigate the influence of different factors (e.g., excavation rate, the thickness of liners, mechanical properties of host rock) on the tunnel's response and stability. Especially, the reliability analysis based on the direct MCS method highlighted the critical effect of the many parameters, notably the uncertainty of the Burgers rock properties on the exceedance probability at the long-term of two failure modes: by an excessive convergence of tunnel beyond an acceptable value or by stress beyond the allowable long-term strength of concrete liner. These results provided by the MCS method were then used as the benchmark to validate the chosen Kriging-based reliability analysis. Indeed, the Kriging metamodel can be established iteratively through an enrichment process of the Design of Experiment (DoE). In order to find out the best candidate points to enrich the DoE, which should be not only close to the limit-state but also far away from the training samples of existing DoE, we proposed a modification of the well-known AK-MCS method by adding a constraint distance on the selection of the new training point of DoE. The validation and significant improvement of this modified AK-MCS were demonstrated in an academic problem and the case of the deep tunnel in the Burgers rock. Comparing the results provided by the MCS method and the modified AK-MCS metamodel revealed a drastic decrease in the required number of structure response evaluations in this last method. Our investigations also highlighted the advantage of using a subset of new training points to enrich the DoE that could reduce the number of iterations to attain the convergence of the exceedance probability. This fact is useful for the study in which the numerical evaluations of the performance function can be conducted through parallel calculations.

The efficiency of the modified AK-MCS was then examined in the more complex problem, which aims at investigating the long-term stability of deep drift support in the context of nuclear waste disposal in the Callovo Oxfordian (COx) claystone. The significant interest in such a context lies in the stability analysis of the gallery's concrete support during the period of exploitation for about 100 years to ensure the functionality of the mechanical system within the tunnel. Especially, a compressible outer liner is considered in the support system of drift. It is a highly porous material that allows this compressible outer layer to absorb the host rock's convergence over-time and reduce the radial stress transmitted to the inner liner. Due to the fact that the study focuses mainly on the effect of the uncertainty of COx rock properties and the contribution of the compressible liner on the long-term stability of the concrete inner lining, the considered problem is limited only to the purely mechanical behavior of the circular drift.

In addition, although the developed procedure was performed in the general context that allows us to study drifts constructed in different directions with respect to the principal stress state, only the case of drift excavated following the major horizontal stress was taken into consideration. Indeed, to reduce the overstress on the support element and to limit the inevitable fracture zone induced by excavation, drift parallel to the major horizontal stress (so a configuration in which cross-section stress is isotropic), is a preferred choice for the long-term safety of nuclear waste management. Furthers, in order to characterize the long-term behavior of COx host rock, the viscoplastic model of Lemaitre is chosen. This model's capacity to catch the essentials of time-dependent behavior of COx rock with a reasonable number of parameters to be calibrated significantly facilitates the uncertainty quantification process. The assessment of uncertainties of parameters of this model is carried out using data from the laboratory tests performed on the samples from the same depth, allows to estimate the epistemic uncertainty of the mechanical properties of COx host rock. To limit the problem on the impact of uncertainties of COx, the behavior of the compressible material described by a tri-linear elastic model being considered as deterministic. The numerical investigations of both the deterministic problem (in the 2D plane strain condition using the open-source Code_Aster) and the reliability analysis, based on the modified AK-MCS method, 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 measured at 100 years of this final support of drift decreases when the thickness of the compressible layer increases or when its compressible potential is higher. These numerical results also confirm the tremendous benefit of compressible material on drift support's long-term stability from a purely mechanical point of view. The result also shows that the exceedance probability of the inner liner is negligible for the actual structural design of the considered underground drift. This demonstrates the robustness of the design approach for the Cigéo project.

The consideration of the spatial variability of COx rock properties on the stability of drift support was carried out in the last part of this dissertation. The uncertainty of each mechanical parameter of the host rock is described by random continuous spatial fields, written as correlation functions with respect to the spatial correlation length. The expansion optimal linear estimation method (EOLE) was applied for discretizing the random fields to handle the continuous uncertainty problems. An adaptation of this EOLE method to take into account the correlation of three random fields corresponding to three parameters of the Lemaitre viscoplastic model of COx rock was presented. The probabilistic assessment by using the modified AK-MCS was then undertaken with respect to the discretized random fields in which the similar correlation lengths of all mechanical parameters of COx rock were assumed. Two configurations, isotropic and anisotropic space variability with a total of 100 random variables generated from the discretization step, were chosen to initiate numerical applications. These two case studies allow us to verify the limits of AK-MCS at handling a significant high dimensional problem. Our research shows that AK-MCS can only be reasonably applied for a maximum number of variables around 50. Beyond that, the number of training points in the DoE and the construction of the Kriging surrogate leads to a drastic increase in calculation time. According to abundant literature in this active ongoing research topic, an alternative approach is then explored, consisting of the use of the dimension reduction methods. Among these methods, the applicability of the high-dimensional model representation (HDMR) approach truncated at the first order (Cut-HDMR1) or at the second-order (Cut-HDMR2) was investigated. In this chosen method, the AK-MCS was used to construct each component (i.e., surrogate) of the Cut-HDMR. The numerical investigations showed that the Cut-HDMR1 is not sufficient, and the Cut-HDMR2 was necessary to improve the obtained results. The accuracy of the AK-HDMR2 was tested in the uncertainty problem of both cases of deep tunnels excavated in the viscoelastic and viscoplastic rocks. However, more improvements to this AK-HDMR2 are required to handle the variability problem.

This work, being the first tentative on establishing a methodology for the assessment of longterm stability of deep underground construction, including the uncertainties and the variability of properties of rock host, open many interesting perspectives for future research. Some of them are indicated hereafter.

Firstly, in order to overpass the limits of the time-consumption for a single numerical simulation, it could be worthy of looking for a quick deterministic response. The closed-form solution of the deep tunnel obtained for viscoelastic rock could be extended to consider more realistic conditions such as the anisotropy of the far-field stress, the shape of the cross-section, the anisotropy of mechanical properties of the host rock, and/or the coupled hydro-mechanical, thermo-hydro-mechanical behavior. In fact, a fully analytic solution considering all these elements could be challenging work, but the recent developments in the literature by using the complex variable/potential approaches in the Laplace domain demonstrated its feasibility using a semi-analytic approach.

On this same idea of considering more realistic behavior for rock masses, the progressive, damaging /healing process could be another direction of ameliorating long-term analysis of underground structure, presenting another challenge for the reliability analysis problem.

Secondly, from the methodological point of view, it is necessary to better consider in the methodology presented here a sensitivity analysis. Indeed whilst the efficiency of the Kriging metamodeling technique as compared to the sampling MCS technique was demonstrated, a sensibility analysis seems necessary. By maintaining only the most sensible parameters/factors, this step allows avoiding the curse of dimensionality in the reliability analysis due to a significant increase of the random variables when different aspects (coupled phenomena, damage mechanism, the anisotropic effect of geometry and/or loading) are handled. This sensibility analysis is crucial to reduce the computational demand, particularly in the case that the deterministic problem is only accessible by the numerical simulation.

Thirdly, it is necessary to deeply ameliorate the uncertainty assessment of the behavior of rock masses. In this work, this assessment is performed using laboratory data and some observations in situ, which were taken into account during the transfer of some properties from laboratory ones to the infield. The ad-hoc choice of the constitutive model leads in itself to high uncertainties. It is important for the reliability analysis of drift support in COx claystone, the developments/improvements of previous propositions (i.e., the consideration of the coupled problem, damage mechanism, and the anisotropic effect) to deeply investigate different types of drifts and incorporate all huge data from the in-situ measurements, in order to assess in one hand, as large as possible the behavior of rock masses and, in the other hand, an objective criterion for evaluation of exceedance probability. Indeed, due to the interaction of these complex mechanisms, the uncertainty quantification of COx rock properties by exploiting the measurements in an underground laboratory can only be operated through the numerical inverse procedure. Consequently, it requires in this task a powerful inverse analysis technique (such as the Bayesian or Artificial Neuron Network methods) regarding the considerable timeconsuming of each simulation. Through this inverse analysis procedure, we can expect the possibility of characterizing the spatial distribution (i.e., the estimation of spatial correlations) of the COx rock properties, which remains an unknown value.

Fourthly, the improvement/development of the Kriging metamodel in combination with the HDMR or the other methods to handle the high-dimensional problem is extremely important to deeply investigate the spatial variability effect on the stability of the deep underground structure. In this context, some other well-known techniques to reduce the dimension of the problem, such as the Principal Component Analysis (PCA) or Partial Least Square (PLS), can be contemplated to combine with the HDMR method.

Lastly, further future works can be consist of dealing with spatial-temporal variability. The developments/extensions of the metamodeling technique with respect to time variable, as proposed in some recent contributions in the literature, present an interesting tool to tackle more complex problems, which incorporate, for instance, the unsaturated flow, chemical interaction, and degradation mechanism of liners. Moreover, the spatial variability, limited until now in 2D conditions for the sake of simplicity, must be revisited in the real 3D context. This problem seems only to be overpassed with the help of a high performance computing cluster. This developed procedure of reliability analysis can be extended in the long term for the real-time risk analysis and optimization design of deep underground structure thanks to the capability of uncertainty quantification of the continuous data measurement. This stochastic procedure will also provide a useful and performant tool for designing and implementing the observation/surveillance in situ and/or in the laboratory.

Table A1-1. Rock mass classification and guidelines for supports of 10m span rock tunnels with RMR system [START_REF] Bieniawski | Engineering rock mass classifications: a complete manual for engineers and geologists in mining, civil, and petroleum engineering[END_REF] Table A1-2. Correlation equations between RMR and Q-systems [START_REF] Zhang | Multivariate adaptive regression splines for analysis of geotechnical engineering systems[END_REF] In brief, empirical methods can be widely applied in rock engineering due to the discontinuous, heterogeneous, and anisotropic nature of the rock mass. They are simple to apply and usually provide safe designs of excavation problems, which have been validated by numerous case studies in many countries. Nevertheless, the support selection relies on a finite number of parameters which do not ensure to characterize the rock mass properties entirely.

A.1.1.2. Analytical solutions Due to the limited applications of empirical methods beyond the preliminary design phase, especially the lack of considering the mode of failure and excavation geometry, analytical approaches could be used to get a more exact understanding of the ground response. In underground excavation engineering, the induced stresses caused by the digging and displacements in the massif are obvious the most interesting aspects. Thus, there various theory and solutions have been proposed and investigated for them. These are based on the universal rules of mechanics as well as on strong hypotheses concerning the geometry of the problem (circular tunnel with the horizontal axis at great depth), stratigraphy (a single homogeneous layer), the law of behavior (Tresca or Mohr-Coulomb) and the initial stress state (non-heavy, isotropic and homogeneous soil). Based on the assumptions above, and back to a solution proposed by Kirsch [START_REF] Kirsch | Die theorie der elastizitat und die bedurfnisse der festigkeitslehre[END_REF] and [START_REF] Poulos | Elastic solutions for soil and rock mechanics[END_REF], the tangential/radial stresses (r) and the radial displacements (Ur) of an unlined tunnel are defined in a polar coordinate system, for the case of linear elasticity, as follows: Similarly, for the lined tunnel, the relevant equations for these stresses (r) and displacements (Ur) are found by [START_REF] Timoshenko | Theory of Elasticity[END_REF] as follows:

2 , 2 2 0 1 1 ( / ) i i r i r r r r              (A1.3) 2 2 0 1
1 2 1 ( / )

s i i r s s i r U r E r r r                        (A1.4)
where Es and s are Young's modulus and Poisson's ratio of the liner respectively, r0 is the interior radius of the liner, ri stands for the exterior radius, and i is the magnitude of the radial compression applied to the liner-rock contact [START_REF] Bobet | Characteristic curves for deep circular tunnels in poroplastic rock[END_REF].

For more complex problem, such convergence-confinement method (CMM) (Panet & Guenot, 1982;Panet, 1993;[START_REF] Panet | Le calcul des tunnels par la méthode convergence-confinement[END_REF]Carranza-Torres & Fairhurst, 2000;Unlu & Gercek, 2003), which correctly simulate the final surface settlement and widely used for preliminary assessment of tunnel deformation and support acceptability for simple cases [START_REF] Langford | Application of Reliability Methods To the Design of Underground Structures[END_REF].

This method requires a priori estimation of the stress release coefficient  (or the rate of deconfinement or the confinement loss), which represents the relative position of the considered tunnel section compared to the tunnel face. Selecting the rate of deconfinement corresponding to the convergence occurring before the support starts interacting with the ground is the most vital point in the CMM [START_REF] Panet | The convergence-confinement method AFTES-recommendations des Groupes de Travait[END_REF]. It requires to have a series of pressure-displacement curves, namely: the longitudinal displacement profile (LDP), the support reaction curve (SRC), and the ground reaction curve (GRC). The LDP shows the radial displacement that occurs along the axis of an unlined excavation both ahead of and behind the face [START_REF] Langford | Application of Reliability Methods To the Design of Underground Structures[END_REF]. The solutions for LDP calculation are synthesized by [START_REF] Paraskevopoulou | Analysis of time-dependent deformation in tunnels using the Convergence-Confinement Method[END_REF] as in Table 1-4. Other two curves can be figured out with several assumptions to simplify the analysis, namely: hydrostatic stress field (0); circular and deep tunnels (boundary conditions of the problem to infinity), continuous, homogeneous and isotropic rock mass behavior; bidimensional problem and plan stress field [START_REF] Oreste | The convergence-confinement method: roles and limits in modern geomechanical tunnel design[END_REF].

Based on the analytical solutions above, for an elastic model of tunnel excavation, the released deformation at the tunnel face is suggested to be 25% of the total deformation [START_REF] Panet | Le calcul des tunnels par la méthode convergence-confinement[END_REF]. Specifically, considering the rate of deconfinement, for unlined tunnels, equations (2-9) and (2-10) could be modified as follow: For the lined tunnel, one can refer to Panet's paper [START_REF] Panet | Le calcul des tunnels par la méthode convergence-confinement[END_REF] to obtain further details. Besides, [START_REF] Sagaseta | Evaluation of surface movements above tunnels. A new approach[END_REF] and [START_REF] Verruijt | A complex variable solution for a deforming circular tunnel in an elastic half-plane[END_REF] proposed formulations for surface settlements and horizontal displacements within the framework of an elastic medium. [START_REF] Lee | Subsidence owing to tunnelling. I. Estimating the gap parameter[END_REF], by combining the elastoplastic deformations around the front with the effects of pre-convergence and the vacuum of the ground-support contract, has defined an annular vacuum parameter called GAP. Based on the GAP parameter, surface settlement can be calculated through empirical relationships [START_REF] Janin | Tunnels en milieu urbain: Prévisions des tassements avec prise en compte des effets des pré-soutènements (renforcement du front de taille et voûteparapluie[END_REF]. Also, from this parameter, one can predict the resulting ground deformations using 2D finite element methods [START_REF] Lee | Subsidence owing to tunnelling. I. Estimating the gap parameter[END_REF]. This notion of the annular vacuum was then taken up by Loganathan [START_REF] Loganathan | Analytical prediction for tunneling-induced ground movements in clays[END_REF]. Besides, for circular excavation problems, Mohr-Coulomb (M-C) failure criterion subjected to the hydrostatic in situ stress was given by Duncan [START_REF] Fama | Numerical modelling of yield zones in weak rocks[END_REF]. Hoek-Brown failure criterion can be solved using the formulae by Carranza-Torres (2004). For circular liners with particular reference to composite supports in both elastic and elastoplastic porous media can be found in Carranza-Torres (2009a, 2009b) and Bobet and Einstein (2011).

Briefly, there are several advantages and disadvantages to analytical methods. The analytical methods make it possible to quantitatively determine the design parameters from a model intended to schematize the behavior of the structure under the effect of the stresses applied to it. By analytical solutions, one can quickly obtain orders of magnitude of the calculation parameters and to estimate the influence of certain parameters on the response of the whole supporting ground. However, their direct field of application is so limited due to many restrictive calculation assumptions on which they are based. In this thesis, analytical methods refer to closed-form solutions for the viscoelastic problem are illustrated in Chapter 3.

A.1.1.3. Numerical method Numerical methods at present are widely used to simulate underground excavations due to the complicated geometries, support sequences, irregular shapes, geological conditions. They have the advantage of directly taking into account the response of the terrain and the support, as well as a large part of the specific features of the project (both geometry and geotechnics). Numerical methods in modeling of geomaterials include finite element method (FEM), finite difference method (FDM), finite volume method (FVM), boundary element method (BEM), and discrete element method (DEM). All of them can be used for treating engineering problems numerically. Thus, no matter what type of them is used, the constitutive model, the relevant hypothesis, and simplifying assumptions must be included in the code. In this thesis, we consider only the finite element method because of its convenient and quite common available methods to interact with the users. Thus, hereafter, the numerical method in our work also means the finite element method. It is worth pointing that, while the time consuming is an essential issue in numerical modeling, it needs to base on appropriate constitutive behavior (i.e., Tresca, Mohr-Coulomb, Hoek-Brown, etc.) to determine the resultant support loads [START_REF] Langford | Application of Reliability Methods To the Design of Underground Structures[END_REF]. That is also a key reason for a new proposed model that is carried out in Chapter 4 to fasten the computation in the numerical method.

A.1.1.4. Observational method

The observational method in rock and underground engineering, which dates back to the "learn as you go" or "design as you go", was first formally introduced by Peck [START_REF] Peck | Advantages and limitations of the observational method in applied soil mechanics[END_REF] for minimizing risks. This method is briefly mentioned in Eurocode 7, the new European geotechnical design code: "When prediction of geotechnical behavior is difficult, it can be appropriate to apply the approach known as "the observational method", in which the design is reviewed during construction" [START_REF] Frank | Designers' guide to EN 1997-1 Eurocode 7: Geotechnical design-General rules[END_REF]. In other words, this code states that the observational method is a viable alternative to conventional design when the geotechnical engineer is confronted with knowledge-based uncertainties of ground behavior [START_REF] Nossan | Observations on the Observational Method[END_REF].

Unlike the three mentioned methods, which predict field behavior in advance, the observational one deals with the actual behavior such as stress, strain, or displacement, so forth. If the acceptable limit during construction is reached by the modernized measuring instrument, one must stop the construction and trigger the pre-defined contingency plan by implying that there is no or only a very basic design available prior to the starting of construction [START_REF] Palmstrom | Ground behaviour and rock engineering tools for underground excavations[END_REF]. It can be highly occurred as subproblems of larger problems, through as the actual behavior of the rock masses with excavations frequently differs from the prediction of numerical or analytical calculations. Moreover, observations on the behavior served as an additional input for the choice of the construction method and determination of lining type and thickness [START_REF] Schubert | The Development of the Observational Method[END_REF].

A.1.1.5. Conclusions of deterministic analysis methods Four mentioned approaches (empirical, analytical, numerical, and observational method) can be used while they still lack a standardized method for underground support design. The approach for a deterministic problem is successful whenever the model analyses effectively and considers both the problem's geological and geometric characteristics [START_REF] Langford | Application of Reliability Methods To the Design of Underground Structures[END_REF].

In geotechnical engineering, deterministic methods give only a partial representation of the true margin of safety [START_REF] Goh | Reliability assessment of stability of underground rock caverns[END_REF] and often lead to an extremely conservative design. Thus, it is unable to account for uncertainties in material and load properties and cannot reflect how safe a geotechnical structure really is. Non-deterministic methods provide a means to explicitly treat the uncertainties, and it will be discussed in the next section.

Before addressing non-deterministic methods, it is useful to revisit the more global vision of the mentioned traditional approaches in the process of rock engineering as the Figure A1-1, proposed by [START_REF] Stille | Classification as a tool in rock engineering[END_REF]. A.1.2. Appendix for reliability analysis methods and stochastic stability analysis AK + PCEs AK-Polynomial Chaos Expansions [START_REF] Schöbi | Combining Polynomial Chaos Expansions and Kriging[END_REF] Table A1-5: Two-dimensional autocorrelation functions [START_REF] Castillo | Probabilistic analysis of structures using stochastic finite elements[END_REF])

                                                                               h t t B B B K K K t t K M K K M K h h t t K K K K K K K K e p
  1 1 2 1 0 1 0 1 1 1 2 0 1 2 1 0 1 0 2 0 1 ( ) ( ) ( ) exp exp ( ) ( ) ( )exp ( ) ( )exp exp ( ) ,                                                               h h t B K K K M K K K K h h t t
    0 0 2 2 2 1 1 0 0 ( ) ( ) exp 1 exp exp t t l K K ini r K ini r K m v G G A R d m R v d R v                                                 1 1 0 0 2 2 2 2 1 ( ) ( )exp 1 exp exp , t t l K K fin t t K fin K m v G G A R d m R d R                                             (A2- 10)                                                         t t B B B B K K t t K M K K t B B K t K M K e e G G
                    K K M K e e G t
whereas the corresponding parameters 1 2 1 2 , , , E E   are respectively equal to with: Therefore, 2 ( ) B f t can be expressed in the form:

                                                                 h h t B B K K K M K K K K h t t K K K t K K K K K e p
                               h t t
    2 1 2 0 2 0 2 1 1 1 2 2 1 2 2 0 2 1 2 3 1 2 1 2 0 2 2 3 ( ) ( ) ( ) exp exp ( ) .
( ) exp exp exp ( ) exp In which the solution of   22 p t can be determined using the Eq. (A2-15) as shown in Eq.(2.27)

                                                                         h h B K K M K K K h t K fin t K K K h t K K K t K M K
For simplifying the presentation, the expressions of the supporting pressures, as well as the displacements in the rock mass, are written in integral form, but it is worth noting that their explicitly analytical expressions can be found without difficulty. The numerical integration is only needed for the functions containing the term ( ) t

 determined in the range 0 t t  (i.e., during the excavation stage) such as ones expressed in Eqs. (2.30) and (2.31)) due to the presence of time in both the numerator and denominator of the exponential function (see also Eq. A2-10).

higher than er, and the partition value is selected at 48, which has low amounts of both the mean and the maximum variance error.

(a) (b) Ce travail de thèse vise à évaluer la probabilité de défaillance à long terme du revêtement en béton des tunnels profonds en considérant l'incertitude des paramètres du comportement différé de la roche hôte. Pour cela, une extension/modification de la technique d'analyse fiabiliste par le métamodèle de Krigeage a été d'abord mise en place. La performance et l'efficacité de ce métamodèle modifié en comparant avec la méthode classique MCS a été démontrée dans le cadre d'une construction séquentielle d'un tunnel profond dans une roche viscoélastique. Ce métamodèle de Krigeage est ensuite appliqué pour analyser la stabilité pendant la période d'exploitation de 100 ans du revêtement en béton d'une galerie construite dans la couche argileuse de Callovo-Oxfordien (COx), la formation géologique choisie pour le stockage profond des déchets radioactifs en France. La quantification des incertitudes et des corrélations des paramètres du modèle de Lemaitre pour la roche hôte a été réalisée en utilisant les résultats des essais de fluage au laboratoire. Le soutènement de la galerie étudiée est compressible et son comportement décrit par un modèle élastique tri-linéaire. Les résultats montrent l'impact des incertitudes des propriétés de la roche hôte et un grand avantage du voussoir compressible sur la stabilité à long terme du revêtement en béton. Selon ces études, la stabilité à 100 ans du revêtement est considérablement augmentée en fonction de la compressibilité et/ou de l'épaisseur de la couche compressible du système de support de la galerie. Une considération de la variabilité spatiale des propriétés mécaniques de la roche hôte a aussi été effectuée. Pour traiter le problème de grandes dimensions, associé à la discrétisation des champs aléatoires, une adaptation de la méthode Cut-HDMR combinée avec la metamodélisation de Krigeage a été proposée dans cette étude et a été ensuite testée à travers des applications numériques.

Mots-clés: incertitude, variabilité spatiale, analyse fiabiliste, métamodèle de Krigeage, tunnel profond, matériau compressible.

LONG-TERM STABILITY EVALUATION OF UNDERGROUND CONSTRUCTIONS BY CONSIDERING UNCERTAINTIES AND VARIABILITY OF ROCK MASSES

Summary: This dissertation aims to evaluate the exceedance probability at the long-term of the concrete lining of deep tunnels by considering the uncertainty of the time-dependent behavior of the host rock. To this purpose, an extension/modification of the reliability analysis technique using the Kriging-based metamodel was carried out. The performance and the efficiency of this metamodel compared to the classical Monte-Carlo Simulation (MCS) method were demonstrated in the context of a sequential construction of a deep tunnel in a linear viscoelastic rock. The metamodel is then applied to analyze the stability of the concrete lining of a gallery during a 100-years period of exploitation that is built in the layer of Callovo-Oxfordian (COx) claystone -the geological formation chosen for the deep nuclear waste disposal in France (Andra). By analyzing the raw database from creep tests performed at Andra, we not only found a significant uncertainty and crosscorrelation of the Lemaitre model's parameters but also developed an adopted viscoplastic model to characterize the creep behavior of COx rock. The reliability analysis exhibited that the long-term stability of concrete inner liner has been affected by the host rock properties' uncertainty and the great benefit of the compressible layer. According to these studies, the stability can be determined as a function of the compressibility and/or the thickness of the compressible layer. Moreover, the spatial variability of the mechanical properties of the host rock was also undertaken. To deal with the high dimensional problem associated with the discretization of the random fields, an adaptation of the Cut-HDMR method combined with the Kriging-based metamodeling technique was proposed in this study. The applicability of this method was then tested through numerical applications.

Keywords: uncertainty, spatial variability, reliability analysis, Kriging metamodel, deep tunnel, compressible material.

Figure 1 - 1 :

 11 Figure 1-1: Categories of uncertainty in soil properties ..........................................................

Figure 1 - 4 :

 14 Figure 1-4: Design point (or MPP) in (a) X-space and (b) U-space ........................................

Figure 1 - 5 :

 15 Figure 1-5: Transformation to Standard Normal distribution (Schulze, 2006) .......................

Figure 1 - 6 :

 16 Figure 1-6: Two approaches to creating a sample of size n = 10 ............................................

Figure 1 - 7 :

 17 Figure 1-7: The principle of the Kriging approach (Hawchar, 2017) ......................................

Figure 1 - 8 :

 18 Figure 1-8: Constant stress-strain-time coordinates (ASTM Standard D2990-09, 2011) .......

Figure 1 - 9 :

 19 Figure 1-9: Strain versus time behavior for four viscoelastic models (Dowling, 2012)..........

Figure 2 - 1 :

 21 Figure 2-1: The sequential excavation of the double-lined tunnel and the installation of two linings .......................................................................................................................................

Figure 2 - 2 :

 22 Figure 2-2: Influence of the excavation rate on the deconfinement rate .................................

Figure 2 - 3 :

 23 Figure 2-3: Validation of our solution to (Nomikos et al., 2011) ............................................

Figure 2 - 4 :

 24 Figure 2-4: Comparison of the solution to Wang's ones of (Wang et al., 2014) .....................

Figure 2 - 5 :

 25 Figure 2-5: Problem statement with numerical solutions ........................................................

Figure 2 - 6 :

 26 Figure 2-6: Validation of our solution to numerical solutions with radial displacements .......

Figure 2 - 7 :

 27 Figure 2-7: Validation of our solution to numerical solutions with stresses versus distances

Figure 2 - 9 :

 29 Figure 2-9: Sensitive analysis for the radial displacement and stresses on the liners (2) ........

Figure 2 -

 2 Figure 2-10: Sensitive analysis for the radial displacement and stresses on the liners (3) ......

Figure 2 -

 2 Figure 2-11: Time-dependent exceedance probability of the tunnel convergence ..................

Figure 2 -

 2 Figure 2-12: Exceedance probability of the tunnel convergence at 100 years ........................

Figure 2 -

 2 Figure 2-13: Exceedance probability at 100 years by considering vl, t1, t2 ..............................

Figure 2 -

 2 Figure 2-15: Exceedance probability (on the second liner) ..................................................... 80

Figure 4

 4 Figure 4-2: Point-wise estimator for variance errors at all central cells of the mesh ............

Figure 4 - 3 :

 43 Figure 4-3: Spatial distribution of the mechanical properties of COx rock of case study 1 ..

Figure 4 - 4 :

 44 Figure 4-4: Spatial distribution of the mechanical properties of COx rock of case study 2 ..

Figure 4

 4 Figure 4-5: Iso-values of equivalent stress in the concrete liner of a realization in variability problems: case study 1 (a), case study 2 (b) .........................................................................

Figure 4 - 6 :

 46 Figure 4-6: Influence of the concrete liner thickness on the maximum equivalent stress .....

Figure 4 - 7 :

 47 Figure 4-7: Influence of the concrete liner thickness on the maximum equivalent stress .....

Figure 4 - 8 :

 48 Figure 4-8: Influence of the compressibility of the compressible liner and inner liner thickness ................................................................................................................................................

Figure 4 - 9 :

 49 Figure 4-9: Probability of exceedance at 100 years in the concrete liner (case study 1): ......

Figure 4 -

 4 Figure 4-10: Probability of exceedance at 100 years in the concrete liner (case study 2):....

Figure 4 -

 4 Figure 4-11: Probability of exceedance at 100 years in the concrete liner versus Y (fixed X=50m) .................................................................................................................................

Figure 4 -

 4 Figure 4-12: Flowchart of AK-HDMR ..................................................................................

Figure 4 -

 4 Figure 4-13: Probability of exceedance in the concrete liner in the viscoelastic rock ..........

Figure 4 -

 4 Figure 4-14: Comparison AK-MCS to AK-HDMR2 in exceedance probability in the concrete liner ........................................................................................................................................

Figure 1

 1 Figure 1-1: Categories of uncertainty in soil properties

  of the plastic zone at equilibrium w.r.t tunnel radius R0

  7) where f(x) is the joint probability density function (PDF) of the relevant demands and resistance parameters x=[x1, x2,…, xd] T . The function of performance denotes g(x)=g(x1, x2, .., xd), and the LSS defines as g(x)=0. It can be illustrated with the case of two normally distributed parameters, as in Figure 1-2.

Figure 1 - 2 :

 12 Figure 1-2: Probability of failure for a system consisting of (a) one and (b) two variables

Figure 1 - 3 :

 13 Figure 1-3: Categories of structural reliability analysis methods

Figure 1 - 4 :

 14 Figure 1-4: Design point (or MPP) in (a) X-space and (b) U-space

Figure 1 - 5 :

 15 Figure 1-5: Transformation to Standard Normal distribution[START_REF] Schulze | Reliability as a Public[END_REF]

  ) Evaluation the function g(U) at sampling points by a deterministic model; (2) Determination vector a by solving the set of linear equations related to sampling points; (3) Using FORM to get a tentative design point and index ; (4) Repeat three first steps until the index converges; (5) Computing probability of failure by using FORM and/or SORM.

Figure 1

 1 Figure 1-6: Two approaches to creating a sample of size n = 10 (a) MCS scheme, (b) LHS scheme

Figure 1 - 7 :

 17 Figure 1-7: The principle of the Kriging approach[START_REF] Hawchar | Développement de méthodes fiabilistes dépendant du temps pour l'analyse de durabilité des structures: application au problème de conception fiabiliste dépendant du temps[END_REF]

  are the variance and autocorrelation function with specified correlation parameter , respectively; d stands for the size of design variables, , ,  component of xi, xj and , respectively. The correlation function (R) contains the assumptions about the approximation function, and hence the choice of an appropriate correlation function is an essential task of Kriging metamodeling. The following squared exponential autocorrelation (Gaussian) function of the hyperparameter vector θ is used in our work:

  stationary character of the Gaussian process, for each dimension, the correlation function depends only on the relative positions of two inputs characterized by a scale parameter i (i = 1, …, d). It is first required to obtain the optimal values of this hyperparameter  by the Maximum Likelihood Estimation technique as follows:

Figure 1 - 8 :

 18 Figure 1-8: Constant stress-strain-time coordinates (ASTM Standard D2990-09, 2011) 1.4.1. Constitutive models of linear viscoelastic rocks

  (a) Maxwell model, (b) Kelvin model, (c) Boltzmann model (Kelvin-Voigt), and (d) Burgers model. The strains versus time behavior for these models are illustrated as in Figure 1-9.

Figure 1 - 9 :

 19 Figure 1-9: Strain versus time behavior for four viscoelastic models[START_REF] Dowling | Mechanical behavior of materials: engineering methods for deformation, fracture, and fatigue[END_REF]

Figure 2

 2 Figure 2-1: The sequential excavation of the double-lined tunnel and the installation of two linings

  advancement curve over time can be calibrated from either measurement of tunnel convergence in the actual field or simulations with FEM. For instance, (Wang et al.FEM simulations. More specifically, in Figure2-2, the curves of the parameter  are shown with various excavation rates by using the mentioned values of m1 and m2. The figure illustrates the fact that the faster the excavation rate, the higher the deconfinement rate on the surface of the tunnel.

Figure 2

 2 Figure 2-2: Influence of the excavation rate on the deconfinement rateAs mentioned previously, the first-and second-liner are installed at two different instants

Figure 2

 2 Figure 2-3: Validation of our solution to (Nomikos et al., 2011) (a) Radial displacement, (b) Supporting pressure 2.3.3. Case of tunnels in the viscoelastic general Kelvin rock supported by double linersThe closed-form solution here could be considered as an extension of the prior contributions of(Wang et al. 2013) and[START_REF] Wang | An analytical approach for the sequential excavation of axisymmetric lined tunnels in viscoelastic rock[END_REF]) whose general solution was investigated for nliners tunnel excavated in the generalized viscoelastic Kelvin rock masses. Obviously, the Burgers material of this work behaves like a generalized Kelvin material, and our solution converges towards the solution of(Wang et al. 2013;[START_REF] Wang | An analytical approach for the sequential excavation of axisymmetric lined tunnels in viscoelastic rock[END_REF]) by using an infinite dashpot for the Maxwell element (i.e., M   ). Figure2-4shows good agreements between our study and the solution of[START_REF] Wang | An analytical approach for the sequential excavation of axisymmetric lined tunnels in viscoelastic rock[END_REF]) (the two liners are put at two different instants t1= 10(days) and t2=30(days) respectively with t0=1(day)).

Figure 2 - 4 :

 24 Figure 2-4: Comparison of the solution to Wang's ones of (Wang et al., 2014) On the cavity: (a) Radial displacement, (b) Radial stress, (c) Hoop stress Radial support pressure: (d) on the first liner, (e) on the second liners 2.3.4. Comparison to the numerical solutions This section consists of comparing our closed-form solution with numerical approaches attained from the finite element method (FEM) implemented on the Code_Aster (Code_Aster® software). Figure 2-5 illustrates the geometry and boundary and the grid of the modeling with one-quarter of the tunnel with a total dimension of 45.0(m) based on the symmetry problem.Notably, the initial conditions (namely, the boundary conditions and the sequential excavation and installation of liners) are kept as in the mentioned analytical simulation. That means on the top-and right-boundaries. We impose the hydrostatic stress 6.8(MPa) while on the bottom-and left-ones, the fixed normal displacement is applied. The quasi-instantaneous excavation is investigated with the initial time t0=5(minutes) while the two linings are put at t1=3(hours) and t2=2(days), respectively.

Figure 2 - 5 :

 25 Figure 2-5: Problem statement with numerical solutions (a) Geometry and boundary conditions (b) Mesh of the vicinity of the tunnel

Figure 2 - 6 :

 26 Figure 2-6: Validation of our solution to numerical solutions with radial displacements (a) versus time (on the cavity), (b) versus distance (in the surrounding rock mass of tunnel)

Figure 2 - 8 :

 28 Figure 2-8: Sensitive analysis for the radial displacement and stresses on the liners (1) (a) excavation rate vl, (b) coupled analysis: vl and t1, (c) coupled analysis: vl and t2 * The time considered in case studies of (b) and (c) is 100 years.

Figure 2

 2 Figure 2-9: Sensitive analysis for the radial displacement and stresses on the liners (2) (a) coupled analysis: t1 and t2, (b) coupled analysis for liner thickness: l1 and l2* The time considered in these case studies is 100 years.

Figure 2 -

 2 Figure 2-10: Sensitive analysis for the radial displacement and stresses on the liners (3) for viscoelastic input variables: (a) M, (b) K, (c) GM, (d) GK 2.4. Reliability analysis by the direct MCS

Figure 2 -

 2 Figure 2-11 and Figure 2-12 below show the exceedance probability respectively over time and at t=100 (years). Evidently, the monotonous behavior of tunnel deformation over time results in the same trend as the obtained curve of the probability of exceeding as expected.

Figure 2 -

 2 Figure 2-11 illustrates that the probability grows rapidly in the early years before reaching a lower rate. Depending on the speed of the tunnel excavation, the installation time, and the thickness of the linings, one can notice the strong variation in the exceedance probability at

Figure 2 -

 2 Figure 2-11: Time-dependent exceedance probability of the tunnel convergence (a) excavation rate vl, (b) installation time of liners t1, t2 (c) the first liner's thickness l1, (d) the second liner thickness l2

Figure 2 -

 2 Figure 2-13: Exceedance probability at 100 years by considering vl, t1, t2 (a, c, e, g) on the first liner, (b, d, f, h) on the second liner

Figure 2 -

 2 Figure 2-14: Exceedance probability (on the first liner) (a) Versus time by considering l1, (b) Versus time by considering l2, (c)(d): by considering (l1, l2) at 100 years with lim1=40MPa, 35MPa

Figure 2 -Figure 2

 22 Figure 2-16: A series four branches system: Evolution of the exceedance probability Table 2-8: The average results of example a series of four branches system Method Type Distribution of initial DoE

Figure 2

 2 Figure 2-19: The tunnel excavated problem: Comparison of different methods (a) Evolution of the exceedance probability, (b) Evolution of the constraint distance parameter D

Figure 3

 3 Figure 3-1: Possible architecture of the Industrial Centre for Geological Disposal (Cigéo).

Figure 3

 3 Figure 3-2: Meuse/Haute-Marne Underground Research Laboratory layout[START_REF] Delay | Synthesis of 20 years Research, Development and Demonstration in Andra's Underground Research Laboratory in Bure for Cigéo Project-France[END_REF]

Figure 3

 3 Figure 3-3: Convergence-confinement model w.r.t. the tunnel face progression[START_REF] El Matarawi | Investigating the use of reliability-calibrated partial factors in tunnel design[END_REF]

Figure 3

 3 Figure 3-4: Geometry of the studied disposal tunnel

  , the parameter 2 v  represents the compressible potential of the first support layer of drift to adsorb the high timedependent convergence of host rock. Moreover, ,  d v e are the volumetric and deviatoric part of the strain tensor, while the Bulk and shear moduli are calculated from the corresponding Young modulus and Poisson's ratio of each stage:

  Stress-strain curve of compressible material (a) under the compressive test (Billig et al., 2007), (b) our proposed tri-linear elastic model 3.3.3. Uncertainty quantification of the mechanical properties of COx rock

Figure 3

 3 Figure 3-7: Calibrated coefficients of Lemaitre's model from 88 samples with the histograms

Figure 3

 3 Figure 3-11: Surfaces and curves of viscous distortions versus time and deviatoric stress in the creep tests (a) with the approximation using the mean values (surface), (b) with boundary surfaces (1000 samples)

  numerical results show that the equivalent stress in the concrete liner is highest at the inner radius (i.e., at R=R3), as shown in Fig 3-12.a below. This expected result is consistent with different adopted assumptions, notably the hypothesis of isotropic hydrostatic stress at far-field and isotropic homogeneous material of components in the considered problem.

Figure 3 -Figure 3

 33 Figure 3-12: Equivalent stress in the concrete liner (a) and volumetric deformation in the outer compressible liner (b) at 100 years.

Figure 3 Figure 3

 33 Figure 3-14: Equivalent stresses on 2 nd liner versus time (a) Variation of 1 st thickness (fixed l2=0.5m), (b) Variation of 2 nd thickness (fixed l1=0.2m)

Figure 3

 3 Figure 3-17: Effects of the deconfinement rate with using Set-mean for four input variables (a)(b) to Equivalent stress on 2 nd liner without/with the variety of time, (c) to deformations of 1 st liner

Figure 3

 3 Figure 3-18: Effects of the deconfinement rate with using Set-1 for four input variables (a)(b) to Equivalent stress on 2 nd liner without/with the variety of time, (c) to deformations of 1 st liner 3.4.3. Influence of the compressibility of the compressible liner

Figure 3 -

 3 Figure 3-19: Effects of the limitation of the CMC deformation on interested stresses with using Set-mean (a) with a variation of the thickness l2 , (b) to the relation of strain-stress on CMC

Figure 3 Figure 3

 33 Figure 3-22: Maximum equivalent stress on concrete liner versus the parameter interest (a) Coefficient A, (b) Coefficient B, (c) Coefficient C, (d) Modulus Young E

Figure 3 Figure 3 Figure 3 Figure 3

 3333 Figure 3-25: Effects of the coefficient B with the mean values of A, C, E (Case-1) (a)(b) to Equivalent stress on 2 nd liner without/with the variety of time, (c) to deformations of 1 st liner

Figure 3 Figure 3

 33 Figure 3-29: Effects of the Young Modulus E with the mean values of A, B, C (Case-1) (a)(b) to Equivalent stress on 2 nd liner without/with the variety of time, (c) to deformations of 1 st liner

Figure 3

 3 Figure 3-32: Probability of exceedance at 100 years versus the thickness (a) 1 st liner (fixed l2=0.5m), (b) 2 nd liner (fixed l1=0.2m)

Figure 3

 3 Figure 3-33: Probability of exceedance versus the deconfinement rate (a), the second limitation of deformation of CMC (v2)

  for the random nature. Let us denote by χ the random vector {H(x1), ...H(xN)}. By construction, χ is a Gaussian vector whose mean value 

  of the eigenvalue (diagonal matrix) of the cross-correlation matrix     C Superscripts T denotes the transpose of the matrix or vector.It is denoted IN as a unity matrix of the order N (after ordering the vectors by decreasing eigenvalues). We can define the orthonormal eigenvectors of the correlation matrix ( Φ D

Figure 4 - 1 :Figure 4

 414 Figure 4-1: Eigenvalues w.r.t the index of the expansion order (N) (a) Isotropic cases x = y = 15m (case study 1), (b) anisotropic case: x = 68m, y =5.05m (case study 2)

Figure 4 Figure 4

 44 Figure 4-3: Spatial distribution of the mechanical properties of COx rock of case study 1 (a): coefficient A, (b) coefficient B, (c) coefficient C, (d) Young's modulus E.

Figure 4 -

 4 Figure 4-5 below illustrates the iso values of equivalent stress in the concrete liner of a realization in two variability problems (i.e., isotropic and anisotropic correlation lengths). The heterogeneous distribution of these parameters, as expected, present the consequence of the spatial variability of the COx rock properties. In Fig. 4-6 and Fig. 4-7, we capture the maximum equivalent stress in the concrete liner in time and at 100 years in the two cases 1 and 2. These results confirm the important dependence of the stress state in the concrete liner with respect

Figure 4 Figure 4 Figure 4 Figure 4

 4444 Figure 4-5: Iso-values of equivalent stress in the concrete liner of a realization in variability problems: case study 1 (a), case study 2 (b)

  while zero value is applied for all the elements of the other random variables  comparison with the results of the uncertainty problem tackled in the previous chapter (see Figure

Figure 4 Figure 4

 44 Figure 4-9: Probability of exceedance at 100 years in the concrete liner (case study 1): (a) versus threshold stresses of 2 nd liner, (b) versus iteration numbers (lim2=36MPa)

  ) x = y = 50m Variability (Anisotropic) x = 50m, y =45m Variability (Anisotropic) x = 50m, y =35m Variability (Anisotropic) x = 50m, y =25m Variability (Anisotropic) x = 50m,

Figure 4 -

 4 Figure 4-11: Probability of exceedance at 100 years in the concrete liner versus Y (fixed X=50m)

Figure

  Figure 4-12: Flowchart of AK-HDMR

  4.29)In which, both methods AK-HDMR1 and HDMR2 are developed from 20 initial samples for DoE as Step 2 of Figure4-12.

Figure 4 -

 4 Figure 4-13 below is captured the exceedance probability provided by the two methods AK-HDMR1 and AK-HDMR2. In comparison with the results obtained in the second chapter (the first four rows of Table4.17 below taken from Table2-9), it can be stated that the AK-HDMR1 cannot be accurately evaluated the exceedance probability of the deep underground structure. More precisely, in comparison with the MCS, the results of AK-HDMR1 and AK-HDMR2 present a difference of about 16.5% and 0.53%, respectively (Table4.17 below). Although some scholars demonstrated the efficiency of the HDMR1 method (e.g.,[START_REF] Lelièvre | Development of AK-based method for reliability analyses[END_REF][START_REF] Wang | An adaptive high-dimensional model representation method for reliability analysis of geotechnical engineering problems[END_REF] by investigating several academic problems, the present numerical investigation shows the inaccuracy of such a method to study in the underground structure in which the interaction of random variables is important. The AK-HDMR2 method provides an appropriate result compared to the modified AK-MCS and MCS methods. However, because the number of metamodels that need to be constructed in this method is essential, the number of calls (Ncall) to the deterministic problem increases significantly regarding modifying one's AK-MCS method.

Figure 4 -

 4 Figure 4-13: Probability of exceedance in the concrete liner in the viscoelastic rock (a) using HDMR1, (b) using HDMR2

where

  Ncall, HDMR1=N1 + N2,HDMR1 + Nparallel*Niteration,HDMR1 =20+20*4+3*18=154 (evaluations) Ncall, HDMR2=N1 + N2,HDMR2 + Nparallel*Niteration, HDMR2 =20+20*10+3*32=316 (evaluations). 4.4.4. Adaptation of the AK-HDMR2 to the studied problem of deep drift in the COx rock

4. 4

 4 .5. Applicability of the AK-HDMR2 in the case of deep drift in the viscoplastic COx claystoneWe present in Figures 4-14 the result of the exceedance probability in the concrete liner of drift by ignoring or considering the spatial variability of COx rock properties, which are evaluated from the modified AK-MCS and AK-HDMR2 methods. More precisely, the first figure (Figure4-14a) is related to the uncertainty problems (i.e., ignoring the spatial variability and the AK-MCS case was also mentioned in in Chapter 4). The right figure

Figure 4 -

 4 Figure 4-14: Comparison AK-MCS to AK-HDMR2 in exceedance probability in the concrete liner (a) uncertainty problem, (b) isotropic variability problem

  hydrostatic stress field (σ0 = σv = σh); E and  are Young's modulus and Poisson's ratio of the rock mass respectively; r0 stands for tunnel radius.

Figure A1- 1 :

 1 FigureA1-1: Relationships between ground behavior and rock engineering[START_REF] Stille | Classification as a tool in rock engineering[END_REF]

Figure

  Figure A1-2: State of the art of AK-based methods

Figure

  Figure A3-3: Contours of variance errors as functions of correlation lengths and expansion terms (a): N varies in [8, 20], (b) N varies in [20, 40] (both cases use the isotropic lengths)

  Tran, V. T. Mai, D. Hoxha, and M. N. Vu. (2019). Time-dependent reliability analysis of deep tunnel in the viscoelastic Burgers rock with the sequential installation of liners. Rock Mechanics and Rock Engineering 53. Springer: 1259-1285. https://doi.org/10.1007/s00603-019-01975-6.

Do D. P., M. N. Vu, N. T. Tran, and G. Armand. (2021). Closed-form solution and reliability analysis of deep tunnel supported by a concrete liner and a covered compressible layer within the viscoelastic Burger rock. Rock Mechanics and Rock Engineering (Minor revision). Tran N. T., D. P. Do, D. Hoxha, M. N. Vu, and G. Armand. (2021). Kriging-based reliability analysis of the long-term stability of a deep drift constructed in the Callovo-Oxfordian claystone. Journal of Rock Mechanics and Geotechnical Engineering (Article in preparation for the submission in an international journal). Conférences internationales avec article Tran N. T., D. P. Do, D. Hoxha, and M. N. Vu. (2020). Reliability-based design of deep tunnel excavated in the viscoelastic Burgers rocks. In Geotechnics for Sustainable Infrastructure Development, 375-382. Springer. https://doi.org/10.1007/978-981-15-2184-3_48.

Table 1 -

 1 1: Limit state functions commonly used in literature ................................................. 32 Table 1-2: Recommendations for model choice and use (Simpson et al., 2001) ..................... 41 Table 1-3: General steps of AK-based methods in reliability analysis ................................... 43 Table 1-4: Association of methods and types of problems ...................................................... 46 Table 1-5: Classical Orthogonal Polynomials Distribution (Sudret, 2012) ............................. 47 Table 1-6: Linear viscoelastic models with the primary forms of stress-strain relationship ... 54 Table 2-7: Properties of Burgers rock and the two liners ....................................................... 65 Table 2-8: The average results of example a series of four branches system ......................... 83

Table 2 -

 2

9: The average results of the tunnel excavated in the viscoelastic Burgers rock ...... 85 Table 3-10: Confining pressure and deviatoric stress of seven tri-axial creep tests performed on COx rock ............................................................................................................................. 99 Table 3-11: Statistical values of the Lemaitre's coefficients calibrated from seven creep tests on COx rock ............................................................................................................................. 99 Table 3-12: Correlation matrix between the calibrated coefficients of Lemaitre's model ...... 99 Table 3-13: Materials properties of the compressible layer................................................... 101 Table 3-14: Formation of the variable interest ..................................................................... 107 Table 4-15: The expansion order corresponding to the correlation lengths (adaptive to variance error ) .................................................................................................................................... 120 Table 4-16: Exceedance probability of concrete liner regarding the correlation length of four random fields of COx rock properties (concrete liner thickness l2=0.5m) ............................ 126 Table 4-17: Comparison of the AK-HDMR1 and AK-HDMR2 methods with the MCS, EGA, classical and modified AK-MCS methods (Tunnel excavated in viscoelastic Burgers rock)133

Table 4

 4 Isotropic variability problem) ............................................................................................... 134

	-18: Comparison of the AK-HDMR2 methods with the modified AK-MCS methods
	(Uncertainty problem) ............................................................................................................ 134
	Table 4-19: Comparison of the AK-HDMR2 methods with the modified AK-MCS methods
	(

Table 1 -

 1

1: Limit state functions commonly used in literature Limit state Limit state function (LSF) Variable definitions Support capacity

Table 1 -

 1

2: Recommendations for model choice and use (Simpson et al., 2001) 1.2.3.2. Active learning and Kriging-based methods Kriging model (or Gaussian Process model), an exact interpolative metamodeling technique, has increasingly received attention and interest in the last decade

Table 1

 1 NMCS) from the built Kriging model and calculate the failure probability Pf (e.g., using Eq. (1.23)).

		-3: General steps of AK-based methods in reliability analysis
	1	Generate initial DoE (i.e., initial training points) so-called matrix S.
	2	Determine the vector y (i.e., solve the deterministic problem to assess the exact response function for the training points of the initial DoE)
	3	Generate NMCS random samples for the interpolation by MCS.

4

Construct a Kriging model from S and y (by using a Kriging toolbox).

5 Interpolate g(x (i) ) (i=1, ...,

Table 1 -

 1

4: Association of methods and types of problems 1.2.3.3. Polynomial chaos expansions (PCEs) PCEs are a particular form of polynomial trend function that uses polynomials that are orthogonal w.r.t the probability density function (PDF) of the model parameters (Xinfeng Gao

Table 1 -

 1

5: Classical Orthogonal Polynomials Distribution

[START_REF] Sudret | Meta-models for structural reliability and uncertainty quantification[END_REF]

  y1) and (x2, y2) are the spatial coordinates of two points (2-D domain); 𝜃ℎ and 𝜃v are correlation lengths, respectively denote the horizontal-, and the vertical-correlation distance of soil properties. Notably, the relation between the scale of fluctuation

	where (x1, (  i	, ) i h v and the 
	correlation length can be expressed as	(	, )
	x y x y		y y	(1.48)

Table 1 -

 1 6: Linear viscoelastic models with the primary forms of stress-strain relationship

Table 2 -

 2 7: Properties of Burgers rock and the two liners Case of tunnels in the viscoelastic Burgers rock supported by an elastic liner Initially, it needs a specific example of a deep tunnel being dug in the viscoelastic Burgers rock and supported by a single elastic lining for the validation of our proposed analytical solution.

	Parameters	GM (GPa)	GK (GPa)	M (GPa.s) 	K (GPa.s) 	EL1 (GPa)	L1	EL2 (GPa)	L2
	Value	3.447	0.3447	9 4.137 10 	7 2.068 10 	16.547	0.2 39.1	0.2
	2.3.2.

Table 2 -

 2 9: The average results of the tunnel excavated in the viscoelastic Burgers rock

Table 3 -

 3 10: Confining pressure and deviatoric stress of seven tri-axial creep tests performed on COx rockTable 3-11: Statistical values of the Lemaitre's coefficients calibrated from seven creep tests on COx rock

	No. Samples	1	2	3	4	5	6	7
	Confining pressure (MPa)	2	2	2	2	12	12	12
	Deviatoric stress	11.5	17.25	20.7	20.7	18.5	27.75	33.3
	(MPa)
			Coefficient A	Coefficient B	Coefficient C
	Mean		4.88  10 -6		2.048		0.171
	Standard deviation (std)		4.27  10 -6		0.389		0.013
	COV (std/mean)		87.5%			19%		7.6%

Table 3 -

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

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

Table 3

 3

	-4, and the features of the concrete liner C60/75 (l2 =0.5m,
	E2=39.1 GPa, 2 0.2   ), the deconfinement rate (=0.85), coefficients A, B, C get its mean
	values from Table 3-2, Young modulus (Erocks=6 GPa) are fixed.

Table 3 -

 3 13: Materials properties of the compressible layer

  Base on non-accumulation of eigenvalues  j , around a non-zero value, one can order them in a descending series converging to zero. Let us consider the spectral decomposition of the

		j H x				Cov H	j	Cov H	x	H	x	j	(4.4)
	covariance matrix  Σ :
	 Σ j j j j    		1,... ,	N	(4.5)
	This equation allows one to transform the original vector linearly χ :
								N
	χ	( ) 	1      j μ	( )    j j  j	(4.6)
	where  	j	,	j		1,...,	 N are independent standard normal variables. Substituting for Eq. 4.6 in
	(4.3) and solving the OLE problem in Eq. (4.5) yields the EOLE representation of the random
	field:

Table 4 -

 4 15: The expansion order corresponding to the correlation lengths (adaptive to variance error )

		Isotropy cases			Anisotropy cases
	No.	x = y (m)	N (5%)	N (20%)	No.	x (m)	y (m)	N (5%)

Table 4 -

 4 16: Exceedance probability of concrete liner regarding the correlation length of four random fields of COx rock properties (concrete liner thickness l2=0.5m)

	Study case	Number of	Ncall	Pf (%)
		variables (4*N)

  Evaluation of x * on G and update the DOE. Here, instead of update np metamodels, one can identify only one-dimension metamodel to be enriched. This corresponds to the one with the highest variance at the point x * , as follows:

	( ) * argminU( ) ( )  k k x x x			k		1,...,	p n	(4.26)
	(10) Stop criterion on learning: if this criterion is satisfied, the learning is completed (Step 12).
	Otherwise, the method continues with Step 11.
	(11) 2 * arg max ( , ), *   u  i i G i i x	i		1,...,	n	p		(4.27)
	(12) Checking the condition of reliability analysis by MCS with the formula COV Pf < 0.05
	If this criterion is met, the learning is completed (End). Otherwise, the method continues with
	Step 13

(13) X is updated by a new population as in Step 1 and go back to Step 8.

Table 4 -

 4 17: Comparison of the AK-HDMR1 and AK-HDMR2 methods with the MCS, EGA, classical and modified AK-MCS methods (Tunnel excavated in viscoelastic Burgers rock)

	Method	Type Distribution of	Ncall	Pf (%)
		initial DoE

Table 4 -

 4 18: Comparison of the AK-HDMR2 methods with the modified AK-MCS methods (Uncertainty problem)

	Method	Type Distribution of	Ncall	Pf (%)	Pf(%)
		initial DoE
	Modified AK-MCS	Quasi-uniform LHS	81	37.72	-
	AK-HDM2	Quasi-uniform LHS	336	36.71	-1.09

Table 4

 4

	-19: Comparison of the AK-HDMR2 methods with the modified AK-MCS methods (Isotropic
		variability problem)
	Method	Type Distribution of	Ncall	Pf (%)	Pf(%)
		initial DoE
	Modified AK-MCS	Quasi-uniform LHS	112	50.84	-
	AK-HDM2	Quasi-uniform LHS	364	58.77	15.6

Table A1 -

 A1 3. Analytical solutions for LDP calculation depending on the medium[START_REF] Paraskevopoulou | Analysis of time-dependent deformation in tunnels using the Convergence-Confinement Method[END_REF]

  This last equation presents in effect the standard integral equation as shown in (Eq. (2.24)) with the free term 1 ( )

	t	d		t		d	f t
						(A2-12)
	e	e	G t		d	f t
	B f t and the Kernel	1	1 exp	(	)

  A.2.2. Determination of supporting pressure after the installation of the second linerThe determination of the supporting pressure in the liners after the second installation stage consists of calculating the pressures 12 ( ) p t and 22 ( ) p t with 2 t t  .Using the radial displacement in each liner expressed in Eqs.(2.15) and (2.18), the compatibility condition of the displacement at the interface of two liners (see Eq.(2.23)) can be rewritten in the form:

														B t	takes the same form of Eq. (2.27), which means
	that:
	1 B 	1 B t f t ( ) ( ) 		1 t t 		1   ( ) t Fe 1 		2   ( ) t F e 2 		1 B f	( ) , d  	(A2-13)
	Correspondingly, the solution of	  p t can be calculated: 11
	  1 B p t 11  	( )exp t	  		K K G 	t	  	,	(A2-14)
	10 12 ( ) a p t a p t a p t t t 11 22 10 12 2 2 ( ) ( ), ( )   	(A2-15)

  ) of the characteristic equation Eq. (2.25). Thus, the solution of Eq. (A2-26) has the same form as one presented in Eq. (2.28), meaning that:

	which has the Kernel 2        K e e 2 exp M	  	( G t  K K			)	  	   	and the free term 2 ( ) B f t .
	This				integral							equation	with	the	corresponding	parameters
	1 E	2  K   e	,	2 E	2 M    e		,		1 		0,	2 	K   K G	results from a positive discriminant (	0
					2 B 	2 B t f t ( ) ( ) 		2 t t 		1   ( ) t Fe 1 		2   ( ) t F e 2 		f	2 B	( ) , d  	(A2-27)
	From this last solution of	  B t 2  , one can deduce the solution of   12 p t :
					  2 B p t 12  	( )exp t	  		K K G 	t	  	,	(A2-28)
	f t		e p G	t					t			G	t	e p R	G	t t	A A
						e p R	A A A R	G	d	e	B	(A2-25)
						B	G	t					e p	G	t	d	M e	B	  	G K	t	  	,
	We also obtain the following standard integral equation:
		2 B	( ) t		2 t t 	   	2    K e e 2 M	exp	  	 G K K	( t		2     ) ( )        B d		2 B f t ( ),	(A2-26)

a D 

Table A1-6. Fundamental equations for probability theory

years. Also, at this time, I highly appreciate Laboratory Lamé and The French National Radioactive Waste Management Agency (Andra) for their additional backing.

A.1.2. Appendix for reliability analysis methods and stochastic stability analysis ... The following discussions will focus on the general view of four mechanical parameters' effects on the maximum equivalent stresses in the concrete liner. The previous investigations exhibit a significant difference in the behavior of concrete liner of drift by using different mechanical properties of COx rock, respecting its mean values, and utilizing the Set-1. Here, the broader view of these properties' effects can be stated in Figure 3-22, where the variation of the interested stresses associated with each parameter of COx rock is illustrated. These results are extracted from the direct evaluations of the performance function taken from a reliability analysis performed in the next part. These simulations were conducted on about 500 MCS sets following the Normal distribution of rock properties. In which, we fixed values of concrete liner thickness (l2=0.5 m) to confirm the complexity of the considered problem. Due to the nonlinear behavior of host rock with the eventual correlation of its long-term properties, an evident

The following analyses aim at providing the impact of each input material parameters by taking place two-cases for each sub-case-study. The first case (so-called Case-1) is accompanied by one coefficient that varies, and others will be set at its mean values. While the second one (socalled Case-2) also allows only one parameter changes, yet others will be taken from the Set-1, as in section 3.4.2 above, and fixed 2  v = 0.4. The values of a variable interest in each study can be applied from the formation illustrated as in Table 3-14 below, in which the minimum and maximum values are based on 95% CI (Confidence interval), their average (mean) and standard deviation (STD) obtained from section 4.3.2.

where X denotes for one variable in each case study, i.e., X can be coefficient A, B, C, or Young modulus E.

Influence of the coefficient A

In this section, the coefficient A tends to affect the maximum equivalent stress in the concrete liner proportionally. However, in both cases (using the Set-mean set and the Set-1), there again appears two sharp changes (at the beginning for the Set-mean, and at about A=0.8*10 -5 ). The clear differences in the two figures show the significant impact of coefficient A on the interested stress. Notably, in Figure 3-24, with respect to a higher value of coefficient A, the deformations on the compressible zone exceed the limitation 2  v and induce a pronounced increase of the maximum equivalent stress in the concrete liner. However, in this case, the influence of the adopted value of lim2 is also investigated regarding the fact that the degradation mechanism of concrete (which is not explicitly taken into account) may be significantly varied at the time scale of 100 years.

Figure 3-31a illustrates a representative evolution of the drift's final support's exceedance probability during the iterative construction of the metamodel. Corresponding to the fixed values of the concrete liner's thickness (l2=0.5m) and the compressible potential of the compressible layer ( 2 0.515 v   ), 91 iterations (i.e., 388 total number (Ncall) of direct evaluations of the performance function) are carried out. The estimated probability, Pf=0.11%, is obtained by adopting the allowable stress lim2= fcd=36MPa. An investigation with the other values of lim2 as illustrated in Figure 3-31b presents a decrease as expected of the exceedance probability in the concrete liner when its compressive strength (threshold) increases. The probability could reach to about 95% for the case lim2=16MPa while it is smaller than 1% when the allowable stress is higher than 25MPa. The figures show that the number of direct evaluations of the performance function (Ncall) is more critical for the smaller exceedance probability. in the first liner of the drift, we now investigate the influence of the two-liner thicknesses on the long-term stability of the solid support . For the variation of the compressible layer thickness (l1) (Figure 3-32a with fixed l2=0.5m), the probability of exceedance in the concrete liner at 100 years changes in the range 6. 81% to 0.11% when the thickness increase from 0.1m to 0.2m. For the remaining case (Figure 3-32b with fixed l1=0.2m), the exceedance probability varies from 99.3% (for the case of thickness l2=0.2m) to the value of 0.11% (case l2=0.5m). Indeed, it is smaller than 1% when the thickness of this final support is higher than 0.35(m). These figures illustrate that the final thickness has much more impact than the rest.

Numerical application in case of the deep tunnel in the viscoelastic Burgers rock

We investigate firstly the applicability of the AK-HMDR1 and AK-HMDR2 in the simple case of uncertainty problem of tunnel excavated in the viscoelastic Burgers rock. Note that this case was deeply studied in the second chapter by comparing the obtained results provided by the MCS, classical AK-MCS, and modified AK-MCS methods.

Note that, in this problem, the vector of random variables consists of four parameters of Burgers model x = [GM, M, GK, K]. By adopting the AK-HDMR1 for this random vector, one needs to construct four Kriging metamodels as follow:

whilst for the HMDR2, a total of ten metamodels must be taken place:

. Some deterministic analysis methods for underground excavations

This section mentions the deterministic approach often used in underground excavations, based on partial safety factors, i.e., it does not consider qualitative uncertainty. This step is necessary to quickly understand the core of interesting problems before deeply studying the probabilistic approach or reliability analysis: quantitative propagation of uncertainty on the probability of failure as in the later sections.

A.1.1.1. Empirical approach Since the accumulation of the first-hand experience of engineering practice, the empirical method has been developed along with the engineering judgments [START_REF] Huaxin | Analysis of underground rock excavations accounting for uncertainty[END_REF]. The method is widely used in engineering practice with a convenient application, especially for describing the domain's rock mass quality and/or support requirements. There are three common empirical methods in underground excavation problems based on the ability to quality the rock mass and guide the adaptable support selection, namely: Rock Mass Rating (RMR), Q systems, and Geological Strength Index (GSI).

Rock Mass Rating (RMR): The RMR system was first developed by Bieniawski [START_REF] Bieniawski | Engineering classification of jointed rock masses[END_REF] and developed through several revisions until the latest version in 1989 [START_REF] Bieniawski | Engineering rock mass classifications: a complete manual for engineers and geologists in mining, civil, and petroleum engineering[END_REF]. The common applications of RMR in underground excavations are: (1) for Rock mass classification, (2) estimation of rock mass strength, (3) estimation of the stand-up time and the maximum unsupported rock span, (4) preliminary selection of rock supports (see Table A1-1 and Table A1-2). Note that the recommendation for tunnel supports is constrained to 10 meters span only.

Q systems: The Q system was first found by [START_REF] Barton | Engineering classification of rock masses for the design of tunnel support[END_REF] initially intended for the support design for tunnels in hard rock. The effective applications of the Q-system can be listed as: (1) for Rock mass classification, (2) estimation of rock mass strength, (3) selection of underground excavation support [START_REF] Barton | Estimation of support requirements for underground excavations[END_REF] (see Table A1-2). Its advantages and disadvantages are: The Q-system works best in ground conditions where block probably falls; the application of Stress Reduction Factor, a vital parameter of Q-system, is unclear for buckling, rock burst, and squeezing conditions. The Rock Quality Designation, the first parameter of Q-systems provided by [START_REF] Barton | Some new Q-value correlations to assist in site characterisation and tunnel design[END_REF], has several limitations in characterizing the degree of jointing. As reported by Pells [START_REF] Pells | Developments in the design of tunnels and caverns in the Triassic rocks of the Sydney region[END_REF], the Q system predicted significantly less support than actually adopted in several caverns he studied in Australia [START_REF] Cai | Rock mass characterization and rock property variability considerations for tunnel and cavern design[END_REF].

Geological Strength Index (GSI):

The GSI, proposed by Hoek [START_REF] Hoek | Strength of rock and rock masse[END_REF], estimates the rock mass quality or engineering parameters based on the field observation. The obviously main application of GSI lies in its role as an important parameter in the Hoek-Brown (H-B) failure criterion (Hoek et al., 2002), besides other engineering parameters such as Mohr-Coulomb or rock mass modulus [START_REF] Cai | Rock mass characterization and rock property variability considerations for tunnel and cavern design[END_REF]. The fact that GSI also can be used as a classification method, i.e., high-quality rock masses in the range 60 < GSI < 75, average quality rock masses in the range 40 < GSI < 50, low-quality rock masses in the range 25 < GSI < 30 (Alfonso and Veiga, 2012). Where: Fs: objective function (the factor of safety), computed by deterministic analysis;

mrslope2D: An open-source by [START_REF] Fenton | Risk assessment in geotechnical engineering[END_REF];

APPENDIX 2: APPENDIX TO CHAPTER 2

In this appendix, we detail the process to determine the supporting pressures of two elastic liners by using the compatibility conditions of displacement at the interface between the rock mass and the first liner (Eq. (2.22)) and at the interface of two liners (Eq. (2.23)).

A.2.1. Determination of supporting pressures after installation of the first liner

By substituting the Eq. (2.12) and Eq. (2.15) in Eq. (2.22), the incremental radial displacement of rock from time t1 to the generic time t > t1 can be written in the integral form:

with:

Regarding the conditions expressed in Eq.(2.5), the Eq. (A2-1) written for the instant t with

This Eq. (A2-3) results from the second type Volterra equation [START_REF] Polyanin | Handbook of integral equations[END_REF] for

Substituting Eq. (2. into Eq. (A2-4) yields the following equation:

which can be expanded in the form:

, the following integral equation can be obtained:

Then by defining 1

e RG G a R , we can rewrite Eq. (A2-7) in the form:

the Eq. (A2-9) can be rewritten as:

Thus, the function 1 ( ) B t  in Eq. (A2-8) can be simplified as:

The Eq. (A2-1) written for the instant t (with 2 t t  ) has the following form:

Substituting ( ) H t from Eq. (2. into Eq. (A2-18) we have: 2( )

, after some developments, the simplified form of the integral equation 2 ( )

 can be obtained:

where: The coupled mechanical probabilistic model of the tunnel excavation problem is investigated the reliability analysis with some following main steps:

1. Generating the centroid of the rock cells and its ID of the physical mesh with the scripts written by PYTHON. The output will be in the plain-text-form with the 2-D coordinates of the ID and the labels of all cells of the rock mass mesh.

2. Generating the random fields (coefficients: A, B, C, E) by MATLAB software.

2a. Selecting the boundary dimension for the domain (i.e., considering the acceptance of the extension dimension (lac) from the original mesh to artificial mesh where the variance error must satisfy a specified threshold for the accuracy of the approximation (e.g., 5%).

2b. Selecting the suitable value of the partition (Np).

2c. Selecting the combination of truncated order expansion (N) and the correlation lengths (xy).

This step creates the MAT-file including some necessary parameters: probabilistic values of 4 random input fields (i.e., the mean and the Std); covariance matrix, the corresponding eigenvalues, and eigenvectors (acts as the matrix S for the initial DoE or matrix X for NMCS) 3. For the case of initial DoE, we transmit the matrix S to the rock elements by using Eq. (4.15) for coefficients A, B, and C, Eq. (4.7) for modulus E. The output is stored in the text files as the source of input material data. Then Code-Aster can read and assign them to every cell of the rock mass. After each numerical mechanical computation, the post-processing of the mechanical model response is performed using MATLAB. This software can again connect to Code-Aster to call the next computations (if the iterations are needed).

4. The response is stored in MATLAB and used by the modified AK-MCS being mentioned in Chapter 2 for reliability analysis. Notice that the active learning and the stopping conditions are the same as we investigated in Chapter 3. First, we assume some other prerequisites for the discretization: the truncated order expansion Mterm=20, the isotropic correlation length xy=15m, and partition value Np = 40. Table A3-1 above shows that the extension dimension lac=3m can be the best choice. In which the maximum of variance error is lower than the prescribed tolerance er = 5%, and the maximum variance error, although higher than er, but occurs in the points significantly closer to the boundaries (all points in the zone of physical mesh [0.0, 55.0] have the error of variance satisfying the requirement of this condition, i.e., lower than er). The best choice of the number depends on the maximum variance error as well as the timeconsuming for each computation. Although all these maximum variance error values are still

A.3.2. Selecting prerequisites for the discretization

A.3.2.1. Boundary dimensions for the domain