Recognizing the relationship between the groundwater level and hydrological time-series: a case study of the Ardabil Plain Open Access

Groundwater (GW) is an important factor in meeting urban, industrial, and especially agricultural water demands, particularly in tropical and subtropical semiarid regions (Siebert et al. 2010). GW systems have some characteristics, such as complexity, nonlinearity, multi-scale, and randomness, which are influenced by natural and/or anthropogenic factors. Therefore, it is important to detect such changes using an accurate measurement of the variability contained in the time-series of the process (Nourani et al. 2015). Many researchers pointed to the effect of climate parameters on the decrease of the groundwater level (GWL) time-series (Zwolsman & van Bokhoven 2007; Waibel et al. 2013; Chinnasamy & Ganapathy 2018). The impact of human activities such as GW abstraction, recharge, and reservoir construction on the GW level fluctuations (Xue et al. 2014; Singh et al. 2016; Yang et al. 2017; Amaranto et al. 2018; Deng et al. 2018).

Recently, various methods for calculating the complexity of time-series and signals have been presented in different fields of science and engineering. The wavelet-entropy measure (WEM) can be used to identify the effective factor of fluctuation change. Shannon (1948) introduced the entropy concept to obtain additional information about time-series. The Shannon entropy concept has been used in numerous research papers to analyze signals (Bercher & Vignat 2000; Shardt & Huang 2013; Chen & Li 2014; Castillo et al. 2015; Singh & Cui 2015; Varanis & Pederiva 2015).

The combination of entropy and wavelet concepts has been used to develop a new complexity WEM (Rosso et al. 2006). In addition to the aforementioned studies, several methods have been proposed in the past decades to measure complexity and thus detect and model changes in time-series. For example, Fathian et al. (2016) used Seasonal Auto-Regressive Integrated Moving Average to study the changes in the water level of Lake Urmia. This article aims to identify the changes in statistical characteristics in terms of trend, stationarity, linearity/nonlinearity, and change point detection analyses. Vaheddoost & Aksoy (2017) calculated the entropy of each proposed station with respect to the long-term mean precipitation of the basin. Although few studies have been conducted in the field of watershed engineering that addresses complexity changes (Li & Zhang 2008), biomedical studies have found that deep sleep or anesthesia, illness, and ageing in humans lead to a decrease in the complexity of associated physiological signals (Goldberger et al. 2002). It can be concluded that WEM is a new and efficient index for determining the complexity of time-series, especially hydrological time-series. Komasi & Sharghi (2019) used WEM as a criterion for the degree of time-series fluctuations. Their WEM results showed that the decrease of the aquifer water level in the Silakhor Plain indicates the decrease of the natural GWL time-series fluctuations. The results showed that the decrease in the fluctuation of the runoff time-series has a greater influence on the GWL oscillation patterns compared to the rainfall and temperature time-series. Also, it can be concluded that the climatic parameters are not subjected to significant changes; thus, human activities play a main role in the decrease of the GWL in the Silakhor Plain.

Hydrological time-series are often non-stationary, and various factors such as climate change, human activities, etc. can influence patterns (Nourani et al. 2015). The presence of seasonal variability in hydrologic processes leads to accurate calculation of complexity and fluctuations using wavelet transform coherence (WTC) measures. A useful mathematical tool such as the WTC to measure the relationship between rainfall and runoff is an important step in restoration projects in the Ardabil Plain. The main purpose of the WT analysis as a function of time is to decompose a signal into sub-series at different time frequencies (Danandeh Mehr et al. 2014). For example, Grinsted et al. (2004) used the WT analysis to determine physical relationships between geophysical time-series. Holman et al. (2011) applied the WT approach to detect temporal time-space nonlinear relationships between North Atlantic Oscillation (NAO) atmospheric teleconnection and GWLs. They used the continuous WT and the cross WT (XWT) to distinguish cross-wavelet power of time-series. Tremblay et al. (2011) used correlation and wavelet analyses as well as wavelet coherence to investigate effect relationships between four climate indices (the NAO, the Arctic Oscillation, the Pacific North American pattern, and the El Niño-Southern Oscillation (ENSO) represented by the Multivariate ENSO Index), the GWL time-series, and precipitation and temperature time-series in three Canadian regions. The three Canadian regions studied show drastically different patterns of variability trends in hydrogeological records. Yu & Lin (2015) applied the XWT integration to investigate the non-stationary time–frequency relationship between precipitation and the GWL variations. Their results showed nonlinear and non-stationary rainfall–recharge relationships of a GW system, which may be frequency and spatially due to different frequencies. In a similar study, Henderson et al. (2009) applied the XWT and continuous WT to identify minute fluctuations in GW as a result of daily pumping in the submarine. Kuss & Gurdak (2014) use Singular Spectrum Analysis, the WTC, and lag correlation to quantify the effects of the ENSO (2- to 7-year cycle), NAO (3- to 6-year cycle), Pacific Decadal Oscillation (PDO) (15- to 25-year cycle), and Atlantic Multi Decadal Oscillation (AMO) (50- to 70-year cycle) on precipitation and GWLs in the United States. The results suggest that GWLs are driven in part by inter-annual to multi decadal climate variability and do not depend solely on temporal patterns in pumping. ENSO and PDO have a greater influence on the GWL variations in the United States than the NAO and AMO, particularly in the West and Central regions. Yu & Lin (2015) proposed an integration of the XWT and empirical orthogonal functional analysis to analyze the space–time nonlinear relationships between the rainfall and GW changes. Their method revealed three major space–time patterns of GWLs in the results. The XWT also identified the lagged effects between the rainfall and GW changes. Duvert et al. (2015) analyzed the hydrodynamic response of an agricultural watershed in southeast Queensland, Australia, to low- and high-frequency precipitation variations that occurred over a 25-year period. Results showed that strong internal variations in the precipitation input affected surface water flow more than GWLs. Statistically significant episodes of the WTC were found in a 2- to 4-year band between the Niño3.4 index and GWLs for the upstream piezometers, especially during the 1995–2000 period, which may be related to a strong La Niña event. Oh et al. (2017) developed the combined application of dynamic factor analysis and wavelet analysis to identify complex latent factors controlling the GWL fluctuations in a riverside alluvial aquifer influenced by barrage construction. They found that the key driving forces controlling the GWL time-series data in a complex hydrological setting can be identified and quantitatively assessed through the combined use of dynamic factor analysis and WT and the application of the WTC. Nourani et al. (2018) applied the WTC to identify the effects of hydro-climatological time-series on water level variations in lakes (Lake Urmia and Lake Van) in the Middle East. Upon investigation, the relationships between runoff and the water level showed maximum values (0.9–1) in the two lakes. Drewnik et al. (2018) investigated the variability of GWLs and GW temperature in raised bogs in the Bieszczady Mountains in southern Poland. The WTC results show that the most significant response of peat bogs to weather conditions was observed in summer and autumn. Neves et al. (2019) investigated the relationships between important large-scale atmospheric circulation modes and inter-annual to decadal oscillations in GWLs using the WTC in Portugal. Results show non-stationary relationships, but consistent across different period bands. The contributions of the eastern Atlantic and Scandinavia at a relatively high frequency (<5-year period) are difficult to separate, but their combined influence accounts for about 20 and 40% of the total variance in GWLs in the south and north of the country, respectively. Zhang et al. (2019) used the WTC to analyze the response of GWLs to semi-diurnal tides. The results show strong correlations on a time scale of 0.5, 1, and 15 days (resonance periods), which are then used as prediction periods for artificial neural network models. The prediction results also confirm that semi-diurnal tides and precipitation have a great influence on the GWL, with better prediction in the filled layer. Rezaei & Gurdak (2020) used Singular Spectrum Analysis, the WTC, and lag correlation calculations to analyze and quantify the effects of ENSO, NAO, PDO and AMO on hydro-climate variables of precipitation, temperature, lake level, GW fluctuations, soil moisture, vegetation cover, and the insolation clearness index in the Lake Urmia watershed. Moderate coherence between the PDO and GWLs in most adjacent aquifers was found during the >8-year period from ∼1980 to 2015. Malakar et al. (2021) examined the long-term effects of local precipitation, global climate cycles, and human influence on the GWLs at multiple depths using lag correlation analyses, the WTC, and regression-based dominance analyses. They observed intuitive responses, i.e. a rapid response in shallow GWs and relatively delayed responses to global climate patterns with increasing depth.

According to previous research, the WEM and WTC are new indexes and an efficient measure to determine the complexity and frequency relation of hydrological time-series. Therefore, in the present study, the effect of rainfall and runoff changes on the GWL of the Ardabil Plain aquifer was investigated for the first time using the WEM and WTC. To study the reduction of the GWL in the Ardabil Plain, it is important to determine the causes of these changes in the GWL. To demonstrate the effects of climate change, rainfall time-series are studied and to monitor the effects of human activities, such as digging illegal wells, and the time-series of runoff in the Ardabil Plain are used.

In this study, the effect of rainfall and runoff changes on the GWL of the Ardabil Plain aquifer is investigated using the WEM and WTC criteria. The hydrological time-series are very complicated. The use of the WT by decomposing the time-series into sub-signals can analyze the time-series and provide accurate short- and long-term information at different levels of resolution (Rajaee et al. 2010; Komasi & Sharghi 2019). The detailed abbreviations and definitions used in the paper are listed in Table 1.

The discrete wavelet coefficient for scale a=2^m and location b=2^mn is defined by T_m,n. A finite time-series, in Equation (20) is considered where N=2^M. Regarding this concept, the boundaries of m and n are and ⁠, respectively.

W_m(t)(m=1, 2, … , M) provide the detailed signals which are able to catch small features of interpretational value in the data and the residual component, T̅(t), denotes the historic information of data.

As a result of the WT and to capture the temporal relationships among non-stationary time-series, the WTC was proposed to determine the localized correlation coefficients and their phase relationships among non-stationary signals in time–frequency spaces in which a detailed explanation of the WTC can be found in Torrence & Compo (1998).

For two different time-series, X_n and Y_n (n presents time-scale), with different WTs of and ⁠, XWT is ⁠, where (*) presents the complex conjugate. The cross-wavelet power is defined as ⁠. The XWT identifies regions in time–frequency space where two time-series show high common power, and thus, significance (Nourani et al. 2018).

The study zone of this study is located in the plain of Ardabil (38 °03′–38 °27′N and 47 °55′–48 °20′E) which is in the northwest of Iran and covers an area of about 990 km², as shown in Figure 1. The annual rainfall in the Ardabil Plain is about 304 mm. The wettest month is May and the driest month in the region is August. The average temperature in the Ardabil Plain is about 9 °C. It is noteworthy that this plain is known as the coldest region in Iran. The average number of frost days in the Ardabil Plain is 130 per year. The main rivers of this plain are Balikhli-Chay, Qara-Chay, and Qara-Su, which of course do not flow all year round and fill up the aquifers of the Ardabil Plain. The Balikhli-Chay and Qara-Chay flow into the Qara-Su in the northern part of the Ardabil Plain (see Figure 1). In this study, the GWL of 15 piezometric stations (P1, P2, P3, …, P15) in the Ardabil Plain from the period of 1993–2018 were selected to perform the trend analysis. The GWL altitude in the study zone varies from 1,308 to 1,529 m above mean sea level. Figure 1 shows the locations of GW monitoring stations in the study area. Table 2 shows the statistical characteristics of the Ardabil Plain. In this study, the monthly rainfall, runoff, and the GWL time-series for 1993–2018 are used (300 months) (see Figure 2).

The GWL faces a decrease of 11.43 m in the level of the piezometers in the Ardabil Plain from 1993 to 2018. It is possible to study the changes in the runoff and rainfall time-series using the WEM and to examine climate change and human factors and their interactions in the complexity decrease of the GWL. According to Figure 2, the mean GWL in 1993 was approximately 1,341.45 m, maximum and minimum values were 1,346.75 and 1,335.40 m, respectively. In 2018, the GWL in the Ardabil Plain decreased by 11.43 m in a period of 25 years due to human activities and the effects of climate change. In this way, the mean GWL time-series is divided into three sub-series and the measure WE is calculated for each period. The GWL time-series sub-series consist of 100 monthly data. The GWL sub-series, which are decomposed into five sub-series by the db2 parent wavelet, contain several frequent time-series at levels 1–5. The decomposition of the time-series into different scales by the WT leads to a structural interpretation of the series and provides insights into the frequency domains and history of the signal (Rajaee et al. 2010; Nourani et al. 2015). The db2 wavelet function was selected based on the similarity between the db2 signal shape and the fluctuations of the GWL time-series compared to other wavelet functions (Komasi & Sharghi 2019). The WEM as a complexity measure was applied to each of the three sub-series of the GWL time-series. Finally, the energy in normal form (ρn) was calculated for the decomposed sub-series of the GW signal (levels 1–5). The results of the normalized energy for the decomposed GW signal are shown in Table 3. The WEM in the second time period (2001–2002 to 2009–2010) shows a significant decrease, which represents the complexity decrease in the GW signal of the Ardabil level. As shown in Table 3, the WEM of the GWL signal in the third part of the time period (2009–2010 to 2017–2018) shows a decrease of 28.9922%. The decrease of the WEM in a certain period shows an unfavorable trend in the GWL signal. Also, the decrease of the WEM indicates a decrease in the complexity or fluctuations of the GWL signal in the third period. This indicates the existence of undesirable trends in the GWL of the Ardabil Plain, and the main objective is to evaluate the causes of the decrease in the GWL due to human activities and climate change factors. The WEM of the GWL decreased and showed significant changes. To find the main reason for the decrease of the GWL in the Ardabil Plain, runoff and rainfall signals are examined. To investigate the relationships between rainfall (or runoff) and the GWL parameters, the WTC method is used.

To achieve this goal, the runoff time-series as factors of human activities (Komasi & Sharghi 2019) and rainfall time-series of the Ardabil Plain as factors of climate change (Komasi & Sharghi 2019) are divided into three time sub-series with the same number (100 months). Then, using the db4 mother wavelet, each time sub-series was decomposed into multi-frequency time-series at decomposition levels 1–5 as the GWL parameters. In this part, we try to identify human and climatic reasons for the decline of the GWL. Similarly, divided sub-rainfall (or runoff) can be identified as the dominant reason for the GWL decline. The WEM was finally estimated for three sub-series of rainfall and runoff signals. Tables 4 and 5 show the WEM for rainfall and runoff time-series. From Table 4, it can be seen that the rainfall time-series has a WEM reduction of about −2.94016 and −1.360983% in the second and third time periods, respectively. Thus, in the second and third time periods of the rainfall time-series, there was a reduction of −2.94016 and −1.360983% in the fluctuations, respectively. Table 5 also shows that the WEM has an increase of 0.86613% and a decrease of −1.57813% for the runoff time-series in the second and third time periods, respectively, which has no significant effect on the decreasing GW. The runoff vibration rate in three sub-series shows an increase in the second part and a decrease in the WEM value in the third part. It was found that human activities led to decreased runoff in the Ardabil Plain in the third period, i.e. the complexity of runoff decreases in the third period. The construction of dams and irrigation as human activities lead to a decrease in runoff, so the complexity decreases in the third part. The WEM of rainfall and runoff time-series in Tables 4 and 5 can be compared with the WEM of the GWL in Table 3, which shows that both rainfall and runoff signals play a minor role in reducing the WEM for the GWL signal in the third part of the period. The increase of the WEM of the GWL signal in the first period could be due to the increase of the WEM of the runoff signal.

Figure 3 shows the WEM changes in three sub-series for rainfall, runoff, and GWL parameters in three 100-month periods. This figure showed the WE changes in the time-series for rainfall, runoff, and the GWL in the best form. Moreover, this figure shows the WEM decrease in the GWL, indicating fluctuations or decrease in complexity of the GWL signal in the Ardabil Plain, which is much stronger than the runoff and rainfall series in the third period. It can be inferred that the reduction in the complexity of the runoff is stronger than the reduction in the precipitation parameters in the third period, which affects the reduction in the complexity of the GWL time-series, none of which significantly contributes to the WEM reduction in the GWL time-series. Recently, human activities such as over-exploitation of GW have led to a reduction in the GWL in the Ardabil Plain. As a result, exploitation of GW plays an important role in the decrease of the GWL compared to climate parameters such as rainfall time-series.

To evaluate the coherency and also the seasonal relationship between hydrological time-series (runoff and rainfall) and the GWL time-series, the XWT and WTC were applied to time-series of the Ardabil Plain. Coherence works like correlation. A dark color indicates that two time-series are highly correlated, bold black lines 1 (yellow) and remaining 0 (blue) indicate no or low correlation. The areas within the bold black lines indicate the times and periodicities with statistically significant XWT and WTC at the 5% significance level. The WTC calculates the cross-correlation of two time-series as a frequency function (at different wavelet scales).

The XWT and WTC results between the GWL and rainfall signals of the Ardabil Plain are shown in Figure 4(a) and 4(b). Figure 5(a) and 5(b) show the XWT and WTC between runoff time-series and the GWL.

Before applying the WT to the original rainfall, runoff, and the GWL time-series, you need to standardize them (mean = 0 and variance = 1). As for the coherence results, Figures 4(a) and 5(a) show the frequency of 8–16 months between rainfall and the GWL time-series, and there are also common periodicities between runoff and the GWL time-series, which shows the most coherency in most time period. Figures 4(b) and 5(b) show that the frequency of 8–16 months between the hydrological time-series and the GWL time-series in the XWT plots has low correlation.

Moreover, the wavelet coherence results between the runoff and GWL signals showed almost similar behavior between rainfall and the GWL time-series due to the frequency band.

Among the WTC graphs between rainfall and the GWL signals and runoff and the GWL signals, the runoff parameter showed a high coherency value in 8–16 month periods. A high degree of coherence occurs in most months over a 12-month time scale. Uniform rightwards with delay less correlation confirmed high coherence.

The decline of the GWL due to errors in water resources management has led to a significant problem for human society and the environment. In this study, the multi-scale method WEM and WT coherence method are used to identify the link between the GWL decline and rainfall or runoff time-series to accurately determine the impact of human activities and climate change on the GWL time-series in the Ardabil Plain. In this study, runoff as a factor of human activities and rainfall as a factor of climate change play a role in hydrological processes. In this way, the GWL, rainfall, and runoff signals of the Ardabil Plain are divided into three 100-month sub-series with decomposition levels 1–5 in the WEM calculation. The rainfall and runoff sub-series were divided into several frequent sub-series using the db4 mother wavelet, and the GWL sub-series were decomposed using the db2 mother wavelet. Finally, the WEM was calculated for each sub-series of the GWL, rainfall, and runoff. The results show that the WE measure of the mean GWL signal decreases in the third time period. In the third time period, the WE measure of rainfall and runoff signals is decreasing, but the main reason for the WE decrease in the GWL signals is other factors of human activities (since 2009–2010 to 2017–2018).

The results of the WTC and XWT analysis to detect maximum common local multi-scale correlations and phase relationships between rainfall (or runoff) and the GWL time-series in the Ardabil Plain showed that the time scale of 12 months has a high coherency.

To complete this study, some recommendations for future research are suggested. For example, the study of other hydro-climatological parameters that affect the reduction of the GWL time-series of the Ardabil Plain, such as human intervention and climate variability in recent years. It is also proposed to apply the methodology of this study to other hydro-climatological parameters (e.g. temperature, transpiration, etc.).

The authors have no relevant financial or non-financial interests to disclose.

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.