Groundwater modeling of Musi basin Hyderabad, India: a case study

In general, groundwater flow and transport models are being applied to investigate a wide variety of hydrogeological conditions besides to calculate the rate and direction of movement of groundwater through aquifers and confining units in the subsurface. Transport models estimate the concentration of a chemical in groundwater which requires the development of a calibrated groundwater flow model or, at a minimum, an accurate determination of the velocity and direction of groundwater flow that is based on field data. All the available hydrogeological, geophysical and water quality data in Musi basin, Hyderabad, India, were fed as input to the model to obtain the groundwater flow velocities and the interaction of surface water and groundwater and thereby seepage loss was estimated. This in turn paved the way to calculate the capacity of the storage treatment plants (STP) to be established at the inlets of six major lakes of the basin. The total dissolved solid was given as the pollutant load in the mass transport model, and through model simulation, its migration at present and futuristic scenarios was brought out by groundwater flow and mass transport modeling. The average groundwater velocity estimated through the flow model was 0.26 m/day. The capacities of STP of various lakes in the study area were estimated based on the lake seepage and evaporation loss. Based on the groundwater velocity and TDS as pollutant load in the lakes, the likely contamination from lakes at present and for the next 20 years was predicted.


Introduction
The use of groundwater models is prevalent in the field of environmental hydrogeology. Also groundwater flow and transport models are being applied to investigate a wide variety of hydrogeologic conditions. Models are used to calculate the rate and direction of movement of groundwater through aquifers and confining units in the subsurface (Toth 1963;Witherspoon 1966, 1968;Bear 1972Bear , 1979Bear and Verruijt 1979;Anderson and Wang 1982;Kinzelbach 1986;McDonald and Harbaugh 1988;Anderson and Woessner 1992). Transport models estimate the concentration of a chemical in groundwater beginning at its point of introduction to the environment to locations down gradient of the source (Scheidegger 1961;Fried 1975;Javandel et al. 1984, Grove andStollenwork 1984;Konokow 1976Konokow , 1977Konokow , 1978Konikow and Bredehoeft 1978;Konikow and Grove 1977). Transport models require the development of a calibrated groundwater flow model or, at a minimum, an accurate determination of the velocity and direction of groundwater flow that is based on field data. Essentially, mathematical modeling of a system implies obtaining solutions to one or more partial differential equations describing groundwater regime. The partial differential equation describing three-dimensional groundwater flows may be written in a homogeneous aquifer as where K xx , K yy and K zz are the hydraulic conductivity (LT −1 ) along x, y, z directions, h is the piezometric head (L), S s is the specific storage of the medium, W is the volumetric flux per unit area and represents sources and/or sinks, and t is the time.
Usually, it is difficult to find the exact solution of Eq. (1) and one has to resort to numerical techniques for (1) ∕ x K xx h∕ x + ∕ y K yy h∕ y + ∕ z K zz h∕ z = S s h∕ t ± W 1 3 14 Page 2 of 22 obtaining approximate solutions. In the present study, finite difference method was used to solve the above equation. Herein, first a continuous system is discretized (both in space and in time) into a number of node points in a grid pattern. The partial differential equation is then replaced by a set of simultaneous algebraic equations valid at different nodes. Thereafter, using standard methods of matrix inversion, these equations are solved. Models are conceptual descriptions or approximations that describe physical systems using mathematical equations-they are not exact descriptions of physical systems or processes. The applicability or usefulness of a model depends on how closely the mathematical equations approximate the physical system being modeled. A mathematical model, in general, portrays all the characteristic features of a full-scale system and can be studied under controlled conditions. In order to identify these characteristics, it is imperative, at the outset, to define the system frame work, physical/chemical processes to be modeled, and the objective of modeling. It is only then that an appropriate simulation model can be constructed to represent all the vital characteristics. Thus, a groundwater flow model essentially simulates the groundwater movement and its spatiotemporal variations. The main objective of the flow model is to study the variation of groundwater heads to evolve optimal exploitation schemes. A chemical quality model of an aquifer needs to simulate not only the contaminant movement due to fluid velocity but also the migration arising from the variations in concentration of contaminant and density gradients. The process of aquifer modeling, in general, comprises the following activities: • Identification of parameters characterizing the physical framework of the aquifer and stresses acting on it. • Estimation of the relevant hydrogeological parameters at as many points as possible, particularly those at the boundaries • Interpolation/extrapolation of these parameters to characterize the entire area of study • Integration of the entire data to conceptualize and resurrect the natural system. • An appropriate mathematical description of the conceptual model giving spatiotemporal relationship between those parameters constrained by the relevant fundamental laws of hydraulics, mass transport, etc. • A solution of mathematical equations describing the groundwater regime in terms of observables such as groundwater levels or concentrations of pollutants.
• A sensitivity analysis of the model to identify those parameters which need to be estimated more accurately, and also to decipher the error bounds. • Refinement of the model to progressively bring in plausibility and compatibility between field estimates of various geohydrological parameters through the process of model calibration and validation. • Prognosis of the aquifer response to evolve efficient management options for optimal utilization of groundwater resources.
All the parameters in Eq. (1) should be known at maximum number of points. The validity and applicability of an aquifer model directly depend upon the adequacy and reliability of the data. A comprehensive understanding of the study needs a few relevant topics including geology of the area of study, and therefore it is organized as "Location and geology of study area," "Physical framework and stresses on the groundwater regime," "Methodology of study" under which are "Lake water budget," "Groundwater flow modeling," "Capacity of storage treatment plant (STP)," "Results and discussion" followed by "Conclusions."

Location and geology of study area
The study area is greater Hyderabad, a Metropolitan city in India, covers an area of 40 km 2 and lies between the latitudes 17°23′N to 17°29′30″N and longitudes 78°31′E to 78°35′30″E. It falls in the NE Musi basin which consists of chain of surface water bodies (lakes), namely RK Puramcheruvu, Nadimicheruvu, Bandacheruvu, Patelcheruvu, Peddacheruvu and Nallacheruvu (Fig. 1). The twin cities as is known (Hyderabad-Secunderabad) reflect an undulating topography interspersed with many hillocks and knolls. The drainage is dendritic, characterized by irregular branching of tributary streams in many directions. The drainage system slopes toward the south to join river Musi. The twin city area is underlain by coarse porphyritic granite containing large plagioclase phenocryst and abundant biotite. The origin of granite is considered to be either late or post-tectonic. The Hyderabad granites, which form part of the Dharwarian craton, are referred to as the basement complex or unclassified gneissic complex (Krishnan 1960;Sitaramayya 1968). The study of hydrogeology of the watershed is essential, and observation of the stream network, topography and morphological features have significant bearing on the surface water-groundwater interaction in the area. A major portion of the study area and its environs consists of granites wherein the occurrence, distribution, storage and movement of groundwater are complex phenomena. Further details on geology, physiography, drainage, climate, rainfall, hydrogeology, etc., are area available in the literature (Sundararajan et al. 2012(Sundararajan et al. , 2015.

Physical framework and stresses on the groundwater regime
The physical framework of an aquifer model is defined by the hydrogeological parameters, namely the transmissivity (T) and storativity (S), boundary conditions and the hydraulic heads. The hydrogeological parameters T and S are estimated through long duration pumping tests. The boundary conditions include constant head, no flow or constant flow which may exist at the geometrical boundaries of the aquifer under study. The hydraulic heads (water level) are the most important observable parameter which characterize the response of an aquifer system. Since the model response is calibrated with the observed hydraulic heads, care should be taken to collect maximum data on water level from observation wells and these data have to be reduced with reference to a common datum, namely mean sea level (msl). The stresses on groundwater regime include the input to and output from the aquifer system. Recharge due to precipitation, infiltration from surface water bodies like rivers, lakes, tanks, canals, irrigation return seepage and lateral groundwater inflow comprise the input to an aquifer system. All the above data required for groundwater flow modeling are also required for mass transport modeling. In the mass transport model, the pollutant is a major input termed as stresses on the quality of groundwater. Groundwater draft, evapotranspiration, lateral groundwater outflow, effluence to surface water bodies comprise the output from an aquifer.

Lake water budget
The main objective of the entire study was to estimate the capacity of the storage treatment plants (STP) to be established at the inlets of different lakes, and the contribution of lake seepage to groundwater is the key parameter in the lake water budget. Hydrological studies of lake-watershed systems are typically based on quantifying the water balance among precipitation, evaporation, surface runoff, lake discharge and groundwater inflow to and outflow from the lake. Since all components affect the volume of water in lake, the relationship can be expressed as: where S Lake -change in volume of water in the lake, Prainfall directly on the lake, SW in -surface water inflow to lake (e.g., rivers, streams), GW in -groundwater discharge to the lake, SRO-surface runoff directly to the lake, E-evaporation directly from the surface of the lake, SW out -surface discharge from the lake (river, stream), GW out -groundwater discharge from the lake, and Others-other factors affecting lake storage (e.g., irrigation, withdrawal of water for drinking and industry).
Contributions to the volume of water within a lake from the hydrological components other than groundwater are measurable, or in the case of evaporation, calculated by the standard formulae given by Healy and Cook (2002). However, the quantities of groundwater flow into or from a lake cannot be measured accurately due to the geological complexity of subsurface and detailed field methodology required to adequately characterize the hydro-stratigraphy. Therefore, for convenience, typical studies assume that groundwater has little significance in the water budget or have ignored the groundwater component completely. Thus, under these assumptions, the balance equation, which is typically solved, is given as: The hydrological balance given by Eq.
(2) cannot be used to uniquely assess the role of groundwater inflow and outflow in the water budget of a lake, or its associated contaminant or nutrient input. All values other than groundwater inflow and outflow in Eq. (2) can be measured, and the substitution of these known values into Eq. (2) will produce one equation with two unknowns (GW in and GW out ) which is typically included in the balance equation as groundwater flux to the lake (GW Lake ) given as: Thus, there is not a unique solution to estimating groundwater inflow to or outflow from a lake. At best, only the net groundwater flux to the lake (GW Lake ) can be estimated, or one of the two terms in Eq. (4) is set to zero to allow the other to be estimated (Cole and Fisher 1979;Cook et al. 1977). The contribution of the lake to the groundwater was estimated through the flow model by simulating the lake water and groundwater interaction. (2)

Groundwater flow modeling
The initial stage in developing groundwater flow model is to define the region of interest and establish the boundary conditions for flow. The geophysical investigations carried out in the watershed provided the necessary insight into the aquifer geometry for development of a groundwater flow model. The following information was used for model conceptualization: • No flow occurs across the watershed boundaries and these boundaries coincide approximately with groundwater divides. • Groundwater recharge takes place from the top layer of the watershed. • Continuous groundwater pumping is prevalent in the watershed due to urbanization and land cover. • The groundwater withdrawal was estimated based on the well inventory, average running hour of pumps and residential distribution in the area. • As the watershed is a closed one with streamlet, some outflow may take place. • Seepage from the lakes was an additional input of recharge to the watershed.
Maximum information was abstracted from hydrogeologic, geophysical and water quality databases generated (Sundararajan et al. 2012(Sundararajan et al. , 2015. The entire area was divided into 76 × 106 square grid. Each grid is a finite difference block with a node located at its center covering an area of 100 m × 100 m (Fig. 2). The water level data measured at 73 observation wells during June 2003 were used for steady state model calibration. Visual MODFLOW (Guiger and Franz 1996) software package was used for simulating the groundwater flow. MODFLOW river package was used to incorporate the surface water boundary conditions into the groundwater flow model. The lake and aquifer interaction was simulated by assigning the lake water levels and lakebed elevation data collected from the lake bathymetry as shown in Table 1. The permeability distribution was based on the layer parameters derived from the resistivity investigations and transmissivity values derived from the pumping test data (Rao et al. 2001(Rao et al. , 2003. The initial permeability varied from 2 to 4 m/day, and the specific storage was assigned as 0.005. The simulated vertical section varied from 35 to 50 m depth based on the geoelectric cross sections as the two-layer weathered and fractured aquifer system. An initial recharge of 10% of the total rainfall was fed through the recharge package to the aquifer system. The natural recharge due to urbanization was very low varying from 25 to 50 mm/year, and the aquifer was assumed to be replenished from the lake water and groundwater interaction. The outflow toward Musi river at the southern end of the study area was simulated as a constant head of 462 m (amsl).
The river and constant head boundary condition is shown in Fig. 3. The groundwater withdrawal from the basin was simulated appropriately through well package considering the urbanization, land use, etc., varying from 10 to 50 m 3 / day (Fig. 4). Steady state calibration was achieved for the best fit between observed and computed water levels through progressive manipulation of the permeability values (Table 2), and the permeability distribution (Freeze and Witherspoon 1967) is shown in Fig. 5. The calibrated blockwise recharge worked out is shown in Fig. 6. The match between the computed and observed water levels with error bounds is shown in Fig. 7, and the computed water level contours of June 2003 is shown in Fig. 8. The groundwater balance under steady state conditions is shown in Table 3. The computed versus observed head at 24 observation wells was found matching closely. The computed water level contours are following the observed water levels reported earlier for June 2003. Assessment of interaction between lake water and groundwater regime of six lakes has been computed through zone budget by assigning seven zones in the model, of which six zones represent six lakes and zone 1 represents entire watershed 175 (Fig. 9). The groundwater balance in the watershed has been computed through zone budget in zone 1 (Table 3). The lake contribution to groundwater estimated through zone budgeting of six lakes is presented in Table 4. The results of the zone budget indicate that the Peddacheruvu and Nallacheruvu are gaining groundwater. The net contribution of the lakes to groundwater was estimated as 3500 m 3 /day. The lake contribution to the groundwater is the seepage loss, an important parameter for calculating the capacity of STP to be established at the inlets of different lakes. An average rainfall of about 700 mm is being received directly on the lake surface. Annual evaporation of 2400 mm occurs from the lakes in Hyderabad. The net loss from the water surface due to evaporation works out to 1700 mm/year. Water budget has been computed with an average evaporation loss of 4.5 mm/day from lake water surface at full tank level of Patelcheruvu, Peddacheruvu and Nallacheruvu having an area of 11 ha, 17 ha and 17 ha, respectively. The proposed capacity of STP for different lakes given in Table 5 is briefly discussed hereunder.

Capacity of storage treatment plant (STP)
The computed seepage loss from RK Puramcheruvu is 0.6 MLD. An STP of 0.6 MLD for tertiary treatment has been functioning on Nadimicheruvu. The lake surface area is about 10 ha, and evaporation loss at the rate of 4.5 mm/day from the lake water surface works out to be 0.45 MLD. As the lake is situated in recharge area of the watershed, it contributes 1 MLD as seepage to the groundwater. Under lake restoration, it is envisaged to maintain FTL at Nadimicheruvu, and then it would need a supply of 1.5 176 MLD of treated wastewater. On a conservative basis with allowance for some outflow from the lake, an STP of 2 MLD capacity may be established. Thus, it is recommended to enhance the present capacity of the STP to 2 MLD at Nadimicheruvu. Further, it satisfies the growing urbanization and the disposal of sewerage volume in near future. The computed seepage losses from Bandacheruvu are 0.8 MLD. Based on the experience of Nadimicheruvu and computed seepage losses from groundwater flow model, it is recommended that STPs with a capacity of 2 MLD each may be constructed at RK Puramcheruvu and Bandacheruvu (Table 5). In the case of Patelcheruvu, the groundwater and lake water interaction suggested that it would require about 1 MLD capacity for meeting the requirement of evaporation and seepage losses ( Table 5). The surrounding intensive urbanization plays an important role in the downstream lakes. Thus, it would be appropriate to commission an STP of 3 MLD at Patelcheruvu. The seepage loss in Peddacheruvu is about 1.5 MLD and evaporation loss is about 1 MLD (Table 5). The minimum capacity of 2.5 MLD is required for meeting the daily demands of lake water budget. But the ambience of urbanization is very high around Peddacheruvu for generation of higher volumes of domestic sewerage. Some industrial discharges adjoining the lake are letting out their wastewater into the lake. Considering all the above conditions, it is recommended for the construction of an STP with a capacity of 10 MLD for tertiary treatment.
Nallacheruvu contributes the lowest seepage losses to the groundwater regime, as it is located in discharge area. The seepage losses account only 0.25 MLD, and evaporation losses will be about 1 MLD (Table 5). Urbanization is very high surrounding the lake and generation of sewage is about 4 MLD. Since the lake is situated at the end of watershed, the diverted surplus flows from upstream lakes could also pass through it. Taking all the likely diverted and surrounding sewerage flows into consideration, an STP of 10 MLD would be necessary for maintaining the lake water level at FTL under the lake restoration program. The flow measurements carried out at the inlet/outlet of the above three lakes support sustainable influent streams of sewerage to the proposed STPs. The groundwater velocity estimated through the flow model was 0.26 m/day. This parameter was used in the mass transport model. Further, the likely contamination of groundwater through lake water interaction has been analyzed through mass transport modeling by assigning TDS concentration of sewage as input in the lakes. The Mass Transport Model Software (MT3-D, Zheng 1990) has been utilized in conjunction with visual MODFLOW for the basin. The database of 2003 with an overall load of 1000 mg/l of TDS in all the six lakes was simulated to arrive at the existing scenario of the pollution spread for the first and second layers as shown in Figs. 10 and 11. The likely contaminant migration from lakes for the next 20 years shows that the lake water seepage would impact up to 600-800 m in the downstream of lakes (Figs. 12, 13).

Results and discussion
The database generated through hydrogeological, lake bathymetry, geophysical, and geochemical studies was used to build the groundwater flow and mass transport model. The geoelectric section obtained from the resistivity investigations helped in conceptualizing the aquifer system and calibrating the groundwater flow model. The thickness of the subsurface layers as well as the field estimates of transmissivity was used in the sensitivity analysis during the model calibration. The permeability estimates through model calibration varied from 1.8 to 3.6 m/day. The calibrated recharge worked out varied from 45 mm/year in the upstream to 25 mm/year in the downstream. The average groundwater recharge and draft were estimated as 3000 m 3 / day and 5750 m 3 /day, respectively. The net contribution of lake seepage (Fig. 14) to groundwater was 3500 m 3 /day. The zone budget results presented in Table 4 indicate that in general the annual recharge to the lake is more than that of the lake contribution to groundwater. In the case of Nallacheruvu, the groundwater contributes more and hence the lake is gaining water. Further, the percentage field saturation data presented in Table 5 reveal that the moisture percentage of Nallacheruvu and Peddacheruvu at 0.9 m depth below lakebed is more comparable to that at lakebed. This confirms the results of the zone budget wherein there is a rise of water in Nallacheruvu and Peddacheruvu, indicating the lakes gain groundwater. The average groundwater velocity estimated through the flow model was 0.26 m/day. The capacities of STP of various lakes in the study area were estimated based on the lake seepage and evaporation loss (Fig. 15). Based on the groundwater velocity and TDS as pollutant load in the lakes, the likely contamination from lakes at present and for the next 20 years was predicted.

Conclusions
Models are used to calculate the rate and direction of movement of groundwater through aquifers and confining units in the subsurface. Accordingly, the present study has revealed that the average groundwater velocity estimated through the flow model was 0.26 m/day. The capacities of STP of various lakes in the study area were estimated based on the lake seepage and evaporation loss. Based on the groundwater velocity and TDS as pollutant load in the lakes, the likely contamination from lakes at present and for the next 20 years was predicted.