## 1 Introduction  
Advanced techniques have been largely developed to ensure the stability at the short- and long-terms of tunnels. For example, the innovative support system using the double flexible/concrete liner of tunnel drilled by TBMs in the complex geological conditions has been proposed in tunneling projects [1]. However, the complexity of the geological formations' behavior combined with their high uncertainty remains a great challenge on the stability prediction of tunnel that is usually excavated in the anisotropic condition of initial stress state [2; 3; 4; 5; 6]. In the last two decades, the probabilistic approach has been largely accepted as the most performant method to account for explicitly the uncertainty effect instead of using different factors of safety as in the traditional method. This stochastic approach allows quantifying the uncertainty propagation from input parameters (e.g., geological formations' behaviors) to a structure's responses (e.g., maximum stress), or estimating the structure's failure probability (_Pf_) with respect to one or several failure modes. Some well-known methods, such as the Monte Carlo simulation (MCS), the first/second order reliability method (FORM/SORM), or the response surface method (SRM) have become standard and have been widely used in different tunnel projects [3; 4; 7; 8]. So far, the reliability analysis of tunnels was mainly performed in the geological formations whose behaviors are assumed to be time-independent and characterized by an elastoplastic law. For example, ref. [4] proposed an approach to assess the system reliability of a tunnel verifying three failure modes (i.e., support capacity, tunnel convergence, and rock bolt length) in which the uncertainty of the Hoek-Brown parameters was accounted for. In many other contributions, scholars considered the uncertainty of the elastic perfectly plastic soils/rocks using the Mohr Coulomb model to analyse the stability of underground structures [3,8] and of the tunnel face [7]. The investigations conducted in these studies elucidated the important effect of the elastic moduli and strength of geological formations on the failure probability of constructed undergrounds. However, the adopted assumption of time-independent behavior can underestimate the obtained results, especially in the tunnel projects excavated in the host rock exhibiting a significant creep behavior. It is worth to note that in some cases, the time-dependent behavior of rock can contribute for more than 70% of the total convergence whilst the evolution in time of convergence (i.e., the well-known time delay effect) can significantly affect the stability of tunnel supports [9, 10, 11]. Thus, the stability of tunnels must also be evaluated at the long-term (i.e., at the end of their service life of about 100 years), especially for deeply buried tunnels constructed in the rock masses owing a high time-dependent behavior. This issue was extensively considered in the literature, but the uncertainty of creep behavior of the geological formations is usually neglected. For example, adopting some simple rheological models to describe the creep phenomenon of viscoelastic rocks, different scholars derived closed-form solutions to estimate displacements in rock masses and stress states in tunnel supports [12, 13, 14, 15, 16]. These solutions provided an interesting tool for quick analyses and parametric studies of tunnel responses at the long-term, which are usually expensive and time-consuming in numerical simulations. From their studies, these authors demonstrated the significant effects of creep rocks when an evolution in time of convergence on the tunnel surface induces a quick increase of compressive stresses in liners.