Comparison of parametric and non-parametric methods for trend identification in groundwater levels in Sirjan plain aquifer, Iran

The main objective of trend analysis is to identify the existence or non-existence of significant increasing or decreasing trends in a data series either by using parametric method, for example, the linear regression test, and non-parametric tests, such as the Spearman's Rho (SR) and Mann–Kendall (MK) methods, the two well-known and most widely used non-parametric tests in identifying the trends in climatic parameters (Hirsch et al. 1982; Bui et al. 2012; Güçlü 2018; Zamani et al. 2018; Ali et al. 2019; Noori et al. 2019; Meshram et al. 2020). Hirsch et al. (1991) explained the reasons for using the non-parametric tests instead of parametric methods in various studies analyzing the trends. The parametric methods are developed based on assumptions, such as the normality, the stationarity and the independency of time series, while these assumptions are most uncommon to be satisfied in a hydrological time series. Further, the parametric tests are very sensitive to the presence of the outliers in the data series, which is not the case with the non-parametric methods. Khalili et al. (2016) reported that the non-parametric methods can be applied to both the linear time series, as well as the non-linear time series. While applying the non-parametric test for trend analysis, it is assumed that there is no significant autocorrelation in the given time series. Hamed & Rao (1998) investigated the effect of existence of autocorrelation in time series on trend detection by using non-parametric tests. They proposed a methodology for eliminating the effect of autocorrelation from a data series and applied the proposed method for identifying trends in precipitation and stream flow time series. On using the MK test after eliminating the autocorrelation effect from time series instead of the classic MK test, an increase in the accuracy of trend detection was observed (Hamed & Rao 1998).

In order to investigate the water table variations, many parametric and non-parametric tests have been developed by various researchers for investigating trends in the groundwater table in different parts of the world. Gehrels et al. (1994) studied the variations in surface water runoff and groundwater in the Netherlands and found that the groundwater levels declined due to the drainage, occurrence of drought and over-exploitation of groundwater by farmers. Almedeij & Al-Ruwaih (2006) reported that the groundwater level has a direct relationship with precipitation and indirect relationship with temperature based on the investigations pertaining to the water table in a residential area in Kuwait. Panda et al. (2007) studied the impacts of drought and the effects of over-exploitation on groundwater levels in Orissa (India) by identifying the trends in groundwater levels of 1,002 piezometric wells by using the MK test during the pre-monsoon and post-monsoon seasons from 1994 to 2003. The results revealed that the groundwater declined due to rainfall reduction during dry years and high temperatures, and over-exploitation of the groundwater did not recover through the recharge in wet years. In the post-monsoon (pre-monsoon) season, about 51% (59%) of the piezometers showed declining trends in the groundwater levels in Orissa. Lee et al. (2007) considered the changes in the groundwater levels at Daegu, Korea and claimed that construction of subway tunnels might have led to decline in the groundwater levels. Shahid & Hazarika (2010) investigated groundwater drought using the data of 85 wells located in the northwestern districts of Bangladesh from 1998 to 2002. The over-exploitation of groundwater for irrigation in the dry season and decrease in the aquifer recharge during dry periods were two main reasons for declining trends in the groundwater levels in the northwestern districts of Bangladesh. Shamsudduha et al. (2009) analyzed trends in groundwater levels in the Ganges–Brahmaputra–Meghna delta of Bangladesh during the period 1985–2005 by applying the non-parametric seasonal-trend decomposition procedure. The seasonality dominated the observed variance in groundwater levels, but groundwater level declined at the rate of about 1 m/year in urban and peri-urban areas, and at the rate of 0.1–0.5 m/year in the areas where there was over-exploitation of groundwater resources for rice cultivation during the dry season. However, groundwater levels increased at the rate of 0.5–2.5 cm/year in the estuarine and southern coastal regions.

Tabari et al. (2012) studied the trends in monthly, seasonal, and annual time series of groundwater levels in Mazandaran province of Iran during the period of 1985–2007 using the MK test and reported a combination of negative and positive trends in the piezometeric wells. However, increasing trends were much more pronounced than the decreasing trends in the groundwater levels. Spatial analysis showed concentration of significant increasing trends in the central parts of the province having irrigated tracts of paddy cultivation. Vousoughi et al. (2013) investigated the variations in the groundwater levels and hydrochemical variables in the Ardebil plain, Iran, in 32 piezometers using the MK test. The groundwater levels in the Ardebil plain declined at the rate of about 18 cm per year. Abdullahi et al. (2015) studied the variations in groundwater levels in seven stations of Terengganu, Malaysia by using the MK test and reported that the groundwater levels had experienced a significant negative trend with occurrence of the highest decline during 2005–2009. Xing et al. (2018) reported that the groundwater levels in the Jinan spring area had experienced a significant downward trend from 1956 to 2013. Fang et al. (2019) analyzed trends in groundwater levels in the Beijing plain during the period 2005 to 2013 using the MK test and found decreases in groundwater levels, which might have created future water stresses in the region. Anand et al. (2019) analyzed the variation of groundwater level in the Lower Bhavani River basin, India using GIS (geographical information system) and MK test. The results of the MK test illustrated that the groundwater fluctuation is very high in the southeastern and northeastern parts of the basin, and it is moderate in the northern and northwestern parts of the study area. Lasagna et al. (2020) investigated the trend of piezometric levels of a shallow unconfined aquifer in Piedmont plain (Italy) through the analysis of recorded time series in 37 piezometers during the period of 2002–2017. They used the MK test and reported that half of the piezometers showed no significant trend.

In arid and semi-arid regions, groundwater plays an important role in the economic and social development of the region. Lack of accurate knowledge and over-exploitation of groundwater resources lead to irreversible environmental damage, such as declining groundwater table, land subsidence, decreasing discharges from wells and Qantas, deteriorating water quality, changing groundwater flow direction, and salt water intrusion (Nayak et al. 2006; Sattari et al. 2018). In order to assess the groundwater resources, the status and optimized water management of aquifers, it is highly necessary to consider the water table oscillations quite accurately. Accurate investigation of variations in water table provides valuable information for water resources managers and policy makers to design appropriate plans and policies for meeting the required water demands reliably. In recent decades, the groundwater level of Sirjan plain has declined remarkably. Since groundwater is the main source of water for various needs in this region, it is very important to consider the fluctuations of the water table for optimal water management. In the present study, the trends in groundwater levels in the Sirjan plain were identified by different parametric and non-parametric methods. The effect of autocorrelation on the performance of trend identifying tests was also investigated and the best method for analyzing the trends in the groundwater levels in the Sirjan plain was determined. To the best of our knowledge, this is the first study analyzing the trends in groundwater levels in the Sirjan plain. Investigation of the causative factors of variation in groundwater levels in Sirjan plain (Iran) is another important objective of this study.

Sirjan plain (54° 57′ to 56° 26′ E and 28° 47′ to 29° 58′ N) with a geographical area of 3,982 km² is located in the Kerman province, Iran (Figure 1). The western area of this plain is formed of quaternary deposits consisting of river alluvium, alluvial terraces, sand deposits with small sand-dunes in the southern part of the plain, and salty sandy-clay flats with salt coverage over 10 cm thick in the Kafeh Namak. The bedrock adjacent to the Kafeh Namak is formed from Miocene Marl and is approximately horizontal. The average thickness of alluvium is about 150 m with a reported maximum thickness of 250 m based on the geophysical experiments’ results. The depth of the bedrocks in the northern and northeastern parts of the plain is rising, and then is found to be decreasing towards the southern and southwestern portions of the plain. The average annual precipitation in this plain is about 177 mm. The main sources of the plain recharge in Kerman are the three rivers of Tanguyeh, Hosseinabad, and Ostor, which are usually dry in the summer season. About 98% of the extracted water from the Sirjan aquifer is used in the irrigation of gardens, especially for Pistachio.

The region has traditionally been supported by irrigated agriculture, predominantly with groundwater tapped by Qanats and hand-dug wells. In recent decades, the rising population and the agricultural development in this region has led to an increase in the demands of fresh water. Therefore, both quantity and quality of water resources play a dominant role in the availability of drinking water in the region. However, in recent decades, many deep wells have been dug and huge amounts of groundwater extracted to meet the various water demands of the local population. As a result, the shortage of groundwater has intensified, leading to a decline in the groundwater table and in the deterioration of water quality.

In this study, the monthly data of groundwater levels of 60 piezometers distributed uniformly across the Sirjan plain during 2005–2018 were used for trend analysis. The locations of the considered piezometers across Sirjan plain are shown in Figure 1. In Figure 2, the contour lines of average groundwater levels in the Sirjan plain are plotted. The groundwater flow direction in the Sirjan plain is from northwest to east and south (see Figure 2).

The regression method is used for analyzing the relationship between variables. In this method, the dependent variable is defined as a function of independent variables. The regression method is applied usually for two objectives: determination of relationship among variables and prediction of dependent variables based on the given independent variables.

where, and are the average of independent and dependent variables, respectively. After determining the linear relation among variables, the null hypothesis states that the coefficient is zero, while the alternative hypothesis states that the value of coefficient is significantly different from zero.

In the MK test, the null hypothesis H₀, (there is no significant trend in the time series) will be accepted if at the α significance level, ⁠, otherwise the H₀ will be rejected and the alternative hypothesis (existence of significant trend at the significance level of α) will be accepted (Dinpashoh et al. 2014).

The Innovative Trend Method (ITM), proposed by Şen (2012), was used in this study to analyze trend identification in groundwater levels in annual, seasonal, and monthly time scales. In the ITM test, the studied time series should be divided into two equal sub-series. Then, both the sub-series should be separately sorted in the ascending order. A scatter plot is drawn from two sorted sub-series and one sub-series (X_i) lies on the X-axis and the other sub-series (X_j) lies on the Y-axis (Figure 3). If the data points concentrate on the 1:1 (45°) line of the resulting scatter plot, it indicates that there is no significant trend in the studied time series. Otherwise, if data points are accumulated in the upper (below) area of the 1:1 line, it means that there is increasing (decreasing) trend in the considered time series (Şen 2012).

In the present study, five tests used in identifying trends, namely, linear regression (LR), Spearmen's Rho (SR), classic Mann–Kendall (MK), modified Mann–Kendall (MMK), and Innovative Trend Method (ITM), were applied with the main aim of choosing the best method for analyzing the trends in groundwater levels in the Sirjan plain. The results obtained by using these tests in seasonal and annual time scales are given in Tables 1–5. The box and whisker plots of trend statistics for 60 studied stations (W1 to W60) in seasonal and annual time scales are given in Figure 4. The lines inside the boxes represent the median of trend test statistics; while the lower and upper lines of the boxes display the 25 and 75% percentiles, respectively (see Figure 4). In addition, the lower and upper parts of the whiskers and vertical lines display the corresponding minimum and maximum values of the trend tests’ statistics of series, respectively. The range of trends in groundwater levels obtained through the t-statistic of the linear regression test is larger than the other four non-parametric tests (see Figure 4). The range of Z statistics of the MK and the SR tests are almost the same, while the Z statistic values obtained by using the MMK test were found to have the lowest variance. Therefore, test statistics’ values from non-parametric methods are largely close to each other, but the calculated t-statistics of the LR test are somewhat higher than the corresponding statistics of non-parametric tests because the parametric tests are sensitive to outlier values and are affected by the extreme values. The non-parametric tests are more robust against the outliers and the extreme data in comparison to the parametric tests.

The autocorrelation function (ACF) plots of groundwater levels at the W1 station in different time scales: in the month of January and annual timescale are shown in Figure 5. The lag 1 and lag 2 correlation coefficients in all the considered timescales were found to be significant, and both correlation coefficients became zero after reduction (Figure 5). The existence of autocorrelation in any time series may lead to incorrect detection of significant trend by the test (Hamed & Rao 1998; Dinpashoh et al. 2014). The LR method overestimated trends in groundwater levels in the Sirjan plain because parametric tests have high sensitivity to the outliers. This same issue was observed at the W2 station as no significant trends in groundwater level were witnessed by using all four non-parametric tests in the spring season. On the other hand, the LR test detected statistically significant increasing trends in the groundwater levels because of the effect of autocorrelation coefficients.

The results obtained by using the SR and MK tests, both being non-parametric tests, are not affected by the outliers and the distribution of the data. However, it is still a very important determinant of the autocorrelation coefficients in these tests that leads to unrealistic results. The results obtained through these two tests (the MK and the SR) were found to be close to each other in the case of most of the stations and the same trends at the same level of significance (see Tables 1 and 3). Therefore, the power of these two tests was found to be almost identical.

However in the case of the MMK test, the effect of all autocorrelation coefficients was omitted leading to comparatively smaller values of the Z statistic obtained from this test in comparison to the Z statistics values obtained from both the MK and the SR tests. The effect of autocorrelation coefficients on these tests is clearly visible in the case of station W2. Both the MK and the SR tests showed significant increasing trends at 5% significance level in different timescales (annual and seasonal: summer, autumn, and winter). On the other hand, no significant trends were observed by using the MMK test in these timescales after eliminating the effect of autocorrelation coefficients at station W2. In the cases of stations such as W9 and W16, the three non-parametric tests (MK, SR, and MMK) yielded similar results. After examining the autocorrelation coefficients at both stations (W9 and W16), it was observed that the autocorrelation in the data series of these stations was found to be too weak to influence the amount of computational statistics. Hamed & Rao (1998) showed that application of the MK test after eliminating the effect of all autocorrelation coefficients was more accurate than the classical MK test.

In addition to the MK, the SR, the LR, and the MMK tests, the results obtained by using the ITM method are also presented in Table 5. The results using the ITM method also confirmed the existence of decreasing trends in the groundwater levels in the study area, which is similar to the results of the MMK and the SR methods (Table 5). One of the important advantages of this method is to provide good initial results for decision-making so that it can be used for cognitive stage studies. However, the biggest disadvantage of this method is the failure to provide a decision index and to determine the range of significant and non-significant changes, unlike the MMK and the SR methods that have an index for judging the significance of trend. Therefore, according to the explanations presented in this section, it can be concluded that the MMK test has a better performance than other examined tests in analyzing the trends in groundwater levels, which is the main reason to use the MMK test in this study for analyzing changes in groundwater levels in the Sirjan plain.

Table 6 presents the results of trends obtained through the MMK test in the groundwater levels of 60 wells in the Sirjan plain on a monthly timescale. Results indicate that groundwater levels decreased in most parts of the Sirjan plain. Significant decreasing trends were observed in 56 out of 60 studied wells in all 12 months in the Sirjan plain. Two stations (W9 and W58) experienced an insignificant reduction in the groundwater level in Sirjan plain as well. However, increasing trends, although statistically insignificant, were observed in groundwater levels at W2 and W16 stations in the Sirjan plain. Among the four stations that did not witness significant changes in the groundwater levels, three stations W2, W9, and W16 are located in the eastern parts of the plain, where there is no agricultural activity and over-exploitation for irrigation purposes due to soil salinity and desert connectivity (Rahnama & Mirabbasi 2010). The W58 station, a high-elevation site, also has less agricultural activity in the surrounding areas of the well.