3D Interpretation of Resistivity Data for Groundwater Potential Assessment of Pakhli Plain, Mansehra District, Pakistan

Abstract The present research describes a method of combining geostatistical analysis with geophysical inversion of electrical resistivity data conducted in Pakhli Plain, northwestern Himalayas, Pakistan. The raw data has been collected from the Technical Report VII-I on Ground Water Resources in Pakhli Plain, Mansehra District. Subsequently, the data has been deciphered and broadened from one dimensional resistivity data into a 2D model that can be entirely visualized and deduced in a spatial sense. Interpretation and calibration of the electrical resistivity curves with the lithologies and geophysical logs of boreholes suggests possible identification of distinctive sedimentary accumulations occurring within the Pakhli Plain. The 2D and 3D gridding and visualization is imperative to map the extents of the alluvial deposits within the Pakhli Plain formed during the periods of extreme tectonic activity. The coarser sediments are associated with lower levels of resistivity as measured in the electrical surveys, whereas the finer sediments exhibit characteristically lower resistivities. Therefore, the zones of low and high resistivity values are indicative of particles associated with coarser and finer sediments, respectively. It has been mentioned that the Pakhli Plain has remained a lacustrine zone during some time in the geological past as indicated by low resistivities representing finer sediments in the middle of the Plain. Consequently, the overall transmissivity of the sediments is low, which imply poor conditions for commercial groundwater production in the Pakhli Plain. Moreover, high resistivity zones of coarse material could be further investigated for groundwater potential areas. In particular, the prime objectives of the present study include 3D modeling of underground resistivity and its exploration in terms of groundwater potential on the basis of distribution of low resistivity zones.


Introduction:
The diversity of climate and fertile land of Indus Plains, a large variety of crops is grown to support the agricultural economy of Pakistan. The average annual rainfall of the country is about 250 mm. Agriculture under such climate requires irrigation to meet crops water requirement. Irrigation in the country is dependent on conjunctive use of surface water of Indus Basin Irrigation System (IBIS) which delivers 172 BCM and groundwater of underlying Indus Plain Aquifer (IPA). A huge percent of groundwater resources of Pakistan occurs in alluvium aquifer under Indus Plain. Pakhli area is an extensive part of Mansehra District exposed in Khyber PakhtunKhwa province, Pakistan. It is located on both sides of Karakorum highway from Mansehra city approximately 3 kilometers north of Shinkiari Village, between 34o-05' N and 73o-20' E. Geologically, Pakhli Plain is widespread portion of southern flank of the northern Hazara highlands which are reflected as western extension of the Himalayan chain ( Figure 1a). Lithologically, this area is chiefly comprised of metamorphic rocks of the Tanawal Formation as well as Mansehra batholith with special prominence of granites. The Mansehra Granite forms bedrock of alluvial fill present in southeastern part of the plain and the remaining area is underlain by metamorphosed bodies specifically schists and quartzites.
Groundwater in the Pakhli Plain primarily depends upon water table conditions (Figure 1b) where depth to water is about 20 m in most part of the plain. Additionally, groundwater movement of the described area is consistent with surface drainage. The Bhut Katha, Ichhar Nallah and Siran River are targeted as draining points of sub-surface water. Approximately 51 x 106 m3 is the estimated total annual recharge, whereas approximately 1.0 x 106 m3 groundwater usages is restricted for domestic purpose. Around 884 x 106 m3 groundwater stored up to bedrock has been projected by Sheikh, [1]. The exploration of groundwater through resistivity approaches is measured as effective and meaningful. Several researchers [2][3][4] inferred the subsurface hydrogeological and geomorphological features in addition to bedrock position by using this analytical scheme. Clay content in fresh water aquifer systems has a robust relationship with electrical resistivity phenomenon. These contents of clays are considered as good electrical conductors whereas coarser sediments gradually behave as poor conductors [5]. Besides these, subsurface lithology reveals important relation to resistivity data, while the results are generally predicted non-unique and of low-resolution comparative to well logging. The strong correlation between resistivity and aquifer characteristics (porosity, permeability, hydraulic conductivity and transmissivity) [6] is challenging to establish in contrast to correlate with lithology. The current research work highlights the groundwater potential of sedimentary accumulations of the Pakhli Plain based on clarification of electrical resistivity curves and calibration with the lithologies and geophysical boreholes logging.

Materials and Methods
Resistivity method is employed in the present study for determination of subsurface hydrogeological and geomorphological features in addition to bedrock position [2][3][4]. To achieve the objectives of the proposed research, a systematic study carried out using multiple analytical techniques. The groundwater exploration of any area is complex work and it needs integration of conventional hydrogeological and geophysical surveys. The application of techniques depends upon scope and purpose of the investigations of availability of fresh water. In the present study the task of groundwater investigations and evaluation was accomplished in different stages of surveys and investigations; literature review, application of surface geophysical methods, exploratory well drilling and sampling, quantitative aquifer evaluation, water quality evolution for utilization of explore groundwater.
In resistivity method, the distribution of electrical potential in the ground around a charged electrode mainly depends on the resistivity and division of the surrounding soils and rocks. This approach is based on direct current (DC) which is applied between two charged electrodes mounted in the ground and computes the difference of electrical potential between the other two additional electrodes that do not transmit current. Furthermore, the potential electrodes are placed in line between the current electrodes and principally these are flexible to relocate anywhere. Finally, the outcome of resistivity survey data is in the form of apparent resistivity values after interpretation. Apparent resistivity ρa, is an Ohm's-law ratio of measured voltage V to applied current l, multiplied by a geometric constant k which depends on the electrode array ( Figure 2).

Figure 2:
Equipotential lines and bowl-shaped current lines for current electrodes A and B on a uniform half-space.

Integrated Geophysical Surveys
The geophysical methods are the best scientific tools to use for the wide coverage and give rapid results in groundwater investigations. Each geophysical method has limitation and advantage, therefore integrated approach of using several geophysical methods was applied in the Pakhli Plain. The research area was investigated by applying electrical resistivity and induced Polarization surveys in field. The electrical resistivity is most extensively used in groundwater investigations as it gives the signal of water quality in calculation to the lithology of the concerned area. Vertical Electrical Soundings (VES) were conducted at regular grid for systematic and useful investigating evidences of the groundwater. Some factors carefully a grid of 5 km x 5 km spacing was implemented in the project area for the investigated depth of 300 m. In adding together to this a relatively coarse grid of 20 km x 20 km was implemented for investigated depth of 1000 meters / bed rock for mapping of thickness of alluvium aquifer. In this way, ABEM SAS 4000 Terameter was used for VES and IP measurement. The resistivity modeling was accomplished by using IX1D software.
The data was attained from the Technical Report VII-I on Ground Water Resources in Pakhli Plain Mansehra District presented by Sheikh [1]; it was then interpreted and broadened from one-dimensional resistivity data into 2D model and 3D models. The vertical electric resistivity sounding method is based on the variation in subsurface resistivities with depth. Classically, subsurface resistivities are termed as formation resistivities (ρf) which is calculated as the product of Formation Factor (F) and the resistivity of the groundwater (ρw). Formation resistivity can be typified in the form of following equation: Depending upon the physical properties of the formation like porosity, permeability clay content and degree of cementation, the formation factor (F) can be evaluated. More importantly, as the clay content increases, the formation factor declines however, it rises if porosity drops with increasing the degree of cementation. Commonly, for the loose sand, formation factor is almost calculated in the range from 3.5 to 7 which is considered higher for gravel and boulder-type lithologies. The formation factor for sands mixed with clay and silt may be around 2 owing to exertion of reducing effect. From the water samples taken from boreholes or wells, the resistivity of groundwater (ρw) is determined through following equation: This indicates that the value of electrical conductivity (EC) is the reciprocal of groundwater resistivity. Where ρw is calculated in ohm.m and electric conductivity (EC) is measured in μs/cm (micro-siemens per centimeter). Sometimes resultant curves produced from resistivity data may need to match with several realistic resistivity distributions since most geophysical methods generate results which can be tentative to interpret. Therefore, in order to finalize a geological model, it is necessary to calibrate the resistivity data against another individually derived data set. So, in present case, we observed & measured lithological logs, surface geology, water levels and electrical conductivities [1]. Using the IPI2WIN software, 2000 [7][8][9][10][11][12] it is conceivable to attain the inverse modeling of the resistivity after calibration with borehole lithologies. Ω. Figure 3 indicates the Borehole information W-4 location along with VES point No 46, W-11 nearest to VES point No 21 and W-12 is close to VES point No 33, in the study area. L/2 max in VES = 600, [1] .

Pump Test Analysis of the Pakhli Plain
Pump test is significant to measure water level in water wells as a function of time. In pumping test, pumping of water from a well is important to calculate discharge of the well with time. Generally, this measurement is accomplished at a constant or regular interval. Figure 5 reveals the location of pump test in the study area. The drawdown in the well is deliberated through piezometers either in the same well or in other wells at known distances from the well. Hydraulic properties of the aquifer such as Storability (S), Hydraulic conductivity (K) and Transmissivity (T) can be assessed by substitution of these measurements in an applicable way in well flow equation. The cone of depression is formed around the well when a well is pumped with a constant discharge (Q) and this cone can be perceived in pumping well along with observation wells [13]. Pump tests are conducted in order to analyze the flow of groundwater towards the water well. During pump test, the well is pumped with a precise mode at determined rates and the subsequent effects on water levels are restrained in both the pumping well and observation boreholes. Additionally, chemical analyses of ground water, the flow of spring and channel streams are also reflected in measurements [14] ( Table 1). Statistical kriging algorithm is basically used for 3D gridding in order to generate a 3D voxel grid model from 3D data. The data to be gridded must have X, Y, Z and data fields, where X and Y are the real-world Universal Transverse Mercator (UTM) coordinates in meters, Z is the depth in meters and data field is the resistivity in Ω.m. In order to grid the data in 3D, it is critical to state the X, Y and Z cell sizes. These parameters govern the smoothness of resulting 3D grid. The averaging effect increases, by using a larger cell size whereas, a smaller cell size will display more detail, but potentially extra noise in highly anisotropic data. Nodes in the model that have no data within the blank distance of the node will be blank and this blank distance is ruled by the "Cell Size" (Table 2; Figure 6). For example, if blank distance is 4, so grid nodes that are not within 4 cells of observed data will be blank. For the calculation of each node, perceived data will be searched out to evaluate the Maximum radius between the Minimum and Maximum Points. The Maximum radius is designated as the number of cells. If more than the Maximum number of points is found for a particular search radius, only the closest points up to the maximum number of points will be used (GeoSoft Oasis Montaj Tutorial).
These parameters comprise Strike, Dip, Plunge, Strike weight and Dip weight which governs the anisotropic behavior of 3D gridding process. In many geologic settings, there is a natural bias (anisotropy) of the features of interest in the direction of geologic strike. These parameters allow to specify the dominant geologic strike and define preferential weighting in the strike and dipping directions. Strike is in degrees (0-360) clockwise from the coordinate Y axis which is often North and dip is in degrees (0-180) relative to the right perpendicular of strike. Moreover, plunge is the angle of plunge from horizontal and along the strike-dip plane and can be between -90 to +90 degrees. The strike and dip weights (>=1.0) are the elliptical weighting to apply along strike and down dip respectively. The default value of 1.0 in both produces no directional bias in the data. Typically, if data is biased along strike and dip, it needs to set the weightings to 2.0. This doubles the importance of data along strike and dip. If our data has plunge, or is only strike-affected, it needs to set the strike weight to 2.0 and the dip weight to 1.0 respectively. The second and related use of these parameters is to correct anisotropy in the observed sampling pattern. Kriging work best when data is nominally and evenly distributed (GeoSoft Oasis Montaj Tutorial). Weight of Dip 1 Figure 6: (a) Resistivity soundings points presented in 3D, (b) Resultant grid after the application of 3D gridding parameters and (c) Resistivity maps viewed in 3D domain

Induced Polarization
The Induced Polarization (IP) technique in this research assumes the measurement of both resistivity and the polarizability at the same time. This procedure provided a convenient means of examining the dependence of IP response on variables such as grain size, clay type (mineralogy), cation exchange capacity, and water content [16]. Where resistivity is not capable of differentiating between the clay and groundwater formation, the IP can assist in indicating the presence of clay (15). The IP measurements were taken with resistivity for deeper investigations at coarse grid of 20 km x 20 km. The interpret USA software IX1D was used for processing of both VES and IP data.

Interpretation of Resistivity Data
The interpretation of resistivity data to measure the resistivity values which taken from the geophysical survey at the research area. The resistivity values overlap for groundwater because of variations in interaction of contributing factors like porosity, fluid content and salt concentration. The entire range of resistivities in study area is divided into four zones. The VES curves uncovered that the resistivity values of the region are uniform related to similar lithology (Figure 7). The VES results with RMS error are given in table 3.

Interpretation of the pump tests
Depending upon the subsurface geology, the results of the pump test may be simple or complex. All the tests were examined by using Jacob's straight line method ( Figure 8): 1) Test in W-8 well was performed with constant Q of 1412 m3/day which provided the average value of Transmissivity 235 m2/day (Table 1; Figure 8a); 2) In W-9 well, test was operated with constant Q of 250 m3/day which yielded T value to be 20 m2/day (Table 1, Figure 8b); 3) In W-10, test was employed with constant Q 1221 m3/day which generated T value to be 179 m2/day (Table 1; Figure 8c); 4) Finally, W-12 the test was conducted at a constant discharge (Q) of 648.7 m3/day which indicated the average value of Transmissivity 11 m2/day (Table 1; Figure 8d). Furthermore, water conductivity of the aquifer has been determined quite low by the analysis of the pump tests.

Interprertation of Bedrock Analysis
Modeled Resistivity data in 3D ( Figure 6) suggest presence of bedrock at variable depths in all over the area. The map depicts an increase in the resistivity value with depth inferring that the bedrock is occupying more space in the subsurface and decreasing thickness of the alluvium. The resistivity map (i.e. at 300 m and at 500 m) proposes that the bedrock exist at variable depths ( Figure 4). The high resistivity in the eastern and southern portion of the study area implies that the bedrock is composed of igneous rocks (probably granites). Bedrock occurs between 30 m to 60 m in southern and eastern sides whereas deeper in the center of the study area, i.e. more than 200 m depth. The exact depth to the bedrock has not been confirmed in the center of the study area.
The comprehensive interpretation and modeling of resistivity data from both methodologies suggest that the thickness of alluvial cover is more towards center of the project area, and towards the southern and eastern boundaries it becomes comparatively shallow. It contains a huge amount of gravel, boulder and clay layers. The alluvium is relatively coarser along the mountain ranges as indicated by resistivity values (Res > 50 Ohm.m); whereas in the center of project area, the decrease in resistivity values designate the dominance of finer deposits like clay and silt. Moreover, thickness of the alluvium is variable; it is smaller near the mountain ranges and measured more towards the center of the survey area. Schist and granites are underlying mostly at number of places. Although ultimate weathering product of the granite is clay yet the resistant minerals like quartz can be transported to far distances. The coarsegrained granite weathered into pebbles and further movement reduced the sizes which ultimately dumped into the basin to make alluvium. Relatively coarser materials exist near the mountain ranges and the finer materials are downstream. This unit is almost 100 meters thick. However, the alluvium is thick towards the center of the study area and increases above 200 m. In the present research, these boulders and gravel deposits are pondered as worthwhile aquifers. Finer lodges a central part of the plain in the subsurface and is characterized by clay with varying thickness from place to place. During Himalayan orogeny, the uplifting of bed rock as depicted in Figure 8a signify the resistivity maps in a 3D view and consequently, various faults have been induced in the bedrock with the depressions creating a space for sediment fill eroded from the neighboring mountain ranges. Examples of such perceived depression are manifested by blue color. A predominance of clay in various parts of the study area (dominated by blue color) suggests the presence of large standing water bodies at the time of sediments deposition. Borehole confirms the occurrence of gravel patches intermixed with boulders. This is owing to the stream channels which reworked the dumped material through the past and can be regarded as erosional and depositional stage. The gravel patches found in many boreholes imply the presence of stream channels attributing to the flash flooding phases. This clarifies that sediments must have been significantly eroded and replaced by different set of sediments. The alluvial fill is mostly composed of gravel, boulders and having a considerable thickness.
On the other hand, true resistivity is requisite for correlation of resistivity and subsurface geomorphology. This resistivity is attained by modeling of the resistivity soundings. The water table occurs at depths between 13 m to 50 m (obtained from borehole information). By using IPI2Win software, the resistivity data has been modeled and developed into a sequence of homogeneous horizontal layers having individual resistivities. A model resistivity curve designed on this basis will be equivalent to the field curve measured at the surface provided that the assumption of horizontal layers is not violated. The layer model accomplished by this technique due to its large number of layers is too complicated to be used for interpretation. Hence, these models are simplified by exchanging a large number of layers with one single layer so that the equivalent curve remains same. The ultimate curve essentially makes use of existing geology such as maps, lithologic and geophysical data of boreholes, water table depth and subsurface electrical conductivities. Resistivity data has been correlated with the borehole data and diverse lithologies have been allocated different resistivity values. These are presented in Table 4 in following sequence:  Table   The table indicates that the dry sediments above the water table represent comparatively higher value of resistivity commonly in excess of 200ohm.m. However, the increase of fine concentration with the sediments may decrease the resistivity values. The coarse sediments encompassing the gravel and boulders inter-layered with clay reveal resistivity ranging between 50-300 ohm.m. The increase of finer sediments like clay and silt within these sediments decrease the resistivity values whereas the increase of the coarse contents within the sediments will increase the resistivities. The bedrock which contains of metamorphic rocks mainly schist and igneous rocks generally granite is characterized by the resistivity value greater than 300 ohm.m and may reach in excess of 1000 ohm.

Conclusions
Based on resistivity and pump test data (discussed above), it can be concluded that 1. Generally, the thickness of alluvium is more than 200 meters only in a small central area of the Pakhli Plain. Whereas the thickness of the alluvial fill ranges from zero to 120 meters in the southeastern, northwestern and eastern part of the plain. 2. The coarse material like gravels and boulders are represented by low resistivities while clay is characterized by high resistivities. 3. Interpretation of resistivity data (VES) with induce polarization acquired from IX1D revealed that the area is good for groundwater potential as resistivity curves gives 28 to 108 Ohm m values which are indications of fresh water in study area. 4. Bedrock occupies more space with depth and is signified by higher resistivities. 5. The depth of the bedrock increases rapidly from mountains to center of the study area. 6. The pumping test results exhibit that the transmissivities of the alluvial sediments are poor.

Acknowledgement
The present work is a part of M.Phil thesis which was completed by Abid Javed in University of the Punjab, Lahore, Pakistan. The Authors are thankful to Dr. Muhammad Riaz and Dr. Tehseen Zafar for actively participated to modify the manuscript.