Impact of climate variation and human activities on groundwater quality in northwest of Iran

Groundwater, as a water resource, is fundamental to the sustainable development of many societies (Xue et al. 2014). The groundwater science and management community face important challenges due to climate change and its direct/indirect impacts on groundwater (Dettinger & Earman 2007). A lack of major streams in semi-arid regions and rapid growth of population and development has led to over-drafted groundwater resources in recent years. As a result, global groundwater resources may be threatened by uncertain consequences of climate change and anthropogenic activities (UNESCO 2008).

With population growth and economic development, the disparity between water supply and demand is becoming increasingly critical, and people exploit groundwater to meet local water demands (Xia & Chen 2001; Sha et al. 2003; Xia & Zhang 2008). Since 1960, access to pumped wells has caused a rapid worldwide increase in groundwater development for municipal, industrial, and agricultural purposes (Konikow & Kendy 2005). Due to increasing demand for groundwater, several countries, especially in semi-arid and arid zones, are experiencing water shortages resulting from the imbalance between demand and supply (Chen et al. 2010). The combination of climate change and over-exploitation of groundwater has various negative effects such as decline in groundwater level, land subsidence, and saltwater intrusion (Manca et al. 2014) as well as the change in groundwater quality.

Many tools/techniques have been used to study groundwater quantity and quality changes. In regards to groundwater quantity, Shamir et al. (2015) used eight dynamically downscaled global circulation models (GCMs) to measure groundwater response to climate change and concluded that water supply reliability in Arizona region would decrease. Touhami et al. (2015) evaluated the climate change impact on groundwater resources using GCMs and a hydrological model and concluded that average annual recharge will reduce, in the range of 3–17%, in comparison with baseline conditions. Ali et al. (2012) studied impacts of climate change using GCMs and a groundwater model in southwestern Australia and concluded that climate change effects on aquifer systems are expected to be modest.

In terms of groundwater quality, Demirel & Güler (2006) identified anthropogenic factors influencing groundwater quality in the Mersin–Erdemli basin (Turkey) by using hierarchical cluster analysis (HCA), principal component analysis (PCA) and geochemical modeling techniques. Furthermore, Lin et al. (2012) evaluated factors influencing groundwater quality using PCA, HCA and geostatistical techniques in Malaysia. Kim et al. (2014) used model-based cluster analysis to differentiate the contribution of anthropogenic and climatic factors to the groundwater quality condition in South Korea. Their results demonstrated that a bivariate normal mixture model was more robust in comparison with multivariate analysis. Pulido-Velazquez et al. (2015) evaluated the potential impacts of climate and land use change in groundwater resources by coupling (SWAT (Soil and Water Assessment Tool)) with MODFLOW, MT3DMS mass-transport model and Markov chain for quantifying trends in the future. The results showed decreasing groundwater recharge and quality. Machiwal & Jha (2015) explored the spatio-temporal variation of groundwater quality parameters via box-whisker plots and a geographic information system (GIS)-based groundwater quality index. Colombani et al. (2016) used a density-dependent numerical model to quantify the actual and future salinization of a coastal aquifer in the Po Delta (Italy) and showed significant salinization by 2050. Bodrud-Dozaa et al. (2016) applied water evaluation indices and a number of statistical approaches to characterize the water quality in Bangladesh. The study revealed that EC, total dissolved solids (TDS), Ca2+ and Fe values in groundwater samples exceeded thresholds. Islam et al. (2018) established an irrigation water quality index used in south-central Bangladesh and concluded strong to weak spatial dependence of IWQ parameters indicating large-scale natural and small scale anthropogenic effects.

The coupling of hydrological and hydrogeological models (e.g. Holman et al. 2012; Pulido-Velazquez et al. 2015) and sequential coupling of numerical models (e.g. Candela et al. 2009; Narula & Gosain 2013) have been used for groundwater quality simulation. On the other hand, conventional tools/techniques, such as multivariate statistical techniques, PCA and HCA, time series modeling and GIS for graphical and statistical interpretation of groundwater quality have been adopted for the management of groundwater quality (e.g. Freeze & Cherry 1979; Karanth 1987; Sara & Gibbons 1991; Jha et al. 2007; Steube et al. 2009; Machiwal & Jha 2010).

The gene expression programming (GEP) mathematical modeling technique was presented for the first time by Ferreira (2001). This technique is different from some of the other data-driven modeling techniques in that the derived model is not completely a ‘black-box’ and the relationship between the input and the output can be expressed in a mathematical representation (Fernando et al. 2012).

The GEP model has been successfully applied to water resources modeling. Nourani et al. (2012) used GEP for rainfall-runoff prediction. Guven & Aytak (2009) used GEP for the stage-discharge relationship while Hisami et al. (2011) used GEP for statistical downscaling of precipitation in climate change conditions. Nadiri et al. (2016) used GEP, neuro-fuzzy and vector machine methods to assess the vulnerability of groundwater to contamination resulting from anthropogenic activities in Ardabil aquifer, northwest Iran.

Anthropogenic and climatic factors may not be easily separable as far as the groundwater salinity is concerned. Hence, modern approaches and tools such as grey rational analysis (GRA) are necessary for analysis of groundwater quality. GRA is an important part of grey system theory pioneered by Professor Deng in 1982 for handling uncertainties in small data samples with imprecise information. The grey systems have expanded to include GRA, grey modeling (GM), and prediction theory (Xue et al. 2014). In recent years, grey system theory has been widely used in different fields. Wong et al. (2006) used the grey relational method to test for change points in hydrological data and illustrated its applications to the river Shunde in China. Groundwater may be regarded as a grey system due to lack of data and insufficient information (Yo-Ping & Chi-Chang 1996). For instance, Gau et al. (2006) selected suitable groundwater recharge sites using the GRA approach. Xue et al. (2014) used GRA to identify major factors influencing groundwater levels in North China and concluded that human activity was the main factor.

GEP is able to compensate for incomplete or contradictory datasets when studying the vulnerability against groundwater salinity (Nadiri et al. 2016). This technique, without any assumption on the structure of the relationship between dependent and independent variables, identifies appropriate relationships for any given time series (Khatibi et al. 2011). GEP has not been used to predict the EC as a function of groundwater level. Also, GRA has not been used to identify major factors influencing the EC.

In this paper, the Mahabad plain located in northwest Iran is used as a case study in order to differentiate the effect temporal variation of agricultural pattern, construction of dams and groundwater exploitation, representing human activities, on the groundwater system. Moreover, the impact of the change in rainfall and temperature due to climate change on the groundwater was studied. The groundwater system was simulated via a groundwater model while the GEP was adopted for salinity simulations. GRA was also applied to identify the main factors influencing groundwater salinity.

As represented in Figure 1, the Mahabad watershed is one of the Urmia Lake sub-basins while Mahabad aquifer, with an area of 173 km², is located within this watershed.

The Mahabad aquifer is unconfined with a thickness of 30–90 m. The main direction of the flow is from south to north. The aquifer specific yield is 0.03 and hydraulic conductivity varies from 0.3 to 10 m/day. One thousand and six hundred pumping wells with a range of 2–60 L/s are used as a principal water supply in this area. The groundwater plays an important role in the economic and social life of the study region so a significant portion of agricultural and drinking water is supplied from groundwater storage. The annual average pumping volume is 18.82 mm³ and has increased at an annual average rate of 0.88 mm³ over the past 40 years. Moreover, EC has increased to a rate of 7 μmohs/cm over the previous 10 years.

The structural pattern in the region is a result of various tectonic events from the Precambrian to the Quaternary. Qum Formation (Mq) forms the dominant lithology of the region, which has a very important role in recharging the aquifer. Groundwater is recharged by Mahabad River (Qt) as well as through a permeable horizon surrounding the plain boundary, including Quaternary sediments (Qt and Qsm) in some areas. The main structure of Mahabad plain consists of alluvial deposits. The dominant lithology of the region is Qum Formation (Mq) which exerts a major influence on groundwater quality. Together with Mq, Qsm is also influential on groundwater salinity. Figure 2 presents the geological map of the Mahabad aquifer.

The study region is located in a semi-arid climate, receiving an average annual precipitation of 306 mm. Mahabad River is the main surface water drainage system in the study area that flows from south towards north of the plain. Drinking water supply from Mahabad Dam is the first priority and irrigation demand and meeting the environmental demands have second and third priority, respectively. The conjunctive use of surface water and groundwater has been planned to supply for irrigated agriculture. Based on the data supplied by the Regional Water Organization, Western Azerbaijan Province, Table 1 summarizes the water resource condition of the plain while Figure 3 shows the time series of population in the region from 1971 to 2010.

As illustrated in Figure 4, a groundwater model simulates the groundwater level while the EC simulation is conducted by GEP. Major factors influencing the EC are recognized via GRA application.

GEP is an evolutionary algorithm invented by Ferreira in 1999 (Ferreira 2001). The base of GEP is the invention of chromosomes capable of representing any expression tree (ET) (Sattar 2014). GEP involves computer programs of different sizes and shapes encoded in linear chromosomes of fixed length. Due to the linear structure of chromosomes, using genetic operators such as mutation, transposition and recombination always produces accurate and valid structures (Ferreira 2001). In this study, the GeneXproTools 4.0 software package was used to simulate EC based on groundwater level input.

The first step in EC simulation is generation of an initial population via a series of functions and terminals. In brief, the procedure is as follows:

• 1.

  An initial set of control variables such as terminals, function, numerical constants and fitness function are chosen.

• 2.

  Chromosome architecture involving the number of genes, head size, and linking functions are defined.

• 3.

  Contingency of operator and the number of chromosomes as controller GEP are determined.

• 4.

  GEP randomly formulates the chromosomes of the parent program and yields first-generation offspring (Sattar 2014) using implementation genetic operators (mutation, recombination, inversion).

• 5.

  Using Equation (1) as fitness criteria, the fittest offspring is determined by GEP.

• 6.

  With implementation of genetic operators on the selected offspring, many second-generation offspring are generated.

• 7.

  GEP evolution continues as per steps 4–6 until the specified program fitness is met. The final GEP model (the fittest offspring of generation) is scored on a set of performance indicators (Sattar 2014).

• 8.

  Steps 1–7 are repeated with a different set of control variables to produce another GEP model (Sattar 2014).

GRA is a measure of geometric proximity between different discrete sequences, a reference sequence and at least one comparison sequence, within a system. The proximity is described by the GRA, which is a measure of the similarities of discrete data that may be arranged in sequential order (Xue et al. 2014). Detailed information on GRA analysis is available in Wong et al. (2006).

Here X0′ and Xi′ denote the reference and comparison series, k denotes period and I denotes time series.

In Equation (9), ρ ∈ [0, 1] denotes the distinguishing coefficient. Generally, ρ is taken as 0.5 based on the minimum information (Xue et al. 2014).

MODFLOW discretization consists of 100 rows by 50 columns, and one layer, resulting in 4330 active cells. The square cell size is 200 m. Inflow and outflow boundary conditions were defined as specific head (Drichlet) and no flow boundaries were considered as Neumann condition.

Recharge factor (β) is taken as 36.8% based on the results of Nasri et al. (2014).

Figure 5 shows the structure of the developed model as well as the location of observation wells in 2010. Groundwater withdrawal from 1600 wells is estimated at 31.06 mm³ in 2010. The red cells represent specified head boundary, the blue cells represent Mahabad River and the green cells represent observation wells.

The groundwater simulation model was calibrated in steady state in the year 2010 and in transient state over the 12-month stress period from September 2010 to August 2011. The groundwater heads recorded from September 2011 to August 2012 were used to validate the MODFLOW model. The simulated and observed heads in three piezometers in northern, central and southern parts of Mahabad aquifer are presented in Figure 6.

As shown in Figure 6, differences between the observed and simulated values in three piezometers were less than 1 m. Thus, the results are satisfactory. Also, according to Table 2 and Figure 7, there is reasonable agreement between the simulated and observed groundwater level in both the calibration and validation stages.

After calibration and validation of the groundwater model, the groundwater level was simulated over the 1971–2010 period. Variation of annual average groundwater level is presented in Figure 8, based on which it may be concluded that the groundwater level declines at an annual average rate of 0.24 m over the 40-year studied period.

Based on the geological formations in western and eastern parts, the region may be divided into three zones so that the EC could be simulated using GEP in each zone. Some 70% of monthly EC data were used for training and the rest for testing. Two sets of functions were chosen and the following parameter settings (as shown in Table 3) were used for developing the coupled model in GenXProTools^®.

Function set 1 = {+ , −, √, exp, ln, x², x³}

Function set 2 = {+ ,−, × , /}

The results of EC simulation using function set (1) were reasonable in comparison with function set (2). The RMSE, MAE, R² and f_i were considered in evaluation of the performance of GEP model. The results of EC simulation in training and testing stages, as illustrated in Figure 9 and Table 4, indicate a reasonable simulation in both training and testing steps.

The EC temporal variation is illustrated in Figure 10 with an annual average increase of 11.13 μmohs/cm.

According to Figure 11(a)–11(c), groundwater EC has been roughly uniform in the central part of the plain close to the river while in the northeast, EC rises from 989 to 1017 and from 1085 to 1635 during the 1971–2010 period.

In the southern part of the plain, EC is approximately uniform while it varies in the eastern and western parts. EC varies in the east from 1750 to 2800 μmohs/cm and in the west from 850 to 1400 μmohs/cm during 1971–2010.

The climate change impact on water resources in semi-arid regions is expected to be significant. In particular, groundwater quality is affected by recharge due to precipitation as well as phreatic evaporation. Figure 12(a) and 12(b) show the annual average precipitation and temperature during the 1971–2010 period. The result of Mann–Kendall nonparametric tests on annual precipitation and temperature indicated a decreasing and increasing trend in precipitation and temperature, respectively.

In this region, three significant human activities occurred as follows: dam construction in 1969, agricultural development and crop pattern change (since 1990), and groundwater overdraft (since 1990).

During the study period, groundwater abstraction has increased, especially since the 1990s (Figure 13(a)). In Mahabad region, the total groundwater recharge was estimated as 45.5 mm³/y over the 1971–2010 period whereas the total discharge was 71.8 mm³/y. Thus, groundwater has been over-exploited with a net abstraction of 26.3 mm³/y. The annual average abstraction was estimated at 18.82 mm³/y from 1971 to 2010 with agriculture accounting for 95% of the total withdrawal.

The construction of upstream reservoirs impacts groundwater recharge/discharge by altering the surface water regime. As illustrated in Figure 13(b), water released from the reservoir has decreased from 1971 to 2010, due of development in upstream areas combined with decreasing precipitation and increasing temperature. Time series of groundwater withdrawal, reservoir release and Mahabad River discharge are presented in Figure 13(c). A declining groundwater table may be one of the causes for the surface water reduction trend.

Agricultural development and crop pattern change, through groundwater exploitation and irrigation return flow, impacted the groundwater quality. Since 1990, the agricultural pattern in Mahabad region changed such that crops with higher water demands (such as sugar beet, apple, vegetable and alfalfa) were cultivated. Most farmers also planted twice a year.

The GRA was used to rank the impact of climate change and human factors on groundwater salinity. The reference and comparison series are illustrated as follows:

where X0′ is the reference series (representing EC) and Xi′ is the comparison series. River discharge in Equation (14) is one of the water balance components obtained from the groundwater simulations.

After pre-processing data (as illustrated in Table 5), we calculated deviation sequences Δ0i, as presented in Table 6.

The calculated grey relational coefficient ξi(k) and grey relational grade γi of each series using Equations (9) and (10) are presented in Table 7.

According to Table 7, groundwater withdrawal and irrigation return flow are the main factors influencing the EC, whereas effect of river discharge comes last. Thus, we may conclude that the effects of human activities on groundwater salinity are more pronounced than the climate change in Mahabad region.

This study proposes an approach to identify the main factors influencing groundwater salinity in an over-drafted aquifer in northwestern Iran. Groundwater flow was simulated using MODFLOW over the past 40 years. Salinity of groundwater, as a function of groundwater level, was successfully simulated using GEP. Results also revealed that annual average EC increased over the past 40 years, especially in eastern and western parts of the study region.

High EC values in the western and eastern regions indicate that the saltwater–freshwater intrusion has occurred in more than 50% of the study area. The stability of EC values in the central region may be related to the role of Mahabad River.

Climate change in Mahabad region can be characterized by the increasing trend in annual temperature and decreasing trend in annual precipitation. Reduced water recharge and increasing temperature may be responsible for the higher concentration of salt in the aquifer.

Construction of Mahabad Dam as a human intervention impacted the river recharge and thus the quality of groundwater. As illustrated in Figure 13(b), water released from the reservoir has decreased from 1971 to 2010 due to reduced dam inflow as a result of upstream development and climate change. The results showed that upstream development was more effective than the climate change for reducing the dam inflow. The results are in agreement with Wang et al. (2006).

Declining water release exacerbated the exploitation of groundwater because of the declining surface water allocation for irrigation. According to the results of groundwater flow simulation by MODFLOW, net river discharge decreased over the 40 year study period (Figure 13(c)). It appears that declining groundwater level has also been partly responsible for reduction of river discharge.

We used GRA for ranking the effects of climate change and human activities on groundwater salinity. The results revealed that groundwater exploitation and irrigation return flow were the main factors influencing the EC. Thus, we conclude that anthropogenic activities in comparison with climate change was the dominant factor in groundwater quality deterioration. Such a conclusion are in agreement with Xue et al. (2014).

The results of GRA, regarding the change of groundwater quality due to human activities, are in agreement with the results of Machiwal & Jha (2015) and Demirel & Güler (2006). On the other hand, the results of this paper are not matched with Islam et al. (2018) results due to high groundwater exploitation.

Climate change, together with a future increase of the water demand due to further agricultural- and other development plans, require an adaptation plan. Therefore, modernized irrigation technology and new regulations to cover water resources management is urgently needed to achieve sustainable development.

On the basis of this study, appropriate measures to manage groundwater salinity may be adopted in the study aquifer. The findings of this study are useful for policy and decision makers to formulate efficient groundwater use in order to ensure a safe and sustainable groundwater supply.

Sophocleous (2000) noted that the sustainable yield of an aquifer must be considerably less than recharge to sustain both the quantity and quality of streams, and groundwater-dependent ecosystems.