There is an increasing interest in utilizing surrogate models to perform the global sensitivity analysis (GSA) for the uncertainty quantification since they allow to quantify efficiently the relative importance of each input parameter with relatively small computational costs. This paper aims at investigating the feasibility of applying this method to a stochastic soil-tunnel system in seismic conditions. An advanced technique named sparse polynomial chaos expansions (SPCE) is introduced to build a surrogate model that allows performing a GSA combined with the Sobol’ indices. The accuracy and efficiency of the SPCE-GSA method are validated by comparison with the true numerical predictions and Monte Carlo simulations. Parametric sensitivity analyses are performed for a wide range of the soil shear wave velocities, ground motion intensities, probability distribution types, and coefficients of variation, assuming a nonlinear soil behavior. The influences of the sampling size, polynomial degree, and ground motion characteristics on the variability in the sensitivity indices are evaluated. The results show that the soil shear wave velocity and modulus reduction factor are the two variables that significantly affect the seismic deformations of tunnels.