Specific yield and retention percentages (values in percent by volume).
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
- Himalayan watershed
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 [1–4]. 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 . 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 . 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).
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 .
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.
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:
Q = quantity of water discharged, in cfs.
K = hydraulic conductivity (constant factor).
i = hydraulic gradient, in feet.
A = cross-sectional area, in square feet.
The hydraulic gradient (i) determines the direction of groundwater flow following Eq. 2
h = hydraulic head
L = horizontal distance from h1 to h2
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 . 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 , the parameter estimation programs PEST , and UCODE  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) , Direct Solution (DE45) , Density (DEN1) , Horizontal-Flow Barrier (HFB1) , Interbed-Storage (IBS1) , Reservoir (RES1) , and Streamflow-Routing (STR1) package . 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  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 . 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 [30–33].
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|
|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|
|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|
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 [36–38], 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 , 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 8–10.
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).
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 14a–c. 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.
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.
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.
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 .
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).
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.