Investigating spatial and temporal trend of groundwater quality in relation to water balance in 2007–2017: a case study of Chaharmahal va Bakhtiari Province, Iran

The extensive exploitation of water resources in Chaharmahal va Bakhtiari province has led to a destructive impact on the water balance and quality of the region. In order to evaluate water quality of the study area, water quality parameters from 132 wells were analyzed to prepare spatial distribution maps of the IRWQI index. To analyze spatial and temporal rainfall anomalies, the SPI index was spatially interpolated using the ordinary Kriging method. Principal component analysis was used to investigate the relationship between water quality parameters. The suitability of data for PCA was evaluated by the Kaiser–Meyer–Olkin and Bartlett tests. Additionally, water balance components of the study area, including surface runoff and ground water, were simulated using the WetSpass-M model. According to the results of the IRWQI index, 54 wells, mostly located in built-up and agricultural lands, had poor quality․ Investigation of the average groundwater quality during the years 2007 to 2017 shows that the trend of groundwater quality decreased. Comparison of drought and water quality maps showed similar patterns so that in areas with extreme drought, water quality was bad. The highest and lowest recorded concentrations for nitrates were related to built-up and rangeland lands with concentrations of 35 and 21 mg/l, respectively. Comparison among land use classes showed that in the rangelands, groundwater quality was better and nitrate level was lower compared to agricultural and built-up areas. Overall, the results of this study show that water quality can be affected by land use types and water balance components.


Introduction
Groundwater is the most important natural source of drinking water for many people around the world (Belkhiri et al. 2020). In some countries, especially in arid and semi-arid areas, dependence on groundwater is very high. Groundwater with water supply capability during periods of drought is still the most reliable and stable source of water, especially in cases in which the surface water supply is limited (Dube et al. 2020). Other benefits of these vital resources include supporting a variety of socio-economic developments and maintaining a wide range of ecosystem functions and services (Mustapha et al. 2019;Shakerkhatibi et al. 2019). Since the quality of groundwater is as important as its quantity, it is necessary to monitor its quality (Neisi et al. 2018;Abbasnia et al. 2019). Groundwater quality change is a function of physical and chemical patterns in an area that is determined by geological and human activities (Belkhiri et al. 2020). Recently, severe groundwater pollution and exploitation have become a global concern due to population growth and widespread demand in the agricultural and industrial sectors (Bhattacharya and Bundschuh 2015). Groundwater can be directly or indirectly affected by climate change and human activities, such as excessive pumping of groundwater and the diffusion of pollutants into the groundwater recharge area (Bundschuh et al. 2017). Climate change will intensify anthropogenic pressures by interfering with the frequency and regime of the aquifer recharge due to changes in rainfall, soil characteristics, and land use (Ricolfi et al. 2020). In addition, climate change affects the recharge/ discharge balance of aquifers by potentially reducing aquifer recharge rates. Therefore, deepening the knowledge of water balance on the regional scale is beneficial for better Extended author information available on the last page of the article management of groundwater resources, both qualitatively and quantitatively (Barbieri et al. 2021). Environmental changes and human activities lead to pollution through excessive concentrations of physicochemical elements in groundwater such as chloride, lead, nitrate, and ammonia (Rashid et al. 2019). In general, the groundwater system is controlled and affected by several complex interactions between different heterogeneous factors. The spatial and temporal variability approach can be effective in analyzing the current hydrochemical image of groundwater (Huang et al. 2018)․ In this regard, it is necessary to identify accurate and reliable methods to understand the spatial and temporal variability of groundwater quality․ Identifying potential toxins in groundwater are also essential to implement strategies to reduce the health risks associated with declining groundwater quality, as well as identifying sources of water pollution (Dube et al. 2020). In recent years, with the increase in the number of physicochemical parameters of groundwater, a wide range of statistical methods are now used for proper analysis and interpretation of data. Most researchers have focused on assessing the spatial distribution of groundwater quality using many statistical methods (Guettaf et al. 2017;Kumar et al., 2014;Singh et al. 2013;Belkhiri and Lotfi 2014)․ Spatial techniques and collected field data have been shown to be one of the most accurate and efficient methods of describing and understanding the spatial distribution of groundwater. For example, the Geographic Information System (GIS) has been widely used as a tool for groundwater quality assessment due to its reliability, application in various environments, as well as its ability to predict and visualize complex data (Rashid et al. 2019;Shakerkhatibi et al. 2019). It also provides new and useful insights into groundwater studies through spatial analysis approaches (Slama and Sebei 2020). The strength of GIS further lies in its role in exploring hydrogeological data and providing spatial information needed to make better decisions (Machiwal et al. 2018). In addition, the application of physics-based hydrological models with GIS to estimate water balance components at the catchment scale has now been improved. Such approaches consider the impact of human intervention and climate diversity on hydrological processes (Alemaw and Chaoka 2003). Consequently, it is possible to estimate the average components of long-term water balance such as surface runoff, evapotranspiration, and groundwater recharge using water balance models integrated in a GIS environment (Alemaw and Chaoka 2003;Kahsay 2008;Julius 2010;Aish 2014;Rwanga and Ndambuki 2017)․ Water balance is defined as the net change of water storage considering all inputs and outputs of the hydrological system (Gebru and Tesfahunegn 2020). Water balance is closely related to ecosystem functions and services because it affects abiotic and biological processes in ecosystems (Mercado-Bettin et al. 2019;Ponce Campos et al. 2013).
Conceptualization of these processes has led to the definition of a large number of models whose suitability depends on the specific case study and its spatio-temporal scale (Abdollahi et al. 2019;Wang and Pullar 2005). The WetSpass model in the GIS environment, which integrates all hydrological processes, has become an excellent option for dealing with more complex hydrological factors in small catchments (Gebru and Tesfahunegn 2020). This model has been used successfully in various studies around the world (Graf and Przybylek 2014;Melki et al. 2017;Molla et al. 2019;Ashaolu et al. 2020;Vu et al. 2020). In this study, the monthly water balance is investigated using the WetSpass-M distribution model. Numerous studies have also been conducted on the quality of groundwater resources (Ebrahimi et al. 2016;Yehia et al. 2017;Khalid 2019;Ratolojanahary et al. 2019;Loh et al. 2020;Sargazi et al. 2021;Masoud and Ali 2020) or the relationship between drought and water quality (Hrdinka et al. 2012;Mosley 2015;Li et al. 2018;Kim et al. 2019;Mishra et al. 2021). But the present study intends to investigate the quality of groundwater from several aspects and examine its relationship with water balance components. To our knowledge, few studies have been conducted in this field. Jeihouni et al. (2018) evaluated the groundwater quantity in the eastern plains of Urmia Lake, Iran during 11 years using a new method of estimating the groundwater balance based on water level table data and 3D modeling. Groundwater quality was also monitored during 10 years using GIS and geostatistics. Also, saltwater intrusion was investigated through the generated quality maps and regression analysis. The results showed that the groundwater balance was negative during the study period. Also, saltwater was spread from the west to other areas, which will increase electrical conductivity, chloride and sodium concentration, and will cause many environmental and agricultural problems. In the study of Nedaw (2010), water balance and hydro chemical characteristics of groundwater in Koraro Area, Tigray, Northern Ethiopia were investigated in order to evaluate the water potential and quality. For this purpose, twelve water samples were collected from hand-dug wells and springs in different lithology and were tested for major cations, anions, and some trace metals. According to the findings, the average annual rainfall, actual evapotranspiration, bare land evaporation, runoff, and groundwater recharge were 548.5 mm, 431.3 mm, 372 mm, 71.33 mm, and 56.7 mm, respectively. The results showed that the groundwater recharge accounts for only 10% of the precipitation, which is due to high evapotranspiration rate and low permeability of surface materials related to high temperature and low humidity. Also, the groundwater quality from the hand-dug wells dug in shale is very poor for both drinking and irrigation purposes. Although water quality is often managed separately from water quantity, it has a close and complex relationship with it. Water quantity is the basic factor of water quality. The nature of this relationship is highly dependent on the individual watershed and the type of aquatic ecosystem, because changes in water quality or quantity may cause immediate changes in the structure and function of ecosystems (Merz 2013). Although water balance and water quality have been investigated in the mentioned studies, the relationship of each water quality parameter with the components of water balance has not been investigated independently, so the present study investigates this issue. Since the outputs of the WetSpass-M model are water balance components maps (an important feature present in few models), the contribution of each well from the values of water balance components can be calculated. Therefore, the relationship between the values of qualitative parameters and water balance components was evaluated in spatial and temporal terms. Also, in this study, the effect of various factors such as drought, land use, amount of nitrates and polluting sources on the quality of groundwater has been investigated. The structure of this article is presented in several sections. First, in the introduction, the importance of groundwater quality and water balance estimation and the studies conducted in these fields is described. Then, the study area is introduced. After that, the trend of temporal and spatial changes in the quality of groundwater resources and water balance in Chaharmahal va Bakhtiari province in Iran has been studied, critical areas (points with low water quality) have been identified using water quality zoning maps, the impact of nitrate and drought on groundwater quality has been assessed using 2017 zoning maps, and, finally, the relationship between the values of water quality variables and the water balance of the upstream watershed has been evaluated and the interdependence among the quality variables has been investigated. Finally, the obtained results are presented.

Study area
Chaharmahal va Bakhtiari province is located in Iran between 49° 28′-51° 25′ E longitude and 31° 09′-32° 48′ N latitude with 16328 km 2 area (Fig. 1). In this study, groundwater data of 132 drinking water wells were used to investigate the effect of land use on drinking water quality in 2016. In terms of geology, Chaharmahal va Bakhtiari province has two different territories, Zagros and Central Iran, along with the Sanandaj-Sirjan zone. In terms of general geology in Chaharmahal va Bakhtiari province, the formations of the third geological period (Zagros region) have intermediate layers and facies of the first period as well. These formations, which geologically belong to the Hormuz series, are seen next to salt domes and faults (Dare Rudi et al. 2015). The set of suitable conditions for the event of the karst process has caused the creation of various water springs in Chaharmahal va Bakhtiari province. The existing major karst springs are often located at the base of karst erosion and are permanent. One of the main characteristics of these springs is their direct dependence on atmospheric precipitation and significant changes in their minimum and maximum water levels (Riahipour and Khalili 2012).

Description of WetSpass-M model
WetSpass (Water and Energy Transfer between Soil, Plants and Atmosphere under quasi-steady state model) is a distributive model used to estimate the long-term mean of spatial patterns of surface runoff, evapotranspiration and groundwater recharge (Wang et al. 2015). This model has been developed in the Department of Hydrological and Hydraulic Engineering of the University of Brussels to predict the transfer of water and energy among soil, plants and atmosphere under quasi-steady state. The original version of this model was presented as an extension in ArcView software on seasonal and annual time scales for temperate humid regions such as Belgium . Recent versions of the model have been developed with the ability to simulate interception, runoff, evapotranspiration, soil water balance and groundwater recharge on a monthly time scale (Mustafa et al. 2017). The latest version of the model is independent of ArcView software and has been developed with the aim of simulating water balance components in arid and semi-arid regions (Abdollahi 2015). In this study, a new and monthly version of this model (WetSpass-M 1.3) has been used. Also, in this study, the baseline flow rate was separated from the observed daily flow data using WHAT software (Web-based Hydrograph Analysis Tool) and the recursive digital filter method. The equation of this method is as follows (Eckhardt 2008): where q t and q t−1 are the filtered base flow at time t and t − 1, respectively, Q t is the total flow at time t (in m 3 /s), α is the constant recession curve, and BFI max is the maximum value of the base-flow index (Eckhardt 2008).

Data collection
Required input data of the WetSpass-M distribution model include monthly rain data, temperature, evaporation from the pan, groundwater depth and wind speed at a height of two meters, number of rainy days per month and digital elevation model maps (DEM), slope, land use and soil texture. Meteorological and hydrometric data were collected from the Regional Water Company of Chaharmahal va Bakhtiari province. Water quality data including electrical conductivity (EC), sodium adsorption ratio (SAR), nitrate, phosphate, total hardness (TH) and the hydrogen potential of water (pH) in 132 groundwater samples from 2007 to 2017.

Research methodology
First, the WetSpass-M model inputs were prepared in order to calculate the monthly water balance in the study area during the statistical period of 2000-2013. For this purpose, monthly maps of rain, temperature, evaporation from the pan, groundwater depth and wind speed at a height of two meters and digital elevation model, slope, land use and soil texture were prepared as raster maps with ASCII format and a pixel size of 500 × 500 m in the ArcGIS and ILWIS (Integrated Land and Water Information System) software environment. The number of rainy days in each month was also prepared as a spatial average. The statistical period of 1996 to 2000 was selected for validation and 2001 to 2013 for calibration of the model, and the performance of the model was evaluated using the Nash-Sutcliffe efficiency coefficient. The model was calibrated manually, using the trial-and-error method and considering the Nash-Sutcliffe coefficient as a target function. The value of the Nash-Sutcliffe efficiency coefficient is calculated as follows: where ENS is the Nash-Sutcliffe efficiency coefficient, Q obs is the observed flow, Q sim is the simulated flow, and Q sim is the average of the estimated flow (Nash and Sutcliffe 1979). Also, data on water quality parameters in 132 wells in the region during the years 2007 to 2017 were collected, and the results of the measured samples were analyzed using SPSS software version 24, and the normality of water quality data was evaluated using the Kolmogorov-Smirnov test.
In the next step, IRWQI GC (IRAN Water Quality Index for Groundwater Resources Conventional Parameters) was used to assess the quality of groundwater resources in the study area. This index is calculated as follows: where w i is the weight of the ith parameter, n is the number of parameters, and I i is the index value for the ith parameter of the ranking curve. Table 1 shows the weight of each parameter, and when there is no value for a parameter, it is In order to determine the descriptive equivalent of the calculated index, Table 2 is used (Hashemi 2011).
The Standardized Precipitation Index (SPI) is one of the most widely used indicators of drought meteorology, calculated on short-term (3, 6 and 9 months) and long-term (12, 18, 24 and 48 months) scales. McKee et al. (1993) developed the SPI index to study drought. Thom (1966) found that the gamma distribution fits the climatic rain time series well. The gamma distribution is defined by its frequency or probability density function: where α > 0 is the shape parameter, β > 0 is the scale parameter, X is the amount of precipitation, and Г(a) is the gamma function.
Gamma probability density distribution parameters are estimated from the sample data using the maximum likelihood method for each station, for the selected time scale, and for each month of the year: where n is the number of rain observations, and x is the average cumulative rainfall for a specific month during the statistical period. Since the gamma function is not defined for (X = 0), zero mm of precipitation and the rainfall distribution may have zero values; the total cumulative probability, which also contains zero values, is obtained from the following equation: where q is the probability of zero rainfall. If m is the number of rainfall data of zero in time series n, then, q is calculated from the following equation: After calculating the total cumulative probability, it is possible to calculate H(X) or the value of the standard normal random variable, which is also probability with mentioned probability, with a mean of zero and a standard deviation of one. The following equations represent the SPI with respect to the values of H(X): Also, based on the results of McKee et al. (1995), different classes of SPI are presented in Table 3.
The geostatistical method is a useful tool for analyzing the structure of spatial variability, interpolation between point observations, and mapping of interpolated values with an  associated error map (Zhou et al. 2011;Arslan 2012). Several studies have reported that groundwater quality is generally characterized by significant spatial variation (Alexander et al. 2017;Maroufpoor et al. 2017). Today, various geostatistical methods are widely used to predict spatial changes in groundwater quality (Jasmin and Mallikarjuna 2014;Belkhiri and Narany 2015;Belkhiri et al. 2017). Kriging is a statistical interpolation approach that consists of several methods, including index kriging, simple kriging, ordinary kriging, and co-kriging, which is commonly used to estimate the spatial distribution of variables (Lee et al. 2007;Dindaroglu 2014;Gyamfi et al. 2016). In this study, spatial distribution maps of IRWQI index, SPI index, and groundwater nitrate were performed using the ordinary kriging method. Recent studies have successfully used advanced geographic statistical programs, such as kriging, to map the spatial variations of groundwater toxins based on point observations (Elubid et al. 2019;Talebiniya et al. 2019;Momejian et al. 2019). A major advantage of geostatistical techniques is that they create a spatial variation of interpolation and also estimate the uncertainty associated with each point (Elubid et al. 2019). For example, the results of studies by Elubid et al. (2019) and Masocha et al. (2019) showed the ability of kriging to contribute to proper understanding of groundwater quality for drinking purposes. Evaluation of the interpolation methods was performed using the Jackknife cross validation procedure. This method is based on the temporary deletion of a sample and its estimation and then returning the sample to the data set and repeating this operation for all available samples. The result of using this method is to create two data sets including estimated values and initial observed values that can be evaluated by calculating different statistical indicators, the results of estimating and optimizing the various characteristics of the site estimator (Mohammadi 2006). To select the best model, the coefficient of determination (R 2 ) and root mean square error (RMSE) was used: where Z* is the estimated value at point x i , Z is the observed value at point x i , and N is the number of points (Wilks 2005). After reviewing the interpolation methods, the best method was selected in terms of the higher amount of R 2 and lower amount of RMSE (higher accuracy) and IRWQI index zoning maps for 2007-2017, and nitrate levels for the years 2016, 2017 and SPI for 2016-2017 were prepared using the appropriate obtained variogram in ArcGIS 9.3 software. Spatial correlation was examined, and the variogram was drawn using version 5.1 of GS+ software. Also, the most important physiographic factors affecting the production of water quality parameters were identified using the principal component analysis (PCA) method in SPSS software. PCA is one of the multivariate statistical methods that help better interpret information when faced with a large amount of information by reducing the complexity of analyzing the primary variables of the problem (Camdevyren et al. 2005).
The suitability of the statistical population for PCA was also assessed by the KMO (Kaiser-Meyer-Olkin) test.
Varimax rotation was also used to improve the relationship between inputs and primary agents and to separate them for membership in agents (Ouyang 2005). In principal component analysis, a correlation matrix is formed and then, the components that cover a small part of the changes in society are left out. Then, the number of main factors that can have a comprehensive description of the variables is determined, and the variables are placed under their group according to the degree of correlation with these components, and the main variables in each component are determined based on the maximum factor load (Lausch and Herzog 2002). The value of the KMO coefficient is between zero and one and is obtained from the following relation: Where r ij is a simple correlation coefficient between the variables i and j, and a ij is a partial correlation coefficient between them. Small values of KMO indicate that the correlation between the pairs of variables cannot be explained by other variables, so the application of factor analysis of variables may not be justified.
As shown in Table 4, if the KMO value is less than 0.5, the data are not suitable for factor analysis, and if it is between 0.5 and 0.69, factor analysis can be done with more caution; if its value is greater than 0.7, the correlations between the data will be suitable for factor analysis (Kalantari 2003). In addition, in order to ensure the appropriateness of the data for factor analysis, the correlation matrix that underlies the analysis is not equal to zero in the population, so the Bartlett sphericity test should be used, which is obtained as follows: where n is the number of subjects, p is the number of variables, and R is the absolute value of the correlation matrix determinants. This statistic has a chi-square distribution with 0.5 p (p − 1) degree of freedom. The amount of information available in |R| is checked by examining the relationship between the number of observations and the number of variables and the probability of error to reject the null hypothesis so that there is no difference from the identity matrix. Identity matrix is a matrix in which in which all diagonal elements are 1, and all non-diagonal elements are zero. The Bartlett test examines the hypothesis that the observed correlation matrix belongs to a community with uncorrelated variables. The variables must be correlated in order for a factor model to be significant; otherwise, there is no reason to explain the factor model. If the hypothesis that the variables are not related is not rejected, factor analysis will not be applicable and should be reconsidered. The significant chi-squared test also indicates the necessary minimum conditions for performing factor analysis (Kalantari 2003). In addition, in this study, the ANOVA test was used to compare the level of nitrate and IRWQI index between different land uses, and the Duncan test was used for multiple comparisons. The ANOVA test is used to compare the mean of two or more communities (the effect of an independent grouping variable on a quantitatively dependent variable). This test shows whether there is a difference between the means of the groups or not, but it is not clear between which groups there is a difference. Therefore, post-experimental tests should be used for this purpose; the most important are the Duncan, Tukey and Scheffe tests (Kalantari 2003). In this study, water quality in areas with pollutant sources was investigated.
First, a table related to the location of pollutants and then, its point file in ArcGIS software was prepared, and three buffers were placed around the position of each contaminant with radii of 500 m, 1 and 5 km. Then, the coordinate system of the created polygon maps was determined and converted to raster maps and then, to ASCII format. After preparing the zoned map of IRWQI index in 2017, water quality based on IRWQI index up to radii of 500 m, 1 and 5 km around each source of pollution was determined in ILWIS software using the cross command. Finally, the water quality around each source of pollution was compared in radii of 500 m, 1 and 5 km. In the last step, the monthly relationship between the values of water quality variables and the water balance of the upstream watershed was investigated. For this purpose, first, the table related to the position of the wells was prepared, and the wells inside the basin were numbered. Then, the location file of the wells was prepared using ArcGIS software, and a buffer with a radius of 500 m was placed around each well, such that some wells were piled on top of each other, and so a number of wells were removed. Next, the coordinate system of the created polygon map was determined and converted to raster format and then to ASCII. After that, the share of each well from the amount of water balance components in the model output maps (monthly runoff, groundwater recharge, actual evapotranspiration and interception maps) was determined in ILWIS software using the cross command for each month. Finally, the relationship between the values of quality parameters and water balance components was evaluated in terms of time and place.

Investigation of the water quality parameters and calculation of the IRWQI index
For water quality parameters, the optimum pH is 6.5-8.6, and its allowable value is 6.5-9. pH is an important indicator for evaluating the quality and pollution of the aquifer systems. Its optimal amount in wells may be an indicator for the presence of alkaline in the groundwater (Ram et al. 2021). Total hardness is distinguished as the amount of dissolved calcium and magnesium in the water. Water moving through soil and rock is an exceptional solvent for calcium and magnesium in that nature sometimes dissolves natural minerals and carries them into the groundwater. The high concentration of TH in water may cause kidney stones and heart disease in human beings (Ram et al. 2021). According to Standard No. 1053 of the Institute of Standards and Industrial Research of Iran, the optimum value of total hardness in drinking water is 1000, and the maximum allowable is 1500; the optimal level of sulfate is 250 mg/L, and its maximum allowable is 400 mg/l; the optimal level of ammonia is 1.5 mg/L, the optimum level of sodium is 200 mg/L, and its maximum allowable is 200 mg/l, which is allowed up to 250 mg/L in the absence of a high quality water source in the region. Sulfate is dissolved and leached from iron sulfides, rocks containing gypsum, and other sulfur-bearing compounds. Sodium is also a highly reactive alkali metal that is present in most groundwater. Many rocks and soils have sodium compounds that dissolve easily and liberate sodium into groundwater (Ram et al. 2021). The high concentration of sodium in groundwater may be related to the mechanism of cation exchange (Kim and Yun 2005). These parameters were desirable in all wells. However, the maximum optimal level of magnesium is 30 mg/L, and 7 wells (5 wells with built-up land use, 1 well with rangeland land use, and 1 well with agricultural land use) were found to have more than the allowable level of magnesium. Magnesium is an important parameter accountable for the hardness of the water (Ram et al. 2021). Also, the optimum turbidity in drinking water is 1 NTU according to Standard No. 1053 of the Iranian Institute of Standards and Industrial Research, and its maximum allowable level is 5 NTU. Six wells (1 well with built-up land use, 1 well with rangeland land use, and 4 wells with agricultural land use) had excessive turbidity levels. The optimal level of chlorine in drinking water according to Standard No. 1053 of the Iranian Institute of Standards and Industrial Research is 250 mg/l, and its maximum allowable is 400 mg/l. The amount of chlorine in the two wells with agricultural and range land uses was excessive. Higher chlorine value in groundwater is dangerous for human health (Pius et al. 2012). With respect to nitrate and nitrite levels, Standard No. 1053 states that the total concentration ratio of each recommended value should not be more than one and the maximum allowable nitrate and nitrite are 50 and 3 mg/l, respectively. In general, nitrate concentration exceeding the permissible limit creates health risks (Krishna Kumar et al. 2015). The high nitrate values in drinking water increase the risk of gastric ulcer/ cancer and other health risks to infants and pregnant women (Rao 2006). Nitrogen contained in bedrock contributes to nitrogen saturation, which leads to leaching of nitrogen and the increasing of nitrate concentration in groundwater (Ram et al. 2021). In 2017, nitrate levels were only allowed in wells with built-up and agricultural land uses. This can be attributed to excessive use of nitrogen fertilizers and leaching of soil nitrate due to irrigation in these areas. According to the results obtained from the IRWQI index for each well, 54 wells had poor quality; according to the land use map of the region, 20 wells were related to built-up use and 24 wells were related to agricultural land use․ Reasons for the low quality of water in these wells include the use of chemical fertilizers, the expansion of urban and industrial wastewater and livestock. These results were consistent with the study conducted by Smith et al. (2020) in the Yucatan Peninsula aquifer. The reason for the low quality of groundwater in these land use areas is linked to the infiltrated agrochemicals and manure used in the wet season.
In the case of zoned water quality maps, the higher the IRWQI index, the higher the water quality. The status of water quality in each year should be considered according to the water quality classification table in terms of the value of this index ( Table 2). The average groundwater quality of the region during the years 2007 to 2017 is decreasing and in terms of spatial changes, and groundwater quality is different in different regions. It may be better to study the groundwater quality situation spatially and temporally, because areas whose water quality is declining over time have been identified and studies can be conducted to investigate the factors affecting groundwater quality reduction, in order to manage eliminate pollutants and to prevent further pollution in these areas. Spatial and temporal changes in groundwater quality according to the value of IRWQI index in different places in each year varies and no specific trend can be observed for all cities in the study area. Figure 2 shows the zoned water quality map based on the IRWQI index in the cities of the study area in 2007 and 2015. In 2007, the best water quality was observed in the city of Lordegan, where the water quality was in the relatively good class, and the worst water quality was in the northeastern and eastern parts of Shahrekord, where the water situation was in the medium class. In 2015, the water quality of the city of Ardal was in the good class, while the water quality in the northern and northeastern areas of the city of Borujen was reduced and in the bad class.

Comparison between groundwater quality and nitrate level
In order to compare groundwater quality and nitrate level, nitrate zoning maps were prepared in 2016 and 2017. The Fig. 2 Comparison of water quality in the study area in 2007 and level of nitrate in the groundwater of the study area in 2017 was lower than in 2016. Of course, the level of nitrate varies in different cities. Ardal had less nitrate, and Borujen and Shahrekord had more. Figure 3 presents the comparison of the results of zoned maps of groundwater quality based on the IRWQI index with the level of nitrate, for example, in the zoning map of water quality and the level of nitrate in 2017. This figure shows that the quality of groundwater can be affected by the level of nitrate because the color pattern of these maps is almost identical. As the amount of nitrate increases, the quality of groundwater decreases. In the areas of the maps colored green, the level of nitrate is lower and the IRWQI index is higher, which indicates better groundwater quality. In both years, the groundwater nitrate level in Ardal has been less; on the other hand, water quality has been better in this city.

The results of studying the effect of drought on water quality using the SPI index
In order to investigate the relationship between drought and groundwater quality, the drought zoning map of the study area based on SPI index was prepared using the kriging method in ArcGIS software for 2016-2017. According to Fig. 4, the lower and upper limits of this index are equal to − 2 and 0.39, which according to the classification of the values of SPI index in Table 3 indicates extreme drought and a near-normal situation, respectively. Also, according to Fig. 4, groundwater quality can be affected by drought, as demonstrated by the color pattern of these maps being almost the same. So, with the decrease in SPI index, the IRWQI index has increased. For example, a comparison of these maps shows that in parts of Borujen that are under extreme drought according to the SPI index, water quality was also bad according to the IRWQI index. Also, most of Ardal, in which water quality is in normal condition based on the SPI index, according to the IRWQI index in 2016, had relatively good water quality and also had good water quality in 2017.

Results of factor analysis
In this study, in order to analyze the main components and the main factor analysis, SPSS statistical software was used, and the appropriateness of the statistical population was investigated using the KMO test. Table 5 presents the results of the KMO test. According to this table, since the value of the KMO statistic equals 0.613, factor analysis can be done. The results of the Bartlett sphericity test were also significant, in the sense that there was a significant correlation between the variables. Table 6 shows the share of variables in the factors after rotation. Each variable is placed in a factor that has a high correlation with that factor. According to this table, the first three components explain about 73% of the changes.

Results of comparing land use in terms of nitrate level and IRWQI index
The ANOVA test was used to compare nitrate level and the IRWQI index between different uses. The results of the ANOVA test showed that there was a significant difference between different land use rangelands, and agricultural and built-up use at the level of 0.01. The Duncan test was used to investigate the difference between land use in nitrate level  and IRWQI index (Figs. 5 and 6). Figure 5 shows a difference between rangeland and built-up uses in terms of IRWQI value, but there was no difference between agricultural use and rangeland or between agricultural and built-up use.
In addition, the highest and lowest values of IRWQI index were related to rangeland and built-up uses with the values of 52 and 45, respectively. According to Fig. 6, there was a difference between rangeland and built-up uses and between rangeland and agriculture in terms of nitrate level, but there was no difference between agricultural and built-up uses. Also, the highest and lowest recorded concentrations for nitrate ions were related to built-up and rangeland uses with concentrations of 35 and 21 mg/L, respectively. The high level of nitrate in agricultural and built-up land uses can have several reasons. In another study, Alvi et al. (2005) reported the source of nitrate contamination of groundwater resources due to groundwater contamination with human and animal waste.
The results of the Lee et al. (2003) study investigating nitrate concentration using GIS showed that the agricultural use had the greatest effect on increasing nitrate concentration in low rainy seasons and built-up use also increased the concentration during the rainy season. Determining the nitrate source and reducing it is the most effective way to prevent nitrate contamination of groundwater, improve groundwater quality, and possibly reduce the discharge of nitrate to these sources (Anca 1999). Smith et al. (2020) also showed that agricultural land use has a vast effect on groundwater quality, due to producing a higher concentration of nitrate, ammonium, potassium and electrical conductivity.

Results of water quality study in polluting sources
Water quality was determined according to the IRWQI index up to the radius of 500 m, 1 and 5 km around each source of pollution and water quality around each source of pollutants in these radii. Table 7 shows water quality in areas with polluting sources, and Table 8 shows water quality in landfills. In these tables, it is expected that as the distance from the landfill and the source of pollution increases, the IRWQI index will increase, and the groundwater quality will improve. Although the changes in water quality in these radii were very small, Table 8 shows that only in some areas (numbered 1, 2, 12, 15, 16, 25 and 26) did the water quality improve by moving away from the landfill. This is not true in other areas, and it can be said that other factors have been involved, the most important of which is human interventions. Given that the trend of water quality changes from a radius of 500 m to a radius of 1 km in almost all areas, it is possible that in the selected radii in this study, especially at the radius of 5 km (a long distance), there was another source of pollution. Furthermore, some groundwater quality parameters such as nitrate under the influence of the wastewater collection network and soil physical characteristics

WetSpass-M model calibration and validation results
The outputs of the WetSpass-M model are in the form of raster maps as well as a file containing a summary of the spatial averaging of the parameters. Investigation of model performance using the Nash-Sutcliffe efficiency coefficient showed satisfactory results; the values of this coefficient were 0.62 for runoff and base flow in the calibration period and 0.57 and 0.55 for runoff and base flow in the validation period, respectively.

Water balance simulation results
According to the WetSpass-M model's simulation results of water balance components, the average annual share of runoff, groundwater recharge (groundwater recharge, soil moisture and subsurface flows) and evapotranspiration (evaporation from soil surface, transpiration from plants and interception) from 544 mm of rainfall were 17% (95 mm), 54% (292 mm) and 24% (119 mm), respectively. Therefore, on average, the highest share of water balance components from the annual rainfall of the basin was related to groundwater recharge, and the results obtained were consistent with those of Soleimani Motlagh's (2016) studies in the Sarab Seyed Ali watershed in western Iran and Bayati (2016) in the Vanak watershed located in Chaharmahal va Bakhtiari and Isfahan provinces. The trend of changes in water balance components and the average annual rainfall during 1996 to 2012 can be seen in Fig. 7. According to this figure, the trend of groundwater recharge changes had an almost similar slope to the average annual rainfall of the study basin during the statistical period 1996-2012. Runoff also had a downward trend with a gentler slope, but evapotranspiration and interception had a gentle upward trend. The results of evaluating the relationship between qualitative variables and water balance components The output maps of the WetSpass-M water balance model include actual evaporation, groundwater recharge and runoff maps on a monthly time scale. In this study, the share of each well from the distributed amount of water balance components in the basin for each month was determined using ILWIS software and based on the geographical location of the study wells. Given the available information and in order to investigate the relationship between quality parameters and water balance components, the statistical period 2007 to 2012 was selected. Water quality parameters were measured in different months of each year in the selected statistical period. Therefore, the relationship between the values of qualitative parameters and the equivalent values of actual evapotranspiration, groundwater recharge and runoff estimated by the model in these months during the years 2007 to 2012 were investigated. Therefore, the coefficient of determination between the qualitative parameters of EC, SAR, pH, nitrate, phosphate and total hardness and water balance components including actual evapotranspiration, runoff and groundwater recharge were calculated during the study period for each well. The highest values of correlation coefficients related to the linear relationship between the actual evapotranspiration values and the qualitative parameters EC, SAR, pH, nitrate and total hardness are shown in Table 9.
The results related to the highest correlation coefficients between runoff and quality parameters are presented in Table 10.
Regarding the correlation between groundwater recharge and water quality parameters, the results varied in each well so that the correlation was observed in a small number of the 45 study wells. The results are shown in Table 11.
In general, the results showed that changes in water quality parameters were a function of fluctuations in monthly evapotranspiration, runoff and groundwater recharge in some  wells. Thus, the values of SAR, pH, nitrate and total hardness increased with increasing evapotranspiration in some wells (Table 9). Also, the values of EC and total hardness decreased with increasing runoff in some wells (Table 10). The results of land use and soil texture of the areas related to the existence of these wells also showed that there was correlation among water quality parameters and evapotranspiration and runoff on average in clay and silty-clay loam textures and medium rangeland use. With increasing groundwater recharge, EC values decrease, and SAR and pH values increase in the mentioned wells (Table 11), in which there was on average silty-clay loam soil texture and medium rangeland and built-up uses. However, as mentioned, the correlation between groundwater recharge and water quality parameters was observed in a small number of study wells in the basin. In addition, the results showed no correlation between phosphate levels and water balance components. In addition, a parametric correlation test (Pearson coefficient in SPSS) was performed in SPSS software to investigate the correlation between water quality parameters and water balance components on an annual time scale. The results of this test showed that considering the Pearson coefficient of 0.6 as a strong correlation, there was a correlation between evapotranspiration and pH in 2007 (with increasing evapotranspiration, the pH of water increases) and between groundwater recharge and SAR in 2009 (SAR decreases with increasing groundwater recharge), which were significant at the 0.01 level (Table 12). Groundwater recharges naturally through rivers, streams, and lakes. An increase in surface water supply and the pollutants may be associated with an increase in the hydraulic and contaminant loading. As a result, in the presence of contaminating factors, it might lead to potential for groundwater quality problems (Lee and Jones-Lee 1995). The evaporation process is followed by the dissolution process and slow movement of groundwater from the recharge areas to downstream increases the hydraulic conductivity of water samples in the direction of the groundwater flow (Abdolahi et al. 2015). Also, no strong correlation existed between water balance values and quality variables. Most parameters had poor correlation.

Conclusion
In the process of evaluating water resources, our quantitative knowledge of these resources is formed. This is important for the sustainable management of the water resources of Chaharmahal va Bakhtiari province, which is a major supplier to several other provinces. In order to study the trend of groundwater quantity, the fully distributed WetSpass-M monthly water balance model was used. Fully distributive models consider the spatial variations of hydrological parameters and calculate the parameters for each pixel due to their ability to connect to GIS. Although the input data of these models are high, because of their capabilities they have received more attention. Results of model performance using Nash-Sutcliffe efficiency coefficient showed satisfactory results for runoff and base flow. The results of distributed simulation of water balance components also showed that the highest share of water balance components from the average annual rainfall of the basin was related to groundwater recharge. The quality of water is as important as it is necessary for the health of the environment. Various factors such as droughts and land use are influencing factors on groundwater quality. The similar color pattern in SPI index and nitrate value maps with the obtained water quality index map confirms the impact of drought on the groundwater nitrate. The effect of land use on drinking water quality in the study area showed that among the three land uses of rangelands, agricultural and built-up uses, water quality was better in rangelands and their level of nitrate was less than the other two land uses. Nitrate zoning maps also showed that groundwater quality could be affected by nitrate level, because in areas with high nitrate level, groundwater quality was lower according to the IRWQI index. The reasons for the high level of nitrate in built-up and agricultural uses compared to rangeland can be attributed to sewage, excessive use of nitrogen fertilizers and leaching of soil nitrate due to irrigation in these areas. Another result of this study for 2016-2017 was that drought was an effective factor in reducing water quality; in areas, where drought was more severe, groundwater quality was worse. The results of comparing the water quality of Chaharmahal va Bakhtiari provinces also showed that the Ardel administrative division had predominantly higher quality as it has lower nitrate levels and higher IRWQI values. Regarding the investigation of groundwater quality, water quality zoning maps based on the IRWQI index were prepared for the years 2007 to 2017. In general, the groundwater quality of the study area decreased over time.
In our study, the trend of spatial and temporal changes in the quality of groundwater resources in Iran's Chaharmahal va Bakhtiari province was investigated. Quantity and quality of groundwater face risks such as declining  levels and reducing groundwater recharge due to excessive increase in the drilling of illegal wells and reducing rainfall, increasing population and urban areas, followed by increasing urban and industrial wastewater, development of agricultural lands and increasing the use of chemical fertilizers. As a result, groundwater is polluted, and its quality is reduced. These problems have led to many studies in this field in many different parts of the world.
The results of these studies show the importance of proper management and planning in the field of water resources. In addition, the results related to the relationship between groundwater quality parameters and water balance components showed that this correlation was weak on an annual scale and was strong on a monthly scale in some wells. Consequently, the values of SAR, pH, nitrate and total hardness increased with increasing actual evapotranspiration and runoff in some wells located in clay and siltyclay textures and medium rangelands. Regarding the correlation between groundwater recharge and water quality parameters, the results were different in each well, and a correlation was seen in a low number of the 45 study wells. EC values decreased with increasing groundwater recharge and SAR and pH values increased in some wells located in silty-clay loam soil and in the use of medium rangelands and built-up. Another result of this part of the study showed no correlation between phosphate levels and water balance parameters. Better understanding of water quantity and quality may help to improve water management, lead to proper planning for water exploitation and to applying best management practices.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.