**3. Materials and methods**

**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

Silt and clay dominate the subsurface lithology where gravel deposits are present in discon‐ tinuous 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

December and January [15].

**2.2. Hydrogeologic framework**

8 Groundwater - Contaminant and Resource Management

of the inhabitants of the city.

**Figure 5.** Schematic of lithological section in the study area.

The remote sensing (RS) image data of Landsat-8 were used to analyze the land-cover/landuse 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:


#### **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 ground‐ water (**Figure 6**) and is expressed as Eq. 1 as follows:

$$
\underline{Q} = \text{K}iA\tag{1}
$$

where


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

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

$$\frac{\dot{a} = (h\_1 - h\_2)}{L} \tag{2}$$

where

*h* = hydraulic head

*L* = horizontal distance from *h*1 to *h*2

**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**

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

10 Groundwater - Contaminant and Resource Management

*L* = horizontal distance from *h*1 to *h*2

where

*h* = hydraulic head

of hydrological processes.

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

1 2 *i hh* ( ) *L*

**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

= - (2)

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 standalone codes or can be integrated with MODFLOW. A stand-alone model or program commu‐ nicates 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 fixedhead 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 Finitedifference 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.


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


**Table 2.** Hydraulic characteristics of the aquifer in the study area.


**Table 3.** Brief summary of the tube wells used in the modeling study.

must be the actual values while all other initial heads can be set arbitrarily. For transient-flow

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 Finitedifference 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

simulations, the initial heads must be the actual values [30–33].

**Material Porosity Specific yield Specific retention**

**Permeability coefficient m/s** **Storage coefficient** **Specific capacity m3**

**/day/m**

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).

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

**m2 /s**

**Table 2.** Hydraulic characteristics of the aquifer in the study area.

**3.3. Model conceptualization and setup**

12 Groundwater - Contaminant and Resource Management

mean water table depth was 28 m.

**Test hole Easting Northing Transmissivity**

**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 [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:

	- **◦** 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:


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 8**–**10**.
