This study develops a reliability-based framework for predicting and optimizing tunnel stability in rock masses under surcharge loading while explicitly accounting for both aleatory and epistemic uncertainties. A unified dataset for twin circular and square tunnels is generated using Adaptive Finite Element Limit Analysis under the generalized Hoek–Brown criterion. The results demonstrate that probabilistic predictions obtained using Natural Gradient Boosting provide accurate stability estimates together with well-calibrated uncertainty bounds, consistently outperforming multiple baseline machine-learning models. Validation against more than 300 independent Optum G2 simulations confirms strong agreement with numerical benchmarks. A dedicated uncertainty decomposition analysis further shows that neglecting either input uncertainty or model uncertainty can lead to misleading and potentially unsafe reliability estimates, underscoring the necessity of joint uncertainty propagation. Overall, the proposed framework enables robust, uncertainty-aware tunnel design under reliability constraints and provides a practical decision-support tool for rock engineering applications.