Received 22 March 2016; accepted 23 May 2016; published 27 May 2016
The greater Punjab plain which is divided into five Doabs (land between two rivers) is suffering from agricultural stresses, because of rapid growth of agriculture and population. Groundwater is major source of water for growing agriculture and population in Pakistan. Groundwater is present in voids and it fully saturates it; porosity and permeability of that stratum have greater importance in groundwater recovery and protection  . Along with production, the groundwater protection is also very important, particularly in arid regions of the world. Agricultural Non-point Pollution Sources (NPS) and urban wastes are the major threat to the portability of groundwater in Punjab, Pakistan. An assessment of these NPSs on groundwater is an important environmental concern.
Several methods of aquifer vulnerability have been proposed; among these methods, index methods are the most popular. There are different such models but the most useful is DRASTIC model. Scientific history of the world has witnessed many examples of development and applications of index methods like DRASTIC, GOD, IRIH, pollutants  . Aquifer vulnerability assessment by S-model is a cheaper alternative to index models. In S- model the protective capacity of vadose zone is assessed by determining lithological composition and thickness through electrical resistivity sounding. Longitudinal conductance (S) values are grouped into different vulnerability zones. The resistivities of earth material depend on soil type, porosity and fluid present in pore spaces and can be determined by electrical resistivity surveying, a geophysical technique used for the determination of such resistivities variations  . The use of geophysical methods in groundwater studies has improved greatly because of development in computers  . Electrical resistivity sounding is economical and highly suitable for hydrological, geotechnical, mining and environmental studies. There are different geophysical methods used for groundwater development through determination of formation resistivity, properties of a fluid inside pores and protective capacities of material above the saturated zone   . VES has proved a non-distractive and appropriate method for determination of various subsurface geologic strata and groundwater quality. This technique, VES with Schlumberger electrode configuration is very famous for the groundwater studies in arid regions  . ERS has been used in various groundwater studies for the achievement of multiple objectives in numerous hydrogeological conditions  . Some of important studies are, groundwater quality  , delineation of fresh and saline water boundary  , for imaging the industrial contamination plume  , determination of aquifer parameters   and hydraulic conductivity and transmissivity  .
The present study aimed at determining the pollution potential of the unsaturated layer present above the aquifer layer. Electrical resistivity method is utilized by adopting Schlumberger array. Dar Zarrouk parameter’s longitudinal conductance (S) values, which reflect the permeability of unsaturated layer, provide the vulnerability of aquifer to the surface contaminates.
2. Material and Methods
2.1. Study Area
Study area (Khezarabad Diary Form) lies in district Sargodha of Punjab at geographic coordinates Latitude 32˚14'3.01"N and Longitude 72˚45'30.76"E (Figure 1). District Sargodha is part of lower portion of Chaj Doabs, alluvium between Jhelum and Chenab River. There are two crop seasons as Rabbi and Kharif. The region is flat with well-developed canal system. Area is irrigated by tributaries from Jhelum River and motor driven pumps. This city is renowned for citrus production world-wide as climate conditions are highly suitable for its cultivation. Khizarabad livestock dairy farm has been selected as case study to assess protective capacity assessment of vados material zone. Case study research allows research to test hypothesis empirically. Availability of canal water is perennial study area consists of loamy soil (http://dlfpunjab.webs.com/khizerabad.htm).
The climate of the study area is arid as rainfall is less than evapotranspiration rate. The annual mean rainfall is 426 mm. Mean annual temperature is 24˚C with the mean maximum of 31˚C and mean minimum of 18˚C. Mean Relative Humidity at day time is 71% whereas the night time it is 79%. Study area is characterized by Monsoon showers during July and August contributing significantly in annual precipitation (Figure 2). In October and November, very little precipitation is recorded at this station.
Figure 1. Regional location of study area on Pakistan map with digital elevation model, major roads, rivers, canals, electrical sounding points and tube wells locations.
Figure 2. Average monthly rainfall recorded at Sargodha station for period 2001-2010.
2.3. Geology and Physiography
Geologically study area consists of thick alluvial deposits of Quaternary age over basement rocks. Major geological units are active flood plain and abandoned flood plain. There are also basement rocks exposed to surface  . The alluvium complex of Quaternary age deposited on metamorphic and igneous rocks of Precambrian age in Chaj Doab. The Alluvium consists of fine to medium sand with slit and minor clay lens. The strata lack both horizontal and vertical continuity. There are three important geological units as Potwar Losses, Piedment Deposits and Aeolian Deposits  . The underlain aquifer has good hydrogeological characteristics such as hydraulic conductivity and Transmissivity. Another important factor that influences hydrogeology of region is gradual exposure of bedrocks from central operation of the Doab to Kirana where bed rocks outcrop (Figure 3). This ridge influences greatly the groundwater flow everywhere in the Doab.
2.4. Concept of Vulnerability
It is dimensionless index which is considered as the function of hydrogeological factors, sources of contaminants and anthropogenic effects in a specific site  . The concept of groundwater vulnerability to contamination is based on the assumption that the physical environment may provide some degree of protection to ground- water against natural and human impacts with respect to contaminants in the groundwater. The vulnerability of a certain area can be described by the degree of susceptibility of that area to groundwater pollution  . In 1968, the French Margate was the first one who used the term vulnerability in Hydrogeology, thereafter, the concept was adopted worldwide  . Groundwater vulnerability to contamination is defined as “The sensitivity of groundwater quality to an imposed contaminant load, which is determined by the intrinsic characteristics of the aquifer”  . Groundwater vulnerability to contamination is the tendency or likelihood of contaminants to reach a specified position in the groundwater system after introduction at some location above the uppermost aquifer. As can be inferred from the above definition, groundwater vulnerability is not an absolute or measurable property, but an indication of the relative possibility with which contamination of groundwater resources will occur. This understanding implies a very basic vulnerability concept that all groundwater is vulnerable  . In the pre- sent study vulnerability is taken equal to protective capacity.
2.5. Resistivity Method
2.5.1. Electrical Resistivity Survey
Apparent resistivities in the field are measured by using SAS-300 equipment by GeoWatAssociates, Islamabad (PVT), a consulting company in July 2011  . Vertical Electrical Soundings provide information of different geological section and vertical variations of different water quality zones  . In this survey Schlumberger array is used because of its data acquisition and interpretations are less labor intensive and economical  . In this
Figure 3. Physiographic units of Chaj Doab of which study area a small portion modified from (Hussain et al., 1990)  .
array the central potential electrodes are fixed and the other peripheral current electrodes are moved horizontally in order to target different depths. The two outer electrodes A and B are used for the current, and the resulting potential difference is measured across the two inner electrodes M and N. MN/2 is always kept sufficiently small relative to AB/2. The distance from center of current electrodes is high comparative to potential electrodes (Figure 4). Apparent or field resistivity is calculated by using Equation (1).
where ρa is apparent resistivity of earth material measured in ohm-meter, I is current applied in milliamp ere (mA) and K is geometric constant, its value depends on the electrode configurations. For Schlumberger electrode configuration its value can be calculated as
where AB is current electrode separation and MN is potential electrode separation. Resistivity values of dry rocks are very high while presence of clay and salted water in the pores reduces the resistivity significantly  .
Resistivity data can be interpreted both qualitatively and quantitatively. Qualitative interpretation can be performed by apparent resistivity curves plotted between filed resistivity values and different AB/2 values. This will provide an initial guess about the variation of resistivity in field. Qualitative interpretation of resistivity data is based on two basic assumptions that earth is horizontal and resistivity of last layer has infinite thickness  . Each layer in 1D model is assumed as uniform and isotropic.
Primary data about VES locations was collected from report  . Resistivity values are grouped into three zones as upper unsaturated, lower high resistivity and low resistivity based on the modeling results of resistivity models
IX1D is used for the interpretation of resistivity curves in terms of layered resistivity model by matching measured curves with theoretical master curves using master curve matching techniques  . These models provide resistivity and thickness of each layer (Figure 5). Geoelectrical section from the results of ID models is obtained along profile AB (Figure 6). These layered parameters are used for the calculation of Longitudinal Conductance (S) and resistivity of unsaturated zone. The Dar Zarrouk (DZ) Parameter (Longitudinal Conductance S) is calculated by layered parameters obtained from ID resistivity curve modeling. Following equation is used for the calculation of DZ parameter.
where S is used in for the assessment of protective capacity of aquifer, hi is the thickness of unsaturated layer and is the resistivity of that layer
2.6. Preparation and Generation Thematic Map
2.6.1. Unsaturated Thickness Map
Isopack thematic map is created by using IWD interpolation algorithm of ArcGIS 10 software. IWD works on the principle that value at unmeasured point can be predicted from neighboring measured points by assigning a weight factor, which varies inversely with the distance from measured point  .
Figure 4. Sketch diagram of Schlumberger configuration array, AB are current electrodes and MN are potential electrodes.
Figure 5. ID resistivity model of VES-1 created in IX1D resistivity software.
Figure 6. Sketchof geoelectrical section along profile AB.
2.6.2. Layer Resistivity Map
The first layer resistivity is interpolated. The resultant map gives resistivity of first layer at different locations in the study area which is then used in aquifer protective capacity assessment. According to Braga et al. (2014), lithology of unsaturated zone is difficult to map with ERS, however, litho-facies of saturated zone can be determined with this method. Study is underlain by unconfined alluvial sandy aquifer and in this case lithology of saturated and unsaturated can be considered as same  . The lithologies of unsaturated zone are built by extrapolation of lithology of saturated zone   .
2.6.3. Longitudinal Conductance Map
Earth martial above water table acts as natural filter against contaminants passing through it. The vadose zone material lithology in term of permeability and porosity is very important in determining pollution potential of underlying aquifer. Such vulnerability potential is related directly with longitudinal conductance S for clayey aquifer, that can be used to determine protective capacity which is taken as infiltration time of contaminates. Here PC is related to S which is measured in units of meter ohm  .
Above layer parameters are used for the calculation of protective capacity as given in Equation (3). The input parameters as well as the protective capacity is evaluated in GIS environment by using ArcGIS 10.1. The algorithms utilized by current investigation is IWD present in Spatial Analyst Tools Module of ArcGIS 10.1 software. Raster Calculator is used for the protective capacity evaluation (Equation (3)). Other major tools are Raster to Polygon conversion, Geo-referencing, reclassify and Geometry calculation.
3. Results and Discussions
3.1. Thickness and Resistivity of Unsaturated Zone
The thickness of dry strata over aquifer varies from 14 (minimum value) to 57 meters depth (maximum value). The north western portion has shallow groundwater depth increases vulnerability of surface contaminates to reach aquifer by reducing their travel times. The south eastern and south western parts have very deep groundwater and provide greater protection to underlying aquifer (Figure 7). Higher the resistivity of material of vadose zone lower is the protective capacity (Equation (3)). Resistivity of first geoelectrical layer varies from less than 12 to greater than 15 Ωm (Figure 8). Major portion of area is occupied with resistivity range of 13 - 14 Ωm.
3.2. Protective Capacity Mapping
The aquifer protective capacity or vulnerability map is classified into three classes using quartile classification of ArcGIS 10.1. The resultant vulnerability is divided into three classes as low (>3.8 siemens), moderate (1.6 - 2.5 siemens) and high (<1.6 siemens). Major portion of study area is under moderate vulnerability zone (Figure 9). The high vulnerability is the reflection of shallow water table and low apparent resistivity of material overlie
Figure 7. Spatial distribution of thickness of unsaturated zone.
Figure 8. Spatial distribution of resistivity of 1st geoelectrical layer.
Figure 9. Protective capacity of unsaturated zone divided into different vulnerability zones.
the aquifer while reverse is true for low vulnerability zones.
This study was carried out to provide vulnerability assessment of unsaturated layer by using electrical resistivity surveying technique of geophysical method. Longitudinal conductance (S) is used for the calculation of relative aquifer protection from the surface contaminants. The resultant vulnerability capacity is divided into three zones high, moderate and low. The results indicate that the most part of study is under low to moderate vulnerability while a small portion is under high vulnerability. The areas having high vulnerability should be taken under consideration of Khezarabad Diary Form administration.
Authors would like to thank GeoWat Associates, Islamabad (PVT) for providing necessary data for this study. We would also like to thank Mis. Farzana Ahmad Deen, who made English corrections.