Integrated geophysical technique for groundwater salinity delineation, an approach to agriculture sustainability for Nankana Sahib Area, Pakistan

Abstract Owing to the climatic dryness and shortage of surface water, Pakistani farmers have been forced to rely on groundwater resources to meet the irrigation demand for crops. Pumping an excessive amount of water from the aquifers is turning it brackish. Our study aims to select a suitable depth of freshwater for installing tube wells in a study area, Tehsil Nankana Sahib. Sixteen vertical electrical soundings (VES) were arrayed in this study area, and water samples were also collected to record their electrical conductivity. The interpolation maps for different parameters were generated by using the inverse distance weighted interpolation method. The obtained results depict that the central and southeastern region has good groundwater at depths from 60 to 90 m, the northeastern region has marginally fit water, while the northwestern region is marked as having an unfit groundwater zone. This research article is not only a case study but also can be used for other areas of Punjab, Pakistan, facing similar groundwater problems. However, the comparison of VES results with borehole data is critical to achieve better accuracy. An extensive exploratory borehole study of the area is recommended to further characterize the subsurface lithology.


Introduction
All over Pakistan, groundwater salinity for irrigation purposes has become an important topic of discussion (Anjum et al. 2016). Pakistan's climate mainly comprises of arid to semi-arid conditions and a distinct irrigation system principally dependent on surface water availability, which is primarily supported by the Indus basin (Aftab et al. 2018). The canal system has not been so effective in some areas due to a shortage of surface water. So, the aquifers are the main source of irrigation for crops in arid areas of Pakistan (Raza et al. 2017). The farmers have been contrived to rely on groundwater supplies due to drought conditions and global climate change . In Pakistan, the uncontrolled usage of groundwater resources is increasing, and in turn, it is pushing demand for water supplies. Without any hydrogeological survey, numerous tube wells have been installed in the country.
Overall, the number of private tube wells in the country was estimated at around 957,916, including 824,879 in Punjab, according to data published by the Ministry of Food, Agriculture, and Livestock, Government of Pakistan, Islamabad (Sikandar et al. 2010). In the areas having canals with insufficient surface water supply, minimal irrigation using groundwater results in degradation of groundwater supplies with an increase in salinity (Sikandar et al. 2017). Aside from excessive groundwater extraction, one of the main reasons for this distressing water situation is the lack of planning to develop more water storage (Raza et al. 2017). Approximately 57% of farmers are reported to keep farms falling because of irrigation water scarcity. Approximately 30% of farms have either a problem with salinity or some other restriction and are bound to have less cultivation than that available in Rechna Doab (Jehangir et al. 2002). The current study area has no significant surface irrigation sources in the area due to which farmers are forced to rely on the groundwater for irrigation of crops. The maximum thickness of any freshwater storage zone will always be near a recharge source like a river, canal, or stream or areas having the highest rate of rainfall (Greenman et al. 1967). However, this case is not applicable for the Ravi river and Sutlej River as they carry sewerage water with them, which pollutes the groundwater (Khan & Khan 2020). In addition, areas having the maximum thickness of freshwater where there is no recharge resource might be associated with the position of former stream courses (Greenman et al. 1967).
The methods commonly used to study the subsurface, such as boreholes, are very costly. Since drilling several boreholes to find out the depth of freshwater is not very economical and feasible, there is a need for a cost-effective and easy method to deeply study the groundwater salinity trends in the area (Gee et al. 2020). So, the current study was planned, which aims at marking the variation of groundwater salinity in the area and selecting the targeted depth of installing tube wells in the prospected area. The main point is to develop long-standing sustainability from the aquifers (Todd and Mays 2005). So, the study of groundwater resources becomes important day by day to fulfil the demand for agriculture and industrial use (Dhakate et al. 2015). Geophysical techniques are important tools to study the groundwater resources for being one of the cheap, non-invasive, and non-destructive techniques.
One geophysical technique named vertical electrical sounding (VES) is utilized in our work. Geoelectrical techniques are often used to determine the thickness, depth, and aquifer boundary (Omosuyi et al. 2007;Asfahani and Abou 2013). The VES has been widely used for hydrology studies worldwide (Coker 2012;Asfahani and Abou 2013). The Electrical Conductivity (EC) data is also used in this study to determine the salinity of groundwater. The bulk electrical resistivity can be extracted from the EC, and the resistivity can be used to separate the subsurface layers. The bulk electrical resistivity varies with changes in lithology, structure, fracturing, local topography, weathering thickness, degree of interconnected particles, and water content (Dhakate et al. 2015;Sikandar et al. 2017). A comparison of VES results with EC results of water samples obtained from the already installed nearby tube wells will help to attain better accuracy in terms of interpretation of good ground water zones.
Sometimes it is very difficult to delineate the areas having freshwater, due to intermixing of the true resistivities and their thicknesses over a large area. In resistivity sounding technique, the different resistivity models can produce resistivity curves that are essentially the same, this phenomenon is called the principles of equivalence and suppression (Singh et al. 2004;Sanuade et al. 2019). To mitigate this problem Dar Zarrouk (DZ) parameters are introduced. The DZ parameters can be defined with longitudinal conductance, longitudinal resistivity and transverse resistance (Utom et al. 2012). The DZ parameters computed from the 1 D geoelectrical method aid in correctly identifying the saline and freshwater interfaces without intermixing and confusion. The importance of these parameters in establishing understandable vision of the occurrence and distribution of fresh and saline water aquifers, has been demonstrated by Hasan et al. (2019). The DZ parameters only give an estimation of the saline and freshwater interface of the area, this is the reason its comparison with the EC of groundwater and borehole lithology data is important. To study the aquifer vulnerability towards contamination, integrated electrical conductivity (IEC) is also derived from VES results (Gemail et al. 2017). This will help us to demarcate the areas with a greater chance of getting contaminated in the future.
The VES has been extensively used in other areas of Rachna Doab to study the distribution of groundwater salinity (Sikandar et al. 2017;Khalid et al. 2019). The groundwater quality in the intermediate layers or shallow depth is far superior to that of the deep layers, which have a maximum thickness and poor water quality . Aftab et al. (2018) developed the geological model that contributes to the variation of soil water composition after intermixing to aquifers to study the variation of the water quality and salinity trends in the area. The total dissolved solids (TDS) for groundwater varies all over the area in Rechna Rachna Doab from 300 to 2500 mg/L (Aftab et al. 2018). The study by Gilani (2013) states that there were very undetectable changes in the quality of these water samples within these 6 months.
The lithological explanation will be limited due to the lack of borehole data for the study area. The limitation of EC data is that it only provides approximations for EC and TDS ranges. To further classify the quality of water, detailed laboratory measurements are required. Three data sources, including the VES, EC, and drainage pattern of the area, were used in the current study to reduce errors in the interpretation of the saline and freshwater interface.
Overall, the VES will provides subsurface variation in resistivity ranges that can be used to interpret the quality of groundwater in the area. Furthermore, the IEC calculated from VES and EC values will assist us in identifying the areas most vulnerable to contamination. The drainage pattern analysis will enables us to identify areas with few paleo-stream channels and measure surface runoff. All of these data are important for understanding aquifer recharge and will be used in our work.

Study area
The study area Tehsil Nankana Sahib, as shown in Figure 1, lies in the Rechna Doab, of Punjab Province, Pakistan. The area of Rechna Doab is 28,000 km 2 . Its southeastern boundary is marked by the river Ravi, and the northeastern boundary is marked by the river Chenab (Raza et al. 2017). The research work is dedicated to the northern limb of Nankana Sahib Tehsil ( Figure 1) with an area of 1626 km 2 and is extended to Tehsil with a coverage of 580 km 2 . Nankana Sahib is an important agricultural area, and the rate of cultivation in this area is high (Gul et al. 2020). Wheat, rice, sugar cane, vegetables, tobacco, roses, and marigold are the most dominant and economic crops cultivated in this area (Riaz and Javaid 2012). The area has extreme weather conditions in summer and mild conditions in winter, while the average annual temperature is 24.1 C (Farooq et al. 2017).
The area has a 500 mm annual average total rainfall, which is not enough to support the irrigation requirements of crops cultivated in this area (Yameen et al. 2019). Moreover, the canal irrigation system has not been enough to support agriculture. Hence to justify this water need, agriculture in this area is mainly supported by groundwater (Jehangir et al. 2002). In the study area, farmers were facing a problem with groundwater salinity and the gradual deepening of the tube well water. As stated in a recent study, the area is termed a saline groundwater zone or highly contaminated zone in the Punjab province (Maqsoom et al. 2020). The continuous extraction of groundwater causes salinity issues in the region. In the study area, 85% of the tube wells have been installed between the depth of 30-120 m, out of which 50% electric tube wells are operated at a depth of <45 m and others between 45 and 75 m bore depth (Qureshi et al. 2019).

Drainage pattern
For slope and topography of the study area, an Advanced Land Observation Satellite Digital Elevation Model (ALOS DEM) was downloaded from Alaska Satellite Facility (https://search.asf.alaska.edu/). This model has a spatial resolution of 12.5 m. The drainage network based on Strahler's definition (Strahler 1957) and slope in the area were generated using the hydrological analysis and slope tools, respectively, in ArcGIS 10.5 software ). The ALOS DEM gives an understanding of the topography of the study area and its surroundings, where the elevation from mean sea level ranges from 103 to 240 m with an average elevation of 147 m. The area demonstrates a gentle slope from North to South with inconsiderable variation in elevation, as shown in Figure 2. The drainage network also follows the slope pattern, which gives a clear idea about the surface runoff that flows from the north to the south. The eastern part consists of a canal and mainstream in the study area, which also runs towards the southwestern part. The upper northwestern part consists of fewer streams and no canals or surface streams for irrigation. Therefore, irrigation is mainly dependent on groundwater. Surface information is also necessary to study the depth of underground water. This is also a shred of supporting evidence to study the identification of groundwater aquifer thickness and depth in the upper northwestern part for proper tube well installation.

Vertical electrical sounding
The VES method is the most commonly used method for underground water studies (Asfahani and Abou 2013;Dafalla et al. 2015). In our study, the VES survey data were obtained from the District Agriculture Engineering Department. Data from 16 random locations were collected from the study area. The instrument used for this survey was an ABEM Terrameter SAS 4000. The purpose of using the VES technique was to map the resistivity zones having the presence of good quality and a considerable quantity of groundwater (Dafalla et al. 2015). The instruments give the in-field values of apparent resistivity at the different spacing between the current electrodes (Atat and Ekpo 2018). The apparent resistivity is defined by the current and potential difference, as well as the geometry of the electrodes. The current electrodes are used to inject the current into the ground while potential electrodes measure the voltage (Beck 1981;Niaz et al. 2018;Olaseeni et al. 2019).
The interpretation of resistivity data is based on the underlying geological layers occurring in a sequence of homogeneous, isotropic, and horizontal strata of the finite thickness (Koefoed 1979). Resistivity modelling is used to calculate apparent resistivities at various electrode spacings in terms of true resistivity and underlying layer thicknesses. The widely used Schlumberger array is employed in our study (Olowofela et al. 2005;Szalai et al. 2009). The VES technique is deployed by using the Schlumberger array in which the depth of penetration of current is 20-30% of half of the spacing between the current electrodes (Beck 1981;Telford et al. 2004). However, The depth of investigation specified by Barker (1989) and Kim et al. (2011) is half of the electrode spacing (AB/2). In other words, the depth of investigation can be approximated as AB/4 as suggested by Kim et al. (2011). The maximum spacing between the current electrodes in this study has been kept between 280 and 300 m.
The values of apparent resistivity are used as an input database in an inversion software, IPI2WIN (Dhakate et al. 2015). This software formulates the earth resistivity model in the form of curves. The number of subsurface layers and the true thickness of each layer is automatically solved after inversion in the IPI2WIN software. Based on the true resistivity difference, subsurface layers are differentiated based on lithology type (Table 1).
The resistivity ranges are specified for the respective region where our study area lies. However, these values may differ for different regions by the change in lithology, chemical composition of the rock, structural changes, and water content present. This is the reason to overcome the overlapping phenomenon of the layer's resistivity and their thickness due to suppression and equivalence principles (Singh et al. 2004;Sanuade et al. 2019), thus the DZ parameters are introduced.
The DZ parameters include transverse resistance, longitudinal conductance, and longitudinal resistivity, they can be used to characterise the salinity variation of the area (Yasin et al. 2014;Olaseeni et al. 2019). To approximate these parameters, values of true resistivity and layer thickness were used in various settings (Henriet 1976). The transverse resistance is denoted by T r in Equation (1). The longitudinal conductance is known as the conductance along the bedding plane direction through a 1 m column and is denoted by S C (Parasnis 1979;Nwanko et al. 2011). However, the longitudinal resistivity is defined as the resistivity of rock unit parallel to the plane of the stratification and is denoted by q L (Henriet 1976). The transverse resistance is defined as the total resistance through a l-m column cut perpendicular to the plane (Parasnis 1979;Nwanko et al. 2011). They have been used in computing a distribution of surface potential and the section consists of N fine layers with thicknesses h 1 , h 2 , h 3, :::, h n and resistivity q l , q 2 , q 3 , :::, q n for a block of unit square area and thickness H ¼ P N i¼1 h i (Batayneh 2013). The D-Z parameters, i.e. S C , T r and q L are defined as following: The units for the transverse resistance, longitudinal conductance and longitudinal resistivity are Xm 2 , mho, and Xm, respectively.

IEC
The hydraulic conductivity variation of the Nankana Sahib area is estimated by using the VES and EC data. Using the results of VES a large number of empirical equations that relate electrical resistivity to hydraulic conductivity have been reported in the literature (Niwas and Celik 2012). The hydraulic conductivity values can be used to Very low resistivity zone <15 Xm The admixture of fine sand, clay, and silt study aquifer protection. The clay covering, which has a low permeability, is primarily responsible for the protection of the area's groundwater aquifer. The protective layer's hydraulic resistance is determined by the effective porosity of sandy materials and the clay content of clay materials above the aquifer. The effective porosity and clay content influence the materials' electrical resistivity (Kirsch 2009). The electrical conductivity can be found out by using the following relation (Heaney 2003;Telford et al. 2004).
where q i represents the electrical resistivity. The electrical conductivity is being considered as equivalent to the hydraulic conductivity K i ð Þ (Gemail et al. 2017). The K i depicts the value of hydraulic conductivity which in our case is electrical conductivity value according to above relation. In the above relation, the hydraulic conductivity K i ð Þ is inversely proportional to the resistivity of the layers above the aquifer q i (Chandra et al. 2008;Said Gemail 2012). This value of K i ð Þ can be used to find out the aquifer vulnerability index (AVI).
where h i is the thickness of the layers above the aquifer and C is AVI. According to (R€ ottger et al. 2005), the AVI is called as IEC which can be derived as The estimated IEC can be termed in ohm À1 X À1 ð Þ or siemens (S).

Electrical conductivity and total dissolved solids measurement
The VES result cannot be used alone to interpret specific lithology for the exploration of groundwater because of the varying nature of subsurface material upon the applied current. The different material depicts different resistivity against the injected current. For comparison, the EC data of water samples were also obtained from the District Agriculture Engineering Department. The EC was measured by using a portable EC meter (Sugirtharan et al. 2015). The Electrical Conductivity/TDSs Meter (Hanna Model) was used for this purpose (Arnold et al. 2005). The water samples were shaken thoroughly, and the electrode of the instrument was first washed with distilled water and then with the sample water. With the standard solution of potassium chloride 0.01 m at 25 C, the conductivity of the instrument was standardized, and then the electrode was immersed in the sample water, and EC readings were recorded. This process was repeated three times to get better accuracy (Akaahan et al. 2017). The TDS can be defined as the total amount of dissolved solids in the water, while salinity is the proportion of salts in the water (Hinze 1982). The TDSs and electrical conductivity relationship for groundwater can be estimated with the help of Equations (4) and (5) The electrical conductivity of the water samples has a direct relation with the concentration of dissolved solids in water. Equation (7) is used if the value of EC of the water sample is <5 ds/m or 5000 lS/cm (Fipps 2003). The results show 10% accuracy when correlated with laboratory TDS measurement as explained (Dahaan et al. 2016).
The data for the current work was obtained from the District Agriculture Department; they only provided the VES and EC data for the groundwater salinity studies.

Interpolated maps of area
The interpolated maps of true resistivity and their thickness, transverse resistance, longitudinal conductance, and longitudinal resistivity were generated by using GIS to mark the different salinity zones in the area. All the DZ parameters were calculated from the VES points at the same spatial locations (grid data). The inverse weighted distance (IDW) technique, a tool in ArcGIS10.5 software, was used in this paper to interpolate the VES grid data all over the surface in the study area. The IDW is a deterministic technique that estimates a value for any unsampled location based on the basic assumption that things close to each other are more similar than things farther apart. The IDW method interpolated the surface from the known VES points all over the study area. The IDW approach has been utilized for VES by many other authors Prasanna et al. 2019;Arunbose et al. 2021). We kept the cell size of 12.5 m in the interpolation process, which was similar to the spatial resolution of the ALOS DEM that we used earlier in the analysis. We chose measurement points at random in areas where local farmers had reported high salinity in the ground water. The minimum distance between the points is $329 m between VES 10 and VES 11, while the maximum distance between the points is $29,525 m between VES 2 and VES 4. Our goal was to create a report that defined a water aquifer depth suitable for drinking and irrigation in the study area.

Vertical electrical sounding
The values of true resistivity of the subsurface layers fall between 0.403 and 772 Xm as output from the IPI2Win software. The curves in Figure 3 explain the relationship between the apparent resistivity and the electrode spacing. The change of depth gives the variation in apparent resistivity of subsurface layers. Figure 3 shows the type of curves for 4, 5, and 6 layers, respectively.
In IPI2WIN software, with the help of the process of iteration, the black line is kept approximately equal to the red line by keeping the RMS error between 2% and 3% (Utom et al. 2012;Ur Rashid et al. 2017). The curve in Figure 3(a) is a 4 layers curve in which q1 < q2 > q3 < q4, where the high resistive layer is entrapped between the two low resistivity layer and the 4th layer resistivity is greater than third layer resistivity. The curve in Figure 3(b) is a five-layer curve in which q1 > q2 < q3 > q4 < q5, where the second layer is entrapped between the two high resistive layers and fourth layer resistivity, is less than the fifth layer resistivity and the third layer resistivity. The curve in Figure 3(c) is a six-layer curve in which q1 > q2 < q3 > q4 > q5 < q6, where the second layer is entrapped between the two high resistive layers and the fourth layer has less resistivity than third layer and fifth layer is entrapped between two high resistivity layers as sixth layers resistivity is greater than the fifth layer. The results obtained from these curves are used to calculate maps for the true resistivity subsurface variation and the DZ parameters.

The true resistivity of subsurface layers
The interpolated maps of true resistivity data obtained from the VES results help to identify the interpolated extent of true resistivity for different layers concerning different thicknesses and VES locations. The VES layers of true resistivity values and their thickness true resistivity variation are shown in Figure 4. The criteria for the fitness of groundwater for irrigation purposes are given in Table 2. The surface layer is the aggregate of the first two layers obtained by using raster calculator tool in arc map (Gemail et al. 2017). The area of Rachna Doab on an average has water table depth of Figure 3. Some types of VES curves having 4-, 5-, and 6-layers, In the above curves, the current electrode spacing AB/2 (m) is on the X-axis and true resistivity is denoted by qa (Xm) values on the Y-axis. Red and the Blue curve gives information about the relation between AB/2 and true resistivity value. The blue line gives information about resistivity value variation while the red line is the master curve. more than 18 meters . The aggregate of the first and second layer which is termed as surface layer appears to be less than 18 m except two points which are VES 7 and VES 8 as in Figure 4(c). However, the aggregates of the depths of the third-and fourthlayers ranges between 26 and 208 m. This is the reason, the aquifer layer is considered as the aggregate of the third and fourth layer, respectively, for this study.
The resistivity of the surface layer ranges between 4.42 and 173 Xm and layers thickness ranging between 0.5 m and 27 m as indicated by Figure 4(a). This variation in the resistivity is caused by the surface inhomogeneity, precipitation, surface runoff, weathering, and recharge conditions. The lowest values of true resistivity in surface layers can be observed in the eastern area while the highest values can be seen in the western area as in Figure 4(a). The high values of the true resistivity in these areas indicates good ground water zone but due to small aquifer thickness the tubewell is never recommended. The resistivity of the aquifer layer shows the high values in the central, north-central and NW areas as in Figure 4(c). Almost 60% of the area can be termed as marginally fit groundwater zone to fit groundwater zone with greater aquifer thickness ranging between 60 and 95 m mainly in the central area. The greater aquifer thickness probably gives more yield and longer-term water supply. The NE and NW areas show the lowest values indicating it as an unfit groundwater zone for irrigation as criteria defined in Table 2 by WAPDA (1978) and Sikandar et al. (2017). The lowest thickness of the aquifer is observed in the NW central area of aquifer layer as in Figure 4(c). Because of the low thickness of the aquifer, it cannot be recommended as the best depth for installing tube wells or dug wells because our aim is to sustain long-term water-supply (Shah et al. 2021).

The Dar Zarrouk (DZ) parameter
The DZ parameters, including longitudinal conductance, transverse resistance, and longitudinal resistivity, are presented in Figure 5. Based on these parameters, different aquifer zones can be given by the specific values (Hasan et al. 2020). The freshwater aquifers were marked as having T r > 1700 Xm 2 , q L > 25 Xm, and S c < 2.6 mhos. The saline water aquifer was marked as having T r < 700 Xm 2 , q L < 15 Xm, and S c >3.3 mhos. While brackish water aquifer was marked having T r from 700 to 1700 Xm 2 , q L between 15 and 25 Xm, and S c between 2.6 and 3.3 mho.
As shown in Figure 5(a), the longitudinal resistivity decreases from the NE and SE to the NW and SW. The longitudinal conductance increases from NW to SW to SE, as shown in Figure 5 (b). As shown in Figure 5(c), the transverse resistance decreases from the SE and NE to the NW and SW, indicating an increase in salinity. While on the NW side, the transverse resistance is high at one point and low at another, the rest of the area represents low values. The same trend can be seen in Figure 5(a-c), where salinity increases from the NE and SE towards the NW and SW.

Electrical conductivity and total dissolved solids
The EC variation in the area is depicted in a map with their specific depth in Figure  6. The salinity is increased by moving from the SE towards the NW. The NE, SE, and central areas have better quality water when compared with the remaining area. The suitability standards for the EC and TDS of water samples for irrigation are deduced from (Robert and Westcot 1985;Kwiatkowski and Pittman 1997) as in Table  3. In this table the effect of water on crops have been defined in three general terms as none, moderate and severe effect on the basis of EC and TDS. Table 4 illustrates the EC and TDS values at VES points, the mentioned depth is the depth of the bore for the previously installed tube well from where the water sample has been collected. Table 2. Irrigation water fitness classification (FFI) based on groundwater salinities (EC) derived from aquifer resistivity measured by VES based on criteria (WAPDA 1978) and (Sikandar et al. 2017 The groundwater quality at VES 8, VES 13, VES 14, VES 15 are good for agriculture and drinking purposes as well because of the low yield of EC as shown in Table 4 at these VES locations (WHO 2008). Based on the EC distribution in the area as shown in Figure 6, an increase in EC can be observed from the SE towards the NW, and three major zones can be observed: freshwater in the SE and Southcentral region, brackish water in the central region, and saline water in the NW region.

Integrated electrical conductivity (IEC) variation in the area
The IEC results serve the purpose to study the vulnerability of the aquifers towards contamination. The low IEC values correspond to highly vulnerable zone which provides poor protection against contaminants and pollutants. The low IEC zone reflects the vertical migration of wastewater from the sandy soils inside the clay cap. Fine sand in this zone may be used as windows for leaking of surface waste water from both the surface drains and the fertilizers used and thus increase levels of contamination in underground water (Thorling et al. 1997). The IEC distribution in Nankana Sahib is presented in Figure 7.
From Figure 7 we can see that the highest values for IEC are observed in the Southcentral and southwestern part of the study area. It indicates good protection from leakage and pollutants. However, the rest of the area in Figure 7 shows relatively low values of IEC usually <1. The lowest values of the IEC are indicators of the highest vulnerability to pollutants and contaminants.

Overlay analysis
The parameters including true resistivity, electrical conductivity, longitudinal conductance, transverse resistance, and longitudinal resistivity were reclassified into three classes. The reclassification classes were based on the division of parameters according to ground water quality for irrigation e.g. fresh, brackish, and saline water. An integrated map was prepared using overlay analysis where all parameters were assigned equal weights. The output map consists of three classes that defines the overall ground water quality for irrigation purposes i.e. fit, marginally fit, and unfit. The overall results in the overlay analysis represents the three general zones extend in the area as in Figure 8a. The major part of the SE and NE area can be termed as fit water zone with a marginally fit water presence in the top corner of NE side. The central part of the area consists mainly of marginally fit water while the western side can be termed as the unfit water zone. The depth of the aquifers as in Figure 8b ranges between 21 to 138 m. The fresh water zone and brackish water zone which lies in the NE and north central area has aquifer depth between 30 and 60 m which can be termed as an ideal depth for the installation tube well for getting fit groundwater. However, for south central area the depth of the aquifer should be ranging between 60 and 90 m. This part of the area majorly consists of fit ground water zone and can be termed as the fresh water zone.

Discussion
Pakistan has the world's largest canal irrigation system according to the Indus, Jhelum, Chenab, Ravi, and Sutluj Rivers (WAPDA 1980). The area lying between the rivers is known as Doab, and our study area lies in Rechna Doab. It is situated between the river Ravi and Chenab in Punjab. Punjab is covered by 200 m of alluvium from the Himalayas, which is underlain by Precambrian basement rock (Farooq et al. 2017). The unconfined aquifer of Punjab has an alluvial plain with more than 350 m thick Holocene and Pleistocene sediments transported by the Indus and its tributaries (France-Lanord et al. 2008) that is directly overlain by Precambrian   The homogeneity of lithology is disturbed by the dispersed foundation in the central part of Doab. Bedrock represents a strong decline at the upper doab where no bedrock was found in a well drilled at Sheikhupura to the depth of 460 m (Aslam 1997). The Rechna Doab is divided into two groundwater basins by the buried ridges. Therefore, aquifers are highly porous, they can efficiently store and transmit water. Lenses with thick clay are not predominant in the area (Anjum et al. 2016). A VES survey was conducted to determine the area of saline water and freshwater and vulnerability classification. The drainage pattern study of the area was also done. Area analysis of the drainage pattern ( Figure 2) gives information about the surface runoff that runs from north to south. It also directs us to mark zones with fewer stream channels and irrigation sources, which could indicate the saline water zone. The modelled VES curves were tuned and interpreted using knowledge from the study area's geology and lithology as in Table 1. Depending on moist conditions, sand resistivity is often greater than clay resistivity; similarly, gravel resistivity is greater than sand resistivity (Akhter and Hasan 2016). The areas having values corresponding to red rectangle values can be termed as a good groundwater potential zone having medium to coarse sand with thin silty clay layers (Sikandar et al. 2010).
The NW and SW area is a saline water zone, which shows low true resistivity values. The results from the VES and EC data show that the VES points VES 2, VES 6, VES 9, VES 10, VES 11, VES 12, and VES 7 in the NW part of the study area are in the saline water zone. The groundwater quality at VES 8, VES 13, VES 14, and VES 15 are good for agriculture and drinking purposes as well because of the low yield of EC (Table 4) at these points. The EC values aids for the determination of the quality of water samples, the high EC values greater than 2500 lS/cm tends to have a severe effect on the growth and production of the crop, while the low values <1000 lS/cm appears to have no bad effect on the crops (Ayers & Westcot 1985;Kwiatkowski and Pittman 1997). However, a lot of work has been done to find the suitability of the acceptable EC of the water. The EC should never exceed 750 lS/cm for the seeding process, but for older grown-up crops, it could be accepted up to 1500 lS/cm (Bauder et al. 2011). The TDS values suitability for different types of crops is stated in Table 3.
For some crops, fresh water should be used because of their sensitive nature, but in the areas of Nankana Sahib, the most dominant crops are wheat and rice (Riaz and Javaid 2012). And those are slight to moderate with the salinity effect of water (Bauder et al. 2011). The TDS values obtained by Equation (7) cannot explain the exact and accurate details of salt intrusions in the water, but it provides the approximate amount of TDS in the water within the accuracy of 10% with the laboratory TDS measurement (Dahaan et al. 2016;Akaahan et al. 2017).
By the comparison of the EC and VES results, it can be deduced that by the increase in depth the true resistivity is increasing in the SE and NE areas. There is a freshwater aquifer having a thickness between the ranges 60 and 90 m in central areas. The NE, and SE areas in Figure 3(c,f) are saline water zones. The salinity in these areas may be caused by the excessive usage of groundwater for irrigation purposes. There is no recharge source in nearby of these areas. In addition, rainwater is the only way to recharge the aquifers in this area. The pumping of groundwater through tube wells in this area is high as agriculture mainly depends upon groundwater in this area, which is not only deepening the depth of the top of the aquifer but is also making it more saline. The rest of the area can be marked as a marginally fit groundwater zone for irrigation. Moreover, the result from the IEC (Figure 7) shows the highest vulnerability for almost 70% of the area, while the central and SW areas show the highest values for the IEC indicating good protecting capacity against pollutants and contaminants. The maximum area falls in the lowest values of the IEC indicating weak protection against contamination. According to the IEC results, the areas currently having marginally fit water are also in danger of getting contaminated in the future. These results indicate that in general, almost all of the areas are on the verge to get contaminated except south-central area and SW area. However, by following proper strategy contamination can be avoided. The installation of waste water treatment plants is highly recommended for this area.
The results from the apparent resistivity, EC, and DZ parameters match each other and have shown the same trend. The trend can also be observed in the overlay analysis where the study area can be differentiated into three general groundwater quality zone i.e. fit, marginally fit and unfit. All of them indicate that the salinity of underground water is increasing from the SE towards the NW. The average depth for installing tube wells in the southern area has to be ranging between 30 and 50 m. In the NE and NW areas installing tubes, wells are not recommended. In the NW area, there is no indication of better-quality water even at a deeper depth. However, by drilling a deeper depth tube wellbore, better quality water can be extracted. In the central region, the tube well can be installed between the depths of 60 and 90 m. While in the SE region tube well is recommended between 30 and 60 m approximately to get better yield and marginally fit the water. Due to the limitation of the accessible data, the results of the VES and EC are not correlated with borehole data. This is the reason the detailed discussion on lithology has not been made in this paper, as lithology can only be confirmed by the comparison of VES interpretation and borehole data. Instead, the criteria set by WAPDA (1978) and Sikandar et al. (2017) in neighbouring areas are considered for the interpretations and recommendations for this study.

Conclusions
The findings from the VES and EC show resemblance to each other in the study area. Through the help of interpolated maps, the study area has been divided into three major vicinities using the EC, VES, and IEC results. The SE and central areas are marked as the fresh groundwater zone having the maximum thickness aquifers. The NE area is marked as a brackish groundwater zone, while the NW part is marked as the saline water zone. The SW part contains less saline water while the salinity increases, moving towards the northwestern. The northwestern part of the study area is interpreted to be highly saline, and so no tube wells or dug wells are recommended in that area. On the NE side, between the depths of 30 and 60 m, tube wells can be installed, while in the central and SE areas an average depth for installing tube well should be between 60 and 90 m. The VES and EC techniques are the most common and economic ones to apply in the field and get approximately accurate results. The current study will act as a guide for farmers of the selected area for the installation of tube wells. Moreover, this study is also applicable for similar areas of Punjab, Pakistan, having salinity issues reportedly in the groundwater. However, better accuracy can be achieved through the comparison of the VES results with borehole data. An extensive borehole study in the area is highly recommended to further explore the subsurface lithology and provide comparable results to the VES and EC data.