**Integration of Groundwater Flow Modeling and GIS**

Arshad Ashraf1 and Zulfiqar Ahmad2 *1National Agricultural Research Center, Islamabad 2Department of Earth Sciences, Quaid-i-Azam University, Islamabad Pakistan* 

#### **1. Introduction**

238 Water Resources Management and Modeling

Sugawara, M., 1995. Tank Model, in "Computer Models of Watershed Hydrology". Singh

Tritz, S., Guinot, V., Jourde, H., 2011. Modelling the behaviour of a karst system catchment

using non linear hysteretic conceptual model. Journal of Hydrology, Journal of

Singh, V.P., 1988. Hydrologic systems, rainfall-runoff modeling. Prentice Hall, NJ.

V.P. [Ed.]. Water Resources Publications, Colorado, pp. 165–214.

Hydrology 397(3-4): 250-262.

The development of a sufficient understanding on which to base decisions or make predictions often requires consideration of a multitude of data of different types and with varying levels of uncertainty. The data for the development of numerical groundwater flow model includes time-constant parameters and time-variant parameters. The time-constant parameters were mainly extracted from thematic data layers generated from GIS and image processing of remote sensing data. Integrated approaches in GIS play a rapidly increasing role in the field of hydrology and water resources development. It provides suitable alternatives for efficient management of large and complex databases developed in different model environments. Remote sensing (RS) technology is capable of providing base for quantitative analysis of an environmental process with some degree of accuracy. It provides an economic and efficient tool for landcover mapping and has its advantages in planning and management of water resources (Ashraf and Ahmad 2008 & Ahmad et al., 2011). One of the greatest advantages of using remote sensing data for hydrological investigations and monitoring is its ability to generate information in spatial and temporal domain. Though, the presence of groundwater cannot be directly ascertained from RS surveys, however, satellite data (GRACE- Gravity Recovery and Climate Experiment and GOCE-Gravity field and steady-state Ocean Circulation Explorer) provides quick and useful baseline information on the parameters that control the occurrence and movement of groundwater such as geomorphology, direction of groundwater flows, lineaments, soils, landcover/landuse and hydrology etc.

Reliable investigation of waterlogging problem can be extremely useful in chalking out suitable water management strategies by reclaiming existing waterlogged areas. The problems of water logging and salinity mostly exist in the irrigated areas like in Indus plains of South Asia. This twin menace generally results from over irrigation, seepage losses through canal and distributory system, poor water management practices and inadequate provision of drainage system. Analysis of high watertable in waterlogged areas and drainage of irrigated areas have not been paid adequate attention in the planning process, partly due to lack of requisite data and partly due to resource crunch in the country. In order to develop suitable water management strategies and controlling the extent of waterlogging in the area, reliable investigation and clear apprehension of the groundwater flow system bear great significance. It is important not only to facilitate the reconstruction of the ecological environment but also to accommodate the sustainable development of the water resources and economy of this region.

According to Choubey (1996), a rapid and accurate assessment of the extent of waterlogged areas could be made using remotely sensed data. Remote sensing (RS) and geographical information system (GIS) offers convenient solutions to map the extent and severity of waterlogging and salinity, particularly in large areas (IDNP 2002). Arora and Goyal (2003) highlighted the use of GIS in development of conceptual groundwater model. Using the logical conditions and analytical functions of GIS domain, the recharge/discharge zones can be delineated effectively (Ashraf et al., 2005). Such zones often provide clues of causative factors of the waterlogging problem. It should be realized that by just using last century's schemes no longer solves challenges related to today's groundwater situation (Zaisheng et al., 2006). In recent years, groundwater numerical simulation models like Feflow which is based on finite element method have been widely applied to groundwater dynamics simulation due to its appealing features such as visualization, interaction and diversified solving methods (Yang and Radulescu, 2006; Russo and Civita, 2009; Peleg and Gvirtzman, 2010; Dafny et al., 2010).

In the present study, remote sensing and numerical groundwater flow modeling were integrated in GIS environment to identify and analyze waterlogged areas and associated groundwater behavior spatially and temporally. A rapid and accurate assessment of the extent of waterlogged areas could be made using remotely sensed and groundwater data. The GIS is used for spatial database development, integration with a remote sensing (RS) and numerical groundwater flow modeling capabilities. Remote sensing technique in conjunction with (GIS) provided a quick inventory of waterlogged area and its monitoring. Finite element groundwater flow model (Feflow 5.1) was developed to configure groundwater equipotential surface, hydraulic head gradient and estimation of the groundwater budget of the aquifer. Steady-state simulations were carried out to describe three-dimensional flow field (head distribution) over the entire model domain. The calibrated heads were used as initial conditions in the transient-state modeling. For transient and predictive modeling, management strategies for pre-stress period and post stress periods were developed on the basis of present water use and future requirements. During transient simulations, model was rerun for period of about 15 years i.e. from 2006 to 2020 to simulate groundwater flow behavior in the model domain. Scenarios of impact of extreme climate events (drought/flood condition) and variable groundwater pumpage on groundwater levels were studied. The integration of GIS with groundwater modeling and Remote sensing (RS) provided efficient way of analyzing and monitoring resource status and land conditions of waterlogged areas. The approach would help in providing an effective decision support tool for evaluating better management options to organize schemes of future monitoring of groundwater resource and waterlogging risks on local and regional basis.

#### **1.1 GIS Functions and methods**

The effort to perform analysis of management scenarios will be substantially reduced by an easily accessible database, a convenient interface between database and groundwater models, visualization and utilities for model inputs and results (Pillmann and Jaeschke, 1990). GIS is an application-oriented spatial information system with a variety of powerful functions to handle for decision support of problems related to spatial dimension. All of the data in a GIS are georeferenced, that is, linked to a specific location on the surface of the Earth through a system of coordinates. Geographical information attaches a variety of qualities and characteristics to geographical locations. These qualities may be physical parameters such as ground elevation, soil moisture or groundwater levels, or classifications according to the type of vegetation and landcover, ownership of land, zoning, and so on. Such occurrences as accidents, floods, or landslides may also be included. A general term 'attributes' is often used to refer to the qualities or characteristics of places and is considered as one of the two basic elements of geographical information, along with locations. GIS technology integrates common database operations such as query and statistical analysis with the unique visualization and benefits of geographic analysis offered by maps. GIS is the most appropriate tool for spatial data input and attribute data handling. It is a computer-based system that provides the following four sets of capabilities to handle georeferenced data (Aronoff, 1989):

Data input

240 Water Resources Management and Modeling

flow system bear great significance. It is important not only to facilitate the reconstruction of the ecological environment but also to accommodate the sustainable development of the

According to Choubey (1996), a rapid and accurate assessment of the extent of waterlogged areas could be made using remotely sensed data. Remote sensing (RS) and geographical information system (GIS) offers convenient solutions to map the extent and severity of waterlogging and salinity, particularly in large areas (IDNP 2002). Arora and Goyal (2003) highlighted the use of GIS in development of conceptual groundwater model. Using the logical conditions and analytical functions of GIS domain, the recharge/discharge zones can be delineated effectively (Ashraf et al., 2005). Such zones often provide clues of causative factors of the waterlogging problem. It should be realized that by just using last century's schemes no longer solves challenges related to today's groundwater situation (Zaisheng et al., 2006). In recent years, groundwater numerical simulation models like Feflow which is based on finite element method have been widely applied to groundwater dynamics simulation due to its appealing features such as visualization, interaction and diversified solving methods (Yang and Radulescu, 2006; Russo and Civita, 2009; Peleg and Gvirtzman,

In the present study, remote sensing and numerical groundwater flow modeling were integrated in GIS environment to identify and analyze waterlogged areas and associated groundwater behavior spatially and temporally. A rapid and accurate assessment of the extent of waterlogged areas could be made using remotely sensed and groundwater data. The GIS is used for spatial database development, integration with a remote sensing (RS) and numerical groundwater flow modeling capabilities. Remote sensing technique in conjunction with (GIS) provided a quick inventory of waterlogged area and its monitoring. Finite element groundwater flow model (Feflow 5.1) was developed to configure groundwater equipotential surface, hydraulic head gradient and estimation of the groundwater budget of the aquifer. Steady-state simulations were carried out to describe three-dimensional flow field (head distribution) over the entire model domain. The calibrated heads were used as initial conditions in the transient-state modeling. For transient and predictive modeling, management strategies for pre-stress period and post stress periods were developed on the basis of present water use and future requirements. During transient simulations, model was rerun for period of about 15 years i.e. from 2006 to 2020 to simulate groundwater flow behavior in the model domain. Scenarios of impact of extreme climate events (drought/flood condition) and variable groundwater pumpage on groundwater levels were studied. The integration of GIS with groundwater modeling and Remote sensing (RS) provided efficient way of analyzing and monitoring resource status and land conditions of waterlogged areas. The approach would help in providing an effective decision support tool for evaluating better management options to organize schemes of future monitoring of groundwater resource and waterlogging risks on local and

The effort to perform analysis of management scenarios will be substantially reduced by an easily accessible database, a convenient interface between database and groundwater models, visualization and utilities for model inputs and results (Pillmann and Jaeschke,

water resources and economy of this region.

2010; Dafny et al., 2010).

regional basis.

**1.1 GIS Functions and methods**


All the geographically related information that is available can be input and prepared in GIS such that user can display the specific information of interest or combine data contained within the system to generate further information which might answer or help resolve a specific problem. The hardware and software functions of a GIS include compilation, storage, retrieval, updating and changing, manipulation and integration, analysis and presentation. All of these actions and operations are applied by a GIS to the geographical data that form its database. The data representation phase deals with putting all information together into a format that communicates the result of data analysis in the best possible way. Geographic Information System can be represented with geometric information such as location, shape and distribution, and attribute information such as characteristics and nature. Any spatial features of the earth's surface are represented in GIS by the *Polygons (*features which occupy a certain area, e.g. lakes, landuse, geological units etc.); *Lines (*linear features, e.g. drainage, contours, boundaries etc); *Points (*points define the discrete locations of geographic features like cities/towns, boreholes, wells, spot heights etc. and by *Attribute data* (properties of the spatial entities). GIS stores information about the real world as a collection of thematic layers with characteristics of common coordinate system. The spatial entities of the earth's features can be represented in digital form by two data models: vector or raster models.

#### **1.1.1 Vector data model**

The method of representing geographic features by the basic graphical elements of points, lines and polygon is said to be the *vector method* or *vector data model. The r*elated vector data are always organized by *themes*, which are also referred to as *layers* or *coverages i.e. administrative* boundaries, infrastructure, soil, vegetation cover, landuse, hydrology, land parcel and others. In a vector model the position of each spatial feature i.e. a point (or node), line (or arc) and area (or polygon), is defined by the series of *x* and *y* coordinates. The

interrelationship between these features is called a topological relationship. Any change in a point, line or area will influence other factors through the topological relationship. A *label* is required so that user can load the appropriate attribute record for a given geographic feature.

#### **1.1.2 Raster data model**

The method of representing geographic features by pixels is called the *raster method* or *raster data model.* The object space is divided into a group of regularly spaced grids or pixels to which the attributes are assigned. Pixels (a term derived for a picture element) are the basic units for which information is explicitly recorded. A raster pixel represents the generalized characteristics of an area of specific size on or near the surface of the earth. The actual ground size depicted by a pixel is dependent on the resolution of the data, which may range from less than a square meter to several square kilometers. Raster data are organized by themes, which are also referred to as layers for example; a raster geographic database may contain the following themes: bed rock geology, vegetation cover, landuse, topography, hydrology, rainfall, temperature. Image data utilizes techniques very similar to raster data, however typically lacks the internal formats required for analysis and modeling of the data. The raster data model has advantages of being simple in structure, provision of easy and efficient overlaying and compatibility with RS image data while the vector data model has compact data structure, provision of efficient network analysis and projection transformation, and can generate accurate map output. The usage of the two models depends on the objectives and planning needs of the GIS project.

The most popular form of primary raster data capture is remote sensing - a technique used to derive information about the physical, chemical, and biological properties of objects without direct physical contact. Information is derived from measurements of the amount of electromagnetic radiation reflected, emitted, or scattered from objects. A variety of sensors, operating throughout the electromagnetic spectrum from visible to microwave wavelengths, are commonly employed to obtain measurements (Lillesand and Kiefer, 2004). Remote sensing provides information of natural indicators and landuse features that can be related with the characteristics of sub-surface aquifer system (recharge and discharge sources) and presence of soil moisture/shallow groundwater condition.

By integrating remote sensing with GIS, an even greater potential for environmental applications is achieved (Ehlers et al., 1989; Shelton and Estes 1980). The RS data has its advantages of spatial, spectral and temporal availability of data covering large and inaccessible areas within short time and thus proves to be an effective tool in assessing, monitoring surface and groundwater resources. The conventional techniques have the limitation to study these parameters together because of the non-availability of data, integration tools and modeling techniques. RS data provides an economic and efficient tool for landuse mapping and has its advantages in the development and management of water resources for optimum use.

#### **1.2 Numerical groundwater flow methods**

The range of numerical methods is quite large, obviously being of use to most fields of engineering and science in general. There are two broad categories of numerical methods i.e.

interrelationship between these features is called a topological relationship. Any change in a point, line or area will influence other factors through the topological relationship. A *label* is required so that user can load the appropriate attribute record for a given geographic

The method of representing geographic features by pixels is called the *raster method* or *raster data model.* The object space is divided into a group of regularly spaced grids or pixels to which the attributes are assigned. Pixels (a term derived for a picture element) are the basic units for which information is explicitly recorded. A raster pixel represents the generalized characteristics of an area of specific size on or near the surface of the earth. The actual ground size depicted by a pixel is dependent on the resolution of the data, which may range from less than a square meter to several square kilometers. Raster data are organized by themes, which are also referred to as layers for example; a raster geographic database may contain the following themes: bed rock geology, vegetation cover, landuse, topography, hydrology, rainfall, temperature. Image data utilizes techniques very similar to raster data, however typically lacks the internal formats required for analysis and modeling of the data. The raster data model has advantages of being simple in structure, provision of easy and efficient overlaying and compatibility with RS image data while the vector data model has compact data structure, provision of efficient network analysis and projection transformation, and can generate accurate map output. The usage of the two models

The most popular form of primary raster data capture is remote sensing - a technique used to derive information about the physical, chemical, and biological properties of objects without direct physical contact. Information is derived from measurements of the amount of electromagnetic radiation reflected, emitted, or scattered from objects. A variety of sensors, operating throughout the electromagnetic spectrum from visible to microwave wavelengths, are commonly employed to obtain measurements (Lillesand and Kiefer, 2004). Remote sensing provides information of natural indicators and landuse features that can be related with the characteristics of sub-surface aquifer system (recharge and discharge sources) and

By integrating remote sensing with GIS, an even greater potential for environmental applications is achieved (Ehlers et al., 1989; Shelton and Estes 1980). The RS data has its advantages of spatial, spectral and temporal availability of data covering large and inaccessible areas within short time and thus proves to be an effective tool in assessing, monitoring surface and groundwater resources. The conventional techniques have the limitation to study these parameters together because of the non-availability of data, integration tools and modeling techniques. RS data provides an economic and efficient tool for landuse mapping and has its advantages in the development and management of water

The range of numerical methods is quite large, obviously being of use to most fields of engineering and science in general. There are two broad categories of numerical methods i.e.

depends on the objectives and planning needs of the GIS project.

presence of soil moisture/shallow groundwater condition.

resources for optimum use.

**1.2 Numerical groundwater flow methods**

feature.

**1.1.2 Raster data model**

gridded or discretized methods and non-gridded or mesh-free methods. In the common finite difference method and finite element method (FEM) the domain is completely gridded (form a grid or mesh of small elements). The analytic element method (AEM) and the boundary integral equation method (BIEM - sometimes also called BEM, or Boundary Element Method) are only discretized at boundaries or along flow elements (line sinks, area sources, etc.). Gridded Methods like finite difference and finite element methods solve the groundwater flow equation by breaking the problem area (domain) into many small elements (squares, rectangles, blocks, triangles, tetrahedra, etc.) and solving the flow equation for each element (all material properties are assumed constant or possibly linearly variable within an element), then linking together all the elements using conservation of mass across the boundaries between the elements (http://en.wikipedia.org/wiki/Hydrogeology). This results in a system which overall approximates the groundwater flow equation i.e. approximating the Laplace equations for flow in a porous media and the partial differential equations governing solute transport.

Finite differences are a way of representing continuous differential operators using discrete intervals (*Δx* and *Δt*). MODFLOW is a well-known example of a general finite difference groundwater flow model. It is one of the most widely used and tested software program developed by U.S. Geological Survey for simulating groundwater flow. It is a threedimensional modular finite difference model that uses variable grid spacing to model groundwater flow. Many commercial products have grown up around it, providing graphical user interfaces to its input file based interface, and typically incorporating preand post-processing of user data. Many other models have been developed to work with MODFLOW input and output, making linked models which simulate several hydrologic processes possible (flow and transport models, surface water and groundwater models and chemical reaction models), because of the simple, well documented nature of MODFLOW. The model has been adopted by the GMS (Groundwater Modeling System) and 'Visual MODFLOW' Graphical User Interfaces. GMS has the ability to import borehole data and create a 3D visualization of the geology, the ability to import GIS files directly into the conceptual model using the 'map' module. Visual MODFLOW is also an user friendly software that has ability to generate 3D visualization graphics and import GIS data.

Finite Element programs are more flexible in design (triangular elements vs. the block elements most finite difference models use) and there are some programs available (SUTRA, a 2D or 3D density-dependent flow model by the USGS; Hydrus, a commercial unsaturated flow model; FEFLOW, a commercial modeling environment for subsurface flow, solute and heat transport processes; and COMSOL Multiphysics (FEMLAB) a commercial general modeling environment). Regional model are used for large scale systems e.g. the P.R.A.M.S. (Perth Regional Aquifer Modeling System) model of the Swan coastal plain aquifer system. Packages need to be selected on the basis of suitability for the intended modeling purpose (e.g. flow modeling, transport modeling, salt water interface modeling etc). Once the model is chosen and the preparation of the input dataset is achieved, one can then proceed for modeling the groundwater behavior that represent the user's requirement. The modeling process comprises of steady-state simulation and transient simulation of the prolific groundwater system, and predictive simulation of groundwater flow/solute transport to study different scenarios.

#### **1.3 GIS and groundwater modeling**

As the geographical location of every item of information stored in a GIS is known, GIS technology makes it possible to relate the quality of groundwater at a site with the health of its inhabitants, to predict how the vegetation in an area will change as the irrigation facilities increases, or to compare development proposals with restrictions on land use. This ability to overlay gives GIS unique power to make decisions about places and to predict the outcomes of those decisions. The ability to combine and integrate data is the backbone of GIS. S*patial modeling* technique in GIS can answer locational and quantitative questions involving, *e.g.,*  the particulars of a given location, the distribution of selected phenomena, the changes that have occurred since a previous analysis, the impact of a specific event in the form of "*Where?*  or *How much?* and/or *What if ?"* scenarios-making in a more realistic way. Some of the advantages of GIS application in groundwater modeling include the following:


Application of GIS technology to hydrological modeling requires careful planning and extensive data manipulation work. In general, the following three major steps are required:


GIS provide spatial data handling and graphical environment to analyze and visualize spatial data related to numerical groundwater flow modeling. The frequently applied GIS functions, which support groundwater modeling during its various steps, are summarized in Table 1. In fact the modeling steps are interlinked to each other and can be varied according to the planning requirements.

In preprocessing, the data is transferred from the GIS to the model while in postprocessing, the data is transferred from the model to the GIS. Developing conceptual models, choosing a computer code, and developing model designs are the preliminary steps in a groundwater modeling project after first establishing the purpose of the project (Anderson & Woessner, 1992). These steps define the numerical models representing a given field situation,

As the geographical location of every item of information stored in a GIS is known, GIS technology makes it possible to relate the quality of groundwater at a site with the health of its inhabitants, to predict how the vegetation in an area will change as the irrigation facilities increases, or to compare development proposals with restrictions on land use. This ability to overlay gives GIS unique power to make decisions about places and to predict the outcomes of those decisions. The ability to combine and integrate data is the backbone of GIS. S*patial modeling* technique in GIS can answer locational and quantitative questions involving, *e.g.,*  the particulars of a given location, the distribution of selected phenomena, the changes that have occurred since a previous analysis, the impact of a specific event in the form of "*Where?*  or *How much?* and/or *What if ?"* scenarios-making in a more realistic way. Some of the

 GIS provides decision support for groundwater management i.e. groundwater pumping for domestic, industrial and agricultural supplies and other actions influencing the regional water cycle: infrastructure, tunnels, waste dump-sites,

 GIS saves much time of collecting large number of geographical data required for groundwater modeling for both pre-processing and post-processing stages and to

GIS can handle large datasets through integration with Database Management System

 With GIS, complex maps can be created and edited much faster that would not be possible by hand, and because the data is stored digitally, the maps are produced with

 GIS can utilize a satellite image to extract useful information i.e. landuse, surface hydrology, soil properties etc. which serves as source data for model conceptualization. GIS has capability to integrate with many hydrological models and techniques, and

Application of GIS technology to hydrological modeling requires careful planning and extensive data manipulation work. In general, the following three major steps are required:

GIS provide spatial data handling and graphical environment to analyze and visualize spatial data related to numerical groundwater flow modeling. The frequently applied GIS functions, which support groundwater modeling during its various steps, are summarized in Table 1. In fact the modeling steps are interlinked to each other and can be varied

In preprocessing, the data is transferred from the GIS to the model while in postprocessing, the data is transferred from the model to the GIS. Developing conceptual models, choosing a computer code, and developing model designs are the preliminary steps in a groundwater modeling project after first establishing the purpose of the project (Anderson & Woessner, 1992). These steps define the numerical models representing a given field situation,

(DBMS) component which provides foundation for all analysis techniques.

transform spatial data according to the modeling requirements.

advantages of GIS application in groundwater modeling include the following:

**1.3 GIS and groundwater modeling**

sewerage systems etc.

improve the model results.

 Development of spatial database Extraction of model layers Linkage to computer models

according to the planning requirements.

the same level of accuracy each time.


Table 1. GIS functions involved in different phases of groundwater flow modeling

including: (a) the most important sources and sinks of water in the field system and how they are to be simulated; (b) the available data on the geohydrologic system; (c) the system geometry (generally the number and type of model layers and the areal extent of these layers); (d) the spatial and temporal structure of the hydraulic properties (generally using zones of constant value or deterministic or stochastic interpolation methods); and (e) boundary condition location and type. The groundwater models output may include: water pressure (head) distribution; flow rates, flow directions; plume movement and particle tracking; water chemistry changes and budgeting. Integration of GIS and hydrologic models follows one of the two approaches; a) To develop hydrologic models that operate within a GIS framework, b) To develop GIS techniques that partially define the parameters of existing hydrologic models (Jain et. al., 1997). In most of groundwater modeling softwares such as Feflow, Modflow, GMS there is an interface that links vector data through compatible GIS formats i.e. *.shp, .lin, .dxf* etc. and raster data formats; *.tif, .bmp, .img* etc.

The GIS data assists in identifying the spatial variability of hydraulic conductivity and recharge, the values of which are estimated under steady-state condition. The effect and sensitivity of different parameters values are tested during the modeling process through calibration technique using packages like PEST and UCODE. Calibration optimizes the input parameter values against the historical water data. The simulation output of the models can be integrated with any thematic data in GIS for verification, impact and/or characterization analysis. There are inbuilt presentation tools in the models that can be used to create labeled contour maps of input data and simulate results. One can fill colors to model cells containing different values and report-quality graphics may be saved to a wide variety of file formats i.e. *surfer*, *dxf, hpgl* and *bmp*. The presentation tools can even create and display two dimensional animation sequences using the simulation results (calculated heads, drawdowns or concentration).

#### **2. Spatial analysis of waterlogged areas**

The study area lies between longitudinal range 73-74 5E and latitudinal range 32- 32 45N in Indus basin, Pakistan (Fig. 1). It covers an area of about 3,417 sq. km within Rivers Jhelum and Chenab flowing in the northwest and southeast, and Upper Jhelum Canal (UJC) and Lower Jhelum Canal (LJC) in the northeast and southwest. The area is well connected with other main cities of the country through metalled roads and rail links. There lie three schemes of Salinity Control and Reclamation Project SCARP-II of Water and Power Development Authority (WAPDA) i.e. Sohawa, Phalia and Busal in the area. The relief increases northward with elevation ranges between 200 and 238 m above mean sea level. The climate of the area is mainly semi arid. May, June and July are the hottest months and the mean maximum and minimum temperature during this period are about 39.5C and 25.4C respectively (Qureshi et al. 2003). The coldest months are December, January and February and the maximum and minimum temperatures during this period are about 21.5C and 5.1C respectively. The rains are erratic and are received in two rainy seasons i.e. about two third of the total rain is received during monsoon season (mid June to mid September) while remaining one third, during the period from January to March.

Major landforms of the area include piedmont plain and basin, flood plains, scalloped interfluves and channel levee remnants (Fig. 1). Waterlogging is a problem in the central parts of the scalloped interfluves and in some parts of the channel-levee remnants. The meander flood plain has nearly level surface and contains local depressions which causes drainage problem. The Piedmont basin contains clayey soils characterized by swelling and shrinking. During the rainy season runoff water from the adjoining land collects in the area, making these soils waterlogged (Soil Survey Report 1967). The soils of the major part of study area are formed from alluvial sediments brought by the rivers from the Himalayan mountain ranges in the north. The alluvial complex forms a unified, highly permeable aquifer, in which water is generally unconfined. The groundwater flows in general from northeast to southwest of study area with the hydraulic gradient ranging from 2.14 to 7.15 ft/mile (PPSGDP 2000). Drainage pattern is mainly dendritic in nature. The waterlogging has affected the main irrigated area of the Indus basin (Fig. 2). The main factors involved in causing the problem of waterlogging and salinity in the area include the following:

 Micro relief in the land surface resulting in appearance of waterlogging in patches, which act as sinks


sensitivity of different parameters values are tested during the modeling process through calibration technique using packages like PEST and UCODE. Calibration optimizes the input parameter values against the historical water data. The simulation output of the models can be integrated with any thematic data in GIS for verification, impact and/or characterization analysis. There are inbuilt presentation tools in the models that can be used to create labeled contour maps of input data and simulate results. One can fill colors to model cells containing different values and report-quality graphics may be saved to a wide variety of file formats i.e. *surfer*, *dxf, hpgl* and *bmp*. The presentation tools can even create and display two dimensional animation sequences using the simulation results (calculated

The study area lies between longitudinal range 73-74 5E and latitudinal range 32- 32 45N in Indus basin, Pakistan (Fig. 1). It covers an area of about 3,417 sq. km within Rivers Jhelum and Chenab flowing in the northwest and southeast, and Upper Jhelum Canal (UJC) and Lower Jhelum Canal (LJC) in the northeast and southwest. The area is well connected with other main cities of the country through metalled roads and rail links. There lie three schemes of Salinity Control and Reclamation Project SCARP-II of Water and Power Development Authority (WAPDA) i.e. Sohawa, Phalia and Busal in the area. The relief increases northward with elevation ranges between 200 and 238 m above mean sea level. The climate of the area is mainly semi arid. May, June and July are the hottest months and the mean maximum and minimum temperature during this period are about 39.5C and 25.4C respectively (Qureshi et al. 2003). The coldest months are December, January and February and the maximum and minimum temperatures during this period are about 21.5C and 5.1C respectively. The rains are erratic and are received in two rainy seasons i.e. about two third of the total rain is received during monsoon season (mid June to mid

September) while remaining one third, during the period from January to March.

causing the problem of waterlogging and salinity in the area include the following:

Micro relief in the land surface resulting in appearance of waterlogging in patches,

Major landforms of the area include piedmont plain and basin, flood plains, scalloped interfluves and channel levee remnants (Fig. 1). Waterlogging is a problem in the central parts of the scalloped interfluves and in some parts of the channel-levee remnants. The meander flood plain has nearly level surface and contains local depressions which causes drainage problem. The Piedmont basin contains clayey soils characterized by swelling and shrinking. During the rainy season runoff water from the adjoining land collects in the area, making these soils waterlogged (Soil Survey Report 1967). The soils of the major part of study area are formed from alluvial sediments brought by the rivers from the Himalayan mountain ranges in the north. The alluvial complex forms a unified, highly permeable aquifer, in which water is generally unconfined. The groundwater flows in general from northeast to southwest of study area with the hydraulic gradient ranging from 2.14 to 7.15 ft/mile (PPSGDP 2000). Drainage pattern is mainly dendritic in nature. The waterlogging has affected the main irrigated area of the Indus basin (Fig. 2). The main factors involved in

heads, drawdowns or concentration).

which act as sinks

**2. Spatial analysis of waterlogged areas**

Lower evaporation than recharge resulting in appearance of swamps

Major contribution to the recharge of the groundwater is from seepage of the rivers, canal irrigation network including main and link canals, distributaries, minors and watercourses/fields, precipitation and return flows of the groundwater use, besides subsurface inflows from nearby zones. The presence of large canal network in the area provides the main source of recharge to the groundwater. Aquifer discharge sources are pumpage from public and private tubewells (water wells), evapotranspiration, outflows to the rivers and drains. The groundwater is mainly abstracted to augment where or when the canal supplies are not adequate especially in the Rabi (wet) season.

The wastelands comprising waterlogged areas, swamps and marshy land along Jhelum and Chenab riverbeds demarcated through RS analysis were related with associated factors like landforms, canal system, topography, climate and groundwater etc.

Fig. 1. Location and landforms extent in the study area

The following thematic layers were developed through on-screen digitizing in ILWIS 3.1 GIS software:


The attribute database of each layer was developed and linked to its respective data map for GIS analysis. The spatial functions of GIS were used to derive the following thematic maps for spatial analysis;


Fig. 2. Effect of waterlogging on irrigated area

The landcover/landuse map was developed through image processing of remote sensing data using ERDAS imagine 8.5 software. These maps including the landcover map have been registered together in GIS environment for integration and analysis (Fig. 3). The ancillary data consists of topographic maps of scale 1:50,000, geomorphology & soil maps, hydrogeology and groundwater levels of observation wells and climatic data.

The landcover map in raster form was converted into vector form and a layer of wasteland polygons was developed which comprises of 585 polygons. Major waterlogged areas were selected by eliminating smaller polygons having an area less than 0.02 km2 leaving total of 361 polygons which covered an aggregate area of about 71 sq. km. The polygon data was analyzed with different thematic data through spatial modeling in GIS. Most of the waterlogged areas (about 28%) were found in the active flood plain followed by about 17% and 13% in the scalloped interfluves and young channel levee remnants. The low-lying

The attribute database of each layer was developed and linked to its respective data map for GIS analysis. The spatial functions of GIS were used to derive the following thematic maps

The landcover/landuse map was developed through image processing of remote sensing data using ERDAS imagine 8.5 software. These maps including the landcover map have been registered together in GIS environment for integration and analysis (Fig. 3). The ancillary data consists of topographic maps of scale 1:50,000, geomorphology & soil maps,

The landcover map in raster form was converted into vector form and a layer of wasteland polygons was developed which comprises of 585 polygons. Major waterlogged areas were selected by eliminating smaller polygons having an area less than 0.02 km2 leaving total of 361 polygons which covered an aggregate area of about 71 sq. km. The polygon data was analyzed with different thematic data through spatial modeling in GIS. Most of the waterlogged areas (about 28%) were found in the active flood plain followed by about 17% and 13% in the scalloped interfluves and young channel levee remnants. The low-lying

hydrogeology and groundwater levels of observation wells and climatic data.

Geomorphology

 Infrastructure Surface hydrology Cities and towns Administrative units

for spatial analysis; Surface elevation

Fig. 2. Effect of waterlogging on irrigated area

 Slope Buffer zones Theissen Polygons Watertable depth Recharge zones Landcover/Landuse

Soils

Fig. 3. GIS has capability to integrate different types of spatial data

areas i.e. parts of the piedmont basin and central part of the scalloped interfluves have especially become victims of waterlogging and salinity (Soil Survey Report 1967). The meander flood plain constitutes about 12% of waterlogged areas. About 11% of waterlogged areas occur in the Piedmont basin which is actually the lowest end of the piedmont plain where, due to absence of slope, basin conditions have developed. The braided riverbed constitutes about 9 km2 of waterlogged areas mostly comprising marshy areas stretched along the Chenab River and to some extent along the Jhelum River. The buffer zones (polygons of 500m interval) were created around the main canal system to analyze the influence of seepage on the waterlogged areas. About 21 sq. km of waterlogged area lie within the buffer zone of 0-500m while only 4.4 sq. km within zone of 1000-1500m indicating a decrease in coverage with the increase in distance from the canal network (Fig. 4). In fact the canal irrigation is a major source of recharge to the groundwater in this area. The watertable beneath the waterlogged areas varies temporally and spatially. During the year 1990, about 4% waterlogged areas lies in zone <1.5m depth whereas about 28% in 1.5-3.0m depth. The watertable zone of 3-6m depth stretched over about 65% area, comprising 63% of the waterlogging polygons. This shows a local effect of the topography and lithological conditions in developing of the waterlogged areas. Analysis of the waterlogged areas with mean annual rainfall indicated maximum of 154 polygons under rainfall range of 500- 600 mm followed by 127 polygons in less than 500mm rainfall range.

The development of a sufficient understanding on which to base decisions or make predictions often requires consideration of a multitude of data of different types and with varying levels of uncertainty (Middlemis 2000). The data for the development of numerical groundwater flow model included the following parameters:

Fig. 4. Spatial analysis of the waterlogged areas using buffering technique of GIS


The time-constant parameters were mainly extracted from thematic data layers generated through GIS and image processing of remote sensing data. The remote sensing data of LANDSAT TM & ETM sensors was used of periods 1990 and 2001. The general approach involves processing of RS image data for analysis of different landcovers and features present in the study area and integration of results in GIS system along with results generated through numerical groundwater flow modeling. The image data was geometrically rectified and analyzed through visual and digital interpretation for the study of landcover/landuse types. The former interpretation was used for qualitative analysis while the latter for quantitative analysis of the landcover/landuse in the RS image data.

Fig. 4. Spatial analysis of the waterlogged areas using buffering technique of GIS

Number, distribution and pumpage from irrigation/drainage tubewells

 Aquifer geometry (areal and vertical distribution of subsurface strata, aquifer thickness etc.). It also requires river network, landcover, soil, surface and subsurface hydrological

The time-constant parameters were mainly extracted from thematic data layers generated through GIS and image processing of remote sensing data. The remote sensing data of LANDSAT TM & ETM sensors was used of periods 1990 and 2001. The general approach involves processing of RS image data for analysis of different landcovers and features present in the study area and integration of results in GIS system along with results generated through numerical groundwater flow modeling. The image data was geometrically rectified and analyzed through visual and digital interpretation for the study of landcover/landuse types. The former interpretation was used for qualitative analysis while the latter for quantitative analysis of the landcover/landuse in the RS image data.

*a. Time Constant Parameters* 

 Hydro-meteorological data Water level monitoring data

River and canal flows and hydraulic features

data input. Hydraulic parameters *b. Time Variant Parameters*  Initially, unsupervised classification method was used to analyze the statistical patterns of various clusters in the image. Later, supervised method of classification was used which is based on background knowledge and experience of the interpreter, available ground information and ancillary data. The surface hydrology including canal irrigation network, streams and drains were digitized partially from RS image data and topographic maps. The statistics of this hydrological network is shown in Table 2. Main canals i.e. UJC and LJC extends upto 62 km and the distributaries mainly of UJ Canal extends upto 398 km in the study area. There are several large streams (about 96 km in length) that originate from the Himalayan Mountains in the north and drain into the Chenab River in the Southeast. The analytical functions of ILWIS, Arcview and ArcGIS softwares were mainly used for digitization, analysis and visualization of the spatial data in the present study.


Table 2. Statistical analysis of hydrological network in the study area

Through visual interpretation of the image, distribution of landcover and features like surface drainage and irrigation network, water bodies, waterlogged areas and soil moisture condition were analyzed. Landforms and structures that can reveal subsurface condition of groundwater were also studied. The false color composite (FCC) of Landsat bands 7, 4, 2 (Red, Green, Blue) showed waterlogged areas in true colors i.e. dark bluish-green spots. The dark appearance of these areas is due to low reflectivity of high moisture contents of standing water consisting of evergreen natural vegetation grown in waterlogged areas. Waterlogged areas are differentiated from fresh water in bands 4 and 5 due to presence of natural vegetation in it. Moist soil can be differentiated from the swamps and waterlogged area in bands 3, 5 and 7, which are indicating higher reflection of moist soil in these bands. The image classification of the RS data indicated about 70% of the land under various types of crops followed by 8% grassland and 6% forest cover. About 10% land is bare soil and 4% is waterlogged and swampy area.

#### **2.1 Development of numerical groundwater flow model**

Numerical modeling employs approximate methods to solve the partial differential equation (PDE) which describes the flow in porous medium. The emphasis here is not on obtaining an exact solution but on obtaining reasonably approximate solution (Thangarajan, 2004). Feflow is a fully integrated modeling environment with a full-featured graphical interface and powerful numeric engines that allow the user to perform any flow or contaminant transport modeling. The components ensure an efficient process for building the finite element model, running the simulation, and visualizing the results. The conceptual model of the study area was developed for carrying out finite element groundwater flow modeling. Based on the hydrogeological conceptual model of the study area, the mathematical model was established as follows:

$$\frac{\partial}{\partial \mathbf{x}} \left( \text{Klh} \frac{\partial H}{\partial \mathbf{x}} \right) + \frac{\partial}{\partial y} \left( \text{Klh} \frac{\partial H}{\partial y} \right) + \frac{\partial}{\partial z} \left( \text{Klh} \frac{\partial H}{\partial z} \right) + \varepsilon = \mu \frac{\partial H}{\partial t} \tag{1}$$

Initial conditions:

$$\left. \left. H \left( \mathbf{x}, y, z, t \right) \right|\_{t=0} = H\_0 \left( \mathbf{x}, y, z \right) \text{ (x, y, z)} \in D \tag{2}$$

Boundary conditions:

$$\left. \mathcal{K} \frac{\partial \mathcal{H}}{\partial \mathbf{n}} \right|\_{\Gamma\_2} = q\left(\mathbf{x}, y, z, t\right) \text{ (\$\mathbf{x}\$, \$y\$, z\$)} \in \Gamma\_2 \quad t \ge 0 \tag{3}$$

Where *K* is permeability coefficient (or hydraulic conductivity coefficient) (m/d) (Due to the lack of experimental data, the aquifer is taken as isotropic); *h* is the distance from the impervious bed of the aquifer to the free water surface, i.e. aquifer thickness (m); *H* is water head (or water table) (m); *ε* is the inflow and outflow factors, i.e. the volume of water that vertically flows into or out of the aquifer in unit time and unit area with inflow being positive and outflow being negative (m3/d/m2); *μ* is storativity; *t* is time (d); *H*0 is initial value of water head (m); *D* is study area enclosed by Г2, and Г2 is the second kind boundary; *n* is the direction of outer normal line of the second kind boundary; and *q* is the volume of water that laterally flows into or out of the aquifer in unit time and unit area on the second kind boundary (m/d).

The finite elements may have both different spatial dimensions and shapes. The order of the underlying interpolation scheme may typically be linear, quadratic or cubic. Continuity may be prescribed not only for the variable themselves but also for their derivatives. The procedures to be followed by UNESCO (1999) are given below;


$$h(\mathbf{x}, \mathbf{y}, z) = \sum\_{j=1}^{N} h\_i \mathbf{y}\_j \tag{4}$$

where '*h*' is piezometric head, *(x, y, z)* are Cartesian coordinates, *'N'* is the number of nodal points in the discretized element grid and '' is the interpolation function or otherwise called shape function.


$$\mathbf{[M]} \{ \frac{dh}{dt} \} + \{ \mathbf{K}(\mathbf{h}) \} \begin{bmatrix} \mathbf{h} \end{bmatrix} = \mathbf{[f]} \tag{5}$$

Where *[M]* denotes a matrix.

vi. Time integration.

252 Water Resources Management and Modeling

2004). Feflow is a fully integrated modeling environment with a full-featured graphical interface and powerful numeric engines that allow the user to perform any flow or contaminant transport modeling. The components ensure an efficient process for building the finite element model, running the simulation, and visualizing the results. The conceptual model of the study area was developed for carrying out finite element groundwater flow modeling. Based on the hydrogeological conceptual model of the study area, the

> *H H HH Kh Kh Kh x xy yz z t*

Where *K* is permeability coefficient (or hydraulic conductivity coefficient) (m/d) (Due to the lack of experimental data, the aquifer is taken as isotropic); *h* is the distance from the impervious bed of the aquifer to the free water surface, i.e. aquifer thickness (m); *H* is water head (or water table) (m); *ε* is the inflow and outflow factors, i.e. the volume of water that vertically flows into or out of the aquifer in unit time and unit area with inflow being positive and outflow being negative (m3/d/m2); *μ* is storativity; *t* is time (d); *H*0 is initial value of water head (m); *D* is study area enclosed by Г2, and Г2 is the second kind boundary; *n* is the direction of outer normal line of the second kind boundary; and *q* is the volume of water that laterally flows into or out of the aquifer in unit time and unit area on the second

The finite elements may have both different spatial dimensions and shapes. The order of the underlying interpolation scheme may typically be linear, quadratic or cubic. Continuity may be prescribed not only for the variable themselves but also for their derivatives. The

i. Discretization of the flow domain into a set of elements, where each element is defined by a number of nodes, for instance 3 or 6-node triangles, 4- 8- or 9-node quadrilaterals. ii. Expression of field parameters such as piezometric head, hydraulic conductivity etc., in

(,,)

1

*j hxyz h*

where '*h*' is piezometric head, *(x, y, z)* are Cartesian coordinates, *'N'* is the number of nodal

*i j*

(4)

' is the interpolation function or otherwise

*N*

2 ,,, *<sup>H</sup> <sup>K</sup> <sup>q</sup> <sup>x</sup> <sup>y</sup> z t*

*n*

procedures to be followed by UNESCO (1999) are given below;

 

(1)

t=0 <sup>0</sup> *H xyzt H xyz* ,,, ,, *x*, , *<sup>y</sup> z D* (2)

<sup>2</sup> *xyz* , , <sup>0</sup> *<sup>t</sup>* (3)

mathematical model was established as follows:

Initial conditions:

Boundary conditions:

kind boundary (m/d).

the following form:

called shape function.

points in the discretized element grid and '

vii. Solution of the global system of linear equations.

The hydrogeological system of the site was modeled as multi-layered using Finite element Model - Feflow. A model grid consisting of superelement mesh of five elements was drawn over the model area using the background information of landuse, landforms and drainage/canal network of the area (Fig. 5). The superelement mesh represents the basic structure of the study domain. The Finite element mesh was generated from the superelement mesh using triangulation option of 6-noded prism for 3-dimensional model. The 3-D mesh consists total of 5,343 elements and 3,928 nodes (Fig. 6). The model layers were developed from point data using Akima's bivariate interpolation method. It is difficult to introduce precise distribution of recharge over the area due to complex distribution of various types of soils in the area. However, equal distribution of the recharges from various sources, on macro level, generally gives the results within practical limits (Thangarajan, 2004). A simple water balance model of an area can be carried out on the fact that there exists a balance between the quantity of water entering into the area, amount store and water leaving the same area during certain period of time. For any hydrologic system, general mass conservation equation can be written as:

$$I - O = \Lambda S \tag{6}$$

Where *I* = Total inflows, *O* = Total outflows and *S* = change in groundwater storage.

Since flow in saturated zone is simulated using groundwater flow model, the net recharge to the groundwater reservoir was computed by assuming aggregate water balance for the unsaturated zone and the land surface. The net recharge of the model area was thus estimated as below:

$$\text{Qnet} = \text{Qmc} + \text{Qdm} + \text{Qwf} + \text{Qrv} + \text{Qrf} + \text{Qrtw} - (\text{Qet} - \text{Qtr} - \text{Qdr}) \tag{7}$$

Where *Qne*t = net recharge to the aquifer, *Qmc* = Recharge from main and link canals, *Qdm* = recharge from distributary/minors, *Qwf* = recharge from watercourses and irrigated fields, *Qriv* = recharge/discharge from rivers, *Qrf* = recharge from rainfall, *Qrtw* = recharge from return flow of tubewell pumpage, *Qet*=discharge by evapotranspiration, *Qtw* = discharge by tubewells, *Qdr* = discharge from drains. The recharge from inflows from adjacent areas was considered implicitly in the groundwater modeling (Sarwar, 1999). The losses of the canal system and other recharge sources were maintained according to the specified limits of the irrigation department and WAPDA (PPSGDP 2000 & WAPDA, 1993).

Fig. 5. Delineation of superelements and recharge zones for groundwater flow modeling in Feflow

Fig. 6. Finite element mesh and point data of observation wells drapped over the soil map

#### **2.2 Model calibration**

254 Water Resources Management and Modeling

Fig. 5. Delineation of superelements and recharge zones for groundwater flow modeling in

Fig. 6. Finite element mesh and point data of observation wells drapped over the soil map

Feflow

The groundwater levels of June 1985 (pre-monsoon period) of 28 observation wells were used as initial condition for executing steady-state simulation. Steady-state simulations were carried out to describe three-dimensional flow field (head distribution) over the entire model domain. Automatic parameter estimation (PEST) method was applied for calibration of the steady-state model (Doherty 1995). The hydraulic conductivity values taken from test holes data were used to develop conductivity zones using Theissen polygons for model calibration. In order to estimate the recharge of the model domain, the model area was divided into five recharge zones (shown in Fig. 5) on the basis of hydrological setup, geomorphology and land capability of the area. The recharge of the zones was characterized by variable infiltration rates of different soils and varying pumpage of groundwater. Initially, steady-state calibration was performed which was fully implicit. The hydraulic conductivity and the recharge values estimated previously were used as initial conditions in the steady-state calibration. The values were adjusted during the calibration runs until the calculated head values became close to the observed heads. Similarly, the specific yield zones were developed using its field data for use in transient state-calibration. The model was rerun for six-month period i.e. April-September 1985, for transient-state calibration. The mean residuals of observed and calculated heads in steady-state and transient-state calibrations are 0.06m and 0.002m with variances of 1.46m and 1.86m, respectively. The calibration results indicated a reasonable agreement between the calculated and observed heads (Fig.7). The sensitivity of the model results was evaluated to variations in hydrologic parameters and modeling assumptions. The errors are expressed as equation 8.

$$RMSE = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left( H\_i - H\_i' \right)^2} \tag{8}$$

Where *RMSE* is root mean square error (m); *n* is total number of measurements; *Hi* is simulated value of groundwater table at the end of *i*th month (m); and *Hi*' is observed value of groundwater table at the end of *i*th month (m). The flow diagram of the methodology followed in the present study is shown in Fig. 8.

Fig. 7. a. Steady-state calibration b. Transient-state calibration

Fig. 8. Flow diagram of integrating Groundwater flow modeling and GIS

#### **2.3 Model simulation for waterlogging analysis**

Strategy of management for pre-stress period and post stress period was developed on the basis of availability of observed data until 2005 and projected hypothetical data for selected periods from 2006 to 2020. The steady-state calibrated model was run for pre-stress period of variable time steps until 2005 (Fig. 9). During the calibrated period of 1985-2005, the watertable had shown an average decline of 0.96m at the rate of 0.05 m/year. The water budget of the steady-state model shows total flux-in of the model area of about 4.6036E+05 m3/d including boundary inflows of 1.6931E+05 m3/d and areal recharge of 2.9105E+05 m3/d. The total flux-out from the model area is estimated as -4.6028E+05 m3/d. The watertable indicated a rising trend from 1988 to 1999 followed by a gradual decline onward. The initial rise may be attributed to record rainfalls that had occurred in Muzaffarabad and Jhelum during period 1997-98 (Siddiqi, 1999). Those rainfalls exaggerated the problem of waterlogging in central parts of the study area. A major breakthrough of groundwater depletion was observed during year 1999 when the last drought prevailed for over 3-4 years in that area. The declining trend of groundwater levels continued in the remaining calibration period. During the drought period, there became shortages in canal water supplies that resulted in low recharge from canals seepage and abstraction of groundwater for irrigation use. The situation had affected the extent of waterlogging which was reduced significantly in different parts of the area. The analysis of velocity variations of groundwater flow in three layers of the model indicated high velocity (>0.12 m/day) in the northwestern part near Rasul Barrage in the first layer. The velocity range 0.04-0.08 m/day was dominant in the central and southeastern parts of the model area. The patches of low velocity zone

Fig. 8. Flow diagram of integrating Groundwater flow modeling and GIS

Strategy of management for pre-stress period and post stress period was developed on the basis of availability of observed data until 2005 and projected hypothetical data for selected periods from 2006 to 2020. The steady-state calibrated model was run for pre-stress period of variable time steps until 2005 (Fig. 9). During the calibrated period of 1985-2005, the watertable had shown an average decline of 0.96m at the rate of 0.05 m/year. The water budget of the steady-state model shows total flux-in of the model area of about 4.6036E+05 m3/d including boundary inflows of 1.6931E+05 m3/d and areal recharge of 2.9105E+05 m3/d. The total flux-out from the model area is estimated as -4.6028E+05 m3/d. The watertable indicated a rising trend from 1988 to 1999 followed by a gradual decline onward. The initial rise may be attributed to record rainfalls that had occurred in Muzaffarabad and Jhelum during period 1997-98 (Siddiqi, 1999). Those rainfalls exaggerated the problem of waterlogging in central parts of the study area. A major breakthrough of groundwater depletion was observed during year 1999 when the last drought prevailed for over 3-4 years in that area. The declining trend of groundwater levels continued in the remaining calibration period. During the drought period, there became shortages in canal water supplies that resulted in low recharge from canals seepage and abstraction of groundwater for irrigation use. The situation had affected the extent of waterlogging which was reduced significantly in different parts of the area. The analysis of velocity variations of groundwater flow in three layers of the model indicated high velocity (>0.12 m/day) in the northwestern part near Rasul Barrage in the first layer. The velocity range 0.04-0.08 m/day was dominant in the central and southeastern parts of the model area. The patches of low velocity zone

**2.3 Model simulation for waterlogging analysis**

(<0.02 m/day) were found in the northern, northeastern and western parts which extend downward in other layers also. In the second layer, velocity zone 0.02-0.04 m/day dominated in most of the central parts. In the third layer, the velocity of less than 0.02 m/day was dominated in most of the northeastern parts of the model area. The regional groundwater flow component in the southern part is indicating existence of a potential aquifer zone in this layer of the model area.

Fig. 9. Integration of spatial data layers in GIS for analysis and output visualization

The calibrated model was rerun to predict the future changes in piezometric heads from period 2006 to 2020. Time series records of previously observed data of precipitation, annual recharge and withdrawals from tubewells were examined which formed basis to generate projected data for groundwater flow simulation. The predictive period of 15 years i.e. 2006- 2020, was chosen to perceive long-term impact of droughts/floods on regional groundwater system. Based on the annual incremental increase/decrease in recharge and/or groundwater pumpage, recharge for various stress periods was adjusted for calibration of groundwater flow model. The predictive model showed an average decline of 0.81m per year in watertable. The waterlogging indicated initially an increase in coverage from year 1990 to 2000 and than a gradual decrease from year 2000 to 2020 (Fig. 10). During pre-stress and post-stress periods, variation in head values ranges between 196 and 234 meters above sea level (masl). In upper reaches of the model domain, fluctuation in heads is low due to presence of less extensive alluvial deposits in piedmont plain. Quantitative analysis of the groundwater aquifer was performed for the base year 2005 using geo-processing techniques of GIS software. Figure 11 shows variation in watertable depth during 2005 in different tehsils (sub-districts) of Mandi Bahauddin and Gujrat districts. Results indicated maximum head drop of about 21m in the Phalia tehsil, while minimum head drop of 14m in the Gujrat tehsil. High variation in groundwater velocity was observed at several locations in the model area (Fig. 12). Maximum range of velocity from 0.006 to 0.09 m/day is found in Mandi Bahauddin and minimum range from 0.003 to 0.035 m/day in the Kharian tehsil. Overall there was relatively low variation in heads observed in the Gujrat district. The comparison in coverage of watertable depth of base year 2005 and predictive year 2020 indicated a noticeable reduction in waterlogged areas i.e. a decrease in waterlogged area (<1.5m watertable depth) from 226.9 km2 in the base year 2005 to 74.8 km2 in the year 2020 (Table 3). The changes are significant in Malakwal and Phalia tehsils. Similarly, the coverage of 1.5-3.0m watertable depth range indicated a decrease from 1299.4 km2 in 2005 to 1104.8 km2 in the year 2020. The coverage of less than 1.5m watertable depth and to some extent of 1.5m - 3.0m watertable depth had shown increase in extents during 1990-2000 period and than gradual decline from 2000 onward indicating shrinking of waterlogging extent in the predictive year 2020. This condition may be caused by natural factors i.e. low precipitation, reduce river flows, and/or human factors like over exploitation of groundwater for agriculture and domestic use in future.

Fig. 10. Development of watertable depth maps for change analysis during 1990-2020

of GIS software. Figure 11 shows variation in watertable depth during 2005 in different tehsils (sub-districts) of Mandi Bahauddin and Gujrat districts. Results indicated maximum head drop of about 21m in the Phalia tehsil, while minimum head drop of 14m in the Gujrat tehsil. High variation in groundwater velocity was observed at several locations in the model area (Fig. 12). Maximum range of velocity from 0.006 to 0.09 m/day is found in Mandi Bahauddin and minimum range from 0.003 to 0.035 m/day in the Kharian tehsil. Overall there was relatively low variation in heads observed in the Gujrat district. The comparison in coverage of watertable depth of base year 2005 and predictive year 2020 indicated a noticeable reduction in waterlogged areas i.e. a decrease in waterlogged area (<1.5m watertable depth) from 226.9 km2 in the base year 2005 to 74.8 km2 in the year 2020 (Table 3). The changes are significant in Malakwal and Phalia tehsils. Similarly, the coverage of 1.5-3.0m watertable depth range indicated a decrease from 1299.4 km2 in 2005 to 1104.8 km2 in the year 2020. The coverage of less than 1.5m watertable depth and to some extent of 1.5m - 3.0m watertable depth had shown increase in extents during 1990-2000 period and than gradual decline from 2000 onward indicating shrinking of waterlogging extent in the predictive year 2020. This condition may be caused by natural factors i.e. low precipitation, reduce river flows, and/or human factors like over exploitation of groundwater for

Fig. 10. Development of watertable depth maps for change analysis during 1990-2020

agriculture and domestic use in future.

Fig. 11. Analysis of the groundwater behavior in different tehsils during base year 2005

Fig. 12. 3-D view of velocity variations during 2005


Table 3. Change in area coverage of different watertable depths for 2005 and 2020

#### **2.4 Scenarios of extreme conditions**

In order to observe impact of natural and human induced factors on groundwater levels and ultimately on waterlogging behavior of the study area, various scenarios were developed. In the first scenario, severe drought condition was assumed to prevail for five-year period 2006-2010. The rainfall data of the severe drought that occurred during 1999-2002 was used as reference in developing the scenario. The hypothetical simulation indicated a mean decline of 0.7m in watertable depth from that of the base year 2005 (Table 4). In the second scenario, extreme wet condition was assumed to prevail for a period of five years from 2011 to 2015. It showed a mean decline of only 0.42m in watertable from that of the base year. High rainfall years will have adverse effects especially in the river flood plains and depression areas having poor drainage system. This condition may exaggerate the problem of waterlogging in the central and southern parts of the area. For effective control of water logging and minimizing the risk of this menace, management strategies and measures like installation of tubewells for vertical drainage; construction of subsurface-drains and tiledrains; planning and designing of future canals on proper lines; lining of water channels and optimizing water-use need to be adopted. Scenarios of groundwater pumpage from deep tubewells were developed to see impact of groundwater withdrawals on watertable depth. The increase in wells from 33 to 66 numbers (double) had shown decline in watertable from depth 0.14m to 0.42m (nearly 3 folds). Similarly 60% increase in rate of pumpage (from 5000m3/day to 8000m3/day) indicated increase in watertable depth to 0.25m for case 4 and 0.71m for case 6. The overall situation of groundwater pumpage indicates linear response of watertable depth with increase in numbers of wells and their rates of discharge. This strategy may help in mitigating the risk of waterlogging in the study area.


Table 4. Watertable behavior in different scenarios (WT depth in 2005 =2.95m)

#### **3. Conclusions**

260 Water Resources Management and Modeling

**6.0-**

**M.B.Din** 0.0 119.3 482.2 206.3 0.0 0.0 25.5 473.8 308.5 0.0 **Malakwal** 60.5 114.4 344.9 128.2 10.4 0.0 39.3 222.8 308.6 87.7 **Phalia** 125.1 781.2 200.6 51.5 1.0 39.7 775.4 262.9 75.8 5.6 **Gujrat** 0.0 223.2 74.0 12.0 0.0 0.0 201.4 93.8 14.0 0.0 **Kharian** 41.3 61.2 98.2 79.2 48.2 34.1 63.2 97.7 81.0 52.0 **Total** 226.9 1299.4 1199.8 477.2 59.7 73.8 1104.8 1151.0 788.0 145.4

In order to observe impact of natural and human induced factors on groundwater levels and ultimately on waterlogging behavior of the study area, various scenarios were developed. In the first scenario, severe drought condition was assumed to prevail for five-year period 2006-2010. The rainfall data of the severe drought that occurred during 1999-2002 was used as reference in developing the scenario. The hypothetical simulation indicated a mean decline of 0.7m in watertable depth from that of the base year 2005 (Table 4). In the second scenario, extreme wet condition was assumed to prevail for a period of five years from 2011 to 2015. It showed a mean decline of only 0.42m in watertable from that of the base year. High rainfall years will have adverse effects especially in the river flood plains and depression areas having poor drainage system. This condition may exaggerate the problem of waterlogging in the central and southern parts of the area. For effective control of water logging and minimizing the risk of this menace, management strategies and measures like installation of tubewells for vertical drainage; construction of subsurface-drains and tiledrains; planning and designing of future canals on proper lines; lining of water channels and optimizing water-use need to be adopted. Scenarios of groundwater pumpage from deep tubewells were developed to see impact of groundwater withdrawals on watertable depth. The increase in wells from 33 to 66 numbers (double) had shown decline in watertable from depth 0.14m to 0.42m (nearly 3 folds). Similarly 60% increase in rate of pumpage (from 5000m3/day to 8000m3/day) indicated increase in watertable depth to 0.25m for case 4 and 0.71m for case 6. The overall situation of groundwater pumpage indicates linear response of watertable depth with increase in numbers of wells and their rates of discharge.

**2020 Area (Km2)** 

> **6.0- 12.0m**

**12.0m < 1.5m 1.5-3.0m 3.0-4.5m 4.5-6.0m** 

**Period** 

**Net change in Watertable depth (m)** 

**2005 Area (Km2)** 

> **3.0- 4.5m**

**4.5- 6.0m** 

Table 3. Change in area coverage of different watertable depths for 2005 and 2020

This strategy may help in mitigating the risk of waterlogging in the study area.

Table 4. Watertable behavior in different scenarios (WT depth in 2005 =2.95m)

1 Severe drought prevails for 5-year period 2006-2010 0.70 2 Wet condition prevails for 5-year period 2011-2015 0.42 3 Pumpage from 33 TWs for 3-year period 2006-2008 0.14 4 60% increase in pumpage from 33 TWs for 3-year 2006-2008 0.25 5 Pumpage from 66 TWs for 3-year period 2006-2008 0.41 6 60% increase in pumpage from 66 TWs for 3-year 2006-2008 0.71

**Scenario Description Simulation** 

**< 1.5m 1.5-3.0m**

**2.4 Scenarios of extreme conditions**

**Tehsil** 

Rapidly developing computer technology has continued to improve modeling methods in hydrology and water resource management. GIS application has provided help in an accurate and manageable way of estimating model input parameters, integration of disparate data layers, conceptualizing of model recharge and discharge sources and visualization of the model output. GIS-based modeling, as a side benefit, also provides an updated database that can be used for non-modeling activities such as water resource planning and facilities management. The characteristics and causative factors of waterlogging/salinity can be investigated not only through remote sensing data that provide a rapid and accurate assessment of the land and water resources but also through numerical groundwater modeling which provides insight of the controlling groundwater behavior. The RS data of high spectral and spatial resolutions would be helpful in reliable assessment of landcover/landuse and identification of potential recharge/discharge areas that could ultimately enhance the quality of groundwater modeling results. The developed model would thus provides a decision support tool for evaluating better management options for sustainable development of land, surface and groundwater resources on micro as well as on macro levels in future.

#### **4. Acknowledgments**

We thank Dr. Alan Fryar Associate professor and Director of graduate study, Department of Earth and Environmental Sciences, University of Kentucky, USA and Dr. Gulraiz Akhtar, Department of Earth Sciences of Quiad-i-Azam University, Islamabad for their valuable input for execution of this work. The data support by Water and Power Development Authority (WAPDA), Pakistan Meteorological Department (PMD) and National engineering services Pakistan (NESPAK) for this research work is also highly appreciated.

#### **5. References**

