Rock properties are estimated using objective lab/in-situ testing and subjective judgements invoking different types of uncertainties, i.e., aleatory, and epistemic, along them, often indicated by varying information levels. This study presents a unified reliability method to integrate the spatial variability of inputs modelled via alternate uncertainty models (intervals and probability-boxes (p-boxes)) with those modelled as stochastic variables via probability distributions. The methodology employs advanced Karhunen-Loe`ve decomposition to generate interval and random fields of inputs modelled via alternate and stochastic models, respectively. Input properties are allocated to the zones of the numerical model in Fast Lagrangian Analysis of Continua-2D (FLAC2D) based on their spatial dependency and correlation functions through a developed MATLAB-FLAC coupled code. The methodology is demonstrated for a deep tunnel in Canada to be constructed along a massive rock prone to brittle failures. Intact rock properties are modelled as stochastic variables due to objective estimation, while Geological Strength Index (GSI) and deformation modulus are modelled using alternate models (interval and p-box, respectively) due to subjective and hybrid estimation (double uncertainty propagation algorithm). The results of the methodology are compared with those of traditional deterministic and random field methods. The methodology reduces the subjectivity invoked by including unavailable additional information (e.g., assuming probability distributions of inputs based on literature) and propagates the originally available information of inputs accurately. The final outputs are the p-boxes of response parameters (i.e., displacements and damage zone) instead of their fixed values (deterministic analysis) and probability distributions (traditional reliability analysis), indicating the propagation of both impreciseness and variability of inputs by the method. For this case study, the p-boxes of outputs were bounding their values/distributions estimated via traditional analyses, verifying the accuracy of the methodology. Further, the impreciseness in the outputs, highest in the damage zone extent, was due to imprecision in the estimated GSI.
