Throughfall temporal stability is critical for maintaining spatial consistency in soil water amount and nutrient inputs, serving as a fundamental pillar for the structural and functional stability of vegetated ecosystems. However, traditional time stability plots rely on subjective visual interpretation. Quantitative alternatives like the relative difference and Spearman rank correlation coefficient methods fail to differentiate the specific contributions of distinct throughfall regimes (dry vs. wet zones) to overall stability. This study developed an integrated framework coupling Gaussian Mixture Model (GMM) with geostatistical analysis and statistical metrics, including the coefficient of variation in ranks (CVr) and global Moran's I. Validated using field data from two dominant shrubs (Caragana korshinskii Kom. and Salixpsammophila C. Wang & C. Y. Yang) in the semi-arid Loess Plateau of China, this framework accurately quantified species-specific throughfall spatiotemporal patterns: S. psammophila exhibited greater spatial heterogeneity and temporal stability along direction and distance gradients. Spatial autocorrelation examination via Moran's I confirmed significant clustering in throughfall distribution (p < 0.05). GMM-identified dry zones, receiving significantly less throughfall (65.0% vs. 77.1%, p < 0.001), maintained higher stability evidenced by smaller CVr (23.4% vs. 68.5%, p < 0.01). They functioned as potential markers for the persistence of overall throughfall patterns, distinct from wet zones acting as biogeochemical hotspots. Therefore, this framework provides a mechanistic tool to decipher throughfall spatiotemporal pattern, enhancing the predictability of throughfall fluxes and offering new insights into ecosystem stability particularly in water-limited environments.