Large-scale climate variability footprint in water levels of alluvial aquifers across Iran

The ability to predict the 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 serves as a potentially effective source for relatively long-term forecasting. In this study, as a first national-scale assessment, we use the continuous wavelet transform, global power spectrum, and wavelet coherence analyses to quantify the linkage between the Atlantic Multidecadal Oscillation (AMO), Pacific Decadal Oscillation (PDO), North Atlantic Oscillation (NAO), and El Niño Southern Oscillation (ENSO) and the representative groundwater levels of the 24 principal aquifers, scattered across 14 different climate zones of Iran. Aquifer storage variations are found to be partially derived from annual to interdecadal climate variability and not solely a function of pumping variability. Moreover, teleconnections are observed to be both frequency and time-specific. The significant coherence patterns between the climate indices and groundwater levels are concentrated at five frequency bands of the annual (~ 1 year), short-interannual (2–4 years), medium-interannual (4–6 years), decadal (8–12 years), and interdecadal (14–18 years), consistent with the dominant oscillations of climate indices. AMO’s strong linkage to groundwater variability is found to be at interdecadal and annual modes of groundwater levels while PDO’s highest imprint is concentrated in interannual, decadal, and interdecadal periodicities. Unlike ENSO for which the highest modulating influence is found to be across the decadal and interannual modes, the NAO’s footprint in aquifers is marked at annual and interdecadal frequency bands. Findings further show that the groundwater variability is driven primarily by a combination of multiple large-scale atmospheric circulations rather than a single atmospheric circulation index. The decadal and interdecadal oscillation modes are 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 so that groundwater variability is more correlated with those climate circulations that have lower frequency such as PDO and AMO. Finally, the data length is found to have a significant effect on the teleconnections if the time series are not contemporaneous and for each particular series, only one value of coherence/correlation is computed instead of separate computations for different frequency bands and different time spans.


Introduction
Since climate variability plays an essential role in the sustainability of water resources (Sadoff and Muller, 2009), the potential for the climate to impact water resource variations on local to global scales has been an increasingly crucial water management subject (Velasco et al. 2015). Periodic control of ocean-atmospheric circulations over groundwater variability offers a potentially important source for relatively long-term forecasting that is of critical importance to drought management (Rust et al. 2018). Despite the effect of interannual to multidecadal climate variability on the spatiotemporal distribution of precipitation, drought, groundwater, and surface-water storage variations (Ropelewski and Halpert, 1986;Cayan and Webb, 1992;Hurkmans et al. 2009;Sadoff and Muller, 2009;Vicente-Serrano et al. 2011;Kuss and Gurdak, 2014;Rezaei and Gurdak, 2020;Rezaei, 2020), it is still poorly understood how climate variability modulates subsurface hydrologic processes and groundwater storage worldwide (Hanson et al. 2006;Green et al. 2011;Kuss and Gurdak, 2014;Velasco et al. 2015;Rust et al. 2018). As the largest accessible freshwater resource, groundwater has a critical role in providing adequate water supplies for human consumption, agricultural and industrial purposes, and ecosystem function (Treidel et al. 2012). How subsurface hydrologic responds to natural climate variability is of particular interest in semiarid and arid climates in which groundwater resource availability and sustainability are crucial (Hanson et al. 2004). Groundwater supplies the primary source of freshwater in a broad part of Iran, where the groundwater pumping from production wells provides solely ~ 60% of the entire water supply across the nation (Mirzaei et al. 2019). Therefore, establishing a better understanding of the influence that climate variability has on groundwater resources in Iran with (semi) arid climate could help in improving sustainable development and management of groundwater and surface-water resources (Hanson et al. 2004;Kuss and Gurdak, 2014;Velasco et al. 2015;Rezaei and Gurdak, 2020).
In spite of the acknowledged importance of groundwater resources, few investigations have documented the linkage between large-scale climate indices and natural variability in Iranian aquifers. Rezaei and Gurdak's (2020) study was the first to quantify the effects of large-scale climatic variability on groundwater resources in Iran. They found that the Pacific-based climate oscillations (PDO and ENSO) have a more powerful influence than the Atlanticbased oscillations (NAO and AMO) on the variability in the groundwater storage surrounding Lake Urmia. Rezaei (2020) further explored the remote teleconnection between ENSO and PDO climate indices and a karst Sarabkalan Spring's discharge in the southwest of Iran with a fastslow flow system which is common in karst aquifers with a coupled conduit-diffuse flow regime. He observed that the enhanced precipitation and spring discharge of the Sarabkalan Spring occurred in connection with a positive PDO phase that overlapped with El Niño.
This work quantifies the linkage between AMO, PDO, NAO, and ENSO and the groundwater levels of 24 principal aquifers scattered over the country. Iran, as a country with a highly variable climate in both time and space, is categorized into 14 different climate zones (Alizadeh-Choobari and Najafi, 2018). The primary aim of this study, as the first national-scale study, is to advance the understanding of how climate variability driven by ocean-atmosphere circulations impacts long-term groundwater storage (16-52 years) across Iran by analyzing groundwater level records from 24 principal aquifers. The influence that data length can have on the teleconnections between climate indices and groundwater levels is further assessed.

Study area and its teleconnection with ocean-atmosphere circulations
Iran with ~ 82 million inhabitants and 1,648,195 km 2 in the area is the second-largest country in the Middle East. Iran has been classified into 14 climate zones by Alizadeh-Choobari and Najafi (2018) based on principal component analysis with a varimax rotation that was accomplished through the 21 meteorological variables such as mean annual precipitation, temperature, and humidity ( Fig. 1). Table 1 lists the winter (DJF), summer (JJA), and annual mean surface temperatures, annual mean precipitation, and annual mean relative humidity of 14 different climate zones identified across Iran. The country's climate varies from mild and humid in the southern coasts of the Caspian Sea with a mean annual precipitation of > 1000 mm to warm and arid in the central to east parts with a mean annual precipitation of < 125 mm. The northwest portions are characterized by a cold and semi-arid climate with annual mean precipitation and temperature of 314 mm and 11.4 °C, respectively. Western Iran that contains climate zones of 2 (490 mm) and 6 (409 mm) receives the secondhighest annual precipitation of the country. The northeastern portions of the nation have a cold and semi-arid climate with annual mean precipitation and temperature of 280 mm and 13.9 °C, respectively. Southern Iran (climate zones of 11, 13, and 14) is characterized by a marine hot and (semi-)arid climate with a mean annual humidity of 54-75% and an annual mean rainfall of 121-243 mm. Figure 1 shows the selected aquifers mapped on the 14 climate zones across the country. The full description of these 14 climate zones is presented by Alizadeh-Choobari and Najafi (2018 in Supporting Information (SI)). This highly variable climate results from the climate impacts imposed by the geographical position of sea/gulf and land in the continent and the mountain ranges (Alborz and Zagros) that border two extremely arid deserts in central Iran (Kavir and Lut). The Zagros and Alborz act as a low-level air mass barrier for the flow of moisture that comes from west to southwestern Iran and the Caspian Sea towards central Iran (Alijani, 2000). Based on the length and continuity of representative water-level records, we selected 24 different principal 1 3 aquifers ( Fig. 1 and Table 2), at least one from each climate zone, to scrutinize their links with large-scale oceanatmospheric circulations. The selected aquifers have various mean elevations from ~ 2 (aquifer 3, at the southern coast of the Caspian Sea) to 2268 m.a.s.l (aquifer 6b, in the Zagros Mountains). The mean water table depth of the aquifers lies in a highly variable range from less than 10 m (aquifers 3 (2.8 m), 2 (3 m), 4 (6.3 m), 5 (8.3 m), and 13 (7.2 m)) to more than 50 m (aquifers 9c (68 m), 7c (58.5 m), and 9b (52.7 m)). Overall, the coasts of the Caspian Sea (with the highest precipitation) and the Lake Urmia catchment encompass the shallowest aquifers. In contrast, the deepest ones belong to central Iran, with the lowest rainfall.
Iran's overall precipitation increases during El Niño and decreases over the La Niña phase. Alizadeh-Choobari and Najafi (2018) pointed out that the El Niño has the highest impacts on the northern coasts of the Persian Gulf, with a warm and arid climate, and the least on southern coats of the Caspian Sea, with a mild and humid climate. On the contrary, the decrease rate in annual precipitation associated with La Niña is particularly significant across the marine warm and arid climate zone to the northern coast of the Oman Sea (Table 1). Rezaei (2021) likewise observed that the ENSO and PDO have the relatively lowest influence on southeastern Iran and the Caspian Sea coasts. In addition to ENSO (mean significant coherence of 0.57), other indices of PDO (0.49), NAO (0.43), and AMO (0.36) were found to be coherently linked to the integrated meteorological and agricultural drought (i.e., precipitation and soil moisture) across Iran (Rezaei, 2021). He also argued that drought across the most upper latitudes of northwestern Iran shows the highest coherence with the Atlantic-based indices (i.e., NAO and AMO). The selected aquifers mapped on the 14 climate zones of Iran presented by Alizadeh-Choobari and Najafi (2018). The black labels refer to climate zones while the blue ones to aquifers Based on previous works on exploring teleconnections between large-scale climate indices and Iran's climate, the indices of ENSO, PDO, NAO, and AMO have been identified as the most important drivers modulating hydroclimate variability in the country (e.g., Nazemosadat and Ghasemi, 2004;Salahi et al. 2005;Raziei et al. 2009;   Dezfuli et al. 2010;Nikzad et al. 2013;Darand, 2014;Nazemosadat et al. 2015;Mohammadrezaei et al. 2020;Nouri and Homaee, 2020;Rezaei, 2021). Much of the country is wetter (drier) than normal conditions during the El Niño (La Niña) phase of ENSO (e.g., Nazemosadat and Ghasemi 2004;Raziei et al. 2009;Dezfuli et al. 2010;Nazemosadat et al. 2015;Nouri and Homaee, 2020;Rezaei, 2021). In addition to the precipitation, the ocean-atmospheric oscillations can also modulate variables of temperature, runoff, and spring flowrate across the nation (e.g., Nazemosadat and Ghasemi, 2004;Nazemosadat et al. 2015;Abbasi and Maleki, 2017;Alizadeh-Choobari and Najafi, 2018;Ahmadi et al. 2019;AlizadehChoobari and Adibi, 2019;Dehghani et al. 2020;Rezaei and Gurdak, 2020). These works, on the whole, reported that (1) Iran is commonly wetter during the El Niño and drier during the La Niña phase of ENSO; (2) precipitation across the country tends to be above normal (below normal) during the negative (positive) phase of NAO; (3) the above (below) normal precipitation generally coincides with the positive (negative) phase of PDO, typically in winter; and (4) precipitation and soil moisture, on the whole, are correlated positively with Pacific-based circulations (ENSO and PDO) and negatively with Atlantic-based indices (AMO and NAO). Recently, Rezaei and Gurdak (2020) and Rezaei (2021) argued that Pacific-based climate indices have a higher impact than Atlantic-based indices on the variability in climate and water storage across Iran. Rezaei (2021) demonstrated that three-coupled indices of ENSO, PDO, and NAO have stronger coherence with Iran's drought than either two-coupled or single indices.

Methodology
We use the monthly time series from (1) ocean-atmosphere indices of AMO, PDO, NAO, and ENSO ( Fig. S2), and (2) the monthly representative groundwater levels from 24 aquifers listed in Table 2. Here, the Southern Oscillation Index (SOI), as the oldest indicator, represents the ENSO index. SOI's positive and negative values are associated with La Niña and El Niño, respectively. The monthly time series of the climate indices of AMO, PDO, NAO, and SOI are freely available online from National Oceanic and Atmospheric Administration (NOAA). Here, the principal unconfined aquifers that have the longest data record were chosen for each climate zone. The representative water table time series for each aquifer were taken from Iran Water Resources Management Company, Tehran. It is necessary to acknowledge that some of the variability in the groundwater level records could be related to human activities and not solely a function of natural climate patterns, owing in part to dry season pumping from aquifers (Hanson et al. 2004;Rezaei and Gurdak, 2020). The simple flowchart of the methodology is displayed in Fig. 2. Prior to analysis, a computer program of the USGS Hydrologic and Climate Analysis Toolkit, HydroCLimATe (Dickinson et al. 2014) was used to remove the long-term trend from each groundwater level series by subtracting a regression-fitted low-order polynomial. The detrending can partially eliminate (1) the anthropogenic signals (such as groundwater extraction and other water resource development), (2) the red noise from the original records to capture the temporal structure of the data series successfully (Kuss and Gurdak, 2014), and (3) much of the long-term multidecadal anthropogenic signals (e.g., groundwater extraction and crop irrigation) from the groundwater level series (Hanson et al. 2004;Gurdak et al. 2007). Although the pumping from groundwater is usually dominant during the dry season of each year and has a short-term effect, it can also significantly influence low-frequency signals in aquifer storage. In practice, human-induced over-extraction from aquifers in each dry season does not fully compensate through the natural recharge that occurred in the following wet season (Panda et al. 2007), so that, the anthropogenic stress on aquifers accumulates year by year. According to the literature review performed by Mohammed et al. (2018), anthropogenic activities such as ongoing humaninduced groundwater over-extraction act as one of the main drivers of the long-term trend of groundwater levels across the world. This study further focuses on interannual to multidecadal signals than seasonal variability. Similarly, the detrending procedure also has been extensively used by previous researchers to remove the anthropogenic impacts from groundwater level time series (e.g., Barco et al. 2010;Hanson et al. 2004;Kuss and Gurdak, 2014;Velasco et al. 2015;Rezaei and Gurdak, 2020). The wavelet analysis conducted here has four methodological steps, including the local continuous wavelet transform (CWT), the cross-wavelet transform (XWT), the wavelet coherence (WTC), and multiple wavelet coherence (MWTC). The local CWT ( W X n (s) for a time series ( x n , n = 1,..., N) with uniform time steps t has been defined as the convolution of x n with the scaled and normalized wavelet (Grinsted et al. 2004): where 0 is the Morlet wavelet, (*) the complex conjugate, s the wavelet scale, 0 dimensionless frequency, and dimensionless time. The significance of CWT is estimated using 1000 randomly constructed synthetic series generated by Monte Carlo methods (Grinsted et al. 2004). Furthermore, the global power spectrum, GPS (as a function of time), was also computed for each time series using a MATLAB script written by Schulte (2019). The degree to which the power spectrum exceeds the background noise was assessed by cumulative arcwise significant test (95% confidence bound). The cross wavelet transform (XWT) for each pair of time series (i.e., x n and y n ) was defined as W XY = W X W Y * , where * is complex conjugation. Finally, WTC for each pair of time series was calculated as follows (Torrence and Webster, 1999;Grinsted et al. 2004): where S(W) is a smoothing operator, S scale smoothing along the wavelet scale axis, and S time smoothing in time. The CWT was employed to define the dominant oscillation modes of each individual series (Grinsted et al. 2004). Next, the XWT was used for each time series to depict the cross-wavelet power of a pair series against the background power spectra (Torrence and Compo, 1998;Grinsted et al. 2004). A cone of influence (COI) characterizes the zone of the wavelet spectrum in which these edge impacts need to be excluded (Torrence and Compo, 1998;Grinsted et al. 2004). The cross-correlation (ranges from 0 to 1) as a function of frequency was computed between a pair series using the WTC. Finally, 1000 randomly constructed synthetic series were generated by Monte Carlo methods to test the significance of the WTC (Grinsted et al. 2004). The multiple wavelet coherence (MWTC) at scale s and location can be formulated as (Hu and Si, 2016): where ⃖�� ⃗ W Y,X (s, ) is the smoothed cross-wavelet power spectrum between multiple predictor variables X and response variables Y, ⃖�� ⃗ W X,X (s, ) signifies either the smoothed auto-wavelet power spectra for the main diagonal of the matrix (i.e., i = j ) or cross-wavelet power spectra when . The MWTC at the 95% confidence level was also computed by the Monte Carlo method (Grinsted et al. 2004). Figure 3 shows the CWT and GPS results for the climate indices (AMO, PDO, NAO, and ENSO) and the groundwater level series. AMO has significant power regions at three frequency bands of annual (particularly after 1998 in the local CWT), decadal (8-10 years, particularly after 1980 in both the local CWT and GPS), and interdecadal (> 16 years as a peak signal in GPS). In both the CWT and GPS graphs, PDO also has significant power regions across the three frequency bounds of interannual (4-6 years), decadal (8-12 years), and, to lesser magnitude, interdecadal (> 16-years). Significant power regions for NAO are observed at annual, interannual (6-10 years), and, to smaller magnitude, interdecadal (> 16 years) modes. For ENSO, there are two significant power patterns at interannual (2-7 years) and decadal (10-13 years) frequency bands. Notably, the power patterns from the Pacific-based indices (PDO and ENSO) are stronger than those from the Atlantic-based ones (AMO and NAO).

Continuous wavelet transform (CWT)
The local CWT and GPS analyses on the representative groundwater levels from the 24 aquifers ( Fig. 3) indicate that significant natural variability exists at annual (1 year), interannual (2-7 years), and decadal-interdecadal (8-18 years) frequency bands, consistent with climate indices. It seems that annual and interdecadal oscillation signals are the The local continuous wavelet transform (local CWT) and its associated global power spectrum (GPS) for both the climate indices (AMO, PDO, NAO, and SOI (i.e., ENSO)) and the aquifers' water levels. The thick-black and dotted lines on the local CWT and GPS are the 95% confidence level. For brevity, only one aquifer from each aquifer group is shown ◂ dominant modes in the majority of the aquifers except for aquifers 11-14 with short data length that we could not compute the interdecadal periodicities for them. Likewise, the dominant power regions in the integrated precipitation and soil moisture time series over Iran have occurred from the interannual to interdecadal periodicity (Rezaei, 2021). For the groundwater in the Lake Urmia catchment, Rezaei and Gurdak (2020) reported that the interannual oscillations are more related to ENSO and lesser to NAO while the decadal-interdecadal modes are typically originating from PDO. They removed the signals with periodicities of less than 2 years, leading to a significant loss of the annual signals. In contrast, the local CWTs here demonstrate that the annual signals from most cases emerge over the whole transect, signifying the importance of annual cycles. The commonality in results of CWT analyses from groundwater here and both precipitation and surface soil moisture by Rezaei (2021) suggests that these hydroclimate variables are responding to a common coherent signal, being driven by the same mechanism or drivers (i.e., large-scale oceanicatmospheric teleconnections). This further confirms by broadly similar patterns existing between the interannual signals from hydroclimate variables and dominant frequencies from ENSO, PDO, and NAO as well as by similarity between the interdecadal signals from hydroclimate variables and dominant modes in PDO and AMO.

Wavelet coherence (WTC)
This study examines the links between groundwater level and climate indices of AMO, PDO, NAO, and ENSO by computing WTC to visually illustrate both those frequency bands and time intervals at which climate indices and groundwater level series are covarying (Torrence and Webster, 1999). Significant coherence is depicted in each run of the WTC analysis in which contours enclose statistically significant periods (P < 0.05), based on a red-noise process as determined by a Monte Carlo experiment (Grinsted et al. 2004). Figure S3 displays the 96 WTC's runs computed for teleconnections between 24 aquifers and four climate indices of AMO, PDO, NAO, and ENSO. For brevity, the results for all 14 aquifer groups are elaborated in section "WTC analyses in more details" in SI and here we only present the general findings. Consistent with the dominant modes obtained from the CWT modes (Fig. 3), the significant coherence patterns between the climate indices and groundwater levels likely occur at five frequency bands of the annual (⁓1 year), short-interannual (2-4 years), mediuminterannual (4-6 years), decadal (8-12 years), and interdecadal (14-18 years). Hereafter, for brevity, we use these terms for dominant modes without referring to the associated year bands. To statistically assess and easily compare the results, the coherence values from all the 24 aquifers are averaged across each of the significant frequency bands (Fig. 4).
The annual variability in the groundwater levels from all the 24 aquifers is coherently linked to AMO (mean coherence value of 0.58) and, to lesser magnitude, PDO (0.47) and NAO (0.45). ENSO (0.32) shows the lowest coherence with annual frequency signals (Fig. 4a). The mean coherence value between AMO and water level across the annual mode is computed to be ≥ 0.70 for aquifers 8b (0.80), 2 (0.80), 6c (0.79), 7a (0.77), 1c (0.76), 6a (0.73), 13 (0.72), 7c (0.72), 6b (0.71), and 1b (0.70), suggesting the critical importance of the annual variability. The AMO's coherence with annual signals of the water levels is highly stronger from 1998 to 2019 compared to the latest 30 years of the twentieth century. As a case example, Fig. 5 illustrates the strong coherence patterns across annual bands of the WTCs for aquifers 1c, 6c, and 8b.
For interannual frequency bands, the significant coherence patterns have occurred at 4-6-year (in most cases) and 2-4-year (in some cases) periods (Fig. S3). For all aquifers, the mean coherence value at the 2-4-year frequency bands for ENSO (i.e., SOI), AMO, PDO, and NAO equals to 0.38, 0.35, 0.33, and 0.33, respectively (Fig. 4). The significant coherence patterns across the 2-4-year modes are more evident in the groundwater levels of the Caspian Sea's southern coast (i.e., aquifers 3 and 4), western Zagros (aquifers 6a, 6b, and 6c), and northern coasts of the Persian Gulf (aquifers 13 and 8b) and the Oman Sea (aquifer 14) (Fig. S3). Figure 4a shows the mean values of coherence over all the aquifers from PDO (the mean coherence value of 0.44), ENSO (0.43), NAO (0.43), and AMO (0.40). Although the coherence means for both the PDO and ENSO are close to each other, based on the extent and significant level of the observed patterns in the WTCs (Fig. S3), the interannual 4-6-year signals of groundwater levels, in decreasing order, are likely controlled by PDO, ENSO, NAO, and AMO. PDO and ENSO usually have the same coherence patterns at the interannual 4-6-year periods in the WTCs (Fig. S3), particularly in groundwater levels from northwestern Iran (aquifers 1b, 1c, 3, and 7b), western Zagros (aquifers 6a, 6b, 6c, and 8a), and northeastern Iran (aquifers 7a and 7c). Consistent with the significant interannual power regions in the CWT for ENSO (Fig. 3), the coherence patterns of ENSO and PDO at the WTCs are concentrated at the time-spans before 2000 and after 2010; the former is more apparent in the majority of aquifers such as 1a, 1b, 2, 4, 6b, 7a, 7c, 8a, 9c, and 12 and the latter is stronger in the aquifers of 1b, 1c, 2, 3, 6c, 7b, 7c, 8b, 9c, 11, 12, and 14 ( Fig. S3). At the interannual 4-6year signals, NAO has the highest coherence with aquifers of 6b, 7a, 7b, 8a, and 9c (Fig. S3); among them, aquifers 7a and 7b belong to upper latitudes to northeastern Iran. Rezaei (2021) likewise reported that NAO has the highest effect on the precipitation and soil moisture at the upper latitudes of the country. AMO relatively has higher coherence with water levels of the northern (aquifers 1b, 1c, 6c, 7a, and 7c) and the central (aquifers 9a and 9c) parts of the country.
Consistent with the significant power regions across the decadal signals shown in the CWTs (Fig. 3), the decadal signals from the groundwater levels shown in Fig. 4 are more correlated with the Pacific-based indices of ENSO (mean coherence value of 0.64) and PDO (0.63) rather than the Atlantic-based indices of NAO (0.41) and AMO (0.37).

Fig. 5
The wavelet coherence (WTC) between climate indices of AMO, PDO, NAO, and SOI and the representative groundwater levels from the aquifers 1, 6, and 8. The thick black curve is the 5% significance level, and the less intense colors reveal the cone of influ-ence. Black arrows are the phase angle, identifying the phase relation between two series for which right-and left-directed arrows representing in-phase and anti-phase relations, respectively (Grinsted et al. 2004) From the above descriptions, each dominant mode in groundwater levels is significantly affected by multiple indices more than a single climate index. Furthermore, Wyatt et al. (2012) found significant co-variances between all North Atlantic indices, and Rezaei (2021) revealed that the interacting effects of ENSO + PDO + NAO on integrated precipitation and soil moisture across Iran are likely larger than each of the individual indices. Therefore, we further examined the multiple interacting effects from some combinations (two-and three-coupled) of the most effective indices by computing the MWTC. Figure 4b and c shows the results from multiple interacting indices through analyzing the MWTCs. Here we select those two-coupled indices which allow to compare the integrated effects of the pacific-based (PDO and ENSO) with Atlantic-based (AMO and NAO) indices as well as to compare the indices of the dominant interannual periodicities (ENSO and NAO) with those of dominant interdecadal modes (AMO and PDO). For the three-coupled indices, Fig. 4c only displays the three most important multiple indices of AMO + NAO + PDO, AMO + PDO + SOI, and NAO + PDO + SOI, displaying the highest and broadest coherence patterns with groundwater levels. At all the dominant annual, decadal, and interannual periodicities, the multiple indices ( Fig. 4b and c) have a higher coherence with groundwater levels compared to single indices (Fig. 4a). Figures 4 and S3 show that the interdecadal signals are more correlated with the lower-frequency indices of AMO (mean coherence value of 0.71) and PDO (0.62) rather than the higher-frequency indices of ENSO (0.29) and NAO (0.29). Given the limit of data length, the interdecadal patterns are only captured in the three aquifers of 1a, 5, and 7b in the northern parts of the country. The significant coherence patterns between AMO and PDO indices and groundwater levels in both aquifers 1a and 7b (in northwestern Iran) are stronger than aquifer 5 (in the southeastern coasts of the Caspian Sea). Similarly, the Caspian Sea coast is among those regions in the country that show the lowest precipitation and soil moisture correlation with climate indices (Nazemosadat and Ghasemi, 2004;Rezaei, 2021), most probably attributed to the contribution of smaller-scale climate regimes triggered by the Caspian Sea and the Siberian high (located to the north of the Caspian Sea). In practice, the Siberian high modulates the intensity of easterly oversea winds (heading to the southwestern coasts, i.e., climate zone 3 in Fig. 1) during late autumn and winter (Nazemosadat and Ghasemi, 2004). El Niño phase is usually associated with a weaker surface pressure of the Siberian high and its associated 500-hPa trough and 200-hPa jet stream (Zhang et al. 1996). Although AMO and PDO have a stronger footprint in the interdecadal signals of the water levels, ENSO can be considered as a weaker counterpart since there appear moderate to strong coherence patterns between ENSO and water levels, particularly before 2000.

Teleconnections between climate indices and groundwater levels
Understanding the nature of annual to interdecadal variability in groundwater storage and their coupling and responding to climate patterns could help in improving our insight from aquifer response to climate variability and their potential susceptibility to change. To see in which part of the country groundwater is more affected by climate indices, in addition to the WTC analysis, the mean coherence values are computed over the five significant frequency bands (annual to interdecadal) from all the WTCs between climate indices and aquifers (Fig. 6).
The groundwater levels are found to be coherently linked to natural climate variability and not solely a function of temporal patterns in pumping since there is significant coherence between AMO, PDO, NAO, and ENSO indices and water levels in many aquifers. This partial control on groundwater storage from natural climate variability has also been reported from different parts of the world, such as the USA, Canada, and Europe (e.g., Hanson et al. 2004;Perez-Valdivia et al. 2012;Kuss and Gurdak, 2014;Velasco et al. 2015;Rust et al. 2019). The aquifers of the western Zagros (mean coherence value of 0.54 for group G6), northeastern (0.48 for group G1), northern (0.47 for group G3 and 0.45 for group G4), and northeastern (0.46 for group G7) regions of the country show the highest coherence with the ocean-atmospheric circulations. These regions, on the whole, have a shallower water table, higher precipitation, and lower temperature. On the contrary, the lowest teleconnections are observed in the aquifers Fig. 6 The climate indices' coherence values averaged over the five dominant frequency bounds from annual to interdecadal across the whole country of the southeastern (0.29 in group G12), and the central to eastern (from central Iran, G10 (0.34), to the northern coasts of the Persian Gulf G11 (0.40)) portions. Furthermore, Alizadeh-Choobari and Najafi (2018) reported that ENSO has the lowest correlation with precipitation in southeastern Iran. Rezaei (2021) likewise observed that the ENSO and PDO have the relatively lowest influence on southeastern Iran. Possible reasons for this are (1) northern Iran and western Zagros receive the highest annual precipitation across the country while the eastern and southeastern portions are the driest regions with a precipitation of < 150 mm/year, and (2) the shallowest aquifers almost belong to northern Iran and western Zagros while the aquifers across the eastern and southeastern Iran tend to have thicker unsaturated zones with a larger damping and lagging contributions. In fact, the arid climate may create strong upward total potential gradients, resulting in a decrease in downward water flux in the unsaturated zone (Walvoord et al. 2003;McMahon et al. 2006). Consequently, the linkage between atmospheric circulations and Iranian aquifers generally weakens from the upper (typically the northwest and western Zagros regions) to lower latitudes. Rezaei (2021) likewise argued the impact that Atlantic-based (i.e., NAO and AMO) indices have on soil moisture and precipitation over the nation, decreasing southward.
Of interest in this study is the observation that the teleconnection between the climate indices and each groundwater level is frequency specific. In other words, an individual index may have a significant coherence with a specific frequency band while it shows no strong linkage with other frequencies. In Iranian aquifers, there are five frequency bands at the WTCs characterized by significant coherence patterns including annual, short-interannual, medium-interannual, decadal, and interdecadal, consistent with the results from the CWTs (Fig. 3). In each significant frequency band, only one or two climate indices are largely effective. The annual signals largely correlate with AMO and, to lesser magnitude, PDO and NAO.  2018) reported that both the ENSO and AMO largely correlate with Iran's precipitation. It has also been observed that AMO has a significant effect on temperature across northwestern Iran (Abbasi and Maleki 2017; Rezaei and Gurdak, 2020). Nonetheless, we observed that ENSO has the lowest correlation with annual signals of groundwater levels while the considerable imprint of ENSO on hydroclimate variables (not groundwater levels) has been previously reported in the literature (e.g., Alizadeh-Choobari and Najafi, 2018; Ahmadi et al. 2019; Alizadeh-Choobari and Adibi, 2019). The possible reason for this is the damping and lagging effects of the unsaturated zone elaborated in the last part of this section.
For the interannual frequency bounds, the significant coherence patterns between climate indices and aquifer's water level have occurred at modes of short-interannual (in most cases) and medium-interannual (in some cases) (Fig. S3). At the short-interannual mode, although the WTCs for AMO, PDO, and NAO indices show relative coherence with groundwater levels, ENSO has a higher coherence with aquifers because ENSO has the strongest power spectrum at short-interannual modes (Fig. 3). At the short-interannual signals, both the PDO and ENSO have relatively stronger influences on groundwater levels compared to NAO and AMO since, beyond the larger mean coherence values (Fig. 4), they show the larger extent and stronger significant coherence patterns at the WTCs (Fig. S3). The stronger coherence patterns of the medium-interannual oscillations can be attributed to either the higher damping of climate signals through the unsaturated zone (Velasco et al. 2015;Rust et al. 2018) or the presence of stronger medium-interannual modes in both the PDO and ENSO compared to the shortinterannual signals (Fig. 3). In the short-interannual signals, both the PDO and ENSO usually have the same coherence patterns with water levels, particularly in northwestern Iran (aquifers 1b, 1c, 3, and 7b), western Zagros (aquifers 6a, 6b, 6c, and 8a), and northeastern Iran (aquifers 7a and 7c) (Fig. S3). Consistent with the significant interannual power spectrums in the CWT analyses for SOI (Fig. 3), the coherence patterns between SOI and groundwater are stronger before 2000 in most aquifers including 1a, 1b, 2, 4, 6b, 7a, 7c, 8a, 9c, and 12 and after 2010 in the aquifers of 1b, 1c, 2, 3, 6c, 7b, 7c, 8b, 9c, 11, 12, and 14. This broad similarity between ENSO and PDO's coherence patterns derives from the synergy between ENSO and PDO, discussed by Newman et al. (2003) and Juanxiong et al. (2004). This appears to be in good agreement with the integrated strong footprints of ENSO and PDO on Iran's drought (Rezaei, 2021) and on soil moisture and dust across the Fertile Crescent, an area that consists of Iraq, Syria, Jordan, Israel, and western Iran (Notaro et al. 2015). Furthermore, the significant correlation between ENSO and hydroclimate variables across Iran has been previously postulated (Alizadeh-Choobari and Najafi, 2018; Ahmadi et al. 2019; Alizadeh-Choobari and Adibi, 2019; Rezaei and Gurdak, 2020). Similarly, ENSO has also been proposed to modulate variations in river flow (Barlow and Tippett, 2008), soil moisture (Sheffield and Wood, 2008), and NDVI (Kariyeva and van Leeuwen, 2012) in the Middle East and Southwest Asia.
Consistent with the significant power spectrums at the decadal oscillations in the CWTs (Fig. 3), the decadal signals from the groundwater levels shown in Fig. 4 are more correlated with the Pacific-based indices of ENSO (0.64) and PDO (0.63) rather than NAO (0.41) and AMO (0.37).
Likewise, the highest coherence between ENSO and PDO indices and Iran's integrated precipitation and soil moisture occurs in two bands of 4-7 years and ~> 10 years (Rezaei 2021). Similar to the 4-6-year signals, both the ENSO and PDO also show similar significant coherence patterns at the decadal oscillation signals for almost aquifers such as 1a, 1b, 1c, 3, 4, 5, 6a, 6b, 7b, 7c, 8b, 9a, 9c, 10b, 11, and 13 (Fig. S3). These aquifers are roughly distributed at a broad portion of the country, indicating the decadal oscillations are one of the dominant signals in groundwater systems that are highly controlled by natural ocean-atmospheric circulations, particularly those from the Pacific Ocean, compatible with the results for the precipitation and soil moisture reported by Rezaei (2021). The CWT results (Fig. 3) also show that the decadal oscillations are the dominant signals observed in both the climate indices of ENSO and PDO as well as in most of the above aquifers. ENSO and PDO have stronger control on the decadal signals of the water levels than the Atlantic-based indices of AMO and NAO since the decadal oscillation, as a dominant mode, in the original ENSO and PDO time series has stronger power than those in AMO and NAO.
Both the mean coherence values (Fig. 4a) and the significant coherence patterns at the WTCs (Fig. S3) indicate that beyond the decadal modes, the interdecadal oscillations are also dominant modes in Iranian aquifers. The interdecadal signals in the water level measurements are reasonably more correlated with the lower-frequency indices of AMO (0.71) and PDO (0.62) compared to the higher-frequency indices of ENSO (0.29) and NAO (0.29) (Figs. 4 and S3), because of low-frequency signals tend to be preserved better than highfrequency signals in groundwater fluctuations (Velasco et al. 2015). Likewise, the higher influence of the lower-frequency index of PDO on groundwater systems is found for the US West Coast (Velasco et al. 2015).
Given the above evidence, each climate index tends to differently impact different oscillation modes of groundwater levels. AMO's significant effect is observed at interdecadal and annual modes while PDO's highest effect is related to decadal and interdecadal oscillations. The highest imprint is observed to be at the decadal and medium-interannual periodicities for ENSO and at annual and interdecadal frequency bands for NAO. Notably, the results from each index are in good agreement with its dominant power spectrums depicted in Fig. 3. As an example, the medium-interannual and decadal modes are the dominant signals in the local and global CWTs for ENSO (i.e., SOI). Overall, the decadal and interdecadal signals are the dominant modes in the climate indices, particularly PDO, AMO, and ENSO (Figs. 4 and S3), compatible with the dominant power regions in the CWT for the majority of aquifers (Fig. 3). It is found that PDO (with the mean coherence of 0.51) has the highest influence on groundwater fluctuations over Iran compared to the other three indices of AMO (0.45), ENSO (0.42), and NAO (0.34). This well agrees with those results obtained by Rezaei and Gurdak (2020) for the aquifers in the Lake Urmia watershed.
Despite the similarity between the teleconnections observed in the majority of aquifers in this study and those in the surface hydroclimate (precipitation, drought, and rivers) in the literature, ENSO has the lowest effect on the annual oscillation modes from the groundwater levels while most of the previous results found that ENSO is one of the most important indices modulating Iran hydroclimate variability (but not groundwater) (e.g., Alizadeh-Choobari and Najafi, 2018;Ahmadi et al. 2019;Alizadeh-Choobari and Adibi, 2019;Rezaei and Gurdak, 2020). We further observed that the groundwater variability, particularly those with lower frequency, is more correlated with lower frequent climate circulations such as PDO and AMO, rather than ENSO. These observations most probably arise from the damping and lagging role that the unsaturated zone has on the climate variability signals. In fact, the damping rate of the unsaturated zone on higher-frequency signals of ENSO (i.e., SOI) and NAO is greater than that on the lower-frequency indices (Bloomfield & Marchant, 2013;Van Loon, 2015;Van Loon et al. 2014;Rust et al. 2018). While this phenomenon is absent in other hydroclimate variables (e.g., precipitation, temperature, surface soil moisture, and drought), so that, the aquifer teleconnections are relatively different from those observed in precipitation, drought, and temperature (Velasco et al. 2015). To clarify this issue in our case, Fig. 7 shows the coherence values averaged across the 1-, 2-4-, 4-6-, and 8-12-year frequency bands versus mean water table depth for each aquifer. The coherence values are generally smaller in the deeper aquifers (i.e., regression lines have a negative slope), particularly for the annual signals which have a higher frequency (Fig. 7a). Furthermore, the highest coherence patterns between AMO and groundwater levels across 1-year periodicity have occurred at the shallowest aquifers of 2, 4, 8b, and 5 with water table depths of < 10 m. In contrast, the lowest coherence between AMO and annual oscillation modes is observed in aquifers 14, 9a, 10a, and 11 to the central and eastern parts of the country with deeper aquifers and drier climates (Fig. S1). In fact, the thicker unsaturated (i.e., deeper water table) zones, particularly those of lower hydraulic conductivities, have a stronger damping and longer lagging effect on the climate variability signals (Rust et al. 2018). Likewise, it is found that shallower aquifers receive greater amplitude of NAO signals across Europe (Rust et al. 2018). Furthermore, as the frequency of the signals increases, the slope of the fitted lines also tends to increase, that is, the slope of the 1-year oscillation modes (Fig. 7a) is greater than that of the lower frequency signals of interannual and decadal modes (Fig. 7b-d). We therefore concluded that the higher-frequency signals are damped more than the lower-frequency modes across the Iranian aquifers.

Multiple interacting controls of climate indices
Comparing Figs. 4a, b, c; S3; and S4 indicates that the coupled indices have stronger and broader significant coherence patterns than each individual index with the groundwater levels, suggesting that multiple interacting controls of indices can better predict the water table behavior across the country. As an example, comparison of the WTCs and MWTCs results respectively from each individual index (Fig. S3) and three-coupled (Fig. S4) indices for each particular aquifer shows that coherence patterns in Fig. S4 are largely stronger and broader than those in Fig. S3, where the colors tend to be more in yellow (higher coherence values) in Fig. S4. Figure 4a,b, and c also reveal that the threecoupled indices, on the whole, have larger mean coherence values than both the two-coupled indices and each individual index across the country. However, among the twocoupled indices, the PDO + SOI index has the highest effect on groundwater levels at the short-interannual, decadal, and interdecadal oscillation modes while the AMO + PDO index tends to highly modulate the annual and mediuminterannual frequency signals (Fig. 4b). Figure 4b also shows that the Pacific-based indices (PDO + SOI with mean coherence value of 0.74) have a stronger effect on Iran's aquifers than the Atlantic-based (AMO + NAO with mean coherence value of 0.70), consistent with the findings made from the integrated precipitation and soil moisture across Iran (Rezaei, 2021). Among the three-coupled ones, except for annual signals in groundwater levels which are correlated with AMO + NAO + PDO, all the other dominant signals from interannual to interdecadal are further correlated to the three-coupled AMO + PDO + SOI index (Fig. 4c). Consequently, the combination of AMO + PDO + SOI is highly capable of explaining the natural groundwater variability in Iran. Wang et al. (2014) also stated that the La Niña coupled with cold PDO highly affects the severe droughts across the Fertile Crescent owing to anomalous ridging across Eastern Europe and northern Asia and subsidence over Fertile Crescent itself.
The multiple interacting impacts from the three-coupled AMO, PDO, and ENSO are mostly attributed to the fact that the ocean-atmosphere circulations usually arise from physically interconnected mechanisms. ENSO, in addition to the north Pacific (i.e., PDO) and tropics, also affects the climate in the North Atlantic (i.e., AMO and NAO) (Bronnimann, 2007). Remind that AMO computes from the sea surface temperature measured across the northern Atlantic. The PDO is dependent upon ENSO on all timescales, so that, it can be taken into account as the reddened response to both atmospheric noise and ENSO . Newman et al. (2016) further argued that PDO is not a single independent event, but rather arises from physical processes through which the ENSO cycle interacts with extratropical variability to derive similar PDO-like SST anomaly patterns. ENSO, on the whole, modulates the PDO throughout the year. As an example, during boreal winter overlapping with El Niño events peak, the Aleutian low deepens and the variations in the surface heat fluxes, wind-driven mixing, and Ekman transport in the upper ocean all interact, leading to a positive PDO phase (Alexander, 1990). Various observations and models have diagnosed that El Niño (La Niña) events are usually associated with negative (positive) NAO-like atmospheric anomaly patterns, especially in winter times (e.g., Mathieu et al. 2004;Bronnimann, 2007), possibly because the North Pacific's atmosphere provides a pathway for ENSO-related diabatic heating in the tropical Pacific to modulate the atmospheric circulation patterns across the North Atlantic (e.g., Graf and Zanchettin 2012;Zhang et al. 2019a). Scaife (2010) argued that vertically propagating Rossby waves, which weaken the stratospheric polar night jet, provides the teleconnection pathway from the tropical Pacific (i.e., ENSO) to the Northern Pacific, and into the stratosphere. Descending wave mean flow interaction then results in a few weeks delay while the ENSO anomaly spreads to the lower stratosphere. Finally, on reaching the tropopause, there is a quick response, likely happening through baroclinic eddies in the storm tracks which equates the negative NAO response observed during El Niño phases across the Atlantic-European region (Scaife, 2010). Zhang et al. (2019b) further argued when ENSO and AMO are in phase, the negative ENSO-NAO correlation is substantial during late winter since the sea surface temperature (SST) anomaly across the tropical North Atlantic responses to both ENSO and AMO. They reported that a warm (cold) SST anomaly in the tropical North Atlantic emerges when El Niño (La Niña) is synchronous with a positive (negative) AMO phase, leading to a substantially weakened NAO. Furthermore, Levine et al. (2017) found that variations in AMO-related tropical Atlantic SSTs lead to variations in the Walker circulation across the tropical Pacific Ocean, modifying ENSO stability on both annual and multidecadal time scales.

Temporal changes in teleconnections
Since the ocean-atmospheric circulations, human demand for groundwater, and groundwater-dependent ecosystems are temporarily variable (Klove et al. 2014), the teleconnection patterns shown in the WTCs are not the same at different times. As a case example, Fig. 5 exhibits the coherence between AMO and groundwater levels of aquifers 1c, 6b, and 8b that intensified after 2003. Likewise, Rust et al. (2018) pointed out that the low-frequency signals in the groundwater levels have both spatial and temporal variability. We therefore scrutinize the coherence variations over time by averaging the coherence value across all the frequency bands to obtain coherence means at different times (Fig. 8). In Fig. 8, the mean standard deviation for ENSO, PDO, AMO, and NAO is estimated to be 0.09 (0.01-0.15), 0.10 (0.02-0.14), 0.10 (0.01-0.15), and 0.07 (0.03-0.12), respectively. Some important observations that can be made from Fig. 8 are as follows: First, although the Pacific-based indices of PDO and ENSO, on the whole, have a higher correlation with Iran's groundwater compared to Atlantic-based indices of AMO and NAO, their teleconnections have significantly weakened after 2000. PDO and ENSO's lowest coherence values are observed from 2003 to 2012, particularly in aquifers 2, 3, 7c, 8a, 8b, 9a, 10b, and 12 as shown in Fig. S3. Likewise, Rajagopalan et al. (1997) argued the influence that ENSO has on Middle Eastern precipitation which intensified during Fig. 8 The mean coherence values between the indices and aquifers' water level series in the different time bands across the country. The mean standard deviation for SOI, PDO, AMO, and NAO is estimated to be 0.09 (0.01-0.15), 0.10 (0.02-0.14), 0.10 (0.01-0.15), and 0.07 (0.03-0.12), respectively. The orange boxes highlighted the severe historical 1998-2000, 2007-2009, and 2010-2011 droughts across Iran the two last decades of the twentieth century partly inferred from variations in the frequency and intensity of ENSO occurrences since the mid-1970s. Although the teleconnections tend to gradually strengthen after 2010 again, it is still weaker than that observed before 2000 for PDO and ENSO indices. Notably, a relative decrease in coherence between indices and water table over the period 1976-1980 is also evident in Figs. 8 and S3 (especially, in the interannual and decadal oscillation signals at the WTCs for aquifers 1a and 7b). A weakness in coherence observed over [2003][2004][2005][2006][2007][2008][2009][2010][2011][2012] leaves open the possibility of either the effect of natural fluctuations of aquifers or the thickening of the unsaturated zone during the two last recent decades corresponding with groundwater over-extraction, severe water table dropping, and more frequent occurrences of severe drought events across Iran (Rezaei and Mohammadi, 2017;Ahmadi et al. 2018), or the global warming effect on teleconnections since, as an example, some researchers argued the wintertime atmospheric teleconnection of ENSO in the North Pacific tends to shift northeastward (e.g., Michel et al. 2020). However, this issue deserves further investigations that are out of the scope of this paper.
Second, unlike the Pacific-based indices, relative effects of the Atlantic-based ones (AMO and, to lesser magnitude, NAO) tend to strengthen on Iran's groundwater after 2000. The impact that AMO has on groundwater storage reached its maximum during the period from 2010 to 2019 (Fig. 8). Figure S3 clearly shows that this reinforcement in AMO's coherence with water levels is more apparent in the annual frequency bands at the WTCs of the majority of aquifers such as 1a, 1b, 1c, 2, 3, 4, 5, 6b, and 6c. Figure S2 demonstrates that unlike the latest 30 years of the twentieth century, AMO turned to its positive (warm) phase after 2000 and remains near the constant value on average, signifying that the positive AMO has possibly stronger effects on Iranian aquifers. Ahmadi et al. (2018) pointed out that the drought's frequency and the temperature over Iran have increased since 2000, overlapping with the start of the warm (positive) phase of AMO.
Third, it seems that the Iranian aquifers respond to droughts sooner than the wet periods. Rezaei (2021Rezaei ( ) determined that 1998Rezaei ( -2000Rezaei ( , 2007Rezaei ( -2009Rezaei ( , and 2011Rezaei ( -2012 Fig. 8) were the most severe droughts across the country since 1979 that have typically occurred in association with the strong La Niña phases. Interestingly, the SOI (i.e., ENSO) curve in Fig. 8 shows relative peaks coincided with these severe drought years, displaying quick aquifers' response to droughts in comparison with the wet periods. As an example, during the severe 1998-2000, 2007-2009, and 2010-11 droughts that are associated with the historical strong La Niña in 2007-2008(Rezaei, 2021, moderate to strong coherence patterns around the 1-2-year periods at WTCs are observed between ENSO and water levels for aquifers of 1a, 1c, 3, 4, 5, 6a, 6b, 6c, 7a, 7b, 8b, 9a, 9c, 10a, 11, 13, and 14. There appears a significant power region in the CWTs of ENSO across ~ 2-year periods from 2007-2012, confirming the strong coherence patterns in those aquifers. This relatively quick aquifer's response to droughts is mostly attributed to the over-extraction from groundwater in the following dry season exacerbating the drought effect. Notably, pumping from groundwater is the primary annual discharge component for Iranian aquifers (Rezaei and Mohammadi, 2017;Mirzaei et al. 2019). Furthermore, unlike the recharge phenomenon, which is significantly affected by the unsaturated zone, the discharge by pumping circumvents the unsaturated zone roles in damping and lagging the discharge signals. The recharge dominantly occurs through the unsaturated zone by natural phenomena. In contrast, the discharge from Iranian aquifers is dominantly performed by pumping without the significant effect of the unsaturated zone. Notably, detrending the groundwater level series by subtracting a regression-fitted low-order polynomial, explained in the "Methodology" section, may fail to remove the seasonal/annual anthropogenic signals completely. Otherwise, the recharge shock during the wet years may percolate to water table with longer time delay and more damping compared to the drought periods, particularly in the deeper aquifers with higher potential in signal attenuation (Bloomfield & Marchant, 2013;Van Loon, 2015;Van Loon et al. 2014;Rust et al. 2018).

Data length effect on teleconnections
The data length can significantly affect the total coherence patterns and our eventual assessment of teleconnections since the above evidence implies that the climate index teleconnections with groundwater series undergo changes over time. The findings show that the data length can significantly affect the results if the time series are not contemporaneous and only one value of coherence/correlation is obtained from a particular series instead of separate computation for different frequency bands and different time spans. As an example, Fig. 9 shows the relationship between the coherence means and data length for the four aquifer groups G1 to G4 (from northwest Iran to southern coasts of the Caspian Sea) and G8 to G11 (from east-center to the northern coasts of the Persian Gulf). Here, adjacent aquifers are used for plotting the curves of data length versus coherence means, owing in part to similar teleconnections of adjacent aquifers. The regression-fitted lines with R 2 = 0.95 for G1 to G4 and R 2 = 0.80 for G8 to G11 (Fig. 9) demonstrate that the coherence means are substantially affected by data length: the larger values obtained from longer data lengths. Consequently, those teleconnections that are only assessed by the results from either coherence means or correlation values (i.e., ignoring the time specific) may not be reasonable except for the contemporaneous time series. Therefore, wavelet analyses performed in this study are reasonable since different coherence values are computed across the different frequency and time bands for each water level series. Finally, this work suggests that the wavelet analysis can assess the teleconnections better than the simple correlation analysis, particularly when the data lengths from different locations/stations are not contemporaneous. The correlations computed over the whole time span would not provide effective information (Yuan et al. 2016).

Physical mechanism behind the teleconnections across Iran
To physically understand the teleconnections between NIÑO, PDO, AMO, and NAO indices and precipitation, and, in turn, the groundwater levels across Iran, it is to be noted that precipitation over the Middle East is linked to the sea surface temperature (SST) in the three particular regions of the equatorial Pacific, North Atlantic, and equatorial Indian Ocean (Niranjan Kumar et al. 2016). The equatorial and Atlantic SSTs have a strong connection through the extension of the Pacific ENSO signal into the Atlantic (e.g., Lee et al. 2013). Niranjan Kumar et al. (2016) argued that the remote impact of ENSO on the Middle East and Iran precipitation is to be through the Rossby wave trains from the equatorial Pacific SSTs. In practice, the phases (positive or negative) of the Rossby wave pattern vary during El Niño and La Niña events which additionally lead to a displacement in the Middle East jet stream positions. During the El Niño (La Niña) phase, the jet streams across the Middle East shift equatorward (poleward), leading to strengthening (weakening) of the jet stream that could bring stronger (weaker) winter storms from the Mediterranean region to Iran (Niranjan Kumar et al. 2016). Alizadeh-Choobari et al.
(2018) also pointed out that the southward shift in the subtropical jet stream over Southwest Asia associated with the El Niño further pushes the Mediterranean storm track equatorward and increases precipitation across Iran, while the intensified jet during La Niña weakens the quasistationary mid-tropospheric planetary waves, resulting in a dry condition in Iran. Remarkably, a positive (negative) PDO leads to strengthening (weakening) and weakening (strengthening) the El Niño and La Niña impacts on Iran's climate, respectively. In practice, both jet stream and storm track (i.e., moist air) tend to shift equatorward (poleward) during positive (negative) PDO, typically in winter times (Sung et al. 2014;Wang et al. 2014), contributing to wetter (drier) conditions in lower latitudes such as Iran. Moreover, Wang et al. (2014) pointed out that the El Niño (La Niña) overlapped with positive PDO deepened (weakened) trough of low pressure that accelerates (impedes) the Middle East jet stream, providing a synoptic condition to increase (decrease) precipitation and the associated recharge rate into aquifers across Iran. Although NAO's imprint on Iran's precipitation is not strong as much as ENSO and PDO (Salahi et al. 2005;Rezaei, 2021), the Atlantic storm track during its positive phase tends to be northward and reinforcing precipitation over western and northwestern Europe and deriving a drier climate across the Middle East (Cullen et al. 2002;Donat et al. 2014). Therefore, higher historical precipitations and, in turn, larger percolation rates to groundwater over Iran are usually mostly accompanied by a strong El Niño overlapped with warm PDO and negative NAO.

Conclusions
The main findings of this study are as follows: -The significant coherence patterns between the climate indices and groundwater levels have occurred at five frequency bounds of the annual (~ 1 year), short-interannual (2-4 years), medium-interannual (4-6 years), decadal (8-12 years), and interdecadal (14-18 years), consistent with the dominant modes obtained from the CWT analyses.
-The annual variability largely correlates to AMO (mean coherence value of 0.58) and, to lesser magnitude, PDO (0.47) and NAO (0.45). Both the interannual Fig. 9 The mean total coherence value versus the mean data length for aquifer groups a G1 to G4 and b G8 to G9 and decadal oscillation signals in the groundwater levels are coherently linked to the Pacific-based indices of PDO and ENSO rather than the Atlantic-based indices of AMO and NAO. The interdecadal modes are more affected by the lower-frequency indices of AMO (0.71) and PDO (0.62) rather than the higher-frequency indices of ENSO (0.29) and NAO (0.29).
-The ocean-atmospheric circulations have the highest footprints on the Zagros and upper latitude regions of the country characterized by the shallowest water table, highest precipitation, and lowest temperature. The lowest teleconnections are rather observed in the aquifers of the southeastern, central, and eastern parts with thicker unsaturated zones and drier conditions. -The influence of each climate index varies much over different modes of groundwater levels. AMO's significant effect is observed at interdecadal and annual modes while PDO's highest effect is related to decadal and interdecadal. The highest influence of ENSO is observed across the decadal and interannual whereas the NAO's highest coherence patterns are related to annual and interdecadal frequency bands.
-Given the larger damping and lagging effect of the unsaturated zone on the higher-frequency modes, the groundwater variability is more correlated with lower frequent climate circulations such as PDO and AMO, rather than ENSO.
-The Pacific-based indices of PDO and ENSO generally have a higher correlation with Iranian aquifers compared to Atlantic-based indices of AMO and NAO.
Notably, the influence that the Pacific-based indices have on Iranian aquifers weakened after 2000 while relative effects of the Atlantic-based ones intensified after 2000, typically after 2010.
-The aquifers respond to droughts sooner than the wet periods recharge, most probably resulting from the groundwater over-extraction that immediately onsets at the dry season of severe droughts exacerbating the drought's impact and circumventing the damping and lagging effects of the unsaturated zones.
-Multiple interacting effects of AMO + PDO + SOI can better modulate the annual to interdecadal groundwater level variability. This enhanced understanding of the influence of multiple interactions from atmospheric processes observed here would be immensely beneficial for improving the state-of-the-art aquifer behavior prediction.
-Finally, this study shows that the data length can significantly affect the teleconnection results; either the time series are not contemporaneous or only one value of coherence/correlation is computed for all the data length instead of separate computations for different frequency bands and different time spans.
Notably, the wavelet analysis results would be strengthened with longer data records as they, in addition to maintaining higher-frequency signals, also allow to capture the multidecadal modes. However, this study reveals that understanding annual to interdecadal behavior in the natural largescale climate is a universally valuable factor in the longterm predictability of groundwater storage variations and conducting effective groundwater management strategies.