Groundwater resources in arid and semi-arid regions of Iran, including the Qazvin aquifer, are under increasing pressure due to overexploitation and climate change. The present study aimed to conduct a multivariate assessment of pollution in the Qazvin aquifer using the Groundwater Quality Index (GWQI) based on eleven parameters: Na⁺, K⁺, Mg²⁺, Ca²⁺, SO₄²⁻, Cl⁻, HCO₃⁻, pH, TDS, EC, and TH. Six copula families were employed, including Gaussian, Student’s t, Clayton, Frank, Joe, and AMH copulas. The data used in this study consisted of records from 28 piezometers across the Qazvin Plain for the period 2006–2022. Given the availability of data from different piezometers, copula functions were applied to analyze the spatial-temporal behavior of the pollution index across the Qazvin Plain. The results showed that the AMH, Clayton, and Frank copulas were the best-performing models for simulating the dependence structure of the GWQI among the piezometers. The joint return period TAND was estimated at 11.11 years, compared with TOR, which was 1.73 years. This represents a 6.4-fold difference and demonstrates that the univariate approach systematically underestimates risk. The findings indicate that conditional copula functions can be used for probabilistic spatial-temporal analysis of groundwater pollution risk in aquifers.