The ability to predict future variability of groundwater resources in time and space is of critical importance in society’s adaptation to climate variability and change. Periodic control of large scale ocean-atmospheric circulations on groundwater levels proposes a potentially effective source of longer term forecasting capability. In this study, as a first national-scale assessment, we use the continues wavelet transform, global power spectrum, and wavelet coherence analyses to quantify the controls of the Atlantic Multidecadal Oscillation (AMO), Pacific Decadal Oscillation (PDO), North Atlantic Oscillation (NAO), and El Niño Southern Oscillation (ENSO) over the representative groundwater levels of the 24 principal aquifers, scattered across different 14 climate zones of Iran. The results demonstrate that aquifer storage variations are partially controlled by annual to interdecadal climate variability and are not solely a function of pumping variations. Moreover, teleconnections are observed to be both frequency and time specific. The significant coherence patterns between the climate indices and groundwater levels are observed at five frequency bands of the annual (~1-yr), interannual (2-4- and 4-6-yr), decadal (8-12-yr), and interdecadal (14-18yr), consistent with the dominant modes of climate indices. AMO’s strong footprint is observed at interdecadal and annual modes of groundwater levels while PDO’s highest imprint is seen in interannual, decadal, and interdecadal modes. The highest controlling influence of ENSO is observed across the decadal and interannual modes whereas the NAO’s footprint is marked at annual and interdecadal frequency bands. Further, it is observed that the groundwater variability being higher modulated by a combination of large-scale atmospheric circulations rather than each individual index. The decadal and interdecadal oscillation modes constitute the dominant modes in Iranian aquifers. Findings also mark the unsaturated zone contribution in damping and lagging of the climate variability modes, particularly for the higher frequency indices of ENSO and NAO where the groundwater variability is observed to be more correlated with lower frequent climate circulations such as PDO and AMO, rather than ENSO and NAO. Finally, it is found that the data length can significantly affect the teleconnections if the time series are not contemporaneous and only one value of coherence/correlation is computed for each particular series instead of separate computations for different frequency bands and different time spans.