Integration of Groundwater Flow Modeling and GIS

Hydrology is the science that deals with the processes governing the depletion and replenishment of water resources of the earth's land areas. The purpose of this book is to put together recent developments on hydrology and water resources engineering. First section covers surface water modeling and second section deals with groundwater modeling. The aim of this book is to focus attention on the management of surface water and groundwater resources. Meeting the challenges and the impact of climate change on water resources is also discussed in the book. Most chapters give insights into the interpretation of field information, development of models, the use of computational models based on analytical and numerical techniques, assessment of model performance and the use of these models for predictive purposes. It is written for the practicing professionals and students, mathematical modelers, hydrogeologists and water resources specialists.


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

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, www.intechopen.com Integration of Groundwater Flow Modeling and GIS 241 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 management (data storage and retrieval)  Data manipulation and analysis  Data output 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.

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

www.intechopen.com
Water Resources Management and Modeling 242 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.

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.

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

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. Spatial 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:  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, sewerage systems etc.  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 improve the model results.  GIS can handle large datasets through integration with Database Management System (DBMS) component which provides foundation for all analysis techniques.  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 the same level of accuracy each time.  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 transform spatial data according to the modeling requirements.
Application of GIS technology to hydrological modeling requires careful planning and extensive data manipulation work. In general, the following three major steps are required:  Development of spatial database  Extraction of model layers  Linkage to computer models 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,

Phase GIS functions Modeling Steps
Pre processing Data collection, Conceptualization 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).

Spatial analysis of waterlogged areas
The study area lies between longitudinal range 73-74 5E and latitudinal range 32-32 45N 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.5C and 25.4C 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.5C and 5.1C 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  Shallow watertable resulting from lenses of clay material present at shallow depth creating perched watertable conditions  High seepage rate from canal system resulting from non-lining of irrigation channels  Lack of water management practices/absence of complementary drainage system  High precipitation 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. The following thematic layers were developed through on-screen digitizing in ILWIS 3.1 GIS software: 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 km 2 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 braided riverbed constitutes about 9 km 2 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: 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.
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 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.

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, Boundary conditions: 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 (m 3 /d/m 2 ); μ 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; 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 the following form: 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.
iii. Formulation of the groundwater flow equation (partial differential equation or PDE) in integral form. iv. Element-wise integration of the integral form of the groundwater equation. v. Assembly of the algebraic matrix equations that result from the integration step into global system of linear equation of the form Where [M] denotes a matrix.
vi. Time integration. 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: 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:

Model calibration
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 Where RMSE is root mean square error (m); n is total number of measurements; H i is simulated value of groundwater table at the end of ith month (m); and H i ' 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. 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 (<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 km 2 in the base year 2005 to 74.8 km 2 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 km 2 in 2005 to 1104.8 km 2 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. 4.5-6.0m 6.0-12.0m < 1.5m 1.5-3.0m 3.0-4.5m 4.5-6.0m

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 5000m 3 /day to 8000m 3 /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.

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