Abstract

Groundwater plays a significant role in domestic, agricultural, and industrial water supply in semi-arid regions. In such areas, rapid growth of population and water demand in conjunction with climate change negatively impacts groundwater quantity and quality. In this research, human activities and climate change effects on groundwater quality in a semi-arid region was studied. First, a numerical groundwater model was calibrated as a tool for simulating an aquifer system. Then, groundwater salinity, as a measure of water quality, was simulated using the gene expression programming ‘GEP’ algorithm. In order to identify major factors influencing the salinity, the grey relational analysis (GR) technique was adopted. Furthermore, in the case study of Mahabad aquifer in northwestern Iran, decreasing precipitation has reduced river flow and aquifer recharge. Crop pattern change, through groundwater exploitation and irrigation return flow, impacted groundwater quality. Changes in temperature, evaporation and crop water requirement has also influenced the surface water-groundwater budget in the region. Groundwater simulations showed a decreasing trend in groundwater level while GEP analysis demonstrated that groundwater electrical conductivity (EC) increased over the past 40 years. Finally, GRA application showed that groundwater withdrawal and agricultural return flow had the highest correlation with increasing EC, compared with the effect of precipitation and temperature as climatic factors.

INTRODUCTION

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.

MATERIALS AND METHODS

Study area

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

Figure 1

Location of the study area in Iran.

Figure 1

Location of the study area in Iran.

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 mm3 and has increased at an annual average rate of 0.88 mm3 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.

Figure 2

Geological map of Mahabad aquifer.

Figure 2

Geological map of 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.

Table 1

The water situation of Mahabad plain during 1971–2000 and 2001–2010

 Dam input Dam release Aquifer withdrawal Precipitation Temperature Evaporation 
 MCM MCM MCM mm °C mm 
Average of first 30 years (1971–2000) 274 328 15.04 314 12.5 1,677 
Average of last 10 years (2001–2010) 223 238 30.1 299 13.1 1,940 
 Dam input Dam release Aquifer withdrawal Precipitation Temperature Evaporation 
 MCM MCM MCM mm °C mm 
Average of first 30 years (1971–2000) 274 328 15.04 314 12.5 1,677 
Average of last 10 years (2001–2010) 223 238 30.1 299 13.1 1,940 
Figure 3

Population in the study region over the 1971–2010 period.

Figure 3

Population in the study region over the 1971–2010 period.

Methodology flowchart

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.

Figure 4

Research flowchart.

Figure 4

Research flowchart.

Groundwater flow modeling

The MODFLOW (Modular Three-dimensional Finite-difference Ground-water Flow Model) (McDonald & Harbaugh 1984) was used to simulate the groundwater flow in this study. The governing partial differential equation for an unconfined aquifer as used in MODFLOW is: 
formula
(1)
where is aquifer transmissivity (L2/T); h is the potentiometric head (L); w is volumetric flux per unit volume representing water sources and/or sinks (T−1); is the specific yield of aquifer (L−1); and is time (T).

Gene expression programming

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).

Indicators used in this study are coefficient of determination (R2), the model fitness function (fi) as represented in Equation (2), and root-mean-square error (RMSE): 
formula
(2)
where T is the observed value, P is the predicted value and R is the selection range.
  • 8.

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

In this study, the relationship between groundwater head and EC is defined by: 
formula
(3)
where w= 1, 2, …, and = 1, 2, …
In order to choose the best relationship between the groundwater level and EC, input variables in different time scales were introduced to the model. Considering R2 and RMSE fitness functions in both training and testing steps, the best pattern is chosen by: 
formula
(4)
where , represent the groundwater level in the last three consecutive months.

Grey relational analysis

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).

In brief, the GRA procedure is as follows. First, according to Equations (5) and (6), reference and comparison series are determined: 
formula
(5)
 
formula
(6)

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

After determining reference and comparison series, data pre-processing is performed. This is a way of transferring the original X0, Xi′ to comparable sequences X0, Xi. Various data pre-processing methods have been proposed for GRA (Kung & Wen 2007). In this study, the following relationship was used: 
formula
(7)
where Xi’ is the original sequences; and Xi is the sequences after data pre-processing.
Determining deviation sequences Δ0i is the next step in which, according to Equation (8), deviation sequences Δ0i between the corresponding values in the reference and comparison series is specified: 
formula
(8)
Calculation of the grey relational coefficient ξi(k) for each series and grey relational grade γi using Equations (9) and (10) is performed: 
formula
(9)
 
formula
(10)

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

RESULTS

Groundwater flow simulation

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.

The increase in groundwater storage due to the naturally occurring precipitation infiltration and percolation is termed ‘precipitation recharge’ (Nasri et al. 2014). In this study, precipitation recharge was determined through the Thornthwait equation (Walton 1970). Required weather data (monthly precipitation, temperature, and evaporation) were obtained from two synoptic stations in Mahabad plain. Irrigation recharge (Ri) was determined based on the total monthly irrigation (I) as follows (Nasri et al. 2014): 
formula
(11)

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 mm3 in 2010. The red cells represent specified head boundary, the blue cells represent Mahabad River and the green cells represent observation wells.

Figure 5

Structure of the MODFLOW model tailored for Mahabad aquifer. Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/aqua.2019.064.

Figure 5

Structure of the MODFLOW model tailored for Mahabad aquifer. Please refer to the online version of this paper to see this figure in color: http://dx.doi.10.2166/aqua.2019.064.

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.

Figure 6

Recorded and simulated heads in (a) central (b) southern and (c) northern parts of Mahabad aquifer.

Figure 6

Recorded and simulated heads in (a) central (b) southern and (c) northern parts of Mahabad aquifer.

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.

Table 2

Measures of error in calibration and validation stages

Piezometer RMSE (m) MAE RMSE (m) MAE 
Central 0.14 0.11 0.14 0.33 
Southern 0.66 0.57 0.68 0.63 
Northern 0.18 0.14 0.28 0.25 
Piezometer RMSE (m) MAE RMSE (m) MAE 
Central 0.14 0.11 0.14 0.33 
Southern 0.66 0.57 0.68 0.63 
Northern 0.18 0.14 0.28 0.25 
Figure 7

Comparison of simulated and observed groundwater head in transient state.

Figure 7

Comparison of simulated and observed groundwater head in transient state.

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.

Figure 8

Groundwater level variations during 1971–2010.

Figure 8

Groundwater level variations during 1971–2010.

Electrical conductivity

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®.

Table 3

GEP model parameters

Parameter Number 
Head size 
Chromosomes 40 
Number of genes 
Mutation rate 0.044 
One-point recombination rate 0.5 
Two-point recombination rate 0.3 
Gene recombination rate 0.1 
IS transposition rate 0.1 
Linking function Addition (+) 
Parameter Number 
Head size 
Chromosomes 40 
Number of genes 
Mutation rate 0.044 
One-point recombination rate 0.5 
Two-point recombination rate 0.3 
Gene recombination rate 0.1 
IS transposition rate 0.1 
Linking function Addition (+) 

Function set 1 = {+ , −, √, exp, ln, x2, x3}

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

The results of EC simulation using function set (1) were reasonable in comparison with function set (2). The RMSE, MAE, R2 and fi 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.

Table 4

GEP model performance measures

EC Test
 
Train
 
fi (test) 
RMSE (μ mohs/cm) MAE MAE MAE 
Central 0.14 0.11 0.11 0.33 664 
Eastern 1.14 1.1 1.1 0.23 56 
Western 0.18 0.16 0.16 0.15 625 
EC Test
 
Train
 
fi (test) 
RMSE (μ mohs/cm) MAE MAE MAE 
Central 0.14 0.11 0.11 0.33 664 
Eastern 1.14 1.1 1.1 0.23 56 
Western 0.18 0.16 0.16 0.15 625 
Figure 9

Simulated-observed EC plots in (a) training stage and (b) testing stage.

Figure 9

Simulated-observed EC plots in (a) training stage and (b) testing stage.

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

Figure 10

Temporal variation of EC during the 1971–2010 period.

Figure 10

Temporal variation of EC during the 1971–2010 period.

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.

Figure 11

Spatial and temporal pattern of simulated EC (μmohs/cm). (a) 1971; (b) 1990; (c) 2010.

Figure 11

Spatial and temporal pattern of simulated EC (μmohs/cm). (a) 1971; (b) 1990; (c) 2010.

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.

Factors influencing EC

Climate change

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.

Figure 12

(a) Annual average precipitation and (b) annual average temperature during the 1971–2010 period.

Figure 12

(a) Annual average precipitation and (b) annual average temperature during the 1971–2010 period.

Human activities

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 mm3/y over the 1971–2010 period whereas the total discharge was 71.8 mm3/y. Thus, groundwater has been over-exploited with a net abstraction of 26.3 mm3/y. The annual average abstraction was estimated at 18.82 mm3/y from 1971 to 2010 with agriculture accounting for 95% of the total withdrawal.

Figure 13

Annual groundwater withdrawal (a), reservoir release (b), and river discharge variations (c) during the 1971–2010 period.

Figure 13

Annual groundwater withdrawal (a), reservoir release (b), and river discharge variations (c) during the 1971–2010 period.

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.

Gray rational analysis

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:

X0′ = {X0′(1), X0′(2), …….. X0′(40)} Electrical conductivity (10) 
X1′ = {X1′(1), X1′(2), …….. X1′(40)} Aquifer recharge due to precipitation (11) 
X2′ = {X2′(1), X2′(2), …….. X2′(40)} Evaporation (12) 
X3′ = {X3′(1), X3′(2), …….. X3′(40)} Groundwater withdrawal (13) 
X4′ = {X4′(1), X4′(2), …….. X4′(40)} River discharge (14) 
X5′ = {X5′(1), X5′(2), …….. X5′(40)} Irrigation return flow (15) 
X0′ = {X0′(1), X0′(2), …….. X0′(40)} Electrical conductivity (10) 
X1′ = {X1′(1), X1′(2), …….. X1′(40)} Aquifer recharge due to precipitation (11) 
X2′ = {X2′(1), X2′(2), …….. X2′(40)} Evaporation (12) 
X3′ = {X3′(1), X3′(2), …….. X3′(40)} Groundwater withdrawal (13) 
X4′ = {X4′(1), X4′(2), …….. X4′(40)} River discharge (14) 
X5′ = {X5′(1), X5′(2), …….. X5′(40)} Irrigation return flow (15) 

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.

Table 5

Pre-processing data

X0 X2 X3 X4 X5 
0.075 0.365 0.000 0.172 0.000 0.322 
0.009 0.313 0.016 0.143 0.004 0.368 
0.000 0.336 0.025 0.565 0.006 0.397 
0.028 0.131 0.016 0.432 0.004 0.187 
0.025 0.219 0.037 0.375 0.139 0.263 
… … … … … … … 
… … … … … … … 
… … … … … … … 
35 0.811 0.302 0.753 0.323 0.888 0.763 
36 0.836 0.107 0.828 0.228 0.863 0.669 
37 0.884 0.347 0.927 0.129 0.553 0.827 
38 0.880 0.380 0.823 0.180 0.637 0.904 
39 0.869 0.370 1.000 0.151 0.86 0.944 
40 0.901 0.179 0.859 0.025 1.000 1.000 
X0 X2 X3 X4 X5 
0.075 0.365 0.000 0.172 0.000 0.322 
0.009 0.313 0.016 0.143 0.004 0.368 
0.000 0.336 0.025 0.565 0.006 0.397 
0.028 0.131 0.016 0.432 0.004 0.187 
0.025 0.219 0.037 0.375 0.139 0.263 
… … … … … … … 
… … … … … … … 
… … … … … … … 
35 0.811 0.302 0.753 0.323 0.888 0.763 
36 0.836 0.107 0.828 0.228 0.863 0.669 
37 0.884 0.347 0.927 0.129 0.553 0.827 
38 0.880 0.380 0.823 0.180 0.637 0.904 
39 0.869 0.370 1.000 0.151 0.86 0.944 
40 0.901 0.179 0.859 0.025 1.000 1.000 
Table 6

The deviation sequences

X1 X2 X3 X4 X5 
0.289 0.075 0.097 0.075 0.247 
0.303 0.006 0.013 0.006 0.358 
0.333 0.025 0.565 0.006 0.039 
0.103 0.012 0.404 0.024 0.159 
0.194 0.012 0.349 0.109 0.238 
… … … … … … 
… … … … … … 
… … … … … … 
35 0.508 0.058 0.488 0.077 0.099 
36 0.729 0.008 0.608 0.028 0.229 
37 0.537 0.043 0.756 0.331 0.094 
38 0.499 0.057 0.7 0. 243 0.011 
39 0.499 0.131 0.718 0.006 0.053 
40 0.722 0.042 0.875 0.099 0.099 
X1 X2 X3 X4 X5 
0.289 0.075 0.097 0.075 0.247 
0.303 0.006 0.013 0.006 0.358 
0.333 0.025 0.565 0.006 0.039 
0.103 0.012 0.404 0.024 0.159 
0.194 0.012 0.349 0.109 0.238 
… … … … … … 
… … … … … … 
… … … … … … 
35 0.508 0.058 0.488 0.077 0.099 
36 0.729 0.008 0.608 0.028 0.229 
37 0.537 0.043 0.756 0.331 0.094 
38 0.499 0.057 0.7 0. 243 0.011 
39 0.499 0.131 0.718 0.006 0.053 
40 0.722 0.042 0.875 0.099 0.099 

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.

Table 7

Calculated grey relational coefficients and grey relational grades for five comparison sequences

X1 X2 X3 X4 X5 
0.582 0.674 0.751 0.853 0.622 
0.570 0.987 0.789 0.802 0.529 
0.543 0.877 0.476 0.474 0.502 
0.814 0.948 0.993 0.559 0.722 
0.681 0.950 0.829 0.595 0.631 
… … … … … … 
… … … … … … 
… … … … … … 
35 0.436 0.978 0.792 0.511 0.813 
36 0.348 0.793 0.770 0.455 0.640 
37 0.442 0.736 0.436 0.401 0.821 
38 0.440 0.537 0.508 0.420 0.995 
39 0.441 0.450 0.535 0.413 0.9 
40 0.350 0.796 0.385 0.336 0.813 
γ0.6 0.705 0.54 0.67 0.62 
X1 X2 X3 X4 X5 
0.582 0.674 0.751 0.853 0.622 
0.570 0.987 0.789 0.802 0.529 
0.543 0.877 0.476 0.474 0.502 
0.814 0.948 0.993 0.559 0.722 
0.681 0.950 0.829 0.595 0.631 
… … … … … … 
… … … … … … 
… … … … … … 
35 0.436 0.978 0.792 0.511 0.813 
36 0.348 0.793 0.770 0.455 0.640 
37 0.442 0.736 0.436 0.401 0.821 
38 0.440 0.537 0.508 0.420 0.995 
39 0.441 0.450 0.535 0.413 0.9 
40 0.350 0.796 0.385 0.336 0.813 
γ0.6 0.705 0.54 0.67 0.62 

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.

DISCUSSION AND CONCLUSIONS

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.

REFERENCES

REFERENCES
Ali
R.
,
McForlane
D.
,
Varma
S.
,
Dawes
W.
,
Emelyanova
I.
,
Hodgson
G.
&
Charles
S.
2012
Potential climate change impacts on groundwater resources of south- western Australia
.
J. Hydrol.
475
,
456
472
.
Bodrud-Dozaa
M. D.
,
Towfiqul Islamb
A. R. M.
,
Ahmeda
F.
,
Dasd
S.
,
Sahae
N.
&
Rahman
M. S.
2016
Characterization of groundwater quality using water evaluation indices, multivariate statistics and geostatistics in central Bangladesh
.
J. Water Sci.
30
,
19
40
.
Chen
H. S.
,
Wang
Z.
&
Gao
Y. H.
2010
Artificial Neural Network Approach for Quantifying Climate Change and Human Activities Impacts on Shallow Groundwater Level – A Case Study of Wuqiao in North China Plain
.
IEEE
,
New York
,
USA
.
Colombani
N.
,
Osti
A.
,
Volta
G.
&
Mastrocicco
M
, .
2016
Impact of climate change on salinization of coastal water resources
.
J. Water Resour. Manage.
30
(
7
),
2483
2496
.
Dettinger
M.
&
Earman
S.
2007
Western groundwater and climate change – pivotal to supply sustainability or vulnerable in its own right?
J. Groundwater
4
,
4
5
.
Fernando
A. K.
,
Shamseldin
A. Y.
&
Abrahart
B. J.
2012
River flow forecasting using gene expression programming
. In
Conference on Hydroinformatics HIC 2012
,
Hamburg, Germany
.
Ferreira
C.
2001
Gene expression programming: a new adaptive algorithm for solving problems
.
Complex Syst.
13
,
87
129
.
Freeze
R. A.
&
Cherry
J. A.
1979
Groundwater
.
Prentice-Hall, Inc.
,
Englewood Cliffs, NJ
.
Gau
H. S.
,
Hsieh
C. Y.
&
Liu
C. W.
2006
Application of grey correlation method to evaluate potential groundwater recharge sites
.
Stoch. Environ. Res. Risk Assess.
20
,
407
421
.
Hisami
M. Z.
,
Shamseldin
A. Y.
&
Melville
B. W.
2011
Statistical downscaling of watershed precipitation using Gene Expressing Programming (GEP)
.
J. Environ. Model. Softw.
26
(
12
),
1639
1646
.
Holman
I. P.
,
Allen
D. M.
,
Cuthbert
M. O.
&
Goderniaux
P.
2012
Towards best practice for assessing the impacts of climate change on groundwater
.
J. Hydrogeol.
20
,
1
4
.
Islam
M. A.
,
Rahman
M. M.
,
Bodrud-Doza
M.
,
Muhib
M. I.
,
Shammi
M.
,
Zahid
A.
,
Akter
Y.
&
Kurasaki
M.
2018
A study of groundwater irrigation water quality in south-central Bangladesh: a geo-statistical model approach using GIS and multivariate statistics
.
Acta Geochim.
37
,
193
214
.
Karanth
K. R.
1987
Ground Water Assessment: Development and Management
.
Tata McGraw-Hill Publishing Company Limited
,
New Delhi
, p.
720
.
Khatibi
R.
,
Ghorbani
A.
,
Hasanpour Kashani
M.
&
Kisi
O.
2011
Comparison of three artificial intelligence techniques for discharge routing
.
J. Hydrol.
403
,
201
212
.
Konikow
L. F.
&
Kendy
E.
2005
Groundwater depletion: a global problem
.
J. Hydrogeol.
13
,
317
320
.
Machiwal
D.
,
Jha
M. K.
2010
Tools and techniques for water quality interpretation
. In:
Advances in Water Quality Control
(
Krantzberg
G.
,
Tanik
A.
,
Antunes do Carmo
J. S.
,
Indarto
A.
&
Ekdal
A.
, eds).
Scientific Research Publishing, Inc.
,
California
,
USA
, pp.
211
252
.
Manca
F.
,
Capelli
G.
,
La Vigna
F.
,
Mazza
R.
&
Pascarella
A.
2014
Wind-induced salt-wedge intrusion in the Tiber river mouth (Rome, Central Italy)
.
J. Environ. Earth Sci.
72
,
1083
1095
.
McDonald
M. G.
&
Harbaugh
A. W.
1984
A Modular Three-Dimensional Finite-Difference Ground-Water Flow Model
.
Scientific Publications Company
,
Reston, VA
,
USA
.
Nadiri
A.
,
Gharekhani
M.
,
Sadeghfam
S.
&
Asghari
M. A.
2016
Groundwater vulnerability indices conditioned by Supervised Intelligence Committee Machine (SICM)
.
J. Sci. Total Environ.
574
,
691
706
.
Nasri
B.
,
Dadmehr
R.
&
Fouche
O.
2014
Water table rising consecutive to surface irrigation in alluvial aquifers: predictive use of numerical modeling. Conference: Chap. 79
. In:
Proceedings of the 12th Congress of IAEG, Engineering Geology for Society and Territory
,
Torino, Italy
,
Volume: 3 (River Basins, Reservoir Sedimentation and Water Resources)
, pp.
379
382
.
Nourani
V. N.
,
Komasi
M.
&
Alami
M. T.
2012
Geomorphology-based genetic programming approach for rainfall-runoff modeling
.
J. Hydroinform.
15
(
2
),
427
445
.
Pulido-Velazquez
M.
,
Peña-Haro
S.
,
García-Prats
A.
,
Mocholi-Almudever
A. F.
,
Henriquez-Dole
L.
,
Macian-Sorribes
H.
&
Lopez-Nicolas
A.
2015
Integrated assessment of the impact of climate and land use changes on groundwater quantity and quality in the Mancha Oriental system (Spain)
.
J. Hydrol. Earth Syst. Sci.
19
,
1677
1693
.
doi:10.5194/hess-19-1677-2015
.
Sara
M. N.
&
Gibbons
R.
1991
Organization and analysis of water quality data
. In:
Practical Handbook of Ground-Water Monitoring
(
Nielsen
D. M.
, ed.).
Lewis Publishers
,
Michigan
,
USA
, pp.
541
588
.
Sha
T.
,
Roy
A. D.
,
Qureshi
A. S.
&
Wang
J.
2003
Sustaining Asia's groundwater boom: an overview of issues and evidence
.
Nat. Resour. Forum
27
,
130
141
.
Shamir
E.
,
Megdal
S. B.
,
Carrillo
C.
,
Castro
C. L.
,
Chang
H. I.
,
Chief
K.
,
Corkhill
F. E.
,
Eden
S.
,
Georgakakos
K. P.
,
Nelson
K. M.
&
Prietto
J.
2015
Climate change and water resources management in the Upper Santa Cruz River, Arizona
.
J. Hydrol.
521
,
18
33
.
Touhami
I.
,
Chirino
E.
,
Andreu
J. M.
,
Sánchez
J. R.
,
Moutahir
H.
&
Bellot
J.
2015
Assessment of climate change impacts on soil water balance and aquifer recharge in a semiarid region in south east Spain
.
J. Hydrol.
527
,
619
629
.
UNESCO & IHP
2008
Groundwater Resources Assessment Under the Pressures of Humanity and Climate Change, a Framework Document, GRAPHIC
.
United Nations Educational, Scientific and Cultural Organization (UNESCO)
,
Paris
,
07 SP. Series No. 2, pp. 6–7
.
Walton
W. C.
1970
Groundwater Resource Evalution
.
McGraw-Hill Book Co.
,
New York
.
Yo-Ping
H.
&
Chi-Chang
H.
1996
The integration and application of fuzzy and grey modeling methods
.
Fuzzy Sets Syst.
78
,
107
119
.