**An Integrated Use of GIS, Geostatistical and Map Overlay Techniques for Spatio-Temporal Variability Analysis of Groundwater Quality and Level in the Punjab Province of Pakistan, South Asia**

**Huzaifa Shahzad <sup>1</sup> , Hafiz Umar Farid <sup>1</sup> , Zahid Mahmood Khan <sup>1</sup> , Muhammad Naveed Anjum <sup>2</sup> , Ijaz Ahmad <sup>3</sup> , Xi Chen <sup>4</sup> , Perviaz Sakindar <sup>5</sup> , Muhammad Mubeen <sup>6</sup> , Matlob Ahmad <sup>7</sup> and Aminjon Gulakhmadov 4,8,9,\***


Received: 24 October 2020; Accepted: 14 December 2020; Published: 17 December 2020

**Abstract:** The rapidly changing climatic scenario is demanding periodic evaluation of groundwater quality at the temporal and spatial scale in any region for its effectual management. The statistical, geographic information system (GIS), geostatistical, and map overlay approaches were applied for investigating the spatio-temporal variation in groundwater quality and level data of 242 monitoring wells in Punjab, Pakistan during pre-monsoon and post-monsoon seasons of the years 2015 and 2016. The analysis indicated the higher variation in data for both the seasons (pre-monsoon and post-monsoon) as coefficient of variation (CV) values were found in the range of 84–175% for groundwater quality parameters. Based on the *t*-test values, the marginal improvement in groundwater electrical conductivity (EC), sodium absorption ratio (SAR) and residual sodium carbonate (RSC) and decrease in groundwater level (GWL) were observed in 2016 as compared to 2015 (*p* = 0.05). The spatial distribution analysis of groundwater EC, SAR and RSC indicated that the groundwater quality was unfit for irrigation in the lower south-east part of the study area. The groundwater level (GWL) was also higher in that part of the study area during the pre-monsoon and post-monsoon seasons in 2015 and 2016. The overlay analysis also indicated that the groundwater EC, RSC and GWL values were higher in south-east parts of the study area during pre-monsoon and post-monsoon seasons of 2015 and 2016. Hence, there is an instant need to apply groundwater management practices in the rest of the region (especially in the lower south-east part) to overcome the future degradation of groundwater quality.

**Keywords:** groundwater quality; groundwater level; geostatistics; *t*-test; spatial distribution modeling

#### **1. Introduction**

Groundwater is becoming a basic requirement for human utilization and crop production [1,2]. Approximately 2.5% of water worldwide is available as freshwater from all the global water resources. The amount of water accessible to humans is available in the form of rivers, lakes, reservoirs, and underground water [3]. Groundwater's function is gaining more importance because a crisis is arising regarding surface water resources with time [4]. In Pakistan, about 50% of water is supplied for irrigation through groundwater sources [5]. Presently, the annual abstraction from groundwater is 60 billion cubic meters [6]. Consequently, the groundwater quality has deteriorated, while the water table has increased in many irrigated areas [7]. The groundwater quality is also rapidly declining worldwide, particularly in developing countries, due to larger dependency on groundwater resources for irrigation [2,8]. Additionally, the deterioration of groundwater quality due to annual and seasonal climate variation may impose pressure on hydrologic and hydrogeologic systems. These seasonal changes were attributed to groundwater quality as a result of monsoonal-driven surface–groundwater interaction [9]. The variation in groundwater quality parameters at a site was mostly related to local conditions and hazards. Similarly, the groundwater level also increased after the precipitation events and then decreased gradually with evaporation. The groundwater level varied within the period from wet to dry seasons and showed seasonal variations because of the seasonal distribution of precipitation and evaporation. Estimation of seasonal flooding impacts on groundwater quality and level due to substantial rainfall is critical to the management of this precious resource [10]. Therefore, it is essential to know the seasonal variation in groundwater quality [11].

Moreover, a key role is played by the various natural processes and anthropogenic activities in degrading groundwater quality [12–14]. To ensure the sustainable safe use of these resources, it is very important to access the groundwater quality as well as its other resources [3]. The groundwater quality is affected by various factors such as regional topography, characteristics of soil, discharge and flow, groundwater circulation through different types of rocks, groundwater recharge, saline water intrusion and hydro-meteorological surroundings of the area [8]. An understanding of the spatial and temporal variation in groundwater quality is essential for sustainable water supplies under changing climate and local environmental pressure [2,15]. For monitoring water quality, traditional approaches have been found to be unreliable due to errors in sampling. Various researchers have developed different graphical and statistical techniques to assess the trends in groundwater quality which vary from simple linear regression to more advanced parametric and non-parametric methods [16–19]. To address the potential groundwater quality problems, a geographic information system (GIS) and trend detection techniques would be useful for examining the long-term water quality variations [20]. The non-parametric methods such as Mann–Kendall (MK), Spearman's rho (SR), Sen's slope estimator (SSE), innovative trend analysis (ITA), the Theil Sen approach and sequential MK have extensively used to detect trends in time series environmental data [21,22]. Similarly, Agca [23] used the *t*-test for the temporal analysis of groundwater quality parameters in Amik plain (South Turkey), which indicated the increasing or decreasing trends of groundwater quality parameters. The management of natural resources can also be achieved using the geographic information system (GIS) at temporal and spatial scale [20,24]. Many authors investigated the GIS contribution in analyzing the spatial distribution of groundwater [25,26].

Therefore, assessment of groundwater spatial distribution is an easy way to analyze the groundwater quality in a matter of its suitability for irrigation [24,27]. The distribution of concentration over space and time that derives from the relation between sample points can be assessed by geostatistical techniques [28,29]. The weighted overlay approach in GIS has also been used to identify the groundwater potential zones in Killinochi, northern Sri Lanka [30]. However, such past studies

did not report the integrated use of statistical, map overlay, and geostatistical techniques for the spatio-temporal variation of groundwater quality and level. The aims of the present research are to use the combine approaches (statistical, geostatistical, and map overlay analysis) for identifying and investigating the spatio-temporal variation in groundwater quality and level in the Haveli Canal Circle (HCC) Command area.

#### **2. Materials and Methods**

### *2.1. Appearances of Study Area*

The present study was conducted in the region of Haveli Canal Circle, Multan irrigation zone (MIZ), Punjab, Pakistan. The region lies between the longitude of 71.12◦ and 72.19◦ and latitude of 29.51◦ and 31.65◦ (Figure 1). The MIZ is situated in an arid and semi-arid regions where the climate is hot in summer and cold in the winter [31]. The highest recorded temperature in the region was 54 ◦C and the lowest recorded temperature was approximately −1 ◦C [27]. It has a flat topography and is suitable for agriculture purposes but receives very little rainfall throughout the year. The average annual rainfall in the study region varies from 100 to 300 mm. The rainfall is unreliable and can be distributed in two seasons (pre-monsoon and post-monsoon). About 60% of the total annual rainfall occurs during the summer season (monsoon rain) and the remaining rainfall is received during the rest of the year [32].

**Figure 1.** Location of study area and sampling points.

Wheat, cotton, sugarcane, and corn are the dominant crops grown in MIZ. The availability of irrigation water from surface sources is a crucial factor for the growth of crops. As the availability of surface water is highly variable, groundwater is pumped to fulfil crop water needs during prolonged dry periods. However, the quality of groundwater is not suitable for irrigation use in many parts of the study area [27]. The aquifers of the study region mostly consist of alluvial deposits which were transported by rivers from the Himalayan mountainous ranges. The alluvium generally contains sands and gravels, or mixtures of sands and gravels. Groundwater in the area is recharged through infiltration from precipitation and seepage from the Chenab River and its associated canal network.

#### *2.2. Data Collection and Analysis*

Two years (2015–2016) of groundwater quality data for 242 monitoring sites (wells) of Haveli Canal Circle, Multan Irrigation Zone (MIZ), Punjab for the pre-monsoon and post-monsoon seasons were collected from the Land Reclamation Department of Multan, Govt. of Punjab, Pakistan. Groundwater samples from every well installed in the study area were collected for pre-monsoon (May–June) and post-monsoon (October-November) meeting the standard protocols. Each sample was examined chemically to identify the groundwater properties of electrical conductivity (EC), sodium adsorption ratio (SAR), and residual sodium carbonate (RSC). The groundwater EC, SAR, and RSC are commonly used to assess the suitability of groundwater for irrigation in Pakistan [33]. Furthermore, the Punjab Irrigation Department also collected data on groundwater EC, SAR, and RSC for groundwater quality monitoring in the study area. The groundwater EC was measured with a portable multi-meter which was calibrated before its use. The SAR and RSC were calculated using the standard laboratory protocol as reported in the literature [34,35]. Similarly, the groundwater level (GWL) data were also recorded using the water level recorder.

The Shapiro–Wilk normality test was conducted to check the data distribution in the Statistical Package for Social Sciences (SPSS) software package. To visualize the normality of the data, normal *Q-Q* plots were used. The *p* values were used to confirm the parameters that showed normal distribution. The normality test showed that the data of all groundwater parameters did not conform to a normal distribution (*p* = 0.05). For exploratory analysis of the data, statistical parameters, i.e., minimum, maximum, mean, median, coefficient of variation (CV) and skewness, were determined using the Statistix 10.0 (Analytical Software, Tallahassee, FL, USA). Spearman's rank correlation analysis was determined in both the seasons for a better assessment of the relationship among groundwater parameters. The Spearman's rank correlation method is a non-parametric test method which can express the level of association between two groundwater parameters [36]. For the temporal evaluation of the 242 datasets during the sampling periods, an independent sample *t*-test was also carried out [23].

#### *2.3. Spatial Variability Analysis*

Spatial variability analysis of groundwater parameters was carried out using the ordinary kriging interpolation technique. Ordinary Kriging is a one of the geostatistical methods that can be used to interpolate a random variable at an unknown location considering its value at the nearby location [37]. Before applying the ordinary kriging interpolation technique, the data were log-transformed to conform to the necessary assumptions for ordinary kriging. The histogram and normal Q-Q plot were drawn using the log-transformed data in the SPSS software package to assess the normality of data and to ensure that the data conformed to the normal distribution. Then, that central tool of geostatistical methods such as semi-variance was applied to quantify the spatial autocorrelation of groundwater parameters based on the log-transformed data [38,39]. All the geostatistical analyses were carried using the GS+ software (RockWare Inc., Golden, CO, USA). Using the binned values fitness method, the best fit model was selected before applying the ordinary kriging interpolations based on the coefficient of determination (R<sup>2</sup> ) and, for scatter plots, each groundwater parameter. Moreover, the ratio of ((C0)/(C0 + C1)) was used to determine the spatial variation of groundwater EC, SAR, RSC, and GWL values. For the given parameters of groundwater, the ratio of ((C0)/(C0 + C1)), i.e., (<25%), (25–75%) and (>75%), indicates strong, moderate, and weak spatial dependence, respectively [40]. The leave one-out cross validation was performed by hiding the values of groundwater samples and estimating its values from the remaining data set. The same procedure was repeated for different semi-variogram models and found the statistics of cross validation error for selecting optimal model. After selecting the suitable model, the kriging method was found to be suitable for the interpolation and management of groundwater parameters using either GS+ software or ArcGIS software package. A study suggested that using the interpolation methods of GS+ software and ArcGIS is efficient in the prediction of unsampled data and gave almost the same result [41]. However, the ArcGIS 10.1 software package was used to interpolate groundwater EC, SAR, RSC, and GWL for both the pre-monsoon and post-monsoon seasons in each year based on the variogram parameters calculated using the GS+ software. A similar approach has also been used for geostatistical and interpolation of top soil properties by Zhang et al. [42]. All groundwater parameters were classified using the manual classification technique based on the groundwater quality standards for irrigation.

#### **3. Results and Discussions**

The minimum, maximum, mean, median, standard deviation (SD), coefficient of variance (CV) and skewness values of groundwater EC, SAR, RSC, GWL for pre-monsoon and post-monsoon seasons in 2015 and 2016 are shown in Tables 1 and 2, respectively. The variability in data can be addressed by the value of CV. A value of CV less than 10% indicates low variability, moderate variability is indicated when 10% ≤ CV ≤ 100%, and high variability is indicated when CV > 100%, respectively [21]. The values of CV for EC during the pre-monsoon and post-monsoon seasons in 2015 were found to be 115.77% and 135.34%, respectively, which indicated the high variability in the data. Similarly, the values of CV for RSC during the pre-monsoon and post-monsoon seasons in 2015 were found as 170.43% and 175.22%, respectively, which also indicated the high variability in the RSC data. Moreover, the values of CV for SAR and GWL during the pre-monsoon and post-monsoon seasons were found to be 84.97%, 92.97%, 79.74%, and 98.29%, respectively, which indicated the moderate variability in the data for both the seasons in 2015 (Table 1). The values of CV for EC during the pre-monsoon and post-monsoon seasons in 2016 were found to be 115.75 and 110.71%, respectively, which also indicated the high variability in the data for groundwater EC. Similarly, the values of CV for RSC and GWL during pre-monsoon and post-monsoon seasons in 2016 were found to be 174.80, 224.19%, 103.35% and 112.47%, respectively, which indicated the high variability in the data for groundwater RSC and GWL (Table 2).

**Table 1.** The statistical summary of groundwater properties in 2015 (*n* = 242).


SD standard deviation, CV coefficient of variation.

**Table 2.** The statistical summary of groundwater properties in 2016 (*n* = 242).


SD standard deviation, CV coefficient of variation.

The paired *t*-test was performed to analyze the difference in mean values of groundwater parameters for pre-monsoon and post-monsoon seasons [43]. The mean EC of 1.49 dS/m during the pre-monsoon season was significantly lower than that of mean EC of 1.65 dS/m during the post-monsoon season in 2015 (*p* = 0.05). The mean values for SAR, RSC, and GWL during the pre-monsoon season were higher than those of post-monsoon season in 2015 (Table 3). Similarly, significant negative difference was determined for groundwater EC between the pre-monsoon and post-monsoon seasons in 2015, having a t-state of −0.887. On the other hand, significant positive difference was observed between the pre-monsoon and post-monsoon seasons for groundwater SAR, RSC, and GWL in 2015. However, the significant positive difference was observed for groundwater EC between the pre-monsoon and post-monsoon seasons, having a t-state value of 0.901 and the significant negative difference was observed for GWL between the pre-monsoon and post-monsoon seasons in 2016 (*p* = 0.05). The difference in results for both the years (2015 and 2016) for seasonal analysis may be due to the variation in rainfall.


**Table 3.** Inequality analysis of means for groundwater quality and level.

\* Significant at *p* = 0.05; electrical conductivity (EC) (dS/m); residual sodium carbonate (RSC) (meq/L); groundwater level (GWL) (m); pre (pre-monsoon season); post (post-monsoon season).

The seasonal variation in groundwater level mainly depends on the annual cycle of rainfall and river water levels [9–11]. The GWL increases due to over-pumping in the pre-monsoon season, which tends to degrade groundwater quality [44]. Annual comparison of the mean value for all groundwater parameters was also observed for both the years (2015 and 2016), which indicated positive significant difference with t-values of 0.291, 1.324, and 0.504 for annual groundwater SAR, RSC and GWL, respectively, at 0.05 significant level (Table 3). The marginal improvement in groundwater SAR and RSC and decrease in GWL were also observed in 2016 as compared to 2015.

## *3.1. Spearman's Rank Corelation Coe*ffi*cient*

The Spearman's rank correlation coefficients (rs) among all groundwater parameters were calculated for correlation analysis in both the seasons. The values of (rs) for each groundwater parameter were shown in Table 4. Interpretation of correlation analysis indicates a quick quality monitoring for groundwater parameters [45]. A significant positive relationship was observed between boring depth and discharge and boring depth and screen length, with r<sup>s</sup> = 0.372 and 0.473, respectively, at 0.05 significant level. The analysis for groundwater EC in the pre-monsoon season of 2015 indicated the significant positive relationship with groundwater EC and SAR in the post-monsoon season with r<sup>s</sup> = 0.930 and 0.729, respectively. Similarly, the analysis of groundwater SAR in pre-monsoon season of 2015 indicated a significant positive relationship with groundwater EC, SAR, and RSC for the post-monsoon season, with r<sup>s</sup> = 0.744, 0.810, and 0.360, respectively. Moreover, groundwater RSC has a significant positive connection in pre-monsoon and post-monsoon seasons with r<sup>s</sup> = 0.842. The significant positive connection between groundwater EC and SAR was also found, with r<sup>s</sup> values varying from 0.297 to 0.783 for both the seasons (pre-monsoon and post-monsoon) in 2015 and 2016. These correlations may indicate some of the impacts of agricultural activities in the study area. This may also be related to the general process governing the groundwater formation [36,46].


**Table 4.** Spearman's rank correlation matrix for groundwater quality parameters.

EC (dS/m); RSC (meq/L); GWL (m); \*\* Significant at = 0.05; \* Significant at = 0.01; pre (pre-monsoon season); post (post-monsoon season).

#### *3.2. Spatial Modelling of Groundwater Parameters*

The spatial dependence of each groundwater parameter for pre-monsoon and the post-monsoon seasons was determined with the help of semi-variograms models for both the years [27]. The ratio of nugget variance (C0) to sill variance (C<sup>0</sup> + C1) was used to determine the spatial variation of groundwater parameters i.e., EC, RSC, SAR and GWL.

The data are considered "strongly" spatially dependent if the ratio is less than or equal to 25%, "moderately" spatially dependent if the ratio ranges from 25 to 75%, and "weakly" spatially dependent if the ratio is greater than 75% [40]. The results indicate that the data of groundwater quality parameters were spatially autocorrelated over the years. From the analysis of the developed semi-variogram, the best fit model was selected based on the coefficient of determination (R<sup>2</sup> ) before the kriging interpolations to analyze the spatial distribution of groundwater parameters (Table 5). The nugget to sill ratios for EC, SAR, RSC, and GWL during the post-monsoon season in 2015 were found to be 13.86%, 12.52%, 15.57%, 13.40% and 8.10%, respectively, which indicated the strong spatial dependence. The nugget to sill ratios during the post-monsoon season in 2016 for EC, SAR, RSC and GWL were found to be 0.03%, 21.60%, 6.84%, 6.33%, 6.33% and 6.38%, respectively, which also indicated the strong spatial dependence. The moderate spatial dependence was analyzed for EC during the pre-monsoon season with nugget to sill ratios of 49.98% (in 2015) and 43.09% (in 2016) for both years. The spatial variability of groundwater properties was affected due to the naturally occurring factors such as topography, regional climate, groundwater flow pattern, groundwater runoff and hydrological conditions. When these naturally occurring factors have greater influence on groundwater, then the chances of strong spatial dependence become maximum for a given area [47]. Similarly, the detailed outcomes about the fitness and selection of different models for the interpolation of groundwater quality parameters are given in Table 5.


**Table 5.** Fitted semi-variogram model of groundwater properties in 2015 and 2016.

C0 Nugget Variance; C0 + C Sill Variance; EC (dS/m); RSC (meq/L); GWL (m).

#### *3.3. Spatial Distribution of Groundwater Quality Parameters*

The spatial distribution of groundwater EC for pre-monsoon season in 2015 indicated that the concentration of EC was higher in the south-east zone of the study area (Figure 2a). The values of EC in this zone were ranged from 2.50 to 4.20 dS/m. A higher concentration of groundwater EC was also found in the northern side for pre-monsoon season in 2015. Similarly, the groundwater EC for the post-monsoon season in 2015 indicated that the concentration of EC was higher in the lower south-east and upper north-west zones of the study area. The lower concentration of EC was observed in the south-west, central south-east, north-west, and north-east parts of the study area and the values ranged from 0 to 1.50 dS/m, respectively (Figure 2b). The groundwater EC for the pre-monsoon season in 2016 also indicated that the concentration of EC was higher in the lower south-east zone and the northern part of the study area (Figure 2c). The groundwater EC concentration showed minor changes in the area, as a higher concentration of groundwater EC was observed in the lower south-east and upper north-east zones (Figure 2d). However, the lower south-east part of the study area continuously showed higher concentration of groundwater EC for both seasons and years. The higher concentration of groundwater EC in that part may be due to the over-abstraction of groundwater. As a result of excessive groundwater extraction, fractured rocks and shales dissolve in the water, which increases the salinity level [48]. Furthermore, the salts tend to accumulate due to strong evaporation in arid and semi-arid areas and consequently, high salinity groundwater forms so that the roots of plants cannot absorb enough water to meet their metabolic requirements [49–51].

Apart from the spatial distribution analysis, overlay analysis was also conducted to analyze the combined effect of the groundwater properties for both seasons. In 2015, the overlay analysis for EC indicated that the lower portion of the south-east zone had a higher concentration of EC in the pre-monsoon and post monsoon seasons. The growth of the high regulation spatial model, GIS, and remote sensing has increased the need for map comparison. The importance of map comparison methods has been recognized and has stimulated growing interest among different researchers [52–54].

The value of EC in this zone was greater than 4.20 dS/m for both the seasons of 2015. The lower concentration of groundwater EC was analyzed in the central south-east and south-west parts of the study area for both seasons of 2015 (Figure 3a). In 2016, the overlay analysis for EC indicated that the lower portion in the south-east and south-west sides had a higher concentration of EC in both seasons (Figure 3b). The area of lower concentration of EC slightly increased in the post-monsoon season as compared to the pre-monsoon season. This indicated that the recharge of freshwater increases in the

**Figure 2.** Spatial distribution of groundwater EC (dS/m) (**a**) for pre-monsoon season in 2015 (**b**) for post-monsoon season in 2015 (**c**) for pre-monsoon season in 2016 (**d**) for pre-monsoon season in 2016.

The spatial distribution of groundwater RSC for pre-monsoon season in 2015 indicated that the concentration of RSC was higher in the upper north-east zone of the study area (Figure 4a). Similarly, the groundwater RSC for the post-monsoon season in 2015 indicated that concentration of RSC was higher in the upper north-east and central south-east and south-west parts of the study area (Figure 4b). The groundwater RSC for pre-monsoon season in 2016 also indicated that the concentration of RSC in groundwater was higher in the upper north-east and central south-east zones of the study area (Figure 4c). Similarly, the groundwater RSC for the post-monsoon season in 2016 indicated that the concentration of RSC was higher in the upper north-east and central south-east parts of the study area

and the values of RSC in that zone ranged from 1.50 to 4.20 meq/l (Figure 4d). In 2015, the overlay analysis for RSC indicated that the upper part of the north-east zone had a higher concentration of RSC in the pre-monsoon and post-monsoon season. The value of RSC in this zone was greater than 4.20 meq/l for both seasons (Figure 5a). Similarly, in 2016, the overlay analysis for RSC indicated that the upper north-east zone had a higher concentration of RSC in both seasons (Figure 5b).

**Figure 3.** Overlay analysis of groundwater EC in pre-monsoon and post-monsoon season. (**a**) EC for 2015 and (**b**) EC for 2016.

The SAR in groundwater evaluates the effect of the sodium hazard with calcium and magnesium concentrations [56]. The spatial distribution of groundwater SAR for the post-monsoon season in 2015 indicated that concentration of SAR, ranging from 5.0 to 10.0, was marginally higher in the upper north-east and lower south-east parts of the study area (Figure 6a). The spatial distribution of groundwater SAR for pre-monsoon season in 2015 indicated that concentration of SAR was lower in the central south-east and south-west parts of the study area and the values ranged from 0 to 5.00 (Figure 6b). It has been reported that seasonal variation affects the groundwater SAR, which can reduce the soil permeability and thus inhibit the absorption of water by crops [57].

The lower concentration of SAR was observed mainly in the central south-east and south-west parts of the study area (Figure 6a). The spatial distribution of groundwater SAR for the pre-monsoon season in 2016 also indicated that the concentration of SAR was higher in the upper north-east and lower south-east parts of the study area (Figure 6c). Higher SAR values in groundwater make soil unfit for plant growth due to a loss of soil permeability. Sodium reduces soil permeability and encourages hardening of the soil. On the other hand, the lower groundwater SAR concentration was observed for the post-monsoon season in 2016 in the upper north-west, north-west, and lower south-west parts of the study area (Figure 6d).

**Figure 4.** Spatial distribution of groundwater RSC (meql/l) (**a**) for pre-monsoon season in 2015 (**b**) for post-monsoon season in 2015 (**c**) for pre-monsoon season in 2016 (**d**) for pre-monsoon season in 2016.

**Figure 5.** Overlay analysis of groundwater RSC. (**a**) RSC for 2015, (**b**) RSC for 2016.

**Figure 6.** Spatial distribution of groundwater SAR (**a**) for post-monsoon season in 2015 (**b**) for pre-monsoon season in 2015 (**c**) for pre-monsoon season in 2016 (**d**) for pre-monsoon season in 2016.

In 2015, the overlay analysis for SAR indicated that the upper part of the north-east and lower south-east and south-west zones had a higher concentration of SAR mainly in pre-monsoon season (Figure 7a). In 2016, the overlay analysis for SAR indicated that the upper north-east zone had a higher concentration of SAR in both the seasons, and the value of SAR in this zone was greater than 15 (Figure 7b). The spatial distribution of groundwater level (GWL) for the pre-monsoon season in 2015 indicated that GWL was high over the lower south-east part of the study area (Figure 8a,b).

**Figure 7.** Overlay analysis of groundwater SAR. (**a**) SAR for 2015, (**b**) SAR for 2016.

It has been reported that with the increase in groundwater level, the osmotic and buoyancy pressure produced by the soil increases to a certain extent. This means that the soil bears a certain additional load within a certain range, resulting in the discharge of pore water from the soil, a reduction in pore water pressure, increase in effective stress, and consolidation and compaction of strata [58]. The GWL for pre-monsoon and post-monsoon seasons in 2016 also indicated that GWL was high in the lower south-east part of the study area. Moreover, the results for the upper north-east part of the study area also indicate that GWL was lower in this zone for both the seasons of 2016 (Figure 8c,d).

In 2015, the overlay analysis for GWL indicated that in the lower south-east part, the groundwater had greater depth in the pre-monsoon and post-monsoon seasons (Figure 9a). Similarly, the GWL results in 2016 indicated that in the lower south-east part, the groundwater had greater depth in pre-monsoon and post-monsoon seasons (Figure 9b). The overall analysis indicated that water quality in the lower south-east part of the study area is deteriorating and unfit for irrigation without possible treatment before its direct application to agricultural lands. The results indicate that potential management measures are needed to ameliorate these impacts, including decreasing groundwater pumping or implementing the activities of groundwater recharge in the region. In order to further reduce the deterioration of groundwater quality and to protect groundwater resources, urgent groundwater management practices are required in the rest of the region.

**Figure 8.** Spatial distribution of groundwater level (GWL) (**a**,**b**) in 2015 and (**c**,**d**) in 2016 for pre monsoon and post monsoon.

**Figure 9.** Overlay analysis of groundwater level. (**a**) GWL for 2015 (**b**) GWL for 2016.

#### **4. Conclusions**

The present study described the integrated use of statistical, GIS, geostatistical, and map overlay approaches for investigating the spatio-temporal variation in groundwater quality and level using the data of 242 monitoring wells in Punjab, Pakistan during the pre-monsoon and post-monsoon seasons of the years 2015 and 2016. The results show that CV values for all the groundwater parameters were greater than 84% for both the seasons (pre-monsoon and post-monsoon) and over the years. The annual comparison for groundwater parameters indicated that the groundwater EC followed a decreasing trend over the years with a 0.05 level of significance. However, groundwater RSC, SAR, and GWL indicated increasing trends over the year, with t-state values of 1.324, 0.291 and 0.504, respectively. However, marginal improvements in groundwater electrical conductivity (EC), sodium absorption ratio (SAR), and residual sodium carbonate (RSC), and a decrease in groundwater level (GWL), were observed in 2016 as compared to 2015 (*p* = 0.05). The results also indicate that the data of groundwater parameters were spatially autocorrelated over the years. The nugget to sill ratios in 2015 for post EC, post-SAR, pre-RSC, post-RSC, pre-GWL, and post-GWL were 13.86%, 12.52%, 15.57%, 13.40%, 7.71% and 8.10%, respectively, which indicates strong spatial dependence. The nugget to sill ratios in 2016 for post-EC, pre-SAR, post-SAR, post-RSC, pre-GWL and post-GWL were 0.03%, 21.60%, 6.84%, 6.33%, 6.61% and 6.38%, respectively, which also indicates strong spatial dependence. The spatial distribution analysis of groundwater EC, SAR, and RSC indicated that the groundwater quality was unfit for irrigation in the lower south-east part of the study area. The groundwater level (GWL) was also higher in that part of the study area during the pre-monsoon and post-monsoon seasons in 2015 and 2016. The overlay analysis also indicated that the groundwater EC, RSC and GWL values were higher in the south-east parts of the study area during pre-monsoon and post-monsoon seasons of 2015 and 2016. The higher values of these parameters in groundwater seem to be due to the simultaneous contribution of natural mineralization, the use of unpurified irrigation water, and anthropogenic inputs. Hence, there is an instant need to apply groundwater management practices in the rest of the region (especially in the lower south-east part) to overcome the future degradation of groundwater quality. The research findings of the present research study provide guidelines for the potential management of groundwater resources in relation to seasonal variation in other regions.

**Author Contributions:** Conceptualization, H.U.F., M.N.A., I.A., and A.G.; methodology, H.S., and H.U.F.; software, H.S., H.U.F., and M.N.A.; formal analysis, H.S., H.U.F., and M.N.A.; data curation, I.A., M.A., and P.S. writing—original draft preparation, H.S., and H.U.F.; writing—review and editing, Z.M.K., M.N.A., I.A., X.C.,

M.M., M.A., and A.G.; funding acquisition, A.G., and X.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was financially supported by the Strategic Priority Research Program of the Chinese -Academy of Sciences, the Pan-Third Pole Environment Study for a Green Silk Road (Grant No. XDA20060303), the National Natural Science Foundation of China (Grant No. 41950410575), the International Cooperation Project of National Natural Science Foundation of China (Grant No. 41761144079), the Xinjiang Tianchi Hundred Talents Program (Grant No.Y848041), the project of the research Center of Ecology and Environment in Central Asia (Grant No. Y934031), and the CAS PIFI fellowship (Grant No. 2021PC0002).

**Acknowledgments:** The authors are grateful to the Department of Land and Reclamation, Multan, for providing groundwater quality data to accomplish this research.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
