In this study, a modified pseudo-dynamic Bishop method (MPDBM) is developed incorporating spatial variability in peak ground acceleration (PGA) and initial phase by way of random field theory. The proposed MPDBM is validated against existing results and illustrated through a large-scale soil slope emphasising the effect of scale of fluctuation on the failure probability of seismic slope stability. The numerical simulations reveal that neglecting the spatial variability of PGA and the initial phase results in a larger failure probability. The maximum increment in the failure probability can be up to almost two orders of magnitude. Ignoring the spatial variability of the PGA and the initial phase tends to overestimate the failure probability or underestimate the seismic slope stability, resulting in a conservative design. The failure probability is more sensitive to exposure time as the spatial variability of the PGA and initial phase becomes more significant.