Two-Dimensional Flood Inundation Modeling in the Godavari River Basin, India—Insights on Model Output Uncertainty

: Most ﬂood inundation models do not come with an uncertainty analysis component chieﬂy because of the complexity associated with model calibration. Additionally, the fact that the models are both data- and compute-intensive, and since uncertainty results from multiple sources, adds another layer of complexity for model use. In the present study, ﬂood inundation modeling was performed in the Godavari River Basin using the Hydrologic Engineering Center—River Analysis System 2D (HEC-RAS 2D) model. The model simulations were generated for six different scenarios that resulted from combinations of different geometric, hydraulic and hydrologic conditions. Thus, the resulted simulations account for multiple sources of uncertainty. The SRTM-30 m and MERIT-90 m Digital elevation Model (DEM), two sets of Manning’s roughness coefﬁcient (Manning’s n) and observed and estimated boundary conditions, were used to reﬂect geometric, hydraulic and hydrologic uncertainties, respectively. The HEC-RAS 2D model ran in an unsteady state mode for the abovementioned six scenarios for the selected three ﬂood events that were observed in three different years, i.e., 1986, 2005 and 2015. The water surface elevation (H) was compared in all scenarios as well as with the observed values at selected locations. In addition, ‘H’ values were analyzed for two different structures of the computational model. The average correlation coefﬁcient (r) between the observed and simulated H values is greater than 0.85, and the highest r, i.e., 0.95, was observed for the combination of MERIT-90 m DEM and optimized (obtained via trial and error) Manning’s n. The analysis shows uncertainty in the river geometry information, and the results highlight the varying role of geometric, hydraulic and hydrologic conditions in the water surface elevation estimates. In addition to the role of the abovementioned, the study recommends a systematic model calibration and river junction modeling to understand the hydrodynamics upstream and downstream of the junction.


Introduction
Flooding is one of the natural hazards that is observed globally; however, flood events have different aspects and generating mechanisms, and localized geomorphological processes such as erosion and sediment deposition along the river play a key role in flood inundation, and consequently, flood impacts [1]. Flood inundation modeling, which typically consists of modeling of flow aspects due to the above complex processes, requires detailed information so that the models can run at the required scale [2]. With the advent of efficient numerical methods in combination with high-performance computing resources in recent times, flood inundation models have made a significant leap in their modeling capabilities, including their capability for producing reliable and accurate estimates of various flow aspects [3]. However, relatively, less research has been pursued in the context of uncertainties that associated with the model estimates [4][5][6][7]. A potential review on the outlet of the Sabari sub-basin, to Kunavaram, where the river Sabari confluences with the Godavari River. The 35 km stretch on the Sabari River considered for the followings reasons, i.e., the Sabari River is the only major tributary that joins the Godavari on its 240 km stretch and its confluence at the river junction Kunavaram generates turbulence and forms backwaters on the Sabari River during very high flows [personal communication with the Central Water Commission (CWC) of India personnel] [80].
The research presented in this paper aims to provide insights on the uncertainty in the model output from six scenarios, which are formed with Shuttle Radar Topographic Mission (SRTM)-30 m-and Multi Error Removed Improved Terrain (MERIT)-90 m digital elevation models (DEMs), initial and optimized Manning's n values and observed and estimated flow boundary conditions. In this study, DEM based cross-sections compared with the measured cross-section. For the six scenarios, we analyzed the error in the water surface elevation at each computational cell and compared stage hydrographs of selected events for locations of upstream and downstream of the river junction. Also, the calculations were done for different model structures. The results obtained from the above analysis assisted in providing the model output uncertainty.

Study Area
The Godavari River flows for 1465 km from west to east and drains an area of 312,812 km 2 , which is about 10% of the total geographical area of India. The Godavari River Basin (GRB) is categorized into upper, middle and lower river basins, and the Lower Godavari Basin (LGB) occupies nearly 28% of the area of the GRB. Two stretches of rivers of 240 km and 35 km length on the Godavari River and Sabari River, respectively, were considered for this study, and these vary in their topographic and hydrologic aspects. The daily average flows (peak water surface elevations) are 465 m 3 Figure 1 shows the location of the study area.

Data
Two digital elevation models (DEMs), i.e., the Multi Error Removed Improved Terrain-90 m (MERIT-90 m) [81] and the Shuttle Radar Topographic Mission-30 m (SRTM-30 m) [82], were used to develop the two-dimensional geometrical mesh of the study area. Land use and land cover (LULC) data with a spatial resolution of 100 m for the years 1985 and 2005 were obtained from the Oak Ridge National Laboratory Distributed Active Archive Center (ORNL DAAC), whereas LULC data with a spatial resolution of 56 m for the year 2011 were obtained from the Indian Space Research Organization (ISRO)'s geoportal Bhuvan [83]. The study area can be categorized into 9 and 12 different land use and land cover classes for 1985 and 2005, and 2011, respectively, based on their land use. Manning's n values were obtained from both the National Land Cover Dataset (NLCD) [84] and the HEC-RAS manual [85]. The observed data of hydraulic variables, such as maximum velocity (V), wetted perimeter (P), Manning's n and bottom frictional slope (S O ), were obtained for each event on a daily time scale from the Krishna and Godavari Board Organization (KGBO), Hyderabad, India. The flow hydrographs at Perur and Konta were used as upstream boundary conditions, and the stage hydrograph at Polavaram was used as a downstream boundary condition for the HEC-RAS 2D model. For the Konta GD location, both the observed and estimated flow hydrographs were considered. The estimated streamflows were obtained from a calibrated hydrologic model, i.e., the Hydrologic Engineering Center-Hydrologic Modeling System (HEC-HMS), in a semi-distributed mode for the Sabari River Basin. We considered three large flood events, of which each event comes from a separate year, and the characteristics of the flood events are shown in the Table 1. Note that the start and end of the event is subjective, and the dates were selected after a detailed analysis of the flow hydrographs of the events. Table 1. Duration of the event (t), Peak flow (Q P ), peak water surface elevation (H p ) and time to peak (T p ) for the selected three flood events (E1, E2 and E3) at a boundary condition gauge location. Q*p reflects the estimated values from the calibrated HEC-HMS model. An unstructured computational mesh of grid size 250 m × 250 m was developed for the flow area ( Figure 2). The 2D flow area consists of 52,241 computational cells, of which many cells, i.e.,~95%, are of 250 m resolution; the maximum, minimum and average area of all the cells is approximately 0.14 km 2 , 0.009 km 2 and 0.062 km 2 , respectively. The 2D area break lines are applied to the cells of high elevation zones over which water does not flow to avoid a numerical error; a high elevation zone can be a structure or an object that acts as a barrier to the flow or that controls the flow direction. The boundary condition lines are drawn at Perur and Polavaram on the Godavari River and at Konta on the Sabari River.

Manning's n Values
Manning's roughness coefficients (Manning's n) for each grid cell of the flow area are assigned based on the land use classes of the corresponding grid cell. The LULC maps corresponding to the years 1985, 2005 and 2011 were used for three events, i.e., 1986, 2005 and 2015, respectively. Note that the LULC maps, which require significant effort, are developed for every 5 or 10 years, assuming that the LULC data do not change drastically in a short period of the time except for in a few situations, e.g., deforestation and conversion of agricultural land into non-agricultural land. Therefore , the existing LULC maps of the  years 1985 and 2011 are used for 1986 and 2015, respectively. The study area is represented  with 09-, 09-and 12-land use classes for the years 1985, 2005

Governing Equations
The 2D St. Venant equations, i.e., both the continuity and momentum equations were solved to estimate various flow aspects of unsteady flow in the river. The equations are given below.
Continuity equation: ∂H ∂t where H is the water surface elevation (m), t is the time (s), h is the water depth in the cross-section (m), u and v are the velocity components (m/s) in the x-and y-directions, respectively, and q is the source or sink flux (m/s). Momentum conservation equation: where Equations (2) and (3) represents the conservation of momentum in the Cartesian coordinates [85], g g is the acceleration due to gravity, v t is the horizontal viscosity coefficient, C f is the bottom friction and f is the Coriolis parameter. Parameters of the right side of the equations, v t and C f , can be estimated using the equations below: where D is a non-dimensional empirical constant and has a range from 0.33 to 0.77 [86]. The value 0.33 adopted in this study, and the value corresponds to moderate transversal mixing of turbulence for gentle meanders with moderate surface irregularities; 'h' is water depth at the cross-section and u * is the shear velocity is given by: where g is the acceleration due to gravity, R is the hydraulic radius (m), S is the energy grade line (EGL) slope, |V| is the magnitude of the velocity vector (m/s) and n is Manning's roughness coefficient (s/m 1/3 ). Model stability condition: The Courant Friedrichs Lewy (CFL) condition [87] was used to maintain the model stability: where C is the Courant number, V w is the wave speed (m/s), ∆T is the numerical computation time interval (or time step) (minutes) and ∆X is the distance interval (m). In the above equation, ∆T was used to estimate the wave velocity and the courant number. Values of 'C' provide information on the convergence of the numerical solution and model stability. The C values greater than 3 and 5 for the full-momentum equation (FME) and diffusion wave equation (DWE), respectively, suggest that the model is unstable. Therefore, the selection of ∆T is an important aspect for CFL condition. Initial values for ∆T are assigned as 10 min. The value of ∆X is 250 m. Coriolis force effect: The Coriolis force acts on the motion of all objects that on earth, which rotates around its own axis. Note that the vertical component of the Coriolis force is ignored, and the horizontal component of the force is proportional to the Coriolis parameter (f ). The Coriolis parameter, f, is calculated as below: where ω = 0.00007292115855306587/s, which is the sidereal angular velocity [86] of the earth, and ϕ is the latitude, for which a value of 17 0 is considered.

Model Parameters
The HEC-RAS 2D model needs to be supplied a set of 12 parameters, including the parameters that are mentioned in the above equations. While values of many parameters are initially assigned to a value based on the expertise or guidance such as the HEC-RAS user manual, they need to be modified so that estimated values are in agreement with observed stage values. Horizontal viscosity coefficient, v t , is computed assuming a D value of 0.33, and the EGL slope is calculated separately for the three boundary condition locations and for the three different events using Equations (4) and (6). The EGL slope was computed as 0.000321538, 0.000226667 and 0.000192592 at Perur and 0.000156154, 0.000180769 and 0.000180769 at Konta for E1, E2 and E3, respectively. The initial conditions regarding the time for the model simulations are 3 h, 2.5 h and 2 h for E1, E2 and E3, respectively. The values of V w , a parameter in the courant stability condition, are 2.33 m/s, 1.55 m/s and 2.1 m/s for E1, E2 and E3, respectively. ∆T, a numerical computation time interval, varied from 10 min to 1 min. The implicit weighting factor (θ) corresponds to weights of the spatial derivatives at current and previous time steps, and the value of θ varied from 0.6 to 1 to achieve the convergence to higher values of C. Note that values of various parameters were the same for all the three events.

Model Simulations
The model ran for six different scenarios that resulted from two different DEM sources, two sets of Manning's n and two sets of boundary conditions. The scenario S1 ( Removal of the local and convective acceleration terms including advection and viscous terms from the FME corresponds to the diffusion wave model and the new form labeled as the diffusion wave equation (DWE). The water surface elevation (H) time series for each grid cells of the flow area and inundation boundary (IB) for the event E2 are generated. The model calibration was done by optimizing the implicit weighting factor (θ), time step (∆T) and Manning's n in a trial-and-error mode so that simulations and observed water surface elevation (H) are in close agreement at the selected three locations, i.e., Dummugudem, Bhadrachalam and Koida on the Godavari River.

Model Performance
The model performance was assessed by comparing the observed and simulated water surface elevations (H) at specific locations for different events as mentioned below (Table 2), i.e., Dummugudem and Koida for E1 and E2, and Bhadrachalam for E3. The performance metrics consist of the correlation coefficient (r), Index of agreement (d), Mean Absolute Error (MAE), Root Mean Squared Error (RMSE), Percentage Bias (PBIAS) [88] and Kling-Gupta Efficiency (KGE) [89]. Table 2. Performance metrics used for the model assessment.

Metric Equation Range
Correlation coefficient, r In Table 2, N indicates the number of observations or the duration of the event, O and P are the observed and simulated values, O and P correspond to average observed and simulated values, r is the correlation coefficient, α is a measure of relative variability in the simulated and observed values and β is the ratio of P and O; the square root term in the KGE equation denotes Euclidian distance (ED) from the ideal point.

Model Calibration
The gridded Manning's n values, extracted from NLCD, exhibited a decreased pattern for three events in a chronological order (Figure 3a The model found to be stable for values of ∆T are 1-, 1-and 2-min for E1, E3 and E2 events, respectively. Figure 6 shows the number of computational grid cells for different values of implicit weighting factor (θ) for the M2 scenario of E1, and it assists in understanding the effect of θ on the model stability. The computational cells are the total number of grid cells minus the cells with an error below the specified threshold. These cells varied at 423, 390, 400, 480 and 553 for the theta values from 0.6 to 1 at an interval of 0.1. An increase in values θ first resulted in a decreased (θ from 0.6 to 0.7) and increased (θ from 0.8 to 1) number of computational cells, which means the error in each cell decreased below the specified threshold in the case of θ at 0.6 and 0.7 and increased in the remaining cases. Approximately, 24-and 10 computational cells with solution converged at the Courant number 3 and 4, respectively. For the values of θ = 0.6 and 0.7, approximately 5 and less than 10 cells converged respectively, with a solution. Convergence of the solution for the number of computational cells with the Courant number 2 increased linearly from θ, 0.6 to 1. In addition, the cells correspond to C4, which refers to that the non-convergence of the model solution decreased in the case of θ = 1; the value of θ = 1 yielded a stable model, whereas θ = 0.6 also exhibited similar stability with more values converged at 1. The scenario with θ = 0.6 yield relatively more number of computational cells with acceptable error as compared to θ = 1. However, the scenario is unstable due to the convergence of the solution with C4. Here, C1, C2, C3 and C4 correspond to the solutions that converged with the Courant number 1, 2, 3 and 4, respectively, for the simulations obtained by using the full-momentum equation (FME).

Analysis of H Error
The water surface elevation tolerance is set to its default values, i.e., 0.003 m, after which the error in water surface elevation (H Error ), i.e., the difference between the computed value and assumed value, for each computational cell is calculated, which is then pooled for all computational cells and plotted in a boxplot (Figure 7). The plots are developed for all six scenarios and for each of the three events, and repeated for all three events. Small H Error values (~<0.1 m) were observed for most of the cells (only two outside values were observed) for the M1, S2, M2 and F2 for (E1) and (E2) events. Relatively higher H Error values were observed for S1 and M1 scenarios as compared to rest of the scenarios (median of H Error is 0.05 m). Boxplots of smaller length and without a clear median suggests smaller and left-skewed H Error values for E2 as compared to E1. The values are smaller for E1 and E2, and relatively higher for E3, for all scenarios. Nevertheless, H Error for 75% of the cells is well below 0.2 m, and the upper whisker, which corresponds to~99.5 percentile, is below 0.4 m. The largest errors (outsiders of the boxplots, which correspond to less than 1% of the cells) vary with each configuration and event, and the values are smaller for event E2 as compared to E1 and E3. MERIT

Comparison of Stage Hydrographs
Both observed and estimated stage hydrographs were compared for locations Dummugudem and Koida for E1 and E2, and for Bhadrachalam for E3. Figure 8 has five plots, of which each plot corresponds to a location and an event. Observed and estimated stage values from different scenario are denoted in different colors. Estimated stage hydrographs from all scenarios suggest realistic simulated hydrographs for all events. Nevertheless, overestimated stage values were observed for two events, E1 and E2, at both locations, Dummugudem and Koida, for all scenarios except for four at Dummugudem for E2; these four scenarios consist of optimized Manning's n values but differ in other sources of uncertainty. Irrespective of whether stages were over-or underestimated, the estimated values exhibited a systematic bias for the entire flood event, whereas the performance of event E3 differs slightly, i.e., under-and overestimated values can be observed on the rising and falling limb, respectively. Additionally, it is observed that the scenario with NLCD-based Manning's n values produced hydrographs of higher magnitude as compared to scenarios with optimized Manning's n values.

Comparison of Flood Inundation Maps for the Event E2
Map of inundation boundary (IB) generated for all the six scenarios, for 22 September 2005, and compared with the flood extent estimated from the RADARSAT-1 imagery shown in the [67] for a downstream stretch of the Sabari River. Figure 9 shows maps of IB obtained from this study for six scenarios. While the maps of IB from all scenarios appear to be not deviating significantly, closer look of the same suggests the contradictory. High-and low-areas of inundation observed for scenarios that use NLCD-based Manning's n values and combination of SRTM-30 m and optimized Manning's n values, respectively. The disagreement in the area can also be due to the different time stamps between satellite image coverage and model simulation.

Comparison of Flood Inundation Maps for the Event E2
Map of inundation boundary (IB) generated for all the six scenarios, for 22 September 2005, and compared with the flood extent estimated from the RADARSAT-1 imagery shown in the [67] for a downstream stretch of the Sabari River. Figure 9 shows maps of IB obtained from this study for six scenarios. While the maps of IB from all scenarios appear to be not deviating significantly, closer look of the same suggests the contradictory. Highand low-areas of inundation observed for scenarios that use NLCD-based Manning's n values and combination of SRTM-30 m and optimized Manning's n values, respectively. The disagreement in the area can also be due to the different time stamps between satellite image coverage and model simulation.

Model Performance
The model performance analyzed by calculating various performance metrics for all six different scenarios for two locations, i.e., Dummugudem and Koida, for events E1 and E2, and Bhadrachalam for event E3 (Tables 3-5). Table 3. Model performance metrics estimated by comparing simulated and observed water surface Figure 1 with six scenarios. The results are presented for the upstream location-Dummugudemand the downstream location-Koida. For event E1 and for location Dummugudem, both scenarios M2 and F1 exhibited relatively low error, and consequently, small error verification metrics. However, only the F1 scenario exhibited higher values of goodness of fit measures. The scenarios F1 and F2, where 50% of the peak flow of E1 was underestimated showed a similar performance to the other scenarios, where observed flow boundary conditions were used (S1 to M2). This suggesting that the model structure is unable to simulate the stages obtained from the estimated flow boundary conditions. In addition, the estimated flow boundary conditions used for E2 and E3 are in better agreement with the observed flows (Table 1). For event E2, all scenarios exhibited similar performance; however, S1 configuration exhibited lower error verification metrics. For both events E1 and E2, and for location Koida, F1 exhibited better performance among all in terms of error and goodness of fit measures. In addition, relatively increased error, and consequently, high values of error verification metrics, were observed for Koida. In general, configuration M1 exhibited higher error verification metrics, and better performance was observed for the configuration that used optimal Manning's n values. For event E3 and for location Bhadrachalam, all scenarios exhibited a similar performance, and scenarios M2 and F2 exhibited slightly smaller error verification metrics. In general, simulations of the event E3 at the Bhadrachalam exhibited better performance as compared to the events E1 and E2 at two other locations; note that the event E3 has a relatively standard hydrograph without much variation.

Event Hydrodynamics
To understand the role of model structure in the simulation of the water surface elevation, the model used the FME and DWE for two scenarios, i.e., S2 and M2. The two scenarios, S2 and M2, consist of optimized Manning's n values and observed flow boundary conditions in combination with the two DEM sources. Estimated hydrographs from two different model structures were compared at six locations that were randomly chosen on the main course of the Godavari River starting from upstream to downstream as well as at three GD locations that were mentioned earlier; six locations, i.e., L1, L2, L3, L4, L5 and L6, are shown in Figure 10. Three locations, L1, L2 and L3, are on the upstream of the river junction, and the remaining three locations are on the downstream of the river junction. Figure 10 has stage hydrographs for six locations and for three events, i.e., a total 18 plots. Each plot has four stage hydrographs and each set of two scenarios corresponds to simulations generated from FME and DWE. The plots suggest the following: (i) simulations from each set of the scenarios are of similar values for most of the locations for all events. Thus, no significant role of DEM is found except for location L1 for event E1 and locations L1 and L2 for event E2; (ii) for the event E1, the simulations derived from the DWE are relatively high for the downstream locations, whereas the simulations from the FME are relatively high for the upstream locations; (iii) for the events E2 and E3 and for all the locations, the simulations from FME obtain high values as compared to the simulations from DWE. Thus, the results suggest the importance of the model structure FME vs. DWE and the source of DEM for a few locations.
Furthermore, the estimated stage hydrographs based on FME and DWE were compared with the observed stage hydrograph at Dummugudem and Konta (for E1 and E2) and Bhadrachalam (for E3) GD stations ( Figure 11). For the event E1 and for both locations, peak stage values were overestimated by two scenarios, i.e., DWE and MERIT-90 m and FME and SRTM-30 m, and underestimated by the other two scenarios. However, the estimated stage hydrographs are realistic for location Koida. For the event E2 and for the location Dummugudem, over-and underestimated stage hydrographs were observed for models runs of DWE and FME, respectively, whereas for location Koida, estimated stage hydrographs of all locations are overestimated. For the event E3 and for location Bhadrachalam, stages based on FME closely match with the observed values, and the remaining locations underestimated stages. The performance of the simulations varies with each event and scenarios as well indicating no systematic influence of model structure and DEM source. Additionally, note that these results differ with the results of the comparison of the streamflows at six locations along the main river, i.e., the estimated stage values are of higher values for the model run with FME for many locations for the event E2, and the opposite was observed for the location Koida for E2. Thus, the findings highlight the differences in FME-and DWE-derived estimates and hint at a detailed study.

Insights on the Uncertainty
The HEC-RAS 2D model has no direct uncertainty estimation feature, and this study attempts to calculate the uncertainty in the model output considering different input sources of uncertainty as well as model structure. The two different DEMs, two sets of Manning's n and boundary conditions consist of both observed and estimated flows resulted six different scenarios, whereas the model structure varied treating the conservation of momentum based on either FME or DWE. In this section, we comment on these different sources.
While the SRTM-30 m DEM exhibited a relatively high H Error as compared to the H Error from MERIT-90 m DEM, it is interesting to know that these two DEMs exhibited a value of H Error , <0.1, <0.25, <0.5 and >0.5 m for 85, 95, 99.8 and 0.2 percentage of the HEC-RAS 2D computational grid cells. While the difference is not significant for most of the cells, it is to be understood that the time stamp associated with the DEM also plays a role; for example, the DEMs may not be representative for the year 1986, as these DEMs and the surveyed cross-section profiles correspond to post-major flood event in 1986. A comparison of the DEMs with the measured values also showed the difference; however, instrument-based uncertainty should not be ignored, i.e., the measurements taken in the river course might be affected by the boat movement, which is not in the desired line relative to the point on the river bank. In addition, constant erosion, sediment deposition and changes in the morphology affect the DEM accuracy, and consequently, the calculations.
The scenarios with optimized Manning's n values yield stages of low error as compared to the values that read from the NLCD database, which has values as per the LULC category. This observation suggests that the land use category might not be truly representative in terms of either spatial resolution or, importantly, in terms of temporal dimension. Note that the LULC maps of a spatial resolution 100, 100 and 56 m correspond to the years 1985, 2005 and 2011 used for three events, 1986, 2005 and 2015, respectively. In addition, the season also play a role; for example, the river stream in the LULC map of 2005 (for E2) is invisible, which might be due to the data acquisition done in the dry period. This resulted in considering the high roughness values instead of lower values, which remained same even after optimization. Additionally, the uncertainty associated with Manning's n value for a land use or land cover class, i.e., multiple values or ranges, brings additional uncertainty, particularly in the modeling of extreme events. Estimated stages using Manning's n based on the NLCD data set are of high values for events E1 and E3, whereas the estimated values are in decent agreement with observed stage values for event E2.
Using of either observed or estimated streamflows at Konta on the Sabari River, which is one of the flow boundary conditions, resulted in different scenarios. The model performance in terms of goodness of fit and error metrics suggests the following: (i) relatively, there is a large influence of flows from the main river as compared to the tributary; (ii) there are changes in the velocity with respect to the time and along the channel; the rise in H in the downstream indicates that the increased turbulence is the velocity head added to any downstream change in bed level; (iii) large values of H influenced by rainfall event type, downstream hydraulic structures and river junction.
In addition to the above, uncertainties in estimates of various model parameters need to be explored systematically. These parameters are energy gradient line (EGL) slope, eddy viscosity transverse mixing coefficient (D), implicit weighing factor (θ) and water surface tolerance parameter. The EGL slope can be computed with much less data. The EGL slope and D were explored for optimized values in the trial-and-error mode, but no significant change in the simulations were found. These two parameters are responsible for the velocity or energy head losses; however, the no significant effect suggests a systematic approach to explore the role of the parameters. In addition, convergence of the solution at the Courant number 2 when the error in cells is greater than the prescribed tolerance (0.003 m) impacts the analysis of error.
Similar estimated stages for the two scenarios of a model structure in combination with differences in the estimated stage values correspond to two different model structures, which emphasizes the role of local and convective acceleration terms. The simulated H values when DWE was used for the event E1 showed high peak flow values (>10 m) compared to the simulations using FME at L4 to L6 (Figure 10). This clearly highlights the role of the acceleration terms in simulating stage values. Furthermore, for locations L4 to L6, downstream of the river junction, which have relatively gentle river meanders, exhibited relatively high H values when modeled with DWE; local and convective acceleration terms were ignored. It is interesting to note that the DEM source dominated these terms for the event E2. The estimated H reflects hydrodynamics at six random locations as well as at GD locations, and suggest an increased error at Koida (downstream) as compared to Dummugudem (upstream).

Discussion
It was observed that the flood inundation in the Lower Godavari Basin (LGB) is a function of the event magnitude and topographic aspects, such as elevation and roughness. Changes in the geometric aspects of river, for example contraction and expansion of the river, river meandering as well as river junction, play a significant role in the simulation of water surface elevation (H). In addition, the model structure that calculates the flow aspects also has a significant influence on the model output (stage values). A comparison of crosssections also shows that DEM-based values are significantly missing. Cross-section profile suggests the following: (i) elevation values from two different DEMs are missing for a few locations; (ii) disagreement among the elevation values estimated from DEM sources and the measured values; (iii) large variations in the elevation values in the SRTM-30 m DEM.
A comparison of stage hydrographs from six scenarios indicates the high error for the event E1, which is the largest flood event, and the high error can be attributed to not providing the event rainfall and lack of detailed hydrodynamics. The median of H Error is more for the event E3, which can be due to the land use and elevation changes, causing changes in Manning's n values. Furthermore, for event E3, H values estimated from FME are higher compared to DWE. The estimated values of H for E2 showed the differences with respect to the elevation models ( Figure 11). The error in the estimation of H is higher in the downstream location. The uncertainty can be attributed to the failure of capturing the detailed hydrodynamics for the flood events that resulted from heavy rainfall (e.g., 206.13 mm for E1, 395.46 for E2 and 359.88 mm for E3 in the Sabari sub-basin), and higher streamflow upstream (20,000 m 3 /s to 60,000 m 3 /s in the Godavari River and 9000 m 3 /s to 20,000 m 3 /s in the Sabari River). The error in the H values is influenced by the mixing at the river junction as well as its energy-based junction method used in HEC-RAS 2D. Currently, the model assumes a 1D subcritical flow; however, parameters such as the angle of junction, zone of separation and flow ratio between the main and tributary are not considered. The river junction has its own role that needs to be explored [90][91][92][93][94].
The calibration of three parameters (θ, ∆T and Manning's roughness coefficient) showed an overall improvement in the estimation of H through six scenarios. The calibration of θ for the suggested range (0.6 to 1) for computational cells (error for the cells > tolerance error) exhibited values of error exceeding 1 m. Optimization of ∆T for three events helped in estimating the consistent wave velocity. For the two events (E1 and E3), the optimized Manning's n values showed a greater sensitivity due to the land use change and the uncertainty from the data source and resolution. The LULC data were obtained at different resolutions from two sources (see Section 2.3.2) and the downscaling or aggregation of the data to the DEM's resolution caused a mismatch of land characteristics. Manning's n values from the different sources for the same land use or land cover feature were different; thus, the choice of roughness values and their optimization was also considered as one of the sources of uncertainty.

Conclusions
In this study, a two-dimensional flood inundation modeling using HEC-RAS 2D was performed for six different scenarios that formed based on different sources of error in the model output. The study was carried out for a stretch of the Godavari River from Perur to Polavaram and a stretch of the Sabari River, from Konta, which is the outlet of the Sabari River Basin, to Kunavaram, where the Sabari River meets the Godavari River. The model ran in unsteady mode for the selected three floods events, namely E1, E2 and E3, observed in three different years, namely 1986, 2005 and 2015, respectively. The model output, i.e., water surface elevation was analyzed for six scenarios for the three events, for selected locations. The configuration with SRTM-30 m DEM and Manning's n from NLCD data set yielded high error values and the scenarios that used optimized Manning's n values yielded small errors (E1 and E3). However, all scenarios exhibited a similar performance for event E3 at Bhadrachalam, and the estimated values of the event E3 are in better agreement with the observed stage values as compared to the events E1 and E2. Between the two events, E1 and E2, the estimated values for event E2 are in better agreement with the observed values; moreover, note that the event E1 is the largest event that was observed over the 50-year period of the record. The low performance of event E1 can be attributed to the uncertainty in the DEM, consequently land use classification or corresponding Manning's n values. In addition, the model ran for two different model structures developed using either the full momentum equation or the diffusion wave equation. The estimated values highlighted the differences between the two model structures, i.e., the model that used FME generated relatively high values for events E1 and E3 for the locations along the main course of the Godavari River. However, the opposite was observed for E2 at the GD location Koida.
The study allows to understand the role of the different input sources and model structures on the model output. This study recommends the following: (i) identification of a DEM that more accurately represents the river geometry; (ii) optimization of Manning's n values to reflect the LULC; (iii) analysis of the model with different model structures. While definite conclusions were drawn, a few questions remain to be answered, for example, the different influence of the model structure on the estimated stages at randomly selected locations and GD locations for events E1 and E2; role of the events in the context of varying flow magnitude and initial conditions. In addition, it is a known fact that the river junction plays a significant role, and its role in terms of mixing of flows, turbulence and backwater effects vary with various aspects of the flood event. In this context, the future work may include the following, i.e., increasing the sample of flood events, accounting for various hydraulic aspects such as the eddy viscosity and EGL slope, addressing the backwater effect in the Sabari tributary and separate the modeling of the river stretches before and after the river junction. In addition, combining the model output that is derived from multiple models, for example, LISFLOOD-FP and MIKE is recommended so that uncertainty in the model output will be accounted. In this context, this study provides insights on modeling uncertainties, with the key initiative to develop an ensemble flood inundation modeling for the study area.  Data Availability Statement: Rainfall and flow data can be obtained from the Central Water Commission (CWC) and India Meteorological Department (IMD), respectively. Land use and land cover data for 2011 can be ordered from ISRO-Bhuvan website. The remaining data that is mentioned in the Section 2.2 is available online.