Flood Inundation Modeling by Integrating HEC–RAS and Satellite Imagery: A Case Study of the Indus River Basin

Abstract

Floods are brutal, catastrophic natural hazards which affect most human beings in terms of economy and life loss, especially in the large river basins worldwide. The Indus River basin is considered as one of the world’s large river basins, comprising several major tributaries, and has experienced severe floods in its history. There is currently no proper early flood warning system for the Indus River which can help administrative authorities cope with such natural hazards. Hence, it is necessary to develop an early flood warning system by integrating a hydrodynamic model, in situ information, and satellite imagery. This study used Hydrologic Engineering Center–River Analysis System (HEC–RAS) to predict river dynamics under extreme flow events and inundation modeling. The calibration and validation of the HEC–RAS v5 model was performed for 2010 and 2015 flood events, respectively. Manning’s roughness coefficient (n) values were extracted using the land use information of the rivers and floodplains. Multiple combinations of n values were used and optimized in the simulation process for the rivers and floodplains. The Landsat 5 Thematic Mapper (TM), Landsat 8 Operational Land Imager (OLI), Moderate Resolution Imaging Spectroradiometer (MODIS) MOD09A1, and MOD09GA products were used in the analysis. The Normalized Difference Water Index (NDWI), Modified NDWI1 (MNDWI1), and MNDWI2, were applied for the delineation of water bodies, and the output of all indices were blended to produce standard flood maps for accurate assessment of the HEC–RAS-based simulated flood extent. The optimized n values for rivers and floodplains were 0.055 and 0.06, respectively, with significant satisfaction of statistical parameters, indicating good agreement between simulated and observed flood extents. The HEC–RAS v5 model integrated with satellite imagery can be further used for early flood warnings in the central part of the Indus River basin.

1. Introduction

Human existence and survival on earth are vulnerable to innumerable damages caused by natural disasters [1,2]. Among these, floods are prominent as the most severe of disasters [3,4,5], affecting several million people annually [6], and claim more human lives compared with other brutal natural catastrophes worldwide [7,8]. River and coastal floods caused a total damage of about USD 46 trillion in 2010, which is expected to further project to USD 158 trillion in 2050 [9]. Floods have wrecked the global economy [10], and ruthlessly affected billions of people across the globe [11]. Floods also cause widespread damage to agriculture, transportation, environmental ecosystems, cultural heritage sites, and an overall disturbance to global economic activities [12,13]. Two primary types of flood damage are extensive damage and indistinct damage. Extensive damage can be assessed and translated into economic terms [14], while indistinct damage such as life loss cannot be assessed as financial losses [15]. Developing countries have suffered both extensive and indistinct damages, and compared with developed countries, suffer more per unit of their GDP due to the frequent and intensive floods under the changing climate [16], as they lack material and financial resources to mitigate the impacts [17].

Among developing countries, Pakistan has suffered a total flooded area of 616,598 km² [18], with an economic loss of about USD 38.171 billion, and lost 13,262 lives during the 25 major catastrophic floods in its history [19,20]. The huge economic and infrastructure losses are unavoidable in extreme floods due to the lack of dams, while the life loss is directly related to the lack of an early flood warning system, flood preparedness and management [21]. Pakistan receives water for its survival from the mighty Indus River basin and its tributaries [22]. The situation of the Indus basin is not only due to a lack of water storage facilities on the main Indus River, but also due to the lack of dams on its tributaries which also receive extreme flows [23,24]. Researchers are predicting more severe floods owing to climate variability [25], as the changing climate proves to be the leading cause of extreme floods worldwide [26,27,28,29,30]. In climate change scenarios, the timely warning of floods using conventional means is a challenge for better management, which has given rise to the scientific awareness of flood forecasting, preparedness and mitigation strategies among people [31,32].

Floods are also among the most commonly occurring disasters that cause geomorphologic flow channel deviations [33]. Unpredicted floods instigate physical effects on the flow paths and channel beds, i.e., changes in width, variations in bed elevations, flow path locus, flow behavior, and development of islands [15,34]. Sedimentation in heavy flood water further decreases the flood-carrying capacity of streams, hence increasing the possibility of bank overflows [35], and sedimentation usually happens in bars and flood-affected areas [36,37]. When the water level in the river exceeds the river’s banks, water moves towards the floodplain, and the flows can be monitored and analyzed by computer simulations by performing flood inundation mapping. The aforementioned floods, associated damages, and management issues indicate that there is a need to properly accept the flood challenges and develop mitigation strategies [38,39,40,41], which can be helpful to minimize flood damage losses [42]. Therefore, hydrodynamic modeling and flood mapping of the rivers’ basins are required to help the administrative authorities in the case of a flood emergency [43,44,45,46]. Hydrodynamic modeling of rivers and floodplains are performed to determine the possible flood extent and water depths that could be experienced due to a hydrological response, for different areas. All processes of hydrodynamic modeling require an accurate representation of the water channel, banks, flow paths and topography of the floodplain [47].

A hydrodynamic model illustrates the processes that occur throughout a flood event. There are usually two modeling approaches for the assessment of flood inundation and flood depths, namely, a one-dimensional approach and a two-dimensional assessment approach. In one-dimensional modeling, Saint Venant equations are utilized, and in the two-dimensional modeling approach, Saint Venant equations are used in two dimensions [48]. In one-dimensional modeling, the average values of water velocity and terrain cross-section are used. The lateral flow of water, in the channel and floodplain, is not considered in one-dimension. Longitudinal flow is considered for the main channel and flood plains. A one-dimensional model is a suitable choice for assessing the direction of the flow path, but it does not consider lateral flow and the topography of cross-sections. On the other hand, in two-dimensional flood inundation mapping, the lateral and longitudinal flow of the channel and floodplains are considered. Two-dimensional models also represent the topography of the flow channel and floodplain [49]. Hydraulic modeling helps in timely warning, and evacuating, the possible affected areas downstream. Floodplain mapping requires online or offline applications [47], and offline hydraulic analysis is carried out using two-dimensional modeling in those areas where the geometry of the river flow is complex [50].

Both one-dimensional and two-dimensional analyses are carried out using different models combined with a geographical information system [51]. One-dimensional models show the flow through the river path and flooded area only in the horizontal direction. Two-dimensional models depend on incorporating the water height to acquire the average elevation of the water, and are analyzed by utilizing a suitable numerical methodology. One-dimensional modeling does not need comprehensive data, but two-dimensional modeling requires complex and detailed data of the river channel and floodplain [52]. Flow momentum, of the main channel and floodplain, is also included for steady and unsteady flow estimations in hydrodynamic modeling. The HEC–RAS v5 model has the capability of studying river dynamics under unsteady and steady state conditions [53]. HEC–RAS works in coordination with the flow data and results, as one of the productive flood inundation systems. Numerical exploration is considered one of the strategies available to carry out floodplain mapping. Therefore, in this study, critical methods were performed that helped evaluate the height of the water surface during overflow seasons. The well-known flood water levels were transformed into a topographic map, which outlined the regions of the floodplain that will be determined for flood forecasting [54].

A geographic information system (GIS) is an important tool to estimate the extent of floodwater and is also used for further analysis by creating maps to highlight damaged areas [55]. The integration of HEC-GeoRAS in ArcGIS software is a hybrid methodology for flood inundation mapping, for gathering, examining, and controlling spatial information [56]. Advanced floodplain mapping is obtained by consolidating all spatial and hydrological information [57]. The flood extent based on HEC–RAS simulations can be exported to ArcGIS software, and can be validated using satellite-based flood mapping [58,59,60] to further fine-tune the n values for the rivers and floodplain, and check the accuracy of the flood inundation modeling [61,62]. The mapping, monitoring, and variation of flood propagation can be addressed properly using a GIS [63]. The propagation of flood extent over floodplain mainly depends on the surface roughness of the rivers and floodplain [64,65]. The accurate assessment of Manning’s roughness coefficient is of prime importance in the simulation of flood extent over the floodplain [66].

Inundation modeling for flood mapping over the Indus River basin is highly complicated and a big challenge; it is considered to be one of the largest basins on the globe [67], expanding from the higher mountains of the Himalaya to the Arabian Sea [68,69,70]. Researchers from different organizations of Pakistan are trying to develop an Early Flood Warning System using advanced hydrologic and hydrodynamic modeling of the rivers’ basins. The United Nations Educational, Scientific and Cultural Organization is also trying to use different hydrological and hydrodynamic models in Pakistan and other parts of world to cope with the challenges of floods. The HEC–RAS v5 model is a popular inundation model and is being applied for inundation modeling and flood mapping worldwide [64,71,72,73,74,75,76]. The HEC–RAS model has been applied in different parts of the Indus River basin by researchers, such as in the Tori levee breach [77], Chashma–Taunsa reach [78], Kabul River reach [79], lower Indus River reach [80], and Hunza and Shyok rivers’ reaches [81,82]. However, a few other models are also being applied by researchers in different parts of the Indus basin for flood mapping and inundation modeling [83,84,85]. The central parts of the Indus River basin have still not been addressed; this area experiences severe floods, such as in 2010, which resulted in the flooding of a large area [86,87,88], depression and post-traumatic stress disorder [89,90,91,92,93], diseases and public health crisis [94,95,96,97], massive economic loss [21,98,99], environmental degradation [100,101,102], groundwater deterioration [103,104,105], and huge life loss [18,106,107]. The increased flood events and intensity, population growth, and climate change challenges have drawn the attention of scientists toward advanced spatiotemporal hydrodynamics modeling, instead of traditional flood modeling inundation. The present study focuses on advanced flood mapping by integrating an inundation model and satellite imagery in the central part of the Indus River basin using high resolution terrain information, soil type, and land use information. The calibrated and validated HEC–RAS v5 model can be further used for early flood warning in extreme future flood events.

2. Materials and Methods

2.1. Study Area

The Indus basin has a drainage area of more than 1,165,000 km², shared by China (10.7%), Afghanistan (6.7%), Pakistan (56%), and India (26.6%) [108]. The origin of the Indus River is Mansarovar Lake in the Tibetan Plateau. Before drainage into the Arabian Sea, it passes through Kashmir and Pakistan. The upper Indus basin is located in the range of 32.48° N and 67.33° E as presented in Figure 1, in the mountainous ranges of Karakoram, Hindu Kush, and Himalaya [109].

These mountainous ranges comprise 11,000 glaciers [110], making it one of the most glaciated regions in the world, with a 22,000 km² surface area of glaciers [110]. The altitude varies from north to south and ranges between 200 m and 8500 m above sea level, with an average elevation of 3750 m above sea level. The study area is one of the major sources of fresh water for Pakistan. The study was conducted in the Indus River basin, from Taunsa Barrage to Kot Mithun, Punjab, a central part of the Indus River basin. The geographic boundary of the study area lies in 70.38°–71.10° E and 28.87°–30.64° N, while the elevation ranges from 80 to 140 m. Muzaffargarh, Dera Ghazi Khan, and Rajan Pur districts are subjected to frequent floods from the Indus and Chenab rivers. The Indus River basin has experienced high magnitude floods in 1956, 1973, 1976, 1992, 1994, 2010, 2015, and 2022, in its history [111]. The floods were of medium to high magnitude, which not only caused huge economic loss to the residents, but also claimed several lives. When the flood discharge reaches up to 17,000 m³/s, water starts overbank flow towards the floodplain of the Indus River basin. When the flood discharge reaches up to 14,000 m³/s, water starts to overbank flow towards its floodplain. The Chenab River diverts water into the lower reach of the Indus River, which increases the flood magnitude of the Indus River and causes significant damage downstream.

2.2. Datasets Used in the Study

Inundation modeling through hydrodynamic simulation requires detailed and finer resolution information regarding soil type, land use, and topography [112,113,114,115]. Due to the lack of availability of cloud-free Landsat images, the daily MOD09GA product, as presented in Table 1, was used for land use classification of the study area. The normalized difference vegetation index (NDVI) initially proposed by [116], as presented in Equation (1), was used for land use classification through supervised classification. The land use information was further used to assign the n values, representing resistance to the flow of water over the floodplain, with the help of literature [117,118]. Detailed soil type information drawn from the Harmonized World Soil Database v 1.2 with a spatial resolution of 1 km [119] was used in the inundation modeling, which also plays an important role in hydrologic and hydrodynamic modeling of rivers’ basins [120,121].

$$\mathrm{NDVI}=\frac{\mathrm{NIR}-\mathrm{RED}}{\mathrm{NIR}+\mathrm{RED}}$$

(1)

2.3. Numerical Simulation in the HEC–RAS Model

The two-dimensional inundation modeling needs detailed topographic information of rivers and floodplains [52] to simulate flood propagation over the floodplain, systematically [122,123]. In this study, the radiometrically terrain corrected (RTC) digital elevation model (DEM) [124,125] of phased array type L-band synthetic aperture radar (PALSAR), mounted on the Advanced Land Observing Satellite-1, was used with a spatial resolution of 12.5 m. The DEM was further converted to a continuous surface triangulated irregular network (TIN) to extract rivers’ cross sections and floodplain geometry. The TIN was further processed in HEC-GeoRAS in ArcGIS software for the extraction of rivers’ parameters such as flow paths, rivers’ banks, cross sections, etc. All the parameters were then exported from HEC-GeoRAS to be used in the HEC–RAS model. All the values of Manning’s roughness coefficient were applied to each cross section. The inflow hydrograph at both flow gauges was used as boundary conditions to assign the flow to cells covering the complete river width at 1st cross-sections of the Indus and Chenab rivers. HEC–RAS v5 performs inundation simulation by solving two-dimensional Saint-Venant / diffusive wave equations [126,127]:

$$\frac{\partial \mathsf{\zeta }}{\partial \mathrm{t}}+\frac{\partial \mathrm{p}}{\partial \mathrm{x}}+\frac{\partial \mathrm{q}}{\partial \mathrm{x}}=0$$

(2)

$$\frac{\partial \mathrm{p}}{\partial \mathrm{t}}+\frac{\partial }{\partial \mathrm{x}}\left(\frac{{\mathrm{p}}^2}{\mathrm{h}}\right)+\frac{\partial }{\partial \mathrm{y}}\left(\frac{\mathrm{pq}}{\mathrm{h}}\right)=-\frac{{\mathrm{n}}^2\mathrm{pg}\sqrt{{\mathrm{p}}^2+{\mathrm{q}}^2}}{{\mathrm{h}}^2}-\mathrm{gh}\frac{\partial \mathsf{\zeta }}{\partial \mathrm{x}}+\mathrm{pf}+\frac{\partial }{\mathsf{\rho }\partial \mathrm{x}}\left({\mathrm{h}\mathsf{\tau }}_{\mathrm{xx}}\right)+\frac{\partial }{\mathsf{\rho }\partial \mathrm{y}}\left({\mathrm{h}\mathsf{\tau }}_{\mathrm{xy}}\right)$$

(3)

$$\frac{\partial \mathrm{q}}{\partial \mathrm{t}}+\frac{\partial }{\partial \mathrm{y}}\left(\frac{{\mathrm{q}}^2}{\mathrm{h}}\right)+\frac{\partial }{\partial \mathrm{x}}\left(\frac{\mathrm{pq}}{\mathrm{h}}\right)=-\frac{{\mathrm{n}}^2\mathrm{qg}\sqrt{{\mathrm{p}}^2+{\mathrm{q}}^2}}{{\mathrm{h}}^2}-\mathrm{gh}\frac{\partial \mathsf{\zeta }}{\partial \mathrm{x}}+\mathrm{qf}+\frac{\partial }{\mathsf{\rho }\partial \mathrm{y}}\left({\mathrm{h}\mathsf{\tau }}_{\mathrm{yy}}\right)+\frac{\partial }{\mathsf{\rho }\partial \mathrm{x}}\left({\mathrm{h}\mathsf{\tau }}_{\mathrm{xy}}\right)$$

(4)

where h is the depth of water (m), $\mathrm{g}$ is the acceleration due to the gravity (m s⁻²), p and q are the specific flow (m² s⁻¹) in cartesian directions, n is the Manning resistance, $\mathsf{\zeta }$ is the surface elevation (m), $\mathsf{\rho }$ is the density of water (kg m⁻³), $\mathrm{f}$ is Coriolis (s⁻¹), and ${\mathsf{\tau }}_{\mathrm{xx}}$, ${\mathsf{\tau }}_{\mathrm{yy}}$ and ${\mathsf{\tau }}_{\mathrm{xy}}$ are components of effective shear stress [127]. The initial terms of the Saint Venant Equations (3) and (4) were ignored for the diffusive wave equations. Both equations were applied with all possible combinations of n values for rivers (main channels) and floodplains, for simulation purposes.

2.4. Satellite-Based Flood Extent Mapping

The satellite datasets used for water bodies delineation in the study are also presented in Table 1. The geographic projection of MODIS and Landsat was transformed into a universal transverse Mercator projection. The NDWI equation, initially proposed by [128], was applied as presented below:

$$\mathrm{NDWI}=\frac{\mathrm{Green}-\mathrm{NIR}}{\mathrm{Green}+\mathrm{NIR}}$$

(5)

Flood water contains a higher amount of sediments, therefore, a different modified NDWI (MNDWI) equation, initially proposed by [129], was also needed for flood mapping. MNDWI and MNDWI2 have also been used for flood mapping by researchers [130,131,132], as described below:

$$\mathrm{MNDWI}1=\frac{\mathrm{Green}-\mathrm{SWIR}1}{\mathrm{Green}+\mathrm{SWIR}1}$$

(6)

$$\mathrm{MNDWI}2=\frac{\mathrm{Green}-\mathrm{SWIR}2}{\mathrm{Green}+\mathrm{SWIR}2}$$

(7)

The details and standard referenced (literature-based) threshold values of NDWI, MNDWI 1, and MNDWI 2 used for flood mapping are described in Table 1. The theoretical consideration of threshold indicates flood cells with a value greater than 0 [133]; however, this threshold value of 0 is not practical in all extreme cases [134,135]. By considering the characteristics of water, an adjustment in threshold values is needed according to local conditions [130,136], therefore, the literature-based thresholds were used for the satellite-based flood mapping. All algorithms of water body delineation were applied, and a union method was applied to produce a blended observed flood extent for comparing the HEC–RAS-based simulated flood extent.

Table 1. Satellite datasets and threshold values of NDWI, MNDWI 1, and MNDWI 2.

Satellite Datasets	Method	Threshold Values	References
Landsat 5 TM	NDWI (2,4) MNDWI 1 (2,5)	0.234, 0.205 0.35, 0.45, 0.33	[137,138] [135,137,139]
Landsat 8 OLI	NDWI (3,5) MNDWI 1 (3,6) MNDWI 2 (3,7)	0.113, 0.09 0.286, 0.33, 0.25–0.31 0.462	[137,139] [139,140,141] [140]
MODIS (MOD09GA/MOD09A1)	NDWI (4,2) MNDWI 1 (4,6)	0.0 0.44, 0.34	[142] [143,144]

Note: The satellite captures of floods are presented as supplementary figures (Figures S1 and S2).

Table 1. Satellite datasets and threshold values of NDWI, MNDWI 1, and MNDWI 2.

Satellite Datasets	Method	Threshold Values	References
Landsat 5 TM	NDWI (2,4) MNDWI 1 (2,5)	0.234, 0.205 0.35, 0.45, 0.33	[137,138] [135,137,139]
Landsat 8 OLI	NDWI (3,5) MNDWI 1 (3,6) MNDWI 2 (3,7)	0.113, 0.09 0.286, 0.33, 0.25–0.31 0.462	[137,139] [139,140,141] [140]
MODIS (MOD09GA/MOD09A1)	NDWI (4,2) MNDWI 1 (4,6)	0.0 0.44, 0.34	[142] [143,144]

Note: The satellite captures of floods are presented as supplementary figures (Figures S1 and S2).

2.5. Calibration and Validation of HEC–RAS Model

The calibration and validation of the HEC_RAS v5 model were performed for 2010 and 2015 flood events, respectively. Flow hydrographs were used for the simulation process for both calibration and validation periods. The flow hydrographs had a peak flow of 27,184 and 17,124 m³/s on 2 August 2010, and 5 August 2015 at Taunsa Barrage, respectively. Approximately 950 and 3398 m³/s of flood water was diverted from the Chenab River into the lower reach of the Indus River. Daily MOD09GA (2 August 2010), 8-day MOD09A1 (5 August 2010), and Landsat 5 TM (12 August 2010) were used to delineate the flood extent for 2010. The daily MOD09GA (6 August 2015), 8-day MOD09A1 (13 August 2015), and Landsat 8 OLI (10 August 2015) were used for delineation of the flood extent for 2015. The n values were optimized in the simulation process by comparing the simulated and observed flood extents for both 2010 and 2015. The river routing accuracy was assessed by comparing the simulated and observed magnitudes and times of peak flows at the barrages downstream.

2.6. Accuracy Assessment of HEC–RAS Model

An uncertainty analysis is necessary in flood inundation modeling, for accuracy assessment [145]. The simulated flood cells were compared with the observed flood cells obtained from satellite imagery, for accuracy assessment of the inundation modeling. The accuracy of the simulation was tested using F1 and F2 indicators recommended by [146]:

$$\mathrm{F}1=\frac{\mathrm{A}}{\mathrm{A}+\mathrm{B}+\mathrm{C}}$$

(8)

$$\mathrm{F}2=\frac{\mathrm{A}-\mathrm{B}}{\mathrm{A}+\mathrm{B}+\mathrm{C}}$$

(9)

The Sørensen–Dice coefficient, also known as Dice similarity coefficient (DSC), is the quotient of similarity, and was also used to check the similarity between the simulated and observed flood areas.

$$\mathrm{DSC}=\frac{2\mathrm{A}}{2\mathrm{A}+\mathrm{B}+\mathrm{C}}$$

(10)

The Jaccard distance (JD) was also used to measure dissimilarity between simulated flood areas and satellite-based (observed) flood areas.

$$\mathrm{JD}=\frac{\left(\mathrm{A}+\mathrm{B}+\mathrm{C}\right)-\mathrm{A}}{\mathrm{A}+\mathrm{B}+\mathrm{C}}$$

(11)

where, A is the accurately simulated flooded area (km²) in the satellite imagery, B is the simulated flood area (km²) which is not present (observed) in the satellite imagery, and C is the area (km²) which is not simulated as flood but present (observed) as a flood area in the satellite imagery. The JD ranges from 0 to ∞, and values close to 0 indicate less dissimilarity between simulated and observed flood areas. The F1, F2, and DSC values range from 0 to 1, −1 to 1, and 0 to 1, respectively, where the values close to 1 indicate the best agreement between simulated and the observed flood areas.

3. Results

3.1. Land Use of the Study Area

The major land uses of the floodplains include the production of cotton, rice, sugarcane, maize, and other fodders. The overall accuracy of the land use map of 2010 and 2015, was 64.06 and 61.64%, with 66 and 73 sampling points, respectively. The n values for floodplains against land use information used in the HEC–RAS model are presented in Figure 2. The estimation of n values was mainly dependent on the remote-sensing-based land use information of the floodplains. Propagation of the flood extent in inundation models mainly depends upon the surface roughness offered by the features of the floodplains.

3.2. Performance Evaluation of HEC–RAS Model

Flood inundation simulation was performed for the 2010 flood against the input hydrograph at the very first HEC_RAS station of the study area. The simulation was carried out by assigning n values for the floodplains and the main channels (Indus River and its tributary, the Chenab River). The floodplains’ and main channels’ n values varied from 0.3 to 0.7 for simulation purposes, as presented in

Table 2. Combinations of Manning’s roughness coefficient (n) values of floodplains and main channels for simulation during the calibration period 2010.

Parameter	Floodplain n	Main Channels (Rivers) n
		0.03	0.035	0.04	0.045	0.05	0.055	0.06	0.065	0.07
A	0.03	3626	3716	3878	4029	4129	4167	4172	4172	4172
	0.04	3699	3745	3955	4108	4138	4159	4165	4165	4165
	0.05	3727	3809	3992	4148	4253	4291	4329	4329	4329
	0.06	3817	3886	4034	4217	4329	4461	4461	4461	4461
	0.07	3847	3963	4087	4265	4388	4461	4461	4461	4461
B	0.03	325	412	438	512	496	588	601	601	601
	0.04	447	511	541	601	638	622	624	624	624
	0.05	487	539	561	649	678	699	701	701	701
	0.06	512	586	635	697	714	758	816	816	816
	0.07	529	629	667	758	801	808	816	816	816
C	0.03	1412	1322	1160	1009	909	871	866	866	866
	0.04	1339	1293	1083	930	900	879	873	873	873
	0.05	1311	1229	1046	890	785	747	709	709	709
	0.06	1221	1152	1004	821	709	577	577	577	577
	0.07	1191	1075	951	773	650	577	577	577	577
F1	0.03	0.68	0.68	0.71	0.73	0.75	0.74	0.74	0.74	0.74
	0.04	0.67	0.67	0.71	0.73	0.73	0.73	0.74	0.74	0.74
	0.05	0.67	0.69	0.72	0.73	0.75	0.76	0.76	0.76	0.76
	0.06	0.69	0.70	0.71	0.74	0.76	0.77	0.76	0.76	0.76
	0.07	0.69	0.70	0.72	0.74	0.75	0.76	0.76	0.76	0.76
F2	0.03	0.62	0.61	0.63	0.63	0.66	0.64	0.63	0.63	0.63
	0.04	0.59	0.58	0.61	0.62	0.62	0.62	0.63	0.63	0.63
	0.05	0.59	0.59	0.61	0.62	0.63	0.63	0.63	0.63	0.63
	0.06	0.60	0.59	0.60	0.61	0.63	0.64	0.62	0.62	0.62
	0.07	0.60	0.59	0.60	0.61	0.61	0.62	0.62	0.62	0.62
DSC	0.03	0.81	0.81	0.83	0.84	0.85	0.85	0.85	0.85	0.85
	0.04	0.81	0.81	0.83	0.84	0.84	0.85	0.85	0.85	0.85
	0.05	0.81	0.81	0.83	0.84	0.85	0.86	0.86	0.86	0.86
	0.06	0.81	0.82	0.83	0.85	0.86	0.87	0.86	0.86	0.86
	0.07	0.82	0.82	0.83	0.85	0.86	0.87	0.86	0.86	0.86
JD	0.03	0.32	0.32	0.29	0.27	0.25	0.26	0.26	0.26	0.26
	0.04	0.33	0.33	0.29	0.27	0.27	0.27	0.26	0.26	0.26
	0.05	0.33	0.32	0.29	0.27	0.26	0.25	0.25	0.25	0.25
	0.06	0.31	0.31	0.29	0.26	0.25	0.23	0.24	0.24	0.24
	0.07	0.31	0.30	0.28	0.26	0.25	0.24	0.24	0.24	0.24

. The simulation performance of the HEC–RAS model for the calibration period indicated that the simulated flood area was increased by the increasing n values of the main rivers and floodplains, while the effect of the n value of the main channels was dominant. It is also evident from parameter A in Table 2 that the small increase in n value of the main channel caused a higher increase in the simulated flood extent, while increasing more floodplains’ n values caused a small increase in the simulated flood extent. The observed flood area from the blended satellite-based flood extent was 5038 km². The lowest correctly simulated flood extent was recorded as 3626 km² with an n value of 0.03 for the main channels and floodplains. The highest correctly simulated flood extent was recorded as 4461 km² with n values of 0.055 to 0.07, and 0.06 to 0.07, for main channels and floodplains, respectively. The JD values indicated that the dissimilarities were smaller throughout the simulation process. The highest value of JD was 0.33 at a few smaller combinations of n values of main channels and floodplains, while the lowest value of JD was 0.23, which corresponded to the n values of 0.055 and 0.06 for the main channels and floodplains, respectively. It is also evident from Table 1 that the decrease in JD was higher with a smaller increase in n values of the main channels, while the decrease in JD was smaller even with a large increase in n value of the floodplains. The lowest values of F1, F2, and DSC corresponded to the lower n values of the floodplains and main channels. The highest values (nearest to 1) of F1, F2, and DSC, were 0.77, 0.64, and 0.87, which corresponded to the n values of 0.055 and 0.06, for main channels and floodplains, respectively. Higher values of F1 and F2 for inundation modeling were demonstrated by the researchers [147].

The numerical simulation of 2D diffusive wave equations to predict simulated flood cells were compared with the observed flood extents, as presented in Figure 3a,b. There was a good relationship between the simulated and observed flood extents near the river and in the lower part of the study area. There were some dissimilarities in the simulated and observed flood extents on the floodplains at the left bank of the river. It is evident from Figure 3a that there was a water depth of less than 3 m in most parts of the study area, while there were some areas where the water depth ranged from 3 to 4.5 m. The most difficult/problematic situation was in region of Rajan Pur, and the upper parts of Rahimyar Khan and Muzaffargarh districts, where water depth ranged from 3 to 6 m.

The simulation started with an input hydrograph at the first station for the validation period 2015, and the results are presented in

Table 3. Combinations of Manning’s roughness coefficient (n) values of floodplains and main channels for simulation during the validation period 2015.

Parameter	Floodplain n	Main Channels (Rivers) n
		0.03	0.035	0.04	0.045	0.05	0.055	0.06	0.065	0.07
A	0.03	2472	2584	2628	2653	2715	2766	2791	2791	2791
	0.04	2585	2643	2679	2721	2773	2816	2828	2828	2828
	0.05	2642	2666	2757	2768	2809	2866	2907	2907	2907
	0.06	2709	2761	2804	2817	2848	2955	2996	2996	2996
	0.07	2751	2816	2862	2885	2909	2979	2996	2996	2996
B	0.03	178	197	215	239	269	305	321	321	321
	0.04	191	210	243	296	325	366	406	406	406
	0.05	212	247	281	313	347	394	416	416	416
	0.06	248	286	335	368	386	416	416	416	416
	0.07	279	323	347	381	402	426	446	446	446
C	0.03	830	718	674	649	587	536	511	511	511
	0.04	717	659	623	581	529	486	474	474	474
	0.05	660	636	545	534	493	436	395	395	395
	0.06	593	541	498	485	454	347	306	306	306
	0.07	551	486	440	417	393	323	306	306	306
F1	0.03	0.71	0.74	0.75	0.75	0.76	0.77	0.77	0.77	0.77
	0.04	0.74	0.75	0.76	0.76	0.76	0.77	0.76	0.76	0.76
	0.05	0.75	0.76	0.77	0.77	0.77	0.78	0.79	0.79	0.79
	0.06	0.76	0.77	0.77	0.77	0.78	0.80	0.81	0.81	0.81
	0.07	0.77	0.78	0.78	0.78	0.79	0.80	0.80	0.80	0.80
F2	0.03	0.66	0.68	0.69	0.68	0.68	0.68	0.68	0.68	0.68
	0.04	0.69	0.69	0.69	0.67	0.67	0.67	0.65	0.65	0.65
	0.05	0.69	0.68	0.69	0.68	0.67	0.67	0.67	0.67	0.67
	0.06	0.69	0.69	0.68	0.67	0.67	0.68	0.69	0.69	0.69
	0.07	0.69	0.69	0.69	0.68	0.68	0.68	0.68	0.68	0.68
DSC	0.03	0.83	0.85	0.86	0.86	0.86	0.87	0.87	0.87	0.87
	0.04	0.85	0.86	0.86	0.86	0.87	0.87	0.87	0.87	0.87
	0.05	0.86	0.86	0.87	0.87	0.87	0.87	0.88	0.88	0.88
	0.06	0.87	0.87	0.87	0.87	0.87	0.89	0.89	0.89	0.89
	0.07	0.87	0.87	0.88	0.88	0.88	0.89	0.89	0.89	0.89
JD	0.03	0.29	0.26	0.25	0.25	0.24	0.23	0.23	0.23	0.23
	0.04	0.26	0.25	0.24	0.24	0.24	0.23	0.24	0.24	0.24
	0.05	0.25	0.25	0.23	0.23	0.23	0.22	0.22	0.22	0.22
	0.06	0.24	0.23	0.23	0.23	0.23	0.21	0.19	0.19	0.19
	0.07	0.23	0.22	0.22	0.22	0.21	0.20	0.20	0.20	0.20

. The satellite-based (observed) flood extent was 3302 km², and the lowest and highest correctly simulated flood extents were recorded as 2472 and 2996 km², respectively, with the same corresponding n value combinations of main channels and floodplains. The lowest values of F1, F2, and DSC were 0.71, 0.66, and 0.86, respectively, which also corresponded to the lower n values of floodplains and main channels and calibration period. The highest values (nearest to 1) of F1, F2, and DSC, were 0.81, 0.69, and 0.89, which also corresponded to n values of 0.06, for the main channels and floodplains. The higher values (close to 1) of F1, F2, and DSC indicated good validation results. The highest value of JD was 0.29 at a few smaller combinations of n values of main channels and floodplains, while the lowest value of JD was 0.19, which corresponded to the n values of 0.06 for the main channel and floodplains. The validation results were closely related to the calibration results; it is evident from Table 1 that the decrease in JD was also higher with a slight increase in the n values of the main channels. In contrast, the decrease in JD was smaller even with a large increase in the n values of the floodplains.

The comparison of simulated and observed flood areas is presented in Figure 4a,b. As discharge was not much higher than the calibration period, the simulated and observed flood extents were smaller. There was a similarity in the observed and simulated flood extents in almost all parts of the study area, except the upper parts of the Indus and Chenab rivers’ confluence, where the modeled flood area was more, and the observed flood area was less. The water depth in most parts of the study area was less than 3 m, while some areas received a water depth of 3.5–4.5 m. Similarly, for the calibration period, similar parts of the study area received a flood depth of more than 4.5 m, but the extent of those areas was less than the calibration period. The slightly lower values of F1, F2, and DSC, and slightly higher values of DJ for calibration and validation period, were closely related to the resolution of topographic information which was used in the study. For the calibration period, the coarser resolution MOD09GA product was available the very next day, but Landsat 5 TM was only available 10 days after the flood date. For the validation period 2015, Landsat 8 OLI was available 5 days after the flood date and MOD09GA was available the very next day after the flood date. The MOD09A1 product was available after 3 and 8 days, for the calibration and validation period, respectively. The delay in high resolution Landsat images caused the misclassification of flood cells used in calibration and validation periods. The misclassification of flood cells, along with other topographic, hydrometeorological and roughness factors, caused a slight reduction in F1, F2, and DSC for the calibration and validation periods.

4. Discussion

The slightly lower values of F1, F2, and DSC, and slightly higher values of DJ were due to several errors associated with the satellite imagery and simulation in HEC–RAS. The overall accuracy of the land use map for 2010 and 2015 was not higher. The reason may be due to the small number of sampling data and poor resolution of the MOD09GA product, wavelength, polarization, and incident angle. The n values of the floodplain and main channels, estimated using land use information, play a key role in the propagation of floods over the floodplain, by offering resistance to the flood propagation. The n values within the floodplains are not a fixed quantity; however, they are a scale-dependent parameter representing all energy losses, which may also include drag force due to the presence of bridges. This highlights the significance of land use in hydrodynamic modeling, as dense vegetation surface offers significant friction to floodwater, and the least friction is offered by unvegetated areas [148]. Hence, the spatial variation of n values corresponding to the estimated land use information in the floodplains directly affects the parametrization of n values, and the accuracy of inundation modeling [49,149]. The accuracies of flood inundation modeling, water depth and flood extent are directly influenced by TIN, which is developed using topographic information of the rivers and floodplains [150,151,152]. This highlights the significance of river cross sections and floodplain geometry, which are mainly dependent on the TIN, indicating that the use of RTC DEM can improve the accuracy of hydrodynamic modeling. However, the accuracy of flood inundation simulation modeling can be further enhanced by using higher resolution topographic information of the study area [153]. It was also observed, during the simulation, that propagation of the flood extent was more dependent on the water surface elevation in the main channel, compared with the floodplain surface roughness, which was also observed by other researchers [147]. The wetness and dryness patterns of the floodplains play a key role in the propagation of flood extent [154], as the simulated flood extent and its propagation was directly related to parametrized n values, soil moisture conditions, in situ rainfall, and wind conditions of the floodplains [147]. These scientific findings revealed that the use of in situ or satellite-based soil moisture, rainfall, and wind conditions can be further integrated with hydrodynamic models to increase the overall accuracy of inundation modeling.

An accurate assessment of flood extent can be obtained from in situ observation and field surveys, which is not possible during the flood event. Satellite imagery is the only alternative source for flood mapping during and after the hazard [155,156,157], which can be utilized as a standard product for flood inundation modeling [158]. The flood inundation models have been calibrated [159] and validated [160] by the researchers with SAR satellite imagery. The time delay between the flood date and satellite imagery acquisition is one of the major factors associated with slightly lower accuracy of flood inundation modeling [161]. The flood extent delineation from satellite imagery is also sensitive to the date of acquisition, wavelengths of sensors, polarization, incident angle, and prevailing atmospheric conditions [147]. There is a need to further improve the accuracy of flood inundation modeling in the central part of the Indus River basin, which can be attained by incorporating advanced three-dimensional topographic information [162,163,164,165,166], and distributed friction parameters based on vegetation heights [167,168]. Significant advancements and higher accuracy can be further attained in hydrodynamic modeling by treating topographic features explicitly at grid cell scale, and homogeneous features at sub-grid scale [167].

The low-lying areas of the floodplain exhibited good performance compared with the high-lying areas, and reproduced the extreme event successfully. It is recommended that inundation modeling should always be conducted in conjunction with detailed field survey data, for proper validation of results. Probabilistic flood mapping methodologies [145,169], hybrid stochastic approaches [170], and coupled hydrological–hydrodynamic modeling [171] are also recommended to minimize the uncertainties associated with inundation modeling.

5. Conclusions

This hydrodynamic modeling of the central part of the Indus River basin proved to be helpful for indicating the most flood-vulnerable areas. The overall accuracy of the land use map for 2010 and 2015, was 64.06 and 61.64%, with 66 and 73 sampling points, respectively. The calibration (2010) and validation (2015) of the HEC–RAS model resulted in 0.055 and 0.06 optimized n values for the main channels and floodplains, respectively. The model simulation accuracy at optimized n values was reliable. The DSC, F1, and F2, were 0.77, 0.64, and 0.87 for the calibration period, and 0.71, 0.66, and 0.86 for the validation period, respectively, The JD values were 0.23, and 0.19 for the calibration and validation periods, respectively. The increase in simulated flood area was higher, with a small increase in n values of the main channels (rivers), while an increase in the modeled flood area was less, with a higher increase in floodplains’ n values. From the results, it can be concluded that HEC–RAS inundation modeling needs proper calibration and validation with in situ flood extent, which may not be possible during a catastrophic flood event. The MODIS and Landsat-based flood extent proved helpful for the comparison of simulated flood extent at different combinations of n values of the main channels and floodplains, and ultimately for the optimization of n values. The reliable results of the calibration and validation periods indicate a good agreement exists between simulated and observed flood extents at most of the places of the floodplains. The HEC–RAS v5 model can be helpful for early flood warning if properly calibrated and validated in the large river basins of the world which have extreme hydrometeorological and complex topographic conditions.