The research on the dynamic response of slopes under seismic action is currently focused on the action of a single mainshock without considering the effect of aftershocks, and the spatial variability of material parameters is always ignored. In this paper, the spatial variability of parameters is fully considered, and a reliability analysis framework based on the Newmark method and probabilistic density evolution method (PDEM) is proposed to quantify the effects of aftershocks and spatial variability on the dynamic reliability. First, the physical random function model, Copula function and narrowband harmonic group superposition method are combined to generate the mainshock-aftershock sequence (MAS). In addition, the random field is generated based on the spectral representation and the parameters are assigned to the finite element model based on the corresponding coordinates. Then, the permanent displacement of the slope considering the spatial variability of parameters subjected to the MAS were batch calculated based on the Newmark method, and the effects of the coefficient of variation of cohesion and friction angle (COVC and COVF), aftershock, and peak ground acceleration (PGA) on the permanent displacement of the slope were analyzed by the mean of the displacement. Finally, based on PDEM, the effects of COV and aftershocks on the dynamic reliability of the slope are explained from a probabilistic point of view. The results of the study show that the mean of displacement shows a gradual increase with the increase of coefficient of variation (COV). In contrast, the COVF has a more pronounced effect on the displacement of the slope. In addition, the mean of displacement of the slope subjected to the MAS are larger compared to the single mainshock. If the spatial variability of parameters and the influence of aftershocks are ignored, the seismic performance of the slope will be overestimated.