Flood modelling and its impacts on groundwater vulnerability in sub-Himalayan region of Pakistan: integration between HEC-RAS and geophysical techniques

Abstract Hydropower projects play a pivot role in the development of a country. Constructions of reservoirs create job opportunities and provide cheap energy but at the same time cause several environmental issues. The current study utilizes geoelectric and flood modelling data to develop a relationship between flood scenarios and their effect on river flows in association with the vulnerability of groundwater due to hydroelectric projects on Neelum River, Pakistan. The resistivity data delineated local aquifer systems that comprised both confined and unconfined aquifers ranging from 05 to 48 m depth, having poor to weak protective capacity with good groundwater development potential. The flood zonation models indicate a decline in the flow rate of the Neelum River from 317 to 39 m3/s with a drop in stage and flow velocity that contributes to a high risk of leachate penetration in the poorly protected shallow aquifer. The aquifer systems that mostly lie near the banks of the river face a serious threat of contamination due to low river flow. The flood modelling revealed that in case of dam burst, maximum probable flood will affect the land cover of 30,43,250 m2 and 33,64,433 m2 in Muzaffarabad and Patikka areas, respectively, affecting major population.


Introduction
Pakistan is among the most water-stressed countries in the world, and the growing groundwater quality problems are of grave concern for the increasing population (Nisar et al. 2021).During the summer time, when the stage is low, the dumped effluents get settled near the river bank thus contaminating the subsurface aquifer systems by infiltration.In case of floods or unfortunate dam bursts, the contaminated water finds its way outside the riverbank into the open wells contaminating the groundwater (Nisar et al. 2018;Khair et al. 2019).
Neelum Jehlum Hydro Electric Power (NJHEP) project is located in Pakistan on Neelum River and generates 969 MW of energy for domestic and industrial consumption.This River serves as the major water resource for Neelum and Muzaffarabad districts as it provides water to a population of about 0.84 million, drains 61,776.3ha of land, and is a source of water for hydropower production (Salik 2015).The storage of water and the diversions from the power projects have decreased the River stage from 10 to 2 m in the District Muzaffarabad.
Due to urbanization and growing population, deposition of domestic waste has increased in the River (Siddique et al. 2017).The heavy rainfall during monsoon season and melting of the snow and glaciers in the upper catchment area cause flood in the river (Khan and Rahman 2005).Although the construction of dams will reduce the flood risks, the low river stages will allow localized contamination associated with industrialization and urbanization to dominate the river water resulting in contamination of the aquifer system close to the river bed.
Geophysical techniques make use of the basic physics principles and integrate them with the geology to get the true picture of subsurface abnormalities (Araffa 2013).Among all the available approaches vertical electrical sounding (VES) is rated as the best practice for the vulnerability assessment of such floods/low flows in terms of surrounding lithologies (Nisar et al. 2018).The variations in acquired resistivity also define the contaminants and proper correlation provides pathways for the contaminants (Tellam and Barker 2006;Mo� scicki et al. 2014;Sentenac et al. 2017;Medici and West 2022).
Among several flood modelling software, HEC-RAS is a potential tool for the prediction of flood extent (Kute et al. 2014;Khattak et al. 2016;Thaileng et al. 2016).This software enables the user to carry out one-dimensional steady flow and twodimensional unsteady flow calculations with several input parameters that makes it a first choice among its competitors.
Due to complex geology and active fault system, NJHEP faces continuous threat.The dam is constructed on the hanging wall of the Main Boundary Thrust (MBT).MBT lies between Murree Formation and Panjal Formation.Panjal Formation (consists of marble-quartzitic, agglomeratic slates, dolomitic limestone, sandstone to volcanic rocks) is acting as a hanging wall and is highly sheared due to two active fault systems i.e.MBT and Panjal Thrust (PT).The dam is also located at the eastern limb of the Hazara Kashmir Syntaxes (HKS).The Hazara Kashmir seismic zone is also a very active region.The dam leakage was reported by Khan et al. (2018) and temporarily treated by grouting.MBT is the main cause to produce the severe seismic hazard (Parajuli et al. 2021, Bhusal et al. 2022) for subduction region formed by convergence of Indian and Eurasian plates and in recent revised seismic code of pakistan PGA for maximum credible earthquake is 1.08g (BCP, 2021).The seismic activity along MBT or PT can cause dam failure.The other problem arises when the low river flow mixed with local sewage and industrial effluents cause river water to contaminate.Later on, such water infiltrates and contaminates the subsurface aquifers.This research focuses on the VES and flood-based datasets collected across Patikka and Muzaffarabad regions for assessing the vulnerability of groundwater in flood and drought conditions.The VES data will identify the subsurface aquifer system and calculate its protective capacity by using the Dar-Zarouk parameters.Similar studies have been conducted for the determination of potential and vulnerability of aquifer by Niaz et al. (2016Niaz et al. ( , 2018Niaz et al. ( , 2022) ) in Bhimber, Rawalakot and Bagh districts of Azad Kashmir (Table 1).The HEC-RAS provides information about different flood scenarios to predict the behaviour and extent of the river during floods by generating the flood zonation maps.Thaileng et al. (2016) applied HEC-RAS technique to study the flood in a River Reach in Cambodia.The current study is aimed to demarcate the type of aquifers in the area and their protective capacity against the leachate infiltration due to decrease in the Neelum River flow by the dam construction in the upstream.The study is also aimed to analyze the flood scenarios in different seasons and in case of dam failure due to seismic activity along active faults in the area.The following objectives are followed to perform the current research work: 1. Demarcation of the aquifer system along the Neelum River using geoelectrical survey.2. Vulnerability assessment of the aquifers especially in low stages.3. To analyze the flood scenarios in case of dam failure.

Study area
The study area is a part of HKS.The MBT, PT, Kashmir Boundary Thrust (KBT) and Jhelum fault (JF) are the major constituent of HKS (Figure 1) (Khan et al. 2018).The surface land cover of the study area comprises grass lands followed by urban land that mostly resides near the river.Forests cover a considerable area apart of the mountains as well (Figure 2).
The KBT is an active seismic body (Mw 7.6) that runs across the Neelum River in Muzaffarabad city.The site has major seepage issues along with the dam site because  16,894.53m 2 /day (Mughal et al. 2018).There is a constant threat of tectonic activity which may lead to a dam failure resulting in serious flood in the Neelum River (Report of Neelum Jhelum Hydropower Company 2019).

Hydrogeology of the area
The hydrogeology of the study area is highly variable due to structural complexity.Murree Formation serves as possible aquifer formation (Niaz et al. 2016).The Murree Formation consists of cyclic deposition of clays and sandstone with subordinate siltstone and mudstone (Shah 2009).The confined aquifers in the area are formed due to the cyclic deposition of clays and sandstone, and unconfined aquifers are developed due to the deposition of river sediments as a top layer.The area is highly rough, due to which availability of the groundwater is a big problem in areas away from streams.The aquifers in the region are mostly unconfined in nature and lie mostly in the vicinity of Neelum River.The recharge of the aquifers is associated with Neelum River and seasonal rainfall.These aquifers discharge as topographic springs in low lying areas and are important source of drinking water.The general groundwater flow direction is northeast to southwest.The average rainfall in the area is 1308 mm and temperature ranges between 18 � C and 39 � C in winter and summer seasons, respectively.

VES survey
The methodology is based on four main stages (Figure 3).A total of 16 stations were used to carry out the VES survey (Figure 1).Data were acquired by adopting Schlumberger array technique in ABEM SAS 4000 Terrameter.The acquired data were processed in IPI2WIN software to develop the relationship between the apparent resistivity data and electrodes spacing in the form of a VES curve (Nisar et al. 2018).To interpret the data, partial curve matching technique was used (Figure 4).The lithological column was established on the basis of true resistivity values (Table 1).Dar-Zarouk parameters (i.e.longitudinal conductance and coefficient of anisotropy) were computed from acquired resistivity data.Longitudinal conductance (Eqs.( 1) and ( 2)) is the conductance of current along the direction of the bedding plane through the column of 1 m.It is represented by 'S' and its units are 'mhos' (Niaz et al. 2016).
where 'h' is thickness, 'q' is the true resistive and 'n' is the number of the layers in the section.
The anisotropy coefficient (Eq.( 3)) (Nisar et al. 2018) is defined as the square root of the resistivity ratio measured in two standard directions.
where 'k' is the anisotropy and 'qt' and 'ql' are the transverse and longitudinal resistivity, respectively.The higher values of anisotropy reveal near-surface materials and inhomogeneity within the subsurface (such as topsoil and partially weathered or weathered layers) which may prevent groundwater development.Low values that raised as a result of structural features such as faults, joints and fractures are important in the groundwater development of the area (Olasehinde and Bayewu 2011;Bayewu et al. 2014).

Flood zonation data
The high-resolution digital elevation model (DEM) of 12.5 m resolution was downloaded from the Alaska satellite facility.HEC-RAS software was used for data processing and flood model development.Data preparation was done using the Arc Catalog, Arc Map and the Spatial Analyst tools of ArcGis 10.4.1.All the maps created in this study were projected to UTM (Zone 43 N) coordinate system with datum of WGS-1984.The resulting raster maps of the elevation files created in ILWIS were then converted into a mosaic raster in order to obtain a single DEM raster map of the study area.
In hydraulics engineering, mainly physical and empirical parameters are practiced.Due to complex problems of hydraulic engineering, the empirical parameters are specified by uncertain values with Manning's coefficient.Determination of these coefficients presents an influential and creative task of open channel flows.Several methods have been proposed (e.g.Ramesh et al. 2000;Ding et al. 2004;Abdul Hameed 2010;Ali et al. 2010;Parhi 2013;Shamkhi and Attab 2018) for the determination of manning coefficients.In this study, the HEC-RAS software estimated the normal depth by using the manning formula (Eq.( 4)) (Chow 1959;Kute et al. 2014).
where Q is the flow rate, A is the cross-sectional area of flow, R is hydraulic radius, S is the slope of the channel, n is the surface roughness and K is the constant that depends upon these units.

Geoelectric sections
Geoelectric subsurface lithologies for the Muzaffarabad (Figure 5(a)) and Patikka profile (Figure 5(b)) reveal four major layers.The first layer comprised the sandy soil with gravel and limestone beds at some stations.The second layer is composed of gravel, sandstone with mudstone and clay.The third layer is mostly clay and sandy clay.The fourth layer represents the bedrock that is sandstone throughout the study area.The shallow layers represented the unconfined aquifer system with coarser lithologies as the topmost layer recharging the aquifer.The fourth layer is identified as the confined aquifer that has massive sandstone bed overlain by clay beds.Mughal et al. (2018) also mapped the Murree Formation and reported similar lithological units supporting water accumulation.

Dar-Zarouk parameters and shallow aquifer thickness
Dar-Zarouk parameters assess the protective capacity of the aquifer for surficial contamination; additionally, these parameters provide an idea regarding the hydraulic conductivity of the aquifer and groundwater development (Khan et al. 2021).The classification of Oladapo and Akintorinwa (2007) was used for the protective rating of the aquifer based on longitudinal unit conductance (Tables 2 and 3).The unconfined aquifer thickness maps were generated to have an idea regarding variations in storage of unconfined aquifer.These maps will also serve as validating tools for assessing aquifer protective capacity.GEOMATICS, NATURAL HAZARDS AND RISK

Muzaffarabad area
The value of longitudinal conductance in the Muzaffarabad area ranged between 0.0021 and 0.014 mhos (Figure 6

Patikka area
Across the Pattika profile, the value of the longitudinal conductance ranges between 0.0029 and 0.04 mhos (Figure 7(a)).The western part of the Patikka profile includes  conductance maps (Figure 7(a)) also indicate higher protective capacity towards the north side indicating maximum bearing capacity aquifer system that has good protective capacity.This area is less vulnerable to contamination.

GIS-based flood zonation
The Neelum River profile for the past 4 years from January 2015 to January 2019 revealed that the flow rate has decreased from 317 to 39 m 3 /s (Figure 8(a, b)).The decrease in the flow is due to the maintenance of the reservoir's water level.The mean stage in the Neelum River before the construction of the NJHEP and Kishanganga hydropower project (KHP) was 4 m (Figure 9) at Muzaffarabad station (RS 1040) with a flow velocity of 1.3-1.8m/s (Figure 10 (a, b)).The mean water level     at Muzaffarabad station (RS 1040) has dropped to 2 m with a flow velocity of 0.4-0.8m 3 /s after the construction of projects (Figure 10(c, d)).The remarkable decrease in area covered by the Neelum River suggests flood is a threat to the population living near the banks.
The HEC-RAS simulation was run on the flow rate of the Neelum River during the flood of 2015 and 2016 to construct the flood model.The model shows the river stage at the Muzaffarabad station (RS 1040) rises up to 03 m high relative to the mean stage before the project.The flow velocity at the station (RS 1040) was 03 m 3 /s during the flood.The flood covered an area of 25,28,931.8m 2 , which will affect the land covered area of 17,25,296.8m 2 (Figure 11(a, b)).The dam reservoir of the NJHEP lies on the MBT which is one of the most active faults of Himalayas.The major seismic activity may cause a catastrophic event in the form of a dam burst (Niaz et al. 2018).The model for maximum possible flood in case of dam failure shows that the stage rises up to 10 m at the Muzaffarabad station (RS 1040).This flood water indicates a maximum velocity up to 11 m/s (Figure 10(d)).A gradual change in the flow path of the river is observed despite the construction of hydropower projects due to the presence of competent lithologies (Mughal et al. 2018) surrounding the river path (Figure 11(a)).The inundation map (Figure 11(b)) indicates that during the maximum flood or dam burst, areas of Dahani, Chahla Bandi and Domail might be severely affected.However, Makri and Bala Noor Shah will not be severely affected attributing to their topographically high location.Domail, Chahla Banda, Makri and Bala Noor shah are among the populous areas near the river bank, and flood will effect maximum population in these areas (Figure 2).The Neelum River mean stage, before the construction of the hydroelectric project, was 2 m high (Figure 12) at Patikka station (RS 21050) where the water flows with a velocity of 4 m/s (Figure 13).The mean stage at Patikka station (RS 21050) drops to 1 m and flow velocity decreases to 2 m/s after the project (Figure 13(a)).The Neelum River covered an average area of 6,76,198 m 2 before the dam project that is reduced to 3,75,283 m 2 after the NJHEP and KHP projects.This indicates that floodwater was 2 m higher relative to the average mean stage before the projects (Figure 13(a, b)).The flow velocity of the flood water will be 5 m/s at the station (RS 21050).The maximum flooding over the Patikka area covered an area of 21,49,921.8m 2 .The flood affects the land cover area of 14,73,723.8m 2 that includes vegetation, urban and barren land (Figure 2).A river velocity variation is shown in Figure 14 for Neelum River stages before, after and maximum construction of hydroelectric project.
The model for maximum possible flood in case of dam failure in Patikka shows that the water level rises up to 13 m (Figure 13(d)).This flood water has a maximum velocity up to 10 m/s.The flood zonation map (Figure 15(a)) shows a slight change in river path after the construction of hydropower projects.

Discussion
The aquifers identified in the Muzaffarabad area are having a thickness between 3.5 and 46 m, while in Patikka area the identified aquifers had thickness between 2.5 and 21 m.These aquifers represent the extension of Murree formation in the subsurface (Baig and Munir 2007).The distribution of lithologies indicates fluvial deposition with episodes of high and low stream velocities.The shallow coarse lithologies serve as carriers for the rain/surficial water to the underlying aquifer (Nisar et al. 2023).The shallow and comparatively less thick aquifer system in Pattika area indicates high stream velocity that has resulted in the deposition of coarser lithologies at shallow level as compared to finer lithologies.The surface topography plays a vital role in sediment deposition, with steep slopes and tectonic instability rock masses keep on sliding into the river (Baig and Munir 2007).That is the reason why coarser sediments in the vicinity of the Muzaffarabad city are presented despite the low stream velocity.
The coefficient of the anisotropy across the Muzaffarabad profile ranges from 1.013 to 4.298 (Figure 6(b)).The relatively high value of 4.3 observed in the central part of Muzaffarabad indicates the near-surface inhomogeneity, suggesting the presence of structural features (Nisar et al. 2021) that support the groundwater development (Figure 6(b)).In the north and south, relatively low value of 1.06 is presented that indicates homogeneity such as topsoil and weathered surfaces that do not support groundwater development conditions (Niaz et al. 2016(Niaz et al. , 2018)).The thickness map for Muzaffarabad area indicates an increase in aquifer thicknesses towards the north and east of the area (Figure 6(c)).Higher thicknesses indicate higher water bearing capacity and vice versa.The longitudinal conductance map (Figure 7(a)) also indicates higher protective capacity towards the north side and indicates maximum bearing capacity aquifer system that has good protective capacity (Niaz et al. 2018).These areas are less vulnerable to contamination and leaching due to the presence of good overburden protective capacity.
The HEC-RAS flood modelling indicates alarming results in case of a dam burst.The majority of population in Muzaffarabad and Patikka resides near the banks of the Neelum River.In case of a dam burst, maximum populated areas might be under water    the areas of Pattika, Ghori, Parsaoncha will be under water.The inundation results indicate that during the flood/dam burst situation, these areas will be under more depth of water due to their location (Figure 15(b)).The population of the Ghori and Patikka area is mostly settled near the Neelum River bank.The maximum possible flood might cover the majority of the populated areas of Ghori and Patikka (Figures 2 and 15).Shallow aquifer system identified in Muzaffarabad and Pattika overlain by coarser sediments are a threat for surficial contamination.These shallow conductive systems are also identified by Dar-Zarouk parameters.In case of high flood, these conductive zones will be under flood water with conductive lithologies (as suggested by flood models); the flood water will quickly move into the shallow aquifer system thus contaminating it.

Conclusions
The geoelectrical results reveal the presence of subsurface confined and unconfined aquifers at variable depths.These shallow aquifer systems lie near the bank of the Neelum River which indicates that they are under direct impact of river flow variations.The low values of longitudinal conductance suggest that the majority of the aquifers present have the poor protective capacity and are highly vulnerable to contamination.The anisotropy values are relatively high suggesting the presence of the structural features that support the groundwater development conditions.These results are also supported by previous researchers (i.e.Mogaji et al. 2007;Niaz et al. 2016).
The Neelum River flow rate has dropped from 317 to 39 m 3 /s after the NJHEP and KHP.The high decline in the flow rate aggravates the risk of leachate penetration from sewage water and wastes along the river, thus threatening the shallow aquifer system.Furthermore, the hydropower project lies in earthquake-prone area and therefore faces constant threat of breakout.The flood modelling carried out with projections of dam break reveals that the maximum flood rises up to 11-13 m with the flow rates exceeding 10 m 3 /s in Muzaffarabad and Patikka areas.The Neelum River covers an area of about 38,46,885.2m 2 , which is 30,43,250.2m 2 land cover of Muzaffarabad that would be under devastating flood.The same river covers about 40,40,631.1 m 2 , which is the 33,64,433 m 2 land cover area of Ghori and Patikka that would be under the havoc flood.
Figure 1.(a) Geological map of the study area, locations of VES points, Muzaffarabad and Patikka profiles, are also highlighted.(b) Tectonic map of the Hazara Kashmir Syntaxes (HKS), rectangle highlights the study area (after Javed et al. 2021).

Figure 2 .
Figure 2. (a) Land cover distribution map.(b) Bar graph of land distribution classes of the study area.

Figure 3 .
Figure 3. Flow chart of the research methodology.

Figure 4 .
Figure 4. Resistivity curves showing subsurface resistivity responses across the study area.

Figure 5 .
Figure 5. Lithologies interpreted based on the VES points across the (a) Muzaffarabad and (b) Patikka profiles.
(a)).The longitudinal conductance shows a low value (0.007 mhos) in the southern part near Domail and Bala Noor Shah Towns (Figures1 and 6(a)).The value decreases (0.002 mhos) in the Chala Banda, Makri area and Dhani area (0.014 mhos).The low values indicate that the protective capacity is poor in the southern part of the area due to the presence of coarse lithologies(Nisar et al. 2018).These areas are highly vulnerable to the infiltration of leachate or contamination due to the presence of less thick overburden capacity.The shallow aquifers are directly recharged by the river water, and overburden protective capacity thickness plays a vital role in surface and groundwater interaction during recharging.The quick interaction between surface and groundwater occurs in areas of less thick overburden protective capacity and are prone to contamination.Shafi et al. (2018) also confirmed the bacterial contamination in tap water of the village Killah, Muzaffarabad.The contamination is possibly due to the reduction of flow of Neelum River and as a result addition of sewerage lines directly into the river.Ali et al. (2019) delineated that the heavy metals like lead and chromium are also exceeded from the WHO guidelines in the District Muzaffarabad.Zahoor et al. (2022) delineated high concentration of electrical conductivity, turbidity, calcium, magnesium, and a similar level of coliform in the District Muzaffarabad.

Figure 8 .
Figure 8.(a) Flow rate of the Neelum River from January 2015 to January 2019.(b) Location map of river stations used in the study.

Figure 9 .
Figure 9. Elevation diagram of Neelum River at Muzaffarabad station (RS 1040) (a) before the dam construction, (b) after the dam construction, (c) during the flood conditions and (d) during maximum probable flood condition.

Figure 10 .
Figure 10.Velocity diagram of Neelum River at Muzaffarabad station (RS 1040) (a) before the dam construction, (b) after the dam construction, (c) during the flood conditions and (d) during maximum probable flood condition.

Figure 11 .
Figure 11.(a) Flood zonation map of Muzaffarabad and (b) inundation map of Muzaffarabad area.

Figure 12 .
Figure 12.Elevation diagram of the Neelum River at Patikka station (RS 21050) (a) before the dam construction, (b) after the dam construction, (c) during the flood conditions and (d) during maximum probable flood condition.

Figure 13 .
Figure 13.Velocity diagram of the Neelum River at Patikka station (RS 21050) (a) before the dam construction, (b) after the dam construction, (c) during the flood conditions and (d) during maximum probable flood condition.

Figure 14 .
Figure 14.Velocity profiles of the Neelum River (a) before the dam construction, (b) after the dam construction and (c) during the flood conditions.

Figure 15 .
Figure 15.(a) Flood zonation map of Patikka area and (b) inundation map of Patikka area.
and 15).The flow rate of the Neelum River varies between 229.429 and 340.341 m 3 /day in years 2002-2009 as shown in Figure 16 (Munir 2013).The flow rate decreases from 2002 to 2004 and then elevated in years 2005 and 2006 (Dad et al. 2021).The decline in flow rate is observed after 2006 associated with construction of different power projects on Neelum River.It is worth mentioning that during the maximum flood,

Figure 16 .
Figure 16.Average runoff fluctuation of Neelum River between years 2002 and 2009.

Table 1 .
List of researchers that carried studies related to the techniques used in current research.

Table 3 .
Longitudinal unit conductance against protective capacity rating.