Open access peer-reviewed chapter

Appraisal of Groundwater Flow Simulation in the Sub- Himalayan Watershed of Pakistan

Written By

Zulfiqar Ahmad, Arshad Ashraf and Mohsin Hafeez

Submitted: 05 March 2016 Reviewed: 24 March 2016 Published: 27 July 2016

DOI: 10.5772/63324

From the Edited Volume

Groundwater - Contaminant and Resource Management

Edited by Muhammad Salik Javaid

Chapter metrics overview

2,181 Chapter Downloads

View Full Metrics


Numerical modeling of an aquifer is increasingly used as a power tool for monitoring and management of groundwater. This paper focuses on conceptualizing hydrogeological condition and establishing numerical simulation model using Visual MODFLOW to simulate the continuous depletion of groundwater in the southwestern part of the Soan watershed in Pakistan. An integrated groundwater modeling and management approach was adopted to provide suitable alternatives for water management in different hydro-environments. Geospatial techniques were employed for spatial database development, integration with a remote sensing (RS), and numerical groundwater flow modeling capabilities to simulate groundwater flow behavior. The calibration results indicated a reasonable agreement between the calculated and observed heads. The calibrated heads were used as initial conditions in the transient-state modeling. The modeling approach facilitated in identifying potential groundwater regime besides providing artificial recharge options for sustainable groundwater development.


  • groundwater flow modeling
  • aquifer recharge
  • geoinformatics
  • Himalayan watershed

1. Introduction

The growing trend in global warming has resulted in the occurrence of extreme climate events (drought/floods) that impact the hydrological system of the Himalayan region in South Asia. The response of groundwater behavior to environmental changes has been studied worldwide by numerous researchers, for example, see references [14]. The changing pattern of precipitation affects the recharge of the groundwater, which controls the behavior of water table in the system [5,6]. The Indus groundwater system is also facing the temporal as well as spatial impacts of changing environment. A detailed appraisal of the subsurface aquifer system is thus vital in the context of prevailing conditions of drought/flood and the growing demand of water for sustainable development in the region. In recent years, groundwater numerical simulation models have been widely applied to groundwater dynamic simulation and as a management tool. The development and application of a groundwater model is a common practice for the management of groundwater resources [7]. Many scientific efforts have been made to develop more comprehensive and computationally efficient models involving complex hydrogeologic processes [8,9]. In many numerical models of a groundwater aquifer, for example, Feflow and MODFLOW, the continuous domain of groundwater system is replaced by a discretized grid network and the governing groundwater flow equation is solved at the network nodes. The inverse modeling (also called parameter optimization) is capable of assigning approximate values to the hydrological parameters by employing approximate methods to solve the partial differential equation (PDE), which describes the flow in a porous medium.

Pakistan is in the grip of a water crisis. This involves a hydroelectric power shortfall, per capita water availability less than 1100 m3, falling levels of groundwater, and limited water supplies in metropolitan areas, including the twin cities of Rawalpindi and Islamabad. Because pumping exceeds recharge, groundwater reserves are becoming significantly depleted. Groundwater overdraft has caused the groundwater table to decline remarkably and resulted in a series of ecological issues such as deterioration of water quality, soil aridity, deterioration of vegetation, and land desertification [10,11].

In the present study, a three-dimensional numerical groundwater flow modeling using finite-difference Visual MODFLOW coupled with decision support tools of geoinformatics was applied to analyze the spatial and temporal behavior of groundwater in the valley plain of the sub-Himalayan watershed in Pakistan. The multilayered aquifer system of the watershed is recharged mainly by the surrounding dissected land and rocky mountainous terrain. The behavior of the groundwater system was predicted in response to probable hydrogeological stresses.


2. Case study

The growing need of water supplies for both domestic and irrigation purposes has led to the exploitation of surface and groundwater resources in the Soan watershed of Indus basin. In the last drought of 1999–2002, the groundwater abstractions have been increased due to rapid increase in the number of private tube wells for irrigation purpose. Excessive pumpage of groundwater from these wells causes a gradual decline in water table. The decline in groundwater levels in the Rawalpindi area has imposed not only a threat to the livelihood of the residing communities but also a significant influence on the sustainable development of the watershed as a whole. Therefore, understanding and investigation of the groundwater dynamics in the target area have great significance in facilitating the socioeconomic and ecological development, as well as the sustainable development of the water resources. The watershed contains metropolitan areas of major cities of Islamabad and Rawalpindi, the combined population of which is about 2.8 million, with a density of 880 persons/km2. Surface water supplies are maintained from Rawal and Simly reservoirs, which provide 21 million gallons per day (MGD) to Rawalpindi city and 17 MGD to Islamabad city. Groundwater supplies are 24 MGD from about 200 public tube wells in Islamabad, and 27 MGD from about 300 public tube wells in Rawalpindi. The average needs of Islamabad and Rawalpindi are 65 and 175 MGD, respectively, which exceed the capacity of available water resources.

2.1. Description of the study area

The Soan watershed covers an area of about 1684 km2 with its longitude 72° 54′–73° 34′ E and latitude 33° 26′–33° 56′ N in the Pothwar plateau of Pakistan (Figure 1). In the north and northeast, it is bounded by the Margalla and Murree hills, which are covered with permanent mixed scrub and coniferous forest. The altitude in the watershed increases gradually from the southwest (<500 m) toward the northeast, near the Murree hills (>2000 m). The area is characterized by gentle to steep slopes. The dominant formations (e.g., the Murree and Kamlial belonging to the Rawalpindi Group of Miocene age) are composed of sandstone, shale, and lenses of conglomerates. The Lei Nullah conglomerates of Quaternary age consist of poorly sorted pebbles and boulders of mostly Eocene limestone strata [12]. The most important aquifers are composed of gravels and boulders in the unconsolidated sediments of Pleistocene and Recent age. Alluvium (the channel-fill deposits) consists of dominantly silt and clay with subordinate amount of gravel and sand. The soils are mainly gravelly, medium to fine textured over calcareous material in the north and northeast and medium to coarse textured over sandstone in the south (Figure 2). The main perennial stream in the watershed is the Soan, whose primary tributaries are the Ling, Gumrah Kas, Korang, and Lei Nullah. The Korang and Soan rivers are dammed at Rawal and Simly reservoirs, respectively. Due to high urban sprawl, a major area has been converted into various degrees of built-up category in the model domain (Figure 3).

Figure 1.

Location of Islamabad watershed and study area.

Figure 2.

Dominant soil classes in the Soan watershed.

Figure 3.

Major land cover/land use in the study area.

Figure 4a.

Trend in annual precipitation at Islamabad Airport (1960–2013).

Figure 4b.

The mean monthly precipitation at Islamabad airport.

Lei Nullah extends approximately 30 km from the Margalla hills in the northwest to the Soan River at the southeastern edge. The Lei catchment area is composed of 239.8 km2 (169.0 km2 in Islamabad and 70.7 km2 in Rawalpindi) [13,14]. The total rainfall during the monsoon rainy season is about 60% of the annual average rainfall, that is, about 1169 mm. The annual and mean monthly rainfall trends at Islamabad for 1960–2013 periods are shown in Figure 4a and b, respectively. Floods in the Lei Nullah basin occur during the monsoon season, which usually starts near the end of June with peaks in August and finishes by September. In June, the daily maximum temperature reaches 40 °C, while the daily minimum temperature falls near 0 °C in December and January [15].

2.2. Hydrogeologic framework

Silt and clay dominate the subsurface lithology where gravel deposits are present in discontinuous layers with silty clay. The gravel beds are generally 1–20 m thick and are composed of limestone and sandstone pebbles mixed with sand (Figure 5). The thickness of the gravel beds decreases in the south and west. In contrast to the gravels, the bedrock is usually considered to act as an aquitard rather than an aquifer since well yields are much lower than encountered in the gravels. The maximum thickness of the alluvium is more than 200 m, as encountered in the test holes RWP-6 and 8 in Dhok Khabba and Dhok Ratta, respectively. Individual beds range in thickness to 30 m or more. The thickness of alluvium probably exceeds 300 m, as indicated by a deep resistivity survey. They are grayish brown to reddish in color and appear to have been deposited from reworked loess, Siwalik, and Murree rocks. It can become difficult to differentiate rotary cuttings of bedrock formations from the alluvial clays and silts. Sand formations are comparatively rare. Generally, the sand is disseminated within the gravel lenses or constitutes a minor fraction of the clays and silts. The aquifers are composed of seven permeable geological horizons with a total average depth of 137 m. About 300 water wells are being pumped for about 18–22 h per day year round to meet the growing demand of the inhabitants of the city.

Figure 5.

Schematic of lithological section in the study area.


3. Materials and methods

The remote sensing (RS) image data of Landsat-8 were used to analyze the land-cover/land-use status and surface hydrological conditions that were found to be helpful in conceptualizing the recharge/discharge sources of aquifer. Discharge data for the Soan and Korang rivers were obtained from the Water and Power Development Authority (WAPDA) and Surface Water Hydrology Project (SWHP). Both rivers recharge the aquifer system during rainy months and are sustained by base flow at other times. Digital elevation model (DEM) of SRTM 90 m was used to generate the watershed boundary and elevation classes of the watershed area.

The study focuses on conceptualizing hydrogeological condition and establishing numerical simulation model using Visual MODFLOW to simulate the groundwater flow and predict the future response of the aquifer to the growing pumping in Rawal Town (a major part of the Rawalpindi city). The major steps of modeling process include as follows:

  • Collection and review of available hydrogeological data.

  • Collection and review of details of existing tube wells, including their location, depth, and extraction rates.

  • Quantification of discharge rate, with assessment of the rates and periods of extraction.

  • Prediction of the effect of increased extraction or changes in usage patterns and location of the sustainable yield and quality of supply.

  • Preparation of water balance including the assessment of recharge potential in the existing aquifer.

  • Outline recommendation for future exploitation of groundwater aquifer and recharge to groundwater in the study area.

3.1. Modeling groundwater flow

Darcy’s law describes the flow of groundwater and is applied to evaluate aquifer and aquifer material hydraulic characteristics. Darcy’s experiments show that the flow of water through a column of saturated sand is proportional to the difference in the hydraulic head at the ends of the column. Darcy’s law is still used as the basic principle that describes the flow of groundwater (Figure 6) and is expressed as Eq. 1 as follows:



  1. Q = quantity of water discharged, in cfs.

  2. K = hydraulic conductivity (constant factor).

  3. i = hydraulic gradient, in feet.

  4. A = cross-sectional area, in square feet.

The hydraulic gradient (i) determines the direction of groundwater flow following Eq. 2



  1. h = hydraulic head

  2. L = horizontal distance from h1 to h2

Figure 6.

Hydraulic gradient determining the direction of groundwater flow under Darcy’s law.

Hydrological models can be classified based on the application of several criteria, for example, according to the degree of conceptualization of the represented processes as physically based (white-box), conceptual (gray-box), and black-box models. Physically based models use PDEs in space and time to accurately represent in a deterministic way all the processes occurring in the physical system. Black-box models often include stochastic components, for example, relating outputs to inputs through a set of empirical functions, such as simple mathematical expressions, time series equations, autoregressive moving average (ARMA) models, and artificial neural networks (ANNs) disregarding any physical insight. Somewhere in between, conceptual models often represent parts of the system as series of tanks that exchange water with one another. Through a quite simplified description, they represent all the relevant parts of hydrological processes.

3.2. Numerical simulation of groundwater flow

The applications of MODFLOW, a modular three-dimensional finite-difference groundwater model of the US Geological Survey, to the description and prediction of the behavior of groundwater systems have increased significantly over the last few years. It has become the worldwide standard groundwater flow model that is used to simulate a wide variety of systems like water supply, containment remediation, and mine dewatering [16]. Several packages are integrated with MODFLOW that deals with a specific feature of the hydrologic system to be simulated, such as wells, recharge, or river. Models or programs can be stand-alone codes or can be integrated with MODFLOW. A stand-alone model or program communicates with MODFLOW through data files, for example, the advective transport model PMPATH [17], the parameter estimation programs PEST [18], and UCODE [19] use this approach.

PMWIN comes with a professional graphical user interface, the supported models and programs, and several other useful modeling tools, for example, Time-Variant Specified-Head (CHD1) [20], Direct Solution (DE45) [21], Density (DEN1) [22], Horizontal-Flow Barrier (HFB1) [23], Interbed-Storage (IBS1) [20], Reservoir (RES1) [24], and Streamflow-Routing (STR1) package [25]. Simulation results include hydraulic heads, drawdowns, cell-by-cell flow terms, compaction, subsidence, Darcy velocities, concentrations, and mass terms. The graphical user interface allows one to create and simulate models and displays temporal development curves of simulation results including hydraulic heads, drawdowns, subsidence, compaction, and concentrations. In many cases, the development of effective and efficient automatic calibration procedures, based on numerical optimization methods, has replaced the calibration of hydrological models conducted manually by “trial-and-error” parameter adjustment. The purpose of PEST and UCODE is to assist in data interpretation and in model calibration. If there are field or laboratory measurements, PEST and UCODE can adjust model parameters and/or excitation data in order that the discrepancies between the pertinent model-generated numbers and the corresponding measurements are reduced to a minimum. Both codes do this by taking control of the model (MODFLOW) and running it as many times as is necessary in order to determine this optimal set of parameters and/or excitations.

The software package Visual MODFLOW [26] has been used for three-dimensional numerical simulations of groundwater flow and drawdown within the study area. A finite-difference rectangular grid was constructed by establishing a network of nodal points. The model computes drawdown, direction of flow with vector lines, and hydraulic heads on each nodal point [27]. The flow model requires a boundary array for each cell. A positive value in that array defines an active cell (the hydraulic head is computed), a negative value defines a fixed-head cell (the hydraulic head is kept constant throughout the simulation), and the zero value defines an inactive cell (no flow takes place within the cell). A fixed-head boundary exists whenever an aquifer is in direct hydraulic contact with a river, a lake, or a reservoir in which the water level is known. Such a boundary provides an inexhaustible supply of water, which in some situations may be unrealistic [4, 28,29]. Visual MODFLOW requires initial hydraulic heads at the beginning of a flow simulation. For steady-state simulations, the initial heads are starting guessed values for the iterative equation solvers. The heads at the fixed-head cells must be the actual values while all other initial heads can be set arbitrarily. For transient-flow simulations, the initial heads must be the actual values [3033].

3.3. Model conceptualization and setup

The modeling area is bounded by the Korang River in the southeast and the Lei Nullah in the southwest. The hydrogeological system of the area was modeled as multilayered using Finite-difference Visual MODFLOW. The major steps involved during the modeling process were as follows: (a) expression of field parameters such as piezometric head, hydraulic conductivity; (b) formulation of the groundwater flow equation in an integral form; (c) integration of the integral form of the groundwater equation; (d) assembly of the algebraic matrix equations that result from the integration step into global system of linear equation; (e) time integration; and (f) solution of the global system of linear equations. The data from test holes and wells drilled by WAPDA [34,35] had been used in this study. Table 1 shows specific yield and retention percentages for various lithologies. Hydraulic characteristics obtained from the analysis of pumping test data are provided in Table 2. We visited 278 tube wells in Rawal and Pothwar Towns for water-level measurements, as summarized in Table 3 and shown in Figure 7. The mean water table depth was 28 m.

Material Porosity Specific yield Specific retention
Soil 55 40 15
Clay 50 2 48
Sand 25 22 3
Gravel 20 19 1
Limestone 20 18 2
Sandstone 11 6 5

Table 1.

Specific yield and retention percentages (values in percent by volume).

Test hole Easting Northing Transmissivity
coefficient m/s
Specific capacity m3/day/m
TH-1 3214782 1053687 1.5 × 10−2 5.6 × 10−4 7 × 10−2 607
TH-6 3211368 1048417 3.6 × 10−3 1.5 × 10−4 2 × 10−3 116
TH-8 3212866 1047817 1.4 × 10−3 6.9 × 10−5 2 × 10−4 116
TH-9 3212299 1049052 1.3 × 10−2 4.6 × 10−4 5 × 10−2 840

Table 2.

Hydraulic characteristics of the aquifer in the study area.

Category No. of tube wells
East zone West zone
Water supply tube wells in Rawal Town 120 100
Water supply tube wells in Pothwar Town 24
New tube wells installed 26
Private tube wells visited 8
Total 278

Table 3.

Brief summary of the tube wells used in the modeling study.

Figure 7.

Location of the water supply public tube wells surveyed in the study area.

There are two major sources of recharge to groundwater: precipitation and infiltration from surface streams. Direct infiltration to the water table from precipitation is likely to occur especially in July and August, when rainfall is highest, although evaporation and soil moisture deficit in these months are also high. Recharge is also possible during the winter rains in February and March. We assume that 30% of rainfall is lost through runoff into the Soan and Korang rivers, 40% is lost through evapotranspiration from the surface and soil zone, and 15% is lost through the capillary and deep roots. Therefore, about 15% of rainfall, that is, 5.58 × 10−9 m/s [3638], is available to recharge the shallow groundwater. In the Margalla hills, 7–15% of the precipitation tends to recharge the aquifer system through deep fractures of variable sizes in the limestone. Recharge within the urban areas of Islamabad and Rawalpindi is limited by impervious cover. The setting of the flow model is given below:

  • Grids 95 columns × 89 rows. The size of grid is 91.4 m (300 ft) in x and y directions

  • The total modeling area is 70.7 km2

  • Number of layers = 7

  • Layer 1 (25 m) type = unconfined/confined (transmissivity [T] constant)

    • Layer 2 (12 m) type = aquiclude

    • Layer 3 (20 m) type = confined (T constant)

    • Layer 4 (12 m) type = aquiclude

    • Layer 5 (25 m) type = confined (T constant)

    • Layer 6 (4 m) type = aquiclude

    • Layer 7 (22 m) type = confined (T constant)

T-values are used from the results of pumping tests at 35 tube wells. Overall, T is set at 2.8 × 10−3 m2/s for layer 1, 4.6 × 10−4 m2/s for layer 3, 6.85 × 10−4 m2/s for layer 5, and 2.78 × 10−4 m2/s for layer 7. Boundary conditions consider all the cells active with the Korang River as a constant head boundary and the Lei Nullah as a recharge boundary. Horizontal hydraulic conductivity varies depending upon the type of subsurface material in layer 1 and is computed from T. The initial hydraulic heads are set at 28 m from initial guess and observed heads in the field. The top of the layer 1 is set at 426.7 m above the mean sea level (masl), and the bottom of the layer 1 is set at 47 m from the ground surface including 12 m layer of aquiclude. The top of the layer 2 is set at 47 m from the ground surface. The bottom of the layer is set at 106.7 m. Therefore, the thickness of layer 2 (saturated) is considered to be 57.9 m including 12-m layer of aquiclude. The top of the layer 3 is set at 106.7 m from the ground surface. The bottom of the layer is set at 137 m. Therefore, the thickness of layer 3 (saturated) is considered to be 30.4 m including 4-m layer of aquiclude. The specific yield (required for flow-path calculations) is set at 0.06 for layer 1, 0.009 for layer 3, 0.005 for layer 5, and 0.007 for layer 7. The Korang River is considered as the constant head boundary. Lei Nullah was considered to be influent (losing) or effluent (gaining) at places that would automatically adjust the situation during the simulations of different stress period in the steady-state and non-steady state conditions. The River Package was used to assign the following values of hydraulic properties of the rivers to the model cells:

  • Hydraulic conductance of the river bed = 1.5 × 10−4 m2/s

  • Head in the river = 27 m from the top of the layer = 0.0

  • Elevation of the river bed bottom = 27.5 m from the top of the layer = 0.0

Observed hydraulic heads were obtained for 1998 and 2003 from the previous literature review [14], while for 2007 data were collected physically in the field. When data were collected, only one well was switched off while the others were pumping, as it was not actually possible to shut off all the wells simultaneously. Consequently, spikes were apparent for hydraulic heads measured in 2007. Three-dimensional projections of the hydraulic heads are shown in Figures 810.


4. Results and discussion

4.1. Simulation of the groundwater flow

First, the model was run for steady-state condition with one time step of 20 years duration. The model was calibrated by comparing observed and computed heads using the UCODE automatic calibration program [19,39]. The hydraulic conductivity and recharge values were estimated during the process of steady-state calibration. Initially, the model was run to achieve the steady-state condition without pumping the wells and results obtained were calibrated with the available hydraulic head of 1998, as the available heads in 2003 and 2007 were very fluctuating due to the pumping effects. Simulated steady-state hydraulic heads for different aquifer layers are shown in Figure 11. Vector plots indicate groundwater flows toward Korang River in the southeast. Mass balance analysis showed a good balance between the inflow and the outflow components during the steady-state condition (Figure 12).

Figure 8.

Three-dimensional projection of observed hydraulic heads in non-steady-state condition (1998).

For transient-flow simulations, the initial heads were taken from the steady-state model. Withdrawals were set to a realistic range of pumping rates (in the range of 0.0094–0015 m3/s). The non-steady-state model was run for 1, 3, and 5 years. Transient-state simulations carried out up to year 2012 showed that the flow field had changed in the underlying aquifer layers because of pumping, as indicated by the arrowheads in Figure 13. Cone of depressions of layers 1, 3, and 7 for the southwestern area during 2008 are shown in Figure 14ac. In the southwest (around Rawal Town), the flow is diverted toward a composite depression trough, but the flow in the southeastern area (around Pothwar Town) still moves into the Korang River. In the extreme south, the flow tends to rise where several wells are not pumping. A maximum of 20 m drawdown was observed at the end of a 3-year simulation in 2010. As a result, groundwater in the upper horizon (47 m) has already been depleted and pumping continues in the remaining permeable layer 3 through layer 7.

Figure 9.

Three-dimensional projection of observed hydraulic heads with reference to water supply tube wells (2003).

Figure 10.

Three-dimensional projection of observed hydraulic heads with reference to water supply tube wells (2007).

Figure 11.

Steady-state hydraulic head in layers 1, 3, and 7 (after 20 years).

Figure 12.

Mass balance during steady-state calibration.

Figure 13.

Hydraulic head in layers 3, 5, and 7 during 2010 (non-steady state).

The physical properties of the aquifer were found favorable for the development of groundwater with moderate values of transmissivity and specific yield. There is a continuous drop in hydraulic heads in the vicinity of pumped tube wells in Rawal Town and adjacent areas. Simulations indicated a maximum drawdown of 20 m. Existing tube wells could be safely pumped for a period of 5 years at the existing rate of discharge. However, drawdown could reach as deep as the fifth layer of the aquifer. Velocity vectors indicate that the groundwater in Rawalpindi city moves southeast toward the Korang River. Groundwater flow pattern in the layers 3, 5, and 7 is almost identical, but the velocity of flow is different for different layers.

Figure 14a.

Cone of depression in layer 3 (1-year simulation in 2008, non-steady state).

Figure 14b.

Cone of depression in layer 5 (1-year simulation in 2008, non-steady state).

Figure 14c.

Cone of depression in layer 7 (1-year simulation in 2008, non-steady state).

4.2. Management options for groundwater development

Under the present water scarcity conditions, there is a need to initiate integrated water resources management programs, with site-specific interventions, especially to harness the available rainwater. This would not only contribute to groundwater recharge in the basin but also supplement the water supplies to meet future water demand for various uses. Surface water drainage system/network (Figure 1) carries rainwater and is a potential source of recharge to groundwater. The Soan and Korang rivers recharge the aquifer system during rainy months and are sustained by base flow at other times. Therefore, more water wells could be constructed along the Korang River where the underlying aquifer layers have good potential yield. Pumping of water wells in the western part of the study area needs to be monitored strictly in order to enable the groundwater levels to recover. Necessary bye-laws are established to restrict the overexploitation of groundwater and installation of unplanned tube wells. Similarly, bye-laws are established for restricted clearance of vegetation cover for unplanned urban development, as vegetative cover on the free catchment has proven to be an aid to groundwater recharge. The recharging system should be so designed that the amount of recharge during a year is at least equal to the amount of groundwater extracted during the year. Some of the recharging techniques may include as follows:

4.2.1. Dams and pounds

Water-charging pounds in open spaces, dams, and check dams need to be constructed in nullahs and distributaries. Such ponds would serve as water storage for emergency as well as for recharging of the groundwater (Figure 15). Main parks and green belts astride main roads of the twin cities and the dissected valleys around the built-up lands are ideal places to harvest rainwater by making water ponds recharging wells, etc.

Figure 15.

Water ponds provide potential source not only for storing the surface runoff but also for recharging the groundwater.

4.2.2. Infiltration using wells and boreholes

Water can be infiltrated by injection using wells or boreholes in the areas where low-permeability strata overlies target aquifers as in the parts of watershed containing medium- to fine-textured soils over calcareous material (Figure 2). This technique is suitable for deep-seated aquifers that form a source of groundwater for urban areas lying in the valley plains of the watershed. The water-injecting wells have advantage that recharge water can bypass thick impervious layers to be introduced to the most permeable portions of the aquifer [40].

4.2.3. Water spreading

Water is diverted from surface water runoff into infiltration basins, ditches, or low-lying areas where the aquifer to be recharged is at or near to the ground surface as in the areas with flat to gentle topography. The water-spreading method could be practiced in medium- to coarse-textured soils in the gently sloping southeastern and some central parts of the watershed (shown in Figure 2). The exposed rocks of limestone and sandstone in the vicinity of the urban areas also provide potential recharge environment for subsurface water (Figure 16).

Figure 16.

Joints and openings in the rocks facilitate recharge to the subsurface water.



The authors gratefully acknowledge the Australian Department of Education, Employment, and Workplace Relations (DEEWR) for providing post-selection support services and funding the research under the Executive Endeavour Award 2010. The Water and Power Sanitation Authority (WASA), Asian Development Bank (ADB), Osmani & Company (Pvt) Ltd, Department of Earth Sciences, Quaid-i-Azam University (QAU), and College of Earth and Environmental Sciences, Punjab University, are acknowledged for providing access to reliable data. Dr. Alan Fryar and Terencen Hamilton from the University of Kentucky, USA, are gratefully acknowledged to review and edit the manuscript.


  1. 1. Allen DM. Historical trends and future projections of groundwater levels and recharge in coastal British Columbia, Canada, SWIM21 – 21st Salt Water Intrusion Meeting, 21 26 June 2010, Azores, Portugal. 2010: 267–270.
  2. 2. Crosbie, Russell S, McCallum, James L, Walker, Glen R and Chiew, Francis HS. Episodic recharge and climate change in the Murray-Darling Basin, Australia, Hydrogeology Journal. 2012; 20(2): 245–261.
  3. 3. Okkonen, Jarkko, Jyrkama, Mikko and Kløve, Bjørn. A conceptual approach for assessing the impact of climate change on groundwater and related surface waters in cold regions (Finland), Hydrogeology Journal. 2010; 18: 429–439.
  4. 4. Ashraf A and Ahmad Z. Regional groundwater flow modeling of upper Chaj Doab, Indus Basin, Geophysical Journal International. 2008; 173(1): 17–24.
  5. 5. McCallum JL, Crosbie RS, Walker GR and Dawes WR. Impacts of climate change on groundwater in Australia: a sensitivity analysis of recharge, Hydrogeology Journal. 2010; 18: 1625–1638.
  6. 6. Kumar CP. Impact of Climate Change on Groundwater Resources, Munich, GRIN Verlag, .2014. (Accessed on May 2, 2016)
  7. 7. Mohammadi Z. Variogram-based approach for selection of grid size in groundwater modeling, Applied Water Science. 2013; 3: 597–602.
  8. 8. Noronha A and Lee J. On the use of information theory to quantify parameter uncertainty in groundwater modeling, Entropy. 2013; 15: 2398–2414.
  9. 9. Dams Jef, Salvadore, Elga and Batelaan, Okke. Predicting impact of climate change on groundwater dependent ecosystems, 2nd International Interdisciplinary Conference on Predictions for Hydrology, Ecology and Water Resources Management: Changes and Hazards caused by Direct Human Interventions and Climate Change, 20–23 September 2010, Prague, Czech Republic.
  10. 10. Wei XM, Kang SZ, Su XL, et al. Impact of oasis agricultural development on the transforming relationship between surface water and groundwater in the Shiyang River Basin, Transactions of the Chinese Society of Agricultural Engineering. 2005; 21(5): 38–41 (in Chinese).
  11. 11. Ma L, Wei XM. Analysis on aberrance point of annual runoff serials in the downstream of Shiyang River, Agricultural Research in the Arid Areas. 2006; 24(2): 174–177 (in Chinese).
  12. 12. Allchin B. The palaeolithic of the Pothwar plateau Punjab, Pakistan: a fresh approach, Paléorient. 1981; 7(1): 123–134. doi:10.3406/paleo.1981.4291.
  13. 13. Bender FK and Raza HA (Eds.). Geology of Pakistan Tutte Druckerei, GmbH, Salzweg-passau. 1995.
  14. 14. JICA. Compilation of report of Islamabad–Rawalpindi ground water study by the Japan International Cooperation Agency, Capital Development of Authority, Pakistan. 1988.
  15. 15. Ahmad B, Kaleem MS, Butt MJ and Dahri ZH. Hydrological modelling and flood hazard mapping of Nullah Lai, Proceedings of the Pakistan Academy of Sciences. 2010a; 47(4): 215–226.
  16. 16. Kumar CP. Numerical modelling of ground water flow using MODFLOW, Indian Journal of Science. 2013; 2(4): 86–92.
  17. 17. Chiang WH and Kinzelbach W. 3D-Groundwater modeling with PMWIN, A simulation system for modeling groundwater flow and pollution. Springer. ISBN: 978-3-662-05551-9. 2003.
  18. 18. Doherty J. PEST ver. 2.04, Watermark computing, Corinda, Australia. 1995.
  19. 19. Poeter EP and Hill MC. Documentation of UCODE, a computer code for universal inverse modelling, U.S. Geological Survey, Water-Resources Investigations Report 98-4080. Denver, Colorado. 1998.
  20. 20. Leake SA and Prudic DE. Documentation of a computer program to simulate aquifer system compaction using the modular finite-difference ground-water flow model. Open-File Report 88-482, U.S. Geological Survey. 1991
  21. 21. Harbaugh AW. Direct solution package based on alternating diagonal ordering for the U.S. Geological Survey modular finite-difference ground-water flow model: U.S. Geological Survey Open-File Report 95-288; 1995: 46 pp.
  22. 22. Schaars FW and van Gerven MW. Density package, simulation of density driven flow in MODFLOW. KIWA-report SWS 97.511, ISBN 90-74741-42-8. KIWA Research & Consultancy, Nieuwegein, The Netherlands. 1997.
  23. 23. Hsieh and Freckleton. Documentation of a computer program to simulate horizontal-flow barriers using the U.S. Geological Survey’s modular three-dimensional finite-difference groundwater flow model. U.S. Geological Survey, Open-File Report 92-477. 1993.
  24. 24. Fenske JP, Leake S A and Prudic DE. Documentation of a computer program (RES1) to simulate leakage from reservoirs using the modular finite-difference ground-water flow model (MODFLOW), U.S. Geological Survey, Open-File Report 96-364. 1996.
  25. 25. Prudic DE. Documentation of a computer program to simulate stream-aquifer relations using a modular, finite-difference, ground-water flow model, U.S. Geological Survey, Open-File Report 88-729, Carson City, Nevada. 1988.
  26. 26. Schlumberger Water Services (SWS)., 2010 Waterloo, ON N2L 5J2, (Accessed on May 5, 2016)
  27. 27. Pollock DW. Semi-analytical computation of path lines for finite difference models, Ground Water. 1988; 26(6): 743–750.
  28. 28. Ahmad Z and Sergio Serrano. Numerical groundwater flow modeling of the Rechna Doab Aquifer, Punjab, Pakistan, Proceeding of International Groundwater Modeling Conference (IGWMC), MODFLOW 2001 and other Modeling Odysseys, Colorado School of Mines and Geophysics, Denver (12–14 September), USA. 2001.
  29. 29. Ahmad Z and Khawaja AA. Three-dimension flow and solute transport modeling of the Eolian Aquifer in a uniform flow field of a Chemical Plant, Pakistan Journal of Hydrocarbon Research. 1999; 11: 59–74.
  30. 30. McDonald MG and Harbaugh AW. A modular three-dimensional finite difference groundwater flow model. USGS open-file report 83-875: Modeling Techniques, Book 6. 1996.
  31. 31. Ahmad Z, Ashraf M and Siddiqui RF. Groundwater computer modeling of the National Park area in the vicinity of Rawal Lake, Islamabad, Journal of Drainage and Water Management. 1997; 1: 8–21.
  32. 32. Elci A, Molz FJ and Waldrop WR. Implications of observed and simulated ambient flow in monitoring wells, Ground Water. 2001; 39(6): 853–862.
  33. 33. Ghulam M and Ahmad Z. Management of groundwater resources in Punjab, Pakistan, using a groundwater flow model, Journal of Environmental Hydrology. 2007; 15: 1–14.
  34. 34. WAPDA, Lower Indus Report. Physical Resources, groundwater, supl.6.1.6, Tubewells and Boreholes. Lahore, Pakistan. 1965.
  35. 35. WAPDA, Regional groundwater investigations in Islamabad and Rawalpindi areas, supl.12. Ground water studies. Lahore, Pakistan. 1980.
  36. 36. Ahmad K. Pakistan: Lei Nullah basin flood problem Islamabad-Rawalpindi cities, WMO/GWP Associated programme on flood management. (Accessed on May 5, 2016)
  37. 37. Ahmad Z, Kausar R and Ahmad I. Implications of depletion of groundwater levels in three layered aquifers and its management to optimize the supply demand in the urban settlement near Kahota Industrial Triangle area, Islamabad, Pakistan, Environmental Monitoring and Assessment. 2010b; 166(1–4): 41–55.
  38. 38. Ahmad Z, Kausar R and Ahmad I. Highlighting the implications of selenium dispersion from disposal of Kahota industrial triangle area, Islamabad, Pakistan using three-dimension solute transport model, Environmental Monitoring and Assessment. 2010c; 167(1–4): 565–579.
  39. 39. Poeter EP and Hill MC. UCODE, a computer code for universal inverse modelling, Computers and Geosciences. 1999; 25: 457–462.
  40. 40. Gale I, et al. The effectiveness of artificial recharge of groundwater: a review, British Geological Survey, Commercial Report CR/02/108N. 2002.

Written By

Zulfiqar Ahmad, Arshad Ashraf and Mohsin Hafeez

Submitted: 05 March 2016 Reviewed: 24 March 2016 Published: 27 July 2016