In drought frequency analysis, as the number of drought variables increases, the joint behavior between these variables needs to be studied. Therefore, this study aims to develop a flexible four-variate joint distribution function of the regional stochastic nature of drought. Using run theory, drought duration, severity, peak, and inter-arrival time were abstracted from the Standardized Precipitation Evapotranspiration Index (SPEI) aggregated at six months, observed in mainland China between 1961 and 2013. As these drought variables showed significant dependence properties and followed different marginal distributions, we employed and compared six four-variate symmetric and asymmetric Archimedean copulas (i.e., Frank, Clayton, Gumbel–Hougaard). The best-fitting model for each region was carefully selected using RMSE, AIC, and BIAS goodness-of-fit tests. Results revealed that the empirical and theoretical probabilities of the symmetric Clayton in regions NE (Northeast), CS (Central and Southern China), EMC (Entire China), and symmetric Frank in regions NC (North China), SC (South China), IM (Inner Mongolia), NW (Northwest), TP (Tibet Plateau) agreed well. Symmetric Frank copula was considered the best-fit for station-based drought analysis in EMC. Based on these copulas, the drought probabilities and return periods for the occurrence of drought events over the next 5, 10, 20, 50, and 100 years in each region were hereby comprehensively explained, and the results shown here could be helpful in the appraisal of the adequacies of water supply systems under drought conditions in all regions. This study showed that a four-variate copula approach is a vital tool for probabilistic interpretation of hydrological and meteorological data in the different climatic region of mainland China.