The traditional deterministic analysis for tunnel face stability neglects the uncertainties of geotechnical parameters, while the simplified reliability analysis which models the potential uncertainties by means of random variables usually fails to account for soil spatial variability. To overcome these limitations, this study proposes an efficient framework for conducting reliability analysis and reliability-based design (RBD) of tunnel face stability in spatially variable soil strata. The three-dimensional (3D) rotational failure mechanism of the tunnel face is extended to account for the soil spatial variability, and a probabilistic framework is established by coupling the extended mechanism with the improved Hasofer-Lind-Rackwits-Fiessler recursive algorithm (iHLRF) as well as its inverse analysis formulation. The proposed framework allows for rapid and precise reliability analysis and RBD of tunnel face stability. To demonstrate the feasibility and efficacy of the proposed framework, an illustrative case of tunnelling in frictional soils is presented, where the soil's cohesion and friction angle are modelled as two anisotropic cross-correlated lognormal random fields. The results show that the proposed method can accurately estimate the failure probability (or reliability index) regarding the tunnel face stability and can efficiently determine the required supporting pressure for a target reliability index with soil spatial variability being taken into account. Furthermore, this study reveals the impact of various factors on the support pressure, including coefficient of variation, cross-correlation between cohesion and friction angle, as well as autocorrelation distance of spatially variable soil strata. The results also demonstrate the feasibility of using the forward and/or inverse first-order reliability method (FORM) in high-dimensional stochastic problems. It is hoped that this study may provide a practical and reliable framework for determining the stability of tunnels in complex soil strata.