[Objective] Accurate prediction of tunnelling-induced vertical responses in pile foundations remains a critical challenge in urban underground construction. Traditional deterministic analyses of the complex tunnel-soil-pile interactions often assume homogeneous soil properties, neglecting the inherent spatial variability of soil properties. Such simplifications may result in underestimating or misrepresenting pile responses. To bridge this gap, the primary objective of this study is to develop and implement a sophisticated probabilistic model capturing the vertical spatial variability of soil properties. This model aims to facilitate a comprehensive stochastic analysis and provide more realistic and reliable predictions of pile behavior due to adjacent shield tunnel excavation. [Methods] The core deterministic framework employed a well-established two-stage analytical procedure: first, tunnelling-induced free-field ground movements were modeled using the Loganathan-Poulos solution, which accounted for volume loss and tunnel geometry effects on surrounding soils; second, pile foundation responses to these soil displacements were evaluated through Load Transfer Analysis to calculate pile head settlements and axial force distributions along the pile shaft. Undrained shear strength was modeled as a random field to account for the vertical spatial variability. The two-stage deterministic procedure and the vertical random field model for undrained shear strength were integrated within an automated Monte-Carlo simulation framework, which constituted the developed stochastic Two-Stage Analysis model. After a sufficiently large number of iterations, the recorded pile responses were statistically analyzed, yielding the probability distributions, mean values, and standard deviations for the pile head additional settlement and the maximum additional axial force along the pile shaft. [Results] Both the coefficient of variation and the correlation length of the undrained shear strength exerted a significant influence on the statistical moments of the pile responses. An increase in the coefficient of variation generally led to higher mean values of pile head additional settlement and the mean maximum additional axial force. The correlation length caused significant variations in both the mean and, particularly, the standard deviations of pile settlements and axial forces compared to the homogeneous case. The probability of pile head additional settlement exceeding a critical serviceability limit state showed a strong dependence on the vertical spatial variability of soil strength. As the overall variability increased, the computed probability of settlement failure rose significantly. The detrimental impact of soil spatial variability on pile reliability amplified by higher levels of deterministic loading factors. Increased ground deformation caused by larger tunnel volume loss intensified the negative effects of soil variability. In addition, existing structural loads applied to the pile head further magnified the sensitivity of pile responses and the associated failure probability to the underlying soil uncertainty. Essentially, the combined loading conditions made the pile foundation more vulnerable to the adverse consequences of spatially variable soil properties. [Conclusion] (1) The proposed stochastic two-stage model successfully bridges the gap between conventional analytical methods and real-world soil heterogeneity. (2) The developed Monte-Carlo automation program can provide practical guidelines for prioritizing vertical variability characterization in geotechnical investigations and for adopting probabilistic design methods instead of conventional safety factors. (3) The current model neglects horizontal soil variability and construction disturbances. Future work should integrate 3D random fields with machine learning techniques for enhanced prediction.